一、原理

原文章：The Sequence alignment/map (SAM) format and SAMtools

摘要

1、SAM格式是一种通用的，用于储存比对后的信息，可以支持来自不同平台的read的比对结果
2、SAM文件在格式上很灵活，易于压缩、可以高效获取以及是千人基因组计划中使用的比对格式
3、SAMtools可以用于处理储存为SAM格式的比对结果文件，可以做indexing、variant caller 、alignment viewer等操作

背景

1、高通量测序技术的诞生衍生了很多的测序方法和原理不同的比对工具，此时需要一种公认的比对格式，可以支持所有的测序方法和比对工具，这样可以作为比对与下游分析沟通的桥梁。
2、SAM格式应运而生，它不仅支持单双端测序，而且支持组合的reads包括color space reads fromAB/SOLiD ，它可以支持10的11次方甚至更多的碱基对的比对结果。
3、同时开发了用于高效便捷处理sam格式的samtools工具

SAM格式

SAM格式简介

1、SAM格式包括头部和比对部两部分，其中头部每行以@开头，比对部则没有@。所有的行都以tab分隔
2、比对部中包括11个必选项和1个可选项，如果必选项的信息没有获取到，则用*或0表示
3、可选项是由键值对组成，格式为TAG:TYPE:VALUE，这边储存了额外的信息
4、例如"RG"这个tag，保存了read的read group信息，若与@RG头部联用，则每条read都必须有origin、sequencing center以及library信息
5、SAM格式中的每一部分信息及tag信息都有具体的描述

SAM的基本格式

SAM的11中必须项

SAM格式的必选项之一----CIGAR简介

1、传统的CIGAR包括三种选项，反应碱基比对的三种结果：M表示匹配或错配，I表示相对参考基因组是插入的，D表示相对于参考基因组是删除的
2、extend CIGAR又添加了四种格式：N表示相对于参考序列，该碱基skipped，S表示soft clipping ，H表示hard clipping ，p表示padding，这些格式支持剪切、clipping、multi——part和padded

extend CIGAR格式

BAM格式简介

1、BAM格式是SAM的二进制格式，其中的信息于SAM相同
2、BAM格式是通过作者开发的BGZF library工具进行压缩的

sort和idexing

1、一个BAM或SAM文件往往是unsort的
2、根据坐标进行sort往往可以简化数据处理流程同时避免将额外的比对结果加载到内存中，同时一个根据位置进行sort的BAM文件可以进一步index
3、作者使用了UCSC binning scheme 和j简单的线性indexl来实现快速随机的对一个染色质特定区域的read的比对结果进行检索，在很多情况下，要在一个区域中检索比对信息，只要一个seek call 就可以了。

SAMtools软件

1、SAMtools是一个可以解析和处理SAM或BAM格式的比对结果的工具包
2、它可以对比对结果进行转化，sort、merge，去除PCR重复，call SNP以及小的indel突变，以及可以用text格式的浏览器查看比对结果
3、本身有两种版本，基于C的和基于JAVA的

二、软件操作

使用说明：http://www.htslib.org/doc/samtools.html
https://www.jianshu.com/p/6b7a442d293f

1、用法汇总

samtools view -bt ref_list.txt -o aln.bam aln.sam.gz

samtools sort -T /tmp/aln.sorted -o aln.sorted.bam aln.bam

samtools index aln.sorted.bam

samtools idxstats aln.sorted.bam

samtools flagstat aln.sorted.bam

samtools stats aln.sorted.bam

samtools bedcov aln.sorted.bam

samtools depth aln.sorted.bam

samtools view aln.sorted.bam chr2:20,100,000-20,200,000

samtools merge out.bam in1.bam in2.bam in3.bam

samtools faidx ref.fasta

samtools fqidx ref.fastq

samtools tview aln.sorted.bam ref.fasta

samtools split merged.bam

samtools quickcheck in1.bam in2.cram

samtools dict -a GRCh38 -s "Homo sapiens" ref.fasta

samtools fixmate in.namesorted.sam out.bam

samtools mpileup -C50 -f ref.fasta -r chr3:1,000-2,000 in1.bam in2.bam

samtools coverage aln.sorted.bam

samtools flags PAIRED,UNMAP,MUNMAP

samtools fastq input.bam > output.fastq

samtools fasta input.bam > output.fasta

samtools addreplacerg -r 'ID:fish' -r 'LB:1334' -r 'SM:alpha' -o output.bam input.bam

samtools collate -o aln.name_collated.bam aln.sorted.bam

samtools depad input.bam

samtools markdup in.algnsorted.bam out.bam

