PCRはNGSのライブラリー調製プロトコールにおける重要なステップである。 PCR前のライブラリー中のユニークなDNA鋳型分子の数が少ない場合、またはユニークなDNAフラグメントを減らすライブラリー調製過程がある場合、いくつかのフラグメントは複数回シーケンスされる可能性がある。 これらのいわゆる"PCR duplicates"は、実験の有効なカバレッジを減少させるので望ましくない。 PCR duplicates比率が高い状態は、シーケンスデプスをあげても改善することはできないので、ライブラリー調製過程から変更する必要性があることを示している。 複数のサンプルを含む大規模なシーケンスプロジェクトでは、解析にバイアスをかける高いPCR duplicates比率を持つ異常なサンプルを特定することが重要となる。
重複リードは、リファンレンスへのアラインメント後に同じ位置にマッピングされるリード（ペアエンドの両方が同じ位置）を調べることで検出できる。この方法で検出された重複はtechnical replicates（PCR duplicatesとoptical duplicates [論文中のref 6]）と、自然なリードの重複（ カバレッジが飽和して同じ位置から偶然２度以上リードが読まれている）の両方を含む。ペアエンドリードを用いたWGSおよびエクソームシーケンケンスでは、ほとんどすべての重複はtechnical replicatesに対応している。しかし例えばRNA seqでは、発現量の高いいくつかの遺伝子は他の領域よりもはるかに多くシーケンスされるため、重複の大部分が自然重複となる。このような状況でPCR duplicatesをすべて削除すると、遺伝子発現値の推定などの下流解析にバイアスをかけてしまう。よってPCR増幅の割合を決定することが重要となる（こちらの質問を参照）。 PCR duplicatesの割合を推定するために、wetのmethodsも報告され（論文中のref 9-13）その有用性も報告されているが、ライブラリ調整に変更が必要であまり使われていない（論文中のref 9,10）。
PCRduplicatesはヘテロなバリアントのコール位置で重複するリードを使いPCR dpilicatesおよび自然な重複の相対的割合を推定するツール。ヘテロ接合部位では、同じ染色体由来のリードと異型染色体由来のリードが1:1で検出されることが予想され、自然重複が生じてもその比率は偏らないと期待されるが、PCR duplicatesがあるとそれが偏ることを調べている。 自然重複の割合が高い1000 Genomeのデータとexomeのデータを使用してPCR dplicatesの割合を見積もり、自然重複が多くてもPCR duplicates見積もり精度が損なわれないことを示している。 また、RNA-seqの重複のほんの一部（5〜30％）しかPCR dpilicatesでないことを実証するために、1000 GenomesのRNA-seqデータを分析している。
インストール
cent OSに導入した。
https://github.com/vibansal/PCRduplicates
git clone https://github.com/vibansal/PCRduplicates.git
cd PCRduplicates/
make all
> ./extract_duplicates
$ ./extract_duplicates
program options are missing, please specify at least one BAM file and a VCF file
bamfiles 0 variantfile None
PROGRAM TO estimate read duplicate cluster counts and extract reads overlapping heterozygous variants from sorted BAM files
./extract_duplicates [options] --bam reads.sorted.bam --VCF variants.vcf > output.reads
=============== PROGRAM OPTIONS ========================================
--mbq : minimum base quality to consider a base, default 13
--mmq : minimum read mapping quality to consider a read, default 30
--treatSE 0/1: consider paired-end reads as Single End (first in pair), default 0
--filterdups 0/1: discard reads marked as Duplicates (cigar flag of 2048), default 0
============================================================================
> python estimate_PCRduprate.py -h
$ python estimate_PCRduprate.py -h
Usage: estimate_PCRduprate.py [options]
Options:
-h, --help show this help message and exit
-i INPUTFILE, --input=INPUTFILE
input file with reads overlapping heterozygous
variants derived from BAM file
-f VARFILTER, --filter=VARFILTER
filter variants based on allele counts, 'exome' or
-s SCALING, --scaling=SCALING
scale the cluster counts (0/1), default =1
-r RANDOM, --random=RANDOM
random sampling of reads (fraction)
-o OUTFILE, --output=OUTFILE
file to output statistics
ラン
解析にはアライメント結果の.bamと、そのヘテロガスな変異をGATK UnifiedGenotyperやsamtoolsなどで検出したvcfファイルが必要となる。
ここでは1000genomeのNA12878.bamを分析してみる。重いのでchr1だけ使う。ヘッダーを分析して、chr1のリードだけ取り出す。
samtools view -H NA12878.bam |less
#chr1は"1"だったので、"1"にマッピングされたリードだけ抽出。
samtools view -b NA12878.bam 1 > chr1.bam
# GATK UnifiedGenotyper.vcfで変異をコール。
gatk -T UnifiedGenotyper -R GRCh38.fa -I chr1.bam -o UnifiedGenotyper.vcf
ここからPCRduplicatesのコマンド。duplicatesを抽出。
./extract_duplicates --bam chr1.bam --VCF UnifiedGenotyper.vcf > sample.hetreads
VCF file ../UnifiedGenotyper.vcf has 21765 variants
vcffile ../UnifiedGenotyper.vcf chromosomes 1 hetvariants 16767 21765
reading sorted bamfile to identify PCR duplicates ../R1R2_sorted.bam
processed 2000000 reads, useful fragments 479958, filtered-reads 214800 current read->tid 0
PCR duplicates marked 0 total-reads 2015396 frac 0.0000 discarded 438905
#clusters 1:660240 2:62153 3:17213 4:7507 5:4017 6:2494 7:1553 8:1137 9:859 10:671 11:496 12:381 13:292 14:219 15:214 16:169 17:149 18:127 19:98 20:99 21:86 22:63 23:50 24:52 25:67 26:51 27:34 28:41 29:26 30:35 31:20 32:15 33:20 34:19 35:25 36:22 37:17 38:16 39:13 40:15 41:11 42:20 43:7 44:3 45:14 46:10 47:9 48:11 49:7 50:5 51:11 52:8 53:6 54:9 55:8 56:5 57:7 58:7 59:6 60:8 61:6 62:3 63:4 64:6 65:4 66:2
total reads (PE=1) 996091 unique-reads 760972 duplicates:235119, duplication rate 0.23604
python estimate_PCRduprate.py -i sample.hetreads > sample.PCRdups
------------simple procedure for estimating PCR-rate-----------
1 74468 1.0 PCR-forclus: 0.0 dup-rate 0.0 C-total 66050.0 C-0000 66050.0 0.0
2 7395 1.5733 PCR-forclus: 0.2134 dup-rate 0.0353 C-total 5759.0 C-0000 4108.0 0.2867
3 1944 2.1465 PCR-forclus: 0.2845 dup-rate 0.0494 C-total 1310.0 C-0000 640.0 0.301
4 781 2.7198 PCR-forclus: 0.3201 dup-rate 0.0595 C-total 482.0 C-0000 198.0 0.2566
5 351 3.2931 PCR-forclus: 0.3414 dup-rate 0.0664 C-total 195.0 C-0000 73.0 0.2178
6 249 3.8663 PCR-forclus: 0.3556 dup-rate 0.0719 C-total 127.0 C-0000 33.0 0.2363
7 117 4.4396 PCR-forclus: 0.3658 dup-rate 0.0753 C-total 48.0 C-0000 12.0 0.2063
8 91 5.0128 PCR-forclus: 0.3734 dup-rate 0.0792 C-total 35.0 C-0000 12.0 0.1418
9 58 5.5861 PCR-forclus: 0.3793 dup-rate 0.0817 C-total 24.0 C-0000 5.0 0.1781
total reads-truncated 103532 unique 85454 duplicate fraction 0.174612680138 PCR-dup-est_f1 0.0745135360726 PCR-dup-est_fvar: 0.0817057595829 singletons 66050
最終出力。
grep FINAL_PCR_RATE sample.PCRdups > sample.PCRduprate.estimate
> cat sample.PCRduprate.estimate
$ cat sample.PCRduprate.estimate
FINAL_PCR_RATE 0.0726992266167 READ_DUPLICATION_RATE 0.236
引用
A computational method for estimating the PCR duplication rate in DNA and RNA-seq experiments.
Bansal V.
BMC Bioinformatics. 2017 Mar 14;18(Suppl 3):43.
