中篇：VAF 路線缺什麼，多變異區域能補上哪三格

為什麼重要

本頁的目標是改善全基因體尺度的 subclone 系統發生重建， 而非增加局部陳述的數量。

上篇的結論界定了起點。在固定寬度 20 kb 視窗、每 Mb 5 個變異的條件下， 含變異的固定寬度視窗中有 95% 僅含單一變異；此類視窗內「一條單倍型」與「一個變異」一一對應， 因此捨去標籤後，單倍型頻率譜近乎退化為 VAF 頻率譜。 頻率值本身並非這條路線的主要貢獻。

其餘 5% 的連鎖視窗所提供的資訊性質不同：它們給出的不是更多資料點， 而是資料點之間的關係 —— 哪兩個變異不可能屬於同一群、哪一群是哪一群的祖先。 而這類關係不需要大量：一棵 clone treeclone tree 克隆演化樹描述腫瘤內各群細胞祖先關係的樹：節點是一群帶有相同變異組合的細胞，邊代表在祖先之上又多拿到變異。要注意同一組群集常常有多棵樹同時相容。A tree describing ancestral relationships among cell populations in a tumour: nodes are groups of cells sharing a mutation set, edges represent additional mutations acquired on top of an ancestor. Multiple trees are often compatible with the same clusters.完整條目 → 所要確定的是群與群之間的少數幾個關係， 而非一萬五千個變異各自的歸屬。

「進入同一個混合模型」一語尚未指明任何機制。 本頁要確立的是兩件事：該混合模型的具體形式，以及 5% 的觀測作用於其中哪一個變數。 第一節先建立僅以 VAF 重建的完整流程作為對照基準， 其後三節逐一說明 5% 能補上的三格。 把三者合成一條估計式、以及那條估計式能不能真的建起來，是下篇的內容。

此外，本實驗室現行的實作已經在計算上圖右欄的內容。 它對每一個單倍型連鎖區段（一個連鎖視窗 × 一條 單倍型家族單倍型家族 一條 germline 單倍型，加上由它衍生的 somatic 單倍型把 read 依 HP tag 分成的兩組之一。家族一是 HP1 與從它長出來的 HP1-1，家族二是 HP2 與 HP2-1；歸不到任何一條 germline 單倍型的 HP3 不屬於任何一族。叫「家族」是因為它把一條 germline 單倍型與由它衍生的 somatic 單倍型收在同一組裡 —— 分組看的是 germline 那一層，不是有沒有帶 somatic 突變。

在同一個 phase block 內，一個家族對應一條染色體拷貝，所以「兩族」就是那個位置上的兩條同源染色體。但兩件事不成立：其一，軟體判定不出哪一族來自父親、哪一族來自母親（那需要另外定序父母）；其二，標號只在該 phase block 內有定義，跨 block 的「家族一」並非同一條染色體。One of the two groups reads are split into by HP tag. Family 1 is HP1 plus the somatic haplotype HP1-1 derived from it; family 2 is HP2 plus HP2-1; HP3, which cannot be assigned to either germline haplotype, belongs to neither. It is called a family because it groups a germline haplotype together with the somatic haplotypes descended from it — the split is by the germline layer, not by whether a read carries a somatic mutation. Within one phase block a family corresponds to one chromosome copy, so the two families are the two homologous chromosomes at that locus. Two things do not follow: which family is paternal cannot be determined without sequencing the parents, and the labels are defined only within that phase block.完整條目 →，以下簡稱連鎖區段） 都算過那一批 read 的共現狀態，並判定過屬於串接或分岔； 所需的原料因此已經在手上，但目前停留在連鎖視窗層級，未被傳遞至全基因體尺度。

本頁為研究指引，預設讀者已完成上篇，內容以式子、參數與失效模式為主。 討論範圍限於單倍型。

概念與互動

一、VAF 路線的估計式

以下建立不使用任何長讀關係、僅以一維頻率譜重建的標準流程 —— PyClone、SciClone 與 MOBSTER 一系共同的骨架，共四個步驟。

