动物随随便便几个G,植物动不动十多个G,而细菌基因组一般就几M,基因组结构简单,没有重复序列杂合度这些组装障碍。我现在很多文章细菌基因组开始用三代测序组装了,还有些文章是二代结合三代测序???纳税人的钱不是钱的感觉,这篇文章是测试原核生物的三代数据组装软件及评测。测试软件及数据类型持续更新~
三代reads一般长度是几K到几十K,组装植物动物基因组,三代数据的优势是很明显的。二代短测序平台存在较三代更明显的GC偏好性,以及重复序列的问题。
格式转化软件
#下载两个测试数据 hdf5 文件
wget -c -r -nd -np -k -L -p https://sra-download.ncbi.nlm.nih.gov/traces/era21/ERZ/002695/ERR2695069/ERR2695069_hdf5.tgz
wget -c -r -nd -np -k -L -p https://sra-download.ncbi.nlm.nih.gov/traces/era21/ERZ/002632/ERR2632363/ERR2632363_hdf5.tgz
# md5sum check
软件安装:
## 格式转化软件
pbh5tools
## 组装软件
# canu
数据预处理:
### hdf5文件解压
$ tar -zxf *tgz
#解压后会生成多个文件
m160607_042209_00127_c100967912550000001823217106101682_s1_p0.metadata.xml
m160607_042209_00127_c100967912550000001823217106101682_s1_p0.bas.h5
m160607_042209_00127_c100967912550000001823217106101682_s1_p0.1.bax.h5
m160607_042209_00127_c100967912550000001823217106101682_s1_p0.3.bax.h5
m160607_042209_00127_c100967912550000001823217106101682_s1_p0.2.bax.h5
### 转换格式
# ”--readType“ 可选ccs,subreads,unrolled。不了解可以看测序原理
$ python /XXX/pbh5tools/bin/bash5tools.py m160607_042209_00127_c100967912550000001823217106101682_s1_p0.bas.h5 --outFilePrefix ERR2632363.subreads --readType subreads --outType fasta
$ python /XXX/pbh5tools/bin/bash5tools.py m160706_123219_00127_c100991152550000001823221807191673_s1_p0.bas.h5 --outFilePrefix ERR2695069.subreads --readType subreads --outType fasta
hdf5解压文件:
m160607_042209_00127_c100967912550000001823217106101682_s1_p0.metadata.xml
m160607_042209_00127_c100967912550000001823217106101682_s1_p0.bas.h5
m160607_042209_00127_c100967912550000001823217106101682_s1_p0.1.bax.h5
m160607_042209_00127_c100967912550000001823217106101682_s1_p0.3.bax.h5
m160607_042209_00127_c100967912550000001823217106101682_s1_p0.2.bax.h5
待讲解。。。。。。。。
我简单统计了一下两个fasta文件, 基因组本身大小约1.8M,两个测试数据平均长度就三四K,深度分别是257X, 743X。。。。,既然数据这么多,我简单过滤一下,过滤掉小于2K的reads。
软件测评 ——基因组组装
1.canu
用过这款软件装过基因组的都知道。这款软件留给大家最深刻的印象就是:慢。由于数据量小,这一劣势在细菌基因组上没有任何表现。
组装分为三步correct,trim,assemble
###偷懒批处理
$ cat sample.list | while read a;do echo "canu -correct -p ${a} -d ${a}_cor genomeSize=1.8m -pacbio-raw ${a}.subreads.2k.fasta 1>canu.${a}.cor.log 2>canu.${a}.cor.err && \
canu -trim -p ${a} -d ${a}_trim genomeSize=1.8m -pacbio-corrected ${a}_cor/${a}.correctedReads.fasta.gz 1>canu.${a}.trim.log 2>canu.${a}.trim.err && \
canu -assemble -p ${a} -d ${a}_0.045 correctedErrorRate=0.045 genomeSize=1.8m -pacbio-corrected ${a}_trim/${a}.trimmedReads.fasta.gz 1>canu.${a}.assemble.log 2>canu.${a}.assemble.err";done >run.canu.sh
# 看一行运行脚本
$ cat run.canu.sh | head -n 1
canu -correct -p ERR2632363 -d ERR2632363_cor genomeSize=1.8m -pacbio-raw ERR2632363.subreads.2k.fasta 1>canu.ERR2632363.cor.log 2>canu.ERR2632363.cor.err && \
canu -trim -p ERR2632363 -d ERR2632363_trim genomeSize=1.8m -pacbio-corrected ERR2632363_cor/ERR2632363.correctedReads.fasta.gz 1>canu.ERR2632363.trim.log 2>canu.ERR2632363.trim.err && \
canu -assemble -p ERR2632363 -d ERR2632363_0.045 correctedErrorRate=0.045 genomeSize=1.8m -pacbio-corrected ERR2632363_trim/ERR2632363.trimmedReads.fasta.gz 1>canu.ERR2632363.assemble.log 2>canu.ERR2632363.assemble.err
### cor,trim大大减少了数据量
file format type num_seqs sum_len min_len avg_len max_len
ERR2632363_cor/ERR2632363.correctedReads.fasta.gz FASTA DNA 6,104 69,719,744 205 11,422 29,482
ERR2695069_cor/ERR2695069.correctedReads.fasta.gz FASTA DNA 41,000 146,577,915 3 3,575.1 27,783
ERR2632363_trim/ERR2632363.trimmedReads.fasta.gz FASTA DNA 6,085 63,644,442 1,005 10,459.2 27,739
ERR2695069_trim/ERR2695069.trimmedReads.fasta.gz FASTA DNA 39,487 142,538,939 1,000 3,609.8 27,740
组装结果:
file format type num_seqs sum_len min_len avg_len max_len
ERR2632363_0.045/ERR2632363.contigs.fasta FASTA DNA 1 1,760,152 1,760,152 1,760,152 1,760,152
ERR2695069_0.045/ERR2695069.contigs.fasta FASTA DNA 8 1,888,736 2,055 236,092 1,872,113
错误率:
The default is 0.045 for PacBio reads, and 0.144 for Nanopore reads.
For low coverage datasets (less than 30X), we recommend increasing [correctedErrorRate] slightly, by 1% or so.
For high-coverage datasets (more than 60X), we recommend decreasing [correctedErrorRate] slightly, by 1% or so.
ERR2695069组装结果有点不理想,在79X的情况下,尝试多个correctedErrorRate值,可得到更好组装结果
运行时间:
没有统计运行时间,在我小服务器满载的情况下,两个组装一起跑,三步大概分别耗时:约1h,几分钟,几分钟。。。