上一篇记录了mapping这一步的软件选择，也讲到了对于sam文件如何考虑MAPQ的过滤，这篇我主要想记录一下bam文件在进行call variation之前的处理。
包括MAPQ的筛选等都是这个处理的一部分。
既然要处理bam文件，不得不提bam文件的格式。
处理生物信息的数据时会发现，文件有各种各样的格式，眼花缭乱，这也是我一开始接触课题接触分析流程时的感受。但是多和这些文件打交道以后，会发现，大多数的不同格式的文件其本质都是文本文件，为了某种特殊的处理或记录需要，按一定的规则记录信息。包括之前接触的fasta、fastq文件，也包括sam文件，以后后面步骤会接触到的vcf文件等。
bam文件是sam文件的二进制格式，这里贴一张图，来说明一下为什么要转换成bam文件来处理。
可以看到，不经过处理的sam文件是非常大的，非常不利于存储和处理，而转换格式后的bam文件小很多，所以很多处理软件也是针对bam文件进行开发的。
包括bwa，bowtie2的输出文件，都是sam文件，但是sam文件具体是什么样子的，我这里不会展开来讲，一是因为官网说明书已经说的很清楚了，二是因为你随手一个百度、谷歌（如果你翻得动墙）就有很多很多的人发帖介绍，而内容大体都是类似的，我自认为也没有什么能比他们讲得更好的。但是了解这个文件的内容确实是非常重要的。
以
@开头的行记录了header信息，之后的行记录了每个reads的mapping信息，一般对于这样的sam文件我们要做的处理大体就是先转换为bam文件，进行sort，再进行merge，最后mark duplicates，并建立索引。接下来我将一步步介绍怎么去处理，以及为什么我是这么处理的。在开始之间先使用samtools对参考序列进行索引建立，作用是用于samtools软件快速识别，这一步也是一劳永逸的，只要不把索引文件删除。
$ samtools faidx re.fa
#得到索引文件：re.fa.fai
首先是把sam文件转换成bam文件，我们通过samtools view来实现，代码如下：
$ samtools view -S in.sam -b > out.bam
$ samtools view -b in.bam -S > out.sam
但实际上，在真正的操作中，我们是不会保留sam文件的，甚至是不会产生sam文件的，因此，我们通常这么来写命令。
$ bwa mem -t 12 -M -R "@RG\tID:W2018001\tPL:ILLUMINA\tLB:W2018001\tSM:W2018001" /your/path/of/reference/ucsc.hg19.fasta 1.fq.gz 2.fq.gz | samtools view -q 1 -Shb - > W2018001.bam
# 关于“|”之前的我就不解释了，看我上一篇简书
#-q 1 ：筛选MAPQ用的，意为保留MAPQ >= 1的记录，上一篇简书中讨论过关于MAPQ的由来和分布，这里用“1”也是保险起见
#-h ：表示保留header信息
#-S，-b：表示格式，S指的是sam格式，b指的是bam格式
这样操作的好处是直接可以得到bam文件，不必占用大量的空间存储sam文件。如果你实在是需要sam文件也可以转换回去，但如果你只是想查看一下，也可以通过以下方式实现。还是很方便的。
$ samtools view -h xxx.bam | less
在测序的时候序列是随机打断的，所以reads也是随机测序记录的，进行比对的时候，产生的结果自然也是乱序的，为了后续分析的便利，将bam文件进行排序。事实上，后续很多分析都建立在已经排完序的前提下。关于排序可以通过以下命令完成。
$ samtools sort -@ 6 -m 4G -O bam -o sorted.bam W2018001.bam
# @：线程数
# m：每个线程分配的最大内存
# O：输出文件格式
# o：输出文件的名字
# 输入文件放在最后
接下来要做的是merge，这个不是所有的样本都需要做的！！！
之前有提到，WGS的数据比较大，对于一个样本可能有多对文件，一般有两个思路，一个是先对原始的fastq文件进行merge，一个是对后面的bam文件进行merge。那么我选用的是后者。实现脚本如下：
$ samtools merge -@ 6 -c sorted.merged.bam *.sorted.bam
#@：线程数
#c：当输入bam文件的header一样时，最终输出一个header
#当然也可以用-h file，定义你想要的header，file里的内容就是sam文件header的格式
#第一个bam是输出的文件名，后面排列输入的文件名，我这里用了通配符‘*’，要保证该目录下‘.sorted.bam’结尾的都是你要输入的文件
#当然也可以用-b file，file文件里罗列要merge的文件名，一行一个文件名
下一步就是remove duplicates，为什么要进行这一步呢？先来说一下测序，我们都知道人的基因组是很庞大的（约30亿个碱基对），在测序的时候，先会把基因组的DNA序列通过超声震荡随机打断成150bp的片段，那么从概率上来讲，出现同样的片段（开始和结束位置都一样）的几率是极小的。但是由于PCR对某些位置有偏好的扩增，会导致一些一模一样的reads存在。这些reads其实是一个片段的扩增导致的，多出来的reads，对该区段突变的判断是没什么贡献的，如果不加处理，反而会大大增加那个位置的深度，导致某些结果的不准确。
如下图所示，箭头所指的reads，就是duplicates，我们一般采取的策略是标记或者去除，这样的话，下一步call variation的时候会不考虑这些reads。
关于这一步，有很多软件可以实现，比较多用的是picard和samtools rmdup。我这里用的是GATK包里集成的picard的MarkDuplicates。关于picard和samtools rmdup的效果其实大家可以自己试一下，我很早之前做过的试验是，samtools rmdup速度很快，但是去除的效果稍差。大家用的最多的也是picard。
$ java -jar /you/path/of/gatk/gatk-package-4.0.10.1-local.jar \
MarkDuplicates \
-INPUT sorted.merged.bam \
-OUTPUT sorted.merged.markdup.bam \
-METRICS_FILE markdup_metrics.txt
#也可以加上“REMOVE_DUPLICATES=true”来去除掉这些duplicates，不然就只是标记
到了这一步基本上就处理的差不多了，可以进行call variation了，但是这里还有一步建立索引，这十分的重要！！！！
必须对bam文件进行排序后，才能进行index，否则会报错。建立索引后将产生后缀为.bai的文件，用于快速的随机处理。很多情况下需要有bai文件的存在，特别是显示序列比对情况下。建立索引很简单。
$ samtools index sorted.merged.markdup.bam
到了这一步就基本上完事了，可以通过IGV或tview来查看情况，这都需要建立索引，且索引文件和bam文件在同一个目录下。
$ samtools tview sorted.merged.markdup.bam re.fa
#可以用-p chr:pos（加在tview和sorted.merged.markdup.bam之间）指定要查看的位置
#也可以进去后用敲‘g’输入`chr10:10,000,000' 或 `=10,000,000'查看指定的位置，敲'?'了解更多说明，q退出
下图就是效果了，用“？”，打开左边的帮助界面，其中圆点表示正链比对，逗号表示负链比对，星号表示插入。不同的颜色代表不同的含义，具体怎么调看帮助框。要是觉得不好看的话可以用IGV查看。IGV的效果就是上图duplicates的效果。
同时对于得到的bam文件也可以进行一下查看，对于bam的统计软件就更多了，我这里网罗了一些帖子上的推荐以及我自己日常使用的软件，感兴趣的可以自己去下载来使用一下。
samtools
这是个强大的软件，也自带一些统计工具，上篇简书其实我们就用到了，主要是四个：idxstats，depth，stats，flagstat
$ samtools depth sorted.merged.markdup.bam
示例结果
#chr pos depth
chr1 1 5
结果可以得到染色体名称、位点位置、覆盖深度
-a:输出所有位点，包括零深度的位点
-a -a --aa:完全输出所有位点，包括未使用到的参考序列
-b FILE：计算BED文件中指定位置或区域的深度
-f FILE：指定bam文件
-l INT:忽略小于此INT值的reads
-q INT：只计算碱基质量大于此值的reads
-Q INT：只计算比对质量大于此值的reads
-r CHR:FROM-END：只计算指定区域的reads
$ samtools idxstats sorted.merged.markdup.bam #需要预先进行sort和index
示例结果
#ref sequence_length mapped_reads unmapped_reads
chr1 195471971 6112404 0
结果可看出，分别统计染色体名称、序列长度、mapped数目，以及unmapped数目
$ samtools flagstat sorted.merged.markdup.bam
示例结果：
20607872 + 0 in total (QC-passed reads + QC-failed reads) #总共的reads数
0 + 0 duplicates #重复reads的数目
19372694 + 0 mapped (94.01%:-nan%) #总体上reads的数目以及匹配率；可能有有小偏差
20607872 + 0 paired in sequencing #paired reads的数目
10301155 + 0 read1 #reads1的数目
10306717 + 0 read2 #reads2的数目
11228982 + 0 properly paired (54.49%:-nan%) #完美匹配的reads数：比对到同一条参考序列，并且两条reads之间的距离符合设置的阈值
18965125 + 0 with itself and mate mapped#两条都比对到参考序列上的reads数，但是并不一定比对到同一条染色体上
407569 + 0 singletons (1.98%:-nan%)#只有一条比对到参考序列上的reads数，和上一个相加，则是总的匹配上的reads数。
3059705 + 0 with mate mapped to a different chr#两条分别比对到两条不同的染色体的reads数
1712129 + 0 with mate mapped to a different chr (mapQ>=5)#两条分别比对到不同染色体的且比对质量值大于5的数量
#说明：
1.bwa的mem比对方法生成的bam文件保留sencondly比对的结果。所以使用flagstat给出的结果会有偏差。
2.每一项统计数据都由两部分组成，分别是QC pass和QC failed，表示通过QC的reads数据量和未通过QC的reads数量。以“PASS + FAILED”格式显示
$ samtools stats sorted.merged.markdup.bam
#对bam文件进行详细的统计，也可只统计某一染色体的某一区域chr：from-to
结果包括：
Summary Numbers,raw total sequences,filtered sequences, reads mapped, reads mapped and paired,reads properly paired等信息
Fragment Qualitites #根据cycle统计每个位点上的碱基质量分布
Coverage distribution #深度为1，2，3，，，的碱基数目
ACGT content per cycle #ACGT在每个cycle中的比例
Insert sizes #插入长度的统计
Read lengths #read的长度分布
stats出来的结果可以使用plot-bamstats来画图(samtools目录下misc目录中)
$ plot-bamstats -p outdir/ sorted.merged.markdup.bam.stats
下图就是plot-bamstats操作的输出结果了，可以看到有很多的图。效果还是很好的。
更多关于samtools的工具以及上文提到的工具的其余参数请参考官网：http://www.htslib.org/doc/samtools.html
RSeQC
这也是上篇简书中提到过的，所以安装方式和使用就不赘述了，其实里面还有一些其余实用的工具，当然这款软件的最初使用对象是RNA-seq，但并不影响使用，有些工具是通用的，有一点要注意的是，bam_stat.py里的unique mapping的默认阈值是MAPQ>=30，当然可以通过参数来修改，这个在结果的理解上大家要注意。
bedtools
这是一个经常使用也很实用的软件，我经常会用来截取bam然后在igv上看突变的情况，师姐推荐其中的mutlicov进行覆盖度的统计。我粗略的看了一下bedtools的说明书，用于coverage统计的应该还有coverage，genomecov。感兴趣的可以尝试一下。
bedtools:
$ wget https://github.com/arq5x/bedtools2/releases/download/v2.27.1/bedtools-2.27.1.tar.gz
$ tar -zxvf bedtools-2.27.1.tar.gz
$ cd bedtools2
$ make
#脚本在bin/下的bedtools
#Ubuntu用户也可以使用下述命令来下载：
$ sudo apt-get install bedtools
截取bam文件，查看igv可以用以下命令：
$ bedtools intersect -a sorted.merged.markdup.bam -b region.sorted.bed > target_reads.bam && samtools index target_reads.bam
#其中bed文件的格式可以参考：
#染色体号 起始位置 终止位置
chr1 226173187 226173247
#用"\t"分隔
GATK
GATK不仅可以call variation，里面还包含了很多其他用途的工具包，其中有一个工具叫DepthOfCoverage，可以统计depth和coverage，但是在panel中比较适用，因为有bed范围，且比较小。这个工具的速度是比较慢的，要在全基因组上做不太现实。所以我本人也没去使用。
BAMStats
一款比较早的bam统计软件，没用过，但是看说明使用是极其简单了，说一下怎么安装。感兴趣的可以自己试一下，很简单。
$ wget https://jaist.dl.sourceforge.net/project/bamstats/BAMStats-1.25.zip
$ unzip BAMStats-1.25.zip
bamdst
一款简单好用的深度统计软件。
$ git clone https://github.com/shiquan/bamdst.git
$ cd bamdst/
$ make
安装好后，需要准备.bed文件及.bam文件，以软件提供的MT-RNR1.bed和test.bam为例，使用命令如下：
$ bamdst -p MT-RNR1.bed -o ./ test.bam
输出的结果包含7个文件，为：
-coverage.report
-cumu.plot
-insert.plot
-chromosome.report
-region.tsv.gz
-depth.tsv.gz
-uncover.bed
主要看一下coverage.report文件，里面包含了大量信息。
qualimap
这算是压轴了吧，这个软件是我师姐推荐的，安装使用都比较容易，给出的是PDF或HTML的报告，报告中的信息特别丰富，还有一堆的图，所以在我自己的脚本中也会对每个样本使用该软件统计，简述一下安装和使用。
#第一步：下载
$ mkdir qualimap
$ cd qualimap
$ wget https://bitbucket.org/kokonech/qualimap/downloads/qualimap_v2.2.1.zip
$ unzip qualimap_v2.2.1.zip
$ cd qualimap_v2.2.1
#第二步安装相关的软件
#java
#这个应该都有，这里就是贴一下官网的步骤
$ sudo apt-get install openjdk-6-jre
#R
#官网上也有，我贴的是自己以前安装的记录
$ apt-key adv --keyserver keyserver.ubuntu.com --recv-keys E298A3A825C0D65DFD57CBB651716619E084DAB9
$ apt-get update
$ apt-get install r-base
$ apt install r-base-core
#R包安装，两个方法，第一个听说容易报错，我用的是第二个
$ Rscript scripts/installDependencies.r
#或
$ R
> install.packages("optparse")
> source("http://bioconductor.org/biocLite.R")
> biocLite(c("NOISeq", "Repitools", "Rsamtools", "GenomicFeatures", "rtracklayer"))
> q()
使用也简单，主要分为带gtf文件和不带gtf文件，全基因组的话一般不带gtf文件，然后bamqc其实也只是这个软件中的一个工具，其他的工具感兴趣的也可以看看。
$ qualimap bamqc -bam sorted.merged.markdup.bam --java-mem-size=20G -c -nt 16 -outdir bamqc -outfile bamqc.pdf -outformat PDF:HTML
#参数也没有什么特别要注意的，基本上默认的就可以
找了一个示例结果，发现有23页，我这里就不贴了，大家可以自己尝试一下。
最后贴两张图，是我自己写的脚本得到的深度分布，累积曲线以及覆盖率，因为是自己写的，所以比较糙，横纵坐标标题什么的一律没写。
上图可以看到，深度分布还是比较正态的，最多的30X左右，至于左边0X为什么这么高，是因为参考基因组有些位置就是N，还有一些位置就是我的样本没有覆盖到。
上图可以看到，小于10X的数据的不到2%，超过50%的数据都是高于30X的，还是不错的。
上图我按不同的染色体进行统计覆盖率，去掉了其余的一些未知染色体位置的序列上的信息（这个信息具体要了解参考基因组的特性，比如有些序列目前能明确在人身上，却不知道具体在哪个染色体等，这些信息也是包含在参考基因组中的，仔细看参考基因组会发现，除了22条常染色体，2条性染色体，还有很多其他的序列），这个覆盖率并不是我想想的那般全体高于99%，也没公司给的报告描述的那么好，我不知道是不是因为MAPQ的过滤导致了这样的结果，但是总体的覆盖率还可以。
水平有限，要是存在什么错误请评论指出！请大家多多批评指正，相互交流，共同成长，谢谢！！！
