先建立最小的生物模型

本模組先從基本概念開始，暫不引入癌症、EGFR 或其他縮寫。理解下列五個層次的關係後， 即可建立後續模組所需的資料模型：

人體細胞核內含有染色體。在本教材的簡化模型中，每條染色體可視為一條長 DNA 分子； DNA 以 A、C、G、T 四種鹼基記錄序列資訊。多數正常常染色體具有父源與母源兩份拷貝； 性染色體與腫瘤拷貝數異常將於後續模組另行說明。

基因不是另一種物質，而是 DNA 序列中具有特定功能的區段；可將其類比為程式系統中的模組。 當其中某個位置相對參考序列發生差異，稱為變異；「突變」則常用於描述新取得的變化。 部分變異可能影響蛋白質功能或治療反應，仍需依變異與情境判定。

「定序」到底做了什麼

我們無法把一整條染色體直接貼進文字編輯器。實際流程是：

1. 從檢體取出大量 DNA 分子。
2. 定序儀逐段量測 DNA，把每一段轉成 A／C／G／T 字串。
3. 每一段序列觀測稱為一條 read；依平台與定序深度，一個樣本可產生大量且部分重疊的 reads。
4. 電腦把 reads 對回人類參考序列，找出不同的字母位置。

參考序列：... A C G T T A C G ...
某條 read：... A C A T T A C G ...
                     ↑
                 這裡有一筆差異

所以定序資料不是「病人的完整基因組序列」，而是很多條來自不同實體 DNA 分子的局部觀測。 分析流程先將它們比對至參考序列，再由 variant-calling 步驟產生候選變異清單。

與治療決策相關的示例

EGFR 是一個編碼細胞生長訊號相關蛋白質的基因。 部分帶有活化 EGFR 變異的肺癌可受益於特定 EGFR 抑制劑；另一些變異則可能降低特定藥物的抑制效果。 實際治療判定仍需依藥物、變異與臨床情境評估。

一位非小細胞肺癌病人在治療後，檢測報告同時找到 EGFR 上的 T790M 與 C797S。這兩個名字只是「EGFR 的哪個位置、變成什麼」的標籤； 現在不用記。重要的是：報告只告訴我們兩個突變都存在，卻沒有回答它們是否在同一份 EGFR DNA 上。

報告已知：detected("T790M") = true
          detected("C797S") = true

還不知道：same_DNA_molecule("T790M", "C797S") = ?

下面每一條橫線是一條 read。實心點表示該 read 帶有那個突變， 空心點表示該位置與參考序列相同。先只根據圖上的分子證據作答：

這個案例來自 Niederst 等人的原始研究；後續也有 依排列方式選藥的臨床案例。 這裡用它說明「分子關係會改變決策」，不是提供個別病人的治療建議。

同樣找到兩個突變，read 圖可能給三種答案

前面的 EGFR 報告只說「T790M 有、C797S 也有」。要知道兩者的排列， 需回到產生這張清單的原始 reads。在本簡化情境中，可能的判讀可分為下面三種：

判讀順序先看範圍，再看實心點。如果同一條 read 根本沒有機會同時看見兩個位置， 目前只能判定資訊不足；只有在 read 確實跨過兩點且品質足夠時，點的組合才可支持「同一份」或「分開」的解釋。

長 read 通常可增加同一次觀測涵蓋的距離，但實際可連結範圍仍取決於讀長分布、覆蓋度與品質。 上面的 EGFR 兩點很近，短片段也可能回答；實際分析也可能處理相距數千至數萬 bp 的變異。 定序技術之間的長度與準確度取捨，留到定序單元完整比較。

兩位點問題如何延伸為三種分析情境

前述圖示只包含兩個突變，且假設 reads 品質良好。真實研究處理相同類型的問題， 但資料會逐層變難：先擴大到整個基因體，再加入腫瘤混合，最後還可能拿不到同一病人的正常檢體。 三項工具分別處理這三個層級的分析條件。

這三項工具共同說明如何將 read-level 判讀問題， 擴展為可處理人類基因體與腫瘤資料的分析方法。M11–M13 依相同層級排列： M11 是一般樣本，M12 加入腫瘤與正常對照，M13 再拿掉正常對照。

未觀察到證據，不等於證明不存在

read 1, 2, 3  → 只覆蓋左邊，看到 T790M
read 4, 5, 6  → 只覆蓋右邊，看到 C797S

接下來的路線

此時不需先安裝軟體或另讀生物學教科書。請依序學習，各模組會在需要時引入新的概念與資料物件：

本頁核心重點如下：

為什麼重要

本教材採用同一個簡化模型連結多項概念。本模組先完整建立細胞、DNA 分子與 read 的基本關係， 並在專有名詞首次出現時提供本頁所需的簡要定義。

本模型聚焦於一個核心問題：腫瘤檢體經過標準 bulk 定序後，資料保留了哪些資訊， 又不再包含哪些細胞層級資訊？

概念與互動

基本資料模型

教學用簡化區段模型

圖中元素的意義如下：

三個變異代表不同的細胞譜系階段

三個變異在細胞群中的分布反映其出現的先後與分支關係：

• S1 先出現，因此紅、藍、橘三個腫瘤細胞群皆帶有 S1。
• S2 出現於後續分支，僅存在於藍色細胞群。
• S3 出現於另一後續分支，僅存在於橘色細胞群。

本模型另設定 S1、S2、S3 均位於 HP1。在此簡化情境中，每個新增變異最初發生於其中一份 DNA 拷貝，因此可以進一步分析變異與 HP1 或 HP2 的分子連結。此處的 HP1／HP2 是方向任意的相對標籤， 不代表父系或母系。

由細胞比例推算 VAF

VAF 是支持特定變異的 read 比例。在本模型的假設下，可由細胞比例推算其預期值。 以 S1 為例：

• 帶有 S1 的細胞比例 = 紅色群（20%）+ 藍色群（20%）+ 橘色群（10%）= 50%
• 本模型假設每個細胞在此位置有兩份 DNA 拷貝，而 S1 僅位於其中一份（HP1）
• 因此，支持 S1 的預期 read 比例為 50% × ½ = 25%

因此圖中標示 S1 的預期 VAF 為 0.25。依相同假設，S2 與 S3 的預期 VAF 分別為 0.10 與 0.05。

bulk 定序保留與未保留的資訊

read 內的分子連結

在標準 bulk 定序資料中，read 與來源細胞的對應關係不會被保留；輸出是一組來自多個細胞的 DNA 片段序列。因此，無法僅由這些資料直接辨識每條 read 原先所屬的細胞群。

然而，read 內仍保留分子層級的連結資訊。若某條 read 同時涵蓋兩個變異， 該觀測可直接支持兩個變異位於同一條 DNA 分子。

這項分子連結構成本教材後續分析方法的共同證據基礎：

與 single-cell 定序比較

single-cell 定序會先分離個別細胞再進行定序，因此可追溯 read 的細胞來源， 並直接比較不同細胞的變異組合。

相較之下，single-cell 定序通常需要在成本、細胞數與單細胞覆蓋度之間取捨。 本教材聚焦於 bulk 定序，因此細胞群結構屬於根據分子連結建立的推論模型，而非直接的細胞層級觀測。 後續分析均需維持此證據層級的區分。

判讀練習

練習與自我檢核

VAF 計算練習

依據上圖的細胞比例與本模組的簡化假設，計算三個變異的預期 VAF：

VAF 推論的限制

若三個變異的 VAF 分別為 0.30、0.30、0.08，是否可以據此推論 「前兩個變異位於同一細胞群，且第三個變異較晚出現」？

不可以。VAF 是各位置獨立統計的 read 比例，不包含變異共同出現或發生順序的資訊。 因此，相同 VAF 不代表兩個變異位於同一條分子，較低 VAF 也不能單獨證明變異較晚出現。 可支持的結論僅為：「第三個變異的相對 read 支持較少；其細胞分布與發生順序仍需分子連結及其他證據判定。」

分析與報告時應明確區分三個層級：直接觀測、資料支持的模型推論，以及目前證據無法支持的結論。

為什麼重要

M1 先以 HP1／HP2、founder clone 與 subclone 作為圖示標籤，以呈現整體資料流。 本模組將進一步區分這些標籤，並建立較精確的生物學資料模型。

這些術語形式相近，但所指的資料層級不同。常見混淆是將 variant 與 allele、 或 genotype 與 haplotype 混用，因而難以精確解讀文獻與資料格式。

自學閱讀時，可先將下列五個詞對應至資料結構的不同層級： 地址、該地址的一種值、與參考不同的事件、單一地址的兩份值、跨多個地址的一整串值。 先掌握概念層級，再對應英文術語。

概念與互動

五個容易混用的詞

上表以正常二倍體常染色體為模型，因此同一個 locus 通常有兩份 allele。 這種具有兩份拷貝的狀態稱為 diploiddiploid 二倍體正常人多數常染色體區段通常有兩份拷貝的狀態；性染色體、拷貝數變化與腫瘤細胞可能不同。A state in which a typical autosomal locus has two copies. Sex chromosomes, copy-number-altered regions and tumour cells may differ.完整條目 →；性染色體與腫瘤拷貝數異常將於後續模組另行處理。

關鍵差別在最後兩個：genotype 是單一位置的，haplotype 是跨越多個位置的。 genotype 說「這裡有一份 A、一份 G」，但沒說 A 跟隔壁位置的 C 是不是在同一條染色體上。 補上這個資訊，正是 phasingphasing 定相為可判定的 variants 建立其位於不同實體染色體拷貝上的相位關係；VCF 常以 0|1 等形式表示。長 read 可提供跨位點的分子層級觀測證據。Determining which variants lie on the same physical chromosome, turning 0/1 into 0|1. Long reads provide direct molecular evidence.完整條目 → 要推定的關係之一。

VCF 可用下列符號表示這項差異：

phasing 會以 reads 或統計模型推定等位基因在不同染色體拷貝上的歸屬； VCF 可使用豎線表示已推定的相位結果。

heterozygous 位點提供主要 phasing 資訊

一個 locus 若兩份 allele 不同（heterozygousheterozygous 異型合子在一般二倍體位點上，兩份 allele 不同（例如一份 A、一份 G）。此類位點可作為區分 haplotype 的資訊錨點；複雜拷貝數情況需另行解讀。At a typical diploid locus, carrying two different allele copies (for example A and G). These sites can anchor phasing; copy-number changes require separate interpretation.完整條目 →，例如 0/1）， 在讀取品質足夠的簡化模型中，可用來區分兩條 haplotype —— 看到 A 的 read 可歸到一條，看到 G 的可歸到另一條。