第一步：將每個變異換算為細胞比例

位點 $m$ 的深度為 $d_m$，其中 $a_m$ 條 read 帶變異， ${\mathrm{VAF}}_m=a_m/d_m$。VAF 不可直接互相比較，因為它同時受三個因素稀釋： 正常細胞的混入、該處的拷貝數，以及變異在每個帶原細胞中的拷貝數。換算關係為：

此步驟需要三個外部估計值，且三者的尺度不同：ρ 為全樣本單一值， CN 為每區段一個值，μ 則為每個變異一個值。上篇所述的第一層不確定性即源於此。

第二步：將一萬五千個 c 匯總為一維直方圖並擬合混合模型

設共有 $K$ 群，細胞比例為 $c_1\dots c_K$； 每個變異帶有一個未觀測的群歸屬 $z_m\in \{1\dots K\}$，各群所佔的變異比例為 ${\pi }_k$。 分群即為下式的最大化：

第三步：決定 K

常用做法為 BIC、Dirichlet process，或如 MOBSTER 以冪次分布吸收中性演化的尾端。 此步驟為整條路線最脆弱的環節，因為它屬於模型選擇而非參數估計 —— 它所回答的是「有幾群」，而非「這幾群位於何處」。

第四步：由各群的細胞比例建樹

可用的約束只有一條，且屬純粹的算術： 一個細胞至多屬於某個父節點底下的一個子節點， 所以各子節點的細胞比例加起來不得超過父節點。形式化之後， 節點 $k$ 得為子節點集合 $S$ 的父節點，必要條件為：

關鍵在於它只能排除樹，不能挑出樹： 通過此檢查的候選通常不只一棵，其餘由簡約性parsimony 簡約法在所有與資料相容的解裡，選步數（或成本）最少的那一個。它是一個偏好而不是證據 —— 當多個解並列時，簡約法決定選哪一個，但資料本身沒有排除其他解。Choosing, among all solutions compatible with the data, the one requiring the fewest steps or lowest cost. It is a preference rather than evidence: when several solutions tie, parsimony picks one, but the data has not ruled the others out.完整條目 →等偏好決定。 上篇的測驗即針對此點：通過相容性檢查與被資料證成是兩回事。

三項限制與其對應的用途

二、接合點：多變異連鎖視窗的變異已包含在同一組觀測中

多變異連鎖視窗並非第二批資料。 其中的 1,475 個變異原已包含在前述 15,500 個點之內 —— 它們各自具備 $a_m$ 與 $d_m$，第一節第二步的式子照樣算得出來。 5% 所額外提供的不是新的觀測值。

但也正因為是同一批 read，接上去的方式不能是「再加一項」。 若同一條分子既貢獻了它在某個位點的 alt 計數， 又被換算成一條「這兩個變異不共存」的關係再算一次，那就是把同一個觀測計了兩次； 此時兩項的和不是一個 likelihood，而必須另外校準一個權重去稀釋重複的部分。

正確的接法是把觀測切開，而不是相加：

• 只含一個變異的連鎖視窗（14,025 個變異）：走第一節那條 beta-binomial，一個位點兩格。
• 含兩個以上變異的連鎖區段（1,475 個變異）：其 read 不再逐位點拆成邊際計數， 而是整批記成一張聯合的狀態表 —— 哪幾條分子帶哪一組字母。

兩類的 read 沒有交集，所以兩者相乘即為完整的 likelihood，不需要任何權重。 要點在於第二類已經包含第一類的資訊：一張聯合表沿任一個位點加總， 得到的就是該位點的 alt／ref 計數。所以多變異連鎖區段的變異不再單獨進直方圖 —— 不是被丟掉，是被它自己那張表接手了。

因此兩者的接合點不是某一個變數，而是整組全域參數 （$c$、$T$、$z$）：兩類觀測由同一組參數產生，只是觀測窗口不同 —— 一邊看到的是一個位點的邊際計數，另一邊看到的是好幾個位點的聯合計數。

