最近學(xué)習(xí)多基因組比對時(shí)，看到一則12年前在Biostars發(fā)布的提問, Programming Challange: Pairwise Alignments To Multiple Alignment, 收獲頗多，這里記錄下。

提問者，一開始闡釋了自己的問題，也就是他有10-12個非常近的物種的染色體序列，將這些物種和一個參考染色體比對后，得到了多個結(jié)果。他希望，在生成多序列聯(lián)配結(jié)果的同時(shí)不影響到原本單獨(dú)的比對結(jié)果。那么，他想的就是，在各個序列中加上一些插入，就可以得到相對基因組的全局聯(lián)配。

同時(shí)，他強(qiáng)調(diào)了，自己不是來找序列相似度！他想的是有沒有已有的腳本可以做這些事情，或者提供一些代碼上的建議幫助他完成。

甚至，他還給了一個案例，來說明自己的需求

也就是把下面這段

Ref1: CGACAAT--GCACGACAGAGGAAGCAGAACAGATATTTAGATTGCCTCTCATTTTCTCTCCC
Seq1: CGACAATAAGCACGACAGAGGAAGCAGAACAGATA-----ATTGCCTCTCATTTTC-CTCCC

Ref1: CGACAATGCACGACAGAGGAAGC--AGAACAGATATTTAGATTGCCTCTCATTTTCTCTCCC
Seq2: CGACAAT-CACGACAGAGGAAGCTTAGAACAGATATTTAG---GCCTCTCATTTTCTCTCCC

Ref1: CGACAATGCACGACAGAGGAAG----CAGAACAGATATTTAGATTGCCTCTCA----TTTTCTCTCCC
Seq3: CGACAATGCACGACAGAGGAAGTTTTCAGAACAGATATTTAGATTGCCTCTCAAAAATTTTCTCTCCC

變成下面這段。

Ref1: CGACAAT--GCACGACAGAGGAAG----C--AGAACAGATATTTAGATTGCCTCTCA----TTTTCTCTCCC
Seq1: CGACAATAAGCACGACAGAGGAAG----C--AGAACAGATA-----ATTGCCTCTCA----TTTTC-CTCCC
Seq2: CGACAAT---CACGACAGAGGAAG----CTTAGAACAGATATTTAG---GCCTCTCA----TTTTCTCTCCC
Seq3: CGACAAT--GCACGACAGAGGAAGTTTTC--AGAACAGATATTTAGATTGCCTCTCAAAAATTTTCTCTCCC

我覺得大部分人看到這樣子的提問，就都知道提問者到底需要什么，也就不需要花太多時(shí)間思考題目，問提問者更多細(xì)節(jié)。

排名第一的回答來自于唐海寶老師

他首先回答了，作者需要的工具是TBA/MULTIZ, 可以從Miller Lab下載.

接著補(bǔ)充了一點(diǎn)細(xì)節(jié)，多序列聯(lián)配的潛在原則是，gap和insertion的引入和你比對序列的順序有關(guān)。也就是說，在很多情況下，seq1-seq2-seq3和seq2-seq1-seq3`結(jié)果是不一樣的，“once a gap, always gap”

最后說了TBA軟件的不足，即需要用他們定義的MAF格式作為輸入，也就是用戶得做一些格式轉(zhuǎn)化你工作。并給了一個使用案例

已知參考序列是ref1, 用于比對的序列是seq1, seq2, seq3。比對之后得到ref1.seq1.sing.maf, ref1.seq2.sing.maf, ref1.seq3.sing.maf, 這三個文件。提供一個進(jìn)化樹描述序列的順序，如(((ref1 seq1) seq2) seq3， 表示ref1和seq1近，后面跟seq2近，最后是seq3。

最后運(yùn)行如下命令

tba "(((ref1 seq1) seq2) seq3)" *.*.maf tba.maf

輸出的tba.maf 就是你想要的結(jié)果，Good luck!