使用ALLMAPS進行輔助組裝

簡介

在從頭組裝過程中，確定基因組的scaffolds/contig的順序和朝向是重建染色體非常關鍵的一步。這一步可以由多種輔助組裝策略完成：例如遺傳圖譜, Hi-C, BioNano光學圖譜，10X Chicago 。

一個物種可能會有多個遺傳圖譜绳军，可以是不同項目中的不同定位群體結果，可以是不同軟件如R/QTL, MSTMAP和JOINMAP的分析結果矢腻。遺傳圖譜會因重組率门驾，偏分離(segregation distortion) , PAV(presence-absence variation)和染色體比對多態(tài)位點不同而發(fā)生變化。每一種圖譜都能夠提供不同的證據(jù)多柑，舉個例子奶是，一個scaffold可能在一個圖譜中無法被錨定，但是在另一個圖譜中可以進行錨定，將這些圖譜進行整合就能提最后染色體組裝的精確度聂沙。

如果只用一個圖譜秆麸，對scaffold進行排序只是計算量大一點而已，你需要根據(jù)圖譜中分子標記在每個scaffold的平均距離進行排序就行及汉。ALLMAPS, 正如名字說的那樣沮趣，就是能夠使用所有的圖譜證據(jù)的工具，它能夠計算scaffold的朝向坷随，使其和已有的圖譜的共線性關系最大化兔毒。他有如下亮點：

可重復性：清晰的可計算目標使得讓多種輸入圖譜的共線性關系能夠最大化
靈活性：允許為輸入的圖譜設置權重，更好的處理沖突
強大性：能使用多種遺傳圖譜甸箱，只需要做最小的轉換
易用性：讓基因組構建(FASTA和AGP輸出)和基因組升級(CHAIN輸出)流程化

通常育叁，解決scaffold順序和朝向(OO)是一種NP問題。因為基因組組裝和遺傳圖譜中都可能存在一些錯誤芍殖，我們的目標只能是找到一個近似解豪嗽。ALLMAPS將該問題轉換成旅行商問題，然后用遺傳算法優(yōu)化scaffold OO. 使用遺傳算法優(yōu)化是為了避免在局部最優(yōu)上出現(xiàn)瓶頸豌骏。遺傳圖譜最常見的錯誤是倒置和異位(inversion and translocation)

安裝可以用CONDA龟梦。

conda install jcvi

實際演示

先現(xiàn)在一個測試數(shù)據(jù)集，用于練習了解軟件的總體流程窃躲。比較尷尬的是计贰，這個數(shù)據(jù)集放在了Dropbox上，我們都知道這是一個國內不存在的網站蒂窒，所以我把他轉存到了我的騰訊微云上,鏈接：https://share.weiyun.com/5nwjljN .

第一步：準備輸入文件躁倒。首先你得要提供物理圖譜和遺傳圖譜的對應關系，格式為

Scaffold ID, scaffold position, LG, genetic position

你可以在excel表格中進行操作洒琢，然后另存為csv為文件秧秉，如下所示。你可以發(fā)現(xiàn)有一些scaffold里就只有一個遺傳標記進行對應衰抑，你可以思考下該scaffold最后會如何處理

EXCEL編輯

當然象迎，我們解壓縮后的zip文件里就有這些內容，你可以用一些命令工具(less, head, tail)看下數(shù)據(jù)是如何存放的呛踊。

第二步: 將兩個圖譜進行合并砾淌，最后會得到一個權重文件(weights.txt)和輸入的bed文件

python -m jcvi.assembly.allmaps merge JMMale.csv JMFemale.csv -o JM-2-test.bed

第三步：對權重文件"weights.txt" 進行調整。weights.txt默認每個輸入的圖譜的權重都是1谭网。作者有一個建議就是汪厨，你通過檢查最后的報告和診斷圖，有監(jiān)督地來重新對每個遺傳圖譜進行權重賦值蜻底。

