BAM/SAM文件的一些小知识
前言
如果不是在陈老师这读博，然后开始折腾BAM/SAM文件，我估计这辈子都不会了解到这么多东西吧
SAM/BAM简介
Sequence Alignment Map (SAM) is a text-based format for storing biological sequences aligned to a reference sequence developed by Heng Li and Bob Handsaker et al.
- SAM是李恒大佬在冷泉港鼓捣出来的一种文件格式，存储了测序获得的信息，map到基因组后的各种信息。SAM格式为纯文本格式，字里行间压缩了极大的信息。
- BAM则是SAM的二进制版，在SAM的基础上运用二进制编码，又极大的压缩了SAM文件的体积。
在这个文件基础上，大家常用的各种alignment软件基本都支持SAM/BAM格式，围绕着这个文件格式，李恒等还开发了多个软件和不同语言版本的依赖库，以供大家使用该文本格式储存信息，并且在该文件基础上进行开发。
- samtools: 对SAM/BAM等格式进行各种相关操作的软件，功能最丰富。
- htslib: samtools的C接口，可以通过该库，在C语言中完成samtools的同等功能。
- htsjdk: Java接口
- pysam: Python接口
具体用法，推荐看各自的文档，比如javadoc或者python doc
格式
以下为SAM格式示例（BAM格式参照SAM即可）
@HD VN:1.6 SO:coordinate
@SQ SN:ref LN:45
r001 99 ref 7 30 8M2I4M1D3M = 37 39 TTAGATAAAGGATACTG *
r002 0 ref 9 30 3S6M1P1I4M * 0 0 AAAAGATAAGGATA *
r003 0 ref 9 30 5S6M * 0 0 GCCTAAGCTAA * SA:Z:ref,29,-,6H5M,17,0;
r004 0 ref 16 30 6M14N5M * 0 0 ATAGCTTCAGC *
r003 2064 ref 29 17 6H5M * 0 0 TAGGC * SA:Z:ref,9,+,5S6M,30,1;
r001 147 ref 37 30 9M = 7 -39 CAGCGGCAT * NM:i:1
SAM文件主要由两个部分构成
- header：标记了该SAM文件的一些基本信息，比如版本、按照什么方式排序的、Reference信息等等。
- 本体：每行为一个reads，不同列记录了不同的信息，列与列之间通过tab分隔。
每列的含义：
- QNAME：测序的reads的名字。
- FLAG：二进制数字之和，不同数字代表了不同的意义；比如正负链，R1/R2（双端测序的哪一端）等。
- RNAME：map到参考基因组后的染色体名称。
- POS：1-based 基因组起始位点。
- MAPQ：map的质量。
- CIGAR：一个数字与字母交替构成的字符串，标记了这段reads不同位置的match情况。不同字母的含义后边介绍。
- RNEXT：如果是pair-end测序，这个为mate（另一端中对应的）的read的染色体名称；否则为下一条read的染色体名称。
- PNEXT：同上，read对应的起始位点。
- TLEN：模板的长度。
- SEQ：序列。
- QUAL：序列的质量打分（fasta文件中的那个）。
flag的含义：
flag tests 查询flag中具体数字含义，以及测试flag内容的网页
flag比较特殊的是一点在于flag最后的数字中包含了这个数即为true（包含该特性），不包含这个数字即为false（不包含该特性）。
检测方法也比较巧妙，通过二进制与进行计算。比如：16 and 16 -> 16。结果与测试的数字一致，则表明flag中包含该数字。否则，不包含。
当然，不同语言版本的接口也直接提供了api来直接获取这些特性，不必在进行人工的计算。
cigar的含义
cigar中会包含数字，代表了特定match持续了多少nt；以及不同的字符，代表了不同的match情况。
cigar string examples：
30S512M216N12S （30nt soft clip -> 512nt exact match -> 216nt skipped region -> 12nt soft clip）
30S (30nt soft clip)
40M (40nt exact match)
其中不同的字符及其含义如下：
另一个版本的说明来源
注意
双端测序，正负链的问题
双端测序很麻烦，mate和read两条序列的flag信息基本都是完全相反的，包括strand等信息。那么如何判断这一对reads测到的到底是哪条链就成了问题。
- 最稳妥的方法是，通过GT/AG规则来确定。缺点就是你很难提取到GT/AG这个序列，最起码在我的测试中如此。
- 在我的理解中，cigar中的I、D、N三个字符代表的区域不计入位点序列中。那么N区域的第一个位点周边的序列即为内含子周边的序列，然而在这个位点周边，很难获取到标准的GT/AG序列及其互补序列。可能由junctions的突变造成的，即N区域并不一定是标注的内含子区域
- STAR应该是通过这个途径识别链方向的
alignSJstitchMismatchNmax 0 -1 0 0 4*int>=0: maximum number of mismatches for stitching of the splice junctions (-1: no limit). (1) non-canonical motifs, (2) GT/AG and CT/AC motif, (3) GC/AG and CT/GC motif, (4) AT/AC and GT/AT motif.
- 另外就是根据readNegativeStrand和mateNegativeStrand，配合以First in pair以及second in pair进行判断。
- 这个办法，依赖于是否确定R1和R2的建库方式和建库方向，如果我们确定R1是某个特定方向，那么我们就能够通过这四个参数的组合获取到正确的方向。
- 目前我已知的regtools就是通过这种途径进行链的识别的。
//Get the strand from the bitwise flag void JunctionsExtractor::set_junction_strand_flag(bam1_t *aln, Junction& j1) { uint32_t flag = (aln->core).flag; // 从flag中提取出readNegativeReversed，mateNegativeReversed, first in pair和second in pair的信息 int reversed = (flag >> 4) % 2; int mate_reversed = (flag >> 5) % 2; int first_in_pair = (flag >> 6) % 2; int second_in_pair = (flag >> 7) % 2; // 下面就是判断正负链的过程 // strandness_ is 0 for unstranded, 1 for RF, and 2 for FR int bool_strandness = strandness_ - 1; int first_strand = !bool_strandness ^ first_in_pair ^ reversed; int second_strand = !bool_strandness ^ second_in_pair ^ mate_reversed; char strand; if (first_strand){ strand = '+'; } else { strand = '-'; } //cerr << "flag is " << flag << endl; // if strand inferences from first and second in pair don't agree, we've got a problem if (first_strand == second_strand){ j1.strand = string(1, strand); } else { j1.strand = string(1, '?'); } //cerr <<"flag strand is " << j1.strand << endl; return; } - 根据参数推测，hisat2可能也是这个方式。
--fr/--rf/--ff -1, -2 mates align fw/rev, rev/fw, fw/fw (--fr)
