模組 4 · Data literacy

資料格式讀寫能力

從 FASTA、FASTQ、BAM 到 VCF，先建立檔案層級與範例，再拆解 BAM 的 CIGAR、HP／MM/ML 與 VCF 的 GT／PS。

本模組學習目標

說出 FASTA、FASTQ、BAM、VCF 在 NGS 資料鏈中的角色，以及每種格式的一筆 record 代表什麼
用一個小型例子，從 FASTQ 的 read 追到 BAM 的 alignment record，再追到 VCF 的位點紀錄
區分檔案格式、record、欄位、BAM tag 與 VCF FORMAT，認出 GT、PS、HP 與 MM/ML
在知道 BAM 的角色後，判讀 CIGAR 字串並計算它消耗的參考序列與 read 長度
區分一般沒有 HP 的 untagged read，與帶有 somatic signal 但來源未定的 HP3

為什麼重要

做 NGS 分析時，你不會只拿到一個「結果檔」。同一批 DNA 觀測會在流程中留下不同格式：有的檔案保存參考序列，有的保存原始 read，有的保存每條 read 如何對齊，有的則把許多 read 的證據整理成變異位點。

所以本模組的第一個問題不是「CIGAR 怎麼解碼」，而是：我現在看到的檔案正在描述哪一層資料？先掌握 FASTA、FASTQ、BAM、VCF 的角色，再往下讀 record、欄位與 tag，最後才把 CIGAR 放回 BAM 的位置。

從左到右是參考序列、原始 read、read-level 的對齊紀錄與 site-level 的變異紀錄。下方的閱讀順序是本模組的主線：先判斷檔案格式，再判斷一筆 record，接著讀欄位，最後才讀 BAM tag 或 VCF FORMAT。

本節會用一個很小的 toy locus 貫穿範例。它不是可直接拿去分析的完整檔案，而是把同一個位置在不同資料層級的樣子並排，讓你能回答「這一列在記錄什麼」以及「我要回哪個檔案找證據」。

概念與互動

一、先建立檔案地圖

格式／檔案	它保存什麼	一筆 record 通常代表什麼	它回答的問題
`FASTAFASTA以文字保存 DNA 或 RNA 序列的格式；在本模組中，FASTA 主要用來提供比對所需的參考基因體序列。A text format for DNA or RNA sequences; here it mainly supplies the reference genome used for alignment.完整條目 →`	參考基因體序列	一條序列或 contig 的 header 與鹼基	參考座標上的標準序列是什麼？
`FASTQFASTQ存放 read 序列與每個鹼基品質分數的文字格式；每筆紀錄由四行組成。Text format holding read sequences plus per-base quality scores; four lines per record.完整條目 →`	定序儀產生的 read 序列與品質	四行 read record：名稱、序列、分隔線、品質	定序儀觀察到了哪些片段？
`BAMBAM已比對到參考基因體的 read 的二進位格式（SAM 的壓縮版）。每一列是一條 read，帶著它的位置、CIGAR、品質與各種 tag。Binary format for reads aligned to the reference. One record per read, carrying position, CIGAR, quality and tags.完整條目 →`	已對齊至參考的 read；它是 SAM 文字格式的二進位版本	一筆 alignment record，通常對應一條 read	這條 read 被放在參考的哪裡？對齊得怎樣？
`.bai`	BAM 的索引，不是另一批生物學觀測	由座標建立的查詢索引	如何快速取出某個區間的 BAM records？
`VCFVCF存放 variant 的文字格式。每一列是一個位點，記錄位置、ref/alt allele、品質、FILTER 與各種 FORMAT 欄位。Text format for variants: one line per site with position, ref/alt alleles, quality, FILTER and FORMAT fields.完整條目 →`	由 read 證據提出的變異紀錄	一列通常代表一個基因體位點	哪個位置被提出為 REF／ALT 變異候選？

FASTA 與 FASTQ 不是同一份資料的兩種副本。FASTA 在這裡提供比對基準；FASTQ 保存尚未有參考座標的 read。把 FASTQ 與 FASTA 比對後，才會得到 BAM；再由 BAM 的多條 read 證據整理出 VCF。VCF 的一列也不等於一條 read，而且通常是 caller 提出的候選，不是脫離原始證據的事實宣告。

