使用的數(shù)據(jù)來(lái)自于一篇發(fā)在NC的擬南芥的基因組文章，文章用了minimap/miniasm 進(jìn)行組裝，然后用racon和Pilon進(jìn)行polish, 最后拼接處62 contigs 且N50 = 12.3?Mb。
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR217/003/ERR2173373/ERR2173373.fastq.gz
seqkit seqkit fq2fa ERR2173373.fastq.gz | gzip -c > ERR2173373.fasta
我用的是Canu進(jìn)行組裝，參數(shù)如下
canu \
-p ath -d Athaliana\
useGrid=true \
gridOptions="-S /bin/sh -q wangjw" \
gridEngineArrayMaxJobs=20 \
gridEngineThreadsOption="-pe openmpi THREADS" gridEngineMemoryOption="-l mem_free=MEMORY" \
minReadLength=2000 maxThreads=15 maxMemory=60G \
genomeSize=100m \
rawErrorRate=0.300 \
correctedErrorRate=0.100 \
-nanopore-raw ERR2173373.fasta.gz
Canu默認(rèn)Pacbi的rawErroRate是0.300， Nanopore是0.500。但是根據(jù)我在自己建立的基因組學(xué)群里的討論，目前nanopore的單條read的錯(cuò)誤率大概是12%，所以兩條read在overlap的時(shí)候，最差估計(jì)會(huì)有24%以上的序列差異，于是我嘗試設(shè)置了0.300. 但是由于Nanopore的錯(cuò)誤率不是完全隨機(jī)（經(jīng)群里的小伙伴告知），所以糾錯(cuò)后正確率低于Pacbio, 所以我設(shè)置了0.100. 其他參數(shù)沒(méi)有修改， 最終我拼出了360條contig，N50=4.45M。
我檢查了下最后輸出的report文件. 第一部分表明，大部分的reads都是能夠overlap。
Part II 關(guān)于多少數(shù)據(jù)用于糾錯(cuò)，以及預(yù)期留下多少數(shù)據(jù)。默認(rèn)Canu只選擇最長(zhǎng)的40X進(jìn)行糾錯(cuò)，可以用corOutCoverage=100調(diào)整成100X. 注： rescued 表示的是剩下的沒(méi)有用于糾錯(cuò)的read，他們可能是質(zhì)粒、線粒體等。Canu保留的目的是為了避免在組裝時(shí)缺失序列信息。
Part III: 省下的就是由于太短，不能用于糾錯(cuò)的部分。
最終結(jié)果，我還用MUMMER分析了以下共線性，代碼如下，
nucmer -t 20 --prefix ont2ath Athaliana.fa ath.contigs.fasta
mummerplot -p ont2ath ont2ath.delta --png --filter
基本上每條contig都主要和一條染色體存在很好的共線性，不存在contig的mis-assembly(錯(cuò)誤組裝）現(xiàn)象。
下一步的計(jì)劃
- 只Correction 不Trim 直接組裝，比較組裝效果
- 提高糾錯(cuò)前的錯(cuò)誤率，保持糾錯(cuò)后的0.1錯(cuò)誤，比較組裝效果
- 保持糾錯(cuò)前的錯(cuò)誤率，提高糾錯(cuò)后的錯(cuò)誤率，比較組裝效果
