研究指引 · Subclone 系統發生重建 · Supplement — a local estimand and features for somatic callers

補篇:把目標換成局部系譜,以及它能餵給 caller 什麼

不重建全域細胞層級的 clone tree,改為重建每個單位內部的分子系譜:先用證據預算說明局部結構為何走不進以邊際 VAF 為觀測的全域似然,再把前處理寫成契約(工具鏈門檻、HP3、視窗為讀序連鎖的傳遞閉包、逐樣本連鎖率),然後是局部生成模型、全域頻率譜如何以單向先驗進入而不重複計數、可計數的輸出旗標與其虛無分布,最後是給 ClairS/DeepSomatic 的三層特徵、留一法與分層驗收。

建議先修:估計式與落地

本模組學習目標

  • 把估計目標由全域細胞層級的 clone tree 換成連鎖區段內的分子系譜,並寫出這個目標的形式定義
  • 用證據預算說明局部結構為何走不進以邊際 VAF 為觀測的全域似然,而不是「全域方法不夠準」
  • 寫出局部生成模型:狀態集合、四配子相容性、節點分子比例與兩個錯誤通道
  • 把前處理寫成契約:工具鏈與門檻、HP3 的三種處理、連鎖視窗為讀序連鎖的傳遞閉包、以及逐樣本先算連鎖率
  • 說明全域頻率譜如何以單向先驗進入局部推論,而不把同一批 read 算兩次
  • 定義輸出語義與旗標,並給出 REFINES_GLOBAL 這個可計數、且有虛無分布可比的產出
  • 規格化給 somatic variant caller 的三組特徵、留一法、缺席契約與分層驗收
  • 列出驗收條件、實作順序,以及這條路線不可宣稱的事

為什麼重要

換掉的是目標,不是方法

前一篇把兩類觀測寫成兩個相乘的因子,接上同一組全域參數 (w,T,z),目標仍然是一棵全基因體、細胞層級的 。 這個目標有一個現實問題:以邊際 VAF 與 為觀測的既有方法已經成熟且被系統性評比過, 要在同一個目標上勝過它們,需要的證據量遠大於多變異連鎖視窗能提供的。 中篇的計算已經指出,可用的多變異連鎖區段在低突變負荷樣本只有數百個。

本篇因此換一個估計目標:不重建全域的細胞層級系統發生, 而是重建每一個連鎖區段內部的分子系譜 —— 哪幾種局部狀態真的存在、 各佔多少分子、以及它們之間的先後關係。這個目標的解析度是單倍型連鎖區段 —— 一個連鎖視窗(read 連鎖起來的一段,不跨 )的一條染色體拷貝,不是「全基因體共用的一個 CCF 位置」。

同一批變異,做成一棵全域樹,或做成幾百棵局部小樹 由上而下三個面板,說明估計目標的改動。 最上方是同一批觀測:一條全基因體的長條,上面散布著體細胞變異; 其中四段用虛線框起來,那是有讀序把兩個以上變異連鎖起來的連鎖視窗,依序標為 R1 到 R4; 框外還有許多變異沒有被任何讀序連鎖起來。 中間面板是既有做法:把全基因體所有變異匯集成每個變異的邊際等位頻率, 輸出一棵細胞層級的樹。這棵樹的一個節點是一群細胞, 解析度的連鎖區段是一個全基因體共用的細胞比例位置, 所以細胞比例相同的兩支 lineage 必然被併成同一個節點。 最下方是本規格:每一個框起來的連鎖視窗,各自依兩條染色體拷貝再切成兩個連鎖區段, 每個連鎖區段各輸出一棵小樹,因此圖上有八棵彼此獨立的小樹。 八棵的形狀各不相同:有的只有一條邊,有的是一條鏈,有的分岔, 而且同一個連鎖視窗的兩個家族可以是不同的形狀 —— 這是資料決定的輸出,不是設定。 沒有被連鎖起來的變異不產生任何小樹。 最下方的結論是:兩邊看的是同一批讀序,差別只在把它們摘要成什麼。 左邊是每個變異的邊際計數,右邊是同一條分子上的共現。 而那些小樹不接成一棵大樹,因為局部標籤只在連鎖區段內有定義。 同一批變異:匯集成一棵樹,或留成幾百棵小樹 同一批 read、同一批體細胞變異 R₁ R₂ R₃ R₄ 既有做法:全部匯集成邊際 VAF,輸出一棵樹 全基因體所有變異的邊際計數匯在一起 位置資訊在匯集時就沒了 一個節點 = 一群細胞 解析度:全基因體共用的 CCF CCF 相同的兩支必然併成一個節點 本規格:每個連鎖視窗 × 每條染色體拷貝各一棵 R₁ R₂ R₃ R₄ HP1 HP2 形狀各不相同,同一個連鎖視窗的兩個家族也可以不同;沒有連鎖的變異不產生小樹。 兩邊看的是同一批 read,差別只在摘要成什麼。 左是每個變異的邊際計數;右是同一條分子上的共現。小樹不接成一棵大樹 —— 標籤只在連鎖區段內有定義。
最上方是同一批 read 與同一批變異,下方兩個面板是它的兩種去向。 中間面板把所有變異匯集成邊際 VAF,輸出一棵樹; 最下方讓每個連鎖視窗、每條染色體拷貝各留一棵小樹 —— 八棵形狀各不相同,而且同一個連鎖視窗的兩個家族可以不同。 沒有被連鎖起來的變異不產生小樹,那正是後文「輸入與前處理契約」定義的連鎖率。
 既有做法的目標本規格的目標
估計對象全基因體的細胞群數、比例與樹每個連鎖區段內各局部狀態的比例與先後
觀測進入似然的形式每個變異的邊際 VAF + 拷貝數同一條分子上多個位點的共現
一個節點是什麼一群細胞一類分子狀態
解析度一個全基因體共用的 CCF 位置一個連鎖視窗、一條染色體拷貝
輸出幾棵樹一棵幾百到幾萬棵,彼此不相接
永遠分不出來的CCF 相同的兩支 lineage跨越 read 距離以外的任何關係
多分出一支要多少證據約 20 個變異(兩群差 0.05 時);兩群越近需要越多數十條跨越分子;與兩支靠得多近無關

兩條路線的差別集中在最後一列: 頻率譜要分開兩支,靠的是它們的比例不同;共現要分開兩支,靠的是它們的突變組合不同。 所以「兩支比例相同」對前者是致命的,對後者完全沒有影響 —— 只要它們帶的突變不一樣,分子上就看得出來。

局部這一側的限制在別的地方:它要求每一支各自夠多 (高於該連鎖區段的錯誤地板),而不是要求兩支之間差得夠開。 這個限制由後文「候選模型與輸出旗標」的 FLOOR_LIMITED 旗標承接。

交付物是「成對關係目錄」,不是幾百棵漂亮的樹;而本頁是規格,不是已驗證的方法

既有實作實測到,七個資料集中有六個,最大的一類是「一個連鎖區段只有兩個變異」 (完整分布見「輸入與前處理契約」)。所以這條路線的主要產出是全基因體的成對分子關係目錄 —— 這兩個突變共不共存、誰在前 —— 多層的局部系譜是其中的少數。 輸出格式、資料結構與驗收指標都要照這個比例設計,不要照「幾百棵樹」設計。

另外,本頁是規格,不是已驗證的方法(LiLT,工作名稱)。 文中所有數字都出自既有的原型實作,不是這份規格跑出來的結果 —— 自己實作時應當重新量一次,不要直接引用。

頻率譜為什麼看不到某些局部結構

「全域方法犧牲了局部解析度」如果只是印象,就沒有價值。 它其實是一個機制,而且可以換算成一個可直接與實際資料對照的數字: 要說服頻率譜「這裡是兩群不是一群」,需要幾個變異?

僅以 VAF 頻率譜重建 subclone 的四個步驟,及三項限制各自所在的參數 由左至右共四個步驟,為不使用任何長讀關係、僅以一維頻率譜重建的標準做法。 第一步為換算:單一位點的 read 分為帶變異與不帶變異兩組,其比值為 VAF, 再除以純度的一半換算為癌細胞比例。此步驟需要三個外部估計值, 分別為全樣本單一值的純度、每區段一個值的拷貝數,以及每個變異一個值的 multiplicity。 第二步為將一萬五千個換算後的數值匯為一維直方圖並擬合混合模型, 每一群具有位置、高度與寬度三個參數;每個變異的群歸屬為未觀測的潛在變數, 在該式中被邊際化。 第三步為決定群數,所依據者為模型選擇準則而非參數估計。 第四步為由各群的細胞比例建樹,其條件為父節點的細胞比例不得小於各子節點的比例總和; 此條件通常無法將候選收斂至唯一一棵,其餘由簡約性等偏好決定。 圖的下方列出三項限制,各以箭頭指回其所在的步驟: 群數不可辨識,位於潛在變數與群數; 寬度參數與群數不可區辨,位於寬度; 細胞比例依賴三個外部估計值,位於位置。 此三項限制與後續三個用途一一對應。 僅用 VAF:自 read 至樹的四個步驟 此為改善的對照基準。三項限制列於下方,各對應一個用途。 ① 換算 每個變異 → 一個細胞比例 a = 3 d = 8 VAF = a/d 再換算成 CCF: c = 2·VAF / ρ 須三個外部估計值 ρ 純度 · 全樣本一個 CN 拷貝數 · 每區段一個 μ multiplicity · 每變異 ② 擬合混合模型 15,500 個 c 匯為一維直方圖 CCF → a_m | z_m=k ~ BetaBin(d_m, θ_k, φ) θ_k 位置 · π_k 高度 · φ 寬度 z_m 於此處被邊際化 5% 的作用點即為此變數 ③ 決定 K 群數為何 BIC Dirichlet process MOBSTER 的冪次尾 此為模型選擇, 非參數估計 ④ 建樹 由 {c_k} 定先後 A B C c_A ≥ c_B + c_C 此條件通常無法 收斂至唯一一棵 其餘由簡約性等偏好 決定,非由資料決定 此路線的三項限制 —— 與 5% 的三個用途一一對應 限制一:K 不可辨識 兩群的 c 相等時,似然對 「一群或兩群」完全不變。 增加深度不作用於此方向。 參數:z → K → 由用途一補足 限制二:φ 與 K 不可區辨 單一寬峰與兩個鄰近的 窄峰擬合出近乎相同的圖。 其結果取決於懲罰項。 參數:φ → 由用途二補足 限制三:c 依賴換算 ρ 偏誤使整體平移; CN 或 μ 偏誤為逐變異 位移,將汙染群的組成。 參數:c → 由用途三補足
先看清楚被討論的是哪一系:這四個步驟就是那一系的全部路徑, 而觀測進入模型的唯一形式是第①步那個每變異的邊際計數與該處拷貝數。 PyClone、SciClone、PhyloWGS、Pairtree、DPClust 走的都是這條路徑(圖上只點名了 MOBSTER,因為第③步的冪次尾是它特有的),而它們正是這個領域被系統性評比的那一群。 圖右下角把三項限制各自掛回一個參數,本頁只用得到「限制一:兩群的 c 相等時, 似然對一群或兩群完全不變」—— 那正是下一張表最後一列。 (每項限制底下那句「由用途一/二/三補足」是前一篇補救這三項限制的做法, 本頁走的是另一條路,用不到它們。)

答案取決於兩群靠得多近,而且不是慢慢變多

兩群的期望 VAF頻率譜大約需要幾個變異才分得開
0.20 對 0.10約 4 個
0.20 對 0.15約 20 個
0.20 對 0.18約 120 個
0.20 對 0.20永遠不夠

最後一列不是「很難」,是沒有答案。 兩群的期望 VAF 相同時,把任何一個變異從這一群改指派到另一群, 直方圖完全不變 —— 兩群本來就落在同一根柱子裡。 加深度、加樣本、換更好的演算法都不會動它。 中篇的 K=2 例子即為此情形。

而一支只在一個連鎖視窗裡不同的 lineage,本來就帶不起那麼多變異

這是問題的另一半。一支細胞群若只在某個連鎖視窗內與它的母群體不同, 它與母群體的全部差別就是那個連鎖視窗裡的那幾個突變 —— 三個、五個。 它因此落在上表「開不出來」那一側,而且通常離得很遠。 所以這不是「不容易偵測」,是它在那個模型裡開不出來。

同一批 read,換一種看法

把那三個變異換成「它們在不在同一條分子上」,需要的證據量完全不同。 假設這三個突變其實共存,那麼在 20 條跨越它們的分子裡, 同時帶三個突變的分子預期只有錯誤地板那麼多 —— 1% × 20 條 ≈ 0.2 條

而實際看到 5 條。

這裡不需要門檻,也不需要任何統計檢定:預期 0.2、實際 5,結論就出來了。 而兩條路線看的是同一批 read —— 差別不在資料量,在把它們摘要成什麼。

上面那幾個數字是怎麼算出來的

「需要幾個變異」=多開一群要付的代價 ÷ 每個變異能提供的證據, 兩者都是對數似然,單位是 nat(自然對數底下的資訊量)。 本文不用 nat 當敘述單位,因為換算成「幾個變異」之後, 即可直接與實際資料比較。

  • 代價是模型複雜度罰則 12pKlogMeff: 參數兩個、有效觀測三千時約 8。
  • 每個變異的證據約為 d(θ1θ2)2/2θ(1θ), 其中 d 是深度、θ 是兩群的期望 VAF。 深度 50、0.20 對 0.15 時約 0.43,故 8÷0.4319,表中記為「約 20」。
  • θ1=θ2 時分子為零、商為無限大 —— 那就是兩群期望 VAF 相同的情形。 表中那四列是同一條曲線上的四個點,不是四件不同的事。
  • 共現那一側:5 條對 0.2 條的對數似然比約 12, 遠高於分辨兩個局部結構所需。這裡的錯誤地板(總量 δu,本頁取 1%) 來自後文「局部生成模型」的逐位點錯誤 ε 與家族誤標, 所以 δu算出來的而不是另一個自由參數,而且不是常數,要逐連鎖區段估。

一個必要的但書

上面整段的前提是「似然只吃邊際計數」,而有兩個既有方法滿足這個前提: PairClone 與 TreeClone 已經把「同一條短讀上一對鄰近 SNV 的共現」寫進了全域似然。 所以本規格與它們的差別不是「有沒有用共現」,而是這條通道有多寬 —— 通道寬度完全由讀長決定,而同一扇門,門後的東西差兩個數量級

數量級補充:短讀只有 0.15%–1.5% 的 sSNV 有伴,ONT 是 40.5%–98.1%

兩個 sSNV 要落在同一條讀序上才數得到,所以能不能用共現,取決於讀長與突變密度的比值。 300 bp 的短讀在 5–50/Mb 的突變負荷下只有 0.15%–1.5% 的 sSNV 有伴; ONT 的實測是 40.5%–98.1%(見「輸入與前處理契約」)。 這也是為什麼 PairClone 與 TreeClone 只處理成對位點:在那個產量下, k 大於 2 的視窗幾乎不存在,可變 k 沒有意義。

方法

方法概覽

方法先把互不重疊的 read 分成兩組:單變異連鎖視窗 S1 用於估計全域頻率基準, 多變異連鎖區段 S2 用於逐區段重建局部分子系譜。全域結果只作為局部推論的先驗, 局部結果不回饋全域,因此同一條 read 不會在兩個 likelihood 中重複計數。

整條路線的總覽:從原始 read 到兩份交付物 由上而下四個面板,分別畫出輸入整理、資料切分、全域與局部推論,以及交付物。 第一個面板是前處理,沿用既有工具: 左邊是四條來源未知的原始 read,以虛線外框表示; 中間是定相與體細胞標記之後的同一批 read,已依兩條單倍型分成兩組並帶上變異標記; 右邊是由它們整理出來的稀疏狀態表,畫成一格一格的方陣, 深色格代表帶有該位點的變異,淺色格代表參考型。 第二個面板是把位點切成兩堆:一條全基因體長條上散布著變異, 多數是孤立的,屬於 S1;少數幾處聚在一起並用虛線框起來,屬於 S2。 兩堆的 read 沒有交集,因此下方有兩支箭頭分別把 S1 送進全域基準、把 S2 送進局部推論。 第三個面板左邊是全域基準,畫成一張變異等位頻率的直方圖, 上面有幾個峰,峰頂各標一個圓點代表被定出來的 CCF 原子; 右邊是局部推論,畫出某一個連鎖區段的四條分子, 每條分子上有三個位點,實心代表帶有變異、空心代表參考型, 右方由這些分子推出一棵小小的局部系譜。 兩者之間有一支由左指向右的箭頭,標著單向:全域基準的結果只當作局部推論的先驗,局部不回饋全域。 最下方是兩份交付物: 左邊是逐連鎖區段目錄,畫成數列,每一列是一棵小樹加上幾個顏色不同的旗標; 右邊是給 somatic caller 的特徵,畫成一疊通道, 上面三條是本規格新增的三組特徵,下面三條是 caller 原本就有的。 方法概覽:三個步驟,兩份交付物 前處理 沿用既有工具 原始 read 定相 + 體細胞標記 稀疏狀態表 切成兩堆 read 互斥 S₁ · 單變異連鎖視窗(多數) S₂ · 虛線框起來的多變異連鎖區段 S₁ S₂ 全域基準 頻率譜分群 VAF 輸出:純度、CCF 原子 單向 局部推論 每個連鎖區段一次 該連鎖區段的分子 局部系譜 交付物 A · 逐連鎖區段目錄 每連鎖區段一列:小樹 + 比例 + 旗標 B · 給 caller 的特徵 局部系譜特徵 caller 既有
讀法:虛線外框的 read 是來源未知的原始資料, 上色之後才分得出兩條單倍型;實心圓是變異、空心圓是參考型。 中間那條長條是分水嶺 —— 孤立的變異進 S1,虛線框起來的成群變異進 S2兩堆的 read 沒有交集,所以全域基準算出的原子可作為局部推論的獨立先驗。 右側四條分子推出的小樹,就是方法的產出單元。
步驟資料產出在本方法中的地位
前處理配對 BAM + 候選集合定相、體細胞標記、連鎖區段與稀疏狀態表沿用既有工具;輸入契約必須固定
全域基準S1 的單變異資料純度、CCF 原子、原子不確定度與 overdispersion沿用既有頻率譜方法
局部推論S2 的單一連鎖區段 + 全域先驗狀態集合、比例、先後關係與品質旗標本規格的核心