二、用同一個小例子追過四種格式

假設我們只看 chr7 上一小段 toy reference。為了讓字串短一點，下面的座標從 100 開始；真實檔案會有完整的 chromosome sequence、header 與更多欄位。

1. FASTA：比對的參考序列

reference.fasta
>chr7
ACCTGATC

這段序列的 chr7:100–107 依序是 A C C T G A T C。FASTA 告訴我們「參考是什麼」，但還沒有任何 sample read，也沒有變異結論。

2. FASTQ：尚未定位的 read 觀測

reads.fastq
@read_001
ACCTGATC
+
IIIIIIII
@read_002
ACCTGGTC
+
IIIIIIII

每筆 FASTQ record 有四行：名稱、序列、+ 分隔線、每個鹼基的品質字串。這裡的 I 只是 toy example 的品質字元；重點是 read 還沒有 chr7:105 這種參考座標。

3. BAM：把每條 read 放回參考座標

BAM 是二進位檔，不能像文字檔一樣直接讀出一整列。為了教學，下面用「從 SAM／samtools view 選出的欄位」簡化顯示；這不是 BAM 的二進位內容：

QNAME       RNAME POS  MAPQ CIGAR SEQ       TAGS
read_001   chr7  100  60   8M    ACCTGATC  HP:i:1
read_002   chr7  100  60   8M    ACCTGGTC  HP:i:2

這兩條 read 都從 chr7:100 開始，並各自跨過 chr7:105。第二條 read 在這個位置讀到 G，而 reference 是 A。此時我們仍在讀每條 read 的觀測；還沒有把它們壓縮成一列 VCF。

4. VCF：把同一位置整理成一列變異紀錄

sample.vcf
##fileformat=VCFv4.3
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample
chr7   105 .  A   G   60   PASS   .    GT:PS  0|1:105

VCF 這一列把前面的觀察摘要成：在 chr7:105，reference allele 是 A，alternate allele 是 G。GT:PS 是這個 sample 欄位的格式說明，後面的 0|1:105 才是對應值。真正要確認這一列是否由可靠 read 支持，仍要用 CHROM 與 POS 回到 BAM。

三、從檔案讀到細節：四個層級

層級	先問什麼	本例
檔案格式	這份檔案的資料模型是什麼？	`FASTQ` 是 read；`BAM` 是對齊後的 read；`VCF` 是位點。
record／row	一筆資料代表一個什麼東西？	BAM 的 `read_002` 是一條 alignment record；VCF 的 `chr7:105` 是一個位點 row。
固定欄位	這筆資料的座標、序列或狀態放在哪一欄？	BAM 的 `POS`、`MAPQ`、`CIGAR`；VCF 的 `REF`、`ALT`、`FILTER`。
tag／FORMAT	工具附加了什麼額外判讀？	BAM 的 `HP`、`MM`、`ML`；VCF 的 `GT`、`PS`。

四層是由外往內的，順序不能跳。縮排就是「再鑽進去一層」：先知道這是哪一種檔案，才知道一筆紀錄是什麼；知道一筆是什麼，才知道欄位在描述誰。

BAM tag 與 VCF FORMAT 欄位都像短短的鍵值，但不要把它們混為一談：HP 是寫在 read record 上的 BAM tag；GT 與 PS 是寫在 VCF sample 欄位中的 FORMAT 子欄位。先辨認外層格式，才不會把一個 read 的標籤誤當成整個位點的結論。

同一個 chr7:105，在 BAM 是六筆紀錄、每筆各帶一個 HP，在 VCF 是一列、整列只有一組 GT 與 PS。所以 HP:Z:1 講的是「這條分子來自哪一邊」，GT 0|1 講的是「這個位置的兩份拷貝各是什麼」——它們不是同一個層級的話。

四、先讀 VCF 的關鍵欄位

