一、 基本介紹

(1) 關(guān)于PCR重復(fù)的問題

• RNA-seq一般不去重復(fù)（起始量高，擴(kuò)增數(shù)少；某些RNA表達(dá)量很高，導(dǎo)致該基因的reads的比對(duì)很可能出現(xiàn)一致的情況）
• ChIP-seq一般去重復(fù)（起始量低，擴(kuò)增數(shù)多，PCR偏好性影響定量分析）
• call SNP一般去重復(fù)（PCR擴(kuò)增數(shù)多，出現(xiàn)擴(kuò)增錯(cuò)誤，影響SNP識(shí)別置信度）
• PCR去重工具首選Picard
• 萬事無絕對(duì)，還需參考起始量和PCR擴(kuò)增數(shù)判斷是否去重復(fù)。reads mapping覆蓋均勻度可以判斷是否需要去重復(fù)。
• 根源上解決去重復(fù)問題：起始量高，循環(huán)數(shù)少，reads能長不短，能雙端不單端

(2) 簡化的處理方案

設(shè)一個(gè)基因組有A、B兩個(gè)片段，PCR后得到無論多少條reads，比如n?A + m?B條，在數(shù)據(jù)分析的時(shí)候，理想情況都只保留1條A和1條B（unique reads）用于組裝，而去掉(n-1)條A和(m-1)條B。共有(n-1)條A和(m-1)條B被當(dāng)成duplicated reads看待，盡管它們是正常PCR的正常產(chǎn)物。目前的算法其實(shí)是一個(gè)簡化的處理方案，把所有重復(fù)的reads都去掉了，留下完全不重復(fù)的reads。算法沒有能力區(qū)分“假重復(fù)”（人為造成的重復(fù)序列方面的bias）和“真重復(fù)”（天然存在的重復(fù)序列）。

二、 背景知識(shí)

測(cè)序所得到的reads是由于超聲波或者酶切斷裂得到的，因此這些reads比對(duì)到基因組上的位置是完全隨機(jī)的。那么兩個(gè)reads比對(duì)到相同位置的概率是非常低的。如果兩個(gè)reads比對(duì)情況相同或者極其相似，則很有可能是由于PCR重復(fù)所導(dǎo)致的。

在illumina測(cè)序中，通常有兩種類型的重復(fù)，分別是光學(xué)重復(fù)和PCR重復(fù)。PCR重復(fù)就是在建庫過程中PCR擴(kuò)增導(dǎo)致的重復(fù)。理想情況下，對(duì)打碎的基因組DNA，每個(gè)DNA片段僅測(cè)到一次。但如果這一步擴(kuò)增了6個(gè)cycle，那么每個(gè)DNA片段有了64份拷貝。將擴(kuò)增后所有產(chǎn)物“灑”到flowcell，來自一個(gè)DNA片段的兩個(gè)拷貝，可能會(huì)錨定在兩個(gè)bead上，經(jīng)過測(cè)序得到的這兩條read，就是PCR duplication。光學(xué)重復(fù)是由于illumina測(cè)序時(shí)照相機(jī)錯(cuò)誤地將一個(gè)簇識(shí)別為兩個(gè)或多個(gè)簇，這就會(huì)導(dǎo)致產(chǎn)生完全一樣的reads。光學(xué)重復(fù)可以利用tile的坐標(biāo)進(jìn)行去除。

可以利用reads mapping的均勻程度判斷是否具有重復(fù)。若富集位點(diǎn)周圍的reads均勻覆蓋，那么沒有重復(fù)；若富集位點(diǎn)周圍覆蓋度不均勻，某些區(qū)域猛然升高，那么很有可能需要進(jìn)行PCR去重復(fù)。

也可以從根源上減少這種重復(fù)帶來的影響。起始量很多時(shí)，不需要去重復(fù)。擴(kuò)增數(shù)很少時(shí)，15個(gè)cycle以內(nèi)，不需要去重復(fù)。雙端測(cè)序由于其兩個(gè)reads的位置矯正，有助于去除PCR重復(fù)。reads長度越長，越容易識(shí)別真正的PCR重復(fù)。

