Bioinformatics BAM & SAM Files101

Lewis LovelockLewis Lovelock
12 min read

🎉 欢迎回来！ 🎉 准备好继续探索生物信息学的奇妙世界了吗？ 这次我们要攻克的是 BAM/SAM 文件！ 🚀 如果你已经掌握了VCF文件，那么BAM/SAM文件就像是它的“前传”或者“幕后花絮”。 理解 BAM/SAM 文件，你就能深入了解基因组数据的原始形态，以及变异是如何被检测出来的。 😎

教程大纲：

  1. BAM/SAM 文件是啥？为啥重要？ 🤔 （用大白话解释 BAM/SAM 的用途）

  2. SAM 文件格式解剖 🦴 （拆解 SAM 文件的结构，理解其组成部分）

    • Header (头部): 描述比对过程和数据来源的元数据

    • Alignment Section (比对区段): 每行代表一条 reads 的比对信息

  3. 字段逐个击破！ 🎯 （详细解释 SAM 文件每一列的含义，结合例子，一看就明白）

    • Header 行的类型和常见标签

    • Alignment Section 的 11 个 mandatory 字段

  4. BAM 文件是啥？SAM 的好兄弟！ 👯‍♂️ （解释 BAM 文件和 SAM 文件的关系，以及 BAM 的优势）

  5. SAM 文件示例，实战演练！ 🤩 （用一个简单的 SAM 文件例子，串联前面学的知识）

  6. BAM/SAM 文件处理工具箱 🧰 （介绍常用的 BAM/SAM 工具，并给出简单示例）

    • samtools

    • bamtools

    • bedtools (在 BAM/SAM 文件操作中的应用)

    • IGV (Integrative Genomics Viewer)

  7. 实战演练！常见 BAM/SAM 分析场景 🚀 （告诉你 BAM/SAM 文件在实际分析中怎么用）

  8. 总结与进阶 🎓 （回顾重点，指引你更深入学习的方向）

Let's get started! 🏊‍♀️

1. BAM/SAM 文件是啥？为啥重要？ 🤔

还记得我们之前说的基因组“蓝图” 📘 吗？ 这次，我们有了测序仪，它可以像“复印机”一样，把基因组 “蓝图” 复制成无数份 DNA 片段 (reads)。 BAM/SAM 文件就是用来记录这些 reads 比对回参考基因组 的结果！ 🗺️

想象一下：

  • Reads (DNA 片段): 就像是基因组“蓝图”被撕碎成的小纸条 📜。

  • 比对 (Alignment): 我们要把这些小纸条重新拼回到原始的“蓝图”上，找到它们在基因组上的位置。 🧩

  • BAM/SAM 文件: 记录了每一张小纸条 (read) 被拼回到“蓝图” (参考基因组) 的哪个位置，以及拼对的质量、方向等等信息。 就像是拼图游戏的“结果报告” 📄。

为啥重要？

  • 变异检测的基础 🧬: VCF 文件中的变异信息，就是从 BAM/SAM 文件中 “挖掘” 出来的！ 只有先将 reads 比对回参考基因组，才能发现哪些位置和参考基因组不一样，从而检测到变异。

  • 基因表达分析 (RNA-Seq) 🎤: RNA-Seq 数据也是比对回基因组或转录组，BAM/SAM 文件记录了 RNA 片段的来源和数量，用于定量基因表达水平。

  • ChIP-Seq, ATAC-Seq 等表观遗传学研究 epigenetics: 这些技术产生的 reads 也要比对回基因组，BAM/SAM 文件用于分析蛋白质结合位点、染色质开放区域等信息。

  • 基因组可视化 👁️: BAM/SAM 文件可以加载到基因组浏览器 (例如 IGV, UCSC Genome Browser) 中，直观地查看 reads 在基因组上的分布，帮助我们理解基因组结构和功能。

  • 质量控制 (QC) ✅: BAM/SAM 文件可以用于评估测序质量、比对质量、文库复杂度等，确保实验数据的可靠性。

总之，BAM/SAM 文件是基因组学研究中最基础、最核心的数据格式之一。 理解它们，你才能真正 “看懂” 基因组测序数据！ 💎

2. SAM 文件格式解剖 🦴

