一文看懂如何用QUAST评估组装的基因组
2021-07-12 本文已影响0人
生信阿拉丁
作者：童蒙
编辑：angelica
二代测序和三代测序技术给基因组研究带来了巨大的变革，然而二代技术的读长短、不能很好地处理重复区域；而三代技术读长长、但是准确率低；同时不同的组装软件也有各种不同的算法。为了更好地评测不同的软件、不同的数据在基因组组装方面的优势，作者升级了QUAST为QUAST-LG，增加了对三代测序和大基因组的评测。
一、依赖软件
- 长序列的比对：minimap
- 基因预测和功能比较：GeneMarkS, GeneMark-ES, GlimmerHMM, Barrnap, and BUSCO
- SV检测：BWA, Sambamba, and GRIDSS
- 深度和覆盖度统计：bedtools
- contig结果浏览器：Icarus、Circos
二、安装和使用
安装
conda一键安装，命令如下：
conda install -c bioconda quast
简易使用
./quast.py test_data/contigs_1.fasta \
test_data/contigs_2.fasta \
-r test_data/reference.fasta.gz \
-g test_data/genes.gff
- contig.fasta是组装好的结果，可以同时跑多个文件，进行比较
- -r : 参考基因组
- -g：参考的基因集合
更多的参数说明
更多的参数说明可以参考官网，由于篇幅有限， 大家可以自己去查看一下。
三、结果解释
小编自己拿人的组装结果进行了测试，各个结果统计如下：
表格内容的解释
1 Genome statistics
Genome fraction (%) 94.015
Duplication ratio 1.008
Largest alignment 24186441
Total aligned length 2741819159
NG50 4238925
NG75 692661
NA50 4339909
NA75 1231691
NGA50 3532470
NGA75 603978
LG50 166
LG75 625
LA50 166
LA75 469
LGA50 207
LGA75 756
带大家看下是些什么指标：
- Genome fraction (%) : 基因组中被组装结果覆盖到的碱基数 除以 参考序列的总长度得到的比值；位于重复区域的contig可能会比对到多个位置，因此会被重复计算。
- Duplication ratio：组装结果中可比对到基因组的碱基数 与 基因组中被覆盖到的碱基数的比值；如果组装的结果中重复序列较多，多个contig覆盖同一个基因组区域的话，这个值会大于1。这种情况可能是由于过多的估计了重复序列的拷贝数。
- Largest alignment：将组装结果同基因组进行比对，得到的最大的连续的比对的长度；这个值一般会同largest contig相同，但是如果larget contig有些missassembled，或者部分未比对上去的话，会相对小一些。
- Total aligned length：组装结果中可比对到基因组的碱基数，由于存在未比对和部分未比对，这个数字一般会小于组装的结果的总长度（total length)。
- NG50：类似于N50的概念，指的是将组装结果从大到小排序，当第1到N个contig的总长度大于等于50%的基因组参考序列时的contig的长度，而N50用的则以50%的组装结果作为标准。
- NG75：同NG50类似，将50%换成75%而已。
- NA50：同N50类似，不同的是，不是对组装结果进行排序，而是对比对上的block进行排序，来进行统计大于等于50%的组装结果长度时的block的长度 。A代表alignment。当一个contig比对上多个位置的时候，contig会被打断成多个block。
- NA75：同NA50类似，将50%换成75%而已。
- NGA50：对比对上的block进行排序，使用的是block的长度，来统计 大于等于50%的基因组参考序列时block的长度。
- NGA75：同NGA50类似，将50%换成75%而已。
- LG50：同NG50配套，当达到NG50时，contig的个数。
- LG75：同NG75配套，当达到NG75时，contig的个数。
- LA50：同NA50配套，当达到NA50时，contig的个数。
- LA75：同NA75配套，当达到NA75时，contig的个数。
- LGA50：同NGA50配套，当达到NGA50时contig的个数。
- LGA75：同NGA75配套，当达到NGA75时contig的个数。
2 Reads mapping
# mapped 3254077
Mapped (%) 99.98
# properly paired 0
Properly paired (%) 0
# singletons 0
Singletons (%) 0
# misjoint mates 0
Misjoint mates (%) 0
Avg. coverage depth 14
Coverage >= 1x (%) 99.95
Coverage >= 5x (%) 97.32
Coverage >= 10x (%) 71.41
- mapped: 命令行中，提供的序列比对上的条目上；
- Mapped (%)：比对率；
- properly paired ：能够正确比对（方向和插入片段大小都正确）的reads的个数；
- Properly paired (%)：# properly paired的序列的比例；
- singletons : 单端比对上的reads的个数；
- Singletons (%) ：# singletons 的比例；
- misjoint mates：成对的两条序列，两条比对到不同的contig；
- Misjoint mates (%)：# misjoint mates的比例；
- Avg. coverage depth：平均的覆盖深度；
- Coverage >= 1x (%) ：针对组装的结果，1x覆盖度的比例；
- Coverage >= 5x (%) ：针对组装的结果，5x覆盖度的比例；
- Coverage >= 10x (%) ：针对组装的结果，10x覆盖度的比例；
3 Misassemblies
Misassemblies
# misassemblies 1412
# relocations 1070
# translocations 286
# inversions 56
# misassembled contigs 740
Misassembled contigs length 1038454760
# local misassemblies 5745
# scaffold gap ext. mis. 0
# scaffold gap loc. mis. 0
# possible TEs 296
# unaligned mis. contigs 306
- misassemblies : Quast检测四种misassemblies，分别为：
- relocation：将contig同reference进行比对，左边的block与右边的block的距离超过了1K，或者二者重叠的距离大于1K
- inversion：两个比对的blcok位于不同的strand上
- translocation：两个比对的block位于不同的染色体上
- misassembled contigs: 包含misassemblies事件的contig数
- Misassembled contigs length ：包含misassemblies事件的contig的总得长度
- local misassemblies : 满足以下两个条件的可以称之为local misassemblies，一个是gap或者overlap小于1K且大于indel的长度标准（85bp）的；左右两个block都在同一个strand和染色体上；
- scaffold gap ext. mis.：组装结果中，scaffold比对到参考基因组，开gap长度的错误数，
- scaffold gap loc. mis. : 组装结果中，scaffold比对到参考基因组，符合local misassemblies的要求，开gap长度的错误数
- possible TEs：TE的数目。由于TE导致misassemblies的事件数。
- unaligned mis. contigs : 一个contig存在50%的unaligned的碱基，同时包含一个misassembly事件时，我们认定为unaligned mis. contigs。
4 Unaligned
# fully unaligned contigs 1428
Fully unaligned length 46340452
# partially unaligned contigs 1721
Partially unaligned length 26108910
- fully unaligned contigs：没有比对上参考基因组的contig数
- Fully unaligned length：没有比对上参考基因组的contig总长度
- partially unaligned contigs：部分没有比对上基因组的contig数，长度大于参数值（500bp）
- Partially unaligned length：部分没有比对上基因组的contig中总的没有比对上的碱基数
5 Mismatches
# mismatches 3515608
# indels 969582
Indels length 8739555
# mismatches per 100 kbp 129.06
# indels per 100 kbp 35.6
# indels (<= 5 bp) 835402
# indels (> 5 bp) 134180
# N's 0
# N's per 100 kbp 0
- mismatches ：所有比对上的碱基中，mismatch的个数
- indels ：所有比对上的碱基中，indel的个数
- Indels length ：indel的总长度
- mismatches per 100 kbp：每100kbp中，mismatch的个数
- indels per 100 kbp : 每100kbp中，indel的个数
- indels (<= 5 bp) : 长度小于5的indel个数
- indels (> 5 bp) : 长度大于5bp的indel个数
- N's ：N碱基的个数
- N's per 100 kbp：每100kbp中N碱基的个数
6 Statistics without reference
# contigs 5649
# contigs (>= 0 bp) 5667
# contigs (>= 1000 bp) 5666
# contigs (>= 5000 bp) 5585
# contigs (>= 10000 bp) 4913
# contigs (>= 25000 bp) 3635
# contigs (>= 50000 bp) 2663
Largest contig 25250014
Total length 2817344836
Total length (>= 0 bp) 2817383104
Total length (>= 1000 bp) 2817382316
Total length (>= 5000 bp) 2817082622
Total length (>= 10000 bp) 2812071586
Total length (>= 25000 bp) 2790868742
Total length (>= 50000 bp) 2756006859
N50 5380095
N75 1473207
L50 133
L75 384
GC (%) 40.74
这个就很好理解了，我就简单地挑个写下，比如contigs ：contig的总的个数。
7 Predicted genes
# predicted genes (unique) 201597
# predicted genes (>= 0 bp) 1313830 + 0 part
# predicted genes (>= 300 bp) 99698 + 0 part
# predicted genes (>= 1500 bp) 3122 + 0 part
# predicted genes (>= 3000 bp) 436 + 0 part
- predicted genes(unique)：预测到的基因个数
图片解释
报告最后一部分是紧跟着几个累积分布图，将contig按照大到小排序后，统计不同的指标，可以很直观地看到各个分布下的解释。由于这部分很直观，就不再解释了。
参考资料
- http://quast.sourceforge.net/docs/manual.html#sec2
- Alla Mikheenko, Andrey Prjibelski, Vladislav Saveliev, Dmitry Antipov, Alexey Gurevich, Versatile genome assembly evaluation with QUAST-LG, Bioinformatics, Volume 34, Issue 13, 01 July 2018, Pages i142–i150, https://doi.org/10.1093/bioinformatics/bty266
