蛋白質(zhì)都是由相似的小型結(jié)構(gòu)域組成的。如果我們有若干個(gè)已知的蛋白序列狈定，那我們就可以根據(jù)這些蛋白序列比較其含有的保守域，尋找在蛋白數(shù)據(jù)庫(kù)中上是否也有一樣保守域的蛋白。而后根據(jù)統(tǒng)計(jì)學(xué)模型僵刮，將顯著性較高的蛋白序列預(yù)測(cè)為同一類基因家族蛋白。
隨著蛋白質(zhì)數(shù)據(jù)庫(kù)的日趨完善鹦牛，使用蛋白質(zhì)結(jié)構(gòu)域進(jìn)行序列比對(duì)相比起傳統(tǒng)的序列全長(zhǎng)比對(duì)更具優(yōu)勢(shì)搞糕。對(duì)于每個(gè)蛋白質(zhì)家族，通常有數(shù)千個(gè)已知的同源蛋白可以比對(duì)成較深的多重序列曼追。序列比對(duì)揭示了一種特定于該結(jié)構(gòu)域的結(jié)構(gòu)和功能的進(jìn)化模式(profile)窍仰。這些模式可以用概率模型捕捉到。HMMER能夠從蛋白或核酸序列中提取出域家族從而構(gòu)建隱式馬爾可夫模型（profile hidden Markov models, profile HMMs）礼殊，從而用于同源序列檢索驹吮，注釋新的序列。
1.軟件安裝
需要用到的軟件包括mafft晶伦，hmmer碟狞，seqkit。
$ conda install -c bioconda mafft hmmer seqkit
2.MAFFT多序列比對(duì)蛋白質(zhì)
MAFFT是一款多序列比對(duì)軟件婚陪，相比起多序列比對(duì)的明星軟件ClustalW族沃，MAFFT在準(zhǔn)確性和速度上均具有優(yōu)勢(shì)。準(zhǔn)確性上MAFFT>Muscle>T-Coffee>ClustalW，比對(duì)速度上Muscle>MAFFT>ClustalW>T-Coffee [1]脆淹。因此我們?cè)谶@里采用MAFFT進(jìn)行多序列比對(duì)常空。
將待比對(duì)的序列手動(dòng)收集并存放于ZAR1_new.fas文件后，我們采用MAFFT的方式進(jìn)行多序列比對(duì)盖溺。注意：mafft可調(diào)整的參數(shù)較多漓糙，可根據(jù)需求選擇適當(dāng)?shù)膮?shù)。
$ mafft --localpair --maxiterate 1000 ZAR1_new.fas > ZAR1_aligned.fas
2.1可調(diào)整的比對(duì)算法：
2.1.1 mafft準(zhǔn)確度優(yōu)先的比對(duì)算法（Accuracy-oriented methods）：
#L-INS-i (最準(zhǔn)確的算法;適用于200條序列以下的比對(duì)):
$ mafft --localpair --maxiterate 1000 input [> output]
#G-INS-i (適用于序列長(zhǎng)度相似的比對(duì);200條序列以下為佳):
$ mafft --globalpair --maxiterate 1000 input [> output]
#E-INS-i (適用于包含大范圍非比對(duì)區(qū)的序列;200條序列以下為佳烘嘱；--ep 0選項(xiàng)代表允許超長(zhǎng)gap的出現(xiàn)):
$ mafft --ep 0 --genafpair --maxiterate 1000 input [> output]
2.1.2 mafft速度優(yōu)先的比對(duì)算法（Speed-oriented methods）:
#FFT-NS-i法:
$ mafft --retree 2 --maxiterate 2 input [> output]
#FFT-NS-i法(最高1000次迭代):
$ mafft --retree 2 --maxiterate 1000 input [> output]
#FFT-NS-2法(快):
$ mafft --retree 2 --maxiterate 0 input [> output]
#NW-NS-2法(快昆禽，且不進(jìn)行FFT近似估計(jì)):
$ mafft --retree 2 --maxiterate 0 --nofft input [> output]
#FFT-NS-1法(更快; 推薦在比對(duì)2000條以上序列時(shí)使用):
$ mafft --retree 1 --maxiterate 0 input [> output]
#NW-NS-PartTree-1 (推薦序列數(shù)~10,000到~50,000條時(shí)使用):
$ mafft --retree 1 --maxiterate 0 --nofft --parttree input [> output]
2.1.3 mafft群體間比對(duì)的算法（Group-to-group alignments）:
$ mafft-profile group1 group2 [> output]
以上MAFFT的參數(shù)命令行均有簡(jiǎn)寫形式，詳情請(qǐng)見Mafft Manual 蝇庭。
3.hmmbuild將多序列比對(duì)文件轉(zhuǎn)化為隱馬模型
我們采用HMMER軟件進(jìn)行隱馬模型的建立为狸。
#將多序列比對(duì)文件轉(zhuǎn)化為隱式馬爾可夫模型
$ hmmbuild ZAR1.hmm ZAR1_aligned.fas
4.利用隱馬模型搜索結(jié)構(gòu)域類似的蛋白質(zhì)
通過(guò)隱馬模型搜索蛋白數(shù)據(jù)庫(kù)中符合該結(jié)構(gòu)的蛋白質(zhì)。將剛產(chǎn)生的profile輪廓文件作為輸入遗契，檢索靶向數(shù)據(jù)庫(kù)中符合該輪廓的蛋白序列辐棒，最終按照符合度輸出序列結(jié)果。Hongyang_pep.fa是先行下載好的蛋白質(zhì)組序列文件牍蜂，以fasta格式呈現(xiàn)漾根。
#在Hongyang_pep.fa蛋白質(zhì)組中搜索具有ZAR1.hmm特征的蛋白質(zhì)
$ hmmsearch ZAR1.hmm Hongyang_pep.fa > hmmer_result.out
#設(shè)定bit-score閾值篩選搜索結(jié)果，此處設(shè)定bit-score閾值為15
$ hmmsearch -T 15 ZAR1.hmm Hongyang_pep.fa > hmmer_bit.out
#設(shè)定bit-score閾值篩選搜索結(jié)果鲫竞。默認(rèn)為 10, 表示每個(gè)搜索報(bào)告大約 10 個(gè)錯(cuò)誤結(jié)果辐怕。
$ hmmsearch -E 0.0001 ZAR1.hmm Hongyang_pep.fa > hmmer_e.out
#此處設(shè)定過(guò)濾閾值-E如果是e-100類似形式會(huì)報(bào)錯(cuò)，因此建議比對(duì)后使用awk進(jìn)行過(guò)濾从绘。
此外寄疏，常用到的HMMER命令還包括:
hmmbuild： 用多重比對(duì)序列構(gòu)建HMM模型；
hmmsearch： 使用HMM模型搜索序列庫(kù)僵井；
hmmscan： 使用序列搜索HMM庫(kù)陕截；
hmmalign： 使用HMM為線索，構(gòu)建多重比對(duì)序列批什；
5.獲取匹配的蛋白質(zhì)并進(jìn)行tblastn檢索新的同源蛋白
利用hmm模型在蛋白質(zhì)組序列中尋找相似的蛋白后农曲，可以通過(guò)seqkit提取該序列（新建grep.txt文件并將待提取序列的名稱保存于此）。再通過(guò)tblastn會(huì)將庫(kù)中的核酸翻譯成蛋白序列驻债，在核酸庫(kù)中尋找與該蛋白相似的核酸序列乳规。實(shí)際使用時(shí)可以采用全基因組建立核酸庫(kù)，即可搜索全基因組內(nèi)可能與目標(biāo)蛋白相似的序列合呐。
#seqkit提取隱馬模型預(yù)測(cè)的序列暮的，保存于Actinidia05846.t1.fa文件
$ seqkit grep -f grep.txt Hongyang_pep.fa -o Actinidia05846.t1.fa
#基因組建立核酸庫(kù)并命名為：canu_genome
$ makeblastdb -in 11251AaHscanu.contigs.fasta -dbtype nucl -parse_seqids -input_type fasta -out canu_genome
#核酸庫(kù)中比對(duì)目標(biāo)序列
$ tblastn -db canu_genome -query Actinidia05846.t1.fa -outfmt 6 -out tblastn_canu.result
6.對(duì)新檢索的蛋白序列進(jìn)行HMM結(jié)構(gòu)域的注釋
對(duì)于剛剛找到的蛋白，如果我們希望探究其功能淌实，往往會(huì)對(duì)其結(jié)構(gòu)域進(jìn)行搜索冻辩。畢竟結(jié)構(gòu)決定了功能猖腕。待探究的輸入文件可以是單個(gè)蛋白序列，多蛋白序列微猖，hmm隱馬模型谈息。我們采用pfam網(wǎng)站對(duì)蛋白序列進(jìn)行注釋缘屹。首先需要下載Pfam注釋庫(kù)文件凛剥，Pfam網(wǎng)站中保留的庫(kù)文件目前只有A數(shù)據(jù)庫(kù)，A數(shù)據(jù)庫(kù)代表著經(jīng)過(guò)手工校正的高質(zhì)量數(shù)據(jù)庫(kù)轻姿。此外犁珠，我們還需要對(duì)庫(kù)文件進(jìn)行初步的二進(jìn)制壓縮和引索處理。
#下載Pfam庫(kù)文件
$ wget https://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz
#解壓Pfam庫(kù)文件
$ gzip -d Pfam-A.hmm.gz
#壓縮Pfam庫(kù)文件:此處的Hmm文件以文本形式保存互亮，壓縮為二進(jìn)制有助于加速運(yùn)算犁享，建立成索引數(shù)據(jù)庫(kù)。
$ hmmpress Pfam-A.hmm
而后我們需要用到HMMER軟件中的hmmscan進(jìn)行Pfam注釋豹休。將待比對(duì)的序列放在hongyang_RPM1_like.fa文件中炊昆，使用剛剛壓縮得到的庫(kù)Pfam-A.hmm作為注釋參考。得到的三個(gè)文件result.txt威根，result.tbl凤巨，result.dom分別表示待比對(duì)序列的注釋結(jié)果文件，注釋出的蛋白域信息文件洛搀，帶有起止位置信息的蛋白結(jié)構(gòu)域信息文件敢茁。
#hmmscan搜索蛋白質(zhì)所含的結(jié)構(gòu)域
$ hmmscan -o result.txt --tblout result.tbl --domtblout result.dom --noali -E 1e-5 Pfam-A.hmm hongyang_RPM1_like.fa
#更多可選參數(shù)：
# -h:顯示幫助信息
# -o FILE:將結(jié)果輸出到指定的文件中。默認(rèn)是輸出到標(biāo)準(zhǔn)輸出留美。
# --tblout FILE:將蛋白質(zhì)家族的結(jié)果以表格形式輸出到指定的文件中彰檬。默認(rèn)不輸出該文件。
# --domtblout FILE:將蛋白結(jié)構(gòu)域的比對(duì)結(jié)果以表格形式輸出到指定的文件中谎砾。默認(rèn)不輸出該文件逢倍。該表格中包含query序列起始結(jié)束位點(diǎn)與目標(biāo)序列起始結(jié)束位點(diǎn)的匹配信息。
# --acc:在輸出結(jié)果中包含 PF 的編號(hào)景图，默認(rèn)是蛋白質(zhì)家族的名稱瓶堕。
# --noali:在輸出結(jié)果中不包含比對(duì)信息。輸出文件的大小則會(huì)更小症歇。
# -E FLOAT:設(shè)定 E_value 閾值郎笆，推薦設(shè)置為 1e-5 。
# -T FLOAT:設(shè)定 Score 閾值忘晤。
# --domE FLOAT:設(shè)定domain比對(duì)的E_value閾值宛蚓。類似-E參數(shù)。
# --cpu：多線程運(yùn)行的CPU设塔。默認(rèn)應(yīng)該是大于1的凄吏，表示支持多線程運(yùn)行。但其實(shí)估計(jì)一般一個(gè)hmmscan程序利用150%個(gè)CPU。并且若進(jìn)行并行化調(diào)用hmmscan痕钢，當(dāng)并行數(shù)高于4的時(shí)候图柏，會(huì)報(bào)錯(cuò)：Fatal exception (source file esl_threads.c, line 129)。這時(shí)任连，設(shè)置--cpu的值為1即可蚤吹。
結(jié)果文件示例如下：
Query: Actinidia05846.t1 [L=955]
Scores for complete sequence (score includes all domains):
--- full sequence --- --- best 1 domain --- -#dom-
E-value score bias E-value score bias exp N Model Description
------- ------ ----- ------- ------ ----- ---- -- -------- -----------
1.1e-60 205.0 0.2 1.7e-60 204.5 0.2 1.3 1 NB-ARC NB-ARC domain
7.6e-24 83.9 0.2 2.9e-23 82.1 0.2 2.1 1 Rx_N Rx N-terminal domain
3.4e-06 26.8 16.8 0.0049 16.7 0.5 4.6 5 LRR_8 Leucine rich repeat
參考資料：
- 多序列比對(duì)算法MAFFT以及HMMER和profile文件的使用 CSDN：https://blog.csdn.net/weixin_45429249/article/details/109021162
- HMMER User’s Guide. http://eddylab.org/software/hmmer/Userguide.pdf
- Mafft Manual. https://mafft.cbrc.jp/alignment/software/manual/manual.html
- HMMSCAN使用pfam數(shù)據(jù)庫(kù)對(duì)多序列文件進(jìn)行結(jié)構(gòu)域注釋。http://www.reibang.com/p/f6db8af1e2cb
