原文地址:[protocol]go enrichment analysis作者:生龙活虎背景: 什么是富集分析,自己可以百度.我到目前也没发现一个比较通俗易懂的介绍.直接理解为一种统计学方法就可以了. 用于查看显著性. 富集分析有很多种,最常见的是go富集分析.也有pathway富集分析.[pathway的我目前不会啊 ::_:: ] 工具: 也有很多种,我这里主要是用ontologizer (links:) 需要的文件: 1.gene_ontology_edit.obo (这个类似go的库文件, links:) 2.whole_gene_list (这是所研究物种中,包含的所有gene_list) 通过一下脚本可以从cds文件中提取gene_list grep "^" seq ence.fasta | tr -d ""  seq enceheader.list 3. sample_1_gene_list (其中一个点的gene_list) 4.最后一个也是最关键的.就是要有一个association file .这个通过一下脚本生成 4.1需要2个文件 和一个脚本 文件1:go_ids 文件 格式如下 gene go:, go:,go:,go:,go:,go:, go:,go:,go:, go:00 文件2: goidtermindexfile 这个可以下载到 links: 脚本1: import sys #inp t the gene annotation file (infile) and yo r goidterms index file (gofile)#infile is yo r c stom annotations#gofile is the go ids to go categories index yo download from the go website infile = open(sys.argv[1], r ) gofile = open(sys.argv[2], r ) #parse the goidterms index and store it in the dictionary, gohash gohash = {} for line in gofile: #skip lines that start with ! character as they are comment headers if line[0] != "!": data = line.strip.split( t ) #skip obsolete terms if data[-1] != obs : for info in data: if info[0:3] == "go:": #create dictionary of term aspects gohash[data[0 = data[-1] #here are some col mns that the gaf format wants. #since ontologizer doesn t care abo t this, we can j st make it p db = yo rorganism dbref = pmid: dbevi = iso dbobjecttype = gene dbtaxon = taxon: dbdate =  dbassignedby = pfam #potential obselete goids that yo have in yo r annotation potentialobs = [] #if yo specified to not print o t obsolete goids, then print o t the .gaf if len(sys.argv) == 3: print !gaf-version: 2.0 #loop thro gh the go annotation file and generate assocation lines. for line in infile: data = line.strip.split( t ) #if gene has go annotations if len(data)  1: #gid is the gene id, goids is an array of assocated go ids gid = data[0] goids = data[1].split( , ) #second col mn of the .gaf file that ontologizer doesn t care abo t dbid = "db." gid #third col mn is the name of the gene which ontologizer does se dbobjsym = gid #for each go id in the line, print o t an association line for goid in goids: if gohash.has_key(goid): dbaspect = gohash[goid] dbobjname = myorganism dbid o tarray = [db,dbid,dbobjsym, 0 ,goid,dbref,dbevi, 0 ,dbaspect, 0 , 0 ,dbobjecttype,dbtaxon,dbdate,dbassignedby] #only print o t the .gaf file if yo didn t specify to print o t obsolete goids. if len(sys.argv) == 3: print t .join(o tarray) else: potentialobs.append(goid) #if there is a 4th arg ment, print o t the potential obsolete list if len(sys.argv) == 4: print n .join(set(potentialobs)) 执行: python  go.c stomannotationfile goidtermindexfile  myorganism.gaf 以上4个文件都有了 可以执行ontologizer了 参考: java -jar ontologizer.jar -a myorganism.gaf -g gene_ontology_edit.obo -s moso_target_gene.txt -p seq necheader.list -c parent-child-union -m westfall-yo ng-single-step -d 0.05 -r 1000 5.这个软件是可以画图的.但是需要下载软件.下载步骤如下: (as root ) y m list available graphviz* y m install graphviz* dot -tpng view-4ho rsmind ced-parent-child-westfall-yo  -oexample.png 还可以参考  method: gene enrichment sing ontologizer witho t a predefined association file by damian kao most gene enrichment websites o t there only allow yo to find enrichments for pop lar model organisms sing pre-established gene ontology annotations. i ran into this problem early on d ring my phd when confronted with having to generate enrichment data on schmidtea mediterranea. software and inp ts to get c stom enrichments, i sed ontologizer. the inp ts for ontologizer are: gene ontology obo file a flatfile containing information on all the c rrent go terms and their relationship to other terms. yo can get it here. list of all gene ids if yo have a .fasta file of all yo r genes, yo can se this command to extract the gene id information: grep "^" seq ence.fasta | tr -d ""  seq enceheader.list list of gene ids yo want to find enrichment for the s bset of all yo r genes that yo want the enrichment for, ie. differentially expressed gene. ontologizer can take in an entire folder of gene id lists if yo set a directory as inp t. association file for yo r organism this is a go annotation file (.gaf). the specifications is detailed here. fort nately, ontologizer doesn t really care abo t most of the database information col mns in the .gaf file, so yo can fake most of it. to generate this file, yo will have to first generate go annotations for yo r genes thro gh software like interproscan or blast2go. yo will also need a tab delimited index of go ids to go category which yo can get here. generating a gene association file ass ming yo r go annotations are formatted in two col mns where first col mn is the gene id and second col mn is a comma separated list of yo r go ids: gene go:,  go:,go:,go:,go:,go:,go:, go:,go:, go:00 here is the python script to generate the association file with comments: import sys #inp t the gene annotation file (infile) and yo r goidterms index file (gofile)#infile is yo r c stom annotations#gofile is the go ids to go categories index yo download from the go website infile = open(sys.argv[1], r ) gofile = open(sys.argv[2], r ) #parse the goidterms index and store it in the dictionary, gohash gohash = {} for line in gofile: #skip lines that start with ! character as they are comment headers if line[0] != "!": data = line.strip.split( t ) #skip obsolete terms if data[-1] != obs : for info in data: if info[0:3] == "go:": #create dictionary of term aspects gohash[data[0 = data[-1] #here are some col mns that the gaf format wants. #since ontologizer doesn t care abo t this, we can j st make it p db = yo rorganism dbref = pmid: dbevi = iso dbobjecttype = gene dbtaxon = taxon: dbdate =  dbassignedby = pfam #potential obselete goids that yo have in yo r annotation potentialobs = [] #if yo specified to not print o t obsolete goids, then print o t the .gaf if len(sys.argv) == 3: print !gaf-version: 2.0 #loop thro gh the go annotation file and generate assocation lines. for line in infile: data = line.strip.split( t ) #if gene has go annotations if len(data)  1: #gid is the gene id, goids is an array of assocated go ids gid = data[0] goids = data[1].split( , ) #second col mn of the .gaf file that ontologizer doesn t care abo t dbid = "db." gid #third col mn is the name of the gene which ontologizer does se dbobjsym = gid #for each go id in the line, print o t an association line for goid in goids: if gohash.has_key(goid): dbaspect = gohash[goid] dbobjname = myorganism dbid o tarray = [db,dbid,dbobjsym, 0 ,goid,dbref,dbevi, 0 ,dbaspect, 0 , 0 ,dbobjecttype,dbtaxon,dbdate,dbassignedby] #only print o t the .gaf file if yo didn t specify to print o t obsolete goids. if len(sys.argv) == 3: print t .join(o tarray) else: potentialobs.append(goid) #if there is a 4th arg ment, print o t the potential obsolete list if len(sys.argv) == 4: print n .join(set(potentialobs)) p tting it all together make s re yo have yo r c stom go annotation file generated in the format shown above and the goid to term aspect index file from here (text file version of terms, ids, secondary ids, obsoletes ). save this script and se it by: python  go.c stomannotationfile goidtermindexfile  myorganism.gaf if yo want to check if any of yo r annotation is potentially deprecated or obsolete, se this command: python  go.c stomannotationfile goidtermindexfile y  potentialdeprecated.goids now that yo have the gene assocation file, yo can se ontologizer by: java -jar ontologizer.jar -a myorganism.gaf -g gene_ontology.1_2.obo -o o t -p myorganism.ids -s in the parameters sed are: a - the gaf file g - the obo file o - the o tp t directory p - the list of all gene ids s - the inp t directoryfile of gene s bsets the enrichments with all relevent statistics will be generated per inp t gene list in a tab delimited file.