原核生物基因组三代数据（pacbio/nanopore）组装

动物随随便便几个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，几分钟，几分钟。。。