下篇：完整的估計式，以及它能不能建起來

為什麼重要

中篇確立了三件事：僅以 VAF 頻率譜重建時， 群數 K 不可辨識、峰寬 φ 與 K 不可區辨、 成分位置 c 依賴三個外部估計值； 而多變異連鎖視窗的觀測正好各補一格，且三者作用在同一組待估量上： 突變屬於哪一群的指派 $z$、額外離散度 $\phi$，以及各群的 CCF $c$。

把三者合起來的最直接寫法，是把頻率的對數似然與一串「關係約束」的對數證據 相加，再給每條關係一個權重。這個寫法有一個結構缺陷，而且它不是細節： 多變異連鎖視窗的 read 會被計兩次 —— 一次是它們各自的 VAF 進了頻率項， 一次是同一批分子被換算成關係又進了第二項。 兩項因此不是獨立的資料集，那個和也就不是一個 likelihood； 權重之所以必須存在，正是為了事後稀釋這個重複。

本頁採取另一條路，而它其實更簡單：不要把 read 換算成關係，直接對 read 建模。 把觀測依連鎖視窗切成互不重疊的兩類 —— 只有一個變異的連鎖視窗，與兩個以上變異的連鎖視窗 —— 兩類各寫一個因子，相乘即得完整的 likelihood。 單變異那個因子正是原本的 beta-binomial； 多變異那個因子是同一個模型在多個位點上的一般形式。 兩者不是兩個模型；當視窗只含一個變異時，這個 read 模型的邊界情形就是 beta-binomial。

這個改動的代價與收益都很集中。收益是：權重、關係的可信度、以及為了保護硬約束而設的門檻 全部不需要了，因為證據強弱本來就是似然自己會算的東西。 代價是：觀測通道必須寫完整 —— 逐位點的定序錯誤與整條分子的家族誤標形狀不同， 少寫一種，錯誤造出來的少數幾條 read 就會被讀成一個新的 clone。

十二節分成四段，每一段回答一個不同層次的問題：

閱讀方式：主線全部寫在沒有摺疊的段落、式子與圖裡。 內文中的虛線摺疊區一律是補充 —— 符號的逐項查表、數字怎麼推出來的、 換一種實作寫法會怎樣、與前兩篇的接點。 第一遍可以整批略過，需要時再展開，主線不會因此斷掉。

概念與互動

一、從一個變異走到完整的目標函數

先不要從最後一條長公式讀起。這個 estimator 只是在反覆回答一個問題： 若假設有 $K$ 個突變群，哪一組群比例、clone tree 與突變指派，最能產生實際看到的 read，同時又沒有為了追逐雜訊而多開群？

在拆解之前先講一件會影響整節讀法的事：這一節要建立的是兩個估計問題，不是一個。

所以讀到第一條長公式時不必期待它就是最終目標 —— 它算的是「給定一組參數，資料出現的機率」， 要等到扣掉複雜度代價、加上群數先驗之後才成為可以拿來比較 $K$ 的分數。 兩層各自怎麼跑，第七節會講。以下先把內層拆成四層：

第一層：先分清楚觀測值、未知狀態與校準量

這一層只做一件事：把接下來會反覆出現的符號分成三類 —— 可觀測的（$a_m$、$d_m$）、要推論的（$K$、$z_m$、$c_k$）， 以及由上游餵進來或先校準好的（${\mathrm{CN}}_m$、${\mu }_m$、${\phi }_m$）。 $\rho$ 兩邊都不是 —— 它是推導量，理由在本節後半。

第二層：生物狀態先決定預期 VAF

分子問的是「樣本中預期有多少份帶突變的 DNA」；分母問的是 「這個位置總共有多少份 DNA」。因此 ${\theta }_{\mathrm{mk}}$ 不是新的觀測值，而是 在候選指派 $z_m=k$ 成立時，模型預測的 VAF 中心。

要記住的一句話是：同一群的變異不必落在同一個原始 VAF —— 拷貝數與 multiplicity 會把中心推到別的地方。

第三層：beta-binomial 把「相差多少」換成機率

這一層要的只有一句話：觀測到的 alt count 不會剛好落在 ${\theta }_{\mathrm{mk}}$ 上， 所以要用一個容許波動的分布去問「差這麼多有多不意外」。

因此 ${\phi }_m=0$ 是普通 binomial 的極限；${\phi }_m$ 越大，同一中心周圍可接受的 read-count 波動越大。它不是可任意調寬的「峰寬旋鈕」 —— 進模型的是一個先用獨立資料校準好的 $\widehat{\phi }(d_m,{\mathrm{CN}}_m)$，不是搜尋時可以順手調大的參數。

不要把這裡的 ${\phi }_m$ 與第三節的 ${\tau }_u$ 混在一起。 ${\phi }_m$ 描述跨獨立變異時單位點 alt count 的額外變異；${\tau }_u$ 則是多變異連鎖區段內 整個狀態比例向量的 concentration，而且本頁稍後主張 PCR-free 長讀資料預設不需要額外的 ${\tau }_u$。