VCF 部分	本例	初學者先記住什麼
`CHROM`／`POS`	`chr7`／`105`	用來定位並回到 BAM 查詢；座標必須和 reference build 一致。
`REF`／`ALT`	`A`／`G`	這一列比較的 reference allele 與 alternate allele。
`FILTER`	`PASS`	表示通過該流程設定的篩選；不是「任何情況下都已證明為真」。
`GT`	`0/1` 或 `0\|1`	`0` 是 REF、`1` 是第一個 ALT；斜線代表未定相，直線代表已提供相位順序。
`PSPS tag phase set標示某群 variants 屬於同一 phase block 的編號。同一 PS 內的相位方向可相互比較；不同 PS 的 HP1/HP2 方向彼此獨立。Phase set identifier grouping variants phased together. Haplotype orientation is only consistent within one PS.完整條目 →`	`105`	同一個 PS 的位點屬於同一個 phase set；不同 PS 的 HP1／HP2 方向不可直接當成一致。

0|1 只說明兩個 allele 的相位順序已被表示出來，不代表第一條一定是父系，也不代表這一列不需要回看 BAM。PS 是一個分組編號，不是支持 read 的數量。

五、再讀 BAM：CIGAR 只是其中一個欄位

讀 BAM 時，可以把一筆 alignment record 想成一列「read 如何被放回參考」的資料。POS 告訴你從哪個參考座標開始，MAPQ 是比對位置的信心摘要，SEQ 是 read 序列；CIGARCIGAR描述一條 read 如何對上參考基因體的緊湊字串。例如 5M1I5M 表示兩端各有 5 個對齊欄位，中間有 1 個相對於參考序列的插入；M 可能代表吻合，也可能代表錯配。A compact string describing how read bases and reference positions are consumed by an alignment. For example, 5M1I5M has five alignment columns, one insertion, then five alignment columns; M may contain matches or mismatches.完整條目 → 則用一串數字與字母，描述 read 與參考序列之間各段如何對齊。

先用固定規則讀 CIGAR：每個操作碼分別消耗「參考序列」或「read」。互動元件只做這一個局部任務；它不是另一種檔案格式。

操作碼	意思	消耗參考	消耗 read
`M`	有對齊的欄位；可能吻合，也可能錯配	✓	✓
`I`	insertion：read 有、參考沒有	✗	✓
`D`	deletion：參考有、read 沒有	✓	✗
`S`	soft clip：read 保留，但這段沒有納入對齊	✗	✓
`H`	hard clip：該段序列不保留在這筆 alignment record 中	✗	✗

M 不等於完全吻合。如果要知道 M 中哪些鹼基真的不同，要一起看 read 序列、參考序列或 MD tag；不能只靠 CIGAR 的 M 判斷。

六、soft clip：BAM 對齊資訊的一個判讀例子

當 read 的一端無法可靠地對上參考時，比對器可能把那一端標成 S。soft-clipped 的序列仍保留在 BAM，所以它和 hard clip 不同；只是該段不消耗參考座標。

上方先對照參考與兩條 read 的端部 clipping；下方比較 clipping 分散與集中。集中在同一座標的 soft clip 可列為結構變異斷點候選，但也可能由低品質、重複序列或其他 alignment 問題造成，仍需其他證據確認。

把候選訊號當成候選，不要當成結論

若多條 read 在相近位置、以相同方向出現 soft clip，這比單一 read 的 clipping 更值得檢查，因為它可能表示 sample 與 reference 在該處的結構不同。但「集中」本身不能排除重複序列、低 mapping quality、定序錯誤或比對器選擇；應搭配 coverage、read 序列與其他變異證據。

七、BAM tags：HP、MM/ML tag 與 somatic haplotype

BAM tag	它描述什麼	本模組的判讀邊界
`HPHP tagBAM 裡標示某條 read 屬於哪一條 haplotype 的 tag。LongPhase-S 的 germline haplotag 寫成整數 HP:i:1；somatic haplotag 與 LongPhase-TO 則寫成字串 HP:Z:1-1。BAM tag assigning a read to a haplotype. Integer HP:i:1 for germline haplotag; string HP:Z:1-1 for somatic.完整條目 →`	工具對這條 read 的 haplotype assignment，例如 `HP:i:1` 或其他工具版本的 HP 值	HP1／HP2 是工具標籤，不自動代表父系／母系；沒有 HP 不代表 read 支持相反 allele。
`HP3`	在 somatic haplotype 分層中，帶有 somatic signal，但工具尚未可靠判定它源自哪一條 germline haplotype	HP3 不是「沒有 HP」。一般 untagged read 只是沒有 HP assignment；只有同時有 somatic signal 且來源未定，才使用 HP3 這個概念。
`MM/MLMM/ML tagBAM 中儲存每條 read 甲基化判讀的兩個 tag：MM 記錄位置、ML 記錄機率。在轉換或重新比對流程中可能遺失，需確認 tags 是否保留。BAM tags carrying per-read methylation calls: MM for positions, ML for probabilities. Easily lost during processing.完整條目 →`	`MM` 記錄 read 上修飾的位置與類型；`ML` 提供相應的修飾機率／信心值	這些 tags 可能在 BAM→FASTQ→BAM 的轉換中遺失；看到欄位空白時先檢查流程，不要直接解讀成「沒有甲基化」。

