最近處理的幾組數(shù)據(jù)使我感覺到生信軟件的教程一定要看最新的或官方的，因?yàn)楝F(xiàn)在的技術(shù)和方法發(fā)展太快了，軟件更新?lián)Q代也非?？?，每一個(gè)版本都可能有細(xì)微的區(qū)別，網(wǎng)上找到的幾年前的教程可能在某些方面已經(jīng)不適用于軟件當(dāng)前版本了。

在做一組MeRIP-seq數(shù)據(jù)的分析時(shí)，遇到了bam文件過濾和去重復(fù)的問題，并且再次看了看hisat2和samtools的一些選項(xiàng)和參數(shù)。

先是從網(wǎng)上查到的關(guān)于PCR重復(fù)的問題

• RNA-seq一般不去重復(fù)（起始量高，擴(kuò)增數(shù)少；某些RNA含量很高）
• ChIP-seq一般去重復(fù)（起始量低）
• call SNP一般去重復(fù)（PCR擴(kuò)增錯(cuò)誤影響SNP識別位點(diǎn)置信度）
• PCR去重工具首選Picard
• 萬事無絕對，還需參考起始量和PCR擴(kuò)增數(shù)判斷是否去重復(fù)。reads mapping覆蓋均勻度可以判斷是否需要去重復(fù)。
• 根源上解決去重復(fù)問題：起始量高，循環(huán)數(shù)少，reads能長不短，能雙端不單端

為什么會(huì)出現(xiàn)重復(fù)

在illumina測序中，通常有兩種類型的重復(fù)，分別是光學(xué)重復(fù)和PCR重復(fù)。

建庫中有一步是PCR擴(kuò)增加了接頭的DNA片段。理想情況下，對打碎的基因組DNA，每個(gè)DNA片段僅測到一次。但這一步擴(kuò)增了6個(gè)cycle，那么每個(gè)DNA片段有了64份拷貝。將擴(kuò)增后所有產(chǎn)物“灑”到flowcell，來自一個(gè)DNA片段的兩個(gè)拷貝，可能會(huì)錨定在兩個(gè)bead上，經(jīng)過測序得到的這兩條read，就是PCR duplication。

光學(xué)重復(fù)，是由于illumina測序時(shí)照相機(jī)錯(cuò)誤的的將一個(gè)簇識別為兩個(gè)簇（多個(gè)），這就會(huì)導(dǎo)致產(chǎn)生完全一樣的reads。該重復(fù)可以利用tail的坐標(biāo)進(jìn)行去除。

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偏好性。

設(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ù)序列）。

是否需要去重

若起始量很低，PCR擴(kuò)增次數(shù)很多，那么則需要去除PCR重復(fù)。

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

而ChIP-seq中，由于起始量不高，且沒有那種富集程度很高的位點(diǎn)，因此通常需要考慮去除PCR重復(fù)。

至于call SNP，起始量一般都高（因?yàn)橐ＷC測序深度），此外由于PCR擴(kuò)增會(huì)導(dǎo)致一些序列復(fù)制錯(cuò)誤，這將嚴(yán)重影響SNP位點(diǎn)識別的置信度。因此一般需要去重復(fù)。

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

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

去重工具

samtools

如果多個(gè)reads具有相同的比對位置時(shí)，samtools rmdup將它們標(biāo)記為duplicates，然后直接將識別出來的重復(fù)reads去掉，通常只保留質(zhì)量最高的一條。

該方法對于以下兩種情況，有很好的去除效果：一些reads由于測序錯(cuò)誤導(dǎo)致其不完全相同；比對錯(cuò)誤導(dǎo)致不同的序列比對到相同的位置（可能性不大）。

該方法的缺點(diǎn)：由于samtools去重只考慮reads比對上的起始終止位置，不考慮比對情況，這種去重有時(shí)會(huì)導(dǎo)致測序信息的丟失。

雙端測序數(shù)據(jù)用samtools rmdup效果很差，很多人建議用picard工具的MarkDuplicates功能。samtools的rmdup是直接將這些重復(fù)序列從比對BAM文件中刪除掉，而Picard的MarkDuplicates默認(rèn)情況則只是在BAM的FLAG信息中標(biāo)記出來，而不是刪除，因此這些重復(fù)序列依然會(huì)被留在文件中，只是我們可以在變異檢測的時(shí)候識別到它們，并進(jìn)行忽略。

目前認(rèn)為，samtools rmdup已經(jīng)過時(shí)了，應(yīng)該使用samtools markdup代替。samtools markdup與picard MarkDuplicates采用類似的策略。

Picard

該工具的MarkDuplicates方法也可以識別duplicates。但是與samtools不同的是，該工具僅僅是對duplicates做一個(gè)標(biāo)記，只在需要的時(shí)候?qū)eads進(jìn)行去重。