第四層之一：全資料的 likelihood，讓每批 read 只出現一次

先講清楚「家族」是什麼

這個詞接下來會一直出現，所以先定義。 一條 單倍型家族單倍型家族 一條 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.完整條目 →就是依 HP tag 分出來的兩組 read 之一：

為什麼因子跑在視窗上，而不是跑在家族上

這是很容易寫錯的一步，而且寫錯不會有任何錯誤訊息。 直覺上，既然家族一與家族二的 read 各自形成一張 0／1 狀態表， 把兩張表各算一次機率再相乘似乎最自然。但那樣寫是錯的，有兩個各自獨立的理由：

因此本頁的規格是：資料表依家族分開建立，機率在視窗這一層合併。 兩條家族的狀態表都由同一組 ${\Theta }_K$ 投影出來， 再一起通過第二節那個成對的觀測通道，最後合成一個 $P_v$。

更完整地說，${\Theta }_K$ 包含細胞比例 $w$、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.完整條目 → $T$、 突變指派 $z$、逐拷貝基因型 $G$，以及觀測錯誤參數 $e$。 $c_k$ 由 $w$、$T$ 與 $z$ 決定，${\mu }_m$ 由 $G$ 決定。 $e$ 不是一個錯誤率，而是一組觀測通道參數的縮寫，包括稍後的逐位點錯誤 $\epsilon$ 與 家族誤標率 $\xi$。 外部提供的 copy number 與校準後的 $\widehat{\phi }$ 是條件輸入，為了避免式子過長沒有逐項寫在條件線右側。

$\rho$ 不在這張清單上。本頁的 $w$ 把正常細胞當成第 0 項 一起放進同一個 simplex（${\sum }_{g=0}^K{w}_g=1$，見第二節）， 所以純度是推導量 $\rho =1-w_0$，不是另一個自由參數。 連續自由度因此是 $K$ 個 —— $K=5$ 時是 5 個。

成立的關鍵不是「$S_1$ 與 $S_2$ 兩個名字不同」，而是給定 ${\Theta }_K$ 後，兩者背後的 read 觀測可視為條件獨立且沒有交集。 若一批 read 已進入 $D_v$，它就不能再透過該位點的 $a_m/d_m$ 進入第一個乘積。

要點在於 $S_2$ 的因子已經包含那些變異的邊際 VAF： 若視窗 $v$ 含 $k_v$ 個變異，一張 $2^{k_v}$ 格的聯合表沿任一個位點加總， 得到的就是該位點的 alt／ref 計數。 所以把多變異連鎖視窗的變異再放回第一個乘積，就是把同一批 read 算第二次。

第四層之二：目標函數 = 配適資料 − 複雜度代價 + 先驗偏好

複雜度項不可省略，因為多開一群幾乎總能改善資料配適，卻可能只是在吸收錯誤。 而它扣多少，取決於 $M_{\mathrm{eff}}$ —— 那是考慮同一突變叢與同一連鎖區段內相關性後的 有效獨立資訊量，不是把 read 或變異數直接代入。這個量估錯，模型選擇就整個歪掉。

還有一件尺度上的事要先講清楚：局部的比值只能定出 $c$ 的相對尺度， 不能憑空補出純度的絕對尺度 —— 絕對尺度只能來自 $S_1$ 那一批全基因體的單變異觀測， 或外部的純度估計。

此處寫的是可實作的統計規格，不是已完成或已驗證的 estimator 程式。 以下各項定義的目的，是讓後續實作能逐項測試；在合成資料與真實重複資料通過末節所列驗證之前， 不得把 $\mathrm{Score}$ 的最大值當成已證實的生物學結論。

符號多，不代表待估參數多

到這裡已經出現了二十幾個符號，很容易產生「要估的東西多到不可能估得準」的印象。 但大部分符號不是自由參數：有些是從別的量算出來的，有些是上游分析餵進來的， 有些要先用獨立資料校準好才進模型。把它們分清楚， 是判斷「這個結果可不可信」的第一步 —— 因為風險幾乎都不在自由參數上， 而在那些被當成已知、其實估錯了的量上。

兩個讀法上的要點。其一，真正在搜尋的只有五樣東西： 群數 $K$、clone tree $T$、突變指派 $z$、逐拷貝基因型 $G$， 以及一組 $K$ 個自由度的細胞比例 $w$。其餘不是算出來的，就是進模型前該備妥的。 其二，「不自由估」的那幾個才是這個方法的真正風險所在 —— $\mathrm{CN}$、$\phi$、$\epsilon$、${\xi }_v$、$M_{\mathrm{eff}}$ 估偏時模型不會報錯，它會照樣收斂，只是收斂到錯的 $K$ 與錯的樹。 第四節的錯誤地板與本頁末的問答，講的都是這件事。

