python 学术_Python

前几天实验室一个师兄给我一个质谱结果，让帮忙做下go的功能富集，数据格式大概是这样的：

由于之前做go和kegg时都是跑流程，像这种针对性的go富集还没做过，说到底，还是由于自己手上缺少数据，没有属于自己的项目，很多细节性的问题都没有经历过。但这不妨碍咱一颗求知的心，我们都是在学习中成长。由于没事的时候逛论坛逛的比较频繁，知道数据的第二列是UniPro数据库的accession，然后该怎么办呢？作为生信人，Google是少不了的，看到Google结果，瞬间明了。根据Google的指引我从网上下载了UniProt数据库里的idmapping.tb.gz文件(wget -c -t 10000 ftp://ftp.pir.georgetown.edu/databases/idmapping/idmapping.tb.gz)，大概18G左右，数据结构如下：

一共有22列，依次分别是：UniProtKB accession，UniProtKB ID，EntrezGene，RefSeq，NCBI GI number，PDB，Pfam，GO，PIRSF，IPI，UniRef100，UniRef90，UniRef50，UniParc，PIR-PSD accession，NCBI taxonomy，MIM，UniGene，Ensembl，PubMed ID，EMBL/GenBank/DDBJ，EMBL protein_id；这就有意思了，数据的第八列就是我们想要的go信息。更有意思的是，有了这个数据库信息，我们就可以根据不同数据库的注释信息做go富集啦！

下面要做的是写一个脚本，根据师兄给的结果调出对应的go号，对于会编程的人来说，这点自然不在话下，代码如下：

import sys

USAGE = "\nusage: python %s idmapping.tb.gz blastout outputfile outputfile2\n" % sys.argv[0]

if len(sys.argv) != 5:

print USAGE

sys.exit()

def parseIDmapping(filename):

UniProt_GO = {}

with open(filename, 'r') as f:

for line in f:

lsplit = line.rstrip().split("\t")

if lsplit[7]:

UniProt_GO[lsplit[0]] = lsplit[7]

return UniProt_GO

def parseBlastOut(filename):

tab_res = []

with open(filename, 'r') as f:

for line in f:

lsplit = line.strip('\n').split('\t')

tab_res.append(lsplit[0])

return tab_res

UniProtKB_GO = parseIDmapping(sys.argv[1])

BlastOut = parseBlastOut(sys.argv[2])

OUT = open(sys.argv[3], 'w')

OUT1 = open(sys.argv[4], 'w')

for i in BlastOut:

if i in UniProtKB_GO.keys():

print i

go = UniProtKB_GO[i]

print go

OUT.write(i+"\t"+go+"\n")

else:

OUT1.write(i+"\n")

OUT.close()

OUT1.close()

得到的结果是这样子：

由于使用软件的关系，这种格式貌似还不能达到要求，再写一脚本转换一下：

import re

file1 = open(r"C:\\Users\\wuchangsong\\Desktop\\11.txt")

out_file1 = open(r"C:\\Users\\wuchangsong\\Desktop\\12.txt", "w")

for line1 in file1:

info1 = re.sub('; ','\t',line1)

out_file1.write(info1)

file1.close()

out_file1.close()

结果长这样：

最终结果：

满满的成就感有么有^_^！

版权声明：本文为博主原创文章，未经博主允许不得转载。