Canu簡介
目前Canu已經(jīng)更新到2.0版本匀钧，本文用的是Cau1.6版本，一些參數(shù)可能存在變動(dòng)崖飘，請(qǐng)注意分辨榴捡。
Canu是Celera的繼任者杈女，能用于組裝PacBio和Nanopore兩家公司得到的測序結(jié)果朱浴。
Canu分為三個(gè)步驟，糾錯(cuò)达椰，修整和組裝翰蠢，每一步都差不多是如下幾個(gè)步驟：
- 加載read到read數(shù)據(jù)庫，gkpStore
- 對(duì)k-mer進(jìn)行技術(shù)啰劲，用于計(jì)算序列間的overlap
- 計(jì)算overlap
- 加載overlap到overlap數(shù)據(jù)庫梁沧，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廷支。
這三步可以分開運(yùn)行，既可以用Canu糾錯(cuò)后結(jié)果作為其他組裝軟件的輸入栓辜，也可以將其他軟件的糾錯(cuò)結(jié)果作為Canu的輸入恋拍，因此下面分別運(yùn)行這三步,并介紹重要的參數(shù)。
幾個(gè)全局參數(shù)：genomeSize設(shè)置預(yù)估的基因組大小藕甩，這用于讓Canu估計(jì)測序深度施敢； maxThreads設(shè)置運(yùn)行的最大線程數(shù)；rawErrorRate用來設(shè)置兩個(gè)未糾錯(cuò)read之間最大期望差異堿基數(shù)；correctedErrorRate則是設(shè)置糾錯(cuò)后read之間最大期望差異堿基數(shù)僵娃，這個(gè)參數(shù)需要在 組裝 時(shí)多次調(diào)整概作；minReadLength表示只使用大于閾值的序列，minOverlapLength表示Overlap的最小長度默怨。提高minReadLength可以提高運(yùn)行速度讯榕，增加minOverlapLength可以降低假陽性的overlap。
組裝實(shí)戰(zhàn)
數(shù)據(jù)下載
數(shù)據(jù)來自于發(fā)表在 Nature Communication 的一篇文章 "High contiguity Arabidopsis thaliana genome assembly with a single nanopore flow cell"匙睹。 這篇文章提供了 Arabidopsis thaliana KBS-Mac-74 的30X短片段文庫二代測序瘩扼、PacBio和Nanopore的三代測序以及Bionano測序數(shù)據(jù), 由于擬南芥的基因組被認(rèn)為是植物中的金標(biāo)準(zhǔn)，因此文章提供的數(shù)據(jù)適合非常適合用于練習(xí)垃僚。根據(jù)文章提供的項(xiàng)目編號(hào)"PRJEB21270", 在European Nucleotide Archive上找到下載地址集绰。
## 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ǔ)，可以通過安裝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ò)碍岔。三代測序本身錯(cuò)誤率高，使得原始數(shù)據(jù)充滿了噪音朵夏。這一步就是通過序列之間的相互比較糾錯(cuò)得到高可信度的堿基蔼啦。主要調(diào)整兩個(gè)參數(shù)
- corOutCoverage: 用于控制多少數(shù)據(jù)用于糾錯(cuò)。比如說擬南芥是120M基因組仰猖，100X測序后得到了12G數(shù)據(jù)捏肢，如果只打算使用最長的6G數(shù)據(jù)進(jìn)行糾錯(cuò)，那么參數(shù)就要設(shè)置為50(120m x 50)饥侵。設(shè)置一個(gè)大于測序深度的數(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ù)器沒有安裝gnuplot, gnuplotTested=true 可以跳過檢查辩棒。
第二步：修整。這一步的目的是為了獲取更高質(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)行日志說沒有bubble，其實(shí)只是它沒有檢測到而已染坯。Canu2.0才真正的加上該信息均芽。但作者強(qiáng)烈推薦你先用purge_dups去冗余，避免軟件難以檢測到的冗余序列存在影響后續(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)劲妙，從而更好的分開單倍型。-dr 1 -ca 500 -cp 50會(huì)影響錯(cuò)誤組裝的拆分儒喊，對(duì)于一個(gè)模棱兩可的contig镣奋，如果至少另一條可選路徑的overlap長度至少時(shí)500bp，或者說另一條可選路徑時(shí)在長度上和當(dāng)前最佳路徑存在50%的差異怀愧，那么就將contig進(jìn)行拆分侨颈。
關(guān)于雜合物種組裝的討論，參考https://github.com/marbl/canu/issues/201#issuecomment-233750764
購買SSD避免服務(wù)器IO瓶頸
如果你的服務(wù)器線程數(shù)很多芯义，你在普通的機(jī)械硬盤上運(yùn)行組裝哈垢，而且你的系統(tǒng)還是CentOS，那么你需要調(diào)整一個(gè)參數(shù)扛拨，避免其中一步的IO嚴(yán)重影響服務(wù)器性能耘分。
Canu通過兩個(gè)策略進(jìn)行并行，bucketizing (‘ovb’ 任務(wù)) 和 sorting (‘ovs’ 任務(wù))绑警。 bucketizing會(huì)從1-overlap讀取輸出的overlap求泰，將他們復(fù)制一份作為中間文件。sorting一步將這些文件加載到內(nèi)存中進(jì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 國際許可協(xié)議 (CC BY-NC-ND 4.0) 進(jìn)行許可。
