「BioNano系列」下機原始資料過濾和評估

從這部分開始，就開始涉及一些軟體的操作和資料分析，是以在進入正文之前，我們需要準備好環境。

環境準備

第一步：從

https://bionanogenomics.com/library/datasets/

下載下傳人類測試資料集，以及對應的NA12878人類基因組。

wget http://bnxinstall.com/publicdatasets/DLS/20180413_NA12878_S3_compressed.tar.gz
tar xf 20180413_NA12878_S3_compressed.tar.gz
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/002/077/035/GCA_002077035.3_NA12878_prelim_3.0/GCA_002077035.3_NA12878_prelim_3.0_genomic.fna.gz
gunzip GCA_002077035.3_NA12878_prelim_3.0_genomic.fna.gz
mv GCA_002077035.3_NA12878_prelim_3.0_genomic.fna NA12878.fa           

第二步: 在

https://bionanogenomics.com/support/software-downloads/

下載下傳Solve軟體，

伺服器要滿足如下需求:

• Python=2.7.x
• perl=5.14.x或5.16.x
• R > 3.1.2，并且安裝data.table, igraph, intervals, MASS, parallel, XML, argparser
• glibc >= 2.14 和 gcc庫
• 至少有一個256G節點的記憶體，最好有一些32G記憶體的小節點

tar -zxvf Solve3.3_10252018.tar.gz           

解壓縮後裡有如下幾個檔案夾

• cohortQC: 主要是MQR運作腳本
• HybridScaffold: 單酶系統和雙酶系統混合組裝工具腳本
• Pipeline：從頭組裝的腳本
• RefAligner：用于序列聯配群組裝
• RefGenome：hg19和hg38的cmp檔案
• SVMerge: 用于合并單酶系統得到SV結果
• UTIL: 運作從頭組裝的一些實用shell腳本，可以根據需要進行修改
• VariantAnnotation: 對找到的SV進行注釋
• VCFConverter: 将SMAP和SVMerge的結果輸出成VCF格式

資料過濾

目前的主流Bionano裝置已經是Sapjyr，BNX檔案産生于Saphyr，經由Bionano Access 下載下傳到本地。

從公司拿到的是"RawMolecules.bnx"檔案， 不過我們練習用的資料是"output/all.bnx.gz"，

mkdir test
mv output/all.bnx.gz test
zcat  all.bnx.gz| head -n 20000 | grep "# Run Data" | wc -l
# 281            

我們發現發現資料集來自于281個通道。

Label SNR 過濾: 過濾信噪比較低的分子，信噪比低意味着品質差。 有如下幾個情況，不需要做Label SNR 過濾,或在你BNX檔案的"#rh"部分有"SNRFilterType"定義，就不需要過濾

• 人類樣本不需要過濾。
• Bionanao Access 1.2 以後新的圖像檢測算法得到的BNX檔案不需要SNR過濾.
• 對于AutoDetect或Irysview處理過的資料，預設會進行label SNR過濾，處理之後就是Molecules.bnx

由于我們是人類資料集，是以下面的代碼就不需要運作了，并且絕大部分情況下也用到下面的指令。

perl /opt/biosoft/Solve3.3_10252018/cohortQC/10252018/filter_SNR_dynamic.pl -i RawMolecules.bnx -o Molecules_filter.bnx -P diag_hist.pdf > snr.log &           

分子長度過濾: 過濾短與某個門檻值的分子，公司一般會隻保留100kb或120Kb以上的分子(取決于資料量，資料越多，門檻值越高)

gunzip all.bnx.gz
RefAlinger -i all.bnx -minlen 120 -merge -o output -bnx > run.log &           

在輸出的内容中，注意如下部分

Final maps=1147524, sites=46764680, length= 268845841.028kb (avg= 234.283 kb, label density= 17.395 /100kb)           

表明過濾後，還有113萬條分子，涉及到4676萬标記，總測序量為258G，标記密度是17.395/100kb，平均分子長度大于234Kb.. 測序深度等于總測序量除以基因組大小，這是人類基因組，按照3G計算，那麼深度就是80X.

過濾後平均分子長度應大于200Kb。 标記密度不能過高，過高會因分辨率不夠而無法區分，過低則無法用于比對。一般DLS在10~25 ， NRLS在8~15.

組裝評估: 在正式組裝之前，我們還需要判斷下目前資料是否滿足組裝要求。 為了擷取所需的評估參數，得将過濾後的BNX檔案和基因組模拟模切得到的CMAP進行比對

第一步: 對基因組序列模拟酶切，得到CMAP檔案

perl /opt/biosoft/Solve3.3_10252018/Pipeline/10252018/fa2cmap_multi_color.pl -i ../NA12878.fa -e cttaag 1
mv ../NA12878_CTTAAG_0kb_0labels.cmap .           

CMAP的格式比較簡單，說明如下:

第二步: 用

align_bnx_to_cmap.py

進行比對。 bionano光學圖譜比對的基本原理是基于标記的相對位置。

python /opt/biosoft/Solve3.3_10252018/Pipeline/10252018/align_bnx_to_cmap.py  \
    --prefix  human \
    --mol molecules120k.bnx  \
    --ref NA12878_CTTAAG_0kb_0labels.cmap \
    --ra /opt/biosoft/Solve3.3_10252018/RefAligner/7915.7989rel  \
    --nthreads 80  \
    --output prealign \
    --snrFilter 2 \
    --color 1            

參數說明可自行閱讀幫助說明。我們重點關注輸出結果中如下幾方面内容:

• “contigs/alignmolvref/alignmol_log_simple.txt”裡的“Fraction of aligned molecules”和"Effective coverage of reference (X)". 我們要判斷資料是否符合最低的比對率。比對率和基因組實際情況有關(組裝品質，錯誤率，重複坍縮情況)。對于人類基因可以達到90%以上，對于不怎麼完整度的基因組，即便Bionano的品質很高，比對率可能也隻有30%~40%(僅統計150 kb 的分子）

• "alignments.tar.gz", 裡面包含的三個檔案可以輸入到Bionano提供的另一個工具Access中進行可視化，注意導入要選擇"Anchor to Molecules"。

由于人類基因組足夠的大，是以需要等待一段時間才能處理完成，之後就可以對比對結果有一個更加直覺的了解。

那麼合格後的資料應該如何組裝呢？請等待後續教程！