若兩份相同（homozygoushomozygous 同型合子在一般二倍體位點上，兩份 allele 相同。此位點本身通常無法用來區分兩條 haplotype。At a typical diploid locus, carrying two identical allele copies. The site itself usually cannot distinguish haplotypes.完整條目 →），該位點本身通常不能區分兩條 haplotype。

het 位點密度是影響 phasing 的因素之一；讀長、coverage 與錯誤率也會影響可建立的相位範圍。

腫瘤中的 LOHLOH 異型合子性喪失原本 heterozygous 的區域變成只剩一種 allele。LOH 不等於缺失 —— 也可能是一條 haplotype 遺失後另一條被複製（copy-neutral LOH）。Loss of heterozygosity: a formerly het region retains only one allele. Not necessarily a deletion.完整條目 →（loss of heterozygosity，異型合子性喪失）可能使一個區段失去其中一種 allele， 使原本可用的 heterozygous 位點減少，因而降低可用的 phasing 錨點。

因此，僅依這些錨點的 phasing 在 LOH 區段會更困難，可能需要其他證據或專門方法； 實際影響取決於 LOH 型態、coverage、purity 與演算法。

germline 與 somatic

在特定 somatic calling 流程中，germline calls 常被排除；然而 heterozygous germline SNPSNP 單核苷酸多型性族群中以多型形式存在的單一鹼基差異；本教材主要以可作為 phasing 錨點的 germline SNP 為例。A single-nucleotide polymorphism observed as a population variant. This tutorial mainly uses germline SNPs as phasing anchors.完整條目 → 也可作為建立 haplotype 骨架的錨點之一。錨點不足時，區分 HP1／HP2 及連結 somatic 變異會更困難。

clone 與 subclone

回到 M1 的簡化區段模型。一個細胞取得突變後增生，形成一群共享可辨識突變的細胞，稱為 cloneclone 克隆源自共同祖先細胞，並共享一組可辨識 somatic mutations 的細胞群。A population of tumour cells descended from a common ancestor and sharing the same set of somatic mutations.完整條目 →。 其中某些細胞又取得新突變、再擴增，就形成 subclonesubclone 次克隆clone 內取得額外變異並擴增的子群；治療後可能富集具有抗性的 subclone。A subset of a clone that acquired further mutations and expanded. Treatment-resistant populations are often subclones.完整條目 →。

許多 somatic 變異可能是 passenger；可提高生長或存活優勢的候選變異稱為 driver mutationdriver mutation 驅動突變可提升腫瘤細胞生長或存活優勢的候選變異；passenger mutation 通常不具有已知的選擇優勢。A mutation that confers a growth advantage, as opposed to the far more numerous passenger mutations.完整條目 →， 其判定需依癌別與功能證據。

治療可能改變不同細胞群的相對比例，使具有抗藥優勢的 subclone 富集：

因此，除列出腫瘤中的變異外，也需評估這些變異可能分布於哪些細胞群。 這種細胞組成的不均勻，就叫 intratumour heterogeneityintratumour heterogeneity 腫瘤內異質性同一腫瘤內不同細胞可帶有不同突變組合，因而同時存在多種基因型。Different cells within one tumour carrying different mutation combinations.完整條目 →。

真實證據

評估方法時需有可供對照的高可信度參考標註。 truth set 通常由多種技術與證據建立，但仍可能包含不確定區域；benchmark 還包含評估範圍與程序。

以下列出六個具有參考標註的癌症細胞株作為示例：

這些細胞株具有可用的配對正常樣本與參考標註，可提供區分 germline 與 somatic 的重要對照； 實際判定仍受 coverage、污染、變異類型與參考集品質影響。

HCC1395 將於後續示例中再次出現；可比較兩份資料： HCC1395_HKU 與 HCC1395_NYGC，同一個細胞株、不同的 basecallerbasecalling 鹼基判讀把定序儀的原始訊號（ONT 是電流）轉成 A/C/G/T 字母的步驟。不同 basecaller 版本會產生不同的錯誤特性 —— 所以「同一個細胞株、不同 basecaller」是兩份不同的資料集。Converting raw sequencer signal into A/C/G/T. Different basecallers yield different error profiles.完整條目 → 與分析流程。 兩份資料的組成相似度（0.909）可作為跨流程一致性的觀察，但不能單獨證明方法穩健性。 因此，provenanceprovenance 資料來源履歷記錄資料來源、細胞株、定序平台、basecaller、reference build、caller、benchmark 與混樣方式。在本實驗室中，這些資訊是重要實驗變數，而非附註。The full origin record of a dataset. In this lab it is a first-class experimental variable, not metadata.完整條目 → 應視為重要實驗變數，而非附註。

判讀練習

練習與自我檢核

三項概念檢核

可先自行回答，再核對下方解析：

1. 「這個位置有一份 A、一份 G」—— 這句話講的是 genotype 還是 haplotype？
2. 0|1 比 0/1 多告訴你什麼？
3. 為什麼 heterozygous 位點通常比 homozygous 位點更能提供 phasing 資訊？

答案依序是：genotype（單一位置）／多告訴你「哪一份在哪一條染色體上」／ 因為兩份 allele 不同時，才提供區分兩條染色體的訊號；homozygous 位點通常無法單獨區分。

在實際資料中檢查

若已有一份 VCF，可使用下列指令檢查 genotype 分布：

# 各種 genotype 各有多少個
bcftools query -f '%POS\t[%GT]\n' sample.vcf | cut -f2 | sort | uniq -c | sort -rn

# 只留下 heterozygous 的位點（phasing 能用的錨點）
bcftools view -g het sample.vcf | grep -vc '^#'

het 位點密度會影響 phasing，亦需同時考慮讀長、coverage 與品質。 人類全基因體中的 het 位點通常為數百萬量級，但會因族群、參考基因組、callable region 與過濾條件而異；請以資料集統計為準。

為什麼重要

M1 已建立以下觀察：若 read 同時涵蓋兩個變異位置，且判讀可靠， 可支持兩者位於同一 DNA 分子。本模組將此觀察連結至定序技術，說明讀長差異、長讀的分析優勢與相應限制。

自學閱讀時，可先比較短讀通常能辨識各位置的變異、但較難由單條 read 建立遠距分子關係， 以及長讀可涵蓋較遠位置的特性；再學習平台名稱與錯誤類型。

瞭解這些限制很重要；後續介紹的 somatic-indel 方法， 即聚焦於降低此類錯誤並改善候選判讀。

概念與互動

長讀long read 長讀單條可達數千至數萬鹼基的定序片段（例如 ONT、PacBio）。若可靠地同時覆蓋多個 variant，可提供它們位於同一 DNA 分子上的直接觀測證據。A sequencing read thousands to tens of thousands of bases long. Its length lets one molecule link multiple variants directly.完整條目 →提供的長距離資訊

此圖說明短讀與長讀的主要資訊差異。上半部並非表示短讀無法偵測變異； 在本比較條件下，Illumina 可分別辨識兩個位點，且單鹼基錯誤率通常較低， 但不能由單條 read 建立 3 kb 的分子連結。實際效能依平台、化學版本、basecaller 與變異類型而異。

長讀 indel 判讀的主要限制

ONT 與 Illumina 的錯誤特性不同。homopolymerhomopolymer 同聚物同一個鹼基連續重複的區段，例如 AAAAAA。nanopore 在此類區段較容易發生長度判讀錯誤，是 indel 假陽性的常見來源之一。A run of identical bases. Nanopore miscounts their length, a major source of false indel calls.完整條目 → 是 ONT indel 錯誤的常見高風險區域， 例如 AAAAAA。

原因在判讀原理上：

下表提供效能示例，用於比較 SNV 與 indel 的相對差異。 數值解讀需同時記錄 caller 名稱與版本、資料集、truth set 與指標定義：

在此示例中，indel F1 低於 SNV F1，表示 somatic-indel 的整體偵測效能仍有改善空間； 後續方法會利用更多 read-level 特徵處理此問題。

三項核心數值

以下以簡化數值示例說明低腫瘤比例如何減少可用訊號：

在二倍體、單拷貝、clonal、無偏取樣且 tumor DNA fraction 為 0.2 的簡化模型中， 期望支持數為 $50\times 0.2÷2=5$ 條。訊號數量不足是主要限制之一；背景錯誤、品質與演算法也會影響可靠性。

basecaller 是實驗變數

basecallingbasecalling 鹼基判讀把定序儀的原始訊號（ONT 是電流）轉成 A/C/G/T 字母的步驟。不同 basecaller 版本會產生不同的錯誤特性 —— 所以「同一個細胞株、不同 basecaller」是兩份不同的資料集。Converting raw sequencer signal into A/C/G/T. Different basecallers yield different error profiles.完整條目 → 是把電流訊號轉成 A/C/G/T 的那一步。不同版本的 basecaller 會產生不同的錯誤特性，所以：

真實證據

合成不同 purity 的樣本

若要評估方法在不同 tumor DNA fraction 下的效能，可依計算後的比例混合 tumor 與 normal BAM， 建立預先設定的合成梯度：

# 1. 先量出兩個 BAM 的實際 coverage（本例量到 tumor 50×、normal 25×）
mosdepth --by 500 tumor_cov  tumor.bam
mosdepth --by 500 normal_cov normal.bam

# 2. 抽樣率不是目標 fraction，要從實測深度回推：
#      抽樣率 = （目標 fraction × 合成後總深度）÷ 來源深度
#    總深度取 25×（兩邊都拿得出來），目標 DNA fraction 0.6：
#      tumor  需要 0.6×25 = 15× → 15 ÷ 50 = 0.30
#      normal 需要 0.4×25 = 10× → 10 ÷ 25 = 0.40
samtools view -s 1.30 -b tumor.bam  > tumor_sub.bam
samtools view -s 1.40 -b normal.bam > normal_sub.bam

# 3. 合併成指定 DNA fraction 的樣本
samtools merge dna_frac_0.6.bam tumor_sub.bam normal_sub.bam
samtools index dna_frac_0.6.bam

-s 的小數部分是抽樣率，不是你要的 fraction。 兩個來源深度不同時直接填 0.60／0.40 會得到 30×:10× ＝ 0.75，不是 0.6。 以 tumor 50×、normal 25×、合成後總深度 25× 為例，五個梯度的抽樣率是：

這樣控制的是 tumor DNA/read fraction，五個樣本的總深度也才一致 —— 否則深度與 fraction 兩個變數會同時在動，之後看到效能下降就分不清是哪一個造成的。

判讀練習

實作練習

檢查資料讀長

讀長會影響可連結的變異距離，仍需同時考慮 coverage、品質與演算法。取得新資料時可先檢查：

# 讀長的統計摘要
samtools stats aln.bam | grep -E '^SN' | grep -E 'length|reads mapped'

