使用PfamScan的API对蛋白结构域进行注释

@ZyuanZhang

已于 2023-12-04 10:20:26 修改

阅读量1.4k

点赞数 1

分类专栏：生信文章标签： python

于 2022-07-26 11:24:01 首次发布

本文链接：https://blog.csdn.net/weixin_44065416/article/details/125990336

版权

生信专栏收录该内容

10 篇文章 2 订阅

订阅专栏

0. PfamScan简介：

PfamScan是根据Pfam HMM对蛋白质序列进行蛋白家族及结构域的注释的一个工具。
网页端：https://www.ebi.ac.uk/Tools/pfa/pfamscan/
本地版：ftp://ftp.ebi.ac.uk/pub/databases/Pfam/Tools/

注：因为我只是临时用一下，希望能较快地出结果，所以没有去安装本地版，而是调用了PfamScan提供的REST API，所以本文记录的是调用API来运行PfamScan的过程。

1. PfamScan的API下载：

URL: https://www.ebi.ac.uk/seqdb/confluence/display/JDSAT/PfamScan+Help+and+Documentation#PfamScanHelpandDocumentation-WebServices
在这里插入图片描述
下载pfamscan.py即可，但是运行前要先安装xmltramp2模块(pip install xmltramp2)

2. 运行PfamScan：

运行命令：

python pfamscan.py --email "xxxxx@gmail.com" --database pfam-a --sequence targetSeqs.fasta --evalue 50 --format txt --outfile targetResult

参数含义可通过python pfamscan.py --help进行查看。
此处的--email --database --sequence是必须的参数(--email可以用谷歌邮箱，其他邮箱没尝试过)。

需要注意的是：输入文件(targetSeq.fasta)中序列数目不能超过100条。

运行结果：

运行结束后会生成三个文件：
targetResult.out.txt
targetResult.sequence.txt
targetResult.submission.params

其中targetResult.out.txt就是所需的输出文件(格式如下图所示)
在这里插入图片描述
列名注释：
Seq id : 蛋白/基因的编号
Alignment start : 基因/蛋白序列比对的结构域起始位置
Alignment end : 基因/蛋白序列比对的结构域终止位置
Envelope start : HMM 模型预测的基因/蛋白序列的结构域起始位置
Envelope end : HMM 模型预测的基因/蛋白序列的结构域终止位置
Hmm acc: 基因/蛋白序列对应结构的模型在Pfam中的编号
Hmm name : 基因/蛋白序列对应结构的模型在Pfam中的名称
Type: 基因/蛋白序列匹配到 Pfam 数据库中对应结构的分类水平，蛋白家族或者结构域
Hmm start: 比对上的部分在数据库匹配序列上的起始位置
Hmm end : 比对上的部分在数据库匹配序列上的终止位置
Hmm length：比对上的长度
Bit score：根据比对和 HMM 模型得出的基因/蛋白序列结构的评分，打分越高，可信度越高
E-value：比对的 E值（E值越小，可信度越高）
Significance: 基因/蛋白序列在数据库中匹配结构的数目
Clan : Pfam 数据库中按照蛋白质序列，结构以及 HMM 文件而分成的类群

3. 批量提交蛋白序列

import os
from Bio import SeqIO
import time


def runPfamScan(inFasta, outFile, email="xxx@gmail.com"):
    """
        inFasta中序列数目 <= 100
    """
    cmd = 'python pfamscan.py --email '+email+' --database pfam-a --sequence '+inFasta+' --evalue 50 --format txt --outfile '+outFile
    os.system(cmd)


def splitFasta(infile, threshold=100):
    """
       将fasta文件中的序列拆分成多个fasta文件，每个fasta文件中的序列数目<=100
       threshold: 拆分完的fasta文件中蛋白序列的数目的阈值，默认为100
    """
    id_list, seq_list = [], []
    for record in SeqIO.parse(infile, "fasta"):
        id_list.append(str(record.id))
        seq_list.append(str(record.seq))

    for i in range(0, len(id_list), threshold):
        if i+threshold <= len(id_list):
            outF = open("/".join(infile.split("/")[0:-1])+"/part_"+infile.split("/")[-1].split(".")[0]+"_"+str(i+1)+"_"+str(i+threshold)+".fasta", "w")
            for j in range(i, i+threshold):
                outF.write(">%s\n%s\n" % (id_list[j], seq_list[j]))
            outF.close()
        else:
            outF = open("/".join(infile.split("/")[0:-1])+"/part_"+infile.split("/")[-1].split(".")[0]+"_"+str(i+1)+"_"+str(len(id_list))+".fasta", "w")
            for k in range(i, len(id_list)):
                outF.write(">%s\n%s\n" % (id_list[k], seq_list[k]))
            outF.close()

def main():
	print("Start: ", time.ctime(), flush=True)
    SplitFasta(infile) ## 输入的fasta文件
    f_list = ["/".join(infile.split("/")[0:-1])+f for f in os.listdir("/".join(infile.split("/")[0:-1])) if f.startswith("part_virus")] ## SplitFasta()输出的结果
    for f in f_list:
        runPfamScan(inFasta = vf,
                    outFile = "./results/"+f.split("/")[-1].split(".")[0]) ## 结果保存在result目录下 (可以手动指定其他的目录)
    print("End: ", time.ctime(), flush=True)

if __name__ == "__main__":
	main()