如何去除二代测序数据中的PCR duplication才科学?

最近听闻前公司和另一家公司A因为一个项目中二代测序数据PCR duplication的问题在扯皮。前公司比较两对Paired-end reads的第8~113bp区间内碱基是否完全相同(PE150);而A公司把1~150bp都完全相同的才当作PCR duplication。最终参数的分歧使得两家公司的统计结果相差接近一倍,关乎到数据量是否的满足合同要求。

而前公司这部分QC流程的参数就是我设置的。选择截取reads中间的区域进行比较的原因很简单,因为Illumina测序前面几个碱基和末端的区域通常测序质量较差,测序错误较多。如果要求两对reads完全一致,相当于600bp(150*2*2)不能出现任何测序错误。具体为什么设置成8~113bp我并没有详细去分析和论证,只是因为当时测序都外包给N公司做,而N公司内部QC流程去PCR duplication设置的范围设置的就是这个,我也就顺理成章用了同样的参数。

那么怎样去除二代测序数据中的PCR duplication才是科学的姿势呢?虽然对很多人来说应该是一个无聊的论题,但是我还挺好奇,本文姑且在这里分析一二。

PCR duplication产生的原因和特点

PCR duplication是由二代测序建库过程中的PCR产生的,其特点是插入片段完全相同。因此如果不考虑测序错误,那么测到的序列也应当完全相同。PCR duplication的产生不仅没有提供新的信息,还会影响许多下游的分析。比如通过kmer分析估算基因组的重复率,PCR duplication也会被统计到基因组的重复率中。其他下游分析比如变异检测、利用mate-pair reads进行基因组的scaffolding、Hi-C的数据分析等都会受到影响。

去除二代测序中的PCR duplication的方法

目前有很多现成的软件可以去PCR duplication,我在这里整理了几种常见的软件或方法:

1. 要求全长序列完全相同,以FASTX-Toolkit为代表

FASTX-Toolkit是一个2009年就发布的处理FASTA/FASTQ格式文件的工具包。FASTX-Toolkit中的组件fastx_collapser可用于处理PCR duplication。其原理是将FASTA/FASTQ中的序列进行比较,完全相同的序列只保留一个,和A公司的方式一致。

2. 要求5'端对齐的子序列完全相同,以FastUniq为代表

FastUniq于2012年发表在PLoS ONE上,是一个专门用来快速去除PCR duplication的的软件。相比FASTX-Toolkit仅支持SE reads,这款软件支持PE,并且在运行速度上有较大提高,声称可以在10分钟内处理8700万条reads。

FastUniq支持不同长度的reads序列比较,但是因为是通过归并排序算法(merge sort algorithm)排序后从5'开始逐个位点比较相邻序列的碱基,因此需要两条不同长度的reads在5'端是对齐的。因为需要序列完全相同或者短序列完全被长序列包含,因此和A公司的处理方式也一致。FastUniq的文章中也指出,这种方式由于只基于序列比较,对于杂合或者测序错误引起的SNP是敏感的。

3. 要求截取的部分序列相同,NextClip为代表

NextClip是一款专门针对Nextera MP文库数据的处理工具,于2013年发表在Bioinformatics上。它在处理PCR duplication上使用了kmer-based方法——将一对MP reads的前11个碱基和中间的11碱基共同构成44-mer序列进行比较。

另外有两款比较常用的二代测序数据质控软件会截取部分序列对PCR duplication进行估算——fastp使用了read1的前12bp序列和read2的前32bp序列;FastQC中所有超过75bp的序列会截短成50bp进行比较。两款软件都只是对重复序列进行估算而没有去重复的功能。

另外FastQC的文档中认为用过长的序列会因为测序错误导致低估duplication rate:

Even so, longer reads are more likely to contain sequencing errors which will artificially increase the observed diversity and will tend to underrepresent highly duplicated sequences.

4. Mapping-based method(基于参考基因组比对的方法)

与上面通过reads间比较的de novo法不同,如果有参考基因组,可以通过BWA、Bowtie等比对工具将reads比对到基因组。相同比对位点和方向的reads被认为是可能的duplicates,可以使用Samtools的rmdup或者picard的MarkDuplicates等工具去除或标记重复。

5. 几种方法的总结

这里不讨论基于比对的方法或者unique molecular identifier(UMI)。上面de novo的方法有的用到全长比较,而有的考虑到性能因素只截取少部分序列进行比较(NextClip和fastp最少,一对reads只用到44bp)。序列太短可能会产生假阳性,而序列太长则会对测序错误更敏感。FastQC的文档中就提到,用过长的序列会因为测序错误导致低估duplication rate:

Even so, longer reads are more likely to contain sequencing errors which will artificially increase the observed diversity and will tend to underrepresent highly duplicated sequences.

前面提到两家公司的争论也就在此。一方认为1~150bp过于严格,产生很多假阴性;另一方认为8~113bp过于宽松,会导致假阳性。虽然两条完全随机的序列比较,出现44bp正好相同的概率非常低,只有1/444,但是基因组中并不是完全随机的序列,重复序列会使得问题更复杂。

另外这几个截取部分序列的软件都倾向选择5'端(起始位点为1)。picard MarkDuplicates也是如此,比较的是5'端的坐标。可能是因为越接近3'端,因为InDels导致错位的概率就越高。而reads不同位点测序错误率分布的影响则不是主要关心的问题。

上面这些都只是逻辑上的推理。如何选择合适的序列长度和起始位点,可以使用基于参考基因组比对的结果作为PCR duplication的标准进行评估。

分析方法和结果

1. 测试数据的选取

选取NCBI上公开的一个大肠杆菌E. coli的Illumina HiSeq 4000测序数据和用其组装的基因组,accession number分别为SRR8203406GCA_002507455.1。因为基于参考基因组比对的方法,基因组重复序列会影响比对的准确率,使用像E. coli这种比较简单的基因组有利于提高分析的可靠性。另外这个基因组组装的连续性也还过得去,41条contigs。

2. 序列比对、标记duplicates

Reads比对和标记duplicates分别使用bwapicard进行,为减少比对对结果的影响,只保留properly paired reads:

# mapping
fq1=SRR8203406_1.fastq.gz
fq2=SRR8203406_2.fastq.gz
ref=GCA_002507455.1_ASM250745v1_genomic.fna

bwa index ${ref}
bwa mem -t 8 ${ref} ${fq1} ${fq2} > out.sam
samtools view -@ 8 -bS out.sam -o out.bam

# selection of properly paired reads
samtools flagstat -@ 8 out.bam # properly paired: 98.96%
samtools view -@ 8 out.bam -h -f 83 > out.filter.sam
samtools view -@ 8 out.bam -f 163 >> out.filter.sam
samtools view -@ 8 out.bam -f 99 >> out.filter.sam
samtools view -@ 8 out.bam -f 147 >> out.filter.sam
samtools view -bS -@ 8 out.filter.sam -o out.filter.bam
samtools flagstat -@ 8 out.filter.bam # properly paired: 100%

# filtering FASTQ files
python fastq_filter.py \ 
    out.filter.sam \
    ${fq1} ${fq2} \
    fq1.filter.fastq.gz \
    fq2.filter.fastq.gz
rm *.sam

# mark duplicates
samtools sort -@ 8 out.filter.bam -o out.sort.bam
java -jar ~/software/picard/build/libs/picard.jar MarkDuplicates \
    I=out.sort.bam \
    O=out.markdup.bam \
    M=markdup.metrics

这里只把FLAG为1187、1107、1171、1123的作为duplication(0x400),而2227、2195、2211、2163、2147、2131为嵌合序列的supplementary比对或者比对错误。

3. 真阳性vs假阳性统计

首先统计各指标的值。

分类预测阳性预测阴性
真实阳性真阳性TP假阴性FN
真实阴性假阳性FP真阴性TN

这里我们把通过picard MarkDuplicates得到的重复序列作为set1。把通过de novo序列比较预测的重复序列作为set2;以重复一次为例子,如果read对A和read对B经过比较相同,两者都统计进去;同理,如果read对A被标记为重复,根据picard MarkDuplicates的原理,相同5'坐标和方向的read对B也统计进去。得到set1和set2之后,以set1为标准衡量set2的结果,全集为setU:

TP = set1 & set2
FN = set1 - set2
FP = set2 - set1
TN = setU - (set1 | set2)

然后求真阳性率TPR、假阳性率FPR、精确率Precision和F1 score:

真阳性率TPR = Sensitivity = Recall = TP/(TP+FN)
假阳性率FPR = 1-Specificity = FP/(FP+TN)
精确率Precision = TP/(TP+FP)
F1 = 2*Precision*Recall/(Precision+Recall)

通过坐标从1开始,步长2;reads截取长度从6bp开始(双端共6*2=12bp),步长2进行穷举,计算上述4个指标,最终得到统计结果stat.txt

python result_stat.py out.markdup.bam fq1.filter.fastq.gz fq2.filter.fastq.gz

4. 结果分析

(1) TPR、FPR、Precision、F1随reads截取长度的变化