因此這個 estimator 不是從原始 BAM 無條件地推出 clone tree。 它吃的是一套已經跑過 variant calling、copy-number calling、定相、單倍型標記 與錯誤校準的上游資料，而上述每一項都是它的輸入契約的一部分。 逐項該釘住哪些工具版本與門檻，是另一個層次的問題，本頁不展開。

二、一個多變異連鎖視窗的生成模型

上式的 $P_v(D_v\mid {\Theta }_K)$ 是整份規格的核心。它必須從全域參數一路長到一條 read 上看到的字母，中間不能有任何自由參數 —— 一旦局部比例可以自己亂動， 這個視窗就不再對全域參數施加任何限制，整條路線也就白做了。

先講清楚一件事：兩類觀測是同一條生成鏈的兩個出口

在展開細節之前要先擋掉一個很自然、但會把整件事讀歪的誤解： 單位點 VAF 不是多位點狀態表的先驗。 兩者都是觀測，都由同一組 ${\Theta }_K$ 沿同一條鏈生出來：

這解釋了 CCF 與 clone 比例講的是細胞、 而 HP1／HP2 的狀態表講的是分子，層級不同卻仍能相乘的原因： 它們最後都被換算到同一層 —— 兩個因子都在問「在這組 ${\Theta }_K$ 下，看到眼前這些 read 的機率是多少」。 所以能相乘靠的是兩件事，兩件都必須成立： 其一，細胞 → 分子的投影要正確（這是本節的工作）； 其二，沒有一條 read 同時進入兩個因子（這是第一節的工作）。

這條鏈預設了什麼

投影正確與否，取決於幾條沒有寫在式子裡、但整份規格都靠它們成立的假設。 把它們列出來，是因為違反時的後果各不相同，不列出來就沒辦法判斷某個結果可不可信：

前兩條在多數體細胞 SNV 上成立得相當好，這也是簡約類方法長期沿用它們的理由； 但它們是假設而不是結論，尤其在高突變負荷或 CN 劇烈變動的樣本上要另行檢查。 第四節那個「潛在節點的三種身分」正是這條假設在輸出端的表現形式。

連鎖區段、局部基因型與分子權重

一個連鎖區段 $u$ 是一個單倍型連鎖區段：一個連鎖視窗 × 一條 單倍型家族單倍型家族 一條 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 全部來自同一條染色體拷貝， 另一條拷貝屬於同一個視窗的另一條連鎖區段。因此 不可以在連鎖區段內部再把兩條拷貝混起來寫成各佔一半 —— 那是尚未依家族分組時的寫法，用在已分組的資料上會把細胞比例整體算錯將近一倍。

要同時記住兩件看似相反、其實各管一段的事，第一節那張表已經寫過一次： 狀態表逐家族分開建立（$u$ 這一層），機率在視窗合併（$v$ 這一層）。 以下先在單一家族內把 $q_u$ 建起來，本節末尾再把兩條家族接回同一個 $P_v$。

設 clone $g$ 在拷貝 $b$ 上的全域基因型為 $G_{g,b,m}\in \left\{\mathrm{0,1}\right\}$， 投影到連鎖區段 $u$ 的那幾個位點，即得局部基因型 $H_{\mathrm{ugb}}\in \{\mathrm{0,1}{\}}^{k_v}$。基因型寫在拷貝這一層而不是 clone 這一層， 是為了讓 cis／trans 與 LOH 有地方表達。分子權重則為：

$A_{\mathrm{ugb}}$ 是該拷貝在此處的份數，${\mathrm{CN}}_{\mathrm{ug}}$ 是 clone $g$ 在此處的總拷貝數。 分母是全部分子數，所以 ${\sum }_{g,b}{\eta }_{\mathrm{ugb}}=1$。 拷貝數在這裡進入，不是事後乘上去的；LOH 時兩條同源染色體不再各佔一半， 這條式子自動處理。正常細胞是 $w_0$ 那一項，$\mathrm{CN}=2$， 其局部基因型恆為全參考。

要注意 $q_{u,h}$ 是分子比例而非細胞比例，兩者差一個拷貝數因子。

覆蓋遮罩：沒覆蓋到的位點要邊際化

一條 read 只覆蓋連鎖區段裡的一部分位點。設其覆蓋遮罩為 $M\subseteq J_u$、觀測到的字母為 $y\in \{\mathrm{0,1}{\}}^{\left|M\right|}$，則

也就是把所有「在 $M$ 上與 $y$ 相符」的完整狀態加起來。 一條只看到第一與第三個位點、讀到 $1\text{–}1$ 的 read，同時支持 $101$ 與 $111$。沒覆蓋到的位點不可以填成參考型 —— 那會把一條中立的 read 變成一條反對 $111$ 的證據。

一個連鎖區段一個 p，不是一個遮罩一個

同一個連鎖區段裡不同遮罩的 read 取樣的是同一池分子， 所以它們的比例是相關的。正確的寫法是讓整個連鎖區段共用一個潛在比例向量， 每條 read 各自從它自己的遮罩投影出來：

