Linux下pfam_scan的使用

做基因功能注释都会特别注意基因上有什么功能结构域，通常我们认为，结构域决定了这个基因的功能。随着高通量测序技术的发展，我们完全可以通过一级序列来预测该基因的结构域，pfam数据库是一个不错的选择，本文将记录通过pfam_scan来完成对多序列文件的结构域注释。虽然网页工具也能做这部分内容(https://www.ebi.ac.uk/Tools/hmmer/search/hmmscan)，但是会限制文件大小。所以本地注释方法也需要掌握。

安装

#安装pfam_scan
conda create -n pfam_scan #创建虚拟环境
source activate pfam_scan
conda install pfam_scan

#安装hmmer3
wget http://eddylab.org/software/hmmer/hmmer-3.2.tar.gz
tar -xzvf  hmmer-3.2.1.tar.gz
cd hmmer-3.2
./configure
make
make check
make install 
vim ~/.bashrc#添加环境变量
export PATH=/usr/local/bin:$PATH
source ~/.bashrc#激活环境变量

数据库下载

wget ftp://ftp.ebi.ac.uk:21/pub/databases/Pfam/current_release/Pfam-A.hmm.gz
wget ftp://ftp.ebi.ac.uk:21/pub/databases/Pfam/current_release/Pfam-A.hmm.dat.gz
wget ftp://ftp.ebi.ac.uk:21/pub/databases/Pfam/current_release/active_site.dat.gz
gunzip *.gz

hmmpress Pfam-A.hmm#建库

使用

pfam_scan.pl -fasta /mnt/d/Documents/test.pep.fa -dir /mnt/d/pfam_scan/hmmer-3.2/ -outfile pfam_scan_result.fa -as
#注意：这里-dir后为pfam所在路径非文件

#参数
-dir ：  Pfam_data_file_dir   包含Pfam数据文件的目录[必须] 

-fasta ：  fasta_file   包含序列的输入文件名 [必须]

-e_seq    序列E-value阈值 [不指定则使用默认阈值] 

-e_dom   结构域E-value阈值 [不指定则使用默认阈值]

-b_seq     序列bit score阈值 [不指定则使用默认阈值]

-b_dom    结构域bit score阈值[不指定则使用默认阈值] 

-align       在结果中显示比对片段 [默认关闭] 

-as        预测Pfam-A数据库匹配的active sites[默认关闭] 

-json [pretty]      输出结果使用JSON格式。例如指定值为[pretty]，则输出结果会使用"pretty" JSON格式输出 [默认关闭] 

-cpu     并行工作的CPU数目 [默认全部]

-translate [mode]   将输入序列视为DNA，并在搜索前使用6框翻译的方法进行转换。如果翻译模式[mode]被指定，则必须为"all"或者"orf"。"all"表示完整翻译，包括终止子并且不产生单独的ORFs；"orf"表示只翻译和报告长度大于20的ORFs。
如果使用了翻译参数而没有指定翻译模式，则默认使用"orf"模式。[默认关闭]

结果

pfamscan蛋白结构域部分分析结果说明如下:
(1) seq_id：转录本ID+[0,1,2]，不存在于列表中的转录本为noncoding
(2) hmm start：比对到结构域的起始位置
(3) hmm end：比对到结构域的终止位置
(4) hmm acc：比对到pfam结构域的ID
(5) hmm name：pfam结构域名称
(6) hmm length：pfam结构域的长度
(7) bit score：比对打分分值
(8) E-value：比对的E值，pfam结构域筛选的条件是: Evalue < 0.001

本文参考https://blog.csdn.net/weixin_39886238/article/details/111740890
https://www.jianshu.com/p/fb3bd3de1c38