0%

写了一个siGSEA包

最近花了一周的时间写了一个简化版GSEA分析包,命名为siGSEA

前段时间用GSEA处理了一些数据,发现其官网还在更新的只有GSEA的桌面版和服务器java版,其R版本只有一个2005年写的脚本!!!我那时粗略的看了一遍,感觉使用起来实在有点繁琐,并且GSEA官网的R代码很长,足足有2500+行,个人对其代码的风格很不适应;因此就有了想法自己重新写一个GSEA的R包,这样单从使用的角度来说更加方便点(虽然可能自己以后也用不了几次,而且GSEA桌面版功能更加全面。。。),毕竟现在很多生信分析软件都是以R包的形式,当然也顺便写个R包娱乐下 现在Biocondutor上关于GSEA的包我简单的搜了下,有改善了算法的RGSEA,GSEAlm,npGSEA,GSVA等,有用C写了加快了计算速度的PGSEA(应该不仅如此),但是从引用率来说GSEA还是热门选手

3天写代码,3天完善改进代码,2天写成R包并测试,几乎花了我所有的空余时间。。。

先记录下一下写代码时遇到的一些问题和感想:

  1. 写一个代码其实很容易,但是如何能保持代码简洁前提下减少BUG是个花时间的活
  2. 能用apply就少用for循环,能增加一点效率是一点。能用向量化(如rowmeans)解决的问题就不要用apply和for这种循环的方式,前者比后者效率高的不是一点点,尤其的当矩阵数据量很大的时候
  3. Matrix包的rowMeans的计算效率比base包的rowmeans又能提高不少(一生信技能树小伙伴的提示)
  4. 当矩阵很大时,用矩阵运算比上面的方法都快(GSEA官网就是用矩阵运算)
  5. 循环运算前,记得预先设置好数据结构和输出变量的长度和类型,千万别在循环过程中渐进性地增加数据长度
  6. 利用ifelse()语句可以使你的代码更加简便,提升运算速度
  7. 为了提升速度,想过用并行的方法,查了几种方法,都不太满意
  8. 当有些函数跟其他包冲突时,最好先在函数前面加上指定包,比如Matrix::rowMeans
  9. 将代码写成一个简单的R包形式还是比较容易,但是要其跟biocondutor或者CRAN上的那些正式包来说,还有其他东西需要完善下
  10. 这个R包没有其他很复杂的算法,都是写简单的逻辑关系,所以写起来比较简单
  11. 多用一些合适的函数来处理对应的问题,多搜搜总有最佳的方法
  12. 发现BUG解决BUG真是个体力活!!!

整个R包大约500来行,我留下了主要的几个参数,其他均使用了默认值,并且没写过多的代码来检查和判断输入数据和参数格式的代码,所以整体代码的行数还是比较少的

尽管代码简洁了,但是!!运算速度要比官网的R脚本慢上一点(这还是在数据量较小的情况下),跟GSEA的java版本比,那至少慢了3分之1,数据量越大越明显。。。

有一点遗憾的是,没有对GSEA作图代码做修改,几乎是照搬了其GSEA官网的R代码,等有时间了可能会考虑自己再画一个新的

代码放在了Github上https://github.com/kaigu1990/siGSEA,简单的模仿了下别人发布R包的网页样式,至少看起来是那么一回事了,并配上了简单的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 (expr = data, list = res, gset = NULL, nompval = 0.05, fdr = 0.25)

图片参见github上的demo哦

虽然花了不少时间写了这个R包,但是觉得还是蛮有意思的。。。

本文出自于http://www.bioinfo-scrounger.com转载请注明出处