SAM (Sequence Alignment/Map) 文件是 文本格式，你可以用任何文本编辑器打开它。 它的结构也分为两部分：

  • Header Section (头部区段): 以 @ 开头的行，提供关于比对过程和数据来源的元数据信息，例如：

    • SAM 文件格式版本 (@HD)

    • 参考序列信息 (@SQ)，例如染色体名称和长度

    • Read Group 信息 (@RG)，例如样本、文库、测序平台等

    • Program 信息 (@PG)，例如比对软件和版本

    • Comment 行 (@CO)，注释信息

    • Header 就像是 SAM 文件的 “说明书”，告诉程序如何解析后面的比对数据。

  • Alignment Section (比对区段): 从没有 @ 开头的行开始，每一行代表一条 read 的比对信息。 每一列用制表符 \t 分隔。

    • Alignment Section 才是 BAM/SAM 文件的 “数据主体”，记录了每条 read 的比对结果。

简单来说，SAM 文件也像一个表格，Header 是 “表头说明”，Alignment Section 是 “数据内容”，但 Header 更加结构化，用不同的 @ 开头的行类型来组织信息。 📊

3. 字段逐个击破！ 🎯

我们先来看 Header Section，再详细讲解 Alignment Section 的每一列。

Header Section 行类型和常见标签:

Header Section 的每一行都以 @ 开头，第二个字符表示行类型。 常见的行类型有：

  • @HD (Header): 文件头行，必须是 SAM 文件的第一行。 描述 SAM 文件格式版本和排序顺序。

    • Required tags (必须标签):

      • VN (Version): SAM 格式版本，例如 1.6。

    • Optional tags (可选标签):

      • SO (Sorting order): 比对结果的排序顺序，例如 unsorted, queryname, coordinate, lexicographical, unknown。

示例: @HD\tVN:1.6\tSO:coordinate (SAM 格式版本 1.6，按坐标排序)

  • @SQ (Sequence Dictionary): 参考序列字典行，每条参考序列 (例如每条染色体) 对应一行。 描述参考序列的名称和长度。

    • Required tags (必须标签):

      • SN (Sequence name): 参考序列名称，例如 chr1, chr2, chrX。

      • LN (Sequence length): 参考序列长度，碱基数。

    • Optional tags (可选标签):

      • AS (Assembly ID): 参考基因组 assembly ID。

      • M5 (MD5 checksum): 参考序列 MD5 校验和。

      • UR (URI): 参考序列 URI。

示例: @SQ\tSN:chr1\tLN:248956422 (染色体 chr1，长度 248956422)

  • @RG (Read Group): Read Group 行，每个 Read Group 对应一行。 用于区分来自不同样本、文库或测序批次的 reads。 Read Group 信息会记录在 Alignment Section 的 FLAG 字段中。

    • Required tags (至少需要 ID 标签):

      • ID (Read group ID): Read Group 的唯一标识符，在整个 SAM 文件中必须唯一。

    • Recommended tags (推荐标签):

      • SM (Sample): 样本 ID。

      • LB (Library): 文库 ID。

      • PL (Platform): 测序平台，例如 ILLUMINA, OXFORD_NANOPORE, PACBIO。

    • Optional tags (可选标签):

      • PU (Platform Unit): Platform Unit (例如 flowcell lane)。

      • CN (Sequencing Center): 测序中心。

      • DT (Date): 测序日期。

示例: @RG\tID:sample1_lib1\tSM:sample1\tLB:lib1\tPL:ILLUMINA (Read Group ID: sample1_lib1, 样本: sample1, 文库: lib1, 平台: ILLUMINA)

  • @PG (Program): Program 行，每个运行过的程序 (例如比对软件) 对应一行。 记录程序的名称、版本、命令行参数等信息，用于追溯数据处理流程。

    • Required tags (至少需要 ID 标签):

      • ID (Program ID): 程序 ID，在 @PG 行中必须唯一。

    • Optional tags (可选标签):

      • PN (Program name): 程序名称，例如 bwa, bowtie2。

      • VN (Program version): 程序版本，例如 0.7.17.

      • CL (Command line): 命令行参数。

示例: @PG\tID:bwa\tPN:bwa\tVN:0.7.17\tCL:bwa mem ref.fa reads.fq.gz (程序 ID: bwa, 名称: bwa, 版本: 0.7.17, 命令行: bwa mem ref.fa reads.fq.gz)

  • @CO (Comment): Comment 行，用于添加任意注释信息。

    示例: @CO\tThis SAM file was generated by ... (注释: 这个 SAM 文件由 ... 生成)

Alignment Section 的 11 个 mandatory 字段:

Alignment Section 的每一行代表一条 read 的比对信息，共有 11 个 mandatory 字段 (列)，按照固定顺序排列，用制表符 \t 分隔：

  1. QNAME (Query NAME) 🏷️

    • 含义: Read 的名称 (Query Name)。 通常来自 FASTQ 文件的 read ID。

    • 示例: read_123, ERR123456.1.1

    • 重要性: 标识每条 read 的身份。 对于 paired-end reads, 通常 pair 的两个 reads 有相似的 QNAME (例如只在末尾加 /1 和 /2 区分)。

  2. FLAG (FLAG) 🚩

    • 含义: Bitwise FLAG (标志位)。 用一个数字表示 read 的比对状态和特征。 每个 bit 位代表一种状态 (例如是否 paired, 是否 properly paired, 是否 unmapped, strand 方向等等)。

    • 数值和含义对应关系 (部分):

      • 0x1 (1): template having multiple segments in sequencing (paired-end or mate-pair reads) (是否 paired-end 或 mate-pair reads)

      • 0x2 (2): each segment properly aligned according to the aligner (是否 properly paired)

      • 0x4 (4): segment unmapped (read 未比对到参考基因组)

      • 0x8 (8): next segment in the template unmapped (pair read 的另一端未比对到参考基因组)

      • 0x10 (16): SEQ being reverse complemented (read 序列是否反向互补)

      • 0x20 (32): next segment in the template SEQ being reverse complemented (pair read 的另一端的序列是否反向互补)

      • 0x40 (64): first segment in template (对于 paired-end reads, 是否是 first read)

      • 0x80 (128): last segment in template (对于 paired-end reads, 是否是 last read)

      • 0x100 (256): secondary alignment (secondary alignment, 非 primary alignment)

      • 0x200 (512): not passing filters (read 未通过质量过滤)

      • 0x400 (1024): PCR or optical duplicate (PCR 或光学重复)

      • 0x800 (2048): supplementary alignment (supplementary alignment, 辅助比对)

    • 重要性: 快速了解 read 的比对状态。 可以使用在线 FLAG 解码工具 (例如 SAMtools web site) 或编程工具来解析 FLAG 值。

  3. RNAME (Reference NAME) 🧬

    • 含义: Reference sequence NAME (参考序列名称)。 Read 比对到的染色体名称。 与 @SQ Header 行的 SN 标签对应。 如果 read 未比对到参考基因组 (FLAG 字段中 unmapped bit 为 1)，则 RNAME 为 *。

    • 示例: chr1, chrX, *

    • 重要性: 确定 read 比对到的染色体位置。

  4. POS (Position) 📍

    • 含义: 1-based leftmost mapping POSition (比对起始位置)。 Read 比对到参考序列的起始位置 (1-based 坐标)。 如果 read 未比对到参考基因组，则 POS 为 0。

    • 示例: 10000, 2500000, 0

    • 重要性: 与 RNAME 一起确定 read 在基因组上的位置。

  5. MAPQ (Mapping Quality) 📊

    • 含义: MAPping Quality (比对质量)。 表示比对质量的 Phred-scaled 值。 数值越高，比对质量越好，表示该 read 比对到当前位置是 “真” 的可能性越高。 计算方式与 VCF 的 QUAL 字段类似： -10 * log10(错误概率). MAPQ 值 255 表示 quality score 不可用。 如果 read 未比对到参考基因组，则 MAPQ 为 0 或 255。

    • 示例: 20, 37, 60, 255

    • 重要性: 评估比对的可靠性，用于过滤低质量比对的 reads。

  6. CIGAR (CIGAR string) 🥢

    • 含义: CIGAR string (Compact Idiosyncratic Gapped Alignment Report, 紧凑的特异性间隔比对报告)。 描述 read 与参考序列的比对细节，包括匹配 (match), 插入 (insertion), 缺失 (deletion), 跳过参考序列区域 (skipped region from the reference), soft clipping, hard clipping 等操作。

    • CIGAR 操作符:

      • M (match or mismatch): 比对上的碱基 (match 或 mismatch)。

      • I (insertion to the reference): 插入到参考序列 (read 中有碱基插入)。

      • D (deletion from the reference): 从参考序列中删除 (read 中有碱基缺失)。

      • N (skipped region from the reference): 跳过参考序列的区域 (例如 splice junction)。

      • S (soft clipping): soft clip 的碱基 (read 末端部分碱基未比对上，但序列仍然保留在 BAM/SAM 文件中)。

      • H (hard clipping): hard clip 的碱基 (read 末端部分碱基被完全去除，序列不在 BAM/SAM 文件中)。

      • P (padding): padding (很少用)。

      • = (sequence match): 精确匹配。

      • X (sequence mismatch): mismatch。

    • CIGAR 字符串格式: <length><operator><length><operator>... 例如 50M 表示 50 个碱基匹配, 3M2I5M 表示 3 个碱基匹配，然后 2 个碱基插入，再 5 个碱基匹配。

    • 示例: 50M, 3M2I5M, 100M8S, * (如果 read 未比对到参考基因组)

    • 重要性: 详细描述 read 的比对方式，用于识别插入缺失、剪接位点、soft clipping 等。

  7. RNEXT (Reference name of the Next segment) 🧬

    • 含义: Reference sequence name of the next segment in the template (pair read 的另一端比对到的参考序列名称)。 对于 paired-end reads, 表示 pair read 的另一端比对到的染色体名称。 如果 pair read 的另一端和当前 read 比对到同一条染色体，则 RNEXT 可以是 =. 如果 pair read 的另一端未比对到参考基因组，则 RNEXT 为 *. 对于 single-end reads, RNEXT 通常为 * 或 =.

    • 示例: chr1, =, *

    • 重要性: 用于确定 paired-end reads 的 pair relationship。

  8. PNEXT (Position of the Next segment) 📍

    • 含义: POSition of the next segment in the template (pair read 的另一端比对到的起始位置)。 对于 paired-end reads, 表示 pair read 的另一端比对到的起始位置 (1-based 坐标)。 如果 pair read 的另一端未比对到参考基因组，则 PNEXT 为 0.

    • 示例: 12000, 0

    • 重要性: 与 RNEXT 一起确定 paired-end reads 的 pair relationship 和跨度。

  9. TLEN (Template LENgth) 📏

    • 含义: Template LENgth (template 长度)。 插入片段的长度。 对于 paired-end reads, 如果 reads properly paired 且方向为 forward-reverse (FR), 则 TLEN 为正值，表示插入片段长度。 如果方向为 reverse-forward (RF), 则 TLEN 为负值。 如果 reads not properly paired 或 single-end reads, 则 TLEN 通常为 0.

    • 示例: 150, -150, 0

    • 重要性: 用于评估 paired-end reads 的插入片段大小，以及判断 reads 是否 properly paired。

  10. SEQ (SEQuence) 🧬

    • 含义: SEQuence (read 序列)。 Read 的碱基序列。 如果 read 序列在 FLAG 字段中被标记为 reverse complemented, 则这里的 SEQ 也是反向互补后的序列。 如果 read 未比对到参考基因组，则 SEQ 可以是 *.

    • 示例: ATGCGTACGT..., *

    • 重要性: Read 的原始序列信息。

  11. QUAL (QUALity) 📊

    • 含义: Base QUALity (碱基质量值)。 Read 中每个碱基的质量值，Phred-scaled 格式。 质量值字符串的长度必须与 SEQ 字段的序列长度一致，每个字符对应 SEQ 字段中相应位置碱基的质量值。 质量值字符与数值的对应关系通常是 ASCII 码偏移 33 (Sanger 格式, Illumina 1.8+)。 如果 read 未比对到参考基因组，则 QUAL 可以是 *.

    • 示例: !!!!!!!!!!!!..., *

    • 重要性: Read 中每个碱基的测序质量信息，用于评估测序质量和下游分析的可靠性。