在一些NGS分析流程中（如ChIP-seq、ATAC_seq和call SNP）一般需要考慮去除PCR重復(fù)，而RNA-seq一般不去重復(fù)。

※ PCR本身就是為了產(chǎn)生重復(fù)序列的。理論上來講，不同的序列在進(jìn)行PCR擴(kuò)增時(shí)，擴(kuò)增的倍數(shù)應(yīng)該是相同的。但是由于聚合酶的偏好性，PCR擴(kuò)增次數(shù)過多的情況下，會(huì)導(dǎo)致一些序列持續(xù)擴(kuò)增，而另一些序列擴(kuò)增到一定程度后便不再進(jìn)行，也就是我們常說的PCR偏好性。這種情況對(duì)于定量分析（如ChIP-seq），會(huì)造成嚴(yán)重的影響。（ChIP-seq中，由于起始量不高，且沒有那種富集程度很高的位點(diǎn)，因此通常需要考慮去除PCR重復(fù)。）

※ 此外，PCR擴(kuò)增循環(huán)數(shù)過多，會(huì)出現(xiàn)一些擴(kuò)增偏差，進(jìn)而影響一些突變識(shí)別（比如call SNP）的置信度。PCR重復(fù)也會(huì)增大變異檢測(cè)結(jié)果的假陰和假陽率。（至于call SNP，因?yàn)橐ＷC測(cè)序深度，起始量一般都高，但由于PCR擴(kuò)增會(huì)導(dǎo)致一些序列復(fù)制錯(cuò)誤，這將嚴(yán)重影響SNP位點(diǎn)識(shí)別的置信度，因此一般需要去重復(fù)。）

※ RNA-seq由于其建庫起始量一般都很高，所以不需要去除重復(fù)，而且RNA-seq數(shù)據(jù)中經(jīng)常會(huì)出現(xiàn)某些基因的表達(dá)量十分高，這就導(dǎo)致這些基因打斷的reads的比對(duì)情況有很大概率是一致的。

※ 去除MeRIP-Seq中的重復(fù)片段read應(yīng)該對(duì)低表達(dá)基因有益，但對(duì)高表達(dá)基因卻不友好。我們一般更關(guān)心假陽性錯(cuò)誤而不是假陰性錯(cuò)誤（寧愿錯(cuò)過真正的peak也不將希望將假peak當(dāng)作真的），所以對(duì)duplicate read進(jìn)行過濾，以避免由于PCR 產(chǎn)物引起的假陽性peak，但可能會(huì)錯(cuò)過高表達(dá)轉(zhuǎn)錄本上的peak。

三、 去重工具

(1) samtools

• 如果多個(gè)reads具有相同的比對(duì)位置時(shí)，samtools rmdup將它們標(biāo)記為duplicates，然后直接將識(shí)別出來的重復(fù)reads去掉，通常只保留質(zhì)量最高的一條。
• 該方法對(duì)于以下兩種情況，有很好的去除效果：一些reads由于測(cè)序錯(cuò)誤導(dǎo)致其不完全相同；比對(duì)錯(cuò)誤導(dǎo)致不同的序列比對(duì)到相同的位置（可能性不大）。
• 該方法的缺點(diǎn)：由于samtools去重只考慮reads比對(duì)上的起始終止位置，不考慮比對(duì)情況，這種去重有時(shí)會(huì)導(dǎo)致測(cè)序信息的丟失。
• 雙端測(cè)序數(shù)據(jù)用samtools rmdup效果很差，很多人建議用picard MarkDuplicates。samtools的rmdup是直接將這些重復(fù)序列從比對(duì)BAM文件中刪除掉，而Picard的MarkDuplicates默認(rèn)情況則只是在BAM的FLAG信息中標(biāo)記出來，而不是刪除，因此這些重復(fù)序列依然會(huì)被留在文件中，只是我們可以在變異檢測(cè)的時(shí)候識(shí)別到它們，并進(jìn)行忽略。
• 目前認(rèn)為，samtools rmdup已經(jīng)過時(shí)了，應(yīng)該使用samtools markdup代替。samtools markdup與picard MarkDuplicates采用類似的策略。

