BAM/SAM文件的一些小知識

前言

如果不是在陳老師這讀博五鲫，然后開始折騰BAM/SAM文件，我估計這輩子都不會了解到這么多東西吧

SAM/BAM簡介

See Wiki

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的基礎(chǔ)上運用二進制編碼背镇，又極大的壓縮了SAM文件的體積。

在這個文件基礎(chǔ)上泽裳，大家常用的各種alignment軟件基本都支持SAM/BAM格式瞒斩，圍繞著這個文件格式，李恒等還開發(fā)了多個軟件和不同語言版本的依賴庫涮总，以供大家使用該文本格式儲存信息胸囱，并且在該文件基礎(chǔ)上進行開發(fā)。

• samtools: 對SAM/BAM等格式進行各種相關(guān)操作的軟件瀑梗，功能最豐富烹笔。
• 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文件主要由兩個部分構(gòu)成

• header：標記了該SAM文件的一些基本信息允蜈，比如版本、按照什么方式排序的、Reference信息等等饶套。
• 本體：每行為一個reads漩蟆，不同列記錄了不同的信息，列與列之間通過tab分隔凤跑。

每列的含義：

column

1. QNAME：測序的reads的名字爆安。
2. FLAG：二進制數(shù)字之和，不同數(shù)字代表了不同的意義仔引；比如正負鏈，R1/R2（雙端測序的哪一端）等褐奥。
3. RNAME：map到參考基因組后的染色體名稱咖耘。
4. POS：1-based 基因組起始位點。
5. MAPQ：map的質(zhì)量撬码。
6. CIGAR：一個數(shù)字與字母交替構(gòu)成的字符串儿倒，標記了這段reads不同位置的match情況。不同字母的含義后邊介紹呜笑。
7. RNEXT：如果是pair-end測序夫否，這個為mate（另一端中對應(yīng)的）的read的染色體名稱；否則為下一條read的染色體名稱叫胁。
8. PNEXT：同上凰慈，read對應(yīng)的起始位點。
9. TLEN：模板的長度驼鹅。
10. SEQ：序列微谓。
11. QUAL：序列的質(zhì)量打分（fasta文件中的那個）。

flag的含義：

flag tests 查詢flag中具體數(shù)字含義输钩，以及測試flag內(nèi)容的網(wǎng)頁

flag

flag比較特殊的是一點在于flag最后的數(shù)字中包含了這個數(shù)即為true（包含該特性）豺型，不包含這個數(shù)字即為false（不包含該特性）。

檢測方法也比較巧妙买乃，通過二進制與進行計算姻氨。比如：16 and 16 -> 16。結(jié)果與測試的數(shù)字一致剪验，則表明flag中包含該數(shù)字肴焊。否則，不包含碉咆。

當然抖韩，不同語言版本的接口也直接提供了api來直接獲取這些特性，不必在進行人工的計算疫铜。

cigar的含義

cigar中會包含數(shù)字茂浮，代表了特定match持續(xù)了多少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測到的到底是哪條鏈就成了問題。

• 最穩(wěn)妥的方法是属桦，通過GT/AG規(guī)則來確定熊痴。缺點就是你很難提取到GT/AG這個序列，最起碼在我的測試中如此聂宾。
  • 在我的理解中果善，cigar中的I、D系谐、N三個字符代表的區(qū)域不計入位點序列中巾陕。那么N區(qū)域的第一個位點周邊的序列即為內(nèi)含子周邊的序列，然而在這個位點周邊纪他，很難獲取到標準的GT/AG序列及其互補序列鄙煤。可能由junctions的突變造成的茶袒，即N區(qū)域并不一定是標注的內(nèi)含子區(qū)域
  • STAR應(yīng)該是通過這個途徑識別鏈方向的

    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.
    

• 另外就是根據(jù)readNegativeStrand和mateNegativeStrand梯刚，配合以First in pair以及second in pair進行判斷。
  • 這個辦法弹谁，依賴于是否確定R1和R2的建庫方式和建庫方向乾巧，如果我們確定R1是某個特定方向，那么我們就能夠通過這四個參數(shù)的組合獲取到正確的方向预愤。
  • 目前我已知的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;
    }
    

  • 根據(jù)參數(shù)推測植康，hisat2可能也是這個方式旷太。

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