它不僅考慮reads的比對位置，還會(huì)考慮其中的插入錯(cuò)配等情況（即會(huì)利用sam/bam文件中的CIGAR值），甚至reads的tail、lane以及flowcell。Picard主要考慮reads的5'端的比對位置，以及每個(gè)reads比對上的方向。

因此我們可以從一定程度上認(rèn)為，5' 端的位置、方向、以及堿基比對情況相同，Picard就將這些reads中堿基比對值Q>15的看作是best pair而其他的reads則當(dāng)作是duplicate reads。甚至當(dāng)reads的長度不同時(shí)，Picard依然利用上述原理進(jìn)行去重。

對Picard來說，reads的5' 端信息更為重要。若duplicates是PCR重復(fù)，那么它們的序列不一定完全相同。但是由于PCR擴(kuò)增時(shí)，酶的前進(jìn)方向是5'->3'方向，PCR重復(fù)序列中5' 端的部分相似的可能性更高。

sambamba

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

過濾bam文件和去重復(fù)

SAM（sequence Alignment/mapping）數(shù)據(jù)格式是目前高通量測序中存放比對數(shù)據(jù)的標(biāo)準(zhǔn)格式。BAM是SAM的二進(jìn)制格式。

bam文件優(yōu)點(diǎn)：bam文件為二進(jìn)制文件，占用的磁盤空間比sam文本文件??；利用bam二進(jìn)制文件的運(yùn)算速度快。

舉個(gè)栗子：
序列名稱——E100032134L1C001R0142515642
flag——419
參考序列名稱/比對上的染色體——2L（無法比對為*）
比對的第一個(gè)位置——277（沒有比對上為0）
比對質(zhì)量/MAPQ/比對錯(cuò)誤率的-10log10值——1（255代表質(zhì)量不可用）
比對的具體方式/CIGAR——3S147M
mate比對到的染色體——=（沒有比對上或沒有mate為*，完全比對為=）
mate比對的第一個(gè)位置——2596（0表示該列不可用）
文庫插入片段的長度——2472（無mate則為0）
read序列——......
測序質(zhì)量——......
AS:i:-3 ZS:i:-3 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:147

從網(wǎng)上查到的bam過濾步驟：

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

看上去應(yīng)該這樣過濾：

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

但是實(shí)際上hisat2比對的結(jié)果MAPQ沒有實(shí)際意義，不表示-10log10比對錯(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)告多重比對是可以指定的，由-k參數(shù)控制。

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

我在網(wǎng)上查到，hisat2默認(rèn)會(huì)look for multiple alignments, report best, with MAPQ，但是在實(shí)際結(jié)果中似乎并不是這樣。

-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ù)輸出了正向比對10個(gè)和反向互補(bǔ)比對10個(gè)，正向和反向各有一個(gè)不包含256的flag（例如1個(gè)83，1個(gè)163，9個(gè)339，9個(gè)419），即只有一個(gè)是首要比對位置，其它都是次要的，它們的MAPQ都是1，都包含tag NH:i:10（唯一比對MAPQ都是60，包含NH:i:1）。如果設(shè)置-k 1，這些多重比對就只能輸出一個(gè)結(jié)果了，MAPQ還是1（也有可能是它正向比對和反向比對各有一個(gè)結(jié)果），但是NH:i:1，且在統(tǒng)計(jì)報(bào)告里會(huì)被認(rèn)為aligned concordantly exactly 1 time。

因此調(diào)整過濾方案為：

$ samtools view -h -@ 2 -b -q 1 -F 256 example.bam > example.filtered.bam

過濾后再去重復(fù)，使用samtools需要四步：

$ samtools sort -n xxx.sort.bam -o xxx.sortname.bam
$ samtools fixmate -m xxx.sortname.bam xxx.fixmate.bam
$ samtools sort xxx.fixmate.bam -o xxx.sortposition.bam
$ samtools markdup -r xxx.sortposition.bam xxx.markdup.bam 
 #加上-r就會(huì)直接去掉重復(fù)序列

也許可以用管道把它們合并起來，但是我操作起來會(huì)出錯(cuò)，不知道為什么。上述四步其實(shí)是可行的，不幸的是由于我在自己筆記本電腦上跑這些程序，內(nèi)存占用直接100%，所以只能換軟件了。

我又嘗試了picard MarkDuplicates，這次更加不幸，CPU占用100%了。

于是我換成了sambamba markdup。

$ sambamba markdup  -r -p -t 2  --tmpdir=./ example.filtered.bam  example.markdup.bam

直接跑這個(gè)程序會(huì)報(bào)錯(cuò)顯示“Too many open files”，解決方案是先使用命令ulimit -n 8000。

過濾并去重復(fù)后的數(shù)據(jù)量會(huì)減少很多，接下來就可以進(jìn)行peak calling了。