(2) Picard

• 該工具的MarkDuplicates方法也可以識(shí)別duplicates。但是與samtools rmdup不同的是，該工具僅僅是對(duì)duplicates做一個(gè)標(biāo)記，只在需要的時(shí)候?qū)eads進(jìn)行去重。
• 它不僅考慮reads的比對(duì)位置，還會(huì)考慮其中的插入錯(cuò)配等情況（即會(huì)利用sam/bam文件中的CIGAR值），甚至reads的tile、lane以及flowcell。
• Picard主要考慮reads的5'端的比對(duì)位置，以及每個(gè)reads比對(duì)上的方向。因此我們可以從一定程度上認(rèn)為，5' 端的位置、方向、以及堿基比對(duì)情況相同，Picard就將這些reads中堿基比對(duì)值Q>15的看作是best pair而其他的reads則當(dāng)作是duplicate reads。甚至當(dāng)reads的長度不同時(shí)，Picard依然利用上述原理進(jìn)行去重。
• 對(duì)Picard來說，reads的5' 端信息更為重要。若duplicates是PCR重復(fù)，那么它們的序列不一定完全相同。但是由于PCR擴(kuò)增時(shí)，酶的前進(jìn)方向是5'->3'方向，PCR重復(fù)序列中5' 端的部分相似的可能性更高。

(3) sambamba

• 也可以使用sambamba操作bam文件和去除重復(fù)，據(jù)說該命令運(yùn)行比picard MarkDuplicates快30倍。

四、 使用方法

(1) bam過濾步驟

a. 去除低質(zhì)量比對(duì)（MAPQ<30）
b. 去除多重比對(duì)（一條read比對(duì)到基因組的多個(gè)位置）
c. 去除PCR重復(fù)（不同reads比對(duì)到基因組的同一位置）
d. 去除比對(duì)到黑名單區(qū)域的序列

(2) 查看bam文件

samtools view 04_bam/WTF4.IP.bam | less

第1列：序列名稱—— E100032134L1C001R0030000131
第2列：flag數(shù)字之和—— 163
第3列：參考序列名稱/比對(duì)上的染色體—— 2L（若無法比對(duì)為*）
第4列：比對(duì)的第一個(gè)位置——900488（沒有比對(duì)上為0）
第5列：比對(duì)質(zhì)量/MAPQ/比對(duì)錯(cuò)誤率的-10log10值—— 60（255代表質(zhì)量不可用）
第6列：比對(duì)的具體方式/CIGAR—— 2S145M3S
第7列：mate比對(duì)到的染色體—— =（沒有比對(duì)上或沒有mate為*，完全比對(duì)為=）
第8列：mate比對(duì)的第一個(gè)位置—— 900488（0表示該列不可用）
第9列：文庫插入片段/fragment的長度—— -152（無mate則為0）
第10列：read序列—— CGGCGCCT……
第11列：測(cè)序質(zhì)量—— DDDDDDDD……（D對(duì)應(yīng)ASCII值68，即35+33，Q35）
第12列：可選區(qū)域—— AS:i:-5 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:145  YS:i:-5 YT:Z:CP XS:A:+  NH:i:1

第2列：
4 = 這條read沒有比對(duì)上；8 = mate沒有比對(duì)上；256 = secondary比對(duì)；1024 = read是PCR或光學(xué)copy產(chǎn)生。

第5列：
實(shí)際上hisat2比對(duì)的結(jié)果MAPQ沒有實(shí)際意義，不表示-10log10比對(duì)錯(cuò)誤率，它只有0、1、60三個(gè)值。

• 60：uniquely mapped read, regardless of number of mismatches / indels
• 1：multiply mapped, perfect match or few mismatches / indels
• 0：unmapped, or multiply mapped and with lots of mismatches / indels

hisat2如何報(bào)告多重比對(duì)是可以指定的，由-k參數(shù)控制。

• -k <int>：report up to <int> alignments per read; MAPQ not meaningful
• -a/--all：report all alignments; very slow, MAPQ not meaningful

