前幾天，實驗室的師弟師妹通過本地blast獲取一些沒有基因組注釋物種的蛋白編碼序列。原本以為可以快速地進(jìn)行下一步的選擇壓力分析，沒想到卻在多序列比對這一環(huán)節(jié)出現(xiàn)了棘手的問題。以前，我都是經(jīng)過PRANK軟件進(jìn)行多序列比對，然后再使用Gblocks軟件對數(shù)據(jù)進(jìn)行過濾的?，F(xiàn)在，由于師弟師妹在拼接CDS序列時，有些堿基并不是保留3的倍數(shù)，造成很多編碼序列出現(xiàn)移碼突變，甚至變成了偽假基因（幾百個基因）。雖然，現(xiàn)在可以進(jìn)行多序列比對的軟件有很多，比如Muscle、MAFFT、PRANK等，但它們均不能解決移碼突變的問題。于是，我開始去google搜尋解決方案。功夫不負(fù)有心人，我最終找到了完美的解決辦法，收獲頗豐，所以把這個過程記錄下來。

多序列比對表示不同序列中的氨基酸/核苷酸的位點(diǎn)同源性。將序列比對用于進(jìn)化分析時，處于相同位置的氨基酸/核苷酸位點(diǎn)則被認(rèn)為在進(jìn)化上是同源的，并且具有共同的祖先。通過研究序列比對中同源序列是如何隨時間變化的，可以推斷序列的結(jié)構(gòu)和功能是如何進(jìn)化的。比如，對多個物種同源序列的比較可以發(fā)現(xiàn)發(fā)生替換的位點(diǎn)，進(jìn)而判斷出哪個序列在哪個位置發(fā)生變化。另外，兩個序列之間的長度差異可以通過刪除一個序列中的現(xiàn)有堿基位點(diǎn)或插入另一個序列中的新堿基位點(diǎn)來解釋。兩個序列之間發(fā)生替換差異的進(jìn)化譜系，可以用外群序列來解決，也可以通過推斷兩個序列的祖先的特征狀態(tài)來解決。因此，以上所有情況的結(jié)論準(zhǔn)確性很可能取決于多序列比對（MSA）的準(zhǔn)確性。

最近的研究表明，MSA算法在分析基因組序列時會產(chǎn)生不同的結(jié)果，包括系統(tǒng)發(fā)生樹推斷和適應(yīng)性進(jìn)化的檢測。換句話說，錯誤的MSA將會產(chǎn)生一個非真實的進(jìn)化歷史信號，從而導(dǎo)致錯誤的推斷。為了減少M(fèi)SA錯誤的影響，許多研究人員在優(yōu)化MSA算法和MSA質(zhì)量過濾軟件等方面做出了大量的努力。然而，MSA不僅包含比對錯誤，而且還包含序列自身錯誤（如測序錯誤、組裝錯誤或錯誤的基因注釋等）。所以，正確選擇MSA分析流程可以顯著減少這些錯誤的產(chǎn)生，從而獲得高質(zhì)量的保守同源序列用于下游的分子進(jìn)化分析。

接下來，我主要是介紹一些平時比較常用的主流分析軟件以及它們的優(yōu)缺點(diǎn)。可能還有很多優(yōu)秀的軟件，由于個人精力有限，暫時只講一下我自己比較熟悉的工具，請大家見諒。

01. Multiple Sequence Alignment

2013年，Blackburne 和 Whelan 在Molecular Biology and Evolution期刊發(fā)表他們的研究成果，該成果證明了MSA方法在下游分析中的關(guān)鍵作用，突出了在分析中選擇的不同的MSA軟件對其結(jié)果有明顯的影響。

是最早開始使用的多序列比對工具，由 Feng 和 Doolittle 于1987年提出，該程序有許多版本，可以基于多種平臺，目前引用次數(shù)已近100000次（Thompson et al. 1994; Larkin et al. 2007）。它采用一種漸進(jìn)的比對方法，先將多個序列兩兩比對構(gòu)建距離矩陣，反應(yīng)序列之間兩兩關(guān)系；然后根據(jù)距離矩陣計算產(chǎn)生系統(tǒng)進(jìn)化指導(dǎo)樹，對關(guān)系密切的序列進(jìn)行加權(quán)，然后從最緊密的兩條序列開始，逐步引入臨近的序列并不斷重新構(gòu)建比對，直到所有序列都被加入為止。因此，它是一種試探算法，所以漸進(jìn)比對不能保證能夠得到最優(yōu)的比對。