除了这 11 个 mandatory 字段，SAM 文件还可以包含 optional fields (可选字段)，格式为 TAG:TYPE:VALUE， 用制表符 \t 分隔，放在每一行的末尾。 常见的 optional fields 包括：

  • AS:i:<integer> (Alignment score): 比对得分。

  • XS:i:<integer> (Suboptimal alignment score): 次优比对得分。

  • XF:i:<integer> (Support from forward strand): 正链支持度。

  • XR:i:<integer> (Support from reverse strand): 负链支持度。

  • MD:Z:<string> (Mismatching positions): Mismatching positions to the reference sequence。 用于描述 read 和参考序列之间的 mismatch 位置和碱基变化。 例如 MD:Z:10M5A30M 表示前 10 个碱基匹配，第 11 个碱基是 mismatch (参考序列碱基为 A)，后面 30 个碱基匹配。

  • RG:Z:<string> (Read group): Read Group ID。 与 @RG Header 行的 ID 标签对应。

Alignment Section 的每一行都包含了非常丰富的信息，描述了一条 read 的比对情况。 理解这些字段，你就能够深入分析比对结果，进行后续的基因组学分析。 🤯

4. BAM 文件是啥？SAM 的好兄弟！ 👯‍♂️

BAM (Binary Alignment/Map) 文件是 SAM 文件的 二进制压缩格式！ 它们包含的信息是完全一样的，只是存储方式不同。