# 讀長分布（看有多少條真的夠長）
# -F 0x900 排除 secondary 與 supplementary，否則同一條 read 會被數好幾次
samtools view -F 0x900 aln.bam | awk '{print length($10)}' | sort -n | \
  awk '{a[NR]=$1} END {print "中位數", a[int(NR*0.5)]; print "第 90 百分位", a[int(NR*0.9)]; print "最長", a[NR]}'

# 覆蓋度
mosdepth --by 500 out aln.bam && head out.mosdepth.summary.txt

估算簡化模型中的預期支持數

可使用資料的 coverage，估算在二倍體、單拷貝、clonal 與無偏取樣假設下的預期支持 read 數； 這不是經驗證的 detection limit：

若期望值僅為兩至三條，可靠偵測低頻變異可能較困難。資料量是限制之一， 仍需比較增加定序深度、改善品質與調整模型的相對效益。

為什麼重要

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

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

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

概念與互動

一、先建立檔案地圖

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。

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

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

四、先讀 VCF 的關鍵欄位

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」。互動元件只做這一個局部任務；它不是另一種檔案格式。

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

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

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

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

這裡的 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 開始追證據

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

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

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

練習二：讀 CIGAR 的兩個總和

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

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

為什麼重要

本課程後續工具多以 caller 輸出的候選變異為輸入，進一步進行篩選或校正。 因此需要先理解 caller 的輸出內容與限制。

germline calling 的情境相對單純：一個人只有兩條染色體，變異要嘛在其中一條、要嘛兩條都有。 本節在這個單純的情境下把 caller 的判斷方式講清楚，M07 再處理腫瘤帶來的額外因素。

概念與互動

變異檢測variant calling從已比對 read 中判斷哪些位置具有足夠證據支持變異存在的過程，通常輸出 VCF。Deciding which positions genuinely carry variants, given aligned reads. Outputs a VCF.完整條目 →的判斷內容

某個位置有 30 條 read 支持 A、28 條支持 G。這組數字可以提出 heterozygous 候選， 但不能單獨定案 —— 因為同樣的 30 比 28，畫面可以完全不同：

所以 caller 讀的不是兩個數字，而是這個位置的整張畫面 —— 把覆蓋此處的 read 排成一張表（一列一條 read、一欄一個位置，這種排法叫 pileup）， 再把表上每個格子的資訊編碼成好幾層數值送進模型。Clair3Clair3以深度學習做 germline variant calling 的工具，長 read 上常用。A deep-learning germline variant caller widely used for long reads.完整條目 → 與 DeepVariantDeepVariantGoogle 開發的 germline variant caller，把 pileup 轉成影像再用 CNN 分類。Google's germline variant caller; encodes pileups as images and classifies them with a CNN.完整條目 → 都是這種做法：

兩個工具的內部步驟不同：DeepVariant 把「做數值」「跑模型」「寫 VCF」拆成三個程式； Clair3 則先用輕量模型掃過所有位置，只把信心不足的少數位置交給看完整比對的模型重判。

把 read 先分成兩群，判斷就變簡單了

前面那張畫面還有一個沒用上的線索：這些 read 分別來自哪一條染色體。 一個人的每個位置有兩條染色體，一條來自父親、一條來自母親。 把 read 依這兩條分成兩群（也就是 phasingphasing 定相為可判定的 variants 建立其位於不同實體染色體拷貝上的相位關係；VCF 常以 0|1 等形式表示。長 read 可提供跨位點的分子層級觀測證據。Determining which variants lie on the same physical chromosome, turning 0/1 into 0|1. Long reads provide direct molecular evidence.完整條目 →），同一個畫面會變得好判斷得多：

關鍵在一致性，不在數量。真的 heterozygous 變異只存在於其中一條染色體上， 所以那一群的每一條 read 都應該帶著它；隨機的定序錯誤沒有這個約束，會零星散落在兩群裡。 這個判準對錯誤率較高的平台（例如部分 ONT 資料）特別有用 —— 那正是單看計數最不可靠的情形。

這也解釋了一個看起來矛盾的地方：明明要先有變異才能 phasing，怎麼又說 phasing 幫助 calling？ 因為兩者是來回兩次：先用最有把握的位點做一次 calling，用它們把 read 分群， 再把 HP 標籤HP 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.完整條目 →當成模型的額外輸入重判一次。Clair3 的 full-alignment 階段與 PEPPER-Margin-DeepVariant 都採用這個順序。

同一個「集中或分散」的線索在腫瘤資料裡會再出現一次，只是那時要判斷的不只是真假， 還包括這個變異長在哪一條 haplotypehaplotype 單倍型同一條實體染色體拷貝上，具有一致相位的一組 alleles 或 variants。HP1 與 HP2 是任意的相對標籤。A set of variants that lie on the same physical chromosome copy and are inherited together.完整條目 → 上 —— 需要的證據也更多。以下這張圖先把最重要的陷阱說清楚： 兩種情況的支持 read 總數可以幾乎一樣，所以不能靠數量判斷。

真實證據

本教材示例使用 Clair3Clair3以深度學習做 germline variant calling 的工具，長 read 上常用。A deep-learning germline variant caller widely used for long reads.完整條目 →；部分流程也會使用 PEPPER-Margin-DeepVariant。 以下提供一組示範指令：

docker run \
  -v "${INPUT_DIR}":"${INPUT_DIR}" \
  -v "${OUTPUT_DIR}":"${OUTPUT_DIR}" \
  kishwars/pepper_deepvariant:r0.7 \
  run_pepper_margin_deepvariant call_variant \
  -b "${INPUT_DIR}/alignment.bam" \
  -f "${INPUT_DIR}/reference.fasta" \
  -o "${OUTPUT_DIR}" \
  -t 10 \
  --ont_r9_guppy5_sup

最後的 --ont_r9_guppy5_sup 是與特定 ONT 化學版本／basecaller 資料相容的模型選項。 R10 與 PacBio HiFi 應依所用版本文件選擇相容模型（例如 --ont_r10_q20 或 --hifi）。 模型不匹配可能降低效能，且流程未必報錯；請將平台、basecaller 與模型版本記錄於 provenanceprovenance 資料來源履歷記錄資料來源、細胞株、定序平台、basecaller、reference build、caller、benchmark 與混樣方式。在本實驗室中，這些資訊是重要實驗變數，而非附註。The full origin record of a dataset. In this lab it is a first-class experimental variable, not metadata.完整條目 →。

在本教材流程中，germline VCF 產生後，再由後續 phasing 工具使用其中的 heterozygous 位點作為錨點。 這是本流程的先 call、再 phase 順序；其他工具可能採用不同的整合方式。

判讀練習

實作練習

執行一次 germline calling

# Clair3（Docker 範例；請依環境調整路徑）
docker run -v "${INPUT_DIR}":"${INPUT_DIR}" -v "${OUTPUT_DIR}":"${OUTPUT_DIR}" \
  hkubal/clair3:latest \
  /opt/bin/run_clair3.sh \
  --bam_fn="${INPUT_DIR}/aln.bam" \
  --ref_fn="${INPUT_DIR}/reference.fasta" \
  --threads=10 \
  --platform="ont" \
  --model_path="/opt/models/ont" \
  --output="${OUTPUT_DIR}"

模型應與定序化學版本、basecaller 與工具版本相容。模型不匹配可能降低效能， 且未必觸發錯誤；請在執行前核對模型文件，並將設定記錄於 provenance。

檢查結果的合理性

取得 VCF 後，先進行以下基本合理性檢查：

# 變異總數；典型高品質人類全基因體 callset 約為數百萬筆，實際值依篩選條件而異
bcftools view -H out.vcf | wc -l

# het 與 hom 的比例；實際值依族群、平台與篩選條件而異
bcftools view -H -g het out.vcf | wc -l
bcftools view -H -g hom out.vcf | wc -l

# SNV 與 indel 的比例；實際值依平台與篩選條件而異
bcftools view -H -v snps out.vcf | wc -l
bcftools view -H -v indels out.vcf | wc -l

這些統計若明顯偏離同類資料的預期，可能提示 coverage、模型相容性、參數、過濾條件或樣本品質問題。 請先檢查原因，再進行 phasing，以降低問題在後續流程中傳遞的風險。

原始文獻與程式碼

• Clair3：Zheng Z 等，Symphonizing pileup and full-alignment for deep learning-based long-read variant calling，Nature Computational Science，2022。程式碼在 github.com/HKU-BAL/Clair3。
• DeepVariant：Poplin R 等，A universal SNP and small-indel variant caller using deep neural networks，Nature Biotechnology，2018（doi 10.1038/nbt.4235）。程式碼在 github.com/google/deepvariant。

為什麼重要

M0 的那個問題到這裡才真正有答案。報告說 T790M 有、C797S 也有， 但要回答它們在不在同一份 DNA 上，得先把 read 依來源分成兩群 —— 這件事就是 phasing。

它產出兩個東西，之後每個模組都會用到：

• 變異的相位寫進 VCF —— 哪些變異在同一條 haplotype 上。
• 每一條 read 的 haplotype 寫進 BAM 的 HP tag —— 這條 read 來自哪一條。 這一步叫 haplotagging。

第二個特別關鍵：M7 判斷 somatic 變異真假時用的「支持 read 是否集中在同一條 haplotype」， 以及 M9 依 haplotype 分開統計甲基化，都是在 HP tag 上算的。 沒有這一步，那些判準根本無從計算。

本模組先用一段簡化 pileup 讓你自己分一次群，再說明工具怎麼做同一件事、 以及它在哪些情況下會分錯。

概念與互動

練習一次 phasing 判讀

下面是一段 pileup，包含三個 heterozygousheterozygous 異型合子在一般二倍體位點上，兩份 allele 不同（例如一份 A、一份 G）。此類位點可作為區分 haplotype 的資訊錨點；複雜拷貝數情況需另行解讀。At a typical diploid locus, carrying two different allele copies (for example A and G). These sites can anchor phasing; copy-number changes require separate interpretation.完整條目 → 位點；每條 read 在其覆蓋到的位點上顯示一個鹼基。 請依跨位點的 allele pattern 將 read 分群，使同群 read 的組合盡可能一致，並保留無法可靠判定的 read。

做完之後，注意三件事：

• 有一條 read 在此練習中無法可靠指派 —— 它沒有覆蓋任何 informative het 位點，因此保持未標記即可。 若該 read 另帶有 somatic signal，才可能在後續流程中進入 HP3 類別；一般 untagged read 不等於 HP3。
• 示例中有一個錯誤鹼基（檢查後會被圈出來）。該 read 在其他位點仍可能呈現一致 pattern， 因此可根據多個位點的整體證據進行分群，但此方法不保證修正所有錯誤。
• HP1 與 HP2 整體互換仍可表示相同的相對 phasing。原因見下方說明。

