转录组分析之--GO和pathway富集分析

再也不是任何意义上的闲来无事了,迫不得已要自己重新分析转录组数据,虽然这套流程在学界早已非常成熟,但对于新手来说每一步并不是想象的那么简单,特别是对于非模式生物的各类数据分析。简单记录一下自认为值得注意的地方。

1.GO富集分析
本来按照实验室的传统,用Blast2GO来做GO分析,不巧的是申请到的Blast2GO一周的试用期刚好与春节假期重叠。白白浪费掉了。不得已去找另一个工具WeGo。这里说明一下,我的输入文件是已经被注释好的Gene-GO号。所以能直接使用Wego,如果是原始的fasta格式的序列文件,Wego是无能为力的。

2.pathway富集分析
说到这个,就有点勾起痛苦的回忆了,Kobas无疑是一款运用非常广泛的KEGG富集分析软件,但不巧的是就刚好在哪几天Kobas的官网打不开!据了解,截止到目前(3/5/2018)这个网站还是不能正常使用。

以下全是废话:~~

好在我在实验室的服务器上找到了Kobas2.0本地版,坑爹的事情再次发生,无论如何软件不能调用,环境变量改了又改,都不行。听师兄说大实验室目录下安装了Kobas3.0,久旱逢甘雨,何况我又是一个对于软件用新不用旧的人。等我兴冲冲调用3.0的时候,坑爹的事情还在继续发生,首先是任然不能正常调用,总是报一些有的没得看不懂的错误,一会儿是找不到module,一会儿是rpy2包出问题什么的。然后咨询了管理员,很不幸,他用自己的账号可以正常调用,师兄也可以正常调用,更不幸的是管理员登我的账号也可以调用!我是一个始终不会轻易认输的人,更别说在这种“羞辱”的情况下,不知道这样的性格是好是坏,但至少它让我浪费了不少时间在这种别人看来鸡毛蒜皮的小事上。直到问题解决的时候,前前后后一天半已经过去了,具体的原因是由于我之前安装了miniconda,调用python的时候是默认在miniconda的文件夹下面的,问题就在这个python的路径上,其实之前我已经想到了是环境变量的问题,所以把通过miniconda安装的所有软件和库都注释了,但是每次调用python发现都是在miniconda的文件下,最终我把miniconda的路径索性也注释了,果然 it works。。

废话不多说了,kobas的用法大致是这样

Usage: run_kobas.py [-l] -i infile [-t intype] -s species [-E evalue] [-R rank] [-N nCPUs] [-C coverage] [-Z ortholog] [-b bgfile] [-d database] [-m method] [-n fdr] [-o outfile] [-c cutoff] [-k kobas_home] [-v blast_home] [-y blastdb] [-q kobasdb] [-p blastp] [-x blastx]

其中需要特别准备的是这个输入文件,我们看到输入文件中ID的支持类型是这样的:

-t INTYPE, --intype=INTYPE
                      input type (fasta:pro, fasta:nuc, blastout:xml,
                      blastout:tab, id:refseqpro, id:uniprot, id:ensembl,
                      id:ncbigene), default fasta:pro

所以就着手准备ID,我手上只有棉花在拟南芥中注释的同源基因的信息(TAIR),并没有ENTREz ID(ncbigene)信息,现在要把TAIR转换成ENTREz,各种艰辛在此不表,最终找到了clusterProfiler这个包,这个包中函数bitr()可以执行ID转换。至于为什么不用clusterProfiler做GO和Pathway富集分析,我想那又是一个非常悲伤的故事,大约也和环境配置有关吧。。

接下来看这个包是怎么转换的,

x <- c("GPX3",  "GLRX",   "LBP",   "CRYAB", "DEFB1", "HCLS1",   "SOD2",   "HSPA2",
       "ORM1",  "IGFBP1", "PTHLH", "GPC3",  "IGFBP3","TOB1",    "MITF",   "NDRG1",
       "NR1H4", "FGFR3",  "PVR",   "IL6",   "PTPRM", "ERBB2",   "NID2",   "LAMB1",
       "COMP",  "PLS3",   "MCAM",  "SPP1",  "LAMC1", "COL4A2",  "COL4A1", "MYOC",
       "ANXA4", "TFPI2",  "CST6",  "SLPI",  "TIMP2", "CPM",     "GGT1",   "NNMT",
       "MAL",   "EEF1A2", "HGD",   "TCN2",  "CDA",   "PCCA",    "CRYM",   "PDXK",
       "STC1",  "WARS",  "HMOX1", "FXYD2", "RBP4",   "SLC6A12", "KDELR3", "ITM2B")
eg = bitr(x, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
head(eg)
##   SYMBOL ENTREZID
## 1   GPX3     2878
## 2   GLRX     2745
## 3    LBP     3929
## 4  CRYAB     1410
## 5  DEFB1     1672
## 6  HCLS1     3059

具体用法参见此处

有了ID,一切都好说了

run_kobas.py  -i infile -t intype -s species -d database -o outfile

Kobas到这里的用法基本结束了,下一章节带来灰常赤鸡的WGCNA的使用说明~

Share (本文总阅读量 次)