bam文件过滤和去重复

2022-02-19 本文已影响0人果蝇饲养员的生信笔记

最近处理的几组数据使我感觉到生信软件的教程一定要看最新的或官方的，因为现在的技术和方法发展太快了，软件更新换代也非常快，每一个版本都可能有细微的区别，网上找到的几年前的教程可能在某些方面已经不适用于软件当前版本了。

在做一组MeRIP-seq数据的分析时，遇到了bam文件过滤和去重复的问题，并且再次看了看hisat2和samtools的一些选项和参数。

先是从网上查到的关于PCR重复的问题

RNA-seq一般不去重复（起始量高，扩增数少；某些RNA含量很高）

ChIP-seq一般去重复（起始量低）

call SNP一般去重复（PCR扩增错误影响SNP识别位点置信度）

PCR去重工具首选Picard

万事无绝对，还需参考起始量和PCR扩增数判断是否去重复。reads mapping覆盖均匀度可以判断是否需要去重复。

根源上解决去重复问题：起始量高，循环数少，reads能长不短，能双端不单端

为什么会出现重复

在illumina测序中，通常有两种类型的重复，分别是光学重复和PCR重复。

建库中有一步是PCR扩增加了接头的DNA片段。理想情况下，对打碎的基因组DNA，每个DNA片段仅测到一次。但这一步扩增了6个cycle，那么每个DNA片段有了64份拷贝。将扩增后所有产物“洒”到flowcell，来自一个DNA片段的两个拷贝，可能会锚定在两个bead上，经过测序得到的这两条read，就是PCR duplication。

光学重复，是由于illumina测序时照相机错误的的将一个簇识别为两个簇（多个），这就会导致产生完全一样的reads。该重复可以利用tail的坐标进行去除。

PCR本身就是为了产生重复序列的。理论上来讲，不同的序列在进行PCR扩增时，扩增的倍数应该是相同的。但是由于聚合酶的偏好性，PCR扩增次数过多的情况下，会导致一些序列持续扩增，而另一些序列扩增到一定程度后便不再进行，也就是我们常说的PCR偏好性。

设一个基因组有A、B两个片段，PCR后得到无论多少条reads，比如n・A+m・B条，在数据分析的时候，理想情况都只保留1条A和1条B（unique reads）用于组装，而去掉(n-1)条A和(m-1)条B。共有(n-1)条A和(m-1)条B被当成duplicated reads看待，尽管它们是正常PCR的正常产物。目前的算法其实是一个简化的处理方案，把所有重复的reads都去掉了，留下完全不重复的reads。算法没有能力区分“假重复”（人为造成的重复序列方面的bias）和“真重复”（天然存在的重复序列）。

是否需要去重

若起始量很低，PCR扩增次数很多，那么则需要去除PCR重复。

RNA-seq由于其建库起始量一般都很高所以不需要去除重复，而且RNA-seq数据中经常会出现某些基因的表达量十分高，这就导致这些基因打断的reads的比对情况有很大概率是一致的。

而ChIP-seq中，由于起始量不高，且没有那种富集程度很高的位点，因此通常需要考虑去除PCR重复。

至于call SNP，起始量一般都高（因为要保证测序深度），此外由于PCR扩增会导致一些序列复制错误，这将严重影响SNP位点识别的置信度。因此一般需要去重复。

也可以：利用reads mapping的均匀程度判断是否具有重复。若富集位点周围的reads均匀覆盖，那么没有重复；若富集位点周围覆盖度不均匀，某些区域猛然升高，那么很有可能需要进行PCR去重复。

起始量很多时，不需要去重复。扩增数很少时，15个cycle以内，不需要去重复。双端测序由于其两个reads的位置矫正，有助于去除PCR重复。reads长度越长，越容易识别真正的PCR重复。

去重工具

samtools

如果多个reads具有相同的比对位置时，samtools rmdup将它们标记为duplicates，然后直接将识别出来的重复reads去掉，通常只保留质量最高的一条。

该方法对于以下两种情况，有很好的去除效果：一些reads由于测序错误导致其不完全相同；比对错误导致不同的序列比对到相同的位置（可能性不大）。

该方法的缺点：由于samtools去重只考虑reads比对上的起始终止位置，不考虑比对情况，这种去重有时会导致测序信息的丢失。

双端测序数据用samtools rmdup效果很差，很多人建议用picard工具的MarkDuplicates功能。samtools的rmdup是直接将这些重复序列从比对BAM文件中删除掉，而Picard的MarkDuplicates默认情况则只是在BAM的FLAG信息中标记出来，而不是删除，因此这些重复序列依然会被留在文件中，只是我们可以在变异检测的时候识别到它们，并进行忽略。