cis 與 trans

phasing 的產物就是回答「兩個變異的相對位置」：

• ciscis兩個 variants 在同一條 haplotype 上。Two variants on the same haplotype.完整條目 →：在同一條 haplotype 上；藥物或功能影響需依特定 EGFR 變異與研究證據判定
• transtrans兩個 variants 分別在兩條不同的 haplotype 上。Two variants on opposite haplotypes.完整條目 →：分屬兩條 haplotype

phase set 與 switch error

phasing 可能因缺乏重疊 read 或 informative 位點而中止，形成新的 phase blockphase block一段可建立連續相位關係的區域。read 長度不足、缺少 informative heterozygous 位點或證據不一致時，可能形成不同 phase blocks。A contiguous stretch over which phasing is consistent. It can break when reads are too short, informative heterozygous sites are absent or evidence conflicts.完整條目 →。 VCF 用 PS 欄位標示每個 block。

switch errorswitch errorphasing 結果自某一位置起將兩條 haplotype 的方向對調，是常見的 phasing 錯誤類型。A phasing error in which the two haplotype labels swap from some point onward; a common phasing error mode.完整條目 → 是常見的 phasing 錯誤類型：自某一位置起，演算法將兩條 haplotype 的標籤對調。 它可能發生在 het 位點稀疏、read 重疊不足或證據歧義較高的區域。

haplotagging：將 haplotype assignment 寫入 read

phasingphasing 定相為可判定的 variants 建立其位於不同實體染色體拷貝上的相位關係；VCF 常以 0|1 等形式表示。長 read 可提供跨位點的分子層級觀測證據。Determining which variants lie on the same physical chromosome, turning 0/1 into 0|1. Long reads provide direct molecular evidence.完整條目 → 將 variant 的 haplotype assignment 寫入 VCF； haplotagginghaplotagging根據已 phase 好的 variants，把每一條 read 指派到 HP1 或 HP2，並把結果寫回 BAM 的 HP tag。Assigning each read to HP1 or HP2 using phased variants, and writing the result back as a BAM HP tag.完整條目 → 則將每條 read 的 assignment 寫入 BAM 的 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.完整條目 → tag。 之後可在 IGV 中依 haplotype 分組或著色顯示。

somatic haplotype 的延伸層級

到目前為止介紹的是標準 germline phasing。接著再說明 somatic haplotype 的延伸層級：

若某個 somatic allele 在其他 read 上已有可靠的 HP1-1 assignment， 流程可在條件成立時將攜帶相同 allele 的部分 HP3 read 重新分類為 HP1-1；此步驟可能增加可用 read， 但也需評估相同 allele、mapping 與相位錯誤造成的誤分類風險。

真實證據

以下示範 germline phasing 與 haplotagging 指令：

# 1. 對 germline VCF 執行 phasing
longphase-s phase \
  -s SNP.vcf \
  -b alignment.bam \
  -r reference.fasta \
  -t 8 \
  -o phased_prefix \
  --ont

# 2. 依 phased VCF 將 read 標記為 HP tag
#    注意 -s 吃的是上一步的產物：phase 的 -o 給的是「前綴」，實際檔名是 <前綴>.vcf
longphase-s haplotag \
  -r reference.fasta \
  -s phased_prefix.vcf \
  -b alignment.bam \
  -t 8 \
  -o tagged_prefix

依目前版本，LongPhase-S 的 --ont／--pb 參數要求與 LongPhase-TO 的處理方式不同。 LongPhase-TO 僅支援 ONT，某些不相容參數可能被接受但不產生預期效果；請依版本文件核對， 不要預期流程在所有版本都提示警告。

完成後，可在 IGV 中依 HP tag 分組，觀察 read 是否呈現 haplotype 分層。

判讀練習

實作練習

執行 phasing 與 haplotagging

用上一節那兩道指令跑一次自己的資料（phase 產生 phased_prefix.vcf，haplotag 產生 tagged_prefix.bam）。 兩個步驟都完成後，再依下面三項檢查輸出。

以 IGV 檢查 haplotagging 結果

建議將 tagged_prefix.bam 載入 IGV，依序檢查以下項目：

1. 在 read 上按右鍵 → Group alignments by → phase（或 tag HP）
2. 前往一個 het 位點
3. 在相位品質良好的區域，通常可觀察到 read 分成兩層，且同一層的 read 在該位點上呈現一致 allele pattern

這項可視化有助於理解 haplotype 分層。若分群不清，請檢查輸入 VCF、coverage、mapping quality、 phase set 與工具參數，不能僅憑單一區域判定 phasing 失敗。

量化 phasing 的品質

# phase block 的 contig＋PS 鍵數；需搭配 block 長度與 switch error 評估
bcftools query -f '%CHROM\t[%PS]\n' phased_prefix.vcf | grep -v '\t\.$' | sort -u | wc -l

# 計算帶有 HP tag 的 read 比例
total=$(samtools view -c tagged_prefix.bam)
tagged=$(samtools view tagged_prefix.bam | grep -c 'HP:')
echo "標記率: $tagged / $total"

標記率沒有通用門檻可以照抄：它同時受 read 長度、heterozygous 位點密度、coverage、 平台與工具版本影響，同一份資料的不同區段也會不一樣。判讀方式是自己建立基線 —— 先看同一份資料在 het 位點密度正常的區段拿到多少，再拿可疑區段跟它比； 明顯偏低時才依上方「30% 沒有 HP tag」那份清單逐項排查。

為什麼重要

本模組比較兩種分析設計：有配對正常樣本的 tumor–normal（TN）模式， 以及只有腫瘤樣本的 tumor-only（TO）模式。是否有 matched normal 會改變可用的對照證據與不確定性。

差異不只是少一個輸入檔案；它會改變可辨識的變異類型與後續判讀方式。

概念與互動

同一批 read，兩種對照，兩種輸出

先把整件事放在一張圖裡。下圖的上下兩排腫瘤 read 完全相同—— 同樣六條 read、同樣四個候選位置。唯一改變的是右邊拿什麼當對照：

請注意這張圖真正的重點：位置 ① 與 ③ 都是這個人天生就有的 germline variantgermline variant 生殖系變異在生殖系形成、通常存在於多數細胞的變異；可遺傳給子代，但不代表必然傳遞。A variant arising in the germline and therefore present in essentially every cell of the individual. It may be inherited or arise de novo, and may be passed to offspring.完整條目 →， 在腫瘤樣本裡長得一模一樣（各三條支持、VAF 都是 0.5）。它們的命運不同， 不是因為變異本身有什麼差別，而是因為資料庫知不知道它。

位置 ④ 兩種模式都留得住，那是我們要的 somatic variantsomatic variant 體細胞變異在非生殖系細胞譜系中、通常於受精後取得的變異；可只存在於部分細胞，通常不由親代遺傳給子代。癌症基因體學常分析此類變異。A variant acquired outside the germline, usually after conception. It may be restricted to a subset of cells and is generally not passed to offspring; it is a central focus of cancer genomics.完整條目 →。 caller 產生的是 候選變異candidate 候選變異caller 提出、但尚未被確認的變異位點。本實驗室多數工具位於 caller 之後，進行再校正，處理的就是這些 candidates。A proposed but unconfirmed variant. Most tools in this lab post-process the caller's candidate set.完整條目 →，還不是最終答案；本教材以 ClairSClairS配對 tumor–normal 的 somatic variant caller。tumor-only 版本叫 ClairS-TO。A somatic variant caller for matched tumour–normal data; ClairS-TO is the tumour-only version.完整條目 →、DeepSomaticDeepSomaticGoogle 的 somatic variant caller，同樣有 tumor-only 版本。Google's somatic variant caller, also available in a tumour-only mode.完整條目 → 作為 tumour-normaltumour-normal 配對腫瘤／正常同時定序病人的腫瘤與正常組織；正常樣本提供個體 germline 背景，有助於辨識 somatic 候選。Sequencing a patient's tumour alongside their normal tissue, which supplies the germline background.完整條目 →（TN）calling 的範例，tumour-onlytumour-only 僅腫瘤只有腫瘤樣本、沒有配對正常樣本。在成本或檢體受限時很常見，但少了 germline 對照，分析難度會大幅提高。Having only the tumour sample. Common for cost or specimen reasons, but much harder without a germline baseline.完整條目 →（TO）則以 ClairS-TO 為例， 實際輸入與過濾流程請依版本文件確認。

沒有正常組織時，改用什麼當對照

在部分臨床或研究情境，matched normal 可能因成本、檢體量或取得條件而缺失。 替代方案是公開的族群資料庫：利用 allele frequency 與已知變異紀錄提供 germline 線索。 本教材在 LongPhase-TO／ClairS-TO 情境中將這組集合稱為 PONPON 正常樣本面板／族群資料庫（依流程而異）PON 依分析流程有兩種用法，兩者不可互換：本教材的 tumor-only 流程以它指稱 population germline database set（如 1000G、CoLoRSdb、dbSNP、gnomAD 的聯集）；GATK Mutect2 的 Panel of Normals（PoN）則由多個正常樣本建立，用來標記反覆出現的技術性 artifact。兩者都不能完全取代同一病人的 matched normal。The abbreviation is workflow-dependent and the two uses are not interchangeable. In this tutorial's tumour-only workflows, PON denotes a population germline database set; in GATK Mutect2, a Panel of Normals (PoN) is built from multiple normal samples to flag recurrent technical artifacts. Either can miss a private germline variant or an unseen artifact, so neither replaces a patient's matched normal.完整條目 →（population database set）， 它包含 1000 Genomes、gnomAD、dbSNP、CoLoRSdb 等來源。

上圖下排的橫軸就是這件事的關鍵。資料庫有一個收錄門檻—— 實際下載的檔名裡寫著 af-ge-0.001，意思是只收頻率 0.001 以上的位點。 位置 ① 的族群頻率 0.31，遠高於門檻，查得到；位置 ③ 只有 0.0002，落在門檻左邊， 資料庫裡根本沒有這一筆。「族群常見」與「這個人有」不是同一件事，缺口就從這裡來。

在 GATK Mutect2 語境，PoN 通常指由多個正常樣本建立的 recurrent-artifact panel。 名稱相同、內容不同 —— 而且更重要的是，這兩種資料各守住一個方向，中間會留下一格：

下方元件用同一組候選點數一次數量，看 precision 怎麼掉下來：

TN calling 也不是把兩個 VCF 相減

前面那張圖為了講清楚缺口，把 TN 流程畫得很乾淨：normal 有訊號就排除。 實務上不能這樣做，而且失準的地方在有配對正常樣本時就已經發生：

TINC：normal 裡的訊號強度說明三種不同的事

