转载

写了一个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 转载请注明出处

原文 http://www.bioinfo-scrounger.com/archives/568

正文到此结束

所属分类：编程技术 Java

本文标签： id https 时间 dist IO git 代码 ftp 图片 ip GitHub 参数 ORM tab 数据软件 ioc 测试安装 list 下载服务器 grep 实例 http BIO java App bug js CTO lib
版权声明： 本文为互联网转载文章，出处已在文章中说明(部分除外)。如果侵权，请联系本站长删除，谢谢。
本文海报： 生成海报一生成海报二

其他链接

关于本站

本站定位：个人技术类博客

本站作用：写博客、记日志、闲聊扯淡鼓捣技术。

问题交流

写了一个siGSEA包

热门推荐

相关文章

说给你听

本文目录

随机标签

书籍教程

近期评论

网站信息

其他链接

关于本站

问题交流