转录组分析之--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 (本文总阅读量 次)