第一節的式子把 $z_m$ 邊際化（即 ${\sum }_k{\pi }_k\dots$ 一項）， 所以它對「哪些變異落在同一批分子上」完全沒有話說 —— 直方圖只看得到每個變異各自的位置。 而聯合狀態表直接觀測到的正是 $z$ 的聯合行為，不是逐個變異的邊際。 這就是頻率譜原則上拿不到、而多變異連鎖區段拿得到的東西，也是以下三節的共同來源。

三者因此不是三種要另外加進去的約束，而是同一張表的三個推論。 下篇把這張表的機率逐層寫出來 —— 全域的 clone 如何投影成局部的分子比例， 一條只覆蓋部分位點的 read 如何計入，以及定序錯誤在哪一層進場。

三、用途一：使群數 K 可辨識

不可辨識性的代數形式

設真實情形為兩群 B 與 C，且細胞比例相等：$c_B=c_C=c$， 則兩者在 VAF 軸上的位置亦相等：${\theta }_B={\theta }_C=\rho c/2$。 代入第一節第二步的混合式：

兩者並非近似相等，而是恆等。似然僅透過 ${\pi }_B+{\pi }_C$ 依賴這兩個參數， 故 $({\pi }_B,{\pi }_C)$ 本身不可辨識。 又因 BIC 一類準則對額外參數施以懲罰，模型選擇必然取 K 較小者。

此即「增加定序深度無助於此」的精確意義： 資料量僅透過似然影響結論，而似然沿此方向為常數。 更換 caller、增加變異數目或提高深度，皆不作用於此方向。

數值例

設 ρ ＝ 0.60、二倍體、μ ＝ 1，真實情形為三群：

真實的樹為 A 分支至 B 與 C。但直方圖上只見 $0.30$ 與 $0.12$ 各一根峰， 最大似然解為 $K=2$，樹為 A → B 的線性形式。

此結果有兩處錯誤：群數少計一群，且拓撲由分支誤判為線性。 兩者皆無法由內部檢查發現 —— $K=2$ 對此組資料即為最大似然解， 殘差、擬合優度與後驗預測檢查均無異常。

今加入一個共現觀測：B 中的變異 $i$ 與 C 中的變異 $j$ 落於同一個連鎖區段，而該連鎖區段的狀態表中「同時帶有兩者」那一格是空的。

$K=2$ 的解在 0.12 一峰上僅有單一成分， 必然將 $i$ 與 $j$ 指派至同一群 —— 而同一群的分子本就該同時帶有兩者，於是那一格空著這件事在 $K=2$ 之下 機率極低。機率最高的解因而變成 $K=3$， 其中兩個成分位於同一位置 $\theta =0.12$。 pigeonhole 至此容許分支：$c_B+c_C=0.80\le c_A=1.00$。

因此，單一個共現觀測所改變的是樹的拓撲與群的個數，而非參數的數值。 這也是此路線的論據不在精度的原因：精度指「可區分，但量測是否準確」， 可辨識性指「資料原則上能否區分兩種情形」。前者隨資料量改善，後者不然。

此觀測的作用範圍

加入該觀測之前，該峰的 1,000 個變異全部被指派至同一群 —— 表面上全部有解，實際上是將兩群併為一群。 加入之後該峰須由兩個成分解釋，但兩個成分的 $\theta$ 相等， 故除了被釘住的 $i$ 與 $j$， 其餘 998 個對這兩個成分的似然完全相同，歸屬轉為不可辨識 —— 除非各自另有共現觀測連至 $i$ 或 $j$。

換言之，一個共現觀測所買到的全部位於群層級：

系統發生重建所需者正是左欄三項，而非右欄。 下篇「最稀疏的一對僅有 28 個連鎖區段」之所以足夠，理由即在於此： 28 個遠不足以指派 1,000 個變異，但對「這兩個成分是否相同」 此一是非判定則已然充足。

四、用途二：以量測值取代自由參數 φ