-k <int>：It searches for at most <int> distinct, primary alignments for each read. Primary alignments mean alignments whose alignment score is equal to or higher than any other alignments. The search terminates when it cannot find more distinct valid alignments, or when it finds <int>, whichever happens first. The alignment score for a paired-end alignment equals the sum of the alignment scores of the individual mates. Each reported read or pair alignment beyond the first has the SAM ‘secondary’ bit (which equals 256) set in its FLAGS field. For reads that have more than <int> distinct, valid alignments, hisat2 does not guarantee that the <int> alignments reported are the best possible in terms of alignment score. Default: 5 (linear index) or 10 (graph index). Note: HISAT2 is not designed with large values for -k in mind, and when aligning reads to long, repetitive genomes, large -k could make alignment much slower.

事實(shí)上hisat2默認(rèn)參數(shù)輸出了正向比對(duì)10個(gè)和反向互補(bǔ)比對(duì)10個(gè)，正向和反向各有一個(gè)不包含256的flag（例如1個(gè)83，1個(gè)163，9個(gè)339，9個(gè)419），即只有一個(gè)是首要比對(duì)位置，其它都是次要的，它們的MAPQ都是1，都包含tag NH:i:10（唯一比對(duì)MAPQ都是60，包含NH:i:1。NH: The number of mapped locations for the read or the pair.）。如果設(shè)置-k 1，這些多重比對(duì)就只能輸出一個(gè)結(jié)果了，MAPQ還是1（也有可能是它正向比對(duì)和反向比對(duì)各有一個(gè)結(jié)果），但是NH:i:1，且在統(tǒng)計(jì)報(bào)告里會(huì)被認(rèn)為aligned concordantly exactly 1 time。

samtools view 04_bam/WTF4.IP.bam | cut -f 2,5 | sort | uniq -c
      1   137     1
      3   141     0
      2   147     60
      1   163     0
     15   163     1
      3   163     60
      6   339     0
     89   339     1
      1   355     0
      6   355     1
      3   393     1
      1   403     0
      6   403     1
      1   409     1
      6   419     0
     89   419     1
      1   69      0
      3   77      0
      1   83      0
     15   83      1
      3   83      60
      2   99      60

(3) 用samtools view進(jìn)行過濾

samtools view  -@ 4  -h  -b  -q 30  -F 4  -F 256  example.bam  >  example.filtered.bam
#去掉比對(duì)質(zhì)量<30，沒有比對(duì)到參考序列的，多重比對(duì)的

samtools view  -@ 4  -h  -b  -q 1  -F 256  example.bam  >  example.filtered.bam
#實(shí)際對(duì)hisat2產(chǎn)生的BAM文件可以用-q 1  -F 256過濾
# samtools的-q 1 比 sambamba "not unmapped"過濾掉的序列更多，包括沒有比對(duì)上的reads、錯(cuò)誤過多的多重比對(duì)reads