這裡的 HP3 是一個「帶有 somatic signal、但來源未定」的工具分類；它不等於普通的未標籤 read，也不等於已經證明存在一個獨立的細胞族群。

真實資料與證據

把命令和它讀寫的格式對上

下面不是一串需要盲目複製的指令，而是一張「哪個步驟讀哪個檔案、產生哪個檔案」的對照表。執行前仍需依平台、版本與路徑調整。

# 1. 為參考 FASTA 建立索引；產生 reference.fasta.fai
samtools faidx reference.fasta

# 2. FASTQ + FASTA → SAM；--MD 讓輸出帶有 MD tag（若下游需要）
minimap2 --MD -ax map-ont -t 10 reference.fasta reads.fastq -o alignment.sam

# 3. SAM → 排序後 BAM；再為 BAM 建立 .bai 索引
samtools sort -@ 10 alignment.sam -o alignment.bam
samtools index -@ 10 alignment.bam

.fai 與 .bai 是讓工具快速依座標取資料的索引，不是新的 read 或新的變異結論。比對時使用的 reference FASTA 版本也要和後續查詢座標一致。

把 VCF 位點回溯到 BAM read

假設 VCF 有一列 chr7:105 A>G。先取出它的座標，再用同一個座標區間查 BAM；這個動作把「site-level 的候選」連回「read-level 的證據」。

# 1. 先看 VCF 的非 header rows
bcftools view -H sample.vcf | head -3

# 2. 回到相同座標，看每一條 BAM alignment record
samtools view alignment.bam chr7:105-105

# 3. 用 reference FASTA 產生該位點的 pileup
samtools mpileup -f reference.fasta -r chr7:105-105 alignment.bam

# 4. 把 read 名稱、位置、MAPQ、CIGAR 與 tags 一起留下
samtools view alignment.bam chr7:105-105 | cut -f 1,3,4,5,6,12-

若 VCF 與 BAM 使用不同的 reference build，座標可能看似合法卻查到錯的區域。若 CIGAR 型式、MAPQ 或 pileup 顯示不一致，先檢查 alignment 品質、重複序列與 indel 對齊歧義，再判斷生物學意義。

轉換 BAM 時保留 MM/ML

如果原始 BAM 有甲基化 tags，轉成 FASTQ 再重比對時要明確把 tags 帶過去；否則輸出的 BAM 可能格式正常，卻沒有可用的 MM/ML。

# -T '*' 把所有 tags 帶進 FASTQ header/comment
samtools fastq -T '*' methylcall.raw.bam > methylcall.raw.fastq

# -y 讓 minimap2 把 FASTQ header 裡的 tags 寫回 BAM
minimap2 -ax map-ont -y reference.fasta methylcall.raw.fastq

重比對後請直接檢查 BAM 的 tags 是否仍存在。若漏掉 -T '*' 或 -y，流程可能不報錯，但後續需要 MM/ML 的分析會拿不到原始修飾資訊。

預測與結果檢視

你看到 VCF 的一列：chr7 105 A>G PASS GT:PS 0|1:105。如果你想知道這個候選是否真的有 read 支持，應先回哪一種檔案？

同時，這列中的 0|1 與 PS=105 各自告訴你什麼？

展開答案

先回 BAM。用 VCF 的 CHROM 與 POS 查詢該區間，再查看每條 read 的序列、MAPQ、CIGAR 與可用 tags；需要時再對照 reference FASTA。FASTQ 是尚未放回座標的原始 read，不能直接回答哪條 read 覆蓋 chr7:105。