是一款速度最快的比對軟件之一，在速度和精度上都優(yōu)于ClustalW，目前引用次數(shù)高達(dá)26246次（Edgar 2004）。Muscle采用的是迭代方法進(jìn)行比對運(yùn)算，每一次最優(yōu)化過程就是迭代過程，通過不斷地使用動態(tài)規(guī)劃法重排來糾正這種錯誤，同時對這些亞類群進(jìn)行比對以獲得所有序列的全局比對。一個形象的例子是，10000條長度為350bp的數(shù)據(jù)進(jìn)行比對只需要十幾分鐘，而用ClustalW則可能需要1年。但是，速度快的后果就是準(zhǔn)確度降低。

的第一個版本于2002年發(fā)布，它使用了一種基于漸進(jìn)對齊的算法，利用快速傅里葉變換對序列進(jìn)行聚類（Katoh K et al. 2002）。其后版本的MAFFT增加了其他算法和操作模式，包括更快地對大量序列進(jìn)行比對的選項、更高精度的比對、非編碼RNA序列的比對以及在現(xiàn)有比對中添加新的序列，目前其引用次數(shù)也近達(dá)20000次。它的比對精度要高于Muscle，速度也較快，但相比其他依賴物種進(jìn)化關(guān)系的比對軟件，它的準(zhǔn)確度還有一些差距。

是由 L?ytynoja 于2005年開發(fā)的軟件，目前它的引用數(shù)達(dá)到1662次。它是一種針對DNA、密碼子（codon）和氨基酸序列的概率多重比對程序，基于一種新的算法，能夠確定各序列在與其共同的祖先分化后發(fā)生的堿基替換，插入（insertion）或缺失（deletion）事件。此外，PRANK可以重新構(gòu)建其祖先序列，有DNA翻譯以及回譯選項（ DNA translation/back-translation）。最后，由于PRANK比對準(zhǔn)確度高，相對比較耗時，不太適合基因組數(shù)據(jù)分析。

是一種貝葉斯后驗比對軟件，它利用馬爾可夫鏈蒙特卡羅來探索給定分子序列數(shù)據(jù)的比對和系統(tǒng)發(fā)育的聯(lián)合空間，同時估計消除了對不準(zhǔn)確的比對引導(dǎo)樹的偏差，在比對過程中采用了更復(fù)雜的替換模型，并且自動利用共享插入/刪除中的信息來幫助推斷系統(tǒng)發(fā)育關(guān)系（Suchard et al. 2006）。其準(zhǔn)確度，與PRANK幾乎差不多（我自己并沒有使用過這款軟件）。

Blackburne 的研究結(jié)果顯示，基于進(jìn)化方法的PRANK與BAli-Phy軟件，可以根據(jù)序列的差異性靈活地使用不同的比對標(biāo)準(zhǔn)（評分矩陣與罰分），然而ClustalW、Muscle以及MAFFT等經(jīng)典比對軟件則不能。他們的結(jié)果表明，在構(gòu)建序列比對時把物種的進(jìn)化關(guān)系及距離考慮進(jìn)去能有效提高比對準(zhǔn)確性，因此比起僅應(yīng)用新的計算方法，可能會產(chǎn)生更大的改進(jìn)潛力。

除此以外，還有許多研究也同樣揭示ClustalW、Muscle以及MAFFT等經(jīng)典比對軟件在MSA精確度方面還存在許多缺陷（如下圖所示）。

L?ytynoja A. 2014

Fletcher W & Yang Z. 2010

是第一個可以用于自動調(diào)整含有移碼變異以及假基因的蛋白編碼序列，而不破壞潛在密碼子結(jié)構(gòu)的多序列比對工具（Ranwez et al. 2011）。它的關(guān)鍵特征是在核苷酸水平上對DNA序列進(jìn)行比對，但有可能包括不是三個堿基的倍數(shù)的間隙長度，即產(chǎn)生移碼，同時基于其氨基酸翻譯對產(chǎn)生的核苷酸比對進(jìn)行評分。這使得可以產(chǎn)生保留潛在密碼子結(jié)構(gòu)的核苷酸比對，同時受益于氨基酸序列的更高相似性（Ranwez et al. 2018）?，F(xiàn)在，這個軟件已被許多基因組分析所廣泛使用，其準(zhǔn)確度也相對較好，然而目前并沒有研究去仔細(xì)分析過它的性能。

OrthoMAM v8/v10數(shù)據(jù)庫就是采用這個軟件的分析流程：CDS and exon sequences are aligned at the codon level in two steps. First, the translated amino acids are aligned using MAFFT (Katoh et al. 2005) and gaps are reported onto the nucleotide sequences. This alignment is then refined using MACSE (Ranwez et al. 2011) to obtain a final codon alignment unaffected by frameshifts, misassemblies, and sequencing errors. Amino acid alignments are then filtered to reduce the impact of errors on evolutionary inferences using HMMcleaner (Philippe et al. 2017).