第一項值得單獨看一次，因為它最容易被寫成錯的規則。下圖三欄的腫瘤樣本完全相同 （24 條 read、6 條帶變異、VAF 都是 0.25），差別只在配對正常樣本裡有幾條：

VAF 可作為輔助訊息，但實際值取決於 purity、copy number、zygosity、multiplicity、 clonality 與 depth（完整關係式見 M8）。因此不能只用「有沒有訊號」或單一固定門檻判定。

沒有正常組織時，還剩下什麼證據

回到那張 PONPON 正常樣本面板／族群資料庫（依流程而異）PON 依分析流程有兩種用法，兩者不可互換：本教材的 tumor-only 流程以它指稱 population germline database set（如 1000G、CoLoRSdb、dbSNP、gnomAD 的聯集）；GATK Mutect2 的 Panel of Normals（PoN）則由多個正常樣本建立，用來標記反覆出現的技術性 artifact。兩者都不能完全取代同一病人的 matched normal。The abbreviation is workflow-dependent and the two uses are not interchangeable. In this tutorial's tumour-only workflows, PON denotes a population germline database set; in GATK Mutect2, a Panel of Normals (PoN) is built from multiple normal samples to flag recurrent technical artifacts. Either can miss a private germline variant or an unseen artifact, so neither replaces a patient's matched normal.完整條目 → 涵蓋不到的左下角。要把住在同一格裡的三種東西分開， 必須換一個問法：不再問「這個位置別人有沒有」，而是問「支持它的那些 read，彼此的結構像不像一次真的突變」。

真的變異一定是從某一條既有的 haplotypehaplotype 單倍型同一條實體染色體拷貝上，具有一致相位的一組 alleles 或 variants。HP1 與 HP2 是任意的相對標籤。A set of variants that lie on the same physical chromosome copy and are inherited together.完整條目 → 上長出來的，所以支持它的 read 應該集中在一條 haplotype 上。如果同一個「變異」的支持 read 分散在兩條染色體上， 那等於同一個位置各獨立突變了一次 —— 機率極低，更像比對或定序造成的問題：

一句話總結

真實證據

TN calling 指令示範

docker run \
  -v ${INPUT_BAM_DIR}:${INPUT_BAM_DIR} \
  -v ${OUTPUT_DIR}:${OUTPUT_DIR} \
  -u $(id -u):$(id -g) \
  hkubal/clairs:v0.4.1 \
  /opt/bin/run_clairs \
  --normal_bam_fn ${INPUT_BAM_DIR}/normal.bam \
  --tumor_bam_fn  ${INPUT_BAM_DIR}/tumor.bam \
  --ref_fn ${INPUT_REF_DIR}/reference.fasta \
  --threads 64 \
  --platform ont_r10_dorado_sup_5khz_ssrs \
  --output_dir ${OUTPUT_DIR}

TO 模式：準備 population database set

mkdir PoN
BASE=http://www.bio8.cs.hku.hk/clairs-to/databases
wget -P PoN $BASE/gnomad.r2.1.af-ge-0.001.sites.vcf.gz
wget -P PoN $BASE/dbsnp.b138.non-somatic.sites.vcf.gz
wget -P PoN $BASE/1000g-pon.sites.vcf.gz
wget -P PoN $BASE/CoLoRSdb.GRCh38.v1.1.0.deepvariant.glnexus.af-ge-0.001.vcf.gz

注意：有些文件把資料夾變數寫成 PATH="PoN" 再使用 $PATH；這會覆寫 shell 的可執行檔搜尋路徑。此處改用 BASE， 實際參數仍應依工具版本文件核對。

四個檔案名稱（gnomad、dbsnp、1000g、 CoLoRSdb）都指向族群 germline 資料庫，與本教材的 population database set 定義一致；它不同於 GATK Mutect2 所稱的 recurrent-artifact PoN。 其中兩個檔名帶著 af-ge-0.001 —— 那就是前面那張圖裡的收錄門檻， 低於這個頻率的位點不在檔案裡，所以查不到。

判讀練習

實作練習

比較兩種 calling 設計的結果

此練習比較 paired 與 tumor-only 的輸出差異。若有配對資料，可在相同參考序列與 相近過濾設定下執行兩種模式，再比較結果：

# isec 要求 bgzip 壓縮＋索引，少了索引會直接報錯
bcftools index -f paired.vcf.gz
bcftools index -f tumoronly.vcf.gz

# 配對模式找到、而僅腫瘤模式未找到的候選（數量依資料與設定而異）
bcftools isec -C paired.vcf.gz tumoronly.vcf.gz | wc -l

# 僅腫瘤模式多出的候選（可能包含 private germline、artifact 或 paired caller 漏報）
bcftools isec -C tumoronly.vcf.gz paired.vcf.gz | wc -l

若檔案還是未壓縮的 .vcf，先 bcftools view -Oz -o out.vcf.gz in.vcf 轉成 bgzip 格式再建索引；一般的 gzip 壓出來的檔案 bcftools index 不接受。

檢查多出來的候選是什麼

上面第二條指令的輸出，就是本模組那張圖裡的位置 ③。逐項檢查它們：

下面第 2 步先講一件事再看指令：HP 標籤的型別依工具而異。 longphase-s haplotag 寫的是整數 HP:i:， somatic_haplotag 與整個 longphase-to 寫的是字串 HP:Z:。 只 grep 其中一種，另一種資料會安靜地數出 0，所以下面用 HP:[iZ]: 同時涵蓋兩者。

# 1) 在配對正常樣本裡數這個位置的 ref/alt read —— 接近一半提示 private germline
#    第 5 欄是每條 read 在該位置讀到的字母，可以直接數（. 與 , 代表與參考相同）
samtools mpileup -f reference.fasta -r chr1:1000000-1000000 normal.bam

# 想要精確數字而不是用眼睛數，讓 bcftools 把 AD 欄挖出來
bcftools mpileup -f reference.fasta -r chr1:1000000-1000000 -a AD normal.bam \
  | bcftools query -f '%POS\t%REF\t%ALT\t[%AD]\n'

# 2) 支持 read 落在哪一條 haplotype 上
#    haplotag 過的 BAM 裡直接數 HP 標籤；兩邊都有支持就要提高疑慮
samtools view tumor.haplotagged.bam chr1:1000000-1000000 \
  | grep -oE 'HP:[iZ]:[12]' | sort | uniq -c

第 1 步的兩道指令是同一件事的兩種讀法：samtools mpileup 給的是可以用眼睛看的 pileup 字串，bcftools mpileup 給的是 VCF —— 後者預設不帶 allele depth， 所以要加 -a AD 才數得出來。

兩項要合起來看。VAF 接近 0.5 可提示 private germline， 但仍需排除污染、比對偏差與其他模型不確定性；haplotype 分布可提示 artifact， 但真的兩次獨立突變或 phasing error 都會造成例外。沒有單一證據可以定案。

原始文獻與程式碼

• ClairS（tumour-normal）：ClairS: a deep-learning method for long-read somatic small variant calling，bioRxiv，2023。程式碼在 github.com/HKU-BAL/ClairS。
• ClairS-TO（tumour-only）：ClairS-TO: a deep-learning method for long-read tumor-only somatic small variant calling，Nature Communications，2025。程式碼在 github.com/HKU-BAL/ClairS-TO。
• DeepSomatic：Accurate somatic small variant discovery for multiple sequencing technologies with DeepSomatic，Nature Biotechnology，2025。程式碼在 github.com/google/deepsomatic。

為什麼重要

到目前為止，「腫瘤樣本」在流程裡就是一袋 read。但那袋 read 來自一堆混在一起的細胞： 有癌細胞，也有正常細胞；癌細胞之間彼此也不完全相同；而且同一個位置在不同細胞裡可能有不同的 DNA 份數。

結果就是：同一個觀測到的數字，可以由好幾種不同的細胞組成算出來。 所以這一章的第一件事，是把 purity、copy number、multiplicity、CCF 這幾個名詞各自「在數什麼」講清楚 —— 它們常被混用，而混用會直接讓結論出錯。

概念與互動

先把名詞跟「在數什麼」對起來

這些名詞看起來很多，但它們其實只做一件事：在不同層級上數東西。 puritytumour purity 腫瘤純度樣本中腫瘤細胞所佔的比例。purity 越低，somatic 訊號被正常細胞稀釋得越嚴重，偵測越困難。The fraction of cells in a sample that are tumour cells. Low purity dilutes somatic signal.完整條目 → 與 CCFcancer cell fraction 癌細胞比例帶有某個特定突變的腫瘤細胞佔全部腫瘤細胞的比例。用來區分 clonal（$\approx 1$）與 subclonal（$<1$）突變。不等於 VAF。The fraction of tumour cells carrying a given mutation; distinguishes clonal from subclonal. Not the same as VAF.完整條目 → 數的是細胞； copy numbercopy number 拷貝數某段基因體在細胞內的拷貝數。多數正常常染色體區段為 2，可再分為 major 與 minor allele copy number。How many copies of a genomic segment a cell carries; a typical diploid autosomal segment has two. It can be split into major- and minor-allele copy numbers.完整條目 → 與 multiplicitymutation multiplicity在帶有某突變的細胞中，該突變所占的拷貝數；是 VAF 與 cancer cell fraction 換算時的重要參數。How many mutant-allele copies are present in a cell carrying the mutation; a necessary term when converting VAF to cancer cell fraction.完整條目 → 數的是一個細胞裡有幾份 DNA； VAFVAF 變異等位基因頻率在某個位點上，支持 alt allele 的 read 佔全部 read 的比例。VAF 不等於帶有這個突變的細胞比例。Variant Allele Frequency: the fraction of reads supporting the alternate allele. Not the same as the fraction of cells carrying it.完整條目 → 數的是read。搞混它們，通常就是搞混了分母。

下面用一個玩具檢體把三層對齊：10 個細胞，其中 6 個是腫瘤細胞， 腫瘤細胞裡有 3 個帶著我們要看的那個突變，每個細胞在這個位置都有 2 份 DNA。

這張圖裡有兩件事值得停下來看。

第一，方向是不對稱的。給定上面兩層，VAF 只有一個答案 —— 腫瘤佔 0.6、其中一半帶突變、每個細胞 2 份裡帶 1 份，算出來就是 $0.6\times 0.5÷2=0.15$。 但拿著 $\mathrm{VAF}=0.15$ 往上反推，會有很多組合都成立。整章的難處都在這個不對稱上。

第二，read 不是細胞。一條 read 是一個 DNA 分子，同一個細胞會貢獻很多條 read， 而且哪一條 read 來自哪個細胞在 bulk 資料裡查不到。所以「數 read」永遠只是「數細胞」的間接代理。

