系列文章：
二代测序方法：DNA测序之靶向重测序
NGS数据分析实践：00. 变异识别的基本流程
NGS数据分析实践：01. Conda环境配置及软件安装
NGS数据分析实践：02. 参考基因组及注释库的下载
NGS数据分析实践：03. 涉及的常用数据格式[1] - fasta和fastq格式

测序数据分析中涉及的常用格式：测序得到的是带有质量值的碱基序列(fastq格式)，参考基因组是(fasta格式)，用比对工具把fastq格式的序列比对到对应的fasta格式的参考基因序列，就可以产生sam格式的比对文件。把sam格式的文本文件压缩成二进制bam文件可以节省空间，如果对参考基因组上面的各个区段标记它们的性质，比如哪些区域是外显子、内含子、UTR等等，这就是gtf/gff格式。如果只是为了单纯描述某个基因组区域，就是bed格式文件，记录染色体号以及起始终止坐标，正负链即可。如果是记录某些位点或者区域碱基的变异，就是vcf文件格式。

fasta/fastq(测序数据)→SAM/BAM(比对)→gff/gtf(描述基因组上的结构：坐标&类型)→Bigwig/Wiggle(测序深度)→bed(描述坐标)→vcf(突变信息)

存储序列：fasta/fastq
比对结果显示的文件：sam/bam
展示注释信息：gtf/gff/bed
突变信息：vcf

注：以下绿色PPT截图来自 生信技能树。

2. sam和bam格式

SAM的全称是sequence alignment/map format。而BAM就是SAM的二进制文件(B源自binary)。

SAM格式用来支持高通量测序数据分析：
快速查找与坐标重叠的比对。例如，选择与2号染色体上的坐标323,567,334重叠的比对。（UCSC查看）
根据read的属性进行选择和过滤。例如，我们希望能够快速选择能过比对到反向链上的read。
有效地存储数据。例如，从SAM格式转化成BAM格式，单个压缩文件包含所有样本的数据，每个样本都以某种方式标记。

SAM文件注释及12列的简要说明：

SAM 格式主要包括两大部分：
(1) 标头注释部分（header section）
首先我们看一个例子：

@HD VN:1.0 SO:unsorted
@SQ SN:gi|10141003|gb|AF086833.2|LN:18959
@PG ID:bowtie2 PN:bowtie2 VN:2.2.5 CL:"bowtie2-2.2.5/bowtie2 -align -s --wrapper basic -0 -x 1976.fa -1 SRR1972739_1.fastq -2 SRR1972739_2.fastq"

标头信息可有可无，都是以@开头，用不同的tag表示不同的信息，主要有：
@HD，说明符合标准的版本、对比序列的排列顺序（这里unsorted）
@SQ，参考序列说明 （SN:gi|10141003|gb|AF086833.2|LN 是参考序列的长度）
@PG，使用的比对程序说明（这里是bowtie2）

(2) 比对结果部分（alignment section）
比对结果部分，每一行表示一个read的比对信息，包括11个必须的字段（mandatory fields）和一个可选的字段，字段之间用tag分割。必须的字段的顺序固定，不可用时，根据字段定义，可以为’0‘或者’*‘，这11个字段包括：

第一列: Query Name (QNAME)
进行reads比对时通常表示reads的名字，①如果这条reads比对到多条序列或比对到这条序列的多个位置，相同名字会出现多次。②如果是pair-end reads，相同名字会出现2次，分别表示来自于R1文件的reads和R2文件的reads；如果其mate pair reads也比对2个位置，也会出现2次，则相同名字共出现4次；如果一条reads也比对2个位置，则其mate pair比对1个位置，则共出现3次；如果其mate pair reads没有比对上序列也会出现1次（第三列显示“*”），所以pair-end测序，R1文件和R2文件同时mapping，相同reads的id最少出现2次。

第二列：FLAG
这是一种常用且高效的保存多个布尔特征值的方法。
举个简单的例子：在 SAM 格式中，当 flag 为 1，也即对应的二进制为 01 时，表示该 read 有多个测序数据 ， 一般理解为有双端测序数据 (另一条没被过滤掉)， 而 flag 为 2， 也即二进制 10 时， 表示这条 read 的多个片断都有比对结果， 通常理解为双端 reads 都比对上了， 那么就可以推断出 flag 为 3 时， 也即二进制的 11， 表示该 read 有另一端的 read 并且比对成功， 可以看到， 其实就是 01 加 10。

Flag标识对应的情况说明：

要是所得出的标记不是上述列举的数字，比如83=64+16+2+1，就是这几种情况值和。
这个网站可以直接通过输入你所得的标记数字，直接告诉其对应的信息：https://broadinstitute.github.io/picard/explain-flags.html

第三列： Reference Name (RNAME)
参考序列的编号，实际上就是比对到参考序列上的染色体号，如果注释中对SQ-SN进行了定义，这里必须和其保持一致；若没有mapping上的序列，用“*”表示，可以认为这条read没有比对上的序列，则这一行的第四，五，八，九 列是“0”，第六，七列与该列是相同的表示方法。

第四列： Position (POS)
表示read比对到RNAME这条序列的最左边的位置，如果该read能够完全比对到这条序列（CIGAR string为M）则这个位置是read的第一个碱基比对的位置，如果该read的反向互补序列比对到这条序列，则这个位置是read的反向互补序列的第一个碱基比对的位置，所以无论该read是正向比对到该序列，或是其反向互补序列比对到该序列，比对结果均是最左端的比对位置。注意是从1开始计数，若没有比对上，此处为0。

