分类：Python

python .py -i db.fasta -I blastresult.tab -o selected.txt -O filtered.txt
blastresult.tab是根据P值sort后的文件

from __future__ import division
import re
import sys, getopt
import operator
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import generic_nucleotide

opts, args = getopt.getopt(sys.argv[1:], "hI:i:o:O:")
input_info1= ""
input_info2= ""
out_file1 = ""
out_file2 = ""

for op, value in opts:
    if op =="-o": #ARG_pattern_MIN_7030.fasta or ARG_pattern_MAX_7030.fasta
        out_file1 = open(value, "w")
        filename_pro = str (value)
        name_item = filename_pro.strip().split(".")
        filename = name_item[0]
    elif op == "-i": # All_ARGcontig_remove_S_nr.fasta
        input_info1 = open(value, "r")
    elif op == "-I":
        input_info2 = open(value, "r")
    elif op == "-O":
        out_file2 = open(value, "w")
        filename_pro = str (value)
        name_item = filename_pro.strip().split(".")
        filename = name_item[0]      

dict1 = {} 
dict2 = {}
for record in SeqIO.parse(input_info1,"fasta"):
    h = len(record.seq)
   
    dict1[record.id]=h
    dict2[record.id]=record.description
    
        
for line1 in input_info2:
    info1 = line1.strip("\n").split("\t")
    a = float(info1[2])
    b = float(info1[3])
    c = b/float(dict1[info1[1]])
    if a >= 80.0 and c >= 0.8:
        out_file1.write(str(dict1[info1[1]])+'\t'+str(b)+'\t'+dict2[info1[1]]+'\t'+line1)
    else:
        out_file2.write(str(dict1[info1[1]])+'\t'+str(b)+'\t'+dict2[info1[1]]+'\t'+line1)
input_info1.close()
input_info2.close()
out_file1.close()
out_file2.close()

前几天实验室一个师兄给我一个质谱结果，让帮忙做下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()

结果长这样：
最终结果：
满满的成就感有么有^_^！

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