觀測通道：兩種形狀不同的錯誤

$q_u$ 是「還沒經過任何錯誤」的比例。實際觀測到的比例要再過一層：

兩種錯誤不是同一類參數：一個能先校準，一個不能

既然定序平台固定，錯誤率是不是就可以當成已知常數餵進來？ 對一半，而這一半正好是比較不重要的那一半。兩個通道要分開處理：

這個區別有實際後果，而且方向明確：$k_v\ge 2$ 時錯誤地板幾乎完全由 ${\xi }_v$ 決定，因為 ${\epsilon }^r$ 隨位點數急速掉下去而 ${\xi }_v$ 不會。 所以「把錯誤率當成已知」若指的是 $\epsilon$，大致無害； 若指的是 ${\xi }_v$，那正是第四節錯誤地板訂錯的主要途徑 —— 也就是把幾條誤標的分子讀成一個新 clone 的那條路。 ${\xi }_v$ 尤其要成對估計，因為它同時連結一個視窗的兩條家族。

把兩條家族接回同一個因子

交換閘門的存在，表示一條被標成家族 $u$ 的 read 不一定真的來自家族 $u$。所以真正用來算機率的觀測分布， 必須是「正確標記」與「由隔壁誤標進來」兩個來源的混合。 令 $u_1$、$u_2$ 為視窗 $v$ 的兩條家族， ${\alpha }_{u\text{'}\to u}$ 為「一條最後被標成 $u$ 的 read，實際來自家族 $u\text{'}$」的權重：

HP1／HP2 的方向是 nuisance，而它屬於視窗

家族標號在每個 phase set 內獨立決定， 所以「這個視窗的家族一」對應到全域基因型的哪一條拷貝 $b$ 並不確定。 處理方式是引入一個方向變數 ${\sigma }_v\in \{原樣,對調\}$ 並把它邊際化。 ${\sigma }_v$ 是每個連鎖視窗一個，不是每條家族一個 —— 同一個視窗的兩條家族互為補集，對調就是同時對調，共用同一個方向：

三、預設是 multinomial：一條 read 是一個分子

上式寫了一個 ${\tau }_u$，但它的預設值應該是無窮大，也就是退化成 multinomial。理由是 guardrail 第 7 條的直接後果：

beta-binomial 的 overdispersion 來自「分子池的真實比例偏離期望值」—— 那是一個兩段抽樣的結構：先從細胞抽出分子，再從分子抽出 read。 而 PCR-free 的長讀定序沒有第二段：一條 read 就是一條分子， 不存在同一條分子被讀到兩次的情形。細胞數又以百萬計， 第一段抽樣的變異小到可以忽略。所以在一個連鎖區段內部， read 的計數本來就是 multinomial，額外的離散度沒有來源。

這給出一條明確的實作規則：先用 multinomial，只有在殘差真的顯示額外離散時才加 ${\tau }_u$，而且要說得出它來自哪裡（連鎖視窗內部的 mapping 或覆蓋不均， 而非分子抽樣）。${\tau }_u$ 若要估，也必須跨可比較的連鎖區段分層借力 —— 一張稀疏的計數表沒有能力同時定出 $q_u$ 與它自己的離散度。

四、三種不同的「沒有」

這一節是整份規格裡最容易出事的一節：狀態表上一格是 0， 可能代表三件完全不同的事，而三者的觀測一模一樣。 分不開它們，錯誤造出來的少數幾條 read 就會被讀成一個新的 clone。

第二種與第三種必須分開寫，還有一個純技術但無法迴避的理由： Dirichlet 的參數不可以是零。結構零給出 $q_h=0$， 而 $\tau q_h=0$ 不是一個合法的 Dirichlet 分量。 所以結構零一定要先經過錯誤通道抬到地板之上：

這不是為了數值穩定而加的小常數。 它就是「零不代表不可能」這句話的數學形式 —— ${\delta }_v{r}_{u,h}$ 的大小正是這個連鎖區段能分辨的最小比例， 也就是第十一節那個偵測下界的來源。若 $r_{u,h}$ 由第二節的兩個錯誤通道算出來， 這個地板是可以就地估計的量，不需要另設參數。

這一節的取捨在整份規格裡最容易被低估，所以本頁末的「預測與結果檢視」 用一整個問答把它再走一次：門檻確實不需要了，但錯誤地板變得比以前更重要。

潛在節點的三種身分

觀測到 $000$、$010$、$111$ 而 $011$ 與 $110$ 皆空時， 最小成本解會補進一個中間狀態。這個節點的身分必須標清楚，因為它有三種可能， 而三者的意義完全不同：

五、分類是輸出，不是輸入

現行流程會把每個連鎖區段標成分岔、串接或不相容。 這些標籤在新寫法下照樣產生，但不再是餵進模型的證據 —— 它們只是狀態表本身的形狀。這一節說明為什麼，以及因此有三個問題直接消失。

現行流程對每個連鎖區段所做的分類（分岔／串接／不相容）在新的寫法下照樣產生， 但它們的地位變了：它們是 $q_u$ 的形狀，不是額外的證據。

由此得到一條必須寫進實作的規則： 凡其原始狀態表已經進入似然的連鎖區段，不得再為它加一條約束項。 分岔摘要仍可以印出來給人看、畫成圖、當作 QC，但那是衍生輸出， 不是另一份獨立證據。這條規則就是第一節那個「兩個因子相乘」在實作層的說法。

三個原本需要特別處理的問題，也在這裡一併消失：

六、數百個連鎖區段的局部形狀不同 —— 這是預期，不是矛盾

通過品質檢查的連鎖區段有數百個，而它們的局部拓撲形狀各不相同： 有的是一條線，有的是分岔，有的只有一個節點。此現象本身不構成衝突。

成因是每個連鎖區段只裝得下落在該連鎖視窗內的那幾個變異。 上圖的連鎖區段甲裝到的是 B 群與 C 群的變異，兩群互不包含，故其局部形狀是分岔； 連鎖區段乙裝到的是 B 群與 D 群的變異，後者巢狀於前者，故其局部形狀是一條線。 兩者皆為同一棵全域 clone tree 在不同變異子集上的投影 —— 形狀相異是因為子集相異，而非證據衝突。

合起來為什麼更準、解析度更高

關鍵在於認清頻率譜的輸出並不是一棵樹，而是一棵部分決定的樹。 本例中它換算出三群（$1.00/0.40/0.20$），而這三群在 pigeonholepigeonhole 鴿籠原理由父代與子代的細胞比例限制樹形的算術規則：一個細胞至多屬於父節點底下的一個子節點，所以各子節點的細胞比例加起來不得超過父節點。名稱來自鴿籠原理 —— 東西放進籠子，總量不會憑空變多。在 subclone 重建的文獻中也稱為 sum rule 或 crossing rule。它只能排除樹，不能挑出樹：通過檢查的候選通常仍不只一棵，其餘要靠 parsimony 之類的偏好決定。The arithmetic constraint that limits tree shape from parent and child cell fractions: a cell belongs to at most one child of a given parent, so the children's cell fractions cannot sum to more than the parent's. The name comes from the pigeonhole principle. Also called the sum rule or crossing rule in the subclonal reconstruction literature. It can only rule trees out, never select one: the surviving candidates are usually more than one, and the choice among them falls to a preference such as parsimony.完整條目 → 之下 有兩棵樹同時通過：線性的 A → X → D，與分支的 A → $\{X,D\}$。 取兩棵的交集，才是頻率譜真正確定下來的東西：

兩個連鎖區段補的正是右欄那兩格，而且各補一格： 甲的狀態表在「兩個突變同時出現」那一格只有錯誤地板的量，乙則相反。 候選數因此逐步收斂 —— 這就是「合起來更好」的具體內容：

三項改變的性質並不相同，值得分開看：

第二項尤其值得看清楚：「D 掛在 B 之下還是 C 之下」在頻率上永遠沒有答案。 B 與 C 的期望 VAF 完全相同，所以兩種掛法產生一模一樣的直方圖 —— 這不是精度問題，是中篇第三節那個不可辨識性換一個位置再出現一次。 連鎖區段乙的共現觀測是唯一回答得了的證據。

而三項改變都不需要先把兩個連鎖區段對齊： 兩個連鎖區段各自算出自己那張表的機率，而兩張表都由同一組 $(w,T,z)$ 產生，相乘即完成合併。 形狀不同從頭到尾都不是需要處理的問題，而合併也不需要任何對齊步驟。

由此導出一條實作規則，與第九節那條互為表裡： 可以合併的是機率，不是標籤。 局部標籤（HP1-1）僅在該連鎖視窗內有定義，跨連鎖視窗的同名標籤並非同一條， 故不得將局部拓撲拼接成一棵大樹。所要合併者是它們對 $(w,T,z)$ 的支持度 —— 而這組參數本就只有一份。

真正的矛盾長什麼樣

真正的矛盾並非形狀相異，而是整組連鎖區段找不到任何一組 $(w,T,z)$ 能同時給出夠高的機率。 在似然寫法下它不再讓求解器卡住 —— 因為沒有任何一格的機率是 0， 求解器會取總機率最大的那組參數，少數異常的連鎖區段被多數蓋過。 但它仍然是必須報告的診斷量：擬合後各連鎖區段的狀態表與模型預測的差距， 估計的是整條產線的錯誤率。

惟此處有一個陷阱：殘差小不代表模型正確。 系統性偏誤（例如第八節的排序偏誤）會使錯誤彼此一致 —— 它們朝同一方向偏，因而完全不互相矛盾，卻一併使結論偏移。 殘差只偵測得到隨機錯誤。

七、推論怎麼跑：外層定結構，內層估比例

到這裡為止，寫下來的都是「什麼是好的解」。這一節講「怎麼找到它」。 先講結論，因為它跟多數人的預期不同：真正難的是離散的那一半， 而連續的那一半小到不需要什麼特別的工具。

兩層的分工

內層：$K$ 個自由度，直接解就好

結構一旦固定，$c_k$ 由 $(w,T,z)$ 決定、${\mu }_m$ 由 $G$ 決定、 $\theta$ 與 $q$ 都是 $w$ 的函數，錯誤參數已經校準或就地估好。 剩下的自由參數只有 $w$ 本身 —— 依第一節的清點，那是 $K$ 個自由度 （$K=5$ 時是 5 個）。問題因此化簡成：

這就是一個 constrained maximum-likelihood 問題： 在「每項非負、總和為一」這個定義域（一個 simplex）上找最大值。 「直接解」不是說有閉式解 —— 一般沒有 —— 而是說把目標函數與這兩條限制交給數值最佳化器即可。

這裡不需要 EM，也不需要 MCMC。這一點值得說清楚， 因為看到 likelihood 就想到 EM／MCMC 是很自然的反射，但在這個內層問題上兩者都是殺雞用牛刀： 五個參數的凸定義域最佳化，數值最佳化器幾毫秒就收斂。

但有一個必須誠實標註的保留：這個內層問題不保證是凹的。 定義域確實是凸的，可是 $q_{u,h}$ 是 $w$ 的比值函數 （第二節那條 $\eta$ 的分母也含 $w$），而且狀態表的機率是多個狀態的加權和 —— 兩者都會破壞對數似然的凹性。實務上必須多起點並記錄各起點是否收斂到同一點； 只跑一個起點就宣稱找到最大值，是這一步最容易犯的錯。

外層：不可能窮舉，所以重點是怎麼不列舉

外層才是成本所在。固定 $K$ 之後，樹的形狀、 每個突變掛到哪一群、每個位點的拷貝基因型，全都是離散選擇，數量隨 $K$ 迅速上升 —— 光是有根有標號的樹就已經多到不可能一棵一棵試。所以概念上的骨架雖然是：

for K = 1 ... Kmax:
    列舉／搜尋這個 K 下的離散結構候選
    for each 候選:
        固定結構 → 在 simplex 上最大化 log-likelihood   ← 內層，直接解
        記下這個候選的最佳 log-likelihood
    取這個 K 的最佳值
比較各 K 的 Score(K, Θ̂_K)，取最大者為 K̂

但中間那個「列舉／搜尋」不是一個真的 for 迴圈。 可用的做法包括啟發式搜尋、分支定界、以取樣取代列舉，以及最重要的一項： 用局部證據先把候選砍掉。第八節那個候選集合正是為此而存在 —— 每個連鎖視窗算出來的局部候選整組帶進外層， 不相容的全域結構因而在列舉之前就被排除， 這比列舉完再逐一評分便宜得多。

兩件事要連著記住：其一，$K$ 的那層迴圈是真的迴圈， 而且必須跑完再比 $\mathrm{Score}$，不能中途看 $L_{\mathrm{data}}$ 上升就停 —— $L_{\mathrm{data}}$ 對 $K$ 幾乎單調遞增，會停不下來。 其二，這一整節是本規格中唯一沒有閉式成本估計的部分， 也是第十二節把它排在最後幾階段的理由。

八、現行實作要改的兩處

前面七節寫的是規格。落到現有的程式碼上， 真正要動的只有兩個地方，而且兩處都是「把資訊留下來」而不是「加一道閘門」： 其一，覆蓋遮罩不要併掉；其二，候選排序的分數換成似然。 超立方體列舉、分類輸出、前處理全部沿用。

其一：覆蓋遮罩要留下來，而不是併掉

一個連鎖區段的有效樣本數並非其總 read 數， 而是各位點組合分別有幾條 read 同時覆蓋。 上圖連鎖區段的三個位點對，跨越深度分別為 20、17 與 3。

舊的做法是設一道門檻把跨越深度不足的丟掉 —— 因為在硬約束的寫法下，一條由 $n=3$ 產生的假「不可同群」代價極高。 在似然的寫法下這道門檻不需要了：$n=3$ 的觀測自然只貢獻很小的證據量， $(1-q{)}^3\approx 0.34$ 這個數字就是似然自己會算的東西。

所以這一步的改動不是加一道閘門，而是把遮罩存下來： 狀態表要記成 $(連鎖區段,\mathrm{phase}\mathrm{set},位點集合,覆蓋遮罩,字母樣式,計數)$ 的稀疏列， 而不是先把部分覆蓋的 read 併進某個完整狀態。 門檻退居為報告量：收得幾條、各遮罩各有幾條、哪些連鎖區段的證據其實很薄。

for each unit u:
    for each read r spanning ≥1 site in u:
        emit (u, phase_set, sites(u), mask(r), pattern(r), 1)
    aggregate → sparse rows
report: 每個遮罩的 read 數、連鎖區段的最大跨越深度、
        以及「若某狀態真佔 5%，這個連鎖區段看得到它的機率」

可跨越範圍受限的成因，來自以下三層結構：

其二：候選排序的分數要換成似然

最小成本候選常不只一個，現行流程以 read-AFread-AF read 層級的 ALT 比例在同一個分析區域、同一個單倍型家族的 read 之中，某個 somatic 位點帶 ALT 的比例。分母已限縮到同一條 haplotype 的同一段區域，所以不必經過純度與拷貝數換算；用途是在步數並列的候選樹之間排序。The fraction of reads carrying the ALT allele at a somatic site, computed within one analysis window and one haplotype family. Because the denominator is already restricted, no purity or copy-number correction is needed; it is used to rank candidate trees of equal cost.完整條目 → 的差值總和排序； 問題在於相加的項數由拓撲形狀決定：

其後果為：本應呈現分岔的連鎖區段，將改為呈現串接。 而分岔是拆峰能力唯一的來源，因此該偏誤系統性地移除拆峰的能力， 使中篇第三節所述的 K ＝ 2 得以維持。

修正方式正是第二節那個 $P_v$：在該候選所施加的結構下算出 $q_u$， 再算整張狀態表的似然。如此不同形狀所比較的是對同一組計數的解釋能力， 形狀本身不再決定勝負。而且不必只留一個候選 —— 局部候選集合應整組帶進全域推論，由全域參數決定哪一個站得住。 這是把「先挑一個再往下走」換成「全部留著，讓似然決定」。

現行實作的超立方體列舉照樣沿用，它負責產生候選與必要的潛在節點； 被換掉的只有那個分數。這是整份規格中最小、也最值得先做的一個改動。

九、兩類與相位有關的誤讀

兩者皆與相位有關，方向相反： 其一為相位資訊僅在 block 內有效，跨 block 即無定義； 其二為即使在 block 內相位解得正確，仍可能被解讀為錯誤的生物學。

十、連鎖區段預算

結論先講：在最不利的低突變負荷樣本（5／Mb）上， 全基因體約產出 369 個同家族的多變異連鎖區段； 若腫瘤有三群、質量 0.60／0.25／0.15，最稀疏的那一對群仍分到 28 個。 這個數量看起來很小，但它要回答的是一個是非題（這兩個成分是不是同一個）， 不是要指派一萬五千個點 —— 所以夠用。

依上述模型，若腫瘤含三群、質量分別為 0.60／0.25／0.15， 則跨第 i 群與第 j 群的連鎖區段數為總數乘以 2w_iw_j：

充足性的依據已見於中篇第三節：這些連鎖區段的作用在於確立成分存在且相異， 而非指派 15,500 個點。最稀疏的一對有 28 個獨立連鎖區段， 對「兩個成分是否相同」此一是非判定為充足。

惟這個預算對第十二節的分期很敏感。 若初期只納入 copy-number-neutral 的雜合區域（那是分期的建議做法）， 在 CN 變異廣泛的實體腫瘤中可能只剩一半甚至更少， 最稀疏的那一對會從 28 掉到十幾個 —— 那時充足性的論證就不成立了。 因此分期的第一階段應優先選高突變負荷的樣本； 低負荷樣本要等到允許非中性區域之後才有足夠的連鎖區段。 報告須載明各階段實際納入與排除的連鎖區段數。

十一、連鎖視窗寬度：連鎖區段產量與品質的取捨

視窗開得寬，收得到的位點多，但多數 read 只覆蓋其中一部分，狀態表變得很稀疏； 視窗開得窄則相反。結論是兩種都做，因為在似然的寫法下同時用兩種寬度不需要付任何代價 —— 兩批都進同一個乘積，各自的證據量由各自的計數與遮罩決定。

這是與舊寫法最直接的對照 —— 舊寫法必須替兩批各設一個 ${\omega }_r$， 新寫法什麼都不必設。

十二、實作分期

整份規格不應一次到位。以下順序讓每一階段都能單獨驗證， 而且前一階段的結果可以當作後一階段的初始值：

真實資料與證據

現行實作的量測結果：可用連鎖區段的庫存

以下數值皆來自本實驗室既有實作在七個資料集、六個生物樣本、chr1–22 上的輸出， 在本頁的脈絡下構成可寫成 $P_v$ 的原料清單。

三項須注意的事項

其一，0.909 為目前最強的結果。一個在更換 basecaller 後仍維持、 且能區分癌別的統計量，所承載的是真實訊號。改成似然之後應重新計算此數值： 若維持或提升，表示新的分數確實在解釋同一個訊號；若下降， 則原本的 0.909 有一部分係由共同的覆蓋幾何所致。

其二，「分支類」的比例尚不可直接視為證據產量。 分岔的判定需要兩個單突變狀態皆存在，且雙突變狀態的缺席不是抽樣零也不是錯誤地板 —— 後者正是第四節所檢驗的對象。故 6%–22% 為上限，而非實際可用量。

其三，有一項計數需要釐清。 各資料集的連鎖區段數總和為 71,955，而結論中引用的數值為 85,941。 兩者可能計數不同的對象，此點須在成文之前確認。

預測與結果檢視

尚未解決的問題

• 整條路線尚未經實測驗證。上述各項皆有結構性的理由，但均未在真實資料上檢驗。 最小的驗證即中篇第三節的數值例：合成 ρ ＝ 0.60、CCF 分別為 1.00／0.40／0.40 的三群， 確認僅用 VAF 的分群必然給出 $K=2$ 與線性樹， 再量測需要多少個多變異連鎖區段、錯誤率要估得多準，方能正確分離而不誤拆。
• $\phi$ 的量測方式需要重新設計。第三節指出： 用同一批分子上的變異估 $\phi$ 會系統性偏低，而偏低的方向恰好偏向本路線想要的結論。 改由相距夠遠的同群變異估計之後，$\widehat{\phi }$ 會變大多少是一個必須先回答的問題 —— 它直接決定用途二還剩多少效力。
• 推論的計算成本尚未評估。全域取樣的內層包了每個連鎖區段的候選集合與 ${\sigma }_v$ 加總，而連鎖視窗有數萬個。這是整份規格中唯一沒有閉式估計的部分。
• 內層最佳化的凹性沒有保證。第七節指出， $q$ 對 $w$ 是比值函數、狀態機率又是加權和，兩者都會破壞對數似然的凹性。 因此「找到的是全域最大值」目前只是假設： 需要先用合成資料量出多起點之間收斂到不同解的比例， 才知道多起點要跑幾個、以及區間估計要不要因此加寬。
• 連鎖區段的來源具有偏性。它們僅來自變異密集的連鎖視窗， 而變異密集有時並非偶然：kataegis 為一種局部超突變， 單次事件於數十 kb 內產生一連串變異。 該串變異幾乎必然同時發生且同屬一群，故將貢獻大量高度相關的證據， 在統計上卻被計為數十個獨立觀測。應先偵測並降權為每叢一個有效觀測。
• 一個群未必對應單一細胞群體。本實驗室的甲基化觀測顯示： 單一個「單倍型 ＋ 等位」狀態之下仍可區分出五群甲基化模式。 故 $K$ 的正確讀法恆為「至少此數」。
• 超出 read 跨距的關係仍不可觀測。相距 100 kb 的兩個變異無法落入同一個連鎖區段， 與定序深度無關；第八節的三層結構為此設下硬上限。

原始文獻與程式碼

本頁「現行實作的量測結果」一節的數值均出自： 廖子游，Subclonal reconstruction using somatic haplotagging and methylation profiles with Nanopore sequencing，國立中正大學資訊工程學系碩士學位口試，2026 年 7 月 30 日。

以 read 上的多變異狀態計數（而非單點 VAF）直接建模 subclone， 最接近的既有做法是 PairClone 與 TreeClone：Zhou T, Sengupta S, Müller P, Ji Y, PairClone: a Bayesian subclone caller based on mutation pairs， J R Stat Soc Ser C 2019;68:705–725 （academic.oup.com）； 以及 Zhou T et al., TreeClone: Reconstruction of tumor subclone phylogeny based on genotype pairs， arXiv:1703.03853。 兩者處理的是成對的位點與部分觀測的 read；本頁規格與其相異之處在於 可變的 $k$、逐 phase set 的方向邊際化，以及與全基因體單變異連鎖視窗的相乘接合。 此一相異性尚未經文獻檢索確認，成文前應先查證。

單變異連鎖視窗那個因子所屬路線的代表性實作：SciClone，Miller CA et al., PLoS Comput Biol 2014;10:e1003665；PyClone，Roth A et al., Nat Methods 2014;11:396–398； MOBSTER，Caravagna G et al., Nat Genet 2020;52:898–907 （doi 10.1038/s41588-020-0675-5）， 程式碼在 github.com/caravagnalab/mobster。 換算式與檢定力指標見 Tarabichi M et al., Nat Methods 2021;18:144–155 （doi 10.1038/s41592-020-01013-2）； 評比數據見 Salcedo A et al., Nat Biotechnol 2024 （doi 10.1038/s41587-024-02250-y）。

簡約假設的原始出處為 Camin JH, Sokal RR, A method for deducing branching sequences in phylogeny, Evolution 1965;19:311–326； 有向最小成本 Steiner 樹形圖的標準參考為 Charikar M et al., Approximation algorithms for directed Steiner problems, J Algorithms 1999;33:73–91。 第八節所沿用的候選列舉即屬此系。

前處理所依賴的單倍型標記見 LongPhase-S 預印本 （bioRxiv, 2025，doi 10.1101/2025.11.20.689492）， 程式碼在 github.com/CCU-Bioinformatics-Lab/longphase-s。 參數名稱與欄位一律以原始碼為準。