「三代組裝」使用Canu對(duì)三代測(cè)序進(jìn)行基因組組裝

Canu簡(jiǎn)介

目前Canu已經(jīng)更新到2.0版本，本文用的是Cau1.6版本，一些參數(shù)可能存在變動(dòng)，請(qǐng)注意分辨。

Canu是Celera的繼任者，能用于組裝PacBio和Nanopore兩家公司得到的測(cè)序結(jié)果。

Canu分為三個(gè)步驟，糾錯(cuò)，修整和組裝，每一步都差不多是如下幾個(gè)步驟：

  • 加載read到read數(shù)據(jù)庫(kù)，gkpStore
  • 對(duì)k-mer進(jìn)行技術(shù)，用于計(jì)算序列間的overlap
  • 計(jì)算overlap
  • 加載overlap到overlap數(shù)據(jù)庫(kù)，OvlStore
  • 根據(jù)read和overlap完成特定分析目標(biāo)
    • read糾錯(cuò)時(shí)會(huì)從overlap中挑選一致性序列替換原始的噪聲r(shí)ead
    • read修整時(shí)會(huì)使用overlap確定read哪些區(qū)域是高質(zhì)量區(qū)域，哪些區(qū)域質(zhì)量較低需要修整。最后保留單個(gè)最高質(zhì)量的序列塊
    • 序列組裝時(shí)根據(jù)一致的overlap對(duì)序列進(jìn)行編排(layout), 最后得到contig。

這三步可以分開(kāi)運(yùn)行，既可以用Canu糾錯(cuò)后結(jié)果作為其他組裝軟件的輸入，也可以將其他軟件的糾錯(cuò)結(jié)果作為Canu的輸入，因此下面分別運(yùn)行這三步,并介紹重要的參數(shù)。

幾個(gè)全局參數(shù)：genomeSize設(shè)置預(yù)估的基因組大小，這用于讓Canu估計(jì)測(cè)序深度； maxThreads設(shè)置運(yùn)行的最大線程數(shù)；rawErrorRate用來(lái)設(shè)置兩個(gè)未糾錯(cuò)read之間最大期望差異堿基數(shù)；correctedErrorRate則是設(shè)置糾錯(cuò)后read之間最大期望差異堿基數(shù)，這個(gè)參數(shù)需要在 組裝 時(shí)多次調(diào)整；minReadLength表示只使用大于閾值的序列，minOverlapLength表示Overlap的最小長(zhǎng)度。提高minReadLength可以提高運(yùn)行速度，增加minOverlapLength可以降低假陽(yáng)性的overlap。

組裝實(shí)戰(zhàn)

數(shù)據(jù)下載

數(shù)據(jù)來(lái)自于發(fā)表在 Nature Communication 的一篇文章 "High contiguity Arabidopsis thaliana genome assembly with a single nanopore flow cell"。 這篇文章提供了 Arabidopsis thaliana KBS-Mac-74 的30X短片段文庫(kù)二代測(cè)序、PacBio和Nanopore的三代測(cè)序以及Bionano測(cè)序數(shù)據(jù), 由于擬南芥的基因組被認(rèn)為是植物中的金標(biāo)準(zhǔn)，因此文章提供的數(shù)據(jù)適合非常適合用于練習(xí)。根據(jù)文章提供的項(xiàng)目編號(hào)"PRJEB21270", 在European Nucleotide Archive上找到下載地址。

ENA搜索
## PacBio Sequal
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116568/bam/pb.bam
## MinION
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116595/fastq/ont.fq.gz
# Illuminia MiSeq
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116569/fastq/il_1.fq.gz
wget -c -q ftp://ftp.sra.ebi.ac.uk/vol1/ERA111/ERA1116569/fastq/il_2.fq.gz

下載的PacBio數(shù)據(jù)以BAM格式存儲(chǔ)，可以通過(guò)安裝PacBio的smrtlink工具套裝，使用其中的bam2fasta工具進(jìn)行轉(zhuǎn)換

# build index for convert
~/opt/biosoft/smrtlink/smrtcmds/bin/pbindex pb.bam &
# convert bam to fasta
~/opt/biosoft/smrtlink/smrtcmds/bin/bam2fasta -o pb pb.bam &

PacBio的smrtlink工具套裝大小為1.4G，不但下載速度慢，安裝也要手動(dòng)確認(rèn)各種我不清楚的選項(xiàng), 唯一好處就是工具很全。

運(yùn)行Canu

第一步：糾錯(cuò)。三代測(cè)序本身錯(cuò)誤率高，使得原始數(shù)據(jù)充滿了噪音。這一步就是通過(guò)序列之間的相互比較糾錯(cuò)得到高可信度的堿基。主要調(diào)整兩個(gè)參數(shù)

  • corOutCoverage: 用于控制多少數(shù)據(jù)用于糾錯(cuò)。比如說(shuō)擬南芥是120M基因組，100X測(cè)序后得到了12G數(shù)據(jù)，如果只打算使用最長(zhǎng)的6G數(shù)據(jù)進(jìn)行糾錯(cuò)，那么參數(shù)就要設(shè)置為50(120m x 50)。設(shè)置一個(gè)大于測(cè)序深度的數(shù)值，例如120，表示使用所有數(shù)據(jù)。
