LongQC: A Quality Control Tool for Third Generation Sequencing Long Read Data
安装
conda install h5py conda install -c bioconda pysam conda install -c bioconda edlib conda install -c bioconda python-edlib
cd /path/to/download/ git clone https://github.com/yfukasawa/LongQC.git cd LongQC/minimap2-coverage && make cd .. LongQCPath=`pwd`
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),开启高灵敏度设置。 运行
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.json web_summary.html #可视化结果 结果解读
(Nanopore 左,PacBio 右) 这个图展示是reads长度分布,三代测序数据会显示单独的高峰。如果是转录组测序数据,会有严格的大小选择的数据(如PacBio数据)会向右移动显示更高的峰。 1. Mean read length:样品的reads预测的长度 2. N50:N50值
(图1 nanopore ,图2 PacBio) 上图将显示了per reads的QV分布,y轴为平均QV值,x轴是reads的长度集。
上图为Per read 覆盖度分布图。如果数据没有问题,应该是单峰(图1)。如果是多峰值,可以考虑是否基因组或者材料的杂合情况比较高(宏基因组数据除外)。
上图为覆盖度与reads长度的关系图,用来检查覆盖范围是否有意外波动。在基因组测序数据中,覆盖范围的波动应该在一定范围内(默认为3 σ)。如果超出这个范围,表明数据存在污染,文库质量低,PacBio超载等问题。
(Nanopore 左,PacBio 右) 正常情况下,Normal reads 比non-sense reads的显示更高的值。
这个部分显示的中值覆盖率可能比使用参考基因组进行质控的结果要小。由于将reads比对到未更正的容易出错的整个测序数据集上时,比对是不太敏感的,因此覆盖率会受到这种不太敏感的结果的影响。由于上述影响,估计的基因组/转录组大小往往大于实际大小。一般来说,更好的数据提供更好的尺寸估计。
上图使用了两种不同的方法计算GC含量,其中蓝色是对end to end全长的reads进行计算,而红色是对分块的亚数据集进行计算。 蓝色显示更清晰的分布,因为使用更长的序列使结果受到更小的偏差。由于测序平台不同也会导致相同的数据的reads长度不同,从而导致全长reads不同,而红色的分块亚数据集对长度已经标准化,所以对测序或样本差异更不敏感,更可靠。
可以用来检查特定序列的连接,如接头。如果没有人工序列(如接头),在两个图中峰值应该显示在0的位置。如果观察到类似于接头的序列,则用红色虚线标出平均长度。第一张图表示有接头,第二张图表示没有接头
这里通过DUST算法显示序列的完整性分数 背景知识三代测序的质控问题一个荧光信号被测序软件判读为某一种碱基,存在一定的概率(Pe)发生判断错误。QV就是用来衡量Pe高低的一个定量指标,其定义为: 而三代测序由于其随机分布的碱基错误率,单碱基的准确性不能直接用来衡量数据的质量,所以Pacbio Sequel测序平台的下机数据不会提供单碱基质量值QV。此外,三代测序虽然读长长,但是错误率高(~15%),而且跟二代测序所用的技术原理不同。二代测序的质控软件如FastQC(依赖QV评估测序数据)不适合三代测序的质控。当三代测序数据不提供QV,且没有参考基因组时,要怎么进行质控呢? LongQC的介绍
LongQC中的一些解释
non-sense reads是指没有比对到整个测序数据集的reads,这些reads在整个数据中的分布是任意且分散的,开发者认为这些reads是测序技术或建库中导致的错误reads,或者是污染的reads。通过计算这些non-sense reads在整个测序数据集所占的比例,可以评估测序数据的质量。
由于Sequel测序数据没有提供碱基质量信息,所以在longQC质控结果中关于QV的评估都是0。而Nanopore测序数据有Q-scores,Q-Scores是用来说明一个给定的碱基比它周围的碱基更好还是更差的依据。开发者通过DASCRUBBER软件对数据生成了QV,再进行质控后可以看到QV的解析结果。 |
|