sam格式的結構和意義_BAM/SAM檔案格式的一些小知識

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分隔。

每列的含義：

column

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

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)

其中不同的字元及其含義如下：

cigar from official

另一個版本的說明來源

cigar from blog

注意

雙端測序，正負鍊的問題

雙端測序很麻煩，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 <

return;

}

根據參數推測，hisat2可能也是這個方式。

--fr/--rf/--ff -1, -2 mates align fw/rev, rev/fw, fw/fw (--fr)