BAM 文件的优势：

  • 文件体积更小 📦: 二进制格式比文本格式更紧凑，BAM 文件通常比对应的 SAM 文件小得多，节省存储空间。

  • 读取速度更快 🚀: 二进制格式更易于计算机解析，BAM 文件的读取速度比 SAM 文件快得多，尤其对于大型基因组数据。

  • 更高效的随机访问 ⚡️: BAM 文件可以被索引 (生成 BAM index 文件 .bai)，实现对 BAM 文件特定区域的快速随机访问，这对于基因组浏览器可视化和区域性分析非常重要。

简单来说，BAM 文件就是 SAM 文件的 “升级版”，更高效、更实用。 在实际应用中，我们通常使用 BAM 文件进行存储和分析，只有在需要人工查看比对结果时，才可能将 BAM 文件转换为 SAM 文件。 💻 -> 💾

工具通常可以直接处理 BAM 文件 (例如 samtools, bamtools, IGV)，无需手动转换为 SAM 文件。 当你看到 .bam 文件时，就把它当作 SAM 文件的二进制 “兄弟” 就好！ 😎

5. SAM 文件示例，实战演练！ 🤩

让我们看一个简单的 SAM 文件示例：

@HD    VN:1.6    SO:coordinate
@SQ    SN:chr1    LN:248956422
@SQ    SN:chr2    LN:242193529
@RG    ID:sample1_lib1    SM:sample1    LB:lib1    PL:ILLUMINA
@PG    ID:bwa    PN:bwa    VN:0.7.17    CL:bwa mem ref.fa reads.fq.gz
read_1    99    chr1    10000    60    50M    =    10100    150    ATGCGTACGT...    !!!!!!!!!!!!...    RG:Z:sample1_lib1    AS:i:80    XS:i:75
read_2    147    chr1    10100    60    50M    =    10000    -150    CGTACGTACG...    !!!!!!!!!!!!...    RG:Z:sample1_lib1    AS:i:80    XS:i:75
read_3    0    *    0    0    *    *    0    0    GATTACA...    !!!!!!!!!!!!...    RG:Z:sample1_lib1
read_4    16    chr2    20000    30    100M    *    0    0    ATGCGTACGT...    !!!!!!!!!!!!...    RG:Z:sample1_lib1    MD:Z:50M1A49M

