###go富集
對(duì)populus_simonii_genome_data.emapper.annotations處理
vi chang_go1.py
#!/usr/bin/python
# -*- coding: utf-8 -*-
import sys, os
x = sys.argv[1]
file = open(x, "r")
lines = file.readlines()
for line in lines:
    line=line.strip()
    if line.startswith("#"):
        continue
    else:
        tmp=line.split("\t")
        if tmp[12] == "-":
            continue
        else:
            if(tmp[12].find(',')):
                tmp2=tmp[12].split(',')
                for i in range(0,len(tmp2)):
                    print(tmp[0]+"\t"+tmp2[i])
            else:
                print(tmp[0]+"\t"+tmp[12])
python chang_go1.py populus_simonii_genome_data.emapper.annotations >go.all
然后用R語(yǔ)言畫圖

library(clusterProfiler)
go_anno <- read.delim("go.all", header=FALSE, stringsAsFactors =FALSE)
names(go_anno) <- c('gene_id','ID')

go_class <- read.delim('go_term.txt', header=FALSE, stringsAsFactors =FALSE)
names(go_class) <- c('ID','Description','Ontology')
go_anno <-merge(go_anno, go_class, by = 'ID', all.x = TRUE)

gene_list <- read.delim("genego.list",header=FALSE,stringsAsFactors = FALSE)
names(gene_list) <- c('gene_id')
gene_select <- gene_list$gene_id

go_rich <- enricher(gene = gene_select,
                    TERM2GENE = go_anno[c('ID','gene_id')],
                    TERM2NAME = go_anno[c('ID','Description')],
                    pvalueCutoff = 0.05,
                    qvalueCutoff = 0.1,
                    pAdjustMethod = 'BH',
                    maxGSSize = 200)
df<-as.data.frame(go_rich)
pdf(file="enrichment.pdf",width=9,height=5)
dotplot(go_rich)
dev.off()