生信步驟｜MAFFT結合HMMER進行多序列比對和基于隐馬模型的基因搜尋

蛋白質都是由相似的小型結構域組成的。如果我們有若幹個已知的蛋白序列，那我們就可以根據這些蛋白序列比較其含有的保守域，尋找在蛋白資料庫中上是否也有一樣保守域的蛋白。而後根據統計學模型，将顯著性較高的蛋白序列預測為同一類基因家族蛋白。

随着蛋白質資料庫的日趨完善，使用蛋白質結構域進行序列比對相比起傳統的序列全長比對更具優勢。對于每個蛋白質家族，通常有數千個已知的同源蛋白可以比對成較深的多重序列。序列比對揭示了一種特定于該結構域的結構和功能的進化模式(profile)。這些模式可以用機率模型捕捉到。HMMER能夠從蛋白或核酸序列中提取出域家族進而建構隐式馬爾可夫模型（profile hidden Markov models, profile HMMs），進而用于同源序列檢索，注釋新的序列。

1.軟體安裝

需要用到的軟體包括mafft，hmmer，seqkit。

$ conda install -c bioconda mafft hmmer seqkit
           

2.MAFFT多序列比對蛋白質

MAFFT是一款多序列比對軟體，相比起多序列比對的明星軟體ClustalW，MAFFT在準确性和速度上均具有優勢。準确性上MAFFT>Muscle>T-Coffee>ClustalW，比對速度上Muscle>MAFFT>ClustalW>T-Coffee [1]。是以我們在這裡采用MAFFT進行多序列比對。

将待比對的序列手動收集并存放于

ZAR1_new.fas

檔案後，我們采用MAFFT的方式進行多序列比對。注意：mafft可調整的參數較多，可根據需求選擇适當的參數。

$ mafft --localpair --maxiterate 1000 ZAR1_new.fas > ZAR1_aligned.fas
           

2.1可調整的比對算法：

2.1.1 mafft準确度優先的比對算法（Accuracy-oriented methods）：

#L-INS-i (最準确的算法;适用于200條序列以下的比對):
$ mafft --localpair --maxiterate 1000 input [> output]

#G-INS-i (适用于序列長度相似的比對;200條序列以下為佳):
$ mafft --globalpair --maxiterate 1000 input [> output]

#E-INS-i (适用于包含大範圍非比對區的序列;200條序列以下為佳；--ep 0選項代表允許超長gap的出現):
$ mafft --ep 0 --genafpair --maxiterate 1000 input [> output]
           

2.1.2 mafft速度優先的比對算法（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法(快，且不進行FFT近似估計):
$ mafft --retree 2 --maxiterate 0 --nofft input [> output]

#FFT-NS-1法(更快; 推薦在比對2000條以上序列時使用):
$ mafft --retree 1 --maxiterate 0 input [> output]

#NW-NS-PartTree-1 (推薦序列數~10,000到~50,000條時使用):
$ mafft --retree 1 --maxiterate 0 --nofft --parttree input [> output]
           

2.1.3 mafft群體間比對的算法（Group-to-group alignments）:

以上MAFFT的參數指令行均有簡寫形式，詳情請見Mafft Manual 。

3.hmmbuild将多序列比對檔案轉化為隐馬模型

我們采用HMMER軟體進行隐馬模型的建立。

#将多序列比對檔案轉化為隐式馬爾可夫模型
$ hmmbuild ZAR1.hmm ZAR1_aligned.fas
           

4.利用隐馬模型搜尋結構域類似的蛋白質

通過隐馬模型搜尋蛋白資料庫中符合該結構的蛋白質。将剛産生的profile輪廓檔案作為輸入，檢索靶向資料庫中符合該輪廓的蛋白序列，最終按照符合度輸出序列結果。

Hongyang_pep.fa

是先行下載下傳好的蛋白質組序列檔案，以fasta格式呈現。

#在Hongyang_pep.fa蛋白質組中搜尋具有ZAR1.hmm特征的蛋白質
$ hmmsearch ZAR1.hmm Hongyang_pep.fa > hmmer_result.out

#設定bit-score門檻值篩選搜尋結果，此處設定bit-score門檻值為15
$ hmmsearch -T 15 ZAR1.hmm Hongyang_pep.fa > hmmer_bit.out

#設定bit-score門檻值篩選搜尋結果。預設為 10, 表示每個搜尋報告大約 10 個錯誤結果。
$ hmmsearch -E 0.0001 ZAR1.hmm Hongyang_pep.fa > hmmer_e.out
#此處設定過濾門檻值-E如果是e-100類似形式會報錯，是以建議比對後使用awk進行過濾。
           

此外，常用到的HMMER指令還包括:

hmmbuild： 用多重比對序列建構HMM模型；