解释：

  • Header Section：

    • @HD\tVN:1.6\tSO:coordinate: SAM 格式版本 1.6，按坐标排序。

    • @SQ\tSN:chr1\tLN:248956422: 参考序列 chr1，长度 248956422。

    • @SQ\tSN:chr2\tLN:242193529: 参考序列 chr2，长度 242193529。

    • @RG\tID:sample1_lib1\tSM:sample1\tLB:lib1\tPL:ILLUMINA: Read Group sample1_lib1，样本 sample1，文库 lib1，平台 ILLUMINA。

    • @PG\tID:bwa\tPN:bwa\tVN:0.7.17\tCL:bwa mem ref.fa reads.fq.gz: 比对软件 bwa 0.7.17，命令行 bwa mem ref.fa reads.fq.gz。

  • Alignment Section：

    • read_1 99 chr1 10000 60 50M = 10100 150 ATGCGTACGT... !!!!!!!!!!!!... RG:Z:sample1_lib1 AS:i:80 XS:i:75

      • QNAME: read_1

      • FLAG: 99 (paired-end, properly paired, first in pair, mate reverse strand)

      • RNAME: chr1

      • POS: 10000

      • MAPQ: 60

      • CIGAR: 50M (50 个碱基匹配)

      • RNEXT: = (same chromosome as mate)

      • PNEXT: 10100

      • TLEN: 150

      • SEQ: ATGCGTACGT... (read 序列)

      • QUAL: !!!!!!!!!!!!... (碱基质量值)

      • RG:Z:sample1_lib1: Read Group sample1_lib1

      • AS:i:80: 比对得分 80

      • XS:i:75: 次优比对得分 75

      • 含义: read_1 是 paired-end reads 的 first read, properly paired, 比对到 chr1:10000, 质量 60, 50 个碱基匹配, pair read 比对到同一条染色体 chr1:10100, 插入片段长度 150, Read Group sample1_lib1, 比对得分 80。

    • read_2 147 chr1 10100 60 50M = 10000 -150 CGTACGTACG... !!!!!!!!!!!!... RG:Z:sample1_lib1 AS:i:80 XS:i:75

      • QNAME: read_2

      • FLAG: 147 (paired-end, properly paired, last in pair, mate forward strand)

      • RNAME: chr1

      • POS: 10100

      • MAPQ: 60

      • CIGAR: 50M (50 个碱基匹配)

      • RNEXT: = (same chromosome as mate)

      • PNEXT: 10000

      • TLEN: -150

      • SEQ: CGTACGTACG... (read 序列)

      • QUAL: !!!!!!!!!!!!... (碱基质量值)

      • RG:Z:sample1_lib1: Read Group sample1_lib1

      • AS:i:80: 比对得分 80

      • XS:i:75: 次优比对得分 75

      • 含义: read_2 是 paired-end reads 的 last read, properly paired, 比对到 chr1:10100, 质量 60, 50 个碱基匹配, pair read 比对到同一条染色体 chr1:10000, 插入片段长度 -150 (方向不同), Read Group sample1_lib1, 比对得分 80。 read_1 和 read_2 是 pair reads。

    • read_3 0 * 0 0 * * 0 0 GATTACA... !!!!!!!!!!!!... RG:Z:sample1_lib1

      • QNAME: read_3

      • FLAG: 0 (single-end, unmapped)

      • RNAME: * (unmapped)

      • POS: 0 (unmapped)

      • MAPQ: 0 (unmapped)

      • CIGAR: * (unmapped)

      • RNEXT: * (unmapped)

      • PNEXT: 0 (unmapped)

      • TLEN: 0 (unmapped)

      • SEQ: GATTACA... (read 序列)

      • QUAL: !!!!!!!!!!!!... (碱基质量值)

      • RG:Z:sample1_lib1: Read Group sample1_lib1

      • 含义: read_3 是 single-end read, 未比对到参考基因组。

    • read_4 16 chr2 20000 30 100M * 0 0 ATGCGTACGT... !!!!!!!!!!!!... RG:Z:sample1_lib1 MD:Z:50M1A49M

      • QNAME: read_4

      • FLAG: 16 (single-end, reverse strand)

      • RNAME: chr2

      • POS: 20000

      • MAPQ: 30

      • CIGAR: 100M (100 个碱基匹配)

      • RNEXT: * (single-end)

      • PNEXT: 0 (single-end)

      • TLEN: 0 (single-end)

      • SEQ: ATGCGTACGT... (read 序列)

      • QUAL: !!!!!!!!!!!!... (碱基质量值)

      • RG:Z:sample1_lib1: Read Group sample1_lib1

      • MD:Z:50M1A49M: Mismatching positions, 前 50 个碱基匹配，第 51 个碱基 mismatch (参考序列碱基为 A)，后 49 个碱基匹配。

      • 含义: read_4 是 single-end read, 比对到 chr2:20000, 质量 30, 100 个碱基匹配, reverse strand, Read Group sample1_lib1, 在第 51 个碱基位置有一个 mismatch。

通过这个例子，你应该对 SAM 文件的 Header 和 Alignment Section 的内容有了更清晰的认识！ 🤓

6. BAM/SAM 文件处理工具箱 🧰

