一般来说，我们需要对比对软件产生的SAM/BAM文件进行处理，比如SAM/BAM格式的转换、排序、建立索引或者合并。处理这些任务最主流的软件包有两种：SAMtools和它的Java执行版本Picard [15]。Picard包含更多的工具，在处理数据时也比SAMtools更严格。本节我们将介绍常用的SAMtools工具。

1、 将SAM文件转换为BAM文件。

将文件以BAM格式储存可以节省空间，而且许多下游的分析使用的是BAM格式的文件而不是SAM格式的文件。此处，我们指定输入文件为SAM格式（参数 -S），输出文件是BAM格式（参数 -b），同时输出文件名应当为alignments.bam（参数 -o）。

samtools view -bS -o alignments.bam input.sam

2、将BAM格式文件转换为SAM格式文件，同时包含header信息（-h）。

header行以“@”开始，包含参考序列名和长度（@SQ）、创建这个文件所使用的软件（@PG）以及这个文件是否被排序过以及是如何被排序的（@HD）。

samtools view -h -o alignments.sam input.bam

3、 仅获得header信息（-H）。

samtools view -H alignments.bam

4、使用染色体坐标或者read名对比对文件进行排序（-n）。

坐标排序对于使用基因组浏览器进行分析以及一些分析工具来说是必须的，read名排序对于很多表达定量工具来说是必须的。

samtools sort alignments.bam alignments.sorted

samtools sort -n alignments.bam alignments.namesorted

5、 SAMtools可以使用Unix系统的管道可以将命令行结合起来，来避免产生大的中间文件。

下面这条命令可以将Bowtie2的SAM输出转换为BAM，之后通过染色体坐标进行排序，产生一个文件alignments.sorted.bam：

bowtie2 -q - - phred64 -p 4 -x GRCh37.74 -U reads1.fq | samtools view -bS - | samtools sort - alignments.sorted

6、为排序的BAM文件建立索引。

建立索引可以使我们快速提取特定序列、特定位置的比对结果，对于基因组浏览软件或者一些下游分析工具来说是必须的。下面这条命令可以产生索引文件alignments.sorted.bam.bai

samtools index alignments.sorted.bam

7、 列出有多少reads回帖到每条染色体。这条命令需要在预先为比对文件建立索引条件下才能运行。

samtools idxstats alignments.sorted.bam

8、 基于比对质量筛选比对子集。下面这条命令行会筛选出比对质量大于30的比对子集：

samtools view -b -q 30 -o alignments_MQmin30.bam alignments.bam

9、 基于SAM 标记字段的值筛选比对子集。命令行中的-F选项会将带有指定标记值的reads过滤掉（数字4代表没有回帖的reads），-f选项会将带有指定标记值的reads保留下来（数字2意味着read以合适的配对被回帖到基因组）。更加详细的标记值的介绍，可以看SAM格式的文献 [2]。

samtools view -b -F 4 -o alignments.mapped_only.bam alignments.bam

samtools view -b -f 2 -o properly_paired_reads.bam alignments.bam

10、 基于标记字段获取Read回帖统计信息

samtools flagstat alignment.bam

这份报告包含回帖reads数和正确配对reads的信息，以及有多少配对的reads回帖到不同的染色体上：

52841623 + 0 in total (QC-passed reads + QC-failed reads)

0 + 0 duplicates

52841623 + 0 mapped (100.00%: -nan%)

52841623 + 0 paired in sequencing

28919461 + 0 read1

23922162 + 0 read2

42664064 + 0 properly paired (80.74%: -nan%)

44904884 + 0 with itself and mate mapped

7936739 + 0 singletons (15.02%: -nan%)

999152 + 0 with mate mapped to a different chr

357082 + 0 with mate mapped to a different chr (mapQ>=5)

比对的统计学数据也可以用RseQC包 [16]得到。RseQC包含一些检查比对矩阵的Pythons脚本，比如有多少read被比对、其中有多少比对到基因组上的唯一位置、内部距离分布是多少以及多少对reads回帖到基因组上相同的位置。配对的reads回帖到基因组上同样的位置可以告诉我们reads来源于相同的片段，或许是PCR过度扩增导致的。基本的比对统计信息可以用bam_stat.py工具得到：

