LongQC | 长读长 Pacbio/Nanopore 测序数据质控神器

LongQC: A Quality Control Tool for Third Generation Sequencing Long Read Data

• 由于PacBio RS II 、PacBio sequel 和 Nanopore测序质量欠佳，如果能对长序列进行质控，则也能对测序原始数据具有初步把握。

• 推进LongQC的原因：1. 运行简单 2. 结果展示美观。

安装

• 安装依赖包

conda install h5pyconda install -c bioconda pysamconda install -c bioconda edlibconda install -c bioconda python-edlib

• 安装longQC & minimap

cd /path/to/download/git clone https://github.com/yfukasawa/LongQC.gitcd LongQC/minimap2-coverage && makecd ..LongQCPath=`pwd`

• 由于部分python版本问题，可能还需要额外安装sklearn的python包。

pip install sklearn

参数详解

python $LongQCPath/longQC.py sampleqc -h

#usage: LongQC sampleqc [-h] -o OUT -x preset [-t] [-n NSAMPLE] [-s SUF] [-c TRIM] [--adapter_5 ADP5] [--adapter_3 ADP3] [-f] [-p NCPU] [-d] [-m MEM] [-i INDS] [-b] input #input 输入待质控数据的格式，如fasta、fastq或pbbam#-o OUT 输出文件的路径#-x preset 输入待质控数据的测序类型，以便自动提供数据的接头和olvp参数，测序类型有:pb-rs2, pb-sequel, pb-hifi, ont-ligation, ont-rapid, ont-1dsq#-t 输入转录本、RNA或CDA序列文件（如果有）#-n NSAMPLE 样本序列的数量，默认为5000#-s SUF 输出文件的前缀名字#-c TRIM 输入去掉的接头的保留路径，如果没有输入，接头不会保留。#--adapter_5 ADP5 输入5'接头#--adapter_3 ADP3 输入3'接头#-f 关闭一些sensitive的设置，能运行更快但是准确性下降#-p NCPU 输入线程数，默认4个线程，线程数要求>=4#-d, --db 输入minimap2比对的db库#-m MEM, --mem MEM 亚数据集的内存限制，指定为千兆字节（>0和<2)，默认是0.5千兆字节#-i INDS, --index INDS 给出minimap2 (-I)的索引大小，单位为bp。在小内存的服务器上运行时这个参数要减少。默认是4g。#-b, --short 对于非常短和高错误的reads(500bp)，开启高灵敏度设置。

运行

• 采用hifi测序原始数据试一试

python $LongQCPath/longQC.py sampleqc -x pb-hifi -o out_dir input_reads.bam

• 结果展示

analysis/ #minimap比对结果figs/ #分析图片logs/ #分析日志longqc_sdust.txt #QV的统计情况QC_vals_longQC_sampleqc.jsonweb_summary.html #可视化结果

结果解读

• 质控的结果概要

1. Sample name：测序样品的名字，可通过参数-s设置
2. Yield：碱基的总数量
3. Number of reads：reads数量
4. Q7 bases：碱基质量值（QV）大于等于7的碱基的比例 （由于是Pacbio 测序数据，所以并没有提供QV）
5. Longest read：最长reads的长度
6. Estimated non-sense read fraction：评估出来的non-sense reads所占比例，应该少于30%。比例高表明: 1. 测序数据有问题 2. 覆盖度不高

• 接头信息

用Edlib算法去查找数据中的接头信息。由于Pacbio数据已经自动去掉接头，所以这个部分对Pacbio数据默认是隐藏的。

1. Number of trimmed reads：在reads末端中发现的可能是接头（75%以上的一致性）的序列长度；如果出乎意料地低或者没有检测到接头，可能是在接头连接步骤上有问题。
2. Max seq identity：接头序列与序列之间的最大的一致性值。如果接头仍然存在于数据集中，这个值应该相当高(90%)。
3. Average trimmed length：接头序列在reads中的平均末端位置

• Reads长度

（Nanopore 左，PacBio 右）

这个图展示是reads长度分布，三代测序数据会显示单独的高峰。如果是转录组测序数据，会有严格的大小选择的数据（如PacBio数据）会向右移动显示更高的峰。

1. Mean read length：样品的reads预测的长度

2. N50：N50值

• Per-reads的质量评估

（图1 nanopore ，图2 PacBio）

上图将显示了per reads的QV分布，y轴为平均QV值，x轴是reads的长度集。
理想情况下，短reads和长reads的QV分布应该相似，中值(绿色横线)应该高于7。
图1是开发者用DASCRUBBERR软件计算出QV后，再进行质控的结果。
值得注意的是，图2 由PacBio测序的序列是无法计算出QV。主要是因为PacBio的测序文件质量值是由"!"作占位符，并未有实质的QV意义。

• Per read 覆盖度

上图为Per read 覆盖度分布图。如果数据没有问题，应该是单峰（图1）。如果是多峰值，可以考虑是否基因组或者材料的杂合情况比较高（宏基因组数据除外）。

LongQC使用GMM(用于基因组)或混合Gaussian 和lognorm 分布(用于转录组)自动检测这样的峰值，以便从背景水平中区别真实的峰值。

上图为覆盖度与reads长度的关系图，用来检查覆盖范围是否有意外波动。在基因组测序数据中，覆盖范围的波动应该在一定范围内(默认为3 σ)。如果超出这个范围，表明数据存在污染，文库质量低，PacBio超载等问题。
图1表明存在微弱波动（可忽略），图2是正常的分布图。