$ cat weights.txt 
JMFemale 1
JMMale 1

第四步：對scaffold進行排序骄崩，搭建成準染色體水平

python -m jcvi.assembly.allmaps path -w weights.txt JM-2-test.bed scaffolds.fasta

結果解讀

上面第四步會在當前文件夾下生成許多結果文件，

JM-2-test.agp
JM-2-test.bed
JM-2-test.chain
JM-2-test.chr.agp
JM-2-test.chr.fasta
JM-2-test.fasta
JM-2-test.fasta.sizes
JM-2-test.lifted.bed
JM-2-test.summary.txt
JM-2-test.tour
JM-2-test.unplaced.agp
JM-2-test.unplaced.fasta

可以分為如下幾類

首先是組裝后的基因組序列

"JM-2-test.fasta"
"JM-2-test.agp": 每個scaffold的順序和朝向薄辅，用于上傳Genbank
"JM-2-test.chain": 用于新舊坐標間的轉換要拂，比如說你在之前坐標注釋得到GFF文件就可以用lifOver轉換到新坐標系下

然后可視化報告。每個染色體都會得到一個對應的pdf文件站楚，可視化展示如下：左圖主要關注交叉的線脱惰，表示某些marker存在矛盾。右圖關注是否有單上升/下降趨勢窿春，圖中的斜率反應的是物理距離(x軸)相對遺傳距離(y軸)的變化拉一，可以認為是重組率。低重組率不容易確定在同一個重組區(qū)間內的scaffold的位置和朝向

白色部分為可信區(qū)間旧乞，灰色部分是存疑區(qū)間蔚润。

可視化展示

除了默認的輸出外，你還可以通過movie得到軟件運行過程每次迭代后組裝情況尺栖，不過要額外安裝ffmpeg和parallel

python -m jcvi.assembly.allmaps movie -w weights.txt JM-2-test.bed scaffolds.fasta chr23

最后是總結性報告嫡纠。盡管可視化會給你比較直觀的結果，但是具體的參數(shù)還是要看JM-2-test.summary.txt. 例如遺傳標記的密度延赌，以及有多少序列被錨定到基因組上除盏。

當然ALLMAPS還有一些比較高級的操作：

拆分嵌合contig：嵌合contig指的是一段區(qū)域能夠比對到多個連鎖群中或者染色體中，一個常見的來源就是不同染色體中的重復區(qū)域由于過于相似在組裝的時候坍縮成了一個挫以。ALLMAPS也提供split進行拆分者蠕。
估計gap長度：默認會用100個N在填充兩個scaffold連接的區(qū)域，方便Genbank識別其為未知區(qū)域掐松。你可以根據(jù)遺傳距離在不同物理位置上的對應關系踱侣，預測不同gap的近似大小然后進行填充，子命令是estimategaps大磺。
在ALLMAPS中使用多種遺傳圖譜.

光學圖譜在著絲粒區(qū)間表型不好泻仙，主要是里面的串聯(lián)重復過多，使得限制性內切酶數(shù)目減少量没，遺傳圖譜中在這個區(qū)域依舊有比較多的標記

局限性：

標記數(shù)目和密度玉转。過多計算量大，太小準確性低
基因組組裝和建圖軟件的潛在錯誤殴蹄【孔ィ基因組組裝的原始質量對ALLMAPS的影響非常大，尤其是高度片段化的結果可能會導致錯誤結果袭灯。
基因組結構特征也會有影響刺下，例如染色體倒置、異位和片段重復稽荧。

參考資料：

ALLMAPS: robust scaffold ordering based on multiple maps
https://github.com/tanghaibao/jcvi/wiki/ALLMAPS

使用ALLMAPS進行輔助組裝得到準染色體水平基因組

使用ALLMAPS進行輔助組裝得到準染色體水平基因組

使用ALLMAPS進行輔助組裝

簡介

實際演示

結果解讀