$\phi$ 在第一節中為與其他參數一併擬合的自由參數， 且與 $K$ 不可區辨：一個 φ 較大的寬成分， 與兩個 φ 較小、位置相鄰的窄成分，可擬合出近乎相同的直方圖。 擬合結果取決於懲罰項與初始值，而非取決於生物學。

5% 的連鎖視窗提供$\phi$的直接量測： 同一連鎖區段、同一條單倍型上的數個變異，其真實頻率完全相等 —— 它們位於同一批分子上，生物學離散恰為零， 故其間觀測到的離散度全部來自技術雜訊。 依深度與拷貝數分層估計，即得 $\widehat{\phi }(d,\mathrm{CN})$。

惟這個量測有一個必須寫明的界線：它給的是下界，不是 φ 本身。 理由正是上一句的優點反過來 —— 既然這幾個變異位於同一批分子上， 那麼「這批分子的真實比例偏離期望值多少」對它們而言是共有的， 兩者相減時會消去。剩下的只有各自的 read 抽樣與 base error。 而 $\phi$ 想捕捉的偏離（該變異的真實 CCF 略偏離群心、 該處 CN 或 μ 估錯、局部 mapping 偏誤）恰恰有一大部分是共有的那一種。

實作上的改動僅一處：將 $\phi$ 由待估參數改為代入的值。 其後果是模型不再能以「該群雜訊較大」吸收結構 —— 原先擬合為單一寬峰者，若 $\widehat{\phi }$ 小於擬合所得的 φ，即須改由兩個成分解釋， K 因而增加。此路徑與用途一殊途同歸，惟其作用點在雜訊一端。

但由於代入的是下界，「擬合值大於下界」本身還不足以斷定要多開一群 —— 那個差距也可能來自上表被消去的兩列。 正確的讀法是：下界把「該群雜訊較大」這個藉口的可用範圍限縮， 超出的部分必須被指名（是 CN 估錯、是群心偏離，還是真的多一群），不能默默吸收。 要真正定出 $\phi$，需要的是同群但分子池互相獨立（即相距夠遠）的變異， 而那組變異「真實頻率完全相等」的前提比本節弱得多。這是本路線尚未解決的問題之一。

五、用途三：以比值定出成分間距，使純度退化為整體縮放

第一節第一步的三個估計值，其偏誤的作用方式並不相同：

• ρ 偏誤：所有 $c$ 等比例縮放。樹的拓撲多半不變， 但 pigeonhole 為絕對比較（$c_B+c_C\le c_A$），判定結果仍可能翻轉。
• CN 或 μ 偏誤：其作用為逐變異的位移，受影響的變異將落入其他峰，汙染群的組成。 此類偏誤的後果較前者嚴重。

此項不是一條要另外加進估計式的等式，而是那張聯合狀態表的直接推論： 表裡的每一格本來就是「一組分子數除以全部分子數」，任兩格相除即得上式。 其效果是將 $K$ 個各自自由的位置參數，改寫為一組比值與一個整體縮放； $\rho$ 由「每個變異換算時皆須用到的因子」退化為「最後施加的單一純量」。 換算偏誤僅使各群整體平移，相對關係不變。

這也解釋了為什麼三個用途不需要三套權重去平衡： 它們讀的是同一張表的不同部分，而那張表只有一個機率。

預測與結果檢視

與下篇的銜接

本頁確立的是三個用途各自作用在估計式的哪一個參數上， 以及它們為何是同一張聯合狀態表的三個推論。 下篇把那張表的機率逐層寫出來：全域的 clone 如何投影成局部的分子比例、 拷貝數與純度在哪一層進場、一條只覆蓋部分位點的 read 如何計入、 以及定序錯誤與家族誤標各自造出什麼。

其中一項特別值得先預告，因為它決定整條路線站不站得住： 「這一格是空的」有三種完全不同的成因 —— 那個狀態存在但沒抽到、真的沒有任何 clone 帶它、 或它確實不存在但錯誤仍會造出看起來像它的 read。 三者的觀測一模一樣。混為一談的後果是把錯誤造出來的少數幾條 read 讀成一個新的 clone， 而報告上不會留下任何異常跡象。