canu -correct \
    -p ath -d pb_ath \
    Threads=10 gnuplotTested=true\
    genomeSize=120m minReadLength=2000 minOverlapLength=500\
    corOutCoverage=120 corMinCoverage=2 \
    -pacbio-raw pb.fasta.gz

可以將上述命令保存到shell腳本中進(jìn)行運(yùn)行, nohup bash run_canu.sh 2> correct.log &.

注: 有些服務(wù)器沒(méi)有安裝gnuplot, gnuplotTested=true 可以跳過(guò)檢查。

第二步：修整。這一步的目的是為了獲取更高質(zhì)量的序列，移除可疑區(qū)域（如殘留的SMRTbell接頭).

canu -trim \
        -p ath -d pb_ath
        maxThreads=20 gnuplotTested=true\
        genomeSize=120m minReadLength=2000 minOverlapLength=500\
        -pacbio-corrected ath/pb_ath.correctedReads.fasta.gz

第三步: 組裝。在前兩步獲得高質(zhì)量的序列后，就可以正式進(jìn)行組裝. 這一步主要調(diào)整的就是糾錯(cuò)后的序列的錯(cuò)誤率， correctedErrorRate，它會(huì)影響utgOvlErrorRate。這一步可以嘗試多個(gè)參數(shù)，因?yàn)樗俣缺容^塊。

# error rate 0.035
canu -assemble \
    -p ath -d ath-erate-0.035 \
    maxThreads=20 gnuplotTested=true \
    genomeSize=120m\
    correctedErrorRate=0.035 \
    -pacbio-corrected atg/pb_ath.trimmedReads.fasta.gz
# error rate 0.050
canu -assemble \
    -p ath -d ath-erate-0.050 \
    maxThreads=20 gnuplotTested=true \
    genomeSize=120m\
    correctedErrorRate=0.050 \
    -pacbio-corrected atg/pb_ath.trimmedReads.fasta.gz

最后輸出文件下的ath.contigs.fasta就是結(jié)果文件。

一些寶貴的建議

后續(xù)分析要去冗余

Canu2.0之前的contig盡管運(yùn)行日志說(shuō)沒(méi)有bubble，其實(shí)只是它沒(méi)有檢測(cè)到而已。Canu2.0才真正的加上該信息。但作者強(qiáng)烈推薦你先用purge_dups去冗余，避免軟件難以檢測(cè)到的冗余序列存在影響后續(xù)分析。

高雜合物種組裝

對(duì)于高雜合物種的組裝，Canu建議是用 batOptions=-dg 3 -db 3 -dr 1 -ca 500 -cp 50參數(shù)盡量分出兩套單倍型，然后對(duì)基因組去冗余。

batOptions表示傳遞后續(xù)的參數(shù)給組裝軟件bogart, -dg 3 -db3降低自動(dòng)確定閾值時(shí)的錯(cuò)誤率離差(deviation)，從而更好的分開(kāi)單倍型。-dr 1 -ca 500 -cp 50會(huì)影響錯(cuò)誤組裝的拆分，對(duì)于一個(gè)模棱兩可的contig，如果至少另一條可選路徑的overlap長(zhǎng)度至少時(shí)500bp，或者說(shuō)另一條可選路徑時(shí)在長(zhǎng)度上和當(dāng)前最佳路徑存在50%的差異，那么就將contig進(jìn)行拆分。

關(guān)于雜合物種組裝的討論，參考https://github.com/marbl/canu/issues/201#issuecomment-233750764

購(gòu)買(mǎi)SSD避免服務(wù)器IO瓶頸

如果你的服務(wù)器線程數(shù)很多，你在普通的機(jī)械硬盤(pán)上運(yùn)行組裝，而且你的系統(tǒng)還是CentOS，那么你需要調(diào)整一個(gè)參數(shù)，避免其中一步的IO嚴(yán)重影響服務(wù)器性能。

Canu通過(guò)兩個(gè)策略進(jìn)行并行，bucketizing (‘ovb’ 任務(wù)) 和 sorting (‘ovs’ 任務(wù))。 bucketizing會(huì)從1-overlap讀取輸出的overlap，將他們復(fù)制一份作為中間文件。sorting一步將這些文件加載到內(nèi)存中進(jìn)行排序然后寫(xiě)出到硬盤(pán)上。 如果你的overlap輸出特別多，那么該步驟將會(huì)瞬間擠爆的你的IO.

為了避免悲劇發(fā)生，請(qǐng)?jiān)黾尤缦聟?shù): ovsMemory=16G ovbConcurrency=15 ovsConcurrency=15， 也就是降低這兩步同時(shí)投遞的任務(wù)數(shù)，緩解IO壓力。

參考資料

版權(quán)聲明：本博客所有文章除特別聲明外，均采用 知識(shí)共享署名-非商業(yè)性使用-禁止演繹 4.0 國(guó)際許可協(xié)議 (CC BY-NC-ND 4.0) 進(jìn)行許可。

掃碼即刻交流
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書(shū)系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

友情鏈接更多精彩內(nèi)容