copy number、multiplicity 與 ploidy：一個位置有幾份

中間那層最容易被跳過，但它是後面所有推論的地基。三個詞的差別很小、影響很大： copy number 是「這一段在一個細胞裡有幾份」，multiplicity 是「其中幾份帶著這個突變」， ploidyploidy 倍性細胞內染色體套數。多數正常常染色體為 2（diploid）；腫瘤可呈現 3、4 或其他非整倍體狀態。The number of chromosome sets in a cell. Most normal autosomal regions are diploid (2); tumours may be triploid, tetraploid or otherwise non-diploid.完整條目 → 是「把整個基因體的 copy number 平均起來」。

為什麼要分這麼細？因為同一個突變、同樣「每個腫瘤細胞都有」， 只要拷貝狀態不同，VAF 就完全不同：

把「腫瘤細胞佔多少」「一個細胞有幾份」「突變佔幾份」「多少腫瘤細胞帶著它」四件事寫進同一個式子， 就是下一節的關係式。

VAF 的完整關係式

M1 使用了簡化情境：「purity=0.5、樣本為 diploid，且每個腫瘤細胞帶有一份變異拷貝，因此 VAF 為 0.25」。在一般情況下，關係式為：

VAF 是可由 reads 估計的觀測量；purity、totalCN、multiplicity 與 CCF 通常需由多項證據估計或推論。 這些參數之間存在可識別性限制，因此「估計腫瘤細胞比例」本身就是一個獨立的研究題目。

LOH 的常見機制及其可觀測差異

LOHLOH 異型合子性喪失原本 heterozygous 的區域變成只剩一種 allele。LOH 不等於缺失 —— 也可能是一條 haplotype 遺失後另一條被複製（copy-neutral LOH）。Loss of heterozygosity: a formerly het region retains only one allele. Not necessarily a deletion.完整條目 → 是「原本 heterozygous 的區域只剩下一種 allele」。 但其形成機制不只一種；以下比較兩類常見情形：

低 purity 對 somatic variant detection 的影響

回到最上面那一層。purity 低的時候，突變本身沒有變少，但它在整池 DNA 裡被稀釋， 落到 read 上就只剩幾條支持：

以下把最左邊那一格算清楚：50× coverage、purity 0.2，以及一個 clonal 的 diploid somatic 變異：

因此預期僅有約五條支持 read。ONT 的錯誤率受 basecaller、序列脈絡與 read 品質影響； 在此深度下，少量錯誤 read 可能與變異訊號重疊，降低可分辨性。

這說明 LongPhase-S 為何採用 purity 感知的過濾：低 purity 通常需要重新校準過濾參數， 並以驗證資料在 sensitivity 與 precision 間取得平衡。

在 LongPhase-S 的此實作中，流程先估計 purity，再依 purity 調整過濾參數；其他方法也可能採用聯合或迭代估計。

真實證據

LongPhase-S 以 ASCAT 作為 purity 估計的比較基準。ASCAT 是常用的 allele-specific copy-number、purity 與 ploidy 估計方法； 具體適用範圍與資料平台應依研究設計說明。

該研究報告 LongPhase-S 在低 purity 合成樣本上的估計誤差多數低於 10%， 交叉驗證 MAE 約 4%、$R^2$ 約 0.95；解讀時應同時確認誤差定義、資料集與切分方式。

這項比較僅涵蓋合成樣本，且假設來源 tumor BAM 接近純腫瘤。 細胞株仍可能包含少量正常或異質細胞，因此此假設限定了 purity 結果的適用範圍。

先提出預測，再查看解答

實作練習

合成不同 DNA fraction 的混合樣本

以下是評估混合樣本穩健性的常用設計之一：

# 1. 先量測兩個 BAM 的有效 coverage
mosdepth --by 500 tumor_cov  tumor.bam
mosdepth --by 500 normal_cov normal.bam

# 2. 依實測 coverage 回推抽樣率（-s 的格式是 seed.fraction）
#    抽樣率 =（目標 fraction × 合成後總深度）÷ 來源深度 —— 推導與五個梯度的數值見 M3
#    以 tumor 50×、normal 25×、總深度 25×、目標 DNA fraction 0.6 為例：
samtools view -s 1.30 -b tumor.bam  > t_sub.bam
samtools view -s 1.40 -b normal.bam > n_sub.bam

# 3. 合併
samtools merge dna_frac_0.6.bam t_sub.bam n_sub.bam
samtools index dna_frac_0.6.bam

可建立 0.2 / 0.4 / 0.6 / 0.8 / 1.0 五個 DNA fraction 梯度， 以評估方法在低 DNA fraction 條件下的效能變化；這是一組常見的基礎評估設計。

以 heterozygous 位點密度初步篩查 LOH

# 計算沿染色體分箱的 het 位點密度
# 需與配對 normal、鄰近可呼叫區域及 coverage 比較，以列出 LOH 候選
bcftools view -H -g het sample.vcf | awk '{print $1, int($2/1000000)}' | \
  sort | uniq -c | awk '{print $2, $3, $1}' | head -40

可呼叫區域的 het 密度受樣本、族群、coverage、caller 與過濾條件影響，應以配對 normal 或鄰近區段建立基線。 若某區段密度顯著低於基線，可列為 LOH 候選；purity 混合可能保留部分 het reads，仍需結合 allele balance、coverage 與其他證據。

為什麼重要

一條 ONT read 除了序列，還帶著一串「這個 C 有沒有被加上甲基」的紀錄。 這件事的用處不只是「多一種資料」：同一群細胞的甲基化模式會比較像， 所以甲基化其實是一個關於「這條 read 來自哪一群細胞」的線索 —— 而 bulk 資料本來查不到這件事。

實驗室有兩個題目用到這個線索，但用法剛好相反： 一個把它當成判斷依據（會影響結果），另一個刻意不讓它影響判斷（只當事後註記）。 這一章先講清楚兩個題目各自在做什麼，再談為什麼其中一個要刻意不用。

概念與互動

5mC：序列之外的表觀遺傳資訊

5mC5mC 5-甲基胞嘧啶胞嘧啶第 5 個碳上的甲基化修飾。這是常見的 DNA 甲基化形式，主要見於 CpG；可影響轉錄調控，效應依基因組位置與細胞類型而異。5-methylcytosine, the commonest DNA methylation mark, usually at CpG sites. Affects expression without changing sequence.完整條目 → 是在胞嘧啶（C）的 5 位碳上加上甲基。在序列字母表示上，它不改變 C 的識別； 5mC 可與基因調控及表現相關，但作用取決於基因組脈絡。

在哺乳類，5mC 主要集中於 CpGCpG序列上一個 C 後接一個 G（p 代表兩者之間的磷酸鍵）。哺乳類多數 5mC 位於 CpG，特定細胞或情況亦可見非 CpG 甲基化。A cytosine followed by a guanine. Most mammalian 5mC occurs at CpG sites, although non-CpG methylation is also observed in particular cells and contexts.完整條目 → 位點（一個 C 後面接著一個 G），但也可能出現在非 CpG 脈絡。 在本節的分析中，基本單位是「某個 CpG 位點上，一條 read 的 methylated 或 unmethylated 狀態」。

bisulfite-seq 是測量甲基化的常用方法之一，會將未甲基化的 C 轉換為 U（在定序中通常讀為 T），流程較為繁複。 ONTONT 奈米孔定序Oxford Nanopore Technologies。以電流訊號讀取 DNA，可產生長 read；在 native-DNA 流程與適當模型下，可由訊號推定 5mC，通常不需 bisulfite 轉換。Oxford Nanopore Technologies. Reads DNA through ionic-current signals, producing long reads. With native DNA and a suitable model, 5mC can be inferred without bisulfite conversion.完整條目 → 可利用電流訊號搭配修飾辨識模型推定 5mC，通常不需 bisulfite 前處理； 在完成修飾辨識且保留 MM/ML 後，同一條 read 可同時提供序列變異與甲基化狀態。

這些資料在程式裡是什麼形狀

先把資料結構看清楚，後面兩種用法就都好懂了。 一條 read 在一個區域內會經過好幾個 CpG 位點，每個位點上它是 methylated 或 unmethylated —— 所以一條 read 就是一個 0／1 向量，一個區域就是一張 read × CpG 的表。

用法一：當作證據 —— 判斷一個候選是不是 somatic

第一個題目是 tumor-only 的 somatic 變異判讀：手上有一堆候選位置，要決定哪些是真的 somatic 變異。 甲基化在這裡派上用場的理由很具體 —— somatic 變異只存在腫瘤細胞裡，所以帶著它的那些 read 全部來自同一群細胞， 而不帶它的 read 則來自腫瘤與正常細胞的混合。兩側的甲基化組成因此會不一樣。

實際的特徵有三類：兩側各自的甲基化覆蓋量與其差距、每一側是單一還是混合的比例， 以及把兩側的平均甲基化程度畫成散布圖後的形狀差異。 分類器（XGBoost）本身不懂生物學，它只是學會「這幾個數字長這樣時，候選通常是真的」。

用法二要用到的任務：subclone 重建在做什麼

第二個題目常被術語擋住，但它的核心其實是一個很單純的組合問題。

腫瘤細胞一代一代分裂，過程中只會多拿到突變，不會把突變還回去。 所以拿一個小區域、看其中 k 個 somatic 位置，每個細胞在這 k 個位置上就是一個 0／1 向量； 一個細胞分裂出的後代，向量只會多出幾個 1。

bulk 定序給你的每一條 read，就是這些向量的一次取樣（read 夠長才能同時覆蓋好幾個位置）。 於是問題變成：已知看到了哪些向量，什麼樣的「一次翻一個 0→1」的路徑能最省步數地產生它們？ 那條路徑就是這個區域的候選演化順序。

這樣看的話，subclone 重建的輸出就是「這個區域內、突變發生順序的候選解」， 而不是一棵完整的腫瘤演化樹。也因為是最小成本解，並列的候選會全部保留，不強行選一個。

同一個狀態裡，read 還能再分群

現在把甲基化接回第二個題目。假設一群 read 的 haplotype 相同、在每個 somatic 位置上的 0／1 向量也完全相同 —— 依序列來看它們是「同一個狀態」，無法再分。 但它們的甲基化模式可能明顯分成兩群：

可觀察到：ALT reads 仍可再分成兩群。 它們的 allele 與 haplotype 相同，僅依序列無法區分，但甲基化模式可能不同。

這個現象稱為 within-ALT residual structure；研究依其統計定義，在七個資料集中將其量化為主要殘餘變異來源（51.2%–80.4%）。 換句話說：固定 haplotype 與 allele 後，仍可觀察到額外的 read-level 甲基化分群。