估計目標:一個連鎖區段的三個量

每個連鎖區段只估計三個對象:局部狀態集合、各狀態的專屬分子比例, 以及狀態之間的先後關係。以下先定義估計單位與輸出,再引入似然與先驗。

一個連鎖區段的輸出:稀疏狀態表、唯一的相容樹,以及可以計數的旗標 三格由左到右。第一格是一個連鎖區段的稀疏狀態表: 全參考狀態 62 條、只帶第一個變異的 18 條、同時帶第一與第二個的 5 條、 同時帶第一與第三個的 7 條,而同時帶第二與第三個的是 0 條, 另外還有覆蓋遮罩不完整的分子另行計數。 這個連鎖區段的錯誤地板約為百分之一乘上分子總數,也就是不到一條, 所以 5 條遠高於地板,而 0 條還不足以排除那個狀態存在。 第二格是由這些狀態建出來的樹:全參考在最上,其下是只帶第一個變異的狀態, 再往下分成兩支,分別再加上第二個或第三個變異。 因為狀態集合通過四配子檢定,這棵樹是唯一的; 所以困難的地方不是搜尋樹,而是決定哪些狀態的比例高於地板。 第三格是這個連鎖區段輸出的那一列:狀態集合與各自比例、先後關係與後驗、 有效分子數與遮罩分布,以及四個旗標。 旗標是這條路線可以計數的產出,不是修辭。 一個連鎖區段輸出什麼 ① 稀疏狀態表 000 62 100 18 110 5 101 7 011 0 覆蓋遮罩不完整的另計 地板約為分子總數的百分之一 5 條遠高於地板 0 條還不足以排除 ② 相容則樹唯一 000 100 110 101 難的是選狀態,不是搜尋樹 四配子檢定通過,樹就唯一 ③ 輸出這一列 狀態集合與各自比例 先後關係,附後驗 有效分子數與遮罩分布 旗標 REFINES_GLOBAL FLOOR_LIMITED ORDER_UNRESOLVED INCOMPATIBLE 旗標是可以計數的產出
整節要交代的就是這張圖: 左邊是一個連鎖區段收到的東西(稀疏狀態表),中間是要估出來的東西(哪幾種狀態、誰在誰之前), 右邊是它最後輸出的那一列。中間那棵樹沒有經過搜尋 —— 狀態集合一旦選定,包含關係就把樹決定了,這是下面第三個小節的重點。

估的範圍:一個連鎖區段

沿用前一篇的連鎖區段u 是一個單倍型連鎖區段:一個連鎖視窗 × 一條 , 其位點集合為 Ju、位點數 ku=|Ju|

這樣切有一個具體的理由,不是為了方便: 一個連鎖區段裡的分子全部來自同一條染色體拷貝,另一條屬於另一個連鎖區段。 所以在連鎖區段內部,「兩個變異出現在同一條分子上」才真的等於 「它們在同一條拷貝上先後發生」。 若不先依家族切開,位於兩條拷貝上的一對變異(trans)永遠不會共存, 其計數表會與「兩支互不包含的 subclone」長得一模一樣 —— 連鎖區段這個切法就是為了讓這兩件事落進不同的連鎖區段而不會相遇。

要估的三個量

局部狀態h{0,1}ku, 即該連鎖區段那幾個位點上「帶或不帶變異」的樣式。 一個連鎖區段的估計目標就是下面三個量,沒有第四個

要估的量記號白話
狀態集合Xu{0,1}ku這個連鎖區段裡有哪幾種分子
專屬比例νu,h0,總和為 1每一種各佔多少
先後關係Xu 上的有根樹 Eu誰是誰的祖先

三者合起來記為 Lu=(Xu,Eu)νu, 樹根固定為全參考狀態 0ku

變數補充:ν 是專屬比例,不是累積比例

一條分子只由它所屬那一個節點的狀態產生,不由祖先節點產生; 祖先節點的 ν 因此是「還停留在該狀態的分子」的比例,可以很小。 換句話說,樹的形狀不對 ν 施加任何單調性, 它只規定哪些狀態可以同時存在

先後關係幾乎不用估:四配子相容性

先後關係看起來是三個量裡最難的,其實不是 —— 它不需要搜尋。 兩個位點稱為相容,若四種樣式 00,01,10,11 中至多出現三種。 只要 Xu每一對位點都相容,樹就已經被 Xu 決定了, 不必列舉任何樹形。實作上要做的只有兩件事:逐對檢查相容性, 以及在檢查失敗時把它當成品質訊號而不是新發現