hmmsearch： 使用HMM模型搜尋序列庫；

hmmscan： 使用序列搜尋HMM庫；

hmmalign： 使用HMM為線索，建構多重比對序列；

5.擷取比對的蛋白質并進行tblastn檢索新的同源蛋白

利用hmm模型在蛋白質組序列中尋找相似的蛋白後，可以通過seqkit提取該序列（建立grep.txt檔案并将待提取序列的名稱儲存于此）。再通過tblastn會将庫中的核酸翻譯成蛋白序列，在核酸庫中尋找與該蛋白相似的核酸序列。實際使用時可以采用全基因組建立核酸庫，即可搜尋全基因組内可能與目标蛋白相似的序列。

#seqkit提取隐馬模型預測的序列，儲存于Actinidia05846.t1.fa檔案
$ seqkit grep -f grep.txt Hongyang_pep.fa -o Actinidia05846.t1.fa

#基因組建立核酸庫并命名為：canu_genome
$ makeblastdb -in 11251AaHscanu.contigs.fasta -dbtype nucl -parse_seqids -input_type fasta -out canu_genome

#核酸庫中比對目标序列
$ tblastn -db canu_genome -query Actinidia05846.t1.fa -outfmt 6 -out tblastn_canu.result
           

6.對新檢索的蛋白序列進行HMM結構域的注釋

對于剛剛找到的蛋白，如果我們希望探究其功能，往往會對其結構域進行搜尋。畢竟結構決定了功能。待探究的輸入檔案可以是單個蛋白序列，多蛋白序列，hmm隐馬模型。我們采用pfam網站對蛋白序列進行注釋。首先需要下載下傳Pfam注釋庫檔案，Pfam網站中保留的庫檔案目前隻有A資料庫，A資料庫代表着經過手工校正的高品質資料庫。此外，我們還需要對庫檔案進行初步的二進制壓縮和引索處理。

#下載下傳Pfam庫檔案
$ wget https://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz
#解壓Pfam庫檔案
$ gzip -d Pfam-A.hmm.gz
#壓縮Pfam庫檔案:此處的Hmm檔案以文本形式儲存，壓縮為二進制有助于加速運算，建立成索引資料庫。
$ hmmpress Pfam-A.hmm
           

而後我們需要用到HMMER軟體中的hmmscan進行Pfam注釋。将待比對的序列放在

hongyang_RPM1_like.fa

檔案中，使用剛剛壓縮得到的庫

Pfam-A.hmm

作為注釋參考。得到的三個檔案

result.txt

，

result.tbl

，

result.dom

分别表示待比對序列的注釋結果檔案，注釋出的蛋白域資訊檔案，帶有起止位置資訊的蛋白結構域資訊檔案。

#hmmscan搜尋蛋白質所含的結構域
$ hmmscan -o result.txt --tblout result.tbl --domtblout result.dom --noali -E 1e-5 Pfam-A.hmm hongyang_RPM1_like.fa
#更多可選參數：
# -h:顯示幫助資訊
# -o FILE:将結果輸出到指定的檔案中。預設是輸出到标準輸出。
# --tblout FILE:将蛋白質家族的結果以表格形式輸出到指定的檔案中。預設不輸出該檔案。
# --domtblout FILE:将蛋白結構域的比對結果以表格形式輸出到指定的檔案中。預設不輸出該檔案。該表格中包含query序列起始結束位點與目标序列起始結束位點的比對資訊。
# --acc:在輸出結果中包含 PF 的編号，預設是蛋白質家族的名稱。
# --noali:在輸出結果中不包含比對資訊。輸出檔案的大小則會更小。
# -E FLOAT:設定 E_value 門檻值，推薦設定為 1e-5 。
# -T FLOAT:設定 Score 門檻值。
# --domE FLOAT:設定domain比對的E_value門檻值。類似-E參數。
# --cpu：多線程運作的CPU。預設應該是大于1的，表示支援多線程運作。但其實估計一般一個hmmscan程式利用150%個CPU。并且若進行并行化調用hmmscan，當并行數高于4的時候，會報錯：Fatal exception (source file esl_threads.c, line 129)。這時，設定--cpu的值為1即可。
           

結果檔案示例如下：

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

參考資料：

1. 多序列比對算法MAFFT以及HMMER和profile檔案的使用 CSDN：https://blog.csdn.net/weixin_45429249/article/details/109021162
2. HMMER User’s Guide. http://eddylab.org/software/hmmer/Userguide.pdf
3. Mafft Manual. https://mafft.cbrc.jp/alignment/software/manual/manual.html
4. HMMSCAN使用pfam資料庫對多序列檔案進行結構域注釋。https://www.jianshu.com/p/f6db8af1e2cb