它的意義是：上一節那棵樹的一個節點裡，可能還藏著更細的分群 —— 序列證據已經用完了，但資料還沒說完。這正是甲基化在第二個題目裡的位置： 它標記已經建好的節點，而不參與建樹。

兩種角色的對照

兩個題目用的是同一份 read × CpG 表，但甲基化在流程裡的位置完全不同：

為什麼用法二要刻意不用？因為那個題目最想知道的，正是「這棵樹對不對」。 如果甲基化已經幫忙選過樹，再拿甲基化差異來說「看，這棵樹的分支確實有甲基化差異」， 那只是把自己的輸入念一次而已 —— 這是循環論證。

同一個原則也適用於用法一：那裡的甲基化已經進了模型，所以評估模型效能時 必須用沒有參與訓練的資料，不能拿甲基化特徵本身當成「模型是對的」的證據。

真實證據

實務限制之一是甲基化資料的覆蓋範圍。這裡要留意計數單位：M2 列出的是六個細胞株， 而評估用的是資料集，同一株可能有多份（HCC1395 的 HKU 與 NYGC）。 以細胞株計，六個裡只有五個有可用的甲基化資料；以資料集計則多於五份。 不論用哪個單位，能做甲基化分析的樣本都比能做序列分析的少； 因此甲基化較適合作為區域層級的補充證據，而非唯一判準。

研究統計各資料集中帶有穩定甲基化差異的 sSNV 位點比例， COLO829 為 9.5%、H2009 為 35.4%，資料集間差異明顯。不可評估的主要限制之一是 ALT read 的 coverage 不足； 要估計兩組甲基化差異，兩邊都需要足夠且相對平衡的 reads。

先提出預測，再查看解答

實作練習

甲基化資訊存在 BAM 的哪裡

前面那張 0／1 表不是憑空來的，它從 BAM 的兩個 tag 讀出來：

samtools fastq -T '*' methylcall.raw.bam > raw.fastq
minimap2 -ax map-ont -y reference.fasta raw.fastq

確認 BAM 是否包含甲基化資料

先確認必要的 MM/ML tags，以免後續修飾分析缺少資料：

# 抽查一條 read 是否包含 MM / ML tags
samtools view methyl.bam | head -1 | tr '\t' '\n' | grep -E '^(MM|ML|Mm|Ml):'

# 統計帶有 MM tag 的 reads（抽查多條 reads，不只查看第一條）
total=$(samtools view -c methyl.bam)
withmod=$(samtools view methyl.bam | grep -c 'MM:Z:')
echo "$withmod / $total 條 read 帶甲基化 tag"

若第二個數字為 0，先確認原始 BAM 是否已完成修飾辨識及 tag 命名格式；可能是尚未呼叫修飾，也可能在前處理中遺失。可用以下流程重新保留 tags：

samtools fastq -T '*' original.bam > with_tags.fastq
minimap2 -ax map-ont -y reference.fasta with_tags.fastq | \
  samtools sort -o remapped.bam && samtools index remapped.bam

提取並統計甲基化訊號

# 產生每個 CpG 位點的甲基化比例（bedMethyl 格式）
modkit pileup methyl.bam output.bed --ref reference.fasta

# 依 haplotype 分開統計（需要先 haplotag）
modkit pileup methyl.bam out_hp.bed --ref reference.fasta --partition-tag HP

--partition-tag HP 可將甲基化依 haplotype 分層統計， 以比較兩條 haplotype 的甲基化模式。依 M6 所述，haplotype 方向僅在同一 phase block 內一致， 因此統計應限制於同一 phase block，或先完成方向校正。

在 IGV 中檢視甲基化

在支援 base-modification 且 BAM 含有效 MM/ML、參考序列與相容版本時，IGV 可顯示 ONT 甲基化。載入 BAM 之後右鍵 → Color alignments by → base modification (5mC)。 搭配前面的 Group by phase，即可同時檢視 haplotype 與甲基化資訊。

為什麼重要

後續閱讀的論文通常會包含多張長條圖與折線圖。 判讀圖表未呈現的資訊，是研究與工程實務中的重要能力。

本節的目標不只是計算 F1，也包括培養批判性解讀與驗證能力。

本模組的數字來自兩項工作

下面會反覆引用兩組具體結果。它們是不同的方法、處理不同的問題，先在這裡分清楚， 後面提到時就直接用名字稱呼：

第一項的那句「只會刪掉候選」，是本模組第一個要拆解的東西 —— 它替所有後處理方法設下了一個效能上限。

概念與互動

先釐清三項指標

第三欄的條件限定上述結論。在固定 caller 候選集合且後處理僅移除候選時，recall 至多維持，若誤刪 TP 則可能下降。 caller 未產生的候選不會由此恢復；若後處理包含 re-calling 或 assembly，則需另行評估。

觀察 precision 與 F1 的差異

truth settruth set 標準答案集評估用的參考標註。可靠度依資料來源與建立方法而異，困難區域可能包含標註錯誤。The reference labels used for evaluation. A truth set is not infallible truth: reliability varies, and hard regions can contain label errors.完整條目 → 並非無誤的標準

評估需要參考答案；不同資料來源的標註可靠度可能存在差異：

LongPhase-S 在搜尋過濾參數時， 依不同 cell line 的 benchmark 可靠度乘上相對應的權重。 此設計可能減少 benchmark 差異的影響，但比較不同論文的數字時，仍應確認 truth set 與加權方式是否一致。

三項可能造成評估偏差的因素

最後一項可由具體案例說明。下面並列 indel 精修模型送進 CNN 的兩張影像 —— 也就是兩個 deletion 候選各自的 pileup 編碼：一個是 true positive（模型給 0.997）， 一個是 false positive（0.994）。 兩者的 tumor reads 都顯示集中一致的 deletion，normal reads 未見相應支持。 就這種影像編碼呈現的特徵而言，兩個候選近似。

該工作的結論是：在這種影像輸入下，這兩類候選沒有可辨識的差異。這類位點常落在 homopolymer 或 tandem repeat 區域， benchmark 標註的不確定性可能較高；若要判定 label noise，仍需獨立 truth 或重複實驗支持。

合成 DNA fraction 混合樣本及其適用範圍

常用流程假設來源 tumor BAM 接近純腫瘤，再依指定比例與 normal BAM 混合。 此流程可產生預先設定的 DNA／read fraction，作為近似真值；但來源樣本的純度、輸入深度、每細胞 DNA 量與拷貝數仍需校正， 實際細胞株也可能包含少量正常或異質細胞。

因此「在標示為 purity 0.2（或 DNA fraction 0.2，依研究定義）的樣本上誤差小於 10%」應完整表述為「在以此方式合成的樣本上」； 結果僅適用於該合成方式與來源樣本前提。

真實證據

用對的指標選模型

indel 精修模型有一項值得注意的模型選擇設計：它不取 validation loss 最低的那個 epoch 當最佳模型。

理由是 loss 反映所選損失函數（例如 cross-entropy）的平均目標，不必然等同於分類錯誤率； 若任務目標是「在盡量保留真 indel 的前提下，篩掉更多假的」， 兩者不是同一件事。

因此它改用以下任務導向分數挑 epoch：

λ 由 benchmark 與驗證集的正負比例估計，而非人工指定（該資料約為 1.45）； 應同時說明估計資料、FP 總數為零時的處理，以及權重敏感度。

本例交叉驗證以樣本為分割單位

LongPhase-S 用 leave-k-out 交叉驗證評估腫瘤 DNA 比例的迴歸模型， 每次留下 1 或 2 個 cell line 當測試集。結果 MAE 約 4%、$R^2$ 約 0.95。

這種做法叫 cross-validationcross-validation 交叉驗證輪流保留資料子集作為驗證折，用來估計模型的泛化能力與選擇設定；這些折是 validation，不是最終 test。最終測試集應在模型與設定確定後才使用，並另行保留。Repeatedly holding out validation folds to estimate generalisation and choose model settings. The folds are validation data, not the final test; an untouched final test set is used only after the model is fixed.完整條目 →。此處保留 cell line 作測試，有助於評估跨樣本的泛化能力。 若同一 cell line 的高度相關位點同時出現在訓練與測試中，可能造成 data leakagedata leakage 資料洩漏訓練或模型選擇階段取得了評估資料的資訊，讓效能被高估。實驗室的做法是按染色體切分，確保同一個位點不會同時出現在訓練、驗證與最終測試中。Evaluation information entering training or model selection and inflating performance. The lab splits by chromosome so a locus cannot appear in training, validation and final test partitions at once.完整條目 → 或效能高估。

先提出預測，再查看解答

實作練習

以標準化流程評估效能

除非已處理 variant normalization，否則不建議自行實作比對； 同一個 indel 可能有多種等價表示，多等位與複雜變異也需要額外處理。可使用專門工具：

# SNV / indel 的常用評估工具（hap.py）
docker run -v $(pwd):/data pkrusche/hap.py /opt/hap.py/bin/hap.py \
  /data/truth.vcf.gz \
  /data/query.vcf.gz \
  -f /data/high_confidence_region.bed \
  -r /data/reference.fasta \
  -o /data/result

-f 是關鍵參數之一：它限制在 benchmark 信心較高的區間內評估。 未限制區域時，應分層報告並說明標註信心，否則結果難以直接解讀。

解讀效能數字時的四項檢核

解讀效能數字時，可依序確認：

1. 另一個指標呢？若只報 precision，應同時檢視 recall 的變化。
2. 在哪個資料集、用哪一份標準答案、有沒有限制在高信心區間內？
3. 測試資料是否參與過調參？若門檻在同一批測試資料上挑選，結果可能是對該資料集過度擬合的估計。
4. 絕對數量是多少？precision 從 0.72 到 0.89，如果 TP 只有 20 個，統計上不穩定。

檢查訓練／測試資料是否存在洩漏

# 檢查訓練與測試的染色體名稱是否重疊（僅為初步檢查）
cut -f1 train_sites.tsv | sort -u > /tmp/train_chr
cut -f1 test_sites.tsv  | sort -u > /tmp/test_chr
comm -12 /tmp/train_chr /tmp/test_chr    # 有輸出表示染色體名稱重疊，仍需檢查座標與樣本

# 比較訓練與測試的類別比例
awk '{c[$NF]++} END {for (k in c) print k, c[k]}' train_sites.tsv
awk '{c[$NF]++} END {for (k in c) print k, c[k]}' test_sites.tsv

類別比例差異會影響評估：若訓練集的假陽性佔 4%、測試集佔 19%， 兩組資料的類別分布不同，不宜僅以 loss 判定模型優劣；仍應檢查其他特徵與品質分布。

為什麼重要