• non-sense reads

(Nanopore 左，PacBio 右)

正常情况下，Normal reads 比non-sense reads的显示更高的值。
理想情况下，正常reads的中值（橙色线）应该高于7。

• 有3种情况：

1. 低coverage:当正常reads和non-sense reads的QV中值都在绿色范围，说明数据低coverage。因为当数据覆盖度不够高，正常的reads会因为测序深度不够而比对不上，导致被归为non-sense reads。

2. 噪音数据集：当正常reads和non-sense reads的中值都低于7（都在红色区域），这就推断出数据集的噪声太大，进一步的分析可能会受到严重影响。

3. 好的数据：当non-sense reads的中值低于7（红色区域），正常reads的中值高于7（绿色区域），说明数据是准确的，如图1。而Pacbio测序数据是无法生成QV 图，且无法评估出Non-sense reads（图2）。

• 注意：

这个部分显示的中值覆盖率可能比使用参考基因组进行质控的结果要小。由于将reads比对到未更正的容易出错的整个测序数据集上时，比对是不太敏感的，因此覆盖率会受到这种不太敏感的结果的影响。由于上述影响，估计的基因组/转录组大小往往大于实际大小。一般来说，更好的数据提供更好的尺寸估计。

• GC分布

上图使用了两种不同的方法计算GC含量，其中蓝色是对end to end全长的reads进行计算，而红色是对分块的亚数据集进行计算。

蓝色显示更清晰的分布，因为使用更长的序列使结果受到更小的偏差。由于测序平台不同也会导致相同的数据的reads长度不同，从而导致全长reads不同，而红色的分块亚数据集对长度已经标准化，所以对测序或样本差异更不敏感，更可靠。

• Flanking region分析

可以用来检查特定序列的连接，如接头。如果没有人工序列（如接头），在两个图中峰值应该显示在0的位置。如果观察到类似于接头的序列，则用红色虚线标出平均长度。第一张图表示有接头，第二张图表示没有接头

• 序列完整性

这里通过DUST算法显示序列的完整性分数

背景知识

三代测序的质控问题

二代测序采用边合成边测序的方法，以Illumina为例，进行桥式扩增后，在Flow cell的固相表面上获得上百万条成簇分布的双链待测片段。然后在Flow cell中加入荧光标记的dNTP和酶，由引物起始开始合成子链。由于dNTP存在 3'-OH 以叠氮基团RTG修饰，这会阻碍子链延伸，使得每个循环只能测得一个碱基。测完一个碱基后， Flow cell 通入液体洗掉多余的dNTP和酶，使用显微镜的激光扫描特征荧光信号，根据图像的颜色即可识别该碱基，再用化学试剂将叠氮基团与荧光基团洗去。重复此过程即可识别全部碱基。在进行碱基识别的同时，软件还根据预先设定的规则对每一个碱基进行质量评估，为每一个碱基赋予一个质量值。不同的二代测序平台，如果其测序原理不同，则定义质量值QV所取的参数也不同。

一个荧光信号被测序软件判读为某一种碱基，存在一定的概率(Pe)发生判断错误。QV就是用来衡量Pe高低的一个定量指标，其定义为：
QV=-10 x lg Pe
根据这一定义可知，
如果碱基识别的准确率是99%，则错误率是1%，QV是20；
如果碱基识别的准确率是99.9%，则错误率是0.1%，QV是30；
如果碱基识别的准确率是99.99%，则错误率是0.01%，QV是40
FASTQ格式的每第四行表示这条序列的质量值，用ACSII码表示。所以二代测序的质控软件可以根据质量值QV来评估测序数据的质量。

而三代测序由于其随机分布的碱基错误率，单碱基的准确性不能直接用来衡量数据的质量，所以Pacbio Sequel测序平台的下机数据不会提供单碱基质量值QV。此外，三代测序虽然读长长，但是错误率高（~15%），而且跟二代测序所用的技术原理不同。二代测序的质控软件如FastQC（依赖QV评估测序数据）不适合三代测序的质控。当三代测序数据不提供QV，且没有参考基因组时，要怎么进行质控呢？

LongQC的介绍

• LongQC质控软件可以解决以上的问题。LongQC不依赖QV和参考基因组进行质控，它依靠覆盖度（Coverage)和non-sense reads对测序数据进行评估。具体的原理和步骤如下：

1. LongQC首先将测序数据集随机分成几个亚数据集

2. 用minimap将亚数据集与整个测序数据集比对，计算overlap和coverage

3. 由于在测序文库中核酸序列随机分布，当有足够的coverage,reads跟整个测序数据集之间存在overlap.所以可以通过coverage来评估测序数据。

4. 当overlap少于2个时，就认为这个reads低coverage，将这种reads定义为non-sense reads，这些non-sense reads不能比对到整个测序数据集。

LongQC中的一些解释

• non-sense reads

non-sense reads是指没有比对到整个测序数据集的reads，这些reads在整个数据中的分布是任意且分散的，开发者认为这些reads是测序技术或建库中导致的错误reads，或者是污染的reads。通过计算这些non-sense reads在整个测序数据集所占的比例，可以评估测序数据的质量。

• 碱基质量值QV

由于Sequel测序数据没有提供碱基质量信息，所以在longQC质控结果中关于QV的评估都是0。而Nanopore测序数据有Q-scores，Q-Scores是用来说明一个给定的碱基比它周围的碱基更好还是更差的依据。开发者通过DASCRUBBER软件对数据生成了QV，再进行质控后可以看到QV的解析结果。