GSVA报错

longgu1314

之前一直尝试做GEO数据芯片的GSVA，奈何一直没理解到MSigDB需要下载那些东西以及如何使用Jimmy大神的脚本，这两天再次看的时候觉得忽然好想能看懂一些了，也可以把代码跑到接近最后的位置了，但是很不幸地报错了......
有哪位大神指导这个是什么原因呢？

金恩12

    感觉在使用软件时,还是应该看一下要使用的软件或者R包的官网的document和发表的相关文献.GSVA和GSEA的不同是,GSVA是先将基因对应到gene set上,得到不同的gene set的得分,然后再根据分组,对gene set的得分进行差异分析.GSEA是先根据分组对基因进行差异分析,再将基因富集到gene set中.所以在gsva分析中需要用的文件只有两个:
    (1)基因表达矩阵(gene expression matrix),SYMBOL或者INTERID格式的均可，但要做均一化处理(如TPM值即可)

    (2)geneset:可以直接从MSigDB下载(一般是.gmt格式的文件)，或者自己制作。下载时，要选择合适的geneID，和表达矩阵保持一致即可。（一般都是SYMBOL格式，因为比较容易供使用者观看）
对于楼主说的不太明白MSiDB中的数据,在官网主页(图1)上(http://software.broadinstitute.org/gsea/msigdb/index.jsp),比如我们需要用KEGG相关的的gene set,选择C2,单击C2进入详细介绍页面,就可以看到我们要下载的文件链接了(图2),记得下载和自己的gene ID 格式相同的gmt文件.下载下来之后,可以看到gmt文件的内容也很单纯(图3),每一行表示一个gene set.第一列是gene set的名字,第三列及之后是每个gene set中含有的gene,所以自己制作gene set时,先对gene set命名然后再添加基因就可以了,但要保证其合理性.
    有了这些之后就可以在R中用GSVA包进行处理了.先导入数据

[AppleScript] 纯文本查看 复制代码

gene <- read.csv("XXX.csv",header = TRUE,sep = ",")
genesets <- getGmt("c2.cp.kegg.v6.1.symbols.gmt")

然后添加分组情况

[AppleScript] 纯文本查看 复制代码

nGrp1 <-n
nGrp2 <- m 
#n和m是对照组和处理组的组数。    
design <- cbind(sampleGroup1=1, sampleGroup2vs1=c(rep(0, nGrp1), rep(1, nGrp2)))

开始处理

[AppleScript] 纯文本查看 复制代码

library(limma)
library(GSVA) 
#根据需要调整一下表达矩阵
rownames(gene) <- gene[,1]    #添加gene中的行标题
gene <- gene[,-1]         #删除第一列
topMatrixGSVA <- gsva(data.matrix(gene), genesets, min.sz=1, max.sz=999999, kcdf="Poisson", abs.ranking=FALSE, verbose=TRUE)        #得到GSVA在这些通路上的分值。

到这里GSVA要做的事已经做完了.
后面就是根据分组进行差异分析了,或者不进行差异分组,直接用得分做热图.

longgu1314

报错如图

longgu1314

金恩12 发表于 2018-9-3 22:21
感觉在使用软件时,还是应该看一下要使用的软件或者R包的官网的document和发表的相关文献.GSVA和GSEA的 ...

当时正在入门，感谢耐心回复，虽然现在说谢谢已经过去一年了，还是认真地说句谢谢！