Ranwez 2018

此外，還有幾種基于隱馬爾科夫模型（Profile HMM Methods）的多序列比對工具，這些工具目前我還不清楚它們的具體優(yōu)缺點(diǎn)。

• SEPP（SATé-enabled Phylogenetic Placement）：解決將short reads放入?yún)⒖夹蛄泻蜆涞南到y(tǒng)發(fā)育問題；
• TIPP（Taxonomic Identification and Phylogenetic Profiling）：解決元組數(shù)據(jù)的分類識別和豐度分析問題
• UPP（Ultra-large alignments using Phylogeny-aware Profiles）：解決非常大的數(shù)據(jù)集對齊的問題，這些數(shù)據(jù)集可能包含一些零碎的數(shù)據(jù)，可以將數(shù)據(jù)集多達(dá)1,000,000條序列對齊；
• HIPPI（Highly Accurate Protein Family Classification with Ensembles of HMMs）：解決蛋白質(zhì)家族分類的問題；

02. Trimming

當(dāng)獲得精準(zhǔn)的多序列比對后，接下來的任務(wù)就是要過濾掉一些低質(zhì)量以及高變異度的序列區(qū)域，僅保留進(jìn)化保守的區(qū)域用于后續(xù)分析。其中，產(chǎn)生這些低質(zhì)量區(qū)域的原因主要有兩個方面：

• 生物學(xué)因素：如果比較來自遠(yuǎn)緣物種的蛋白序列，很有可能只有蛋白質(zhì)的功能部分在序列上具有較高的保守性。其他部分，如第一個外顯子或最后一個外顯子區(qū)域，更有可能改變它們的氨基酸序列，無論是在氨基酸含量方面，還是在insertion/deletion存在的情況下。即使當(dāng)多個序列對齊的預(yù)測是正確的，它也可能對系統(tǒng)發(fā)育樹的推斷產(chǎn)生負(fù)面影響；

• 錯誤：測序錯誤，基因組組裝錯誤、基因預(yù)測錯誤以及多序列比對錯誤；

目前主流的多序列比對過濾軟件可分為以下幾個方面：一是，block-filtering：trimAl、Gblocks；二是，segment filtering：HmmCleaner、PREQUA；三是，sliding window analysis：FasParser2。

是較早時間（2000年）就應(yīng)用于MSA過濾的軟件，由Castresana實驗室所開發(fā)，被引用次數(shù)為6194次。它可以將MSA中大片斷非保守性或者是非同源片段給刪除（6-10bp左右的非同源片段則不能很好識別出來），另外它還對Block（即一段連續(xù)的且中間不含Gap的列）的長度進(jìn)行了限制。有研究把這種方法應(yīng)用于構(gòu)建系統(tǒng)發(fā)育樹，發(fā)現(xiàn)通過刪除原始比對中得分比較低的片段以后，Gblocks有效地提高了系統(tǒng)發(fā)育關(guān)系的確定性（Talavera et al. 2007）。但是，Gblocks的問題在于它太武斷地規(guī)定了某個具體的閾值來判斷比對片段的保留或刪除。比較科學(xué)的閾值選擇應(yīng)該要根據(jù)序列的具體情況來確定，而不是對所有基因的比對都采用同一個閾值，因為不同的基因或同一個基因的不同片段，其進(jìn)化速率可能都不相同。

Capella-Gutiérrez等人于2009年開發(fā)的另一款MSA自動調(diào)整過濾軟件，特別適用于大規(guī)模的系統(tǒng)發(fā)育分析，目前引用次數(shù)已有2157次。它的主要優(yōu)點(diǎn)就是速度快，準(zhǔn)確度高，相比于Gblocks，它可以自動選擇在每個特定比對中使用的參數(shù)，主要包括Gap的比例以及氨基酸相似性的水平，從而優(yōu)化 signal-to-noise ratio 。