第五列：Mapping Quality (MAPQ)
比对的质量；比对的质量分数，越高说明该read比对到参考基因组上的位置越准确。

第六列：Compact Idiosyncratic Gapped Alignment Representation (CIGAR)
CIGAR 代表着简要比对信息表达式，其以参考序列为基础，使用数字加字母表示比对结果。 例如 3S6M1P1I4M
前三个碱基被剪切去除了，然后6个比对上了，然后打开了一 个缺口，有一个碱基插入，最后是4个比对上了。

对于mapping状态可分为以下几类：

M：alignment match (can be a sequence match or mismatch)
表示read可mapping到第三列的序列上，则read的碱基序列与第三列的序列碱基相同，表示正常的mapping结果，M表示完全匹配，但是无论reads与序列的正确匹配或是错误匹配该位置都显示为M。
I：insertion to the reference
表示read的碱基序列相对于第三列的RNAME序列，有碱基的插入。
D：deletion from the reference
表示read的碱基序列相对于第三列的RNAME序列，有碱基的删除。
N：skipped region from the reference
表示可变剪接位置。
P：padding (silent deletion from padded reference)
S：soft clipping (clipped sequences present in SEQ)
H：hard clipping (clipped sequences NOT present in SEQ)
clipped均表示一条read的序列被分开，之所以被分开，是因为read的一部分序列能匹配到第三列的RNAME序列上，而被分开的那部分不能匹配到RNAME序列上。
=：表示正确匹配到序列上
X：表示错误匹配到序列上

注意：
1) H只出现在一条read的前端或末端，但不会出现在中间，S一般会和H成对出现，当有H出现时，一定会有一个与之对应的S出现 。
2) S可以单独出现，而H必须有与之对应的S出现时才可能出现，不可在相同第一列的情况下单独出现。
3) N：如果是mRNA-to-genome，N出现的位置代表内含子，其它比对形式出现N时则没有具体解释 。
4) M/I/S/=/X：这些数值的加和等于第10列SEQ的长度。

第七列：MRNM (chr)
这条reads第二次比对的位置，在利用bwa mem产生sam文件时，如果该列是“*”，而第3列RNAME不是“*”则表示该reads比对到第3列显示序列名的序列上，而没有比对到其他位置。在利用bwa aln及bwa sampe比对生成的sam文件，如果和上述情况相同，则第7列为“=”，上述情况均表示该reads只比对到这一个位置 。如果第3列RNAME和第7列MRNM都为“*”，则说明这条reads没有匹配上序列，如果这条reads匹配两个序列，则第一个序列的名称出现在第3列，而第二个序列的名称出现在第7列。

第八列：MPOS (mate position)
该列表示与该reads对应的mate pair reads的比对位置，如果这对pair-end reads比对到同一条reference序列上，在sam文件中reads的id出现2次，Read1比对的第4列等于Read2比对的第8列。同样Read1比对的第8列等于Read2比对的第4列。

第九列：ISIZE
TLEN：signed observed Template LENgth （可以理解为文库插入片段长度）
如果R1端的read和R2端的read能够mapping到同一条Reference序列上（即第三列RNAME相同），则该列的值表示第8列减去第4列加上第6列的值，R1端和R2端相同id的reads其第九列值相同，但该值为一正一负，R1文件的reads和R2文件的reads，相同id的reads要相对来看。在进行该第列值的计算时，如果取第6列的数值，一定要取出现M的值，S或H的值不能取。
详见Illumina中paired end sequencing 和 mate pair sequencing。
ISIZE为正，说明amplicon的start 为POS，amplicon 的end为 pos + isize-1
ISIZE为负，说明amplicon的start为mate position，amplicon的end为mate position - isize-1
ISIZE为0，说明无法计算出amlicon的大小。

第十列：Sequence
序列片段的序列信息，就是read的碱基序列，如果是比对到互补链上则是reverse completed。如果不存储此类信息，此处为’*‘，注意CIGAR中M/I/S/=/X对应数字的和要等于序列长度。

第十一列：ASCII
read质量的ASCII编码。

第十二列：Optional fields
可选的区域
格式如：TAG:TYPE:VALUE，其中TAG有两个大写字母组成，每个TAG代表一类信息，每一行一个TAG只能出现一次，TYPE表示TAG对应值的类型，可以是字符串、整数、字节、数组等。如，AS:i:-10

find . -name *.sam #查找sam格式的文件
bwa mem hg19  MiSeq_SOP/tmp.fq > tmp.sam 
cat tmp.sam | less -S
samtools view -bS tmp.sam > tmp.bam
samtools view tmp.bamsamtools sort tmp.bam > tmp.sorted.bam
samtools view tmp.bam | less -S #无序
samtools view tmp.sorted.bam | less -S #按染色体顺序
samtools view -h tmp.sorted.bam | less -S #-h：查看头文件

其余文件格式见后续。

参考阅读：
http://genome.ucsc.edu/FAQ/FAQformat.html
Sam格式详解手册：http://samtools.github.io/hts-specs/SAMv1.pdf
Sam格式相关文献：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2723002/
https://www.cnblogs.com/yahengwang/p/10695981.html
https://www.jianshu.com/p/a9deb1adb5a6