2、常见用法
2.1 samtools view

samtools view [options] in.sam|in.bam|in.cram [region...]#查看和转化SAM/BAM/CRAM文件
#实例
samtools view -b -1 -o selected.bam -U unselected.bam input.bam chr1:1000-10000
#输出比对到一号染色体1001到10000的比对结果，并保存到selected.bam中，其余结果保存到unselected.bam中
#如果options和region不指定，那就会输出比对的所有内容
#常用options
#-b Output in the BAM format#输出为bam格式
#-1 Enable fast BAM compression (implies -b)
#-h Include the header in the output
#-o FILE Output to FILE [stdout].
#-U FILE Write alignments that are not selected by the various filter options to FILE. When this option is used, all alignments (or all alignments intersecting the regions specified) are written to either the output file or this file, but never both.
#region的表示方法：RNAME[:STARTPOS[-ENDPOS]]
#chr1 Output all alignments mapped to the reference sequence named `chr1' (i.e. @SQ SN:chr1)
#chr2:1000000 The region on chr2 beginning at base position 1,000,000 and ending at the end of the chromosome
#chr3:1000-2000 The 1001bp region on chr3 beginning at base position 1,000 and ending at base position 2,000 (including both end positions).
#'*' Output the unmapped reads at the end of the file. (This does not include any unmapped reads placed on a reference sequence alongside their mapped mates.)
#. Output all alignments. (Mostly unnecessary as not specifying a region at all has the same effect.)
#通过man samtools view 或在samtools说明文档中查看具体参数内容
#其它参数的概要内容
#The -b, -C, -1, -u, -h, -H, and -c options change the output format from the default of headerless SAM, and the -o and -U options set the output file name(s).

#The -t and -T options provide additional reference data. One of these two options is required when SAM input does not contain @SQ headers, and the -T option is required whenever writing CRAM output.

#The -L, -M, -r, -R, -d, -D, -s, -q, -l, -m, -f, -F, and -G options filter the alignments that will be included in the output to only those alignments that match certain criteria.

#The -x and -B options modify the data which is contained in each alignment.

#The -X option can be used to allow user to specify customized index file location(s) if the data folder does not contain any index file. See EXAMPLES section for sample of useage.

#Finally, the -@ option can be used to allocate additional threads to be used for compression, and the -? option requests a long help message.

2.2 samtools sort

#对sam/bam/cram文件排序（根据左边的坐标排序）
samtools sort [-l level] [-m maxMem] [-o out.bam] [-O format] [-n] [-t tag] [-T tmpprefix] [-@ threads] [in.sam|in.bam|in.cram]
#实例
samtools sort -@ 15 -I 5 -m 20M -n -o my_sorted.bam -O bam my_unsorted.bam
#常见option
#-@ INT 设定线程，默认为单线程
#-l INT 设定输出文件的压缩水平，0表示不压缩，1表示快速但低水平的压缩，9表示缓慢但高水平的压缩，范围是0或者1-9，如果不指定-I，则使用默认的压缩水平
#-m INT 设定每个线程需要的最大内存，使用K M G作为后缀，为了避免sort时产生很多中间文件，该参数最小值设为1M
#-n 根据read name进行sort而不是根据染色体的坐标
#-o FILE 最后输出的文件，如果为空，就输出到标准输出中
#-O FORMAT 输出文件的格式，SAM/BAM/CRAM中选一个
#默认情况下的排序，是根据reference（@SQheader recorder），然后根据reference中的位置，之后时REVERSE flag.

2.3 samtools index

#为bam文件加索引
samtools index [-bc] [-m INT] aln.sam|aln.bam|aln.cram [out.index]
#索引经过sort的BAM文件
#实例
samtools index -bc my.bam my.index
#最终索引会以my.index的文件名输出
samtools index -bc my.bam
#索引会输出为my.bai和my.csi两个文件
#常见选项
#-@, --threads INT 用于压缩的线程数目，默认是0
#-b 创建BAI索引，当其它索引类型没有设定时，输出该类索引
#-c 创建CSI索引
#-m INT 创建一个CSI索引，不同于上面的
#如果filename给定了，那么索引文件就会命名为out.index，如果没有给定，那么对于xxx.bam文件作为输入的情况，会输出名为xxx.bai或xxx.csi的文件

2.4 samtools depth：计算每个位置或每个区域的read的深度
2.5 samtools split：通过read group 对一个文件进行分割
2.6 samtools mpileup：从一个比对结果中产生plieup的文本格式的文件
2.7 samtools coverage：产生每个染色体上的覆盖度的柱形图和表格
2.8 samtools markdup：在sort好的bam中标记出重复的比对