k=2 的連鎖區段只會落入三類:串接、分岔、不相容 一個只有兩個位點的連鎖區段,它的分子只可能有四種狀態: 兩個都不帶、只帶第一個、只帶第二個、兩個都帶。 三欄各是一種情形,每一欄把這四種狀態各自的分子數列出來。 三欄的邊際刻意做成完全相同:總數都是八十條,位點 A 的變異分子都是三十九條, 位點 B 的都是十七條 —— 這是整張圖的重點,因為它表示光看每個位點各自的頻率 分不出這三種情形,要看四種狀態怎麼分配。 左欄是串接:沒有任何分子只帶 B,所以帶 B 的分子全部也帶 A, B 只可能長在 A 的後代裡。 中欄是分岔:沒有任何分子兩個都帶,所以 A 與 B 落在互不包含的兩支上。 右欄是不相容:四種狀態都有分子。在「每個變異只發生一次」的前提下這不可能, 所以它不是發現了一支新的 lineage,而是四配子檢定失敗 —— 代表這個連鎖區段的錯誤地板偏高或家族誤標偏多,是一個品質訊號。 讀法是:分類要用似然比檢定,不要對原始計數設門檻, 因為在偵測極限附近,兩種假設只差幾條分子。 同一組邊際,三種分子狀態分配 每一列是一種分子狀態:實心=帶變異、空心=參考型。三欄的邊際完全相同(見最下一列)。 串接 chain A → B 00 41 10 22 01 0 11 17 沒有分子「只帶 B」 → B 長在 A 的後代裡 分岔 fork A / B 00 24 10 39 01 17 11 0 沒有分子「兩個都帶」 → A、B 互不包含 不相容 incompatible 四種都有 00 33 10 30 01 8 11 9 四配子檢定失敗,非新 lineage N 80 · A 的變異 39 · B 的 17 N 80 · A 的變異 39 · B 的 17 N 80 · A 的變異 39 · B 的 17 右欄的出現率可以反過來估計單倍型分型的錯誤率 —— 模型因此自帶一個校正用的觀測量。 分類要用似然比檢定,不是對原始計數設門檻:在偵測極限附近,兩種假設只差幾條分子。
一個 k=2 的連鎖區段,其分子只有四種狀態。 三欄的總數與兩個位點各自的變異分子數完全相同(最下一列)—— 所以光看每個位點的頻率分不出這三種情形,要看四種狀態怎麼分配。 左欄沒有「只帶 B」的分子,故 B 只能長在 A 的後代裡; 中欄沒有「兩個都帶」的分子,故 A、B 互不包含; 右欄四種都有 —— 在「每個變異只發生一次」下不可能, 所以它不是發現了一支新 lineage,是四配子檢定失敗。
為什麼「兩兩相容」就足以決定一棵唯一的樹(含 k3

把每個位點 j 換成一個集合 CjXu帶有該變異的狀態所成的集合。 四種樣式少一種,說的正是 CACB 這兩個集合的關係:

缺哪一種樣式兩個集合的關係樹上的意思
沒有 01CBCAB 在 A 的後代裡
沒有 10CACBA 在 B 的後代裡
沒有 11CACB=A、B 在互不包含的兩支上
四種都有兩集合互相交錯樹上不存在這種關係

所以「相容」=「這兩個集合不交錯」,也就是要嘛一個包含另一個,要嘛不相交

唯一性:一族兩兩不交錯的集合,在數學上就是一棵樹 —— 包含關係本身給出父子關係(Cj 的父親是包含它的最小那一個), 不相交的兩個成為兄弟。沒有任何自由度留給搜尋,所以樹是 Xu函數,不是待估的參數。 剩下的自由度只有一種:一條邊上若有數個位點的集合完全相同, 它們同時發生、在樹上擠在同一條邊,彼此的先後資料沒有講 —— 這就是「本質上唯一」那個「本質上」的全部內容。

k3 為什麼不必額外檢查:關鍵在於兩兩相容不只是必要條件, 它同時是充分條件(perfect phylogeny 定理,Gusfield 1991)。 也就是說不存在「每一對都相容、但三個一起就不相容」的反例 —— 所以三個位點不必檢查 3 組三元組,k 個位點也只要檢查 k(k1)/2 對,成本是 k2 而不是 2k。 這正是「先後關係幾乎不用估」的來源。

實作上因此只需維護一張逐對的相容性表; 它同時就是右欄那個品質訊號的計數來源。

這給出兩個貫穿後面各節的結論:

  • 困難不在搜尋樹,在決定 Xu —— 哪些狀態的比例真的高於錯誤地板。後續的生成模型、先驗與候選選擇都在處理這一件事。
  • 四格全滿不是發現了一支新 lineage,是四配子檢定失敗, 代表這個連鎖區段的錯誤地板偏高或家族誤標偏多。相容性檢定因此同時是品質訊號。
解讀補充:未觀測狀態可能是現存或歷史狀態

Xu 可以含沒有被任何分子觀測到的狀態。 觀測到 110101100 沒看到時, 兩者若不共用一個帶第一個突變的祖先,第一個突變就得發生兩次。 這種狀態只有兩種身分 —— 現存但沒抽到(ν>0), 或已被後代取代的歷史狀態(ν=0)。 前一篇還有第三種:純粹因為「一步一個突變」的表示法而被迫補進來的記帳節點; 本規格允許一條邊帶多個突變,那一種就隨表示法一起消失了。 那條警告仍然適用,但適用範圍窄了一半。

什麼時候算「估完了」

三個核心估計量之外,細胞比例與全域群對應還需要外部資訊。 每項輸出各有成立條件;條件不足時輸出旗標,而不是猜一個值 —— 這條規則是後文輸出旗標的語義來源。

要回答的問題估到什麼需要什麼條件不足時的輸出
有哪幾種、各佔多少Xuνu只需要該連鎖區段的 read 與就地估到的錯誤率FLOOR_LIMITED
誰在誰之前Xu 上的先後關係兩個單一狀態皆可見,且雙狀態的有無高於地板ORDER_UNRESOLVED
換算成細胞比例各節點的 純度與該處 ,由全域基準提供只報分子比例,標 CELL_SCALE_UNAVAILABLE
對得上哪一個全域群局部節點與全域群的對應節點內至少一個變異有可信的全域指派UNLINKED_TO_GLOBAL
局部節點是一類分子狀態,不是一群細胞

上表的「換算成細胞比例」只是尺度轉換,它不會把節點變成細胞群。 一個節點只表示「這些分子在 Ju 上帶有同一組突變」。 兩群基因體其他地方不同、但在這個連鎖視窗內完全相同的細胞,會落進同一個節點。 本實驗室的甲基化觀測正是這件事的實例:單一個「單倍型 + 等位」狀態之下, 仍可再分出數群甲基化模式。因此節點數的正確讀法恆為「至少此數」。 反向也不成立:同一群細胞在不同連鎖區段裡對應到不同的局部狀態, 因為每個連鎖區段看的位點集合不同。

局部生成模型

這一節只回答一個問題:給定某個連鎖區段裡各種真實分子狀態的比例 νu,資料中的一條 read 為什麼會呈現目前看到的 0/1 字母? 由真實分子池到觀測資料共有四步:

  1. 抽出一種真實狀態。一條分子先依 νu,屬於某個完整狀態 Hr=h
  2. 套用這條 read 的覆蓋範圍。read 通常只蓋住位點集合 Mr,未覆蓋的位置保持未知。
  3. 通過觀測錯誤。逐位點錯誤可能改變一個字母,家族誤標則可能把整條分子分到錯的單倍型家族。
  4. 形成實際資料。通過前述步驟後的字母 Yr,才是狀態表裡真正計數的觀測。
局部生成鏈:由這個連鎖區段的狀態比例,長到一條 read 上看到的字母 由左到右四格,說明一條讀序上的字母是怎麼產生的。 第一格是這個連鎖區段的狀態比例:四種局部狀態各佔多少分子,畫成四條長短不同的長條。 這組比例是自由參數,先驗來自全域基準。 第二格是從這組比例抽出一條分子,它有一個真實狀態,圖中畫的是三個位點皆帶變異的 111; 實心圓代表帶有變異,空心圓代表參考型。 第三格是覆蓋遮罩:第二個位點沒有可靠觀測,以虛線空心圓表示未知; 第一與第三個位點讀到 1 和 1。沒觀測到的位點必須邊際化,不可以填成參考型。 第四格是兩個錯誤通道:逐位點錯誤一次只改一個字母, 整條分子的家族誤標則一次把整條分子搬到另一個家族; 通過這兩道之後看到的字母才是真正的資料。 最下方是與前一篇的唯一結構差異,而且只差在第一格。 左邊是前一篇:全域的 clone tree 以一個確定函數決定局部比例,局部沒有任何自由度。 右邊是本規格:同一棵全域樹只以虛線箭頭提供先驗,局部比例是自由參數。 這就是「提高局部解析度」在技術上的全部內容。 一條 read 上的字母是怎麼來的 ① 這個連鎖區段的狀態比例 000 101 111 自由參數 先驗來自全域基準 ② 抽一條分子 它的真實狀態 111 ③ 覆蓋遮罩 第二個位點沒有可靠觀測 要邊際化,不可填成 0 ④ 兩個錯誤通道 逐位點:改一個字母 誤標:整條搬走 過完這兩道之後 看到的字母才是資料 與前一篇的唯一結構差異,就只差在第 ① 格 前一篇 確定函數 ν 沒有自由度 本規格 只當先驗 ν 是自由參數
上半部依序畫出完整生成過程: 先從此連鎖區段的狀態比例 νu 抽出一條具有完整真實狀態的分子, 再只保留該 read 實際覆蓋的位置,最後通過逐位點錯誤與家族誤標,得到資料中觀測到的字母。 下半部才比較與前一篇的差異:前一篇由全域 clone tree 確定 νu; 本規格讓 νu 成為局部自由參數,全域結果只提供先驗。

先看一條只有部分覆蓋的 read

設三位點連鎖區段的狀態集合為 Xu={000,101,111},三種狀態的專屬比例依序為 νu=(0.60,0.30,0.10)。現在有一條 read 只覆蓋第一與第三個位點, 並在這兩處讀到 1–1。在尚未加入錯誤時,這條 read 同時相容於 101111,因此其觀測機率為 0.30+0.10=0.40一條 read 支持的是一組狀態,不是一個狀態 —— 這句話就是下面整節的全部內容。

符號表:HrMrYr 三者的分工(隨時可回來查)
記號在此例中的值意義
Hr000101111這條分子的完整真實狀態,未直接觀測
Mr{1,3}這條 read 實際覆蓋的位點
Hr,Mr0–01–1把完整狀態投影到覆蓋位置後,原本應看見的字母
Yr例如 1–1經過錯誤通道後,資料中真正記錄的字母

三者的關係是一條單向鏈:Hr 經覆蓋遮罩投影成 Hr,Mr, 再經錯誤通道變成資料裡的 Yr。推論要走的是反方向。

先條件於 read 已進入正確的單倍型家族;此時 read-level 的骨架可寫成:

P(Hr=hνu)=νu,h,P(Yr=yMr=M,νu,eu)=hXuνu,hAu,M(yh,eu)

Au,M(yh,eu) 是單一家族內的觀測通道:給定完整真實狀態 h, 它處理部分覆蓋與逐位點錯誤,並給出最後看到 y 的機率。 eu 為本連鎖區段的逐位點錯誤參數。跨家族誤標會連結另一個連鎖區段,於第二個細節另行加入。

這條式子的讀法是:逐一考慮所有可能的真實狀態, 以其分子比例 νu,h 加權,再乘上它經觀測通道變成 y 的機率。 以下各小節先拆開 A,再把跨家族誤標接回來。

第一個細節:未覆蓋位點要邊際化

若暫時忽略所有錯誤,觀測通道只剩下「h 投影至 M 後是否等於 y」這項檢查。上面的完整式子便化為:

qu,M(y)=h:hM=yνu,h

MJu 是該 read 的覆蓋遮罩,y 是它讀到的字母。 沒覆蓋到的位點不可以填成參考型 —— 那會把一條中立的 read 變成反對某個狀態的證據。

一個真實連鎖視窗的 read 覆蓋:八成的分子只看得到一個位點 上半是一個真實分析單位裡 read 覆蓋情況的統計。 這個連鎖區段有三個 somatic 位點、一百九十四條 read,連鎖視窗寬三十九點四 kb。 三條長條由上而下分別是:同時蓋到三個位點的 read 只有三條,佔百分之一點五; 蓋到兩個位點的三十一條,佔百分之十六;只蓋到一個位點的一百六十條,佔百分之八十二點五。 下半把同一批 read 按位點對拆開,看每一對到底有幾條 read 同時蓋到。 第一對與第二對各有十七條與二十條,分別判出分岔與串接; 第三對只有三條,而且三條都是參考型,什麼也判不出來。 最下方是這張圖的結論:這個連鎖區段的演化樹靠的是十七與二十這兩個數字,不是一百九十四; 而第三對實質上缺乏資料,仍被演算法判定了關係。 讀法是:有效樣本數是「每一個位點對各自的跨越深度」,不是這個連鎖區段的總 read 數。 真實資料的覆蓋結構 HCC1395_HKU · chr12:981,725–1,021,146 · HP1 · k=3 · 194 reads · 連鎖視窗 39.4 kb 一條 read 蓋得到幾個位點 蓋到 3 個 3 條 1.5% 蓋到 2 個 31 條 16.0% 只蓋到 1 個 160 條 82.5% 八成的分子只看得到一個位點 —— 它對「兩個位點的共現」完全沒有意見。 把同一批 read 按位點對拆開 A × C n = 3 三條都是參考型 什麼也判不出來 A × B n = 17 RA 12 · AR 5 · 沒有 AA 判為分岔 B × C n = 20 AA 6 · AR 9 · RR 5 判為串接 這個連鎖區段的樹靠的是 17 與 20,不是 194。而 A × C 那一對沒有資料 卻一樣被演算法決定了關係 —— 「沒看到」被當成了「不存在」。
約八成的 read 僅覆蓋三個位點中的一個。 將同一批 read 依位點對拆分後可見:此連鎖區段的判定依據為 17 與 20 兩個數字,而非 194; 其中一對僅有 3 條 read,仍產出了一條關係判定。
證據補充:部分覆蓋為何只能算一個約束

一條有 x 個未覆蓋位點的 read,與 2x 個狀態相容, 那些狀態構成狀態空間的一個子立方體。它主張的是「真實狀態落在其中」, 不是 2x 個觀測 —— 把兩者混為一談會把證據量灌大 2x 倍。 在似然寫法下這是自動的(上式的加總就是那個子立方體),但報告層不自動: 既有實作有一個 k=3 的連鎖視窗,靠一個全跨樣式(3 條 read) 加上十一個部分樣式就「解出唯一解」。 這種連鎖區段的唯一性大部分來自簡約性與先驗,不是來自資料。 因此每個連鎖區段必須分開報全跨分子部分分子各貢獻多少證據; 只有全跨樣式與根是跨候選不變的,被部分 read 見證的狀態不是。 順帶一個反直覺的推論:部分覆蓋越多,相容的候選越多 —— 覆蓋越差的連鎖區段越容易看起來「並列」,而不是越容易被排除。

整個連鎖區段共用同一個分子池

整個連鎖區段共用一組 νu,每條 read 由它自己的遮罩投影出來; 不同遮罩不各配一個獨立的比例向量,因為它們取樣的是同一池分子。 在 PCR-free 長讀下,一條 read 就是一條分子, 所以連鎖區段內部的計數預設為 multinomial,額外離散度沒有來源; 只有殘差真的顯示過度離散時才加,而且要說得出它來自哪裡。

把這個區段的所有 read 記為 Du={(Mr,Yr)}r=1nu。暫時忽略跨家族誤標時,單一連鎖區段的基礎似然只是逐 read 相乘,而每一項都是前述的加權和:

P(Duνu,eu)=r=1nu[hXuνu,hAu,Mr(Yrh,eu)]

不同 read 可以有不同的遮罩 Mr,但全部共用同一組 νu。 同一遮罩的 read 可彙總成 multinomial 計數;兩種寫法是同一個模型。

第二個細節:兩種錯誤通道的形狀不同

前式中的 A 不能只寫成一個無來源的錯誤率。 逐位點錯誤與整條分子的家族誤標會產生不同形狀的假狀態,因此必須分開校準。 前者留在 A 內;後者要把成對的兩個連鎖區段一起寫:

逐位點錯誤與整條分子誤標的形狀不同,而且位點越多,誤標越是主角 左邊是逐位點的定序錯誤。它一次只改一個字母, 每個位點各自獨立,所以要一次改對三個位點才會生出一個完全不同的狀態, 機率是單點錯誤率的三次方,小到可以忽略。 右邊是整條分子的家族誤標。一條本來屬於另一個單倍型家族的分子被標錯, 整條就這樣跳到這個連鎖區段裡,它的三個位點是一起過來的, 所以機率只有一個誤標率,跟位點數完全無關。 中間下方是兩者的比較表。以單點錯誤率千分之五為例, 差一個位點時兩者相當;差兩個位點時誤標大約高一百倍; 差三個位點時高兩萬倍以上。 所以位點越多,能造出假狀態的幾乎只剩誤標這一條路。 最下方是後果:如果模型只寫了逐位點錯誤, 那麼在全參考的背景上冒出來的三位點全變異狀態, 在模型眼中會是「錯誤不可能造出來的」,於是只剩一個解釋 —— 一個新的 clone。 而它真正的來源是隔壁那個家族。 讀法是:一個參數描述不了兩種形狀不同的錯誤,寫漏的那一種會直接變成假的 subclone。 兩種錯誤的形狀不同:一種改一個字母,一種搬走一整條分子 逐位點定序錯誤 每個位點獨立,一次改一個字母 0 0 0 只有中間這個位點被讀錯 0 1 0 要生出一個差三個位點的狀態, 同時錯三次:ε³ 整條分子的家族誤標 三個位點一起過來 另一個單倍型家族 1 1 1 整條被標到這個連鎖區段 1 1 1 機率只有一個誤標率 η, 與位點數無關 差幾個位點,誰才是主角(以 ε = 0.005、η = 0.01、來源比例 0.3 為例) 差 1 個位點 逐位點 5×10⁻³ 誤標 3×10⁻³ 兩者相當 差 2 個位點 逐位點 2.5×10⁻⁵ 誤標 3×10⁻³ 誤標高約 100 倍 差 3 個位點 逐位點 1.3×10⁻⁷ 誤標 3×10⁻³ 誤標高約 20,000 倍 只寫逐位點錯誤會發生什麼 在全 0 的背景上冒出 111:模型認為錯誤造不出來(10⁻⁷), 於是只剩一個解釋 —— 一個新的 clone。而它真正的來源是隔壁那個家族。 讀法:一個參數描述不了兩種形狀不同的錯誤;寫漏的那一種會直接變成假的 subclone。
左右兩種錯誤的形狀不同:一種一次改一個字母,一種一次搬走整條分子。 位點差得越多,逐位點錯誤的機率掉得越快,而誤標完全不受位點數影響 —— 所以差三個位點時,能造出假狀態的幾乎只剩誤標這一條路。

u1=u 為目前的單倍型家族,u2 為同一連鎖視窗的另一個家族, αvu 為「一條最後被標成 u 的 read,實際來自家族 v」的校準權重。 則真正用於 unit u 的觀測機率為:

pu,M(y)=v{u1,u2}αvuhXvνv,hAv,M(yh,ev),vαvu=1

αu1u 是正確標記的來源,αu2u 是由另一家族誤標進來的來源。 這些權重由成對家族的分子數與誤標率導出,不是每個 unit 自由擬合的比例。 因此 production likelihood 以 pu,Mr(Yr) 取代前一個基礎似然括號內的值。

νu,hobs=(1δu)·νu,h+δu·ru,h

這是前述成對觀測通道的摘要寫法。ru,h 由逐位點錯誤與整條分子的家族誤標算出,不是自由參數。 δu地板總量,由逐位點錯誤 ε 與家族誤標算出,不是另一個自由參數; δuru,h 就是這個連鎖區段能分辨的最小比例,也是所有偵測下界的來源。

一格是空的,有三種完全不同的原因 三欄並列,三欄的觀測完全一樣:那一格的計數是零。但成因不同,處理方式也不同。 第一欄是抽樣零。那個狀態真的存在,只是跨越這幾個位點的分子太少,沒抽到。 判準是機率:沒抽到的機率等於一減去它的比例,再取跨越深度次方。 跨越深度三、比例三成時,這個機率大約是三分之一, 所以「沒看到」幾乎沒有排除任何東西。跨越深度五十九時才降到百分之五。 處理方式是讓似然自己算,不需要任何門檻。 第二欄是結構零。沒有任何 clone 帶有那個基因型,所以它的真實比例就是零。 這是我們真正想推論的結論,不是輸入。 第三欄是錯誤地板。那個基因型確實不存在,但定序錯誤與家族誤標仍然會產生看起來像它的 read。 所以觀測到的比例永遠不會真的是零,而是壓在一個地板上。 中間橫跨一條警告:Dirichlet 的參數不允許是零, 所以結構零必須經過錯誤通道抬到地板之上,模型才寫得下去。 最下方是三者的處理方式對照。 讀法是:把這三種零混為一談,就會把錯誤產生的少數幾條 read 讀成一個新的 clone。 同樣是「這一格沒有 read」,成因有三種 ① 抽樣零 狀態存在,只是沒抽到 P(看不到) = (1 − q)ⁿ n = 3、q = 0.30 → 0.34 三分之一的機率誤判成「不存在」 n = 59 才降到 0.05 處理:讓似然自己算 ② 結構零 沒有任何 clone 帶這個基因型 q = 0 這是要推論出來的結論 不是輸入,也不是門檻判出來的 硬性排除等於把結論當成前提 處理:讓資料把 q 壓下去 ③ 錯誤地板 基因型不存在,但錯誤會造出它 q觀測 = (1−ε)·q + ε·r 觀測比例有一個下界 逐位點錯誤 + 家族誤標都在這裡 地板高度決定得出來的解析度下限 處理:寫進觀測通道 為什麼 ② 與 ③ 一定要分開寫:Dirichlet 的參數不可以是零 結構零給出 q = 0,但 Dirichlet 分量必須為正。所以結構零一定要 先經過 ③ 的通道抬到地板之上,式子才寫得下去 —— 這不是技術細節, 它就是「零不代表不可能」的數學形式。 三者的觀測一模一樣,都是「這一格計數為 0」 分得開它們的不是那個 0,是跨越深度、其他連鎖區段的證據,以及就地估到的錯誤率。 混為一談的後果:錯誤造出來的少數幾條 read,會被讀成一個新的 clone。 讀法:門檻只能把三種零一起丟掉;似然可以分別給它們不同的重量。
三欄的觀測一模一樣,都是「這一格計數為 0」。 分得開它們的不是那個 0,而是跨越深度、其他連鎖區段的證據,與就地估到的錯誤率。 第③欄那個地板就是上面那條式子的 δuru,h結構零在觀測上永遠到不了零, 所以「這個狀態不存在」這件事只能由它與地板的距離來講,不能由計數是不是 0 來講。

家族誤標尤其要成對估計,因為它連結同一連鎖視窗的兩個連鎖區段。 在 k2 時它是假狀態的主要來源:只寫逐位點錯誤的話, 全 0 背景上出現的 111 在模型眼中「錯誤造不出來」(約 ε3), 於是只剩一個解釋 —— 一支新的 lineage。而它真正的來源是隔壁那個家族。

觀測要經過幾道閘門:每個位點一道,最後兩個家族之間再交換一次 兩條並排的通道分別是同一個連鎖視窗的兩個單倍型家族。 每條通道由左到右經過相同的幾個階段。 起點是真實的分子組成。 接著是每個位點各一道錯誤閘門,閘門會把參考型讀成變異型,也會把變異型讀成參考型, 兩個方向的機率不一樣,所以閘門不是對稱的。 k 個位點就串 k 道閘門,把它們串起來就是公式裡那個張量積。 最後一道是家族交換閘門:大部分的分子留在自己的通道, 但有一小部分會被標到另一條通道去,圖上用兩條交叉的虛線表示。 終點是實際觀測到的計數。 圖的下方指出這件事的後果:因為交換閘門把兩條通道接在一起, 兩個家族必須一起估。分開估會讓譜系內比例每個連鎖區段都往零點五偏一點, 而且幾萬個連鎖區段偏的方向相同,那不是把譜變寬,是把整條譜平移。 讀法是:公式裡的三個符號各對應圖上一道閘門。 觀測要經過幾道閘門 家族 H1 家族 H2 真實組成 真實組成 位點 1 閘門 ε₁ , η₁ 位點 1 閘門 ε₁ , η₁ 位點 2 閘門 ε₂ , η₂ 位點 2 閘門 ε₂ , η₂ 這兩道串起來就是公式裡的 T = E₁ ⊗ E₂ (k 個位點就串 k 道) 家族交換閘門 δ 1 − δ δ 觀測計數 觀測計數 ε 是「把 REF 讀成 ALT」,η 是「把 ALT 讀成 REF」—— 兩個方向機率不同,所以閘門不對稱。 兩者都可以用同一個連鎖視窗裡的 germline 雜合位點就地估出來,不必外部參數。 交換閘門把兩條通道接在一起,所以兩個家族必須一起估。分開估會讓 ϱ 每個連鎖區段 都往 0.5 偏一點,而且幾萬個連鎖區段偏的方向相同 —— 那不是把譜變寬,是把整條譜平移。
每個位點一道錯誤閘門;最末一道為家族交換閘門, 少部分分子被標至另一條通道。交換閘門連結兩條通道,故兩個家族必須一併估計 —— 分開估的偏誤方向在幾萬個連鎖區段上相同,那不是把估計變寬,是整批平移 (圖下方的 ϱ 是連鎖區段內的分子比例,在本頁就是 νu)。 「輸入與前處理契約」中三種 HP3 處理之所以必須明選,也是這個閘門的問題: 它決定哪些分子進得了這兩條通道、哪些根本不進來。
實作風險:少數錯字為何會改變整個候選集合

既有實作把狀態表當成精確值:沒有逐 read 錯誤模型, 沒有接受一個位點為「已覆蓋」的最低品質門檻,也沒有「一個樣式要幾條 read 才收」的下限。 在頻率路線上,一個誤讀只是把某個 VAF 推偏一點點; 在這條路線上不是 —— 一個被數條 read 共享的誤讀會引入一個樣本從來沒有的狀態, 而那個狀態會改變整個候選集合,連帶改變樹的形狀、旗標與計數。 這就是為什麼本節的錯誤地板不是數值穩定用的小常數, 而是這條路線唯一擋得住這件事的東西;也是為什麼前處理契約的三個門檻 (覆蓋一個位點的最低品質、算連鎖的最低 read 數、保留一個連鎖視窗的最低跨越深度) 必須寫死 —— 它們決定哪些字母有資格進到狀態表裡。

第三個細節:只准 0 → 1 時,LOH 會被誤讀

四配子相容性與 perfect phylogeny 都預設突變只增不減,所以 一段發生 而失去某個突變的區域,在模型眼中會被描述成 「那個狀態從來沒有取得過」,而不是「取得後又失去」。

來源補充:既有實作的 Camin–Sokal 條件在哪裡漏掉了「喪失」

既有實作的 Camin–Sokal 條件同樣只准 01:復發(同一個突變在兩支各發生一次) 不被罰,但喪失完全不在模型的表達範圍內 —— 它不是被罰得很重,是連寫都寫不出來。 既有實作在整條鏈上沒有任何一處修正這件事,所以 LOH 區段的局部系譜會安靜地錯, 而且錯得跟一個正常結果長得一樣。

本規格的處理是把它變成明示的邊界,而不是默默承受: 拷貝數由全域基準提供,故每個連鎖區段都知道自己落在哪一種區段。 非中性區段的連鎖區段一律標 LOSS_UNMODELLED, 其局部系譜只作報告、不進任何樣本層統計量; 第一版並建議只納入 CN-neutral 的雜合區段。 這會減少可用連鎖區段(在 CN 變異廣泛的實體腫瘤可能減少一半以上), 所以納入與排除的連鎖區段數必須逐階段載明 —— 這正好也是連鎖率之外的第二個分層變數。

家族標號的方向不影響局部系譜

家族標號在每個 內獨立決定, 故仍引入 σu{,調} 並邊際化。但與前一篇不同的是, 這裡有一件可以放心的事:局部系譜的形狀完全不依賴 σu —— 連鎖區段內部的狀態集合與包含關係與家族標號無關,所以 σu 不進核心推論,只進註解層。

作用範圍補充:那 σu 到底在哪兩處還是要緊

兩處都在「把局部結果接到外面」的時候: 把節點對應到某一條 germline 單倍型(用於後文 caller 的單倍型特徵), 以及把節點對應到全域群。兩者都屬於註解層 —— 它們錯了,這個連鎖區段的狀態集合、比例與先後關係一個都不會變

與前一篇的唯一結構差異

項目前一篇(全域目標)本規格(局部目標)
局部比例 νu全域 (w,T,z)確定函數,沒有自由度自由參數,帶來自全域基準的先驗
換得的東西每個連鎖區段都對全域參數施加限制連鎖區段內部可以出現全域表達不了的結構
付出的代價局部只能表達全域已經有的群連鎖區段之間不再互相限制,後驗較寬
推論結構全域取樣的內層包住每個連鎖區段連鎖區段彼此獨立,可完全平行

這張表就是整份規格的決定點。 νu 自由,正是「提高局部解析度」的技術內容; 它同時也是「不再限制全域參數」的原因。兩件事是同一個改動的兩面,不能只要一面。

全域頻率基準如何成為獨立先驗

局部比例自由之後,20 條分子撐不起一個穩定的估計。 全域頻率譜在這裡仍然有用 —— 但只能單向進來。 可以這樣做的理由,正是前一篇為了讓似然相乘而建立的那個性質: 單變異連鎖視窗 S1 與多變異連鎖區段 S2 的 read 沒有交集

單向先驗:全域頻率譜只用單變異連鎖視窗擬合 最上方是一條全基因體的長條,上面散布著幾個較深的色塊。 長條的大部分是只含一個變異的連鎖視窗,稱為 S1,佔全部變異的九成; 深色塊是含兩個以上變異的連鎖區段,稱為 S2,佔約一成。 兩堆的 read 沒有交集,這正是前一篇用來把似然寫成兩個因子相乘的同一個性質。 中間是兩個方塊。左邊是全域基準:只用 S1 擬合全域頻率譜, 輸入是 S1 的 alt 計數、深度與拷貝數,輸出是純度、CCF 原子與權重、以及 overdispersion。 全域基準可沿用任一既有實作,不是本規格的貢獻, 唯一的硬性要求是 S2 的位點必須整批排除。 右邊是局部推論:每個 S2 連鎖區段各推論一次局部系譜, 輸入是該連鎖區段的稀疏狀態表,先驗來自全域基準, 輸出是局部狀態集合、比例、先後關係與旗標,連鎖區段之間彼此獨立、可以完全平行。 兩個方塊之間只有一個單向箭頭:局部不回饋全域,這是刻意切斷的回饋。 最下方是警告:如果全域基準用到了 S2 的 read,同一批 read 會先算進頻率項再算進局部項, 那正是前一篇用兩個因子相乘避免掉的重複計數。 另外拷貝數與純度用的是全基因體深度,其中含 S2 的 read, 這個界線有限但必須量出來並報告。 全域訊號怎麼用,才不會把同一批 read 算兩次 S₁:只有一個變異的連鎖視窗 S₂:兩個以上變異的連鎖區段(深色塊) 兩堆的 read 沒有交集 —— 這正是前一篇用來把似然寫成兩個因子相乘的同一個性質。 全域基準 · 只用 S₁ 輸入:S₁ 的 alt/depth/CN 輸出:純度、CCF 原子與權重 用哪一個既有實作都可以 沿用既有方法,不是本規格的貢獻 唯一的硬性要求 S₂ 的位點必須整批排除 單向 局部不回饋全域 (cut posterior) 局部推論 · 每個連鎖區段一次 輸入:該連鎖區段的稀疏狀態表 先驗:全域基準的輸出 輸出:局部狀態集合與比例    先後關係與旗標 連鎖區段之間彼此獨立 可以完全平行,沒有全域取樣 如果全域基準用到了 S₂ 的 read,這條路就白做了 那等於把同一批 read 先算進頻率項,再算進局部項 —— 這正是前一篇指出、 並且用「兩個因子相乘」避免掉的重複計數。 另一個必須量出來的邊界:CN 與純度的估計用的是全基因體深度,其中含 S₂ 的 read。 一個 20 kb 連鎖視窗通常佔其所屬區段不到 1%,這個界線有限,但必須量出來並報告。
上方是同一個切法:兩堆的 read 互不重疊。 因此全域基準只用 S1 擬合,其輸出對局部推論而言是一份獨立的先驗。 中間只有一個單向箭頭 —— 局部不回饋全域,是刻意切斷的。

全域頻率基準只使用 S1

全域基準是一個標準的 分群, 用哪一個既有實作都可以,它不是本規格的貢獻。它必須輸出四樣東西: 純度 ρ^;CCF 原子的位置 c^j 與權重 ω^j; 各原子的後驗不確定度;以及分層估到的 overdispersion φ^(d,CN)

唯一的硬性要求是:S2 的位點必須整批排除。 它是這條路線不重複計數的全部依據,而且可以寫成檢查程式: 全域基準的輸入位點集合與 S2 的交集必須為空。

這個要求的代價比紙上估計大得多,而且逐樣本差很多。 中篇以 Poisson 模型推得「約 9.5% 的變異落在多變異連鎖區段」,若真是如此, 排掉 S2 對頻率譜幾乎沒有影響。但既有實作實測到的是 40.5% 到 98.1%(見「輸入與前處理契約」)—— 在連鎖率高的樣本上, 排掉 S2 等於排掉大半個變異集合。 所以S1 剩幾個變異,是逐樣本必須先算、必須寫進報告的量, 不是一個可以忽略的尾數。

後果補充:S1 被排小之後,該傳什麼、什麼時候該直接棄權

S1 太小時全域基準的原子後驗會變寬,那份不確定度必須在局部推論中逐次傳遞, 不能只傳原子的位置 —— 只傳位置等於把一個很寬的估計當成確定值使用。

而在極端情況下 S1 可能小到撐不起分群。 那時全域基準就該直接棄權,全部連鎖區段標 UNLINKED_TO_GLOBAL, 而不是拿一個估不準的原子集合去當先驗 —— 後者會把全域的雜訊當成局部的先驗知識, 而且在輸出裡看不出來。

兩個邊界必須量出來,不能用講的:深度不是嚴格互斥,先驗不是硬約束

其一,深度不是嚴格互斥的。拷貝數分段與純度用的是全基因體深度, 其中含有 S2 連鎖視窗的 read。一個 20 kb 尺度的連鎖視窗通常佔其所屬 CN 區段不到 1%, 所以這個洩漏有界 —— 但界線要實際算出來寫進報告, 嚴格版本則在 CN 分段時把 S2 連鎖視窗一併排除,並比較兩者的差異。 其二,先驗不是硬約束。全域基準的原子若當成 ν 只能取的值, 局部解析度就歸零了;下面那個逃逸質量存在的唯一理由,就是不讓這件事發生。

先驗:原子加上逃逸質量

全域頻率譜當作局部比例的先驗:原子加上逃逸質量 三個並排的分布圖,橫軸都是某一個局部狀態的分子比例。 左圖是稀疏局部計數自己的似然:只有 20 條跨越分子時,這條曲線又寬又平, 單靠這個連鎖區段什麼都定不下來。 中圖是全域基準的全域頻率譜提供的先驗:幾根很窄的尖峰落在全域算出來的 CCF 位置上, 底下另有一條很低但不為零的水平線,代表逃逸質量, 它保留了「這一支 lineage 只在局部存在、全域頻譜上沒有對應的峰」這個可能。 右圖是兩者合起來的後驗:主峰被拉到某一個原子上而變窄, 同時仍留下一個較小的次峰,那是只有局部看得到的那一支。 最下方是警告:逃逸質量是這條路線唯一一個靠判斷的旋鈕。 設成零,局部結構會被強行貼回全域原子,解析度歸零,等於沒做局部推論; 設成一,等於丟掉全域資訊,20 條分子撐不起任何結論。 它必須預先登錄並做敏感度分析。 全域頻率譜怎麼變成局部比例的先驗 稀疏的局部計數 20 條分子,比例定不下來 似然很寬 單靠這 20 條,什麼都定不下來 全域頻譜給的先驗 原子 + 逃逸質量 尖峰的位置來自全域基準 那條低線是逃逸質量 合起來的後驗 收緊,但沒有把逃逸封死 主峰貼上其中一個原子 次峰是局部才看得到的那一支 逃逸質量是這條路線唯一一個靠判斷的旋鈕 設成 0:局部結構被強行貼回全域原子,解析度歸零 —— 等於沒做局部推論。 設成 1:等於丟掉全域資訊,20 條分子撐不起結論。它必須預先登錄並做敏感度分析。
左:20 條分子自己的似然又寬又平。 中:全域基準給的先驗是幾根窄峰加上一條不為零的底線。 右:後驗被拉緊,但那條底線保留了「這一支只在局部存在」的可能。

先在細胞尺度上寫先驗,再換算成分子比例。設節點 h 的細胞比例為 cu,h

p(cu,h)=(1λ)jω^j·Beta(cu,h;aj,bj)+λ·Beta(cu,h;a0,b0)

(aj,bj) 由全域基準第 j 個原子的後驗均值與變異決定; (a0,b0) 是弱的擴散成分。λ逃逸質量: 先驗上認為這一支不對應任何全域群的機率。 λ0 把局部結構強行貼回全域原子(解析度歸零); λ1 等於丟掉全域資訊。λ 必須預先登錄並做敏感度分析。

由細胞比例換算為連鎖區段內的分子比例,用的是前一篇那條把拷貝數放進分母的式子: 一個 clone 在該處的每一條拷貝各貢獻一份,分母是該連鎖區段全部的分子數。 正常細胞是其中一項,其局部狀態恆為全參考 —— 所以純度不是額外的換算因子。 若該處的 CN 或 multiplicity 未定,這一步就停在分子尺度,輸出標 CELL_SCALE_UNAVAILABLE,而不是套一個預設的二倍體值。

先後關係的先驗:只能軟用

logp(Lu)=(h,h')Eulogωanc(h,h')+κsplit·|Xu|

ωanc 由兩端節點所含變異在全域基準的群指派後驗算出: 若祖先端的群 CCF 不小於後代端,該邊獲得較高的先驗。 κsplit<0 是狀態數的複雜度罰則。 兩者都必須是的:全域基準的群指派本身有不確定度, 把它當成硬性排序會讓局部推論繼承全域的錯誤而且不留痕跡。

先後關係的先驗只能是軟的:全域基準給的是重疊的後驗,不是一條排序 左上格說明全域基準交出來的東西。含節點 h 的那些變異被指派到群 A, 細胞比例的後驗中心是 0.42、區間 0.24 到 0.60; 含節點 h′ 的那些變異被指派到群 B,中心 0.35、區間 0.13 到 0.57。 兩個區間大幅重疊,重疊的那一段被標了出來。 所以「A 在 B 之前」不是一個事實,是一個機率,這裡是 0.68。 右上格把那個機率變成邊的先驗權重:由 h 指向 h′ 的邊得到 log 0.68, 反方向得到 log 0.32,兩者只差 0.75 nat,沒有任何一個方向被封死。 中間一列把 0.75 nat 跟這個連鎖區段自己的觀測放在一起比: 一次共現觀測,也就是預期 0.2 條而實際看到 5 條,值約 12 nat, 長條長了十六倍。所以局部觀測想翻轉先驗給的方向,輕而易舉。 最下方是硬用的後果:把 0.68 當成 1、0.32 當成 0, 就有三成的機率整個連鎖區段的先後從一開始就反了, 而局部後驗看起來一樣自信,輸出裡看不出來。 讀法是:先驗只調初始傾向,它封不死方向,也不該封。 先後關係的先驗,為什麼只能是軟的 ① 全域基準給的不是排序 群 A = 含節點 h 的變異 群 B = 含節點 h′ 的變異 群 A 0.42 群 B 0.35 0 0.4 0.8 重疊這麼多,所以那是機率,不是事實 ② 那個機率就是這條邊的權重 h h′ log 0.68 log 0.32 兩個方向只差 0.75 nat —— 沒有任何一個方向被封死 0.75 nat 有多大?跟這個連鎖區段自己的觀測比一下 先驗給的方向傾向 0.75 nat 一次共現觀測(預期 0.2 條、實際 5 條) ≈ 12 nat 硬用的後果:錯了不會留下痕跡 把 0.68 當成 1、0.32 當成 0,就有三成的機率整個連鎖區段的先後從一開始就反了,而局部後驗看起來一樣自信。
「軟」在這裡是一個可以算出來的量,不是態度。 全域基準交出來的兩個群,細胞比例的後驗大幅重疊, 所以「A 在 B 之前」是機率 0.68 而不是事實; ωanc 把這 0.68 換成邊的權重,兩個方向只差 0.75 nat。 中間那兩條長條說明這有多小:一次共現觀測(預期 0.2 條、實際 5 條)約 12 nat, 是它的十六倍 —— 局部觀測要翻轉先驗給的方向,一次就夠。 反過來,若把 0.68 當成 1 硬用,三成的連鎖區段先後從一開始就反了, 而且局部後驗看起來一樣自信,輸出裡完全看不出來。

為什麼不做完整的聯合模型:那正是前一篇與整合篇的規格。 本規格刻意切斷回饋(cut posterior),換三件事 —— 連鎖區段之間完全獨立因而可以平行、局部比例不會被全域參數壓平、 以及每一個連鎖區段的結論可以單獨檢視與反駁。 代價是這個後驗不是任何聯合模型的邊際分布, 報告中必須如此稱呼,不得寫成「聯合估計」。

候選模型與輸出旗標

每個連鎖區段的推論目標是 Xu 的後驗。 候選狀態集合由「觀測到的樣式 ∪ 相容性所需的中間狀態」產生, 再對每個候選算邊際似然(ν 依上節先驗積掉),最後正規化成後驗。 Xu 不取單一最大值,而是輸出前若干個候選與各自的質量。

每個連鎖區段除了那份後驗,還帶一組旗標。十二個旗標分成三類, 只有第一類是產出,另外兩類都是「這個連鎖區段的話只能講到哪裡」

  • 與全域基準的關係(四個):REFINES_GLOBALGLOBAL_CONSISTENTLOCAL_ONLY_LINEAGECONTRADICTS_GLOBAL —— 第一個是本規格的主要產出。
  • 證據不足以作答(五個):FLOOR_LIMITEDORDER_UNRESOLVEDBRANCH_UNDERPOWEREDPARTIAL_DOMINATEDSTATE_CAP_REACHED —— 每一個都代表「沒看到」,不代表「不存在」
  • 品質訊號(兩個):INCOMPATIBLEKATAEGIS_DOWNWEIGHTED —— 它們指向地板、誤標或取樣偏性,不是新發現。
旗標全表:十二個旗標各自的觸發條件與讀法(實作時逐欄對照)
旗標觸發條件讀法
REFINES_GLOBAL後驗有足量質量落在「至少兩個節點,其所含變異被全域基準高信心指派到同一群」的候選上本規格的主要產出:全域表達不了的結構
GLOBAL_CONSISTENT後驗主質量的節點與群一一對應局部只是把全域結果再確認一次
LOCAL_ONLY_LINEAGE某節點的比例後驗主要由逃逸成分承接該支在全域頻率譜上沒有對應的峰
CONTRADICTS_GLOBAL局部包含關係與全域基準的鴿籠排序相反診斷量:多半是地板、誤標或全域基準的群指派有問題
FLOOR_LIMITED有效分子數不足以把該狀態與地板分開不是「不存在」
ORDER_UNRESOLVED雙狀態格的跨越深度過低節點存在但先後未定
BRANCH_UNDERPOWERED判為無分支,但分子數與地板不足以看到一支佔 5% 的旁支「沒有分岔」在這裡只代表沒看到
PARTIAL_DOMINATED該連鎖區段的證據主要來自部分覆蓋的分子唯一性多半來自簡約性與先驗,不是資料
INCOMPATIBLE四配子檢定失敗品質訊號,不是新 lineage
KATAEGIS_DOWNWEIGHTED該連鎖區段被判定為局部超突變叢整叢降權為一個有效觀測
STATE_CAP_REACHED候選狀態集合觸及列舉上限不得在截斷後仍回報通過
一個連鎖區段定得下來什麼、定不下來什麼,以及每一項各由哪個旗標承接 左格是這個連鎖區段定得下來的三件事,全部由該連鎖區段自己的分子決定: 有哪幾種狀態,這裡是全參考、只帶第一個變異、以及兩個雙變異狀態共四種; 誰在誰之前,只帶第一個變異的那個狀態是另外兩個的祖先; 以及每一種各佔多少分子,依序是 0.66、0.19、0.05、0.07。 三件事都落在「一個連鎖視窗、一條染色體拷貝」之內。 右格是同一個連鎖區段定不下來的四件事,每一項後面掛著承接它的旗標: 某一支對應到哪一個全域群,掛 UNLINKED_TO_GLOBAL; 分子比例換算成細胞比例是多少,缺純度與拷貝數,掛 CELL_SCALE_UNAVAILABLE; 有沒有第四支被錯誤地板吃掉,掛 FLOOR_LIMITED; 某個沒被觀測到的狀態是真的不存在還是沒抽到,掛 ORDER_UNRESOLVED。 最下方是這張圖的重點:旗標不是免責聲明,是欄位。 右格每一項都有名字,所以「這條流程對多少比例的資料拒絕作答」查得出來; 被沉默補齊的東西在統計量上是隱形的。 一個連鎖區段交出去的那一列,說了什麼、沒說什麼 定得下來 —— 用它自己的分子 000 100 110 101 有哪幾種狀態 000 · 100 · 110 · 101 誰在誰之前 100 → 110 · 100 → 101 各佔多少分子 .66 · .19 · .05 · .07 三件事都在「一個連鎖視窗、一條拷貝」之內 —— 沒有用到任何連鎖區段以外的東西 定不下來 —— 每一項各有一個名字 110 那一支對應哪一個全域群? UNLINKED_TO_GLOBAL 0.05 換算成細胞比例是多少? CELL_SCALE_UNAVAILABLE 有沒有第四支被錯誤地板吃掉? FLOOR_LIMITED 沒看到 011,是不存在還是沒抽到? ORDER_UNRESOLVED 旗標不是免責聲明,是欄位 右格每一項都有名字,所以「對多少比例的資料拒絕作答」查得出來;被沉默補齊的東西在統計量上是隱形的。
把十二個旗標讀成一張對照:左格三件事全部由該連鎖區段自己的分子定出來, 右格四件事它定不下來 —— 而每一項各對應一個旗標。 兩格的分界就是上述估計範圍:問題只要越出「一個連鎖視窗、一條拷貝」, 答案就不在左格。這也是為什麼旗標必須是欄位而不是註解: 右格每一項都有名字,「這條流程對多少比例的資料拒絕作答」才查得出來。

REFINES_GLOBAL 是這條路線的操作型定義: 「比全域方法解析度高」這句話在這裡等於「有多少個連鎖區段掛上這個旗標, 而且顯著多於虛無分布」。後文「驗收設計」定義這個虛無分布的產生方式。 沒有這個定義,前面所有論述都不可檢驗。

跨連鎖區段的彙整邊界

前一篇的規則是「可以合併的是機率,不是標籤」。 在本規格下更嚴格:連機率都不合併,因為 νu 是每個連鎖區段自己的參數, 不同連鎖區段的後驗沒有共同的參數可以相乘。局部拓撲當然更不可以拼接成一棵大樹 —— 局部標籤只在該連鎖區段內有定義。

跨連鎖區段的產出可否條件
連鎖區段層級的目錄(每個連鎖區段一列 + 旗標)本規格的主要交付物
樣本層統計量:REFINES_GLOBAL 的個數、局部分岔比例、節點數分布必須同時報告驗收設計的虛無分布
把局部樹拼成一棵全基因體的樹不可標籤跨連鎖區段無定義,且不同連鎖區段看的位點集合不同
把成對的祖先/互斥關係匯出給全域方法有條件使用它的全域擬合必須排除 S2 的 read,否則就是重複計數
把節點數加總當成全基因體的群數不可節點是分子狀態,且各連鎖區段互不對齊
三層資訊:全基因體只有邊際、phase block 知道 cis/trans、read 跨距內才數得到共現 三層由下而上代表資訊愈來愈完整,但涵蓋範圍愈來愈小。 最下層是全基因體,任何兩個位點之間都只有各自的邊際頻率,沒有共現關係。 中間層是 phase block,典型長度一到五個 Mb,在這個範圍內兩個位點是在同一條染色體 還是分開在兩條,是已知的,即使沒有任何一條 read 同時蓋到它們。 最上層是 read 跨距,典型二十個 kb,只有在這個範圍內才數得出「同時帶有兩個變異的分子有幾條」。 右側標出每一層在一個典型樣本裡各有多少筆可用觀測: read 跨距內的聯合計數約數百到數萬筆,phase block 內的 cis/trans 判定約多一百倍, 全基因體邊際則涵蓋全部變異。 讀法是:read 長度限制的是「數得到共現」,phase block 限制的是「知道在不在同一條染色體上」, 這兩件事的範圍差了兩個數量級,常被混為一談。 三層資訊:範圍愈小,資訊愈完整 ③ 全基因體 只有邊際頻率(VAF、β)—— 傳統路線整個住在這一層 全部變異 數百萬筆 ② phase block(1–5 Mb) cis/trans 已知,但數不到共現分子數 cis/trans 判定 約 7 萬筆 ① read 跨距(20 kb) 聯合計數看得到 聯合計數 數百到數萬筆 資訊最完整 涵蓋最廣 最常見的誤解:把 ① 與 ② 當成同一件事。read 長度限制的是「數得到共現」, phase block 限制的是「知不知道在同一條染色體上」—— 兩者範圍差約一百倍。
read 跨距限制的是能否計數共現, 限制的是能否確定兩個位點是否位於同一條染色體。 本規格的每一列都只能自最上層產生,而該層亦為三者中範圍最小者 —— 這就是「不可拼接成一棵大樹」的物理來源,不是保守的選擇。

輸入與前處理契約

局部模型假設連鎖區段與稀疏狀態表已經建立。 前處理沿用既有工具,但其輸出直接決定候選狀態與覆蓋遮罩; 因此工具版本、品質門檻、HP3 處理與連鎖定義都是模型輸入契約的一部分。

下列參數會改變狀態表、連鎖區段與候選集合,因此必須與結果一起記錄, 以保證前處理可重現。

工具鏈與必須釘住的參數

既有的鏈是:正常樣本 germline calling() → 定相 → 腫瘤樣本 somatic calling() → 體細胞標記。兩個步驟各列了兩個可選工具, 但沒有記錄實際用了哪一個、是否取交集,也沒有任何版本號。 這不是文件疏漏,是重跑不出同一份狀態表的直接原因 —— 狀態表換了,連鎖區段、候選集合與每一個下游數字都跟著換。

清單:五項必須釘住的設定,以及各自不釘住會壞在哪
必須釘住的不釘住的後果
germline/somatic caller 的身分、版本,以及是否取交集候選位點集合改變 → 連鎖區段改變 → 全部數字不可比
判定「這條 read 覆蓋這個位點」的最低 base quality 與 mapping quality決定 XRA 的分界,直接改變覆蓋遮罩
沒有 HP 標籤的 read 怎麼處理丟掉與併入是兩種方向相反的偏誤,既有實作未載明
兩個 sSNV 要有幾條共同覆蓋的 read 才算連鎖既有實作沒有下限,所以一條 read 就能開出一個連鎖視窗
一個連鎖視窗要有多少跨越深度才保留同上;下限缺席時,證據極薄的連鎖區段與紮實的連鎖區段在計數上等重
單倍型標籤的五個值,以及連鎖視窗如何由讀序連鎖的傳遞閉包長出來 上排三格說明體細胞單倍型標記輸出的標籤詞彙恰好只有五個值。 左格是 germline 的兩條:HP1 與 HP2,它們是定相軟體給的兩條單倍型, 標號只在這個 phase set 內有效。 中格是可歸因的體細胞改變:HP1-1 與 HP2-1, 表示該分子的體細胞改變可歸因於第一個數字所指的那條單倍型,詞彙到此為止。 右格是無法歸因的 HP3:兩條單倍型都歸不上去, 所以「一個連鎖區段只看得到一條染色體拷貝」這個假設對它並不成立。 中排說明連鎖視窗的定義:在同一個 phase set 內, 只要有某一條讀序同時覆蓋至少兩個體細胞變異,該段就構成一個連鎖視窗; 若另一條讀序疊到已連鎖的位點又碰到新的位點,兩段就合併成更長的連鎖視窗。 圖中左邊三條讀序彼此重疊,把五個位點連成一個連鎖視窗; 右邊兩條之中,一條連起兩個位點形成第二個連鎖視窗,另一條只覆蓋單一位點因而沒有連鎖。 因此位點數 k 數的是位點而不是讀序,一個由數條讀序接起來的連鎖視窗可能一條全跨分子都沒有。 最下方是兩件容易被誤讀的事: 後綴不是獲得序列,1-1 是一個標籤而不是「1 再加一步」,不認得的標籤一律排除; 搜尋成本由 q 決定而不是 k,q 是該家族內真的有變化的位點數,恆有 q 小於等於 k。 標籤只有五個值,連鎖視窗由讀序連鎖長出來 germline 兩條 HP1 HP2 定相軟體給的兩條 單倍型,標號只在 這個 phase set 內有效 可歸因的體細胞改變 HP1-1 HP2-1 體細胞改變可歸因於 第一個數字那條單倍型 詞彙就到此為止 無法歸因 HP3 兩條都歸不上去 「一個連鎖區段只看得到 一條拷貝」對它不成立 連鎖視窗 = 同一個 phase set 內、讀序連鎖的傳遞閉包 W₁ · k=5 W₂ · k=2 沒有連鎖 k 數的是位點不是 read —— 一個由數條 read 接起來的連鎖視窗,可能一條全跨分子都沒有。 兩件會被誤讀的事 後綴不是獲得序列:1-1 是一個標籤,不是「1 再加一步」;不認得的標籤一律排除。 搜尋成本由 q 決定不是 k:q 是該家族內真的有變化的位點數,恆有 q ≤ k。
上排是標籤詞彙恰好五個值,其中 HP3 是本規格的真問題。 下排是連鎖視窗的定義:不是固定寬度,是讀序連鎖的傳遞閉包。 最下方兩件事看似細節,但兩者都會安靜地改變結論。

標籤詞彙只有五個值,而 HP3 是一個真的問題

輸出的詞彙恰好121-12-13 五個值。 1-1 表示這條分子屬於 germline 單倍型 1,且帶有可歸因於它的體細胞改變; 3 表示體細胞改變無法歸因於任何一條 germline 單倍型。 沒有更長的後綴,所以後綴不是一條獲得序列 —— 不可以把它讀成「先 1 再加一步」。不認得的標籤一律排除在所有分組之外。

HP3 破壞的正是「一個連鎖區段只看得到一條拷貝」

估計單位的核心假設是:連鎖區段內部不必處理兩條拷貝的混合。 HP3 的分子不滿足它,因為它們的體細胞改變兩條都歸不上去。 三種處理各有代價,必須明選其一並在輸出載明實際被影響的量: 整批排除(最保守,但在 LOH 與高拷貝區會丟掉大量證據); 當成獨立的第三個連鎖區段(形狀對,但那個連鎖區段的 CCF 換算沒有分母); 當成缺失並在兩個家族之間邊際化(統計上正確,成本最高)。 既有實作採第一種,但沒有報告被排除的量 —— 沒有那個量, ν 的分母就是錯的,而且錯得沒有跡象。

連鎖視窗不是固定寬度,是讀序連鎖的傳遞閉包

既有實作的連鎖視窗定義比「固定寬度 20 kb 視窗」精確得多,也更該沿用: 在同一個 內,凡有某條 read 同時覆蓋至少兩個 sSNV,該段即為一個連鎖視窗; 若另一條 read 疊到已連鎖的位點又碰到新的位點,兩段合併成更長的連鎖視窗。三個推論:

  • k 數的是位點不是 read,所以一個由數條 read 接起來的連鎖視窗 可能一條全跨分子都沒有 —— 這正是覆蓋遮罩必須逐條保留的原因;
  • 連鎖視窗不跨 phase set,因此所有狀態比較都在同一個定相框架內 —— 這是「不同 block 的 H1 不是同一條」那條警告的實作形式;
  • 一個連鎖視窗依家族切成至多兩個連鎖區段,所以連鎖區段數大於連鎖視窗數, 兩者不可混用為同一個分母。
兩種「沒看到」:真的沒有那種細胞,與沒有 read 跨得到 左邊是跨越深度隨距離衰減的曲線。 兩個位點靠得愈近,同時蓋到它們的 read 就愈多;距離接近讀長時,跨越深度掉到接近零。 圖上標出真實例子的三個位點對,跨越深度分別是二十、十七與三。 右邊說明問題所在:演算法看到某個狀態沒有出現,有兩種完全不同的原因。 第一種是真的沒有細胞帶著那個組合,那是生物學結論。 第二種是無 read 同時覆蓋該二位點,屬於未觀測。 目前的做法把兩者當成同一件事,於是把看不到當成不存在。 最下方是修正方式:每一個位點對先用自己的跨越深度算檢定力, 檢定力不足的對標成未定,不允許它去約束樹的形狀,並且在報告裡列出來。 讀法是:缺席要有足夠的觀測撐著,才能當成證據。 兩種「沒看到」 跨越深度隨距離掉得很快 兩個位點的距離 → 跨越深度 n=20 n=17 n=3 檢定力門檻 同樣是「這個狀態沒出現」 ① 真的沒有細胞帶著那個組合 這是生物學結論,可以拿來決定樹的形狀。 ② 無 read 同時覆蓋該二位點 這只是看不到,不能拿來決定任何事。 目前的做法把兩者當成同一件事 修正:每一對先過檢定力這一關 用該位點對自己的跨越深度算檢定力;不足的標成「未定」, 不允許它的「缺席」去約束樹,並且在報告裡逐條列出來。 預期效果:目前被歸進「單層無分支」的連鎖區段,會有一部分正確地移到「未解析」。
「此狀態未出現」有兩種成因:該類細胞確實不存在, 或無 read 跨越該對位點。既有實作沒有設連鎖下限,所以兩者混在一起; 似然寫法讓第二種成因自動只貢獻很少的證據,不需要門檻就分得開 —— 但前提是那條 read 有沒有跨過去,必須記下來而不是先併掉。
定義補充:前處理還要交出一個量 —— 作用位點數 q(搜尋成本要用它,不是 k

前處理交出去的不只是連鎖區段本身,還有一個決定推論成本的量, 而它必須在這裡就算出來,因為它只有看得到狀態表的時候算得出來: 作用位點數 q, 即在該家族內真的有變化的位點數。某個位點若該家族每條 read 都是同一個 allele, 它不產生任何座標,故 qk;用 k 估搜尋成本會系統性高估。 本頁凡涉及狀態空間、頂點、遮罩與上限之處,一律以 q 為準。

連鎖率才是真正的閘門,而它逐樣本差很大

中篇以 Poisson 模型推得「約 9.5% 的變異落在多變異連鎖區段」。 實測不是這樣,而且它的變動幅度是這條路線最重要的一個經驗事實。 既有實作在七個資料集上量到:落入任一連鎖視窗的 sSNV 位點比例由 40.5% 到 98.1%,超過兩倍,而且不依癌別排序 —— 四個乳癌資料集自己就橫跨了幾乎整個範圍。

落入任一連鎖視窗的體細胞變異比例,在七個資料集之間差超過兩倍 七條橫向長條,每條是一個資料集,長度是該資料集中落入任一分析連鎖視窗的 體細胞單核苷酸變異位點比例,右側另列該資料集的變異總數。 由低到高分別是:HCC1395_NYGC 四成零五、HCC1395_HKU 四成五七、 H1437 六成九四、H2009 九成二五、HCC1937 九成七七、HCC1954 九成七七、 COLO829 九成八一。最低與最高之間相差超過兩倍,而且不依癌別排序: 四個乳癌資料集自己就橫跨了幾乎整個範圍。 最前面兩條是同一株細胞株經過兩條不同流程的結果,兩者只差 5.2 個百分點, 所以其他資料集之間五十個百分點以上的落差不是技術雜訊。 因此連鎖率必須逐樣本先算再決定要不要跑,它是分層報告的第一個分層變數, 而不是事後的說明。 連鎖率才是閘門,而它逐樣本差很多 資料集 落入任一連鎖視窗的 sSNV 位點比例 比例 sSNV 總數 HCC1395_NYGC 40.5% 79,739 HCC1395_HKU 45.7% 79,687 H1437 69.4% 77,080 H2009 92.5% 154,465 HCC1937 97.7% 18,690 HCC1954 97.7% 22,400 COLO829 98.1% 37,788 0 100% 同一株細胞株、兩條流程,只差 5.2 個百分點 所以其餘五十個百分點以上的落差不是技術雜訊。連鎖率要逐樣本先算,它是第一個分層變數。
同一株細胞株經兩條流程只差 5.2 個百分點, 所以其餘五十個百分點以上的落差不是技術雜訊。 Poisson 模型給的單一數字在這裡沒有用:真正該做的是逐樣本先算這個比例。

連鎖進來之後,多數連鎖區段只有兩個位點

第二個實測事實同樣改變了規格該長什麼樣:k=2 是七個資料集中六個的最大單一類, 但集中程度由 76.7% 到 16.8%,只有 H2009 把 k 一路填到上限 12。 換句話說,這條路線大部分的產出,是一對位點之間的關係,而不是多層系譜。

連鎖區段實際長什麼樣:多數只有兩個位點,多數形狀是鏈 左欄是每個連鎖視窗含幾個體細胞變異位點。 位點數等於二是七個資料集中六個的最大單一類,但集中程度差很多: HCC1395_NYGC 有七成六七的連鎖視窗只有兩個位點,H2009 只有一成六八, 而且只有 H2009 把位點數一路填到上限十二。 因此多數連鎖區段講的是一對位點之間的關係,要談多層系譜得位點數至少三個。 右欄是重建出來的形狀組成。無分支的類別在每個資料集都佔多數, 由七成一七到八成八五,有分支的類別到處都是少數。 原因寫在下方的方塊裡:一條鏈只要一個觀測到的樣式就建得起來, 一個分岔卻要兩個互不包含的樣式, 所以分岔稀少反映的是這個連鎖區段帶了多少不同的樣式證據,不只是底下的結構。 最下方指出這件事限制了什麼、又沒有限制什麼: 它限制了深度,多數連鎖區段只講得出一對位點的先後; 但它沒有限制可偵測性,前面那個證據預算在兩個位點時就成立, 位點數增加買到的是深度,不是分不分得出來。 連鎖區段實際長什麼樣 一個連鎖視窗有幾個位點 k=2 在七個資料集中有六個是最大類 HCC1395_NYGC 76.7% H2009 16.8% 深色 = 只有兩個位點的連鎖視窗 只有 H2009 把 k 一路填到上限 12 多數連鎖區段講的是一對位點的關係 重建出來的形狀 無分支類在每個資料集都佔多數 無分支 71.7%–88.5% 有分支:少數 為什麼分岔這麼少 一條鏈只要一個觀測到的樣式 一個分岔要兩個互不包含的樣式 稀少反映的是樣式證據量 這限制了什麼,又沒有限制什麼 限制了深度:多數連鎖區段只講得出一對位點的先後,談不上多層系譜。 沒有限制可偵測性:前面那個證據預算在 k=2 就成立 —— k 變大買的是深度。
左欄:多數連鎖視窗只有兩個位點,但集中程度差很多。 右欄:無分支的形狀在每個資料集都佔七到九成,而且原因是機械性的 —— 一條鏈只要一個樣式,一個分岔要兩個互不包含的樣式。

這件事必須誠實地放進定位,但它不動搖前面那個證據預算,兩者常被混在一起: k=2 佔多數限制的是深度(多數連鎖區段只講得出一對突變的關係), 不是可偵測性 ——「預期 0.2 條、實際 5 條」在 k=2 同樣成立, 而分辨兩支等值 CCF 的 lineage 本來就只需要一對位點。

逐項對照:k=2 佔多數改變了什麼、沒有改變什麼
k=2 佔多數影響
限制了深度多數連鎖區段只講得出「這兩個突變共不共存、誰在前」,談不上多層局部系譜
沒有限制可偵測性開頭那個「預期 0.2 條、實際 5 條」在 k=2 同樣成立 —— 它比的是「某個共現狀態存在」與「只有錯誤地板」,與位點數無關
REFINES_GLOBAL不影響:分辨兩支等值 CCF 的 lineage 只需要一對位點的共現或互斥
對交付物的描述要改:主要交付物是「全基因體的成對分子關係目錄」,深層局部系譜是其中的少數

連帶的還有兩個直接後果。其一,caller 的局部系譜特徵覆蓋率不是固定的一成 —— 在高連鎖率的樣本上它可以覆蓋大部分候選,在低連鎖率的樣本上則接近無效。 其二,連鎖率是分層報告的第一個分層變數,必須在跑之前先算、 並用它決定這個樣本值不值得跑,而不是事後拿來解釋結果。

推論演算法與計算預算

  1. 前處理。依連鎖區段輸出稀疏狀態表,每列為 (連鎖區段、phase set、位點集合、覆蓋遮罩、字母樣式、計數)。遮罩不得先併掉。
  2. 就地校準。由 germline-only 連鎖區段、matched normal 與已知重複區估 逐位點錯誤率與家族誤標率,成對估計並分層借力。
  3. 全域基準。只用 S1,輸出原子、權重、純度與 φ^, 並斷言與 S2 的位點交集為空。
  4. 候選狀態集合列舉。先取計數高於地板的樣式,再補相容性所需的中間狀態。 狀態空間是 2qq=2 時 4 格,q=12 時 4,096 格 —— 所以上限必須存在。既有實作用的兩個上限與它們的實際行為見下。
  5. 邊際似然。對每個候選,ν 依上述全域頻率先驗積掉。 先驗是 Beta 混合、似然是 multinomial,非共軛, 故以每個候選數十個 quadrature 或 Monte Carlo 點計算並回報數值誤差。 這一步取代既有實作的 read-AF 排序:既有那個分數不只排序,它會先砍掉 30.03% 的最小成本解(見證據節),而且它沒有任何拷貝數校正 —— 當一個連鎖區段的每個位點都落在最高 read-AF 時,分數退化、全部並列。 換成邊際似然之後,不同形狀比較的是對同一組計數的解釋能力, 而且候選整組留著、由後驗質量表達,不先砍。
  6. 全域基準不確定度的傳遞。由全域基準後驗抽 S(ρ^,c^,ω^),逐組重算再平均。 S 過小會低估區間寬度,須在收斂契約中登錄。
  7. 輸出。每個連鎖區段一列,含前若干個 Xu 候選與質量、 節點比例的區間、旗標與有效分子數。
read-AF 差值總和有形狀偏誤:星狀永遠得零分,串接永遠贏 目前排序候選樹的分數,是把樹上每一組「祖先與後代」的 read-AF 相減再全部加起來。 問題在於相加的項數本身由樹的形狀決定。 圖上三個形狀由左到右:三個突變串成一條線、一條線加一個分支、三個都從根長出來的星狀。 它們的祖先後代配對數分別是三、一、零。 因此串接的分數是第一個減第三個的兩倍,中間的形狀只有一項, 而星狀因為沒有任何祖先後代配對,分數恆等於零。 結論寫在下方:只要 read-AF 不完全相等,串接的分數就是正的, 而星狀永遠是零,所以當最小成本候選裡同時有這兩種形狀時,分支永遠贏不了。 觀測到的最多的類別正好是多層無分支,佔三成八到五成三, 這個偏誤剛好往那個方向推。 最下方是替代方案:改用機率模型比較各個形狀解釋同一組 read-AF 的能力, 而不是加總差值。 讀法是:分數要能跨形狀比較,才不會讓形狀自己決定勝負。 分數的形狀偏誤 分數=把樹上每一組「祖先→後代」的 read-AF 相減,全部加起來。項數由形狀決定。 串接 祖先後代配對:3 組 2(a₁ − a₃) 一線加一分支 祖先後代配對:1 組 a₁ − a₂ 星狀(全分支) 祖先後代配對:0 組 恆等於 0 只要 read-AF 不完全相等,串接就是正分,星狀恆為 0 —— 所以最小成本候選裡同時有這兩種形狀時,分支永遠贏不了 觀測到最多的類別正好是「多層無分支」(38–53%),而這個偏誤剛好往那個方向推。 替代方案:用機率模型比較各形狀解釋同一組 read-AF 的能力,而不是加總差值。
三種形狀的祖先後代配對數分別為 3、1、0。 串接得到正分,星狀恆為 0 —— 所以最小成本候選中同時存在此二形狀時,分支形狀恆不獲選。 而分岔正是拆峰能力唯一的來源,所以這個偏誤系統性地移除的, 剛好就是本規格最想留下來的那一類。
局部推論的七步推論流程、兩個上限與實測成本 左側是局部推論的七個步驟,由上而下:輸出稀疏狀態表、就地校準兩個錯誤率、 跑全域基準、列舉候選狀態集合、算邊際似然、把全域基準的不確定度傳下來、輸出。 其中第四步與第五步用深色標出,因為只有這兩步與既有實作不同; 第三步是沿用,不是本規格的貢獻。 右側三張卡分別說明三件會影響實作的事。 第一張是兩個上限:搜尋節點上限一千,實測有一成二的連鎖區段撞到它; 候選家族上限十萬,實測從來不觸發。兩者都是工程設定不是數學結果, 必須寫進 provenance,而撞到上限的連鎖區段要標記棄權並且仍然留在分母裡。 第二張說明第五步取代掉的是既有實作的 read-AF 排序: 那個分數會先砍掉三成的最小成本解,而且沒有任何拷貝數校正; 改成邊際似然之後候選整組留著,由後驗質量表達。 第三張是既有實作量過的成本:九萬多個連鎖區段兩分鐘、記憶體不到四百 MB, 而成本的變異來自作用位點數的分布,不是連鎖區段的總數 —— 跨資料集差約十九倍。 讀法是:這條流程沒有全域取樣包在外層,所以它是尷尬平行的, 成本隨連鎖區段數線性成長。 局部推論的七步,以及成本落在哪裡 深色的兩步才是與既有實作的差異;③ 是沿用,不是貢獻 局部推論的七步 ① 前處理:稀疏狀態表,覆蓋遮罩不得先併掉 ② 就地校準:逐位點錯誤率與家族誤標率,成對估 ③ 全域基準:只用 S₁,斷言與 S₂ 的位點交集為空 ④ 候選列舉:狀態空間 2^q,q=12 時 4,096 格 ⑤ 邊際似然:ν 依先驗積掉,候選整組留著不先砍 ⑥ 傳遞全域基準的不確定度:抽 S ≈ 50 組逐組重算 ⑦ 輸出:候選與質量、節點比例區間、旗標、分子數 ④⑤ 是方法核心;其餘五步既有實作都已具備。 沒有任何一步把全域取樣包在外層。 兩個上限與棄權 κ_node = 1,000 → 12.47% 的連鎖區段撞到(實測) κ_fam = 100,000 → 實測從不觸發;兩者皆為工程設定 ⑤ 取代舊的 read-AF 排序 舊分數先砍掉 30.03% 的最小成本解 而且沒有任何拷貝數校正 新做法:候選整組留著 由後驗質量表達,不先砍 成本:既有實作已經量過 98,955 段 · 122.03 s · 356.7 MB 13,773,565 個搜尋節點,平均 139.2 跨資料集差 19 倍:292.0 對 15.5 要控成本,該控的是 q 的分布,不是連鎖區段總數 —— 成本隨連鎖區段數線性成長。
左欄七步由上而下。只有 ④⑤ 與既有實作不同, ③ 是沿用而非貢獻,其餘四步既有實作都已具備。 右上是兩個上限的實測行為:節點上限有 12.47% 的連鎖區段撞到,家族上限從不觸發 —— 兩者都是工程設定,不是關於 q 的數學結果。 右中是 ⑤ 取代掉的東西:舊的 read-AF 分數會先砍掉 30.03% 的最小成本解且無拷貝數校正。 右下是實測成本;成本的變異來自 q 的分布,不是連鎖區段總數

上限、棄權,以及棄權必須留在分母裡

既有實作用兩個上限:搜尋節點上限 κnode=1,000 與候選家族大小上限 κfam=100,000。實測結果是 只有節點上限真的會觸發 —— 帶變異的連鎖區段中有 12.47% 撞到它。 兩個上限都是這套輸入與這台機器上的工程設定,不是關於 q 的數學結果; 換一台機器或換一個搜尋順序,數字就會移動,所以它們必須寫進 provenance。

撞到上限的連鎖區段要標記並計數,不可以丟掉

這是既有實作做對、而且值得整條抄過來的一件事: 撞上限的連鎖區段得到一個明示的棄權標記,而且仍然佔據每一個比例的分母。 被無聲丟棄的連鎖區段在完全相同的統計量上是隱形的 —— 報告中將無法判斷流程對多少比例的資料拒絕作答。 另有一條連帶的紀律:同一份報告至少有兩個分母 (全部連鎖區段、帶變異的連鎖區段),兩者回答不同的問題, 任何百分比都必須帶著自己的分母出現。

計數補充:同一個候選底下有幾棵等價的樹,不必抽樣(它沿頂點分解)

上限管的是「候選列不列舉得完」;另一件不同的事是列舉完之後, 同一個候選底下有幾棵等價的樹。 這個數目不需要抽樣,因為它沿頂點分解。 給定一個候選的頂點集合 N,每個非根頂點各自獨立地選一個父節點, 可選數就是它在 N 誘導子圖上的入度:

CT(N)=vN,vrootdN-(v)

既有實作用這條式子精確計算並列數。在本規格下它有第二個用途: 它是簡約先驗在給定狀態集合下的正規化常數, 所以「這個候選有多少種等價寫法」不會被誤算成「這個候選有多少證據」。

成本不是這條路線的瓶頸

既有實作已經量過:整個基因體的局部重建是兩分鐘的工作,記憶體不到 400 MB。 本規格在它之外只加了每個候選的積分與全域基準的 S50 次重算, 仍然遠低於任何全域取樣。要注意的只有一件事: 成本的變異來自 q 的分布而非連鎖區段數 —— q 大的資料集吃掉絕大部分節點,所以要控成本,該控的是 q 的分布, 不是連鎖區段總數。

實測數字:98,955 個連鎖區段的秒數、記憶體與跨資料集 19 倍差異

既有實作處理 98,955 個連鎖區段, solver 122.03 秒、census 58.96 秒、峰值常駐記憶體 356.7 MB、 展開 13,773,565 個搜尋節點 —— 平均每個連鎖區段 139.2 個節點, 但跨資料集差約 19 倍(H2009 每連鎖區段 292.0,HCC1395_NYGC 15.5)。 那個 19 倍就是上面「變異來自 q 的分布」的量: H2009 正是唯一把 k 一路填到上限 12 的資料集。

關鍵仍在於沒有全域取樣包在外層,所以整份工作是尷尬平行的, 成本隨連鎖區段數線性成長。這是放棄全域目標換來的直接好處, 也是它與整合篇那條路線在工程上的主要差別。

外部驗證:Somatic caller

局部系譜沒有完整金標準,因此另以 somatic calling 作為具有 truth set 的外部驗證終點。 此路線不參與局部系譜的核心推論;它只檢驗局部 lineage 一致性是否在 caller 既有訊號之外提供增量。

錨點與候選必須分開

這裡有一個明顯的循環:局部系譜要靠可信的 somatic 變異才建得起來, 而我們想幫的正是那些不可信的候選。解法是把兩者分成兩個角色:

角色來源用途
錨點germline het 位點 + 通過嚴格門檻的 somatic 呼叫建構該連鎖區段的局部系譜
候選caller 在寬鬆門檻下的候選集合被評分的對象,不進錨點集合

因此每個候選 m 的特徵,都是對把它拿掉之後建成的局部系譜 Lu(m) 算出來的。這條留一法不是保險,是必要條件: 否則特徵會直接把「這個候選有幾條 ALT read」抄一遍, 在訓練時看起來很有效,在真正的低頻變異上完全沒有幫助。 錨點也不得取自 truth set —— 那是最典型的

錨點與候選的角色分離,以及留一法算出來的特徵 這張圖解決一個循環:局部系譜要靠可信的體細胞變異才建得起來, 而我們想幫的正是那些不可信的候選。 左欄把兩者分成兩個角色。錨點是 germline 雜合位點加上通過嚴格門檻的體細胞呼叫, 它們負責建樹;候選是 caller 在寬鬆門檻下給出的東西,它是被評分的對象, 絕對不能進入錨點集合,否則循環就回來了。 中欄是把候選拿掉之後才建起來的局部系譜,四個節點分別是全零、一零零、一一零、一零一。 候選不在這棵樹裡,所以算出來的特徵不可能抄它自己的資料。 右欄是把候選的 ALT 分子投到這棵樹的節點上之後看到的兩種形狀。 真的低頻體細胞變異來自某一支既有的 lineage, 所以它的支持分子會集中在一個節點上,巢狀對數似然比因而很高。 定序錯誤與嵌合對局部系譜完全無所謂,分子會散落在各節點之間,對數似然比接近零。 讀法是:不做留一,特徵就退化成把「這個候選有幾條 ALT read」抄一遍 —— 訓練時看起來很有效,在真正的低頻變異上完全沒有幫助。 留一法:候選被拿掉之後,才建局部系譜 錨點建樹、候選被評分。兩個角色混用,循環就回來了 ① 兩個角色分開 錨點 germline 雜合位點 + 通過嚴格門檻的 somatic ↓ 只有錨點建樹 候選 m caller 寬鬆門檻下的候選 被評分的對象 候選進錨點集合, 錨點也不得取自 truth set。 後者是最典型的資料洩漏。 ② 拿掉 m 之後建樹 000 100 110 101 這棵樹叫 L(−m) m 不在裡面 所以特徵抄不到它自己 ③ m 的 ALT 分子落在哪 真的低頻 somatic:集中在一個節點 000 100 110 101 llr_nested 高 → 判為真 定序錯誤或嵌合:散落各節點 000 100 110 101 llr_nested ≈ 0 → 判為錯誤 不做留一,特徵等於把「這個候選有幾條 ALT read」抄一遍 —— 訓練時好看,低頻變異上沒有幫助
左欄是角色分離:錨點建樹,候選只被評分, 候選不進錨點集合,否則那個循環就回來了。 中欄是把候選拿掉之後才建起來的 Lu(m) —— 候選不在裡面,所以特徵抄不到它自己的資料。 右欄是把候選的 ALT 分子投到這棵樹上之後的兩種形狀: 真的低頻 somatic 來自某一支既有 lineage,分子集中在一個節點; 定序錯誤與嵌合對局部系譜無所謂,分子散落各節點llrnested 因而接近零。 這就是 llr_nestedlineage_coherence 量的東西。

特徵群與可用範圍

給 somatic caller 的三組特徵,以及各組的覆蓋率與比較基準 三條橫向的帶狀區塊,由上而下是基礎特徵群、單倍型特徵群與局部系譜特徵群。 每一組左邊寫覆蓋率並畫成長條,中間寫它提供什麼,右邊寫它要打敗的基準。 基礎特徵群所有候選都有:把樣本層的 CCF 原子換算成這個位點的期望 VAF 格點, 再加上就地估到的錯誤地板;要打敗的基準是 caller 自己從訓練資料學到的 VAF 分布。 單倍型特徵群是有 germline 錨點的候選,佔多數:ALT 分子集中在哪一條 germline haplotype 上, 以及該單倍型的深度脈絡;要打敗的基準是 ClairS 已經有的相位通道。 局部系譜特徵群只有連鎖進某個連鎖視窗的候選才有,這個比例逐樣本由四成到九成八: ALT 分子與既有 somatic lineage 的一致性、巢狀似然比與四配子相容性; 這一組沒有現成基準,是本規格新增的。 最下方兩個方塊是兩條必須寫進實作的規則。 左邊是留一法:局部系譜只由錨點建成,被評分的候選本身不得列入錨點, 否則特徵會直接把「這個候選有 ALT read」抄一遍。 右邊是缺席不是零:沒有連鎖進連鎖視窗的候選就沒有局部系譜特徵,每一組都要有「這一組不可用」的指示通道, 而且必須證明在不可用的子集上沒有變差;驗收時三組各報一次 PR,不可只報總體。 局部系譜特徵群那條長條畫成一個範圍,淺色是上界九成八,深色是下界四成, 因為這個覆蓋率逐樣本差超過兩倍,不是一個固定值。 三組特徵:覆蓋率不同,要打敗的基準也不同 基礎特徵群 所有候選都有 樣本層的 CCF 原子換算成 這個位點的期望 VAF 格點 以及就地估到的錯誤地板 要打敗的基準 caller 自己學到的 VAF 分布 單倍型特徵群 有 germline 錨點的候選 ALT 分子集中在哪一條 germline haplotype 上 以及該單倍型的深度脈絡 要打敗的基準 ClairS 已經有的相位通道 局部系譜特徵群 連鎖進連鎖視窗的候選 40.5% – 98.1% ALT 分子與既有 somatic lineage 的一致性、巢狀似然比 與四配子相容性 沒有現成基準 這一組是本規格新增的 留一法:候選不進錨點 局部系譜由錨點建成,被評分的候選 本身不得列入錨點;否則特徵會直接 把「這個候選有 ALT read」抄一遍。 錨點 = germline het + 嚴格門檻的 somatic 缺席不是零 沒連鎖進連鎖視窗就沒有局部系譜特徵。每一組 都要有「不可用」的指示通道,而且 必須證明在不可用的子集上沒有變差。 驗收時三組各報一次 PR,不可只報總體。
三組特徵的覆蓋率與要打敗的基準都不同。 單倍型特徵群是 ClairS 已有的相位通道,是基準而不是貢獻; 局部系譜特徵群的覆蓋率逐樣本由 40.5% 到 98.1%,不是固定的一成; 但它是唯一沒有現成替代的特徵群。

局部系譜特徵群的核心是 llr_nested,而它的邏輯與 ClairS 的相位通道相同、 只是改用更局部的分類單位:真的低頻 somatic 變異來自某一支既有的 lineage, 所以它的支持分子應該整齊地落在該節點的後代位置上; 定序錯誤與嵌合則對局部系譜完全無所謂,會散落在各節點之間。 差別在於相位通道分的是兩條 germline 單倍型, 局部系譜特徵分的則是同一條單倍型內部的數支 somatic lineage。

欄位全表:三組特徵的每一個欄位與定義(實作時逐欄對照)
特徵群欄位意義
基礎vaf_grid_dist觀測 VAF 與「全域基準原子 × 該處 CN/multiplicity」所隱含的期望值格點的最近距離,以抽樣標準差為單位
基礎floor_ratioALT 分子數 ÷ 該連鎖區段錯誤地板的期望分子數
基礎puritycn_localphi_stratum樣本與區段層的脈絡量
單倍型hp_concentrationALT 分子集中在單一 germline 單倍型的比例
單倍型hp_depth_ratiohp_untagged_frac該單倍型的深度脈絡與未標記比例
局部系譜lineage_coherenceALT 分子在 Lu(m) 各節點上的集中度:最大節點佔比
局部系譜llr_nested對數似然比:「ALT 分子構成某一節點的巢狀子群」對「ALT 分子依 νu 隨機散布」
局部系譜four_gamete_viol把候選加進 Xu 後是否產生四配子違反,及其超過地板的程度
局部系譜expected_vaf_local最佳配對節點的 ν 換算出的期望 VAF,與觀測值的差
局部系譜n_spank_ustate_entropy證據量與局部複雜度的脈絡量
全部tier0_oktier1_oktier2_ok該特徵群是否可用的指示通道
缺席不是零

沒有連鎖進任何連鎖視窗的候選就沒有局部系譜特徵群,而那個比例逐樣本由 1.9% 到 59.5% (連鎖率的補數)。把缺席欄位填 0 會讓模型學到「0 代表沒有支持」, 而那正好與「這個位點根本沒有可用的連鎖區段」相反。 每一特徵群都必須有獨立的可用性指示通道,訓練資料必須同時包含可用與不可用的兩種樣本, 而且驗收時三組各報一次 precision/recall,並另報不可用子集。 分層報告在這裡不是謹慎,是必要條件:可用比例本身逐樣本就差兩倍以上, 把 COLO829 那種 98% 與 HCC1395 那種 41% 平均起來, 一邊的增益與另一邊的傷害會同時消失。

三種整合方式

方式做法成本風險
A · 事後重新評分以 caller 分數 + 特徵訓練一個小模型,只改排序低,不需重訓 caller受限於 caller 已經丟掉的候選
B · 加成輸入通道把特徵接進 caller 的張量後重訓高,需要完整訓練資料與算力上限最高,但難以歸因是哪一組特徵有用
C · 前置過濾用特徵先砍候選不建議:不可逆,且會安靜地傷害 recall
三種整合方式各自接在 somatic caller 流程的哪一個位置 上排是 somatic caller 本來的流程:產生候選、算特徵與張量、模型給分、輸出 VCF。 三種整合方式的差別就是接在這條流程的哪一個位置,而位置決定了成本與風險。 方式 A 接在模型給分之後,只重新排序,所以不必重訓 caller,成本最低, 但它只碰得到 caller 已經留下來的候選,被前面丟掉的救不回來。 方式 B 把特徵接進張量之後重訓,位置最靠前,上限最高, 代價是需要完整訓練資料與算力,而且很難歸因是哪一組特徵有用。 方式 C 接在產生候選之前,直接用特徵先砍候選。 它便宜,但它是不可逆的:被砍掉的候選不會出現在任何後續指標裡, 所以它會安靜地傷害 recall,圖上因此標為不建議。 讀法是:建議順序是 A 先做,因為 A 可歸因、每一組特徵可以單獨加減; 而如果 A 沒有效果,B 幾乎不可能有 —— A 的失敗代表這些特徵 在 caller 已有的資訊之外沒有增量。 三種接法接在 caller 流程的哪裡 位置決定成本與風險:越靠前上限越高,也越不可逆 產生候選 特徵/張量 模型給分 輸出 VCF C B A A · 事後重新評分 —— 接在給分之後,只改排序。成本低、不必重訓;但只碰得到 caller 已留下的候選。 B · 加成輸入通道 —— 特徵接進張量後重訓。上限最高;需要完整訓練資料與算力,且難以歸因哪一組特徵有用。 C · 前置過濾(不建議) —— 先砍候選。便宜,但不可逆:被砍掉的不出現在任何後續指標裡,會安靜地傷害 recall。 建議順序 A 先做 —— A 沒有效果,B 幾乎不可能有。
上排是 caller 本來的流程,三個箭頭是三種接法的接入位置 —— 位置決定成本與風險。越靠前上限越高,也越不可逆: A 接在給分之後,只改排序,碰不到已經被丟掉的候選; B 接進張量後重訓,上限最高但難以歸因; C 接在產生候選之前,被它砍掉的候選不會出現在任何後續指標裡, 所以它傷害 recall 的方式是安靜的。

建議順序是 A 先做。它便宜、可歸因(每一組特徵可以單獨加減)、 而且如果 A 沒有效果,B 幾乎不可能有 —— 因為 A 的失敗代表這些特徵 在 caller 已有的資訊之外沒有增量。

不論走哪一種接法,訓練與評估的切分有四條不能違反的規則, 其中最容易犯的是用逐變異的隨機切分 —— 同一個連鎖區段裡的候選共用一份局部系譜,隨機切分會讓指標系統性偏高。

切分紀律:四條規則逐條(切分單位、錨點來源、同一條管線、分層驗收)
  • 不可以用逐變異的隨機切分。同一個連鎖區段裡的候選共用一份局部系譜, 彼此高度相關;隨機切分會讓同一個連鎖區段同時出現在訓練與驗證集, 指標會系統性偏高。切分必須以樣本為單位,連鎖視窗層另做一次巢狀檢查。
  • 錨點必須由 caller 自己產生,不得取自 truth set 或 的標註。
  • 訓練與推論時的特徵必須由同一條管線算出;全域基準的實作、 λ、地板估計方式都要寫進 provenance。
  • 驗收須跨純度、突變負荷、平台與 basecaller 分層報告, 以檢查

驗收設計

核心:局部打散的負對照

連鎖區段層級的結論沒有金標準,但有一個乾淨的負對照。 在每個連鎖區段內部打散分子的 somatic–somatic 配對,同時保留: 每個位點的 somatic 邊際計數、每條分子的覆蓋遮罩、 每個位點的 germline–somatic 列聯表、深度與品質分層。 打散之後,資料在所有邊際統計量上與原始資料相同, 唯一被移除的就是共現

比較量什麼結論的讀法
原始 vs 打散REFINES_GLOBAL 的個數差值即本規格的主要效果量;打散提供虛無分布
原始 vs 打散留一分子的預測似然共現是否真的攜帶邊際以外的資訊
逐連鎖區段可置換分子比例、實際改變比例、獨立打散份數三者任一過低時輸出 SHUFFLE_UNDERPOWERED,不得把「沒有差別」解讀為「沒有共現訊號」

其餘驗收項目

其餘十項分成三段:機率核、相容性、地板檢查程式有沒有寫對; 先驗影響、可辨識性、重現性、半合成檢查結論是不是由資料而非旋鈕決定; 下游、留一法、切分檢查外部驗證有沒有洩漏。 每一項都要有「失敗代表什麼」,否則它只是一個會綠的測試。

驗收全表:十個層級各自必須通過的測試與失敗意義(實作時逐列勾)
層級必須通過的測試失敗意義
機率核每個遮罩下的樣式機率加總為 1;部分覆蓋 read 的邊際化結果與顯式加總一致似然未正規化或把未覆蓋位點當成參考型
相容性對隨機狀態集合,四配子檢定的結果與 perfect phylogeny 建構的成敗一致相容性判定與樹建構不同源
地板以 germline-only 連鎖區段(無 somatic 訊號)估到的假狀態率與模型預測一致地板訂錯 —— 訂低會生出假 lineage,訂高會吃掉真的低頻支
先驗影響λ 掃描下 REFINES_GLOBAL 的個數變化須平滑且可解釋結論由先驗而非資料決定
可辨識性刻意合成兩支 CCF 相同的 lineage:僅用邊際的擬合必須無法分開,加入共現後必須分開程式以隱藏的先驗或 bug 假造了可辨識性
重現性同細胞株跨 basecaller、跨機構的連鎖區段分類一致度;與既有實作的 0.909 對照分類是流程雜訊而非訊號
半合成兩株已知細胞株依已知比例混合,在兩者基因型相異的連鎖視窗回收兩支結構局部推論在已知答案上就不對
下游三組特徵各自的 precision/recall 增益,含不可用子集特徵沒有增量,或在多數子集上造成傷害
留一法候選出現在自己的錨點集合時必須直接失敗特徵抄了候選自己的 ALT 計數
切分S2 與全域基準輸入位點的交集為空;訓練/驗證不共用樣本重複計數或資料洩漏
可分析的頻率連鎖視窗:傳統 VAF 擬合區間與譜系內比例區間的比較 兩條橫軸畫的是同一件事在兩條路線上的可用範圍,範圍都以「佔細胞的比例」為單位。 上排是傳統做法:原始研究刻意只在 VAF 零點一二到零點二四之間擬合, 更低的一端混著假陽性,更高的一端已經進入 clonal 峰。 把它換算成細胞比例,在純度零點六的樣本上大約是零點四到零點八,是一段相當窄的區間。 下排是本頁的做法:下界由每一個連鎖區段自己的深度與地板總量決定,不是一個固定常數, 在每個家族六十條分子、地板總量千分之五時大約是零點零五; 上界則由譜系內的 clonal 堆積決定,暫定零點五。 右側的表列出下界怎麼隨分子數與地板總量變化。 底下標出一個重要提醒:上界零點五是全文最弱的一個數字,必須由模擬重新推定。 讀法是:把比例改成在譜系內部量,動態範圍不但沒有變窄,反而變寬了。 可分析的比例連鎖視窗 傳統 VAF 擬合區間 0 0.5 1.0 換算成細胞比例約 0.4–0.8 HFS 譜系內比例區間 0 0.5 1.0 約 0.05–0.5,且下界逐連鎖區段計算 下界隨深度變化 分子數 δ=0.005 δ=0.01 20 0.10 0.12 30 0.08 0.10 60 0.05 0.07 100 0.04 0.05 下界不是常數,要逐連鎖區段算 把比例改成在譜系內部量,動態範圍不但沒有變窄,反而變寬。 但上界 0.5 是本頁最弱的一個數字:它是從傳統區間的比例類推來的, 必須由模擬重新推定,不要當成已知常數引用。
驗收項目裡「地板」那一項要驗的就是這條下界。 圖上兩條軸是兩種比例的可分析區間:上面是傳統做法在細胞比例上量得動的範圍, 下面是本規格在連鎖區段內分子比例上量得動的範圍(圖上寫「譜系內比例」,就是 νu)。 右下那張小表是下界怎麼算出來的:本頁的地板總量 δu 取 1%, 每個家族 60 條分子時下界約 0.07;分子數加倍就往下掉一截,地板減半也會往下掉一截。 注意下界與地板是兩個量δu輸入, 下界則是它在這個分子數之下換算出來的「一支要佔多少才看得見」。 所以下界必須逐連鎖區段算,不能當成固定常數引用。
實作順序(展開查看)
步驟做什麼可單獨驗證的內容不得宣稱之事
前處理稀疏狀態表(含遮罩)與就地錯誤率校準能重現現行分類的計數;地板與 germline-only 連鎖區段一致任何生物學結論
單區段模型單一連鎖區段的局部後驗,先驗用平坦分布機率核、相容性、邊際化三項測試解析度優於全域
全域先驗S1 建立原子先驗與逃逸質量λ 敏感度;合成的等值 CCF 例子可分開真實樣本的準確度
全樣本目錄輸出目錄與 REFINES_GLOBAL 計數打散負對照的效果量與虛無分布細胞層級的群數
Caller 基準特徵基礎與單倍型特徵 + 事後重新評分是否勝過 caller 既有的相位通道局部系譜特徵的貢獻
Caller 局部特徵加入局部系譜特徵並分層報告可用子集上的增益、不可用子集上的無害總體平均掩蓋分層差異
外部重現重現性與癌別關聯;匯出成對關係跨 basecaller 一致度沒有金標準的準確度宣稱

既有資料與文獻邊界

既有實作的資料規模

以下數值來自本實驗室既有實作在七個資料集、六個生物樣本、chr1–22 上的輸出。 在本規格的脈絡下,它們的意義與前一篇不同:不再是「有幾條約束」, 而是「有幾個連鎖區段可以各自輸出一列」。

已經量到的全基因體譜:七個資料集的連鎖視窗拓撲類別組成 七條堆疊長條,每一條是一個資料集,橫軸是佔全部可分析單位的百分比。 每條分成三段:單層無分支、多層無分支,以及其餘(分支類與未解析加起來)。 單層無分支的範圍是一成九點五到四成六點八, 多層無分支的範圍是三成八點一到五成三點三,在七個資料集裡有六個是最大的一類。 右側標出每個資料集的可分析單位數,從四千二百多到兩萬三千多不等。 下方列出重現性:同一個細胞株在不同機構、不同 basecaller 下的組成相似度是零點九零九, 同為乳癌的樣本之間是零點八四到零點九零, 兩個肺癌樣本之間是零點八六,跨癌別則掉到零點五九到零點七八。 一個換了 basecaller 還能維持零點九零九、卻分得開癌別的統計量,帶著真的訊號。 最下方提醒:這張圖尚未針對覆蓋幾何做校正,所以類別比例還不能直接當成演化結論。 讀法是:這已經是一個可重現的全基因體單倍型頻譜,重新設計要建在它上面。 已經量到的全基因體譜 七個資料集、六個生物樣本、chr1–22,每個連鎖區段是「連鎖視窗 × 單倍型家族」。 HCC1395_HKU 34.5% 52.2% 9,130 HCC1395_NYGC 36.0% 52.3% 5,308 HCC1937 28.8% 53.3% 4,245 HCC1954 46.8% 38.1% 5,647 H1437 29.8% 53.2% 13,740 H2009 19.5% 52.2% 23,128 COLO829 42.8% 45.7% 10,757 單層無分支 多層無分支(七個裡有六個是最大類) 分支類 + 未解析 重現性:同細胞株換機構換 basecaller 0.909;乳癌之間 0.84–0.90; 兩個肺癌之間 0.86;跨癌別掉到 0.59–0.78。換了 basecaller 還守得住,代表訊號是真的。 但這張圖還沒針對覆蓋幾何校正,所以類別比例還不能直接當成演化結論。
七個資料集的連鎖視窗分類組成。 在本規格下這張圖量的是局部系譜的形狀分布:印出來的兩段是單層無分支與多層無分支, 第三段(分支類與未解析合計)是剩下來的 11.5%–28.3%,圖上沒有標數字。 把未解析拆掉之後,真正的分支類是下表那個 6%–22%。 而分岔正是全域頻率譜最無能為力的形狀 —— 也就是說,這條路線最想要的那一類, 在既有的分類裡是最薄的一層。
觀測數值在本規格的意義
切出來的連鎖區段總數98,955 個單倍型連鎖區段分母之一:所有比例的最外層
 其中帶至少一個變異等位85,941 個分母之二:棄權率該用的那一個
  候選集合完整列舉75,224 個(87.53%候選狀態列舉真的跑得完的比例
  撞到搜尋節點上限而棄權10,717 個(12.47%標記並留在分母,不是丟掉
  可排序的連鎖區段71,955 個再往下一層的子集,非前者的替代數字
分支類的比例各資料集 6%–22%最可能掛上 REFINES_GLOBAL 的一類
多層無分支的比例38.1%–53.3%先後關係的來源,數量最多
更換 basecaller 的重現性同細胞株跨機構 0.909驗收設計中重現性測試的既有對照值
單一連鎖區段的覆蓋結構194 條 read,覆蓋 3 個位點者僅 3 條遮罩必須留下來的直接依據

三點必須先講清楚。其一,前四個數字是巢狀的子集,不是四個版本的同一件事 —— 所以每一個百分比都必須帶著自己的分母出現。

算法補充:棄權率該用哪一個分母(12.47% 與 10.83% 回答的是不同問題)

98,955 是全部連鎖區段;85,941 是其中帶變異的;75,224 是其中列舉得完的; 71,955 是再往下可排序的。棄權率該用的分母是 85,941(12.47%), 因為不帶變異的連鎖區段根本沒有東西可列舉、也撞不到任何上限; 用 98,955 算會得到 10.83%,那個數字回答的是另一個問題。

其二,形狀分布是一條啟發式規則作用「之後」的結果。 既有實作用 read-AF 分數排序並列的最小成本解,但它並不只是排序 —— 它只保留分數最高者、丟掉其餘,因而在任何結果被指派之前就縮減了候選集合。 七個資料集合計丟掉 972,592 個最小成本頂點集合中的 292,065 個, 也就是 30.03%,並介入了 71,955 個可排序連鎖區段中的 16,830 個。 所以「解到單一拓撲」的比例不是簡約性本身給的: 單就簡約性可歸因的區間是 64.89% 到 88.26%。 這正是推論演算法要把那個分數換成似然的具體理由 —— 問題不在它排得對不對,在它先砍再排

其三,0.909 是重現性,不是準確度,而且它比較的是 兩份已經過上述縮減的組成。它證明分類不是流程雜訊,不證明分類是對的; 改成本規格的推論之後應重算此值並比較。

分岔比例的檢定力限制

分岔是 REFINES_GLOBAL 最可能的來源,也是最稀缺的一類。 既有實作在七個資料集上的形狀組成,無分支類一律佔 71.7%–88.5% (最高的是 COLO829 的 88.5%,最低的是 H2009 的 71.7%)。

逐資料集:七個資料集的無分支比例與已分類連鎖區段數
資料集無分支合計已分類連鎖區段數
COLO82988.5%10,757
HCC1395_NYGC88.3%5,308
HCC1395_HKU86.7%9,130
HCC195484.9%5,647
H143783.0%13,740
HCC193782.1%4,245
H200971.7%23,128

關鍵是為什麼稀少,而這個理由是機械性的: 一條鏈只要一個觀測到的樣式就建得起來,一個分岔要兩個互不包含的樣式。 所以分岔的稀缺,量到的至少有一部分是「這個連鎖區段帶了多少不同的樣式證據」, 而不是「底下的結構有沒有分岔」 —— 既有實作沒有把這兩者分開,本規格必須分開。

分法是現成的:錯誤地板決定了「第二個樣式若真的存在,看得到它的機率」。 所以每個判為無分支的連鎖區段都要附一個檢定力數字 —— 在該連鎖區段的分子數與地板之下,一支佔 5% 的旁支被看到的機率是多少。 檢定力低的無分支連鎖區段應標 BRANCH_UNDERPOWERED, 而不得計入「無分支佔多少」這個比例的分子。 沒有這一步,形狀組成量到的是覆蓋深度,不是演化。

還有一個落差常被接在這一段後面,但它的來源不是檢定力: 單一拓撲的比例由 35.3%(H2009)到 90.6%(HCC1954),差 55 個百分點。 那個排序與 k 的分布同向 —— k 大的資料集並列的候選多、 解到單一拓撲的少,而不是它們的腫瘤比較單純。 也就是說,單一拓撲比例量的是「有沒有第二個樣式」,k 的分布量的是「有幾個位點」, 兩者都不是演化訊號,而且要分開報。

合起來的結論是:庫存那張表裡 6%–22% 的分支類是上限,不是估計值。 判定分岔需要兩個單一狀態皆存在,雙狀態的缺席不是抽樣零也不是錯誤地板 —— 後者正是局部生成模型中錯誤地板要檢驗的對象。

相關工作與新意邊界

既有能力代表方法本規格不可據此宣稱之處
由邊際 VAF 與拷貝數重建細胞層級 clone treePyClone、SciClone、PhyloWGS、Pairtree、DPClust不可宣稱全域樹準確度提升
把中性演化的尾巴與真正的 subclone 分開MOBSTER不可宣稱首次處理「頻率群集不等於 clone」
成對位點的分子狀態(含相位)建 subclone 與樹PairClone、TreeClone不可宣稱把共現寫進似然為新意,也不可宣稱全域方法只吃邊際 —— 差別在可用位點對的數量級
相位輔助的等位/subclonal 拷貝數Battenberg、Refphase、HATCHet2不可宣稱 phased CN 為新意
somatic caller 使用 germline 相位通道ClairS、DeepSomatic不可宣稱單倍型感知呼叫為新意 —— 單倍型特徵群是基準不是貢獻

候選新意因此只剩兩處,而且都必須先做文獻檢索才能使用「首次」二字: 其一,把可變 k 的多位點共現當成估計目標本身(而不是當成全域樹的約束), 並以逐連鎖區段的旗標與打散負對照給出可計數的效果量; 其二,把同一條 germline 單倍型內部的 somatic lineage 一致性做成 caller 的特徵通道。 PairClone 與 TreeClone 處理的是成對位點且目標仍是全域 subclone 結構; ClairS 的相位通道停在兩條 germline 單倍型這一層。

預期結果與限制

結果必須同時通過兩個獨立終點: REFINES_GLOBAL 的數量應明顯高於局部打散所產生的虛無分布; caller 特徵則必須在留一法下,對既有相位通道提供額外的 precision/recall 增益。 這兩項都不能證明局部節點就是真實細胞群,也不能證明全域 clone tree 因此變準。

預期結果模式

資料情境預期輸出正確結論
單一 clone、突變負荷低幾乎全部 GLOBAL_CONSISTENT,打散後無差異局部沒有可加的資訊,這是正確的負結果
兩支 CCF 相同的 lineage,且其變異落在同一批連鎖區段該批連鎖區段掛 REFINES_GLOBAL,打散後消失本規格的主要成立情境
兩支 CCF 相同,但變異散在不同連鎖區段無旗標,與單一 clone 無法區分read 跨距的硬上限,不是模型失敗
家族誤標率偏高INCOMPATIBLE 比例上升,REFINES_GLOBAL 同步上升,但打散後不下降是地板問題,不得解讀為 subclone
kataegis 叢單一連鎖視窗貢獻大量高度相關的連鎖區段整叢降權為一個有效觀測
下游:局部系譜特徵群可用的子集precision 上升而 recall 持平或微升可歸因於 lineage 一致性
下游:局部系譜特徵群不可用的子集指標不變缺席通道有效;若下降則缺席契約寫錯了

尚未解決的問題

十個問題,每一條都是「目前沒有答案」,不是「還沒寫」。 標題那一句就是問題本身;展開看的是它為什麼還沒被回答,以及要回答它得先做什麼。 其中連鎖率為什麼逐樣本差兩倍以上是第一個該回答的 —— 在弄清楚它之前,無法事先判斷一個新樣本值不值得跑。

整條路線尚未經實測。

本頁各項皆有結構性的理由,但都未在真實資料上檢驗。 最小的驗證是驗收設計中等值 CCF 的合成例子:確認僅用邊際的擬合必然分不開, 再量需要多少條跨越分子、地板要估得多準才分得開而不誤拆。

λ 沒有客觀的訂法。

它是這條路線唯一一個靠判斷的旋鈕。 可以掃描與做敏感度分析,但「掃出來最好看的那個值」不是估計。 是否有辦法由資料本身(例如以留一連鎖區段的預測似然)把它定下來,尚未回答。

全域基準的錯誤會整批傳下去。

切斷回饋的代價是局部無法修正全域; 若全域基準的群數少算一個,所有原子都會偏,而局部只會用逃逸質量吸收, 表現為 LOCAL_ONLY_LINEAGE 大量出現。這個失效模式的診斷還沒設計。

連鎖率為什麼逐樣本差兩倍以上,沒有人知道。

它由 40.5% 到 98.1%, 既有實作沒有把它與突變數、突變密度、深度或讀長分布關聯起來, 而技術重複只解釋得了 5.2 個百分點。 在弄清楚機制之前,無法事先判斷一個新樣本值不值得跑, 也無法判斷低連鎖率是生物學(突變稀疏)還是流程(定相破碎、讀長短)造成的 —— 兩者的補救方式完全不同。這是本規格第一個該回答的經驗問題。

連鎖區段的來源有偏性。

它們只來自變異密集的連鎖視窗, 而密集有時並非偶然(kataegis)。降權規則已寫,但降權的強度沒有依據。

甲基化是一個誘人但被綁住的第四層。

既有實作量到, 在有穩定甲基化差異的位點中,51.2%–80.4% 的差異落在單一單倍型的 ALT read 內部 —— 那正是本規格想分開的地方,所以它看起來是局部系譜特徵群之外的另一條通道。 但等位特異甲基化以順式作用,因此它與界定連鎖區段的單倍型切分並非獨立: 一個甲基化分群與 ALT/REF 軸或單倍型軸對齊,機制上本來就會發生, 不需要任何 lineage 解釋。要把它變成證據,得先有一個能把「順式機制」與 「lineage 差異」分開的設計,而那個設計目前不存在。

LOH 區域怎麼辦,只有一個保守答案。

局部生成模型把非中性區段標記並排除, 但這在 CN 變異廣泛的實體腫瘤可能排掉一半以上的連鎖區段。 要放寬就得讓模型表達突變的喪失,而那會破壞四配子相容性 與「給定狀態集合樹即唯一」這兩個讓局部推論便宜的性質。這個取捨還沒有評估。

與整合篇的關係尚未決定。

本規格刻意切斷回饋,整合篇則走完整聯合。 兩者不是替代關係,但如果聯合模型可行,本規格的局部後驗應該是它的一個近似 —— 這個近似有多差,沒有評估。

「首次」尚不可用,而且門檻比原本以為的高。

PairClone 與 TreeClone 已經把成對共現與相位寫進全域似然,所以候選新意不能再宣稱在「使用共現」這件事上, 只能落在「可用位點對多兩個數量級、k 可變、逐連鎖區段交付」這幾點。 Pairtree、CliP、LICHeE 與任何以 read linkage 輔助 subclone 的後續工作仍待系統檢索。

短讀成對共現的實際產量沒有實測過。

開頭那個 0.15%–1.5% 是 Poisson 估計, 不是量出來的;PairClone/TreeClone 論文用的資料集實際有幾對可用, 以及它們在那個產量下對全域結論的改變幅度有多大,會直接決定 「長讀多兩個數量級」這句話值多少。這是說服審稿人最需要的一個數字,目前沒有。

原始文獻與程式碼

本頁「輸入與前處理契約」、「推論演算法與計算預算」的上限與成本數字, 以及「既有實作的資料規模」的全部計數, 均出自本實驗室既有實作的紀錄:廖子游, Subclonal reconstruction using somatic haplotagging and methylation profiles with Nanopore sequencing,國立中正大學資訊工程學系碩士論文,2026 年 7 月。 該實作與本規格的目標一致(局部、分子層級、不重建全域樹), 但推論方式不同:它以 Camin–Sokal 簡約的最小潛在頂點搜尋產生候選, 再以 read-AF 分數縮減候選集合;本規格把後者換成邊際似然,並補上錯誤地板與全域先驗。 凡本頁引用其數值之處,均為該實作實測到的量,不是本規格已驗證的結果。

以 read 上的多位點狀態計數直接建模 subclone,最接近的既有做法是 Zhou T, Sengupta S, Müller P, Ji Y, PairClone: a Bayesian subclone caller based on mutation pairsJ 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 pairsarXiv:1703.03853

延伸引用:全域重建路線、下游 caller、切斷回饋與四配子檢定的原始出處

作為對照基準的全域重建路線:SciClone,Miller CA et al., PLoS Comput Biol 2014;10:e1003665;PyClone,Roth A et al., Nat Methods 2014;11:396–398; PhyloWGS,Deshwar AG et al., Genome Biol 2015;16:35 (doi 10.1186/s13059-015-0602-8); Pairtree,Wintersinger JA et al., Blood Cancer Discov 2022;3:208–219 (doi 10.1158/2643-3230.BCD-21-0092), 程式碼在 github.com/morrislab/pairtree; MOBSTER,Caravagna G et al., Nat Genet 2020;52:898–907 (doi 10.1038/s41588-020-0675-5)。 評比數據見 Salcedo A et al., Nat Biotechnol 2024 (doi 10.1038/s41587-024-02250-y)。

下游 caller:ClairS,Zheng Z et al., ClairS: a deep-learning method for long-read somatic small variant callingbioRxiv 2023,doi 10.1101/2023.08.17.553778, 程式碼在 github.com/HKU-BAL/ClairS; DeepSomatic,Park J, Cook DE, Chang PC et al., Accurate somatic small variant discovery for multiple sequencing technologies with DeepSomaticNat Biotechnol 2025 (doi 10.1038/s41587-025-02839-x), 程式碼在 github.com/google/deepsomatic

「切斷回饋」這個做法的統計出處為 Plummer M, Cuts in Bayesian graphical models, Stat Comput 2015;25:37–43 (doi 10.1007/s11222-014-9503-z); 該文同時指出這類後驗不是任何聯合模型的邊際分布,也指出其計算上的陷阱, 是「全域頻率基準如何成為獨立先驗」中「不得寫成聯合估計」的依據。 四配子檢定與 perfect phylogeny 的標準參考為 Gusfield D, Efficient algorithms for inferring evolutionary trees, Networks 1991;21:19–28; 簡約假設的原始出處為 Camin JH, Sokal RR, A method for deducing branching sequences in phylogeny, Evolution 1965;19:311–326。

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

本模組術語

Clair3
以深度學習做 germline variant calling 的工具,長 read 上常用。
ClairS
配對 tumor–normal 的 somatic variant caller。tumor-only 版本叫 ClairS-TO。
DeepSomatic
Google 的 somatic variant caller,同樣有 tumor-only 版本。
DeepVariant
Google 開發的 germline variant caller,把 pileup 轉成影像再用 CNN 分類。
LOH(異型合子性喪失)
原本 heterozygous 的區域變成只剩一種 allele。LOH 不等於缺失 —— 也可能是一條 haplotype 遺失後另一條被複製(copy-neutral LOH)。
LongPhase
實驗室開發的長 read phasing 工具,是後續所有工具的共同基礎。輸出 germline haplotype。
LongPhase-S
配對 tumor–normal 版本。做腫瘤 DNA 比例估計(輸出欄位名為 purity)、比例感知的 somatic variant 過濾,以及 somatic haplotagging。
allele-specific copy number(等位特異拷貝數)
把一個區段的總拷貝數拆成兩個親源等位各自的份數,通常記為 major 與 minor。總拷貝數相同而兩個等位不同的狀態(例如 2+0 與 1+1)在深度上完全一致,只有 B-allele frequency 分得開 —— copy-neutral LOH 即為此類。
cancer cell fraction(癌細胞比例)
帶有某個特定突變的腫瘤細胞佔全部腫瘤細胞的比例。用來區分 clonal(1)與 subclonal(<1)突變。不等於 VAF。
clone tree(克隆演化樹)
描述腫瘤內各群細胞祖先關係的樹:節點是一群帶有相同變異組合的細胞,邊代表在祖先之上又多拿到變異。要注意同一組群集常常有多棵樹同時相容。
copy number(拷貝數)
某段基因體在細胞內的拷貝數。多數正常常染色體區段為 2,可再分為 major 與 minor allele copy number。
data leakage(資料洩漏)
訓練或模型選擇階段取得了評估資料的資訊,讓效能被高估。實驗室的做法是按染色體切分,確保同一個位點不會同時出現在訓練、驗證與最終測試中。
distribution shift(分布偏移)
測試資料的組成與訓練資料存在顯著差異(例如正負樣本比例不同),使得模型表現不如預期。
haplotagging
根據已 phase 好的 variants,把每一條 read 指派到 HP1 或 HP2,並把結果寫回 BAM 的 HP tag。
high-confidence region
benchmark 建立者指定的高可信度區間,通常以 BED 檔表示;區間外不宜視為具有相同標註可靠度。
latent node(潛在節點)
建樹時為了讓圖連得起來而補進的中間狀態,沒有被任何 read 直接觀測到。它是模型的產物,不能當成「還沒觀察到的細胞」。
mutation frequency spectrum(突變頻率譜)
把一份樣本裡所有 somatic 變異的 VAF 畫成直方圖後得到的分布。分布上的峰與肩對應不同大小的細胞群,最低頻端的尾巴斜率則被用來判斷有沒有天擇。
phase block
一段可建立連續相位關係的區域。read 長度不足、缺少 informative heterozygous 位點或證據不一致時,可能形成不同 phase blocks。
pigeonhole(鴿籠原理)
由父代與子代的細胞比例限制樹形的算術規則:一個細胞至多屬於父節點底下的一個子節點,所以各子節點的細胞比例加起來不得超過父節點。名稱來自鴿籠原理 —— 東西放進籠子,總量不會憑空變多。在 subclone 重建的文獻中也稱為 sum rule 或 crossing rule。它只能排除樹,不能挑出樹:通過檢查的候選通常仍不只一棵,其餘要靠 parsimony 之類的偏好決定。
單倍型家族(一條 germline 單倍型,加上由它衍生的 somatic 單倍型)
把 read 依 HP tag 分成的兩組之一。家族一HP1 與從它長出來的 HP1-1家族二HP2HP2-1;歸不到任何一條 germline 單倍型的 HP3 不屬於任何一族。叫「家族」是因為它把一條 germline 單倍型與由它衍生的 somatic 單倍型收在同一組裡 —— 分組看的是 germline 那一層,不是有沒有帶 somatic 突變。

在同一個 phase block 內,一個家族對應一條染色體拷貝,所以「兩族」就是那個位置上的兩條同源染色體。但兩件事不成立:其一,軟體判定不出哪一族來自父親、哪一族來自母親(那需要另外定序父母);其二,標號只在該 phase block 內有定義,跨 block 的「家族一」並非同一條染色體。