写了一个siGSEA包
最近花了一周的时间写了一个简化版GSEA分析包，命名为siGSEA
前段时间用GSEA处理了一些数据，发现其官网还在更新的只有GSEA的桌面版和服务器java版，其R版本只有一个2005年写的脚本！！！我那时粗略的看了一遍，感觉使用起来实在有点繁琐，并且GSEA官网的R代码很长，足足有2500+行，个人对其代码的风格很不适应；因此就有了想法自己重新写一个GSEA的R包，这样单从使用的角度来说更加方便点（虽然可能自己以后也用不了几次，而且GSEA桌面版功能更加全面。。。），毕竟现在很多生信分析软件都是以R包的形式，当然也顺便写个R包娱乐下
现在Biocondutor上关于GSEA的包我简单的搜了下，有改善了算法的RGSEA，GSEAlm，npGSEA等，有用C写了加快了计算速度的PGSEA，还有用于评估GSEA结果的GSVA，但是从引用率来说GSEA还是热门选手
3天写代码，3天完善改进代码，2天写成R包并测试，几乎花了我所有的空余时间。。。
先记录下一下写代码时遇到的一些问题和感想：
- 写一个代码其实很容易，但是如何能保持代码简洁前提下减少BUG是个花时间的
- 能用apply就少用for循环，能增加一点效率是一点。能用向量化(如rowmeans)解决的问题就不要用apply和for这种循环的方式，前者和后者效率差的不是一点点，尤其的当矩阵数据量很大的时候
- Matrix包的rowMeans的计算效率比base包的rowmeans又能提高不少（一生信技能树小伙伴的提示）
- 当矩阵很大时，用矩阵运算比上面的方法都快（GSEA官网就是通过标准差的公式以矩阵运算来计算）
- 当有些函数跟其他包冲突时，最好先在函数前面加上指定包，比如Matrix::rowMeans
- 将代码写成一个简单的R包形式还是比较容易，但是要其跟biocondutor或者CRAN上的那些正式包来说，还有其他东西需要完善下
- 这个R包没有其他很复杂的算法，都是写简单的逻辑关系，所以写起来比较简单
- 多用一些合适的函数来处理对应的问题，多搜搜总有最佳的方法
- 发现BUG解决BUG真是个体力活！！！
整个R包大约500来行，我留下了主要的几个参数，其他均使用了默认值，并且没写过多的检查和判断输入数据和参数格式的代码，所以整体代码的行数还是比较少的
尽管代码简洁了，但是！！运算速度要比官网的R脚本慢上一点（这还是在数据量较小的情况下），跟GSEA的java版本比，那至少慢了3分之1，数据量越大越明显。。。
有一点遗憾的是，没有对GSEA作图代码做修改，几乎是照搬了其GSEA官网的R代码，等有时间了可能会考虑自己再画一个新的
代码放在了Github上 https://github.com/kaigu1990/siGSEA ，简单的模仿了下别人的样子，至少看起来是那么一回事了，并配上了简单的README说明以及测试数据
下面以一个TCGA的数据为例，再做一个简单的实例
先安装下siGSEA
install_github("kaigu1990/siGSEA")
读入TCGA的fpkm数据，并整理出有癌旁对照样本的表达谱数据，下载方式还是跟之前的一篇博文GSEA-基因富集分析一样
library(dplyr) library(siGSEA) data <- read.table(file = "TCGA-BRCA.htseq_fpkm.tsv", sep = "/t", header = T, row.names = 1, stringsAsFactors = F, check.names = F) normal_df <- data[, str_sub(names(data), 14, 15) %in% 11:19] tumor_df <- data[, grepl(paste(str_sub(names(normal_df), 1, 12), collapse = "|"), names(data)) & !names(data) %in% names(normal_df)] TCGA_data <- cbind(normal_df, tumor_df)
接着考虑到TCGA的ensembl ID有版本号，因此需要预先去除（如果没版本号，这步可不做）
id <- rownames(TCGA_data)
id <- unlist(lapply(id, function(x){
strsplit(x, "//.")[[1]][1]
}))
然后读入GSEA的chip数据，可以在官方FTP上下载 ftp://ftp.broadinstitute.org/pub/gsea/annotations ，我这次选择下载的是ENSEMBL_human_gene.chip，主要用于将Ensembl id转化为gene symbol，从而得到的expr为siGSEA输入的表达谱数据
chip <- read.table(file = "ENSEMBL_human_gene.chip", sep = "/t", header = T, quote = "", stringsAsFactors = F)
id2 <- unlist(lapply(id, function(x){
if (is.na(match(x, chip[,1]))){
x
}else{
chip[match(x, chip[,1]), 2]
}
}))
tmp <- data.frame(symbol=id2, TCGA_data, stringsAsFactors = F, check.names = F)
tmp <- dplyr::distinct(tmp, symbol, .keep_all = TRUE)
rownames(tmp) <- tmp$symbol
expr <- tmp[,-1]
还需要一个分组信息
group <- c(rep("normal", length(names(normal_df))), rep("tumor", length(names(tumor_df))))
class <- factor(group)
以及一个gene set文件，这个也可以在GSEA官网下载 http://software.broadinstitute.org/gsea/downloads.jsp ，选择研究所需的gene set文件即可
gene_set <- readLines("c2.cp.v6.1.symbols.gmt")
运行 simple_gsea
函数进行GSEA分析
res <- simple_gsea(dataexpr = expr,
group = class,
geneset = gene_set,
gsminsize = 15,
gsmaxsize = 500,
nperm = 1000)
运行 gseaplot
函数，出GSEA结果图
gseaplot <- function(expr = data, list = res, gset = NULL, nompval = 0.05, fdr = 0.25)
图片参见github上的demo哦
虽然花了不少时间写了这个简单版R包，但是觉得还是蛮有意思的。。。
本文出自于 http://www.bioinfo-scrounger.com 转载请注明出处
正文到此结束
热门推荐
相关文章
近期评论
-
666
-
文章和留言都翻到11页了 没有OOM
-
-
朋友，翻页到11页，及以后，会出现OOM，无法访问
-
-
-
-
版本号是多少，你可以下载哪个代码仓库，jdk选1.8 直接跑就行
-
-
Loading...
![[HBLOG]公众号](https://www.liuhaihua.cn/img/qrcode_gzh.jpg)