0|1 表示 REF 與 ALT 的順序已被定相表示；PS=105 表示這個位點屬於編號為 105 的 phase set。它們都不等於支持 read 數量，也不會取代對 BAM 的檢查；HP1／HP2 也不因此自動帶有父系／母系意義。

另一條 BAM record 的 CIGAR 是 60S40M。這條 read 總長多少？它消耗了多少參考序列？

展開答案

read 長 100 bp，參考涵蓋 40 bp。S 仍保留在 read，所以 $60 + 40 = 100$ ；但 S 不消耗參考，只有 40M 消耗 40 個參考位置。

如果許多 read 在相近座標出現相同方向的 S，可把該處列為結構變異斷點候選；仍需結合其他證據，不能由單一 CIGAR 直接下結論。

實作練習

練習一：從一列 VCF 開始追證據

選一列 VCF，先寫下 CHROM、POS、REF、ALT，再逐步回答：

它是一個 site-level row，還是一條 read？
要回 BAM 查哪一個座標區間？
在 BAM record 中，哪些固定欄位告訴你 read 的位置、對齊型式與品質？
有哪些 HP、MM/ML tags？沒有 HP 時，能不能直接寫成 HP3？

最後用 pileup 或 IGV 對照支持與反對的 read；不要只看 VCF 的 FILTER=PASS 就跳過 read-level 檢查。

練習二：讀 CIGAR 的兩個總和

對每個 CIGAR，分別把「消耗 read」與「消耗參考」的操作碼加總，再用上方解碼器核對：

CIGAR	read 長度	參考涵蓋長度	先看哪個規則？
`10M`	10	10	`M` 兩邊都消耗
`5M2I5M`	12	10	`I` 只消耗 read
`5M2D5M`	10	12	`D` 只消耗參考
`3S7M`	10	7	`S` 留在 read、不佔參考

如果你只記得一句話：先辨認檔案的資料粒度，再讀欄位；CIGAR 是 BAM 對齊欄位的解碼規則，不是 NGS 檔案鏈的起點。

學習檢核

本模組術語

BAM: 已比對到參考基因體的 read 的二進位格式（SAM 的壓縮版）。每一列是一條 read，帶著它的位置、CIGAR、品質與各種 tag。
CIGAR: 描述一條 read 如何對上參考基因體的緊湊字串。例如 5M1I5M 表示兩端各有 5 個對齊欄位，中間有 1 個相對於參考序列的插入；M 可能代表吻合，也可能代表錯配。
FASTA: 以文字保存 DNA 或 RNA 序列的格式；在本模組中，FASTA 主要用來提供比對所需的參考基因體序列。
FASTQ: 存放 read 序列與每個鹼基品質分數的文字格式；每筆紀錄由四行組成。
HP tag: BAM 裡標示某條 read 屬於哪一條 haplotype 的 tag。LongPhase-S 的 germline haplotag 寫成整數 HP:i:1；somatic haplotag 與 LongPhase-TO 則寫成字串 HP:Z:1-1。
MM/ML tag: BAM 中儲存每條 read 甲基化判讀的兩個 tag：MM 記錄位置、ML 記錄機率。在轉換或重新比對流程中可能遺失，需確認 tags 是否保留。
PS tag（phase set）: 標示某群 variants 屬於同一 phase block 的編號。同一 PS 內的相位方向可相互比較；不同 PS 的 HP1/HP2 方向彼此獨立。
VCF: 存放 variant 的文字格式。每一列是一個位點，記錄位置、ref/alt allele、品質、FILTER 與各種 FORMAT 欄位。
reference genome（參考基因體）: 用於比對的標準序列（人類資料常用 GRCh38）。所有座標均相對於指定版本；更換 build 會改變座標系，座標不可直接比較。
soft clipping: read 端部未與參考序列對齊，但序列仍保留在 BAM 中（CIGAR 為 S）。大量 clipping 集中於同一位置時，可能提示結構變異斷點，仍需其他證據確認。
somatic haplotype: 在 germline haplotype 之下再細分出來的層級。某條 HP1 的 read 若帶有 somatic 突變，就屬於 HP1-1；HP2 的則是 HP2-1。無法判定來源的歸為 HP3。