TPR、FPR、Precision、F1随reads截取长度变化

  • TPR:从整体上看TPR随着长度增加而逐渐减小,可能的原因是序列越长,测序错误越多,导致真阳性率下降。当使用全长序列进行分析时,TPR仅有0.41。
  • FPR:除了截取长度为12bp(read1和read2各6bp)时FPR高达0.20~0.26,基本上FPR都处于很低的水平,说明只要read1和read2分别取大于等于8bp序列在本例中就能保证非常好的特异性。
  • Precision:同理,除了截取长度为12bp的情况,Precision都相当高。
  • F1 score:F1是平衡Precision和Recall(TPR)的综合评价指标,由于除了截取长度在12bp以外Precision都极高,因此F1的变化趋势和TPR基本保持一致。
  • 较短截取长度(12、16、20bp)会导致FPR和Precision在接近3'端出现很大变化,在图中体现为离群点。

(2) TPR、FPR、Precision、F1随位点的变化

 TPR、FPR、Precision、F1随位点变化

这里因为较短截取长度在接近3'端的FPR/Precision变化太大,所以去除了12、16、20bp的情况。

  • 不同截取长度之间的TPR、FPR、Precision、F1随位点变化趋势非常一致。
  • 由于Precision极高,F1和TPR的变化趋势基本一致。
  • TPR、F1在5'端附近先增加后下降,而在3'端急速下降。这可能和5'端以及3'端较高的测序错误率有关。
  • TPR、F1、Precision都大体呈现随位点下降趋势,FPR则随位点升高。可能是因为越接近3'端因为InDels导致错位的概率越大。
  • FPR在5'端和3'端上升,Precision相反,但总体上变化并不大(Y轴)。

(3) 按照F1排序的Top 5

TPRFPRPrecisionF1RangeLength
0.9593341007180.002122223642360.9910878038140.9749524703313-2016
0.9589724383930.002133514468290.9910374646880.97474132059515-2216
0.9581407986310.002129449770950.9910466920230.97431598912817-2416
0.951843835310.0005920909119240.9974778130270.9741266742715-2420
0.9517061466070.0005970588753340.9974563373180.97404432419213-2220

(4) 计算两家公司的方法以及几种现有软件的指标

MethodTPRFPRPrecisionF1RangeLength
前公司 / N公司0.71033968727.722924938e-050.99955825140.83048921188-113212
A公司 / FastUniq0.41256675612.709798224e-050.99973308180.58409216431-150300
NextClip0.86388645640.00011516642450.99945839170.92674050031-11, 76-8644
fastp0.86983277250.00018787934350.99912277000.93000578061-12; 1-3244

以F1 score作为衡量标准,方法优劣排序如下:

fastp ≈ NextClip > 前公司/N公司 > A公司

这里没有比较FASTX-Toolkit和FastQC,因为这两款软件都只支持单端。

结论 & 讨论

  1. 过长的截取片段的确可能因为测序错误导致重复比例被低估。按照A公司策略选取1-150全长进行分析,TPR只有0.413,按照前公司和N公司的策略选取8-113时TPR为0.71。这可以解释为什么两家公司结算结果有如此大的差异。
  2. 在本例中,仅选取片段>=16bp就能得到不错的特异性。F1值最高时截取的序列范围为13-20,共(20-13+1)*2=16bp。
  3. 5'端和3'端测序错误可能导致重复比例被低估(TPR下降),5'端的影响并不是特别大,但3'端附近的TPR急速下降。
  4. 本例使用picard MarkDuplicates作为标准,为了防止错误比对选择了基因组重复较低的E. coli,并且只用properly paired reads构成的subset进行分析。
  5. 对于复杂基因组测序结果,重复序列比例应该要高很多,因此FPR应该比本例的要高,为了保证方法的稳健性,还是应当提高截取序列的长度,仅使用16bp可能不够。而TPR应该影响不大。

代码获取

本文分析中数据下载、比对、过滤、运行picard MarkDuplicates、以及TPR、FPR、Precision、F1计算的代码可以去我的GitHub查阅:

https://github.com/zengxiaofei/PCR-duplication

参考资料

标签: 高通量测序, 精选

知识共享许可协议 作者: 链接:https://byteofbio.com/archives/8.html
本文采用“署名-非商业性使用-相同方式分享 4.0 国际许可协议”进行许可

已有 2 条评论

  1. wu wu

    您好,文章中FN = set2 - set1的定义错了

    1. 感谢指正!文章里FN和FP写反了,已改正,代码写的是对的。

添加新评论