目前认为，samtools rmdup已经过时了，应该使用samtools markdup代替。samtools markdup与picard MarkDuplicates采用类似的策略。

Picard

该工具的MarkDuplicates方法也可以识别duplicates。但是与samtools不同的是，该工具仅仅是对duplicates做一个标记，只在需要的时候对reads进行去重。

它不仅考虑reads的比对位置，还会考虑其中的插入错配等情况（即会利用sam/bam文件中的CIGAR值），甚至reads的tail、lane以及flowcell。Picard主要考虑reads的5'端的比对位置，以及每个reads比对上的方向。

因此我们可以从一定程度上认为，5' 端的位置、方向、以及碱基比对情况相同，Picard就将这些reads中碱基比对值Q>15的看作是best pair而其他的reads则当作是duplicate reads。甚至当reads的长度不同时，Picard依然利用上述原理进行去重。

对Picard来说，reads的5' 端信息更为重要。若duplicates是PCR重复，那么它们的序列不一定完全相同。但是由于PCR扩增时，酶的前进方向是5'->3'方向，PCR重复序列中5' 端的部分相似的可能性更高。

sambamba

也可以使用sambamba操作bam文件和去除重复，据说该命令运行比picard MarkDuplicates快30倍。

过滤bam文件和去重复

SAM（sequence Alignment/mapping）数据格式是目前高通量测序中存放比对数据的标准格式。BAM是SAM的二进制格式。

bam文件优点：bam文件为二进制文件，占用的磁盘空间比sam文本文件小；利用bam二进制文件的运算速度快。

举个栗子：
序列名称——E100032134L1C001R0142515642
flag——419
参考序列名称/比对上的染色体——2L（无法比对为*）
比对的第一个位置——277（没有比对上为0）
比对质量/MAPQ/比对错误率的-10log10值——1（255代表质量不可用）
比对的具体方式/CIGAR——3S147M
mate比对到的染色体——=（没有比对上或没有mate为*，完全比对为=）
mate比对的第一个位置——2596（0表示该列不可用）
文库插入片段的长度——2472（无mate则为0）
read序列——......
测序质量——......
AS:i:-3 ZS:i:-3 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:147

从网上查到的bam过滤步骤：

1 去除低质量比对（MAPQ<30）
2 去除多重比对（一条read比对到基因组的多个位置）
3 去除PCR重复（不同reads比对到基因组的同一位置）
4 去除比对到黑名单区域的序列

看上去应该这样过滤：

$ samtools view -h -@ 2 -b -q 30 -F 4 -F 256 example.bam > example.filtered.bam  
#去掉比对质量<30，没有比对到参考序列的，多重比对的

但是实际上hisat2比对的结果MAPQ没有实际意义，不表示-10log10比对错误率，它只有0, 1, 60三个值。

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如何报告多重比对是可以指定的，由-k参数控制。

-k <int> report up to <int> alns per read; MAPQ not meaningful

-a/--all report all alignments; very slow, MAPQ not meaningful

我在网上查到，hisat2默认会look for multiple alignments, report best, with MAPQ，但是在实际结果中似乎并不是这样。

-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.

事实上hisat2默认参数输出了正向比对10个和反向互补比对10个，正向和反向各有一个不包含256的flag（例如1个83，1个163，9个339，9个419），即只有一个是首要比对位置，其它都是次要的，它们的MAPQ都是1，都包含tag NH:i:10（唯一比对MAPQ都是60，包含NH:i:1）。如果设置-k 1，这些多重比对就只能输出一个结果了，MAPQ还是1（也有可能是它正向比对和反向比对各有一个结果），但是NH:i:1，且在统计报告里会被认为aligned concordantly exactly 1 time。

因此调整过滤方案为：

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

过滤后再去重复，使用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就会直接去掉重复序列

也许可以用管道把它们合并起来，但是我操作起来会出错，不知道为什么。上述四步其实是可行的，不幸的是由于我在自己笔记本电脑上跑这些程序，内存占用直接100%，所以只能换软件了。

我又尝试了picard MarkDuplicates，这次更加不幸，CPU占用100%了。

于是我换成了sambamba markdup。

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

直接跑这个程序会报错显示“Too many open files”，解决方案是先使用命令ulimit -n 8000。

过滤并去重复后的数据量会减少很多，接下来就可以进行peak calling了。