評(píng)估從頭組裝結(jié)果

Bionano從頭組裝出光學(xué)圖譜CMAP可以和參考序列的CMAP進(jìn)行比對(duì)，通過Access上可視化檢查參考基因組的組裝質(zhì)量，比較兩者間的不同。

這里所用的CMAP圖譜來自于一篇發(fā)表在NC的擬南芥的基因組文章(原本計(jì)劃用他們的bnx文件介紹從頭組裝，但是通訊作者根本不搭理我)，

光學(xué)圖譜的下載方式為:

wget https://submit.ncbi.nlm.nih.gov/ft/byid/w4jcevedkbs-mac-74_bng_contigs2017.cmap

我們可以根據(jù)「三代組裝」使用minimap+miniasm對(duì)nanopore進(jìn)行基因組組裝組裝出對(duì)應(yīng)的物理圖譜，

分析代碼如下:

#模擬酶切
perl /opt/biosoft/Solve3.3_10252018/Pipeline/10252018/fa2cmap_multi_color.pl -i R05C0144.fa -e BspQI 1
# 兩個(gè)圖譜比較
python /opt/biosoft/Solve3.3_10252018/Pipeline/10252018/runCharacterize.py \
    -t /opt/biosoft/Solve3.3_10252018/RefAligner/7915.7989rel/RefAligner \
    -q kbs-mac-74_bng_contigs2017.cmap -r R05C0144_BSPQI_0kb_0labels.cmap \
    -p /opt/biosoft/Solve3.3_10252018/Pipeline/10252018 \
    -a /opt/biosoft/Solve3.3_10252018/RefAligner/7915.7989rel/optArguments_nonhaplotype_noES_irys.xml \
    -n 10

運(yùn)行之后會(huì)在當(dāng)前目錄下生成一個(gè)"alignref"文件夾， 將其中的"q.cmap","r.cmap",".xmap"下載到本地，上傳到access中進(jìn)行可視化

組裝肉眼評(píng)估

上圖中，箭頭指示的部分可能就是光學(xué)圖譜能用于錨定其他contig的部分，這就是下一節(jié)光學(xué)圖譜輔助組裝的原理。

光學(xué)圖譜輔助組裝

NGM(Next-Generation Mapping) Scaffold 流程:

1. 為序列數(shù)據(jù)產(chǎn)生 in silico 圖譜
2. 將序列和Bionano基因組圖譜進(jìn)行比較，找到兩者之間的沖突并嘗試解決
3. 將不沖突的圖譜合并成 hybrid scafold
4. 在序列圖譜和hybrid scaffold之間形成聯(lián)配
5. 得到scaffold的AGP和FASTA文件

整個(gè)流程和Bionano Access完美整合，為使用者提供了方便的操作界面，用于對(duì)scafflod結(jié)果進(jìn)行可視化。流程的腳本在"/opt/biosoft/Solve3.3_10252018/HybridScaffold/10252018"

單酶系統(tǒng)

流程控制腳本為: Solve3.3_版本日期HybridScaffold/版本日期/hybridScaffold.pl, 他接受輸入文件，輸出運(yùn)行過程中的信息，產(chǎn)生輸出文件，最后得到結(jié)果描述。

有四個(gè)必須文件: FASTA格式組裝結(jié)果，CMAP格式的Bionano 基因組圖譜組裝，XML格式的配置文件， RefAligner.

perl hybridScaffold.pl 
    -n FASTA格式序列 (必須)
    -b BIonano CMAP文件 (必須)
    -c  Merge 的XML配置文件 (必須)
    -r RefAligner運(yùn)行工具路徑 (必須)
    -o 輸出文件夾 (必須)
    -B conflict filter level genome maps; 1,2 or3， 決定如何處理沖突，1表示不過濾，2表示在沖突處分割contig，3表示刪除沖突的contig，沒有-M時(shí)一定要加入
    -N conflict filter level for sequences; 1,2 or 3， 決定如何處理沖突，1表示不過濾，2表示在沖突處分割contig，3表示刪除沖突的contig，沒有-M時(shí)一定要加入
    -f 是否覆蓋之前的輸出
    -x 分別進(jìn)行hybrid scaffold 和 genome map的相互比對(duì)
    -y 為輸入的genome maps生成嵌合質(zhì)量分
    -M 輸入手工解決過沖突的文件
    -m: 如果使用了-x或-y參數(shù)，則需要輸入Bionano molecules的BNX文件
    -p 從頭組裝流程的文件路徑，如果使用了-x或, -y 選項(xiàng)，就需要加入這一項(xiàng)
    -q 從頭組裝流程的XML配置文件，如果使用了-x或, -y 選項(xiàng)，就需要加入這一項(xiàng)
    -e 從頭組裝時(shí)的噪音參數(shù)， .errbin或err文件
    -v 輸出流程版本信息

