前些日子，由于各种原因，没想写公众号，其实也不是没有时间，就是一旦一件事情有一天没去做，就觉得其实也可以明天做，明天又觉得做不做都行，再过几天，就可能不去做了。所以，有些时候还是要定个计划和目标的，不管有没有意义，不管做多做少，就像吃饭一样，每天吃多吃少还是要吃点的。
今天我们把转录组的差异基因通过Go富集和kegg分析，看看他们是不是真的是差异表达基因。
1.Go富集分析
由于稻瘟病菌没有完整的Go号和orgdb库，需要自己建一个。
第一步首先需要稻瘟病菌MGG基因号和对应的Go号，这个可以在http://geneontology.org/中获取，但不是很全，还有可以通过序列比对找到其他物种的Go号，总之不管用什么方法，找的号越全，差异分析越全。这里通过以前专业分析的数据爬取了4万+个Go号，应该是很全的了。
比如从这个excel中提取MGG对应的Go号
代码如下
import ref =open('CC.csv','w')for line in open('Expression.csv','r'): BF = line.split(',')[1] CC = line.split(',')[2] MF = line.split(',')[3] BF = CC if BF != '': # print(BF) if '///' in BF: MGG = line.split(',')[0] l = BF.split('///') for s in l: go = re.search(r'd{7}',s.strip()) if go: go = go.group() else: go = 'Null' ev = re.search(r'[[](.+)[]]',s) if ev: ev = ev.group(1) else: ev = 'Null' mgg = re.search(r'MGG_d{5}',s) if mgg: mgg = mgg.group() else: mgg = 'Null' #print(BF) f.write(MGG + ',' +'GO:'+str(go)+','+ 'IEA' +','+'CC'+ 'n') #print(MGG + ',' +'GO:'+str(go)+','+ 'IEA' +','+'MF'+ 'n') else: if '//' in BF: MGG = line.split(',')[0] go = re.search(r'd{7}',BF.strip()) if go: go = go.group() else: go = 'Null' ev = re.search(r'[[](IEA)[]]',BF) if ev: ev = ev.group(1) else: ev = 'Null' mgg = re.search(r'MGG_d{5}',BF) if mgg: mgg = mgg.group() else: mgg = 'Null' # print(BF) f.write(MGG + ',' +'GO:'+str(go)+','+ 'IEA' +','+'CC'+ 'n') # print(MGG + ',' +'GO:'+str(go)+','+ 'IEA' +','+'MF'+ 'n')f.close()
输出是这样的，
一共4万多列。
第二步，有了这个对应的关系，我们就可以用R的Gostats包来进行富集了。
getwd()#获取当前路径library(GOstats)#载入包ago '/Users/zhuxueming/Desktop/R/Mo_goev_forStats.csv',row.names = colnames(ago) 'gene_id',goframeData #去除空数据goFrame #转换成标准模式goAllFrame library('GSEABase')gsc universe as.vector(ago$gene_id)DEGfile '/Users/zhuxueming/Desktop/R/result.csv',header = T,row.names = genes = as.vector(DEGfile$Gene_id)#获取差异基因号#genesparams 'custom', geneSetCollection = gsc, geneIds = genes, universeGeneIds = universe, ontology = 'MF',#选择类型BP，MF，CC pvalueCutoff = 0.05,#p小于0.05 conditional = FALSE, testDirection = 'over')Over write.csv(summary(Over),'/Users/zhuxueming/Desktop/R/GO_MF_hyper.csv') #保存数据
输出结果就是富集的基因类别。
类似的可以把BP，MF，CC都输出
这样，我们就可以根据term和count作图了。
2.kegg分析
和go分析一样，首先需要获取MGG号对应的kegg号，我们可以通过kegg网站爬取相应对应的号码。
import ref = open('all_kegg_path.csv','w')f1 = open('mgr00001.keg','r')dict1 = {}dict2 = {}for i in f1.readlines(): if(i[0] == 'C'): d = re.search('d+',i) if(i[0] == 'D'): h = re.search('Kd+',i) m = re.search('MGG_d+',i) #print(m.group(),h.group(),d.group()) f.write(m.group()+','+h.group()+','+ str(d.group())+'n')f1.close()f.close()
输出结果
有了这个数据，就可以用R包做kegg分析了
library(GOstats)ago '/Users/zhuxueming/Desktop/R/all_kegg_path.csv',header = F,row.names = colnames(ago) 'gene_id',ago$kegg_id "%05d",ago$kegg_id)keggframeData keggframeData keggFrame ' Pyricularia oryzae 70-15')library('GSEABase')gsc universe as.vector(ago$gene_id)DEGfile '/Users/zhuxueming/Desktop/km-k/result.csv',header = T,row.names = genes as.vector(DEGfile$gene_id)params 'custom', geneSetCollection = gsc, geneIds = genes, universeGeneIds = universe, pvalueCutoff = 0.05, testDirection = 'over')Over write.csv(summary(Over),'/Users/zhuxueming/Desktop/km-k/KEGG_hyper.csv')
输出结果
