研究指引 · 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 的三組特徵、留一法、缺席契約與分層驗收
- 列出驗收條件、實作順序,以及這條路線不可宣稱的事
為什麼重要
換掉的是目標,不是方法
前一篇把兩類觀測寫成兩個相乘的因子,接上同一組全域參數 ,目標仍然是一棵全基因體、細胞層級的 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.完整條目 →。 這個目標有一個現實問題:以邊際 VAF 與 等位特異拷貝數allele-specific copy number 等位特異拷貝數把一個區段的總拷貝數拆成兩個親源等位各自的份數,通常記為 major 與 minor。總拷貝數相同而兩個等位不同的狀態(例如 2+0 與 1+1)在深度上完全一致,只有 B-allele frequency 分得開 —— copy-neutral LOH 即為此類。The total copy number of a segment split into the counts contributed by each parental allele, usually reported as major and minor. States with equal total but different split (2+0 versus 1+1) are identical in depth and separable only by B-allele frequency; copy-neutral LOH is exactly such a case.完整條目 → 為觀測的既有方法已經成熟且被系統性評比過, 要在同一個目標上勝過它們,需要的證據量遠大於多變異連鎖視窗能提供的。 中篇的計算已經指出,可用的多變異連鎖區段在低突變負荷樣本只有數百個。
本篇因此換一個估計目標:不重建全域的細胞層級系統發生, 而是重建每一個連鎖區段內部的分子系譜 —— 哪幾種局部狀態真的存在、 各佔多少分子、以及它們之間的先後關係。這個目標的解析度是單倍型連鎖區段 —— 一個連鎖視窗(read 連鎖起來的一段,不跨 phase setphase 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.完整條目 →)的一條染色體拷貝,不是「全基因體共用的一個 CCF 位置」。
| 既有做法的目標 | 本規格的目標 | |
|---|---|---|
| 估計對象 | 全基因體的細胞群數、比例與樹 | 每個連鎖區段內各局部狀態的比例與先後 |
| 觀測進入似然的形式 | 每個變異的邊際 VAF + 拷貝數 | 同一條分子上多個位點的共現 |
| 一個節點是什麼 | 一群細胞 | 一類分子狀態 |
| 解析度 | 一個全基因體共用的 CCF 位置 | 一個連鎖視窗、一條染色體拷貝 |
| 輸出幾棵樹 | 一棵 | 幾百到幾萬棵,彼此不相接 |
| 永遠分不出來的 | CCF 相同的兩支 lineage | 跨越 read 距離以外的任何關係 |
| 多分出一支要多少證據 | 約 20 個變異(兩群差 0.05 時);兩群越近需要越多 | 數十條跨越分子;與兩支靠得多近無關 |
兩條路線的差別集中在最後一列: 頻率譜要分開兩支,靠的是它們的比例不同;共現要分開兩支,靠的是它們的突變組合不同。 所以「兩支比例相同」對前者是致命的,對後者完全沒有影響 —— 只要它們帶的突變不一樣,分子上就看得出來。
局部這一側的限制在別的地方:它要求每一支各自夠多
(高於該連鎖區段的錯誤地板),而不是要求兩支之間差得夠開。
這個限制由後文「候選模型與輸出旗標」的 FLOOR_LIMITED 旗標承接。
交付物是「成對關係目錄」,不是幾百棵漂亮的樹;而本頁是規格,不是已驗證的方法
既有實作實測到,七個資料集中有六個,最大的一類是「一個連鎖區段只有兩個變異」 (完整分布見「輸入與前處理契約」)。所以這條路線的主要產出是全基因體的成對分子關係目錄 —— 這兩個突變共不共存、誰在前 —— 多層的局部系譜是其中的少數。 輸出格式、資料結構與驗收指標都要照這個比例設計,不要照「幾百棵樹」設計。
另外,本頁是規格,不是已驗證的方法(LiLT,工作名稱)。 文中所有數字都出自既有的原型實作,不是這份規格跑出來的結果 —— 自己實作時應當重新量一次,不要直接引用。
頻率譜為什麼看不到某些局部結構
「全域方法犧牲了局部解析度」如果只是印象,就沒有價值。 它其實是一個機制,而且可以換算成一個可直接與實際資料對照的數字: 要說服頻率譜「這裡是兩群不是一群」,需要幾個變異?
答案取決於兩群靠得多近,而且不是慢慢變多
| 兩群的期望 VAF | 頻率譜大約需要幾個變異才分得開 |
|---|---|
| 0.20 對 0.10 | 約 4 個 |
| 0.20 對 0.15 | 約 20 個 |
| 0.20 對 0.18 | 約 120 個 |
| 0.20 對 0.20 | 永遠不夠 |
最後一列不是「很難」,是沒有答案。 兩群的期望 VAF 相同時,把任何一個變異從這一群改指派到另一群, 直方圖完全不變 —— 兩群本來就落在同一根柱子裡。 加深度、加樣本、換更好的演算法都不會動它。 中篇的 例子即為此情形。
而一支只在一個連鎖視窗裡不同的 lineage,本來就帶不起那麼多變異
這是問題的另一半。一支細胞群若只在某個連鎖視窗內與它的母群體不同, 它與母群體的全部差別就是那個連鎖視窗裡的那幾個突變 —— 三個、五個。 它因此落在上表「開不出來」那一側,而且通常離得很遠。 所以這不是「不容易偵測」,是它在那個模型裡開不出來。
同一批 read,換一種看法
把那三個變異換成「它們在不在同一條分子上」,需要的證據量完全不同。 假設這三個突變其實不共存,那麼在 20 條跨越它們的分子裡, 同時帶三個突變的分子預期只有錯誤地板那麼多 —— 1% × 20 條 ≈ 0.2 條。
而實際看到 5 條。
這裡不需要門檻,也不需要任何統計檢定:預期 0.2、實際 5,結論就出來了。 而兩條路線看的是同一批 read —— 差別不在資料量,在把它們摘要成什麼。
上面那幾個數字是怎麼算出來的
「需要幾個變異」=多開一群要付的代價 ÷ 每個變異能提供的證據, 兩者都是對數似然,單位是 nat(自然對數底下的資訊量)。 本文不用 nat 當敘述單位,因為換算成「幾個變異」之後, 即可直接與實際資料比較。
- 代價是模型複雜度罰則 : 參數兩個、有效觀測三千時約 8。
- 每個變異的證據約為 , 其中 是深度、 是兩群的期望 VAF。 深度 50、0.20 對 0.15 時約 0.43,故 ,表中記為「約 20」。
- 時分子為零、商為無限大 —— 那就是兩群期望 VAF 相同的情形。 表中那四列是同一條曲線上的四個點,不是四件不同的事。
- 共現那一側: 條對 條的對數似然比約 12, 遠高於分辨兩個局部結構所需。這裡的錯誤地板(總量 ,本頁取 1%) 來自後文「局部生成模型」的逐位點錯誤 與家族誤標, 所以 是算出來的而不是另一個自由參數,而且不是常數,要逐連鎖區段估。
一個必要的但書
上面整段的前提是「似然只吃邊際計數」,而有兩個既有方法不滿足這個前提: 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 只處理成對位點:在那個產量下, 大於 2 的視窗幾乎不存在,可變 沒有意義。
方法
方法概覽
方法先把互不重疊的 read 分成兩組:單變異連鎖視窗 用於估計全域頻率基準, 多變異連鎖區段 用於逐區段重建局部分子系譜。全域結果只作為局部推論的先驗, 局部結果不回饋全域,因此同一條 read 不會在兩個 likelihood 中重複計數。
| 步驟 | 資料 | 產出 | 在本方法中的地位 |
|---|---|---|---|
| 前處理 | 配對 BAM + 候選集合 | 定相、體細胞標記、連鎖區段與稀疏狀態表 | 沿用既有工具;輸入契約必須固定 |
| 全域基準 | 的單變異資料 | 純度、CCF 原子、原子不確定度與 overdispersion | 沿用既有頻率譜方法 |
| 局部推論 | 的單一連鎖區段 + 全域先驗 | 狀態集合、比例、先後關係與品質旗標 | 本規格的核心 |
估計目標:一個連鎖區段的三個量
每個連鎖區段只估計三個對象:局部狀態集合、各狀態的專屬分子比例, 以及狀態之間的先後關係。以下先定義估計單位與輸出,再引入似然與先驗。
估的範圍:一個連鎖區段
沿用前一篇的連鎖區段: 是一個單倍型連鎖區段:一個連鎖視窗 × 一條 單倍型家族單倍型家族 一條 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.完整條目 →,
其位點集合為 、位點數 。
這樣切有一個具體的理由,不是為了方便: 一個連鎖區段裡的分子全部來自同一條染色體拷貝,另一條屬於另一個連鎖區段。 所以在連鎖區段內部,「兩個變異出現在同一條分子上」才真的等於 「它們在同一條拷貝上先後發生」。 若不先依家族切開,位於兩條拷貝上的一對變異(trans)永遠不會共存, 其計數表會與「兩支互不包含的 subclone」長得一模一樣 —— 連鎖區段這個切法就是為了讓這兩件事落進不同的連鎖區段而不會相遇。
要估的三個量
令局部狀態為 , 即該連鎖區段那幾個位點上「帶或不帶變異」的樣式。 一個連鎖區段的估計目標就是下面三個量,沒有第四個:
| 要估的量 | 記號 | 白話 |
|---|---|---|
| 狀態集合 | 這個連鎖區段裡有哪幾種分子 | |
| 專屬比例 | ,總和為 1 | 每一種各佔多少 |
| 先後關係 | 上的有根樹 | 誰是誰的祖先 |
三者合起來記為 與 , 樹根固定為全參考狀態 。
變數補充: 是專屬比例,不是累積比例
一條分子只由它所屬那一個節點的狀態產生,不由祖先節點產生; 祖先節點的 因此是「還停留在該狀態的分子」的比例,可以很小。 換句話說,樹的形狀不對 施加任何單調性, 它只規定哪些狀態可以同時存在。
先後關係幾乎不用估:四配子相容性
先後關係看起來是三個量裡最難的,其實不是 —— 它不需要搜尋。 兩個位點稱為相容,若四種樣式 中至多出現三種。 只要 裡每一對位點都相容,樹就已經被 決定了, 不必列舉任何樹形。實作上要做的只有兩件事:逐對檢查相容性, 以及在檢查失敗時把它當成品質訊號而不是新發現。
為什麼「兩兩相容」就足以決定一棵唯一的樹(含 )
把每個位點 換成一個集合 : 裡帶有該變異的狀態所成的集合。 四種樣式少一種,說的正是 與 這兩個集合的關係:
| 缺哪一種樣式 | 兩個集合的關係 | 樹上的意思 |
|---|---|---|
| 沒有 | B 在 A 的後代裡 | |
| 沒有 | A 在 B 的後代裡 | |
| 沒有 | A、B 在互不包含的兩支上 | |
| 四種都有 | 兩集合互相交錯 | 樹上不存在這種關係 |
所以「相容」=「這兩個集合不交錯」,也就是要嘛一個包含另一個,要嘛不相交。
唯一性:一族兩兩不交錯的集合,在數學上就是一棵樹 —— 包含關係本身給出父子關係( 的父親是包含它的最小那一個), 不相交的兩個成為兄弟。沒有任何自由度留給搜尋,所以樹是 的函數,不是待估的參數。 剩下的自由度只有一種:一條邊上若有數個位點的集合完全相同, 它們同時發生、在樹上擠在同一條邊,彼此的先後資料沒有講 —— 這就是「本質上唯一」那個「本質上」的全部內容。
為什麼不必額外檢查:關鍵在於兩兩相容不只是必要條件, 它同時是充分條件(perfect phylogeny 定理,Gusfield 1991)。 也就是說不存在「每一對都相容、但三個一起就不相容」的反例 —— 所以三個位點不必檢查 3 組三元組, 個位點也只要檢查 對,成本是 而不是 。 這正是「先後關係幾乎不用估」的來源。
實作上因此只需維護一張逐對的相容性表; 它同時就是右欄那個品質訊號的計數來源。
這給出兩個貫穿後面各節的結論:
- 困難不在搜尋樹,在決定 —— 哪些狀態的比例真的高於錯誤地板。後續的生成模型、先驗與候選選擇都在處理這一件事。
- 四格全滿不是發現了一支新 lineage,是四配子檢定失敗, 代表這個連鎖區段的錯誤地板偏高或家族誤標偏多。相容性檢定因此同時是品質訊號。
解讀補充:未觀測狀態可能是現存或歷史狀態
可以含沒有被任何分子觀測到的狀態。 觀測到 與 而 沒看到時, 兩者若不共用一個帶第一個突變的祖先,第一個突變就得發生兩次。 這種狀態只有兩種身分 —— 現存但沒抽到(), 或已被後代取代的歷史狀態()。 前一篇還有第三種:純粹因為「一步一個突變」的表示法而被迫補進來的記帳節點; 本規格允許一條邊帶多個突變,那一種就隨表示法一起消失了。 latent nodelatent node 潛在節點建樹時為了讓圖連得起來而補進的中間狀態,沒有被任何 read 直接觀測到。它是模型的產物,不能當成「還沒觀察到的細胞」。An intermediate state added during tree construction to keep the graph connected, not directly observed in any read. It is a product of the model and must not be read as an unobserved cell population.完整條目 → 那條警告仍然適用,但適用範圍窄了一半。
什麼時候算「估完了」
三個核心估計量之外,細胞比例與全域群對應還需要外部資訊。 每項輸出各有成立條件;條件不足時輸出旗標,而不是猜一個值 —— 這條規則是後文輸出旗標的語義來源。
| 要回答的問題 | 估到什麼 | 需要什麼 | 條件不足時的輸出 |
|---|---|---|---|
| 有哪幾種、各佔多少 | 與 | 只需要該連鎖區段的 read 與就地估到的錯誤率 | FLOOR_LIMITED |
| 誰在誰之前 | 上的先後關係 | 兩個單一狀態皆可見,且雙狀態的有無高於地板 | ORDER_UNRESOLVED |
| 換算成細胞比例 | 各節點的 CCFcancer cell fraction 癌細胞比例帶有某個特定突變的腫瘤細胞佔全部腫瘤細胞的比例。用來區分 clonal()與 subclonal()突變。不等於 VAF。The fraction of tumour cells carrying a given mutation; distinguishes clonal from subclonal. Not the same as VAF.完整條目 → | 純度與該處 拷貝數copy 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.完整條目 →,由全域基準提供 | 只報分子比例,標 CELL_SCALE_UNAVAILABLE |
| 對得上哪一個全域群 | 局部節點與全域群的對應 | 節點內至少一個變異有可信的全域指派 | UNLINKED_TO_GLOBAL |
局部節點是一類分子狀態,不是一群細胞
上表的「換算成細胞比例」只是尺度轉換,它不會把節點變成細胞群。 一個節點只表示「這些分子在 上帶有同一組突變」。 兩群基因體其他地方不同、但在這個連鎖視窗內完全相同的細胞,會落進同一個節點。 本實驗室的甲基化觀測正是這件事的實例:單一個「單倍型 + 等位」狀態之下, 仍可再分出數群甲基化模式。因此節點數的正確讀法恆為「至少此數」。 反向也不成立:同一群細胞在不同連鎖區段裡對應到不同的局部狀態, 因為每個連鎖區段看的位點集合不同。
局部生成模型
這一節只回答一個問題:給定某個連鎖區段裡各種真實分子狀態的比例 ,資料中的一條 read 為什麼會呈現目前看到的 0/1 字母? 由真實分子池到觀測資料共有四步:
- 抽出一種真實狀態。一條分子先依 ,屬於某個完整狀態 。
- 套用這條 read 的覆蓋範圍。read 通常只蓋住位點集合 ,未覆蓋的位置保持未知。
- 通過觀測錯誤。逐位點錯誤可能改變一個字母,家族誤標則可能把整條分子分到錯的單倍型家族。
- 形成實際資料。通過前述步驟後的字母 ,才是狀態表裡真正計數的觀測。
先看一條只有部分覆蓋的 read
設三位點連鎖區段的狀態集合為
,三種狀態的專屬比例依序為
。現在有一條 read 只覆蓋第一與第三個位點,
並在這兩處讀到 1–1。在尚未加入錯誤時,這條 read 同時相容於
101 與 111,因此其觀測機率為
。一條 read 支持的是一組狀態,不是一個狀態 —— 這句話就是下面整節的全部內容。
符號表:、、 三者的分工(隨時可回來查)
| 記號 | 在此例中的值 | 意義 |
|---|---|---|
000、101 或 111 | 這條分子的完整真實狀態,未直接觀測 | |
| 這條 read 實際覆蓋的位點 | ||
0–0 或 1–1 | 把完整狀態投影到覆蓋位置後,原本應看見的字母 | |
例如 1–1 | 經過錯誤通道後,資料中真正記錄的字母 |
三者的關係是一條單向鏈: 經覆蓋遮罩投影成 , 再經錯誤通道變成資料裡的 。推論要走的是反方向。
先條件於 read 已進入正確的單倍型家族;此時 read-level 的骨架可寫成:
是單一家族內的觀測通道:給定完整真實狀態 , 它處理部分覆蓋與逐位點錯誤,並給出最後看到 的機率。 為本連鎖區段的逐位點錯誤參數。跨家族誤標會連結另一個連鎖區段,於第二個細節另行加入。
這條式子的讀法是:逐一考慮所有可能的真實狀態, 以其分子比例 加權,再乘上它經觀測通道變成 的機率。 以下各小節先拆開 ,再把跨家族誤標接回來。
第一個細節:未覆蓋位點要邊際化
若暫時忽略所有錯誤,觀測通道只剩下「 投影至 後是否等於 」這項檢查。上面的完整式子便化為:
是該 read 的覆蓋遮罩, 是它讀到的字母。 沒覆蓋到的位點不可以填成參考型 —— 那會把一條中立的 read 變成反對某個狀態的證據。
證據補充:部分覆蓋為何只能算一個約束
一條有 個未覆蓋位點的 read,與 個狀態相容, 那些狀態構成狀態空間的一個子立方體。它主張的是「真實狀態落在其中」, 不是 個觀測 —— 把兩者混為一談會把證據量灌大 倍。 在似然寫法下這是自動的(上式的加總就是那個子立方體),但報告層不自動: 既有實作有一個 的連鎖視窗,靠一個全跨樣式(3 條 read) 加上十一個部分樣式就「解出唯一解」。 這種連鎖區段的唯一性大部分來自簡約性與先驗,不是來自資料。 因此每個連鎖區段必須分開報全跨分子與部分分子各貢獻多少證據; 只有全跨樣式與根是跨候選不變的,被部分 read 見證的狀態不是。 順帶一個反直覺的推論:部分覆蓋越多,相容的候選越多 —— 覆蓋越差的連鎖區段越容易看起來「並列」,而不是越容易被排除。
整個連鎖區段共用同一個分子池
整個連鎖區段共用一組 ,每條 read 由它自己的遮罩投影出來; 不同遮罩不各配一個獨立的比例向量,因為它們取樣的是同一池分子。 在 PCR-free 長讀下,一條 read 就是一條分子, 所以連鎖區段內部的計數預設為 multinomial,額外離散度沒有來源; 只有殘差真的顯示過度離散時才加,而且要說得出它來自哪裡。
把這個區段的所有 read 記為 。暫時忽略跨家族誤標時,單一連鎖區段的基礎似然只是逐 read 相乘,而每一項都是前述的加權和:
不同 read 可以有不同的遮罩 ,但全部共用同一組 。 同一遮罩的 read 可彙總成 multinomial 計數;兩種寫法是同一個模型。
第二個細節:兩種錯誤通道的形狀不同
前式中的 不能只寫成一個無來源的錯誤率。 逐位點錯誤與整條分子的家族誤標會產生不同形狀的假狀態,因此必須分開校準。 前者留在 內;後者要把成對的兩個連鎖區段一起寫:
令 為目前的單倍型家族, 為同一連鎖視窗的另一個家族, 為「一條最後被標成 的 read,實際來自家族 」的校準權重。 則真正用於 unit 的觀測機率為:
是正確標記的來源, 是由另一家族誤標進來的來源。 這些權重由成對家族的分子數與誤標率導出,不是每個 unit 自由擬合的比例。 因此 production likelihood 以 取代前一個基礎似然括號內的值。
這是前述成對觀測通道的摘要寫法。 由逐位點錯誤與整條分子的家族誤標算出,不是自由參數。 是地板總量,由逐位點錯誤 與家族誤標算出,不是另一個自由參數; 就是這個連鎖區段能分辨的最小比例,也是所有偵測下界的來源。
家族誤標尤其要成對估計,因為它連結同一連鎖視窗的兩個連鎖區段。 在 時它是假狀態的主要來源:只寫逐位點錯誤的話, 全 背景上出現的 在模型眼中「錯誤造不出來」(約 ), 於是只剩一個解釋 —— 一支新的 lineage。而它真正的來源是隔壁那個家族。
實作風險:少數錯字為何會改變整個候選集合
既有實作把狀態表當成精確值:沒有逐 read 錯誤模型, 沒有接受一個位點為「已覆蓋」的最低品質門檻,也沒有「一個樣式要幾條 read 才收」的下限。 在頻率路線上,一個誤讀只是把某個 VAF 推偏一點點; 在這條路線上不是 —— 一個被數條 read 共享的誤讀會引入一個樣本從來沒有的狀態, 而那個狀態會改變整個候選集合,連帶改變樹的形狀、旗標與計數。 這就是為什麼本節的錯誤地板不是數值穩定用的小常數, 而是這條路線唯一擋得住這件事的東西;也是為什麼前處理契約的三個門檻 (覆蓋一個位點的最低品質、算連鎖的最低 read 數、保留一個連鎖視窗的最低跨越深度) 必須寫死 —— 它們決定哪些字母有資格進到狀態表裡。
第三個細節:只准 0 → 1 時,LOH 會被誤讀
四配子相容性與 perfect phylogeny 都預設突變只增不減,所以 一段發生 LOHLOH 異型合子性喪失原本 heterozygous 的區域變成只剩一種 allele。LOH 不等於缺失 —— 也可能是一條 haplotype 遺失後另一條被複製(copy-neutral LOH)。Loss of heterozygosity: a formerly het region retains only one allele. Not necessarily a deletion.完整條目 → 而失去某個突變的區域,在模型眼中會被描述成 「那個狀態從來沒有取得過」,而不是「取得後又失去」。
來源補充:既有實作的 Camin–Sokal 條件在哪裡漏掉了「喪失」
既有實作的 Camin–Sokal 條件同樣只准 :復發(同一個突變在兩支各發生一次) 不被罰,但喪失完全不在模型的表達範圍內 —— 它不是被罰得很重,是連寫都寫不出來。 既有實作在整條鏈上沒有任何一處修正這件事,所以 LOH 區段的局部系譜會安靜地錯, 而且錯得跟一個正常結果長得一樣。
本規格的處理是把它變成明示的邊界,而不是默默承受:
拷貝數由全域基準提供,故每個連鎖區段都知道自己落在哪一種區段。
非中性區段的連鎖區段一律標 LOSS_UNMODELLED,
其局部系譜只作報告、不進任何樣本層統計量;
第一版並建議只納入 CN-neutral 的雜合區段。
這會減少可用連鎖區段(在 CN 變異廣泛的實體腫瘤可能減少一半以上),
所以納入與排除的連鎖區段數必須逐階段載明 —— 這正好也是連鎖率之外的第二個分層變數。
家族標號的方向不影響局部系譜
家族標號在每個 phase setphase 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.完整條目 → 內獨立決定, 故仍引入 並邊際化。但與前一篇不同的是, 這裡有一件可以放心的事:局部系譜的形狀完全不依賴 —— 連鎖區段內部的狀態集合與包含關係與家族標號無關,所以 不進核心推論,只進註解層。
作用範圍補充:那 到底在哪兩處還是要緊
兩處都在「把局部結果接到外面」的時候: 把節點對應到某一條 germline 單倍型(用於後文 caller 的單倍型特徵), 以及把節點對應到全域群。兩者都屬於註解層 —— 它們錯了,這個連鎖區段的狀態集合、比例與先後關係一個都不會變。
與前一篇的唯一結構差異
| 項目 | 前一篇(全域目標) | 本規格(局部目標) |
|---|---|---|
| 局部比例 | 全域 的確定函數,沒有自由度 | 自由參數,帶來自全域基準的先驗 |
| 換得的東西 | 每個連鎖區段都對全域參數施加限制 | 連鎖區段內部可以出現全域表達不了的結構 |
| 付出的代價 | 局部只能表達全域已經有的群 | 連鎖區段之間不再互相限制,後驗較寬 |
| 推論結構 | 全域取樣的內層包住每個連鎖區段 | 連鎖區段彼此獨立,可完全平行 |
這張表就是整份規格的決定點。 讓 自由,正是「提高局部解析度」的技術內容; 它同時也是「不再限制全域參數」的原因。兩件事是同一個改動的兩面,不能只要一面。
全域頻率基準如何成為獨立先驗
局部比例自由之後,20 條分子撐不起一個穩定的估計。 全域頻率譜在這裡仍然有用 —— 但只能單向進來。 可以這樣做的理由,正是前一篇為了讓似然相乘而建立的那個性質: 單變異連鎖視窗 與多變異連鎖區段 的 read 沒有交集。
全域頻率基準只使用
全域基準是一個標準的 頻率譜mutation frequency spectrum 突變頻率譜把一份樣本裡所有 somatic 變異的 VAF 畫成直方圖後得到的分布。分布上的峰與肩對應不同大小的細胞群,最低頻端的尾巴斜率則被用來判斷有沒有天擇。The distribution obtained by histogramming the VAFs of all somatic variants in a sample. Peaks and shoulders correspond to cell populations of different sizes; the slope of the low-frequency tail is used to test for selection.完整條目 →分群, 用哪一個既有實作都可以,它不是本規格的貢獻。它必須輸出四樣東西: 純度 ;CCF 原子的位置 與權重 ; 各原子的後驗不確定度;以及分層估到的 overdispersion 。
唯一的硬性要求是: 的位點必須整批排除。 它是這條路線不重複計數的全部依據,而且可以寫成檢查程式: 全域基準的輸入位點集合與 的交集必須為空。
這個要求的代價比紙上估計大得多,而且逐樣本差很多。 中篇以 Poisson 模型推得「約 9.5% 的變異落在多變異連鎖區段」,若真是如此, 排掉 對頻率譜幾乎沒有影響。但既有實作實測到的是 40.5% 到 98.1%(見「輸入與前處理契約」)—— 在連鎖率高的樣本上, 排掉 等於排掉大半個變異集合。 所以 剩幾個變異,是逐樣本必須先算、必須寫進報告的量, 不是一個可以忽略的尾數。
後果補充: 被排小之後,該傳什麼、什麼時候該直接棄權
太小時全域基準的原子後驗會變寬,那份不確定度必須在局部推論中逐次傳遞, 不能只傳原子的位置 —— 只傳位置等於把一個很寬的估計當成確定值使用。
而在極端情況下 可能小到撐不起分群。
那時全域基準就該直接棄權,全部連鎖區段標 UNLINKED_TO_GLOBAL,
而不是拿一個估不準的原子集合去當先驗 —— 後者會把全域的雜訊當成局部的先驗知識,
而且在輸出裡看不出來。
兩個邊界必須量出來,不能用講的:深度不是嚴格互斥,先驗不是硬約束
其一,深度不是嚴格互斥的。拷貝數分段與純度用的是全基因體深度, 其中含有 連鎖視窗的 read。一個 20 kb 尺度的連鎖視窗通常佔其所屬 CN 區段不到 1%, 所以這個洩漏有界 —— 但界線要實際算出來寫進報告, 嚴格版本則在 CN 分段時把 連鎖視窗一併排除,並比較兩者的差異。 其二,先驗不是硬約束。全域基準的原子若當成 只能取的值, 局部解析度就歸零了;下面那個逃逸質量存在的唯一理由,就是不讓這件事發生。
先驗:原子加上逃逸質量
先在細胞尺度上寫先驗,再換算成分子比例。設節點 的細胞比例為 :
由全域基準第 個原子的後驗均值與變異決定; 是弱的擴散成分。 是逃逸質量: 先驗上認為這一支不對應任何全域群的機率。 把局部結構強行貼回全域原子(解析度歸零); 等於丟掉全域資訊。 必須預先登錄並做敏感度分析。
由細胞比例換算為連鎖區段內的分子比例,用的是前一篇那條把拷貝數放進分母的式子:
一個 clone 在該處的每一條拷貝各貢獻一份,分母是該連鎖區段全部的分子數。
正常細胞是其中一項,其局部狀態恆為全參考 —— 所以純度不是額外的換算因子。
若該處的 CN 或 multiplicity 未定,這一步就停在分子尺度,輸出標
CELL_SCALE_UNAVAILABLE,而不是套一個預設的二倍體值。
先後關係的先驗:鴿籠pigeonhole 鴿籠原理由父代與子代的細胞比例限制樹形的算術規則:一個細胞至多屬於父節點底下的一個子節點,所以各子節點的細胞比例加起來不得超過父節點。名稱來自鴿籠原理 —— 東西放進籠子,總量不會憑空變多。在 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.完整條目 →只能軟用
由兩端節點所含變異在全域基準的群指派後驗算出: 若祖先端的群 CCF 不小於後代端,該邊獲得較高的先驗。 是狀態數的複雜度罰則。 兩者都必須是軟的:全域基準的群指派本身有不確定度, 把它當成硬性排序會讓局部推論繼承全域的錯誤而且不留痕跡。
為什麼不做完整的聯合模型:那正是前一篇與整合篇的規格。 本規格刻意切斷回饋(cut posterior),換三件事 —— 連鎖區段之間完全獨立因而可以平行、局部比例不會被全域參數壓平、 以及每一個連鎖區段的結論可以單獨檢視與反駁。 代價是這個後驗不是任何聯合模型的邊際分布, 報告中必須如此稱呼,不得寫成「聯合估計」。
候選模型與輸出旗標
每個連鎖區段的推論目標是 的後驗。 候選狀態集合由「觀測到的樣式 ∪ 相容性所需的中間狀態」產生, 再對每個候選算邊際似然( 依上節先驗積掉),最後正規化成後驗。 不取單一最大值,而是輸出前若干個候選與各自的質量。
每個連鎖區段除了那份後驗,還帶一組旗標。十二個旗標分成三類, 只有第一類是產出,另外兩類都是「這個連鎖區段的話只能講到哪裡」:
- 與全域基準的關係(四個):
REFINES_GLOBAL、GLOBAL_CONSISTENT、LOCAL_ONLY_LINEAGE、CONTRADICTS_GLOBAL—— 第一個是本規格的主要產出。 - 證據不足以作答(五個):
FLOOR_LIMITED、ORDER_UNRESOLVED、BRANCH_UNDERPOWERED、PARTIAL_DOMINATED、STATE_CAP_REACHED—— 每一個都代表「沒看到」,不代表「不存在」。 - 品質訊號(兩個):
INCOMPATIBLE、KATAEGIS_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 | 候選狀態集合觸及列舉上限 | 不得在截斷後仍回報通過 |
REFINES_GLOBAL 是這條路線的操作型定義:
「比全域方法解析度高」這句話在這裡等於「有多少個連鎖區段掛上這個旗標,
而且顯著多於虛無分布」。後文「驗收設計」定義這個虛無分布的產生方式。
沒有這個定義,前面所有論述都不可檢驗。
跨連鎖區段的彙整邊界
前一篇的規則是「可以合併的是機率,不是標籤」。 在本規格下更嚴格:連機率都不合併,因為 是每個連鎖區段自己的參數, 不同連鎖區段的後驗沒有共同的參數可以相乘。局部拓撲當然更不可以拼接成一棵大樹 —— 局部標籤只在該連鎖區段內有定義。
| 跨連鎖區段的產出 | 可否 | 條件 |
|---|---|---|
| 連鎖區段層級的目錄(每個連鎖區段一列 + 旗標) | 可 | 本規格的主要交付物 |
樣本層統計量:REFINES_GLOBAL 的個數、局部分岔比例、節點數分布 | 可 | 必須同時報告驗收設計的虛無分布 |
| 把局部樹拼成一棵全基因體的樹 | 不可 | 標籤跨連鎖區段無定義,且不同連鎖區段看的位點集合不同 |
| 把成對的祖先/互斥關係匯出給全域方法 | 有條件 | 使用它的全域擬合必須排除 的 read,否則就是重複計數 |
| 把節點數加總當成全基因體的群數 | 不可 | 節點是分子狀態,且各連鎖區段互不對齊 |
輸入與前處理契約
局部模型假設連鎖區段與稀疏狀態表已經建立。 前處理沿用既有工具,但其輸出直接決定候選狀態與覆蓋遮罩; 因此工具版本、品質門檻、HP3 處理與連鎖定義都是模型輸入契約的一部分。
下列參數會改變狀態表、連鎖區段與候選集合,因此必須與結果一起記錄, 以保證前處理可重現。
工具鏈與必須釘住的參數
既有的鏈是:正常樣本 germline calling(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.完整條目 →) → LongPhaseLongPhase實驗室開發的長 read phasing 工具,是後續所有工具的共同基礎。輸出 germline haplotype。The lab's long-read phasing tool and the common foundation for everything else. Produces germline haplotypes.完整條目 → 定相 → 腫瘤樣本 somatic calling(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.完整條目 →) → 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/somatic caller 的身分、版本,以及是否取交集 | 候選位點集合改變 → 連鎖區段改變 → 全部數字不可比 |
| 判定「這條 read 覆蓋這個位點」的最低 base quality 與 mapping quality | 決定 X 與 R/A 的分界,直接改變覆蓋遮罩 |
| 沒有 HP 標籤的 read 怎麼處理 | 丟掉與併入是兩種方向相反的偏誤,既有實作未載明 |
| 兩個 sSNV 要有幾條共同覆蓋的 read 才算連鎖 | 既有實作沒有下限,所以一條 read 就能開出一個連鎖視窗 |
| 一個連鎖視窗要有多少跨越深度才保留 | 同上;下限缺席時,證據極薄的連鎖區段與紮實的連鎖區段在計數上等重 |
標籤詞彙只有五個值,而 HP3 是一個真的問題
體細胞標記haplotagging根據已 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.完整條目 →輸出的詞彙恰好是
1、2、1-1、2-1、3 五個值。
1-1 表示這條分子屬於 germline 單倍型 1,且帶有可歸因於它的體細胞改變;
3 表示體細胞改變無法歸因於任何一條 germline 單倍型。
沒有更長的後綴,所以後綴不是一條獲得序列 ——
不可以把它讀成「先 1 再加一步」。不認得的標籤一律排除在所有分組之外。
HP3 破壞的正是「一個連鎖區段只看得到一條拷貝」
估計單位的核心假設是:連鎖區段內部不必處理兩條拷貝的混合。 HP3 的分子不滿足它,因為它們的體細胞改變兩條都歸不上去。 三種處理各有代價,必須明選其一並在輸出載明實際被影響的量: 整批排除(最保守,但在 LOH 與高拷貝區會丟掉大量證據); 當成獨立的第三個連鎖區段(形狀對,但那個連鎖區段的 CCF 換算沒有分母); 當成缺失並在兩個家族之間邊際化(統計上正確,成本最高)。 既有實作採第一種,但沒有報告被排除的量 —— 沒有那個量, 的分母就是錯的,而且錯得沒有跡象。
連鎖視窗不是固定寬度,是讀序連鎖的傳遞閉包
既有實作的連鎖視窗定義比「固定寬度 20 kb 視窗」精確得多,也更該沿用: 在同一個 phase setphase 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.完整條目 → 內,凡有某條 read 同時覆蓋至少兩個 sSNV,該段即為一個連鎖視窗; 若另一條 read 疊到已連鎖的位點又碰到新的位點,兩段合併成更長的連鎖視窗。三個推論:
- 數的是位點不是 read,所以一個由數條 read 接起來的連鎖視窗 可能一條全跨分子都沒有 —— 這正是覆蓋遮罩必須逐條保留的原因;
- 連鎖視窗不跨 phase set,因此所有狀態比較都在同一個定相框架內 —— 這是「不同 block 的 H1 不是同一條」那條警告的實作形式;
- 一個連鎖視窗依家族切成至多兩個連鎖區段,所以連鎖區段數大於連鎖視窗數, 兩者不可混用為同一個分母。
定義補充:前處理還要交出一個量 —— 作用位點數 (搜尋成本要用它,不是 )
前處理交出去的不只是連鎖區段本身,還有一個決定推論成本的量, 而它必須在這裡就算出來,因為它只有看得到狀態表的時候算得出來: 作用位點數 , 即在該家族內真的有變化的位點數。某個位點若該家族每條 read 都是同一個 allele, 它不產生任何座標,故 ;用 估搜尋成本會系統性高估。 本頁凡涉及狀態空間、頂點、遮罩與上限之處,一律以 為準。
連鎖率才是真正的閘門,而它逐樣本差很大
中篇以 Poisson 模型推得「約 9.5% 的變異落在多變異連鎖區段」。 實測不是這樣,而且它的變動幅度是這條路線最重要的一個經驗事實。 既有實作在七個資料集上量到:落入任一連鎖視窗的 sSNV 位點比例由 40.5% 到 98.1%,超過兩倍,而且不依癌別排序 —— 四個乳癌資料集自己就橫跨了幾乎整個範圍。
連鎖進來之後,多數連鎖區段只有兩個位點
第二個實測事實同樣改變了規格該長什麼樣: 是七個資料集中六個的最大單一類, 但集中程度由 76.7% 到 16.8%,只有 H2009 把 一路填到上限 12。 換句話說,這條路線大部分的產出,是一對位點之間的關係,而不是多層系譜。
這件事必須誠實地放進定位,但它不動搖前面那個證據預算,兩者常被混在一起: 佔多數限制的是深度(多數連鎖區段只講得出一對突變的關係), 不是可偵測性 ——「預期 0.2 條、實際 5 條」在 同樣成立, 而分辨兩支等值 CCF 的 lineage 本來就只需要一對位點。
逐項對照: 佔多數改變了什麼、沒有改變什麼
| 佔多數 | 影響 |
|---|---|
| 限制了深度 | 多數連鎖區段只講得出「這兩個突變共不共存、誰在前」,談不上多層局部系譜 |
| 沒有限制可偵測性 | 開頭那個「預期 0.2 條、實際 5 條」在 同樣成立 —— 它比的是「某個共現狀態存在」與「只有錯誤地板」,與位點數無關 |
對 REFINES_GLOBAL | 不影響:分辨兩支等值 CCF 的 lineage 只需要一對位點的共現或互斥 |
| 對交付物的描述 | 要改:主要交付物是「全基因體的成對分子關係目錄」,深層局部系譜是其中的少數 |
連帶的還有兩個直接後果。其一,caller 的局部系譜特徵覆蓋率不是固定的一成 —— 在高連鎖率的樣本上它可以覆蓋大部分候選,在低連鎖率的樣本上則接近無效。 其二,連鎖率是分層報告的第一個分層變數,必須在跑之前先算、 並用它決定這個樣本值不值得跑,而不是事後拿來解釋結果。
推論演算法與計算預算
- 前處理。依連鎖區段輸出稀疏狀態表,每列為 (連鎖區段、phase set、位點集合、覆蓋遮罩、字母樣式、計數)。遮罩不得先併掉。
- 就地校準。由 germline-only 連鎖區段、matched normal 與已知重複區估 逐位點錯誤率與家族誤標率,成對估計並分層借力。
- 全域基準。只用 ,輸出原子、權重、純度與 , 並斷言與 的位點交集為空。
- 候選狀態集合列舉。先取計數高於地板的樣式,再補相容性所需的中間狀態。 狀態空間是 : 時 4 格, 時 4,096 格 —— 所以上限必須存在。既有實作用的兩個上限與它們的實際行為見下。
- 邊際似然。對每個候選, 依上述全域頻率先驗積掉。 先驗是 Beta 混合、似然是 multinomial,非共軛, 故以每個候選數十個 quadrature 或 Monte Carlo 點計算並回報數值誤差。 這一步取代既有實作的 read-AF 排序:既有那個分數不只排序,它會先砍掉 30.03% 的最小成本解(見證據節),而且它沒有任何拷貝數校正 —— 當一個連鎖區段的每個位點都落在最高 read-AF 時,分數退化、全部並列。 換成邊際似然之後,不同形狀比較的是對同一組計數的解釋能力, 而且候選整組留著、由後驗質量表達,不先砍。
- 全域基準不確定度的傳遞。由全域基準後驗抽 組 ,逐組重算再平均。 過小會低估區間寬度,須在收斂契約中登錄。
- 輸出。每個連鎖區段一列,含前若干個 候選與質量、 節點比例的區間、旗標與有效分子數。
上限、棄權,以及棄權必須留在分母裡
既有實作用兩個上限:搜尋節點上限 與候選家族大小上限 。實測結果是 只有節點上限真的會觸發 —— 帶變異的連鎖區段中有 12.47% 撞到它。 兩個上限都是這套輸入與這台機器上的工程設定,不是關於 的數學結果; 換一台機器或換一個搜尋順序,數字就會移動,所以它們必須寫進 provenance。
撞到上限的連鎖區段要標記並計數,不可以丟掉
這是既有實作做對、而且值得整條抄過來的一件事: 撞上限的連鎖區段得到一個明示的棄權標記,而且仍然佔據每一個比例的分母。 被無聲丟棄的連鎖區段在完全相同的統計量上是隱形的 —— 報告中將無法判斷流程對多少比例的資料拒絕作答。 另有一條連帶的紀律:同一份報告至少有兩個分母 (全部連鎖區段、帶變異的連鎖區段),兩者回答不同的問題, 任何百分比都必須帶著自己的分母出現。
計數補充:同一個候選底下有幾棵等價的樹,不必抽樣(它沿頂點分解)
上限管的是「候選列不列舉得完」;另一件不同的事是列舉完之後, 同一個候選底下有幾棵等價的樹。 這個數目不需要抽樣,因為它沿頂點分解。 給定一個候選的頂點集合 ,每個非根頂點各自獨立地選一個父節點, 可選數就是它在 誘導子圖上的入度:
既有實作用這條式子精確計算並列數。在本規格下它有第二個用途: 它是簡約先驗在給定狀態集合下的正規化常數, 所以「這個候選有多少種等價寫法」不會被誤算成「這個候選有多少證據」。
成本不是這條路線的瓶頸
既有實作已經量過:整個基因體的局部重建是兩分鐘的工作,記憶體不到 400 MB。 本規格在它之外只加了每個候選的積分與全域基準的 次重算, 仍然遠低於任何全域取樣。要注意的只有一件事: 成本的變異來自 的分布而非連鎖區段數 —— 大的資料集吃掉絕大部分節點,所以要控成本,該控的是 的分布, 不是連鎖區段總數。
實測數字: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 倍就是上面「變異來自 的分布」的量: H2009 正是唯一把 一路填到上限 12 的資料集。
關鍵仍在於沒有全域取樣包在外層,所以整份工作是尷尬平行的, 成本隨連鎖區段數線性成長。這是放棄全域目標換來的直接好處, 也是它與整合篇那條路線在工程上的主要差別。
外部驗證:Somatic caller
局部系譜沒有完整金標準,因此另以 somatic calling 作為具有 truth set 的外部驗證終點。 此路線不參與局部系譜的核心推論;它只檢驗局部 lineage 一致性是否在 caller 既有訊號之外提供增量。
錨點與候選必須分開
這裡有一個明顯的循環:局部系譜要靠可信的 somatic 變異才建得起來, 而我們想幫的正是那些不可信的候選。解法是把兩者分成兩個角色:
| 角色 | 來源 | 用途 |
|---|---|---|
| 錨點 | germline het 位點 + 通過嚴格門檻的 somatic 呼叫 | 建構該連鎖區段的局部系譜 |
| 候選 | caller 在寬鬆門檻下的候選集合 | 被評分的對象,不進錨點集合 |
因此每個候選 的特徵,都是對把它拿掉之後建成的局部系譜 算出來的。這條留一法不是保險,是必要條件: 否則特徵會直接把「這個候選有幾條 ALT read」抄一遍, 在訓練時看起來很有效,在真正的低頻變異上完全沒有幫助。 錨點也不得取自 truth set —— 那是最典型的資料洩漏data 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.完整條目 →。
特徵群與可用範圍
局部系譜特徵群的核心是 llr_nested,而它的邏輯與 ClairS 的相位通道相同、
只是改用更局部的分類單位:真的低頻 somatic 變異來自某一支既有的 lineage,
所以它的支持分子應該整齊地落在該節點的後代位置上;
定序錯誤與嵌合則對局部系譜完全無所謂,會散落在各節點之間。
差別在於相位通道分的是兩條 germline 單倍型,
局部系譜特徵分的則是同一條單倍型內部的數支 somatic lineage。
欄位全表:三組特徵的每一個欄位與定義(實作時逐欄對照)
| 特徵群 | 欄位 | 意義 |
|---|---|---|
| 基礎 | vaf_grid_dist | 觀測 VAF 與「全域基準原子 × 該處 CN/multiplicity」所隱含的期望值格點的最近距離,以抽樣標準差為單位 |
| 基礎 | floor_ratio | ALT 分子數 ÷ 該連鎖區段錯誤地板的期望分子數 |
| 基礎 | purity、cn_local、phi_stratum | 樣本與區段層的脈絡量 |
| 單倍型 | hp_concentration | ALT 分子集中在單一 germline 單倍型的比例 |
| 單倍型 | hp_depth_ratio、hp_untagged_frac | 該單倍型的深度脈絡與未標記比例 |
| 局部系譜 | lineage_coherence | ALT 分子在 各節點上的集中度:最大節點佔比 |
| 局部系譜 | llr_nested | 對數似然比:「ALT 分子構成某一節點的巢狀子群」對「ALT 分子依 隨機散布」 |
| 局部系譜 | four_gamete_viol | 把候選加進 後是否產生四配子違反,及其超過地板的程度 |
| 局部系譜 | expected_vaf_local | 最佳配對節點的 換算出的期望 VAF,與觀測值的差 |
| 局部系譜 | n_span、k_u、state_entropy | 證據量與局部複雜度的脈絡量 |
| 全部 | tier0_ok、tier1_ok、tier2_ok | 該特徵群是否可用的指示通道 |
缺席不是零
沒有連鎖進任何連鎖視窗的候選就沒有局部系譜特徵群,而那個比例逐樣本由 1.9% 到 59.5% (連鎖率的補數)。把缺席欄位填 0 會讓模型學到「0 代表沒有支持」, 而那正好與「這個位點根本沒有可用的連鎖區段」相反。 每一特徵群都必須有獨立的可用性指示通道,訓練資料必須同時包含可用與不可用的兩種樣本, 而且驗收時三組各報一次 precision/recall,並另報不可用子集。 分層報告在這裡不是謹慎,是必要條件:可用比例本身逐樣本就差兩倍以上, 把 COLO829 那種 98% 與 HCC1395 那種 41% 平均起來, 一邊的增益與另一邊的傷害會同時消失。
三種整合方式
| 方式 | 做法 | 成本 | 風險 |
|---|---|---|---|
| A · 事後重新評分 | 以 caller 分數 + 特徵訓練一個小模型,只改排序 | 低,不需重訓 caller | 受限於 caller 已經丟掉的候選 |
| B · 加成輸入通道 | 把特徵接進 caller 的張量後重訓 | 高,需要完整訓練資料與算力 | 上限最高,但難以歸因是哪一組特徵有用 |
| C · 前置過濾 | 用特徵先砍候選 | 低 | 不建議:不可逆,且會安靜地傷害 recall |
建議順序是 A 先做。它便宜、可歸因(每一組特徵可以單獨加減)、 而且如果 A 沒有效果,B 幾乎不可能有 —— 因為 A 的失敗代表這些特徵 在 caller 已有的資訊之外沒有增量。
不論走哪一種接法,訓練與評估的切分有四條不能違反的規則, 其中最容易犯的是用逐變異的隨機切分 —— 同一個連鎖區段裡的候選共用一份局部系譜,隨機切分會讓指標系統性偏高。
切分紀律:四條規則逐條(切分單位、錨點來源、同一條管線、分層驗收)
- 不可以用逐變異的隨機切分。同一個連鎖區段裡的候選共用一份局部系譜, 彼此高度相關;隨機切分會讓同一個連鎖區段同時出現在訓練與驗證集, 指標會系統性偏高。切分必須以樣本為單位,連鎖視窗層另做一次巢狀檢查。
- 錨點必須由 caller 自己產生,不得取自 truth set 或 高信心區域high-confidence regionbenchmark 建立者指定的高可信度區間,通常以 BED 檔表示;區間外不宜視為具有相同標註可靠度。Intervals that benchmark authors consider reliable, usually supplied as a BED file. Results outside those intervals should not be assumed to have the same label reliability.完整條目 →的標註。
- 訓練與推論時的特徵必須由同一條管線算出;全域基準的實作、 、地板估計方式都要寫進 provenance。
- 驗收須跨純度、突變負荷、平台與 basecaller 分層報告, 以檢查分布偏移distribution shift 分布偏移測試資料的組成與訓練資料存在顯著差異(例如正負樣本比例不同),使得模型表現不如預期。Test data differing in composition from training data, so measured performance does not transfer.完整條目 →。
驗收設計
核心:局部打散的負對照
連鎖區段層級的結論沒有金標準,但有一個乾淨的負對照。 在每個連鎖區段內部打散分子的 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 計數 |
| 切分 | 與全域基準輸入位點的交集為空;訓練/驗證不共用樣本 | 重複計數或資料洩漏 |
實作順序(展開查看)
| 步驟 | 做什麼 | 可單獨驗證的內容 | 不得宣稱之事 |
|---|---|---|---|
| 前處理 | 稀疏狀態表(含遮罩)與就地錯誤率校準 | 能重現現行分類的計數;地板與 germline-only 連鎖區段一致 | 任何生物學結論 |
| 單區段模型 | 單一連鎖區段的局部後驗,先驗用平坦分布 | 機率核、相容性、邊際化三項測試 | 解析度優於全域 |
| 全域先驗 | 以 建立原子先驗與逃逸質量 | 敏感度;合成的等值 CCF 例子可分開 | 真實樣本的準確度 |
| 全樣本目錄 | 輸出目錄與 REFINES_GLOBAL 計數 | 打散負對照的效果量與虛無分布 | 細胞層級的群數 |
| Caller 基準特徵 | 基礎與單倍型特徵 + 事後重新評分 | 是否勝過 caller 既有的相位通道 | 局部系譜特徵的貢獻 |
| Caller 局部特徵 | 加入局部系譜特徵並分層報告 | 可用子集上的增益、不可用子集上的無害 | 總體平均掩蓋分層差異 |
| 外部重現 | 重現性與癌別關聯;匯出成對關係 | 跨 basecaller 一致度 | 沒有金標準的準確度宣稱 |
既有資料與文獻邊界
既有實作的資料規模
以下數值來自本實驗室既有實作在七個資料集、六個生物樣本、chr1–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%)。
逐資料集:七個資料集的無分支比例與已分類連鎖區段數
| 資料集 | 無分支合計 | 已分類連鎖區段數 |
|---|---|---|
| COLO829 | 88.5% | 10,757 |
| HCC1395_NYGC | 88.3% | 5,308 |
| HCC1395_HKU | 86.7% | 9,130 |
| HCC1954 | 84.9% | 5,647 |
| H1437 | 83.0% | 13,740 |
| HCC1937 | 82.1% | 4,245 |
| H2009 | 71.7% | 23,128 |
關鍵是為什麼稀少,而這個理由是機械性的: 一條鏈只要一個觀測到的樣式就建得起來,一個分岔要兩個互不包含的樣式。 所以分岔的稀缺,量到的至少有一部分是「這個連鎖區段帶了多少不同的樣式證據」, 而不是「底下的結構有沒有分岔」 —— 既有實作沒有把這兩者分開,本規格必須分開。
分法是現成的:錯誤地板決定了「第二個樣式若真的存在,看得到它的機率」。
所以每個判為無分支的連鎖區段都要附一個檢定力數字 ——
在該連鎖區段的分子數與地板之下,一支佔 5% 的旁支被看到的機率是多少。
檢定力低的無分支連鎖區段應標 BRANCH_UNDERPOWERED,
而不得計入「無分支佔多少」這個比例的分子。
沒有這一步,形狀組成量到的是覆蓋深度,不是演化。
還有一個落差常被接在這一段後面,但它的來源不是檢定力: 單一拓撲的比例由 35.3%(H2009)到 90.6%(HCC1954),差 55 個百分點。 那個排序與 的分布同向 —— 大的資料集並列的候選多、 解到單一拓撲的少,而不是它們的腫瘤比較單純。 也就是說,單一拓撲比例量的是「有沒有第二個樣式」, 的分布量的是「有幾個位點」, 兩者都不是演化訊號,而且要分開報。
合起來的結論是:庫存那張表裡 6%–22% 的分支類是上限,不是估計值。 判定分岔需要兩個單一狀態皆存在,且雙狀態的缺席不是抽樣零也不是錯誤地板 —— 後者正是局部生成模型中錯誤地板要檢驗的對象。
相關工作與新意邊界
| 既有能力 | 代表方法 | 本規格不可據此宣稱之處 |
|---|---|---|
| 由邊際 VAF 與拷貝數重建細胞層級 clone tree | PyClone、SciClone、PhyloWGS、Pairtree、DPClust | 不可宣稱全域樹準確度提升 |
| 把中性演化的尾巴與真正的 subclone 分開 | MOBSTER | 不可宣稱首次處理「頻率群集不等於 clone」 |
| 以成對位點的分子狀態(含相位)建 subclone 與樹 | PairClone、TreeClone | 不可宣稱把共現寫進似然為新意,也不可宣稱全域方法只吃邊際 —— 差別在可用位點對的數量級 |
| 相位輔助的等位/subclonal 拷貝數 | Battenberg、Refphase、HATCHet2 | 不可宣稱 phased CN 為新意 |
| somatic caller 使用 germline 相位通道 | ClairS、DeepSomatic | 不可宣稱單倍型感知呼叫為新意 —— 單倍型特徵群是基準不是貢獻 |
候選新意因此只剩兩處,而且都必須先做文獻檢索才能使用「首次」二字: 其一,把可變 的多位點共現當成估計目標本身(而不是當成全域樹的約束), 並以逐連鎖區段的旗標與打散負對照給出可計數的效果量; 其二,把同一條 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 已經把成對共現與相位寫進全域似然,所以候選新意不能再宣稱在「使用共現」這件事上, 只能落在「可用位點對多兩個數量級、 可變、逐連鎖區段交付」這幾點。 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 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。
延伸引用:全域重建路線、下游 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 calling,bioRxiv 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 DeepSomatic,Nat 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()與 subclonal()突變。不等於 VAF。
- clone tree(克隆演化樹)
- 描述腫瘤內各群細胞祖先關係的樹:節點是一群帶有相同變異組合的細胞,邊代表在祖先之上又多拿到變異。要注意同一組群集常常有多棵樹同時相容。
- copy number(拷貝數)
- 某段基因體在細胞內的拷貝數。多數正常常染色體區段為 2,可再分為 major 與 minor allele copy number。
- data leakage(資料洩漏)
- 訓練或模型選擇階段取得了評估資料的資訊,讓效能被高估。實驗室的做法是按染色體切分,確保同一個位點不會同時出現在訓練、驗證與最終測試中。
- distribution shift(分布偏移)
- 測試資料的組成與訓練資料存在顯著差異(例如正負樣本比例不同),使得模型表現不如預期。
- haplotagging
- 根據已 phase 好的 variants,把每一條 read 指派到 HP1 或 HP2,並把結果寫回 BAM 的
HPtag。 - 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 依
HPtag 分成的兩組之一。家族一是HP1與從它長出來的HP1-1,家族二是HP2與HP2-1;歸不到任何一條 germline 單倍型的HP3不屬於任何一族。叫「家族」是因為它把一條 germline 單倍型與由它衍生的 somatic 單倍型收在同一組裡 —— 分組看的是 germline 那一層,不是有沒有帶 somatic 突變。
在同一個 phase block 內,一個家族對應一條染色體拷貝,所以「兩族」就是那個位置上的兩條同源染色體。但兩件事不成立:其一,軟體判定不出哪一族來自父親、哪一族來自母親(那需要另外定序父母);其二,標號只在該 phase block 內有定義,跨 block 的「家族一」並非同一條染色體。