明確一點(diǎn): -c 要求的XML文件真的不是無腦用，需要修改其中fasta2cmap的enzyme部分

實(shí)際運(yùn)行案例:

cp /opt/biosoft/Solve3.3_10252018/HybridScaffold/10252018/hybridScaffold_config.xml .
# 用vim修改hybridScaffold_config.xml中的enzyme
perl /opt/biosoft/Solve3.3_10252018/HybridScaffold/10252018/hybridScaffold.pl \
    -n R05C0144.fa \
    -b kbs-mac-74_bng_contigs2017.cmap \
    -c hybridScaffold_config.xml \
    -r /opt/biosoft/Solve3.3_10252018/RefAligner/7915.7989rel/RefAligner \
    -o R05C0144 \
    -B 2 -N 2 \
    -f 

運(yùn)行過程中會(huì)輸出scaffold N50等一些參數(shù)。N50僅僅提升了1.1M，估計(jì)是作者bionano數(shù)據(jù)不夠多。

組裝的FASTA在"R05C0144/agp_fasta"文件下，而"R05C0144/hybridScaffold_archive.tar.gz"可以上傳到Access查看組裝效果, 下圖就是一個(gè)典型的混合組裝

典型的混合組裝結(jié)果

當(dāng)然具體分為哪幾步，以及每一步調(diào)用的腳本如下所示:

第一步： 將FASTA轉(zhuǎn)成CMAP格式,

Step 1

用到一個(gè)perl腳本, fa2cmap_multi_color.pl, 通過對(duì)基因組序列進(jìn)行模式搜索尋找可能的酶切位點(diǎn)，默認(rèn)輸出在"fa2cmap"文件夾下

第二步: 識(shí)別并解決沖突。

Step 2

沖突可能來自于真實(shí)的等位基因，或者時(shí)組裝錯(cuò)誤，最終的結(jié)果就是在聯(lián)配中出現(xiàn)過多無法比對(duì)上的標(biāo)記(labels). Hybrid Scaffold流程會(huì)先用RefAligner將第一步得到的cmp去跟Bionano基因組圖譜比，然后用AssignAlignType.pl識(shí)別沖突交界處。輸入文件為RefAligner運(yùn)行后得到的XMAP和CMAP文件，以及原始序列和原始Bionano基因組圖譜。統(tǒng)計(jì)每個(gè)聯(lián)配中比對(duì)和未必對(duì)標(biāo)記數(shù)，根據(jù)XML配置文件中"assignAlignType.max_overhang" 參數(shù)設(shè)置最大可以容忍的無法聯(lián)配的標(biāo)記數(shù)。最后會(huì)輸出"assginAlignType.xmap"(列出沖突位置),以及"assignAlignType_r.cmap"(無沖突序列), "assignAlignType_q.cmap"(無沖突圖譜)。更重要的是"conflicts.txt"，記錄著每個(gè)可能的位置。

之后流程用cut_conflicts.pl解決不一致的位置， 輸出"conflicts_cut_status.txt", 可以手工編輯，有監(jiān)督的進(jìn)行處理。

第三步: 合并兩者的組裝結(jié)果，形成Hybrid scaffold

Step 3

這一步用MergeNGS_BN.pl腳本完成，它會(huì)調(diào)用RefAligner進(jìn)行迭代兩兩配對(duì)合并，輸入文件是下面的其中一個(gè)

• 原始輸入
• 沖突解決后的組裝(cut_conflicts.pl輸出結(jié)果)
• 沒有沖突的組裝(AssignAlignType.pl的結(jié)果)

每一種輸入都是一種選項(xiàng)，我們可以嘗試不同的輸入，最后進(jìn)行比較。

第四步: 將序列圖圖譜和基因組圖譜比對(duì)到hybrid scaffold

Step 4

第五步: 生成hybrid scaffold表征的AGP和FASTA文件

Step 5

一些注意事項(xiàng)：

• Bionano很難處理Hi-C數(shù)據(jù)引起的基因組中朝向/排序的錯(cuò)誤。所以先Bionano混合組裝，然后才是Hi-C
• 覆蓋度: 至少50X，NLRS隨著覆蓋度提高并不會(huì)有明顯增強(qiáng)圖譜連續(xù)性，DLS(例如DLE0-1) 100X以上的覆蓋度能夠明顯提高某些植物和東西的圖譜連續(xù)性。
• 當(dāng)前的Hybrid Scaffold 流程無法很好處理單倍體信息，所以上一步的從頭組裝一定要是nonhaplotype.