#例如
ls ./04_bam/*.bam | while read id ; do (samtools view -@ 4 -h -b -q 1 -F 256 $id > ./06_filtered/$(basename $id ".bam").filtered.bam) ; done  &

-@：線程數(shù)。
-h：設(shè)定輸出sam文件時(shí)帶header信息，默認(rèn)輸出sam不帶header信息。
-b：輸出為bam格式，默認(rèn)輸出SAM格式。
-f：需要的flag。
-F：過濾掉的flag。
-q：have mapping quality >= INT

(4) 用sambamba view進(jìn)行過濾

# sambamba view  -t 4  -h  -f bam  -F "[XS] == null and not unmapped and not duplicate"  example.bam  >  example.filtered.bam   #未經(jīng)實(shí)踐的步驟

-t：線程數(shù)
-h：print header before reads (always done for BAM output)
-H：output only header to stdout (if format=bam, the header is printed as SAM)
-f：--format=sam|bam|cram|json|unpack。輸出文件的格式，默認(rèn)sam
-F：--filter=FILTER。set custom filter for alignments（[XS]是bowtie2多重比對(duì)才會(huì)有的得分）
--num-filter=NUMFILTER： filter flag bits; 'i1/i2' corresponds to -f i1 -F i2 samtools arguments; either of the numbers can be omitted

(5) 用samtools markdup去重復(fù)

按read name排序：samtools sort -n xxx.sort.bam -o xxx.sortname.bam
然后：samtools fixmate -m xxx.sortname.bam xxx.fixmate.bam
按position排序：samtools sort xxx.fixmate.bam -o xxx.sortposition.bam
最后：samtools markdup -r xxx.sortposition.bam xxx.markdup.bam #加上-r就會(huì)直接去掉重復(fù)序列

(6) 用picard MarkDuplicates去重復(fù)

# picard MarkDuplicates  REMOVE_DUPLICATES=true  MAX_FILE_HANDLES=800  VALIDATION_STRINGENCY=LENIENT  I=./06_filtered/2LF4.IP.bam  O=./06_markdup/2LF4.IP.markdup.bam  M=./06_markdup/2LF4.IP.markdup.log.txt  2>>log/picard_MarkDuplicates_log.txt  &

(7) 用sambamba markdup去重復(fù)

sambamba-markdup  [options]  <input.bam>  [<input2.bam> [...]]  <output.bam>
sambamba markdup  -t 4  -r  -p  --tmpdir=./06_markdup/tmp/  ./06_filtered/2LF4.IP.bam  ./06_markdup/2LF4.IP.markdup.bam  2>>log/sambamba_markdup_log.txt  &
# 如果出現(xiàn)Too many open files報(bào)錯(cuò)，需要通過使用ulimit -n 8000或添加--overflow-list-size=600000來解決

-t：線程數(shù)
-r：remove duplicates instead of just marking them
-p：show progressbar in STDERR
--tmpdir=TMPDIR：specify directory for temporary files
--overflow-list-size=OVERFLOW_LIST_SIZE：從哈希表中拋出的reads得到第二次機(jī)會(huì)滿足它們的pairs的溢出列表的大?。J(rèn)為200000 reads）；增加大小會(huì)減少創(chuàng)建的臨時(shí)文件的數(shù)量。

(8) 建立索引

# for i in `ls ./06_markdup/*.markdup.bam` ; do (samtools index $i) ; done  &
# sambamba去重復(fù)時(shí)會(huì)自動(dòng)建立索引

(9) 進(jìn)行統(tǒng)計(jì)

ls ./06_markdup/*.markdup.bam | while read id ; do (samtools flagstat $id > ./06_markdup/$(basename $id ".bam").stat) ; done  &
ls ./04_bam/*.bam | while read id ; do (samtools flagstat $id > ./04_bam/$(basename $id ".bam").stat) ; done  &

五、 舉個(gè)例子

(1) 用samtools view進(jìn)行過濾

ls ./04_bam/*.bam | while read id ; do (samtools view -@ 4 -h -b -q 1 -F 256 $id > ./06_filtered/$(basename $id ".bam").filtered.bam) ; done  &
#去除沒有比對(duì)上的reads、錯(cuò)誤過多的多重比對(duì)reads（-q 1）；去除多重比對(duì)次要位置（secondary，-F 256）

(2) 用sambamba markdup去重復(fù)

ulimit -n 8000
#去除PCR重復(fù)
# sambamba_markdup.sh
for i in `ls ./06_filtered/*.filtered.bam` 
do 

    i=`basename $i`
    i=${i%.filtered.bam}
    echo $i >>log/sambamba_markdup_log.txt
    
    sambamba markdup  -t 4  -r  -p  --tmpdir=./06_markdup/tmp/  ./06_filtered/$i.filtered.bam  ./06_markdup/$i.markdup.bam  2>>log/sambamba_markdup_log.txt

done

(3) 建立索引

# for i in `ls ./06_markdup/*.markdup.bam` ; do (samtools index $i) ; done  &
# sambamba去重復(fù)時(shí)會(huì)自動(dòng)建立索引

(4) 進(jìn)行統(tǒng)計(jì)

ls ./06_markdup/*.markdup.bam | while read id ; do (samtools flagstat $id > ./06_markdup/$(basename $id ".bam").stat) ; done  &
ls ./04_bam/*.bam | while read id ; do (samtools flagstat $id > ./04_bam/$(basename $id ".bam").stat) ; done  &