LongPhaseLongPhase實驗室開發的長 read phasing 工具，是後續所有工具的共同基礎。輸出 germline haplotype。The lab's long-read phasing tool and the common foundation for everything else. Produces germline haplotypes.完整條目 → 是本教材後續工具共用的基礎，處理的是一般（非腫瘤）樣本的 phasing。 開始分析腫瘤資料前，先掌握這個基本流程。

它值得看清楚內部機制，原因是：後面 LongPhase-S 與 LongPhase-TO 的腫瘤專用處理， 都是接在這張圖與這套投票之上。不知道圖長什麼樣，就看不出腫瘤版改了什麼。

概念與互動

輸入與輸出

注意這個順序：先 call variant，再 phase，最後 haplotag。 LongPhase 本身不做 variant calling —— 它需要別人先給它 het 位點當錨點。

四種變異，一張圖

LongPhase 最關鍵的設計選擇，是把四種變異放進同一張圖一起定相， 官方文件稱為 co-phasing。理由來自長讀本身：一條 read 可能同時跨過一個 SNP、 一個甲基化位點、一個小的 indel 與一個 SV —— 這四件事被同一個分子綁在一起， 分開 phasing 就把這個證據丟掉了。

程式內部給每個節點記一個型別，原始碼裡的註解就寫著 0=SNP 1=SV 2=MOD 3=INDEL 4=tandem repeat INDEL。 型別不只是標籤，它會改變後面的判斷 —— 例如串聯重複區的 indel（型別 4）投票權重被降到 0.1， 因為那類位置的定序錯誤率特別高。

一條邊怎麼決定方向，又怎麼被丟掉

邊不是「有」或「沒有」，而是四個加權計數。跨過兩個位置的每一條 read， 依它在兩邊看到的 allele 落進 rr、ra、ar、aa 之一：

兩個細節值得記住，因為它們解釋了很多「為什麼這裡沒接起來」：

• 低品質的鹼基會被打折。邊的任一端 base quality 低於 --baseQuality（預設 12）時， 那條 read 對這條邊只貢獻 --edgeWeight（預設 0.1），而不是 1。
• 甲基化的門檻更嚴。SNP 與甲基化位點之間的邊，門檻從 0.7 收緊到 0.3 —— 甲基化訊號比較吵，要求證據更乾淨才採信。

heuristic phasing：加權投票

圖建好之後，LongPhase 不求全域最佳解。它沿著座標從左往右掃， 一次決定一個變異：前面最多 --connectAdjacent（預設 35）個變異各投一票， 票的權重就是上面那個 20 / 1 / 0.1，兩邊加總大的一邊贏。

這是一個典型的 heuristic：快、在長讀上通常夠準，但代價是它對兩件事敏感 —— 證據剛好平手，以及單一 read 主導。③ 那個特例保護就是為了後者而存在的。 另外，距離超過 --distance（預設 300000，約一個 centromere 的長度）的兩個變異不互相投票。

最後兩道信心門檻

圖定好相之後還有一步，程式把結果反過來檢查一次：先用變異的相位決定每條 read 屬於哪條 haplotype，再用已標好的 read 回頭重算每個變異。兩步各有一道門檻：

三種失敗模式

第一項與第三項在腫瘤樣本上可能更明顯，因為 LOHLOH 異型合子性喪失原本 heterozygous 的區域變成只剩一種 allele。LOH 不等於缺失 —— 也可能是一條 haplotype 遺失後另一條被複製（copy-neutral LOH）。Loss of heterozygosity: a formerly het region retains only one allele. Not necessarily a deletion.完整條目 → 與 copy number 異常會減少可用的 phasing 錨點。 LongPhase-S 與 LongPhase-TO 都針對這些腫瘤資料情境提供額外處理。

真實證據

最小可跑範例

# 1. 準備輸入（見 M4）
samtools faidx reference.fasta
minimap2 --MD -ax map-ont -t 10 reference.fasta reads.fastq -o aln.sam
samtools sort -@ 10 aln.sam -o aln.bam && samtools index aln.bam

# 2. germline calling（Clair3 或 PEPPER-DeepVariant，見 M5）
#    產生 SNP.vcf

# 3. phasing
longphase-s phase \
  -s SNP.vcf \
  -b aln.bam \
  -r reference.fasta \
  -t 8 \
  -o phased \
  --ont

# 4. haplotagging
longphase-s haplotag \
  -r reference.fasta \
  -s phased.vcf \
  -b aln.bam \
  -t 8 \
  -o tagged

co-phasing：把其他三種變異一起帶進去

# 甲基化要先用 modcall 產生 VCF，才能餵給 phase
longphase-s modcall \
  -b aln.bam \
  -r reference.fasta \
  -t 8 \
  -o modcall

# 四種一起 co-phase；每種變異會各自輸出一份 phased VCF
longphase-s phase \
  -s SNP.vcf \
  -b aln.bam \
  -r reference.fasta \
  --indels \
  --sv-file SV.vcf \
  --mod-file modcall.vcf \
  -t 8 \
  -o phased \
  --ont

注意：modcall 的 help 顯示 -b, --bam-file， 但選項表註冊的長選項其實是 --methylbamfile；短選項 -b 兩者都通。 同理 phase 的 help 寫 --sv-window 與 --sv-threshold， 選項表註冊的是 --svWindow 與 --svThreshold。

把圖直接印出來看

前面那三張圖不必只當示意圖 —— --dot 會為每一條染色體輸出一個 dot 檔，裡面就是實際建出來的邊：

longphase-s phase -s SNP.vcf -b aln.bam -r reference.fasta --ont --dot -o phased

# 每一行是一條邊：<位置>.<allele> -> <位置>.<allele>
# allele 相同 = 平行連接，相反 = 交叉連接
head -20 phased.chr20.dot

如何驗證輸出

# 在已定相的 heterozygous records 檢查 | 與 PS
grep -v '^##' phased.vcf | head -5

# BAM 裡應該有 HP tag
samtools view tagged.bam | grep -o 'HP:i:[12]' | sort | uniq -c

# 有多少 read 沒被標到？
samtools view -c tagged.bam
samtools view tagged.bam | grep -c 'HP:i:'

最後兩行的差就是沒過 --readConfidence（0.65）那道門檻的 read。 請將差值除以總 read 數，再依資料品質與分析目標判讀；比例偏高時，回到 M6 的排查清單。

參數對應到哪一步

參數很多，但每一個都落在前面某一張圖上。按階段記比按字母記容易得多：

判讀練習

實作練習

執行最小流程

# 僅執行小區域測試，以縮短重複測試時間
REGION="chr20:1000000-2000000"

longphase-s phase \
  -s SNP.vcf \
  -b aln.bam \
  -r reference.fasta \
  -t 8 \
  -o phased \
  --ont

longphase-s haplotag \
  -r reference.fasta \
  -s phased.vcf \
  -b aln.bam \
  --region "${REGION}" \
  -t 8 \
  -o tagged

把兩道門檻的效果量出來

這個練習把前面那張「兩道信心門檻」的圖變成實際數字。 調鬆 --readConfidence，看有多少 read 從沒有 HP 變成有：

# 預設 0.65
longphase-s phase -s SNP.vcf -b aln.bam -r reference.fasta --ont -o strict
# 調鬆到 0.55
longphase-s phase -s SNP.vcf -b aln.bam -r reference.fasta --ont -m 0.55 -o loose

# 各自 haplotag 之後比較 HP 標記率
for p in strict loose; do
  total=$(samtools view -c ${p}.tagged.bam)
  tagged=$(samtools view ${p}.tagged.bam | grep -c 'HP:i:')
  echo "${p}: ${tagged} / ${total}"
done

標記率會上升 —— 但那不代表結果更好。放寬門檻是把不確定的 read 也貼上標籤， 標記率與正確率是兩件事，要用有 truth 的資料才分得出來（見 M10）。

確認參數名稱

# 看目前版本接受的長選項
longphase-s phase --help

# help 與實際行為不符時，以選項表為準（檔案在 repo 根目錄）
grep -A40 'static const struct option longopts' Phasing.cpp

原始文獻與程式碼

LongPhase：Lin JH、Chen LC、Yu SQ、Huang YT， LongPhase: an ultra-fast chromosome-scale phasing algorithm for small and large variants，Bioinformatics，2022。程式碼在 github.com/twolinin/LongPhase。

為什麼重要

前面幾章的 phasing 都假設「一個人只有兩條染色體拷貝」。腫瘤樣本不成立： 裡面混著正常細胞與好幾群癌細胞，某些區段還多了或少了拷貝。 在這個前提下，把 read 分成 HP1／HP2 兩堆已經不夠用。

LongPhase-SLongPhase-S配對 tumor–normal 版本。做腫瘤 DNA 比例估計（輸出欄位名為 purity）、比例感知的 somatic variant 過濾，以及 somatic haplotagging。The matched tumour–normal version: tumour DNA fraction estimation (the output field is named purity), fraction-aware somatic filtering, and somatic haplotagging.完整條目 → 的做法是多做一層：除了兩條 germline haplotype， 再把「帶著 somatic 變異的那些 read」單獨拉成它們的後代分支。 這一層做出來之後有三個連帶好處 —— 可以直接估腫瘤來源 DNA 佔了多少、 可以用來重新判斷變異真假、也可以看出哪些 read 屬於同一個癌細胞群。

第一項要先把量講清楚，因為工具的參數名稱與輸出欄位都叫 purity， 但它估的是 腫瘤 DNA 比例tumour DNA fraction 腫瘤 DNA 比例樣本 DNA 中源自腫瘤的比例。與 tumor purity（細胞比例）在 aneuploid 或 WGD 的情況下會不一樣。The fraction of DNA originating from tumour cells; diverges from cellular purity under aneuploidy or WGD.完整條目 →（分母是分子）， 不是 cellular puritytumour purity 腫瘤純度樣本中腫瘤細胞所佔的比例。purity 越低，somatic 訊號被正常細胞稀釋得越嚴重，偵測越困難。The fraction of cells in a sample that are tumour cells. Low purity dilutes somatic signal.完整條目 →（分母是細胞）。 M3 與 M8 的重要區分 #2 已建立這兩個量在非二倍體之下並不相等； 本章說明 LongPhase-S 落在哪一邊，以及為什麼它不可能落在另一邊。

概念與互動

先看懂 somatic haplotagging 本身

整章的其他東西都建立在這一步上：somatic haplotaggingsomatic haplotagging把帶有 somatic 突變的 read 指派到它們源自的 germline haplotype。這是 LongPhase-S 的核心產出。Assigning somatic-mutation-bearing reads to the germline haplotype they arose from. The core output of LongPhase-S.完整條目 → 就是給每條 tumor read 貼一個標籤。 判斷只用兩件事 —— 這條 read 上的 germline allele 屬於哪一條，以及它有沒有帶 somatic allele。