該軟件是由中國科學(xué)院昆明動物研究所孫艷波副研究員所開發(fā) link地址，并于2018年在bioinformatics期刊發(fā)表。FasParser2主要是可以在友好界面窗口下對一些常規(guī)序列的操作，尤其是對多基因串聯(lián)合并、序列提取、低質(zhì)量序列鑒定、快速尋找趨同氨基酸變異位點(diǎn)等進(jìn)行批處理化。尤其是，該軟件增加了有效過濾序列比對質(zhì)量的模塊（滑窗分析），可實現(xiàn)序列中對非同源序列的清除，適用于大規(guī)模檢測序列中經(jīng)歷達(dá)爾文正選擇的位點(diǎn)，以及重建系統(tǒng)發(fā)育等分析。軟件作者通過序列模擬測試發(fā)現(xiàn)，其修剪效率非常高，它能夠刪除95%以上的錯誤列，且基本不依賴于比對軟件的選擇，明顯優(yōu)于其他同類程序，如Guidance和Gblocks。另外，經(jīng)過AlignTrim處理后的比對結(jié)果，可顯著降低正選擇檢測中的檢陽性事件（低于5%），且不會導(dǎo)致正選擇的檢測率降低，從而提高了檢測正選擇基因的準(zhǔn)確性。該作者，將其與Gblocks、trimAl軟件進(jìn)行了比較，結(jié)果發(fā)現(xiàn)FasParser2的處理效果要優(yōu)于這兩款軟件（該成果還未正式發(fā)表）。

Yanbo Sun

FasParser

與這兩款軟件也是近兩年新發(fā)表的軟件，它們采取了新的過濾算法（主要基于隱馬爾可夫模型），能夠準(zhǔn)確識別多序列比對中的非同源片段，排除了由于測序質(zhì)量、基因注釋以及可變剪切造成的序列錯誤，其準(zhǔn)確度以及速度相對較快，適用于大規(guī)模數(shù)據(jù)分析（Di Franco A. 2019；Whelan S. 2018）。值得注意的是，PREQUAL針對于非比對的蛋白編碼序列以及氨基酸序列，而HmmCleaner只能操作氨基酸序列（可以是比對后的）。另外，它們對輸入文件要求較高，一般fasta格式，序列中只能存在Gap和序列，不允許存在“？”、“！”、“X”等非編碼字符，還要注意換行符問題。

03. Know Your Limits

• 數(shù)據(jù)類型： DNA vs. RNA, coding vs. non-coding nucleotides (wobble bp), AAs, proteins, etc.
• 數(shù)據(jù)特性：substitution (≠ mutation) rate strength (↑ vs. ↓), indel size and rate (% gap & gap length), pairwise sequence identity (PID), etc.
• 數(shù)據(jù)矩陣特性：# of tips, # of sequences, (alignment length ∝) data matrix weight, e.g., light (K, M) vs. heavy (G, T), etc.
• 計算資源：CPU time and RAM memory

04. Summary

每個人的數(shù)據(jù)可能與其他人的數(shù)據(jù)有很大差別，一定要學(xué)會正確選擇適合自己數(shù)據(jù)的軟件類型，這樣得到的結(jié)果才是準(zhǔn)確、可靠、可重復(fù)的。

生物信息學(xué)與計算機(jī)類似，更新速度很快，幾年前還在使用的經(jīng)典軟件，現(xiàn)在可能已經(jīng)過時了，因而保持與時俱進(jìn)對生物信息人員是十分重要的。

最后，感謝中國科學(xué)院昆明動物所孫艷波老師給與的指導(dǎo)與幫助。

05. 這里放一張?zhí)幚硇蛄兄写嬖谝拼a突變的序列比對流程，已親測（歡迎交流）

圖片發(fā)自簡書App

問題解決了，收獲滿滿，開心?。。?/div>

參考文獻(xiàn)：

1. Blackburne B P, Whelan S. Class of multiple sequence alignment algorithm affects genomic analysis[J]. Molecular biology and evolution, 2012, 30(3): 642-653.
2. L?ytynoja A. Phylogeny-aware alignment with PRANK[M]//Multiple sequence alignment methods. Humana Press, Totowa, NJ, 2014: 155-170.
3. L?ytynoja A, Goldman N. Phylogeny-aware gap placement prevents errors in sequence alignment and evolutionary analysis[J]. Science, 2008, 320(5883): 1632-1635.
4. Fletcher W, Yang Z. The effect of insertions, deletions, and alignment errors on the branch-site test of positive selection[J]. Molecular biology and evolution, 2010, 27(10): 2257-2267.
5. Ranwez V, Douzery E J P, Cambon C, et al. MACSE v2: toolkit for the alignment of coding sequences accounting for frameshifts and stop codons[J]. Molecular biology and evolution, 2018, 35(10): 2582-2584.
6. Sun Y B. FasParser2: a graphical platform for batch manipulation of tremendous amount of sequence data[J]. Bioinformatics, 2018, 34(14): 2493-2495.
7. Whelan S, Irisarri I, Burki F. PREQUAL: detecting non-homologous characters in sets of unaligned homologous sequences[J]. Bioinformatics, 2018, 34(22): 3929-3930.
8. Di Franco A, Poujol R, Baurain D, et al. Evaluating the usefulness of alignment filtering methods to reduce the impact of errors on evolutionary inferences[J]. BMC evolutionary biology, 2019, 19(1): 21.