KeepNotes blog

Stay hungry, Stay Foolish.

0%

写了一个siGSEA包

Posted on Edited on In ,
Symbols count in article: 3.1k Reading time ≈ 3 mins.

最近花了一周的时间写了一个简化版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转载请注明出处

Related Posts