处理 BAM/SAM 文件，也需要一些强大的工具。 这里介绍几个常用的命令行工具和图形界面工具：

  • samtools 🧰: BAM/SAM 文件处理的 “瑞士军刀”！ 功能非常全面，速度快，是 BAM/SAM 分析的必备工具。 常用功能：

    • 查看 BAM/SAM 文件: samtools view input.bam (默认输出 SAM 格式到标准输出) , samtools view -H input.bam (只查看 Header), samtools view -b input.sam -o output.bam (SAM 转 BAM)

    • 排序 BAM 文件: samtools sort input.bam -o sorted.bam (按坐标排序), samtools sort -n input.bam -o sorted_by_name.bam (按 read name 排序)

    • 索引 BAM 文件: samtools index sorted.bam (生成 BAM 索引文件 sorted.bam.bai)

    • 过滤 BAM 文件: samtools view -b -q 20 input.bam > filtered.bam (过滤 MAPQ < 20 的 reads), samtools view -b -f 0x2 input.bam > properly_paired.bam (只保留 properly paired reads)

    • 统计 BAM 文件信息: samtools flagstat input.bam (read 比对状态统计), samtools stats input.bam > stats.txt (生成 BAM 统计报告), samtools depth input.bam > depth.txt (计算基因组 coverage depth)

    • 合并 BAM 文件: samtools merge merged.bam input1.bam input2.bam (合并多个 BAM 文件)

    • 提取 BAM 文件子集: samtools view -b input.bam chr1:10000-20000 > region.bam (提取 chr1:10000-20000 区域的 reads), samtools view -b -r read_list.txt input.bam > subset.bam (提取 read_list.txt 中 read names 对应的 reads)

    • 更多功能: BAM/SAM 格式转换, pileup (碱基堆叠), 变异 calling, 基因型 calling 等等...

  • bamtools 🧰: 另一个流行的 BAM 文件处理工具集，功能与 samtools 有些重叠，但也有一些独特的功能。 常用功能：

    • 过滤 BAM 文件: bamtools filter -in input.bam -out filtered.bam -mapQuality ">20" -isProperPair (过滤 MAPQ > 20 且 properly paired 的 reads)

    • 统计 BAM 文件信息: bamtools stats -in input.bam (生成 BAM 统计报告)

    • 转换为 BED 格式: bamtools convert -in input.bam -format bed > reads.bed (将 reads 转换为 BED 格式)

    • 更多功能: BAM 文件索引, 合并, 拆分, 提取区域, 编辑 Header 等等...

  • bedtools 🧰: 主要用于 BED, GFF/GTF 文件操作，但也可以处理 BAM 文件，进行区域性操作。 常用功能：

    • BAM 文件与 BED 文件 intersect: bedtools intersect -abam input.bam -b regions.bed -bamsufix bam -a -o output.bam (提取与 regions.bed 文件中区域 overlap 的 reads)

    • 计算 BAM 文件在 BED 区域的 coverage: bedtools coverage -abam input.bam -b regions.bed > coverage.txt (计算每个 BED 区域的 coverage)

  • IGV (Integrative Genomics Viewer) 👁️: 图形界面的基因组浏览器，由 Broad Institute 开发。 用于可视化 BAM/SAM 文件，直观地查看 reads 在基因组上的分布，以及变异、基因结构等信息。

    • 加载 BAM/SAM 文件和索引文件 (.bai): 可以直接拖拽 BAM 文件和索引文件到 IGV 窗口，或者使用 "File -> Load from File" 菜单加载。

    • 基因组导航: 可以跳转到特定基因或基因组区域，放大缩小，平移浏览。

    • 查看 reads 比对细节: 可以查看每条 read 的比对位置、CIGAR string、碱基质量值等详细信息。

    • 与其他数据类型整合: 可以同时加载 VCF 文件, BED 文件, GTF 文件, WIG 文件等，进行多组学数据整合可视化。

    • 用于质量控制和结果验证: 通过 IGV 可以快速检查比对质量, 变异 calling 结果, RNA-Seq coverage 等，进行质量控制和结果验证。

示例： 使用 samtools 过滤 BAM 文件

假设你有一个名为 alignments.bam 的 BAM 文件，你想过滤掉比对质量 (MAPQ) 低于 20 的 reads，并保存为新的 BAM 文件 filtered_alignments.bam。

命令：

samtools view -b -q 20 alignments.bam > filtered_alignments.bam

解释：

  • samtools view: 运行 samtools 的 view 子命令，用于查看和转换 BAM/SAM 文件。

  • -b: 输出 BAM 格式。

  • -q 20: -q 参数指定 MAPQ 阈值，20 表示只保留 MAPQ 大于等于 20 的 reads。

  • alignments.bam: 输入 BAM 文件。

  • > filtered_alignments.bam: 将标准输出重定向到文件 filtered_alignments.bam，保存过滤后的 BAM 文件。

运行这个命令后，你就会得到一个过滤后的 BAM 文件 filtered_alignments.bam，其中只包含比对质量较高的 reads。 🎉

7. 实战演练！常见 BAM/SAM 分析场景 🚀

