寫在前面
目前還是在繼續(xù)整一些擦邊“比較基因組”的工具。前段時間停下來，原因只有一個：我發(fā)現(xiàn) TBtools 的Simple Ka/Ks Calculator 速度確實(shí)不夠快。在分析的過程中，我們常常會對兩類基因?qū)线M(jìn)行 Ks 計(jì)算：
- 共線性分析結(jié)果：比如 MCScan 和 I-ADHoRe 的分析結(jié)果，這系列軟件已經(jīng)挖掘了共線性，于是得到的基因?qū)Ρ旧硪膊皇呛芏?，往往?w+ 級別
- blastP 結(jié)果：即相對初步很原始的，具有一定相似性的基因?qū)?，未?jīng)過共線性分析，往往是 10w+ 級別
無論哪一種，要對全部基因?qū)M(jìn)行 Ka/Ks 計(jì)算，那么對于普通 PC 來說，計(jì)算量也不是太小。主要計(jì)算分兩步：
- 雙序列的codon alignment，一般批量分析時，先做 protein alignment ，隨后基于 CDS 轉(zhuǎn)換為 codon alignment
- 進(jìn)行 ka/ks 計(jì)算，其中在兩基因間最快的最方便的最常用的是 NG86
很明顯，看起來都很簡單。但具體實(shí)踐是另一回事。
舊版本的 Simple Ka/Ks Calculator
早前，我在推文已經(jīng)提及，寫這個功能最原始的目的，還是更好的理解大體計(jì)算邏輯。當(dāng)然寫完之后，我還是寫了推文，具體見《批量計(jì)算dn/ds，粗略篩選正選擇？所有人-技能Get !》，自認(rèn)為還可以。后面課題組做相關(guān)分析其實(shí)也用上了。
當(dāng)初的實(shí)現(xiàn)比較簡單，雙序列比對上，直接調(diào)用 muscle （大家都清楚，TBtools 打包了muscle）。整體效果其實(shí)還可以。就是：
- 1w+ 基因?qū)?，需?~20 min
- 10w+ 基因?qū)Γ枰?~4 hour
從某個角度來說，靜置過夜似乎是不錯的辦法。盡管往往只有我們做實(shí)驗(yàn)的時候才會這么干。
增強(qiáng)版本
事實(shí)上，頻繁地通過命令行調(diào)用外部程序，對于程序的開銷極其大。于是，最好的解決辦法是直接用原生的雙序列比對邏輯。為了保證自己寫的邏輯正確，我反反復(fù)復(fù)其實(shí)間斷鼓搗果幾次，整體周期應(yīng)有一年以上，其實(shí)這不能怪我，具體如下：
- EMBOSS Needleman-Wunsch，表現(xiàn)最佳，調(diào)用needle，主要作為準(zhǔn)確參考
- Muscle，多序列比對軟件
- mafft，多序列比對軟件
- clustal w2
- biojava
- 我自己寫的
具體得分如下
事實(shí)上，我是有點(diǎn)意外的。當(dāng)然，我后續(xù)只做了 100 來個基因?qū)?。needle 表現(xiàn)最優(yōu)，這個跟它本身軟件開發(fā)目的有關(guān)。其次是 我寫的（其實(shí)我寫的也是 needle 算法，但是我也發(fā)現(xiàn)了有一點(diǎn)瑕疵），隨后才是 biojava 和 muscle。
各類 bioXXX 庫當(dāng)然很優(yōu)秀，我也相信他的實(shí)現(xiàn)其實(shí)很好。具體沒去看源碼。但這一年來我并不是沒有嘗試在全網(wǎng)檢索，看看是否有表現(xiàn)與 emboss needle 一樣好的。結(jié)果是，沒有。這個完全可以理解。
雙序列比對是很多人的入門算法，尤其是動態(tài)規(guī)劃。大家最多考慮是 gap+extension 的罰分。但是序列足夠多，足夠長的時候，常常會遇到連續(xù)的雙元選擇，這一步最優(yōu)的結(jié)果，可能會導(dǎo)致下一步的局部最優(yōu)。當(dāng)然，這個跟具體打分實(shí)現(xiàn)有關(guān)（后續(xù)我會繼續(xù)優(yōu)化）。Emmm，簡單的實(shí)現(xiàn)，到處都是，畢竟大家目的是理解動態(tài)規(guī)劃的概念，而不是寫一個著實(shí)可用的Code。
既然測試結(jié)果，我自己寫的并不亞于 biojava，更不亞于 原始的 調(diào)用 muscle 版本，那么就可以做增強(qiáng)。
基于測試：
- TBtools 原始版本 單線程 計(jì)算 1.6w 個基因?qū)Γ?0min
- 改進(jìn)版本 單線程計(jì)算 1.6w 個基因?qū)Γ?min
多線程
市面上應(yīng)是幾乎找不到 PC 或者 laptop 還是只有單線程的CPU的... 多線程從某種角度可以加速一些分析任務(wù)。
在使用我自己寫的雙序列比對邏輯之前，我想到的加速方案是 多線程。
但測試結(jié)果如下：
- 單線程 10000個基因?qū)Γ?2 min
- 四線程 10000個基因?qū)Γ?5 min
開了多線程，反而更慢，甚至比 15 min 還久。這個其實(shí)很好理解，因?yàn)殚_辟線程是需要開銷的，另外本身就是一個重IO的過程，調(diào)用muscle中設(shè)計(jì)系列IO操作，尤其是程序之間的通訊。Sad。
不過，現(xiàn)在既然用自己寫的雙序列比對邏輯，那么就會得到
- 單線程 10000個基因?qū)Γ? min
- 四線程 10000個基因?qū)Γ?~20 s
開心！完美，不僅單線程加速了，還可以使用多個線程，從某個角度來說，并行再加速...當(dāng)然，線程不是越多越多....
增強(qiáng)版本的 Simple Ka/Ks Calculator
相比于舊版本，增加了一個 CPU 參數(shù)：1）保持為 0 ，那么仍然會自動使用舊版本的計(jì)算邏輯，調(diào)用muscle，速度慢；2）調(diào)整為 1 或者更高，比如 4，那么會使用我自己寫的邏輯，速度快
使用實(shí)例
既然現(xiàn)在，我們可以進(jìn)行 10w+ 基因?qū)墑e的 Ks 快速計(jì)算（也就10min），那么可以做的事情就很多，比如，看看香蕉基因組的 Ks Plot，（注：直接走一個蛋白全集BlastP，計(jì)算Ks 即可，無需經(jīng)過共線性分析，如MCScan等）
看看多年前的香蕉基因組 Paper
感覺還不錯，于是我們還可以看看擬南芥的
其實(shí)也很好。
還有更多
比如，在擬南芥上，你可以畫一個原始的Blast結(jié)果，同時還畫一個共線性分析之后的結(jié)果
大體參數(shù)
注：前面我們看到了，Ks分布大體就是兩個Peaks，0.6 和 1.8。
總的來說，共線性分析之后確實(shí)看起來清晰了很多，盡管仁者見仁。
榕樹基因組的情況
正好前幾天測試 “Find Best Homology” 的時候，下載了榕樹基因組文章的數(shù)據(jù)。那就跑跑看。
看看榕樹基因組文章原本的點(diǎn)圖
然后再看看 TBtools 目前版本的點(diǎn)圖（整體跑下來，從輸入數(shù)據(jù)到輸出，大概 20 min）
比較一下，就會發(fā)現(xiàn)，基因組文章原圖中，可以非常明顯的看出兩個榕樹物種的共線性情況（應(yīng)是進(jìn)行了過濾，只相似度最高的基因?qū)M(jìn)行可視化，但具體就不清楚了）。而 TBtools 輸出的這個圖，其實(shí)也還可以。
拿到這個圖，我一眼下去是覺得有多倍化事件的痕跡。于是找了課題組的沉睡的小五郎，聊了下這個圖，也聊了下一些基因組分析的雜事。
針對這個圖，我們一致的觀點(diǎn)是：圖中顯示的應(yīng)該似乎雙子葉共有的古三倍化事件?？梢钥吹?，兩個榕樹物種基因組共線性很好（藍(lán)色點(diǎn)線），同時存在一個基因又另外對應(yīng)兩個位置（紅色點(diǎn)先）。
通過 Ks 著色，可以大體判斷基因之間的分化時間。事實(shí)上，我是挺想畫出一個還不錯的圖的。區(qū)分不同時期的倍化事件。不過香蕉基因組似乎不給我機(jī)會。
勉強(qiáng)調(diào)了下出圖參數(shù)
多少可以看得出來，香蕉基因組有兩次獨(dú)立相對近期的全基因組復(fù)制，其中一個比較老（灰色），一個比較近（藍(lán)色），當(dāng)然，還有一次更老的，可以看到都是紅色的點(diǎn)線。
寫在后面
天下武功，唯快不破！
香蕉基因組注釋，確實(shí)比較差。開展相關(guān)研究，或許得先從基因結(jié)構(gòu)注釋矯正做起。
寫到這里，已經(jīng)忘了為啥要寫了。
不過，話說回來，如果沒有完善，簡便，高效的流程，我相信我?guī)缀醪豢赡茉谝粌蓚€小時內(nèi)做完以上的分析，而且還可視化。同時，寫了這個推文。