附：SAM格式的详细说明

文章摘自：http://www.biotrainee.com/thread-2704-1-1.html

1）Sam (Sequence Alignment/Map)

1) SAM 文件产生背景

随着Illumina/Solexa, AB/SOLiD and Roche/454测序技术不断的进步，各种比对工具产生，被用来高效的将reads比对到参考基因组。因为这些比对工具产生不同格式的文件，导致下游分析比较困难，因此一个通用的格式可以提供一个很好的接口用于链接比对与下游分析(组装，变异等，基因分型等)。因此SAM格式应运而生，主要是用来存储测序reads与参考序列比对结果信息的一种文件格式，以TAB为分割符，支持不同平台的短reads及长reads(最长为128Mbp)。

2）格式解读

我们用文献中的例子来详细解释sam格式。

2.1）首先看一个比对事件：

image

ref是参考序列，Read r001/1和 r001/2组成read pair，r003是嵌合体(chimeric read) ,r004表示 split alignment事件
2.2）相应的sam格式是：

image

image

由此我们可以看到，SAM是由两部分组成：分为标头注释信息(header section)和比对结果(alignment section)。标头信息可有可无，都是以@开头，用不同的tag表示不同的信息，主要有:

• @HD，说明符合标准的版本、对比序列的排列顺序（这里为coordinate）
• @SQ，参考序列说明 （SN:ref,LN 是参考序列的长度)
• @PG，使用的比对程序说明（这里没有给出）
• 比对结果部分(alignment section)每一行表示一个片段（segment）的比对信息，包括11个必须的字段（mandatory fields）和一个可选的字段，字段之间用tag分割。必须的字段有11个，顺序固定，根据字段定义，可以为’0‘或者’*‘，这11个字段是：

1）QNAME：比对片段的（template）的编号；
2）FLAG：位标识，template mapping情况的数字表示，每一个数字代表一种比对情况，这里的值是符合情况的数字相加总和；进一步学习可查看https://broadinstitute.github.io/picard/explain-flags.html
3）RNAME：参考序列的编号，如果注释中对SQ-SN进行了定义，这里必须和其保持一致，另外对于没有mapping上的序列；
4）POS：比对上的位置，注意是从1开始计数，没有比对上，此处为0；
5）MAPQ：mappint的质量；
6）CIGAR：简要比对信息表达式（Compact Idiosyncratic Gapped Alignment Report），其以参考序列为基础，使用数字加字母表示比对结果，比如3S6M1P1I4M，前三个碱基被剪切去除了，然后6个比对上了，然后打开了一个缺口，有一个碱基插入，最后是4个比对上了，是按照顺序的；
7）RNEXT：下一个片段比对上的参考序列的编号，没有另外的片段，这里是’‘，同一个片段，用’=‘；
8）PNEXT：下一个片段比对上的位置，如果不可用，此处为0；
9）TLEN：Template的长度，最左边得为正，最右边的为负，中间的不用定义正负，不分区段（single-segment)的比对上，或者不可用时，此处为0；
10）SEQ：序列片段的序列信息，如果不存储此类信息，此处为’‘，注意CIGAR中M/I/S/=/X对应数字的和要等于序列长度；
11）QUAL：序列的质量信息，格式同FASTQ一样

除了上述11列外，可以有额外列：

image

第二列FLAG每个数值含义如下，如果符合下面多种情况，则为以下数字之和：

1 read是pair中的一条（read表示本条read，mate表示pair中的另一条read）
2 pair一正一负完美的比对上
4 这条read没有比对上
8 mate没有比对上
16 这条read反向比对
32 mate反向比对
64 这条read是read1
128 这条read是read2
256 第二次比对
512 比对质量不合格
1024 read是PCR或光学副本产生
2048 辅助比对结果

其中第六列Extended CIGAR ：

M： match/mismatch
I ：插入 insertion(和参考基因组相比）
D: 删除 deletion(和参考基因组相比）
N： 跳跃 skipped(和参考基因组相比）
S： 软剪切 soft clipping ，（表示unaligned,）
H： 硬剪切 hard clipping （被剪切的序列不存在于序列中）
P： 填充 padding(表示参考基因组没有，而reads里面含有位点

2）Bam (Binary Alignment/Map)

bam文件是Sam 文件的二进制压缩格式，保留了与sam 完成相同的内容信息。SAM/BAM 文件可以是未排序的，但是按照坐标(coodinate)排序可以线性的监控数据处理过程。samtools可以用来转化bam/sam文件，可以merg，sort aligment,可以去除duplicate,可以call snp及indels.