最近學(xué)習(xí)多基因組比對(duì)時(shí)嗓节，看到一則12年前在Biostars發(fā)布的提問(wèn), Programming Challange: Pairwise Alignments To Multiple Alignment, 收獲頗多，這里記錄下。
提問(wèn)者傲茄，一開始闡釋了自己的問(wèn)題巷懈，也就是他有10-12個(gè)非常近的物種的染色體序列赂毯，將這些物種和一個(gè)參考染色體比對(duì)后怀跛，得到了多個(gè)結(jié)果。他希望相叁，在生成多序列聯(lián)配結(jié)果的同時(shí)不影響到原本單獨(dú)的比對(duì)結(jié)果遵绰。那么，他想的就是增淹，在各個(gè)序列中加上一些插入椿访，就可以得到相對(duì)基因組的全局聯(lián)配。
同時(shí)虑润，他強(qiáng)調(diào)了成玫，自己不是來(lái)找序列相似度！他想的是有沒(méi)有已有的腳本可以做這些事情拳喻，或者提供一些代碼上的建議幫助他完成哭当。
甚至，他還給了一個(gè)案例冗澈，來(lái)說(shuō)明自己的需求
也就是把下面這段
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
我覺(jué)得大部分人看到這樣子的提問(wèn)，就都知道提問(wèn)者到底需要什么亚亲，也就不需要花太多時(shí)間思考題目彻采，問(wèn)提問(wèn)者更多細(xì)節(jié)腐缤。
排名第一的回答來(lái)自于唐海寶老師
他首先回答了，作者需要的工具是TBA/MULTIZ, 可以從Miller Lab下載.
接著補(bǔ)充了一點(diǎn)細(xì)節(jié)肛响，多序列聯(lián)配的潛在原則是柴梆，gap和insertion的引入和你比對(duì)序列的順序有關(guān)。也就是說(shuō)终惑，在很多情況下，seq1-seq2-seq3和seq2-seq1-seq3`結(jié)果是不一樣的门扇，“once a gap, always gap”
最后說(shuō)了TBA軟件的不足雹有，即需要用他們定義的MAF格式作為輸入，也就是用戶得做一些格式轉(zhuǎn)化你工作臼寄。并給了一個(gè)使用案例
已知參考序列是ref1, 用于比對(duì)的序列是seq1, seq2, seq3霸奕。比對(duì)之后得到ref1.seq1.sing.maf, ref1.seq2.sing.maf, ref1.seq3.sing.maf, 這三個(gè)文件。提供一個(gè)進(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!
