细菌基因组组装初探（以伯克霍尔德菌Burkholderia为例）
随着高通量测序的迅猛发展，越来越多的物种的基因组被测序解读出来。基因组是分子生物学，遗传学，进化生物学等研究的基石，如果你研究的物种没有好的参考基因组序列，就不容易开展基因功能分析，转录组学相关，群体遗传等研究。所在实验室的主要研究物种没有好的基因组，开展实验时深感基因组的重要。不过基因组测序的成本和分析难度是一个门槛，需要大量的人力物力的投入。一般来说基因组小而简单就很容易进行测序组装注释，像细菌真菌都比较例子。所以今天以一个细菌基因组作为例子尝试和了解一下基因组测序组装。
1. 软件安装
需要下列软件
• FastQC v0.11.5
• Trimmomatic v0.36
• SPAdes v3.11.1
• MegaHit v1.1.1
• KmerGenie
• Abyss v2.1.5
• MaSuRCA
• velvet
• QUAST v4.5
• bowtie2 v2.2.5
• anvi’o v3
使用miniconda安装上述软件，miniconda的使用见我的简书教程Miniconda的安装与使用（以转录组分析软件为例）。
czh@ubuntu ~/Desktop/De_novo_genome_assembly
$ conda create -n de_novo_genome fastqc trimmomatic SPAdes MegaHit QUAST bowtie2 anvio -y
2. 获取细菌基因组原始数据
czh@ubuntu ~/Desktop/De_novo_genome_assembly
$ mkdir -p De_novo_genome_assembly & cd De_novo_genome_assembly
czh@ubuntu ~/Desktop/De_novo_genome_assembly
$ wget https://github.com/AstrobioMike/happy_belly_tutorial_data/raw/master/genomics_de_novo_working_dir.tar.gz
Length: 1086929422 (1.0G) [application/octet-stream]
Saving to: 'genomics_de_novo_working_dir.tar.gz'
genomics_de_novo_wo 100%[===================>] 1.01G 2.82MB/s in 10m 36s
2019-05-31 07:06:03 (1.63 MB/s) - 'genomics_de_novo_working_dir.tar.gz' saved [1086929422/1086929422]
czh@ubuntu ~/Desktop/De_novo_genome_assembly
$ wget https://github.com/AstrobioMike/happy_belly_tutorial_data/raw/master/genomics_de_novo_working_dir.tar.gz
Length: 516473079 (493M) [application/octet-stream]
Saving to: 'genomics_de_novo_downloaded_results.tar.gz'
genomics_de_novo_do 100%[===================>] 492.55M 749KB/s in 6m 50s
2019-05-31 07:03:29 (1.20 MB/s) - 'genomics_de_novo_downloaded_results.tar.gz' saved [516473079/516473079]
czh@ubuntu ~/Desktop/De_novo_genome_assembly
$ tar -xzvf genomics_de_novo_working_dir.tar.gz & tar -xzvf genomics_de_novo_downloaded_results.tar.gz
3. 原始数据质控
3.1 使用fastqc对原始数据进行质控
czh@ubuntu ~/Desktop/De_novo_genome_assembly
$ cd working_dir/
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir
$ ls
BCep_R1_QCd_err_cor.fastq.gz B_cepacia_raw_R2.fastq.gz
BCep_R2_QCd_err_cor.fastq.gz processing_commands.txt
B_cepacia_raw_R1.fastq.gz reference_genome
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir
$ fastqc B_cepacia_raw_R1.fastq.gz B_cepacia_raw_R2.fastq.gz -t 8
fastqc对两个双端测序的原始文件的质控结果为html格式文件：B_cepacia_raw_R1_fastqc.html，B_cepacia_raw_R2_fastqc.html。
从结果中可知， 正向序列有4149343条，总共有序列读长为151。
测序每个碱基质量图的解读，Y轴是质量评分，Q = -10log10(e)，e为测错的概率。
Relationship Between Sequencing Quality Score and Base Call Accuracy
| Quality Score | Probability of Incorrect Base Call | Inferred Base Call Accuracy |
|---|---|---|
| 10 (Q10) | 1 in 10 | 90% |
| 20 (Q20) | 1 in 100 | 99% |
| 30 (Q30) | 1 in 1000 | 99.9% |
蓝色线表示每个碱基的质量评分，红线是质量评分中位数，黄色的箱线图表示四分位差。
3.2 使用Trimmomatic切除原始数据中低质量碱基
(de_novo_genome)
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir
$ trimmomatic PE B_cepacia_raw_R1.fastq.gz B_cepacia_raw_R2.fastq.gz BCep_R1_paired.fastq.gz BCep_R1_unpaired.fastq.gz BCep_R2_paired.fastq.gz BCep_R2_unpaired.fastq.gz LEADING:10 TRAILING:10 SLIDINGWINDOW:5:20 MINLEN:151 -threads 8
TrimmomaticPE: Started with arguments:
B_cepacia_raw_R1.fastq.gz B_cepacia_raw_R2.fastq.gz BCep_R1_paired.fastq.gz BCep_R1_unpaired.fastq.gz BCep_R2_paired.fastq.gz BCep_R2_unpaired.fastq.gz LEADING:10 TRAILING:10 SLIDINGWINDOW:5:20 MINLEN:151 -threads 8
Quality encoding detected as phred33
Input Read Pairs: 4149343 Both Surviving: 598814 (14.43%) Forward Only Surviving: 657483 (15.85%) Reverse Only Surviving: 859607 (20.72%) Dropped: 2033439 (49.01%)
TrimmomaticPE: Completed successfully
参数解读
"LEADING:10" 序列开头质量评分低于10的剔除。
"TRAILING:10" 序列末端质量评分低于10的剔除。
"SLIDINGWINDOW:5:20" 滑窗参数，从第一个碱基开始，以5个碱基为滑窗计算平均质量评分低于20的位置剔。
"MINLEN:151" 检查的总序列长度不低于151。
质控结果
最后只有大约14%的双端序列剩下，大约600,000双端序列，也就是(600,000 paired reads) * (302 bps per paired read) = 181.2 Mbps数据，大多数Burkholderia菌基因组大约为 8.5 Mbps，意味着只有20X覆盖度，数据量偏少，理想的覆盖度是50-100X左右，因此需要改变质控策略，放低质控要求。
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir
$ trimmomatic PE B_cepacia_raw_R1.fastq.gz B_cepacia_raw_R2.fastq.gz BCep_R1_paired.fastq.gz BCep_R1_unpaired.fastq.gz BCep_R2_paired.fastq.gz BCep_R2_unpaired.fastq.gz CROP:140 LEADING:10 TRAILING:10 SLIDINGWINDOW:5:20 MINLEN:140 -threads 8
TrimmomaticPE: Started with arguments:
B_cepacia_raw_R1.fastq.gz B_cepacia_raw_R2.fastq.gz BCep_R1_paired.fastq.gz BCep_R1_unpaired.fastq.gz BCep_R2_paired.fastq.gz BCep_R2_unpaired.fastq.gz CROP:140 LEADING:10 TRAILING:10 SLIDINGWINDOW:5:20 MINLEN:140 -threads 8
Quality encoding detected as phred33
Input Read Pairs: 4149343 Both Surviving: 1484739 (35.78%) Forward Only Surviving: 638815 (15.40%) Reverse Only Surviving: 877835 (21.16%) Dropped: 1147954 (27.67%)
TrimmomaticPE: Completed successfully
再次质控得到大约36%的数据，150万个双端序列，大约450Mbps数据，足够50X的覆盖度了。
再用fastqc检测一下切割低质量序列后的数据
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir
$ fastqc BCep_R1_paired.fastq.gz BCep_R2_paired.fastq.gz -t 4
序列被切割成140bp的长度了，低质量序列被剔除，但数据质量依然不够完美。
4. 基因组组装
从纯二代测序从头组装基因组这篇文章和Happy Belly Bioinformatics网站中得知主流的组装软件MaSuRCA, IDBA-UD, SOAPdenovo2, Abyss, velvet，Spades和MegaHit，我逐个去这些软件的主页上去看了看。
• MaSuRCA 在2019年5月份刚更新的版本，支持二代短序列和三代长序列的混合组装。
• IDBA-UD 最近的更新是2016年7月。
• SOAPdenovo2 最近的更新是2017年1月。
• Abyss 最近的更新是2018年12月。
• velvet 最近的更新是2018年7月。
• Spades 在2019年5月份刚更新，该软件更新非常频繁。
• MegaHit 在2019年5月份刚更新，该软件更新比较频繁，软件作者之一是SOAPdenovo2的第一作者，文章强调SOAPdenovo2和IDBA-UD在组装土壤宏基因组的时候内存需求巨大，所以构建了MegaHit软件去解决这个问题。
综上，打算使用Spades, MegaHit, MaSuRCA, Abyss 和velvet这5个软件去分别组装细菌基因组。
4.1 序列错误矫正
使用spades进行序列错误矫正，可以提高基因组组装质量。
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir
$ spades.py -1 BCep_R1_paired.fastq.gz -2 BCep_R2_paired.fastq.gz -o spades_error_corrected_reads -t 12 -m 24 --only-error-correction
4.2 spades组装
尝试kmer分别为21,33,55,77,89去进行组装，--careful模式减少组装错误。
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir
$ spades.py -1 BCep_R1_QCd_err_cor.fastq.gz -2 BCep_R2_QCd_err_cor.fastq.gz -o spades_kmers_set_careful_assembly -t 4 -k 21,33,55,77,89 --careful --only-assembler
4.3 MegaHit组装
默认参数下组装
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir
$ megahit -1 BCep_R1_QCd_err_cor.fastq.gz -2 BCep_R2_QCd_err_cor.fastq.gz -o megahit_default_assembly -t 4
--min-count参数下组装
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir
$ megahit -1 BCep_R1_QCd_err_cor.fastq.gz -2 BCep_R2_QCd_err_cor.fastq.gz -o megahit_min_count_3_assembly -t 4 --min-count 3
4.4 MaSuRCA组装
参考组装工具汇总(占坑用)，这个软件强调双端数据不要进行质控切割序列和误差矫正，直接喂给软件原始数据就行。
#创建config.txt文件
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir
$ masurca -g config.txt
# 编辑config.txt中的内容
$ vim config.txt
主要分为DATA和PARAMETERS，DATA部分用来指定PE（双端illumina普通文库), JUMP（illumina大片段文库), OTHER（其他平台的测序结果)。
#config.txt文件中需要修改的行
PE= pe 350 50 /home/czh/Desktop/De_novo_genome_assembly/working_dir/B_cepacia_raw_R1.fastq /home/czh/Desktop/De_novo_genome_assembly/working_dir/B_cepacia_raw_R2.fastq
# 350为测序时插入片段大小， 50为标准差，没有大片段文库，在JUMP前加#，将其注释掉。
产生组装程序并运行
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir
$ masurca config.txt
Verifying PATHS...
jellyfish OK
runCA OK
createSuperReadsForDirectory.perl OK
creating script file for the actions...done.
execute assemble.sh to run assembly
#组装运行
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir
$ ./assemble.sh
4.5 velvet组装
参考了velvet软件进行基因组组装和使用 velvet 组装细菌基因组这两个教程。
#使用不同的Kmer值进行组装，从Kmer=31开始，以2为间隔设置Kmer值到97。
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir/velvet
$ velveth assemble 31,97,2 -shortPaired -fastq.gz -separate BCep_R1_QCd_err_cor.fastq.gz BCep_R2_QCd_err_cor.fastq.gz
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir/velvet
$ for dir in assemble_*; do velvetg $dir & done
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir/velvet
#查看不同Kmer值下组装的N50, Kmer=89时，N50最大。
$ grep -P "n50 of \d+" assemble_*/Log -o
assemble_33/Log:n50 of 813
assemble_35/Log:n50 of 921
assemble_39/Log:n50 of 1183
assemble_41/Log:n50 of 1346
assemble_43/Log:n50 of 1532
assemble_45/Log:n50 of 1734
assemble_49/Log:n50 of 2283
assemble_51/Log:n50 of 2681
assemble_53/Log:n50 of 3020
assemble_57/Log:n50 of 4125
assemble_59/Log:n50 of 4851
assemble_61/Log:n50 of 5587
assemble_63/Log:n50 of 6538
assemble_65/Log:n50 of 8029
assemble_67/Log:n50 of 9756
assemble_69/Log:n50 of 12173
assemble_71/Log:n50 of 14681
assemble_73/Log:n50 of 16829
assemble_75/Log:n50 of 19888
assemble_77/Log:n50 of 20261
assemble_79/Log:n50 of 23545
assemble_81/Log:n50 of 23545
assemble_83/Log:n50 of 26101
assemble_85/Log:n50 of 25393
assemble_87/Log:n50 of 23581
assemble_89/Log:n50 of 27885
assemble_91/Log:n50 of 25776
assemble_93/Log:n50 of 24638
assemble_95/Log:n50 of 23882
assemble_97/Log:n50 of 23798
4.6 Abyss组装
Abyss组装需要你给定一个kmer值，这里用kmer=89进行组装。实际中你需要通过循环语句，实现多个kmer 梯度组装，得到最优组装结果。KmerGenie这个软件可以帮助估算最优kmer值。
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir/abyss
$ abyss-pe k=89 name=test89 in='/home/czh/Desktop/De_novo_genome_assembly/working_dir/abyss/BCep_R1_QCd_err_cor.fastq.gz /home/czh/Desktop/De_novo_genome_assembly/working_dir/abyss/BCep_R2_QCd_err_cor.fastq.gz'
5. 不同软件组装质量的比较
QUAST是一个优秀的组装质量比较工具。QUAST 系列软件还可以比较使用不同方法组装的转录组或者宏基因组的质量。 从NCBI上下载两个关于Bulkholderia cepacia ATCC 25416 基因组的参考文件.
czh@ubuntu ~/Desktop/De_novo_genome_assembly/working_dir
$ quast.py -o quast_B_cep_out -R reference_genome/BCep_ref.fna -G reference_genome/BCep_ref.gff -l "spades_kmers_careful, megahit_default, megahit_min_count_3, MaSuRCA, velvet_89, Abyss_89" spades_kmers_set_careful_assembly/contigs.fasta megahit_default_assembly/final.contigs.fa megahit_min_count_3_assembly/final.contigs.fa masurca/CA/9-terminator/genome.scf.fasta velvet/assemble_89/contigs.fa abyss/test89-scaffolds.fa -t 8 -m 1000
5个软件中初步看来是spades组装效果最好，但是velvet和abyss是需要给定Kmer值的，我粗略的给了Kmer=89去组装，所以我这个quast组装评估结果并不可靠。此外，只利用了350bp的小片段文库的二代数据去组装，而如今上述软件大多支持二代小片段和大片段共同组装，二代三代混合组装以及纯三代pacbio和nanopore的组装，存在不同的测序策略和不同的组装策略。细节有很多，在此只是简单过一遍分析的大体流程，存在大量问题有待探索。spades这个软件还是非常适合小基因组的组装的，推荐！
参考文献
生信小白学习系列：如何进行基因组组装？(1)
De novo genome assembly Happy Belly Bioinformatics
纯二代测序从头组装基因组
velvet软件进行基因组组装
使用 velvet 组装细菌基因组