了解了 BAM/SAM 格式和工具，我们来看看一些常见的 BAM/SAM 分析场景：

  • 质量控制 (QC): 评估测序质量、比对质量、文库复杂度等。 常用工具： samtools flagstat, samtools stats, bamtools stats, Picard CollectAlignmentSummaryMetrics, FastQC (虽然 FastQC 主要针对 FASTQ 文件，但也可以分析 BAM 文件的部分信息)。

  • 基因组可视化: 使用基因组浏览器 (例如 IGV) 直观查看 reads 在基因组上的分布，用于检查比对结果、变异位点、基因表达情况等。 常用工具： IGV, UCSC Genome Browser, JBrowse.

  • 变异检测 (Variant Calling): 从 BAM/SAM 文件中检测基因组变异 (SNP, Indel, Structural Variation)。 常用工具： GATK HaplotypeCaller, bcftools mpileup + call, FreeBayes, Strelka2 (癌症变异检测)。

  • 基因表达定量 (RNA-Seq): 使用 RNA-Seq 数据 (BAM/SAM 文件) 定量基因和转录本的表达水平。 常用工具： STAR, HISAT2 (RNA-Seq 比对软件), featureCounts, HTSeq-count, RSEM, Salmon (RNA-Seq 定量软件)。

  • ChIP-Seq/ATAC-Seq 分析: 分析 ChIP-Seq/ATAC-Seq 数据 (BAM/SAM 文件) 识别蛋白质结合位点或染色质开放区域。 常用工具： MACS2, SPP, HOMER (peak calling), deepTools (coverage 计算和可视化), DiffBind (差异结合分析)。

  • 拷贝数变异 (CNV) 分析: 从 BAM/SAM 文件中推断拷贝数变异。 常用工具： CNVkit, Control-FREEC, copycat。

  • 甲基化分析 (Bisulfite Sequencing): 分析 Bisulfite Sequencing 数据 (BAM/SAM 文件) 定量 DNA 甲基化水平。 常用工具： Bismark, BS-Seeker2, methylKit.

这些只是 BAM/SAM 文件分析的一些常见应用，实际上 BAM/SAM 文件是基因组学研究的基石，几乎所有基因组数据分析流程都离不开 BAM/SAM 文件。 🚀

8. 总结与进阶 🎓

太棒了！ 🎉 你又完成了一个生物信息学核心知识点的学习：BAM/SAM 文件！ 现在你应该对 BAM/SAM 文件的格式、字段含义、常用工具和应用场景有了全面的了解。 你离生物信息学忍者又近了一步！ 🥷

回顾一下重点：

  • BAM/SAM 文件记录了 reads 比对回参考基因组的结果。

  • SAM 文件是文本格式，BAM 文件是二进制压缩格式，BAM 更常用。

  • SAM 文件由 Header Section 和 Alignment Section 组成。

  • Alignment Section 的每一行代表一条 read 的比对信息，包含 11 个 mandatory 字段。

  • 理解每个字段的含义是解读 BAM/SAM 文件的关键。

  • samtools, bamtools, IGV 是常用的 BAM/SAM 文件处理工具。

  • BAM/SAM 文件广泛应用于变异检测、基因表达分析、表观遗传学研究等领域。

进阶学习方向：

  • 深入了解 SAM/BAM 格式规范: 阅读 SAM/BAM 格式规范文档 (官方文档或维基百科)。

  • 学习更多 BAM/SAM 工具和高级用法: 深入学习 samtools, bamtools, bedtools, GATK 等工具的各种功能和参数，掌握更高级的 BAM/SAM 文件操作技巧。

  • 掌握编程处理 BAM/SAM 文件: 学习使用 Python 或 R 库 (例如 pysam, Rsamtools) 编程处理 BAM/SAM 文件，实现更复杂的分析流程。

  • 学习基因组比对算法和软件: 了解常用的基因组比对算法 (例如 BWA, Bowtie2, STAR, HISAT2) 和比对软件的原理和参数，深入理解比对过程。

  • 实践！实践！再实践！: 下载一些公开的 BAM/SAM 文件 (例如 ENCODE, TCGA)，使用工具进行练习，解决实际生物学问题。

基因组数据的世界充满挑战，也充满机遇。 BAM/SAM 文件是通往基因组奥秘的钥匙。 继续努力，不断探索，你将成为基因组数据的 “解码大师”！ 🔑

祝你学习进步！ Happy BAM-ing! 😄

0
Subscribe to my newsletter

Read articles from Lewis Lovelock directly inside your inbox. Subscribe to the newsletter, and don't miss out.

bioinformatics

Written by

Lewis Lovelock
Lewis Lovelock

I am a developer working in BGI located in Shenzhen. I am familiar with genomics 🧬 and coding. I love 🏀 👩🏻‍💻 and Hiphop🎵