软件6 —— sambamba和picard

一、 基本介绍

(1) 关于PCR重复的问题

• RNA-seq一般不去重复（起始量高，扩增数少；某些RNA表达量很高，导致该基因的reads的比对很可能出现一致的情况）
• ChIP-seq一般去重复（起始量低，扩增数多，PCR偏好性影响定量分析）
• call SNP一般去重复（PCR扩增数多，出现扩增错误，影响SNP识别置信度）
• PCR去重工具首选Picard
• 万事无绝对，还需参考起始量和PCR扩增数判断是否去重复。reads mapping覆盖均匀度可以判断是否需要去重复。
• 根源上解决去重复问题：起始量高，循环数少，reads能长不短，能双端不单端

(2) 简化的处理方案

设一个基因组有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）和“真重复”（天然存在的重复序列）。

二、 背景知识

测序所得到的reads是由于超声波或者酶切断裂得到的，因此这些reads比对到基因组上的位置是完全随机的。那么两个reads比对到相同位置的概率是非常低的。如果两个reads比对情况相同或者极其相似，则很有可能是由于PCR重复所导致的。

在illumina测序中，通常有两种类型的重复，分别是光学重复和PCR重复。PCR重复就是在建库过程中PCR扩增导致的重复。理想情况下，对打碎的基因组DNA，每个DNA片段仅测到一次。但如果这一步扩增了6个cycle，那么每个DNA片段有了64份拷贝。将扩增后所有产物“洒”到flowcell，来自一个DNA片段的两个拷贝，可能会锚定在两个bead上，经过测序得到的这两条read，就是PCR duplication。光学重复是由于illumina测序时照相机错误地将一个簇识别为两个或多个簇，这就会导致产生完全一样的reads。光学重复可以利用tile的坐标进行去除。

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

也可以从根源上减少这种重复带来的影响。起始量很多时，不需要去重复。扩增数很少时，15个cycle以内，不需要去重复。双端测序由于其两个reads的位置矫正，有助于去除PCR重复。reads长度越长，越容易识别真正的PCR重复。

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

※ PCR本身就是为了产生重复序列的。理论上来讲，不同的序列在进行PCR扩增时，扩增的倍数应该是相同的。但是由于聚合酶的偏好性，PCR扩增次数过多的情况下，会导致一些序列持续扩增，而另一些序列扩增到一定程度后便不再进行，也就是我们常说的PCR偏好性。这种情况对于定量分析（如ChIP-seq），会造成严重的影响。（ChIP-seq中，由于起始量不高，且没有那种富集程度很高的位点，因此通常需要考虑去除PCR重复。）

※ 此外，PCR扩增循环数过多，会出现一些扩增偏差，进而影响一些突变识别（比如call SNP）的置信度。PCR重复也会增大变异检测结果的假阴和假阳率。（至于call SNP，因为要保证测序深度，起始量一般都高，但由于PCR扩增会导致一些序列复制错误，这将严重影响SNP位点识别的置信度，因此一般需要去重复。）

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

※ 去除MeRIP-Seq中的重复片段read应该对低表达基因有益，但对高表达基因却不友好。我们一般更关心假阳性错误而不是假阴性错误（宁愿错过真正的peak也不将希望将假peak当作真的），所以对duplicate read进行过滤，以避免由于PCR 产物引起的假阳性peak，但可能会错过高表达转录本上的peak。

三、 去重工具

(1) 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采用类似的策略。

(2) Picard

• 该工具的MarkDuplicates方法也可以识别duplicates。但是与samtools rmdup不同的是，该工具仅仅是对duplicates做一个标记，只在需要的时候对reads进行去重。
• 它不仅考虑reads的比对位置，还会考虑其中的插入错配等情况（即会利用sam/bam文件中的CIGAR值），甚至reads的tile、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' 端的部分相似的可能性更高。

(3) sambamba

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

四、 使用方法

(1) bam过滤步骤

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

(2) 查看bam文件

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

第1列：序列名称—— E100032134L1C001R0030000131
第2列：flag数字之和—— 163
第3列：参考序列名称/比对上的染色体—— 2L（若无法比对为*）
第4列：比对的第一个位置——900488（没有比对上为0）
第5列：比对质量/MAPQ/比对错误率的-10log10值—— 60（255代表质量不可用）
第6列：比对的具体方式/CIGAR—— 2S145M3S
第7列：mate比对到的染色体—— =（没有比对上或没有mate为*，完全比对为=）
第8列：mate比对的第一个位置—— 900488（0表示该列不可用）
第9列：文库插入片段/fragment的长度—— -152（无mate则为0）
第10列：read序列—— CGGCGCCT……
第11列：测序质量—— DDDDDDDD……（D对应ASCII值68，即35+33，Q35）
第12列：可选区域—— 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没有比对上；8 = mate没有比对上；256 = secondary比对；1024 = read是PCR或光学copy产生。

第5列：
实际上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> 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.

事实上hisat2默认参数输出了正向比对10个和反向互补比对10个，正向和反向各有一个不包含256的flag（例如1个83，1个163，9个339，9个419），即只有一个是首要比对位置，其它都是次要的，它们的MAPQ都是1，都包含tag NH:i:10（唯一比对MAPQ都是60，包含NH:i:1。NH: The number of mapped locations for the read or the pair.）。如果设置-k 1，这些多重比对就只能输出一个结果了，MAPQ还是1（也有可能是它正向比对和反向比对各有一个结果），但是NH:i:1，且在统计报告里会被认为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进行过滤

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

samtools view  -@ 4  -h  -b  -q 1  -F 256  example.bam  >  example.filtered.bam
#实际对hisat2产生的BAM文件可以用-q 1  -F 256过滤
# samtools的-q 1 比 sambamba "not unmapped"过滤掉的序列更多，包括没有比对上的reads、错误过多的多重比对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  &

-@：线程数。
-h：设定输出sam文件时带header信息，默认输出sam不带header信息。
-b：输出为bam格式，默认输出SAM格式。
-f：需要的flag。
-F：过滤掉的flag。
-q：have mapping quality >= INT

(4) 用sambamba view进行过滤

# sambamba view  -t 4  -h  -f bam  -F "[XS] == null and not unmapped and not duplicate"  example.bam  >  example.filtered.bam   #未经实践的步骤

-t：线程数
-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。输出文件的格式，默认sam
-F：--filter=FILTER。set custom filter for alignments（[XS]是bowtie2多重比对才会有的得分）
--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去重复

按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就会直接去掉重复序列

(6) 用picard MarkDuplicates去重复

# 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去重复

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  &
# 如果出现Too many open files报错，需要通过使用ulimit -n 8000或添加--overflow-list-size=600000来解决

-t：线程数
-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得到第二次机会满足它们的pairs的溢出列表的大小（默认为200000 reads）；增加大小会减少创建的临时文件的数量。

(8) 建立索引

# for i in `ls ./06_markdup/*.markdup.bam` ; do (samtools index $i) ; done  &
# sambamba去重复时会自动建立索引

(9) 进行统计

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  &

五、 举个例子

(1) 用samtools view进行过滤

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  &
#去除没有比对上的reads、错误过多的多重比对reads（-q 1）；去除多重比对次要位置（secondary，-F 256）

(2) 用sambamba markdup去重复

ulimit -n 8000
#去除PCR重复
# 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去重复时会自动建立索引

(4) 进行统计

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  &