python bam_stat.py -i accepted_hits.bam

它可以产生下表，如果reads的比对质量大于30则认为该reads是唯一的（通过选项-q 可以改变这个阈值）

#=============================================

#All numbers are READ count

#=============================================

Totalrecords: 52841623

QCfailed: 0

Optical/PCR duplicate: 0

Non primary hits 3098468

Unmapped reads: 0

mapq < mapq_cut (non-unique): 1774335

mapq >= mapq_cut (unique): 47968820

Read-1: 26128297

Read-2: 21840523

Reads map to‘+’: 24085239

Reads map to‘-’: 23883581

Non-splice reads: 35970095

Splice reads： 11998725

Reads mapped in proper pairs: 39702036

Proper-paired reads map to different chrom: 0

Chipster中的SAM/BAM操作工具和比对统计工具

Chipster在工具分类中有很多基于SAMtools的工具。它们可以实现SAM和BAM之间的互相转换、排序、建立索引、亚集、合并BAM文件、计算每条染色体上的配对子数量以及配对子的总数量以及从配对子中得到一条一致序列。一些工具需要BAM索引文件。同时选择BAM文件和索引文件，在参数面板上检查这些文件放到了正确的位置上。

在基因组中对reads进行可视化有很多用途，也是强烈推荐的。用户可以对新的转录本的结构进行可视化、判断新的连接，检查不同外显子的覆盖度以及是否有重复的reads、点插入/删除和SNPs的峰出现等。重要的是，用户可以将数据与参考注释比较。

一些基因组浏览器可以将高通量测序数据可视化，包括Integrative Genomics View IGV [17]，JBrowse [18]，Tablet [19]以及UCSC [20]和Chipster基因组浏览器。想了解这部分内容，我们推荐读者阅读Briefings in Bioinformatics中关于第二代测序可视化的特别讨论 [21]，这部分内容包含了一些详细介绍基因组浏览器的文章。当然IGV软件也很强大。

我们在这里简短介绍一下Chipster这个基因组浏览器。Chipster基于Ensembl注释对数据进行可视化，支持包括BAM、BED、GTF、VCF和tsv在内的几种文件格式。用户可以将数据放大到核苷酸水平，高亮数据与参考序列的差异，同时观察自动计算的覆盖度信息（全或链特异）。使用BED文件也可以对这种分数进行可视化。重要的是，不同的数据可以同时可视化。比如，用户可以并排看到RNA-seq数据和微阵列测量的拷贝数异常。因为Chipster基因组浏览器整合到完善的分析环境中，用户不需要导出和导入数据到一个外源的应用中。当然，如果仅仅是为了可视化的目的，用户可以导入BAM文件到Chipster中。这种情况下，导入的文件会被自动排序并建立索引。

使用Chipster在基因组背景下对reads进行可视化

我们用TopHat2的结果文件accepted_hits.bam和deletions.bed进行可视化作为例子。

1、可以使用BED文件进行导航辅助，所以先分开它到另一个窗口：双击文件打开它到一个单独的电子表格并点击“Detach”。

2、在可视化面板上选择BAM和BED文件以及可视化方法“Genome Browser”，同时将可视化面板最小化来换取更大的观察区。

3、 从“Genome”下拉菜单中选择hg19，点击“Go”。用户可以使用鼠标滚轮进行放大和缩小以及如果需要可以改变覆盖度的尺度。

4、 使用分离的BED文件来有效观察删除突变：点击一个删除突变开始的坐标（第一列），浏览器会移动到那个位置。用户也可以通过点击第四列抬头使用相应分数对BED文件进行排序（带有删除突变的reads数）。

比对文件可以用不同的软件操作，比如SAMtools和Picard，这些比对软件可以让我们有效的取回map到特定区域的reads或唯一比对的reads。诸如RseQC之类的工具可以提供关于比对reads的重要的质量信息。一些基因组浏览器可以在基因组背景下可视化比对结果。对数据的可视化是高度推荐的，因为没有什么比人眼更可以识别数据的有趣模式了。