部件分冊 · 第 12 頁 · InterSubMod（inter_sub_mod）

InterSubMod 吃什麼、做什麼、吐什麼

這是目前可編譯、可執行，並產出 per-region 甲基化矩陣與統計輸出的程式。七資料集 exact-PS 漏斗與候選拓樸數字來自另一條 research exact_ps_topology_af C++ solver + Python runner 流程，不是 inter_sub_mod 本體輸出。本頁把它的 3 個必填輸入、8 個處理階段、17 種輸出檔案全部攤開，每個檔案都附實際 head 出來的 header 與資料行。

一句話

給它一個腫瘤 BAM、一份參考基因組、一份體細胞突變清單（VCF），它會對每個突變位點開一個視窗，把跨過該位點的每一條 read 的甲基化狀態排成矩陣，算 read 兩兩之間的距離、做階層分群與統計檢定，全部落成檔案。

為什麼要這樣做：一般 bulk 定序把所有細胞的訊號平均掉了。ONT 長 read 可以在單一分子上同時看到「這條分子帶的是突變型還是正常型」與「這條分子上的 CpG 甲基化型態」，因此能在 read 層級（而非樣本層級）描述距離與 algorithmic clustering，並檢驗指定 mutation／haplotype label 是否與 read-distance variation 或 cluster assignment 關聯；這不證明 unsupervised cluster truth、cellular group 或因果。

實測可跑具名的 2026-08-06 內部 HCC1395 單點 receipt 記錄 2.9 秒、exit 0；含 normal BAM 與 LOH 標註的同一內部資料組合也 exit 0。這是特定機器、輸入與 commit 的歷史收據，不是一般 runtime 承諾。

● 用之前一定要知道的一件事 Git 只發布原始碼與可重建契約，不把 build output 當成受版控的發布物。先 checkout release manifest 指定的 immutable commit，再用 repo 外的新 build 目錄編譯；執行檔 inventory 應由該次 build receipt 動態記錄。

01輸入 — 3 個必填、2 個選填

下表的「實際路徑」都是磁碟上真實存在的檔案，可直接拿去跑。

參數	必填	這是什麼、程式拿它做什麼	限制與實測值
`--tumor-bam` `-t`	必填	腫瘤樣本的長 read 比對檔。程式只抓跨過該突變位點的 read，並從中讀取 MM/ML 甲基化標記與 HP 單倍型標籤。	需 `.bai` 索引（缺了只印 Warning 不擋，但隨機存取實際上會失敗）。實測：`data/bam/HCC1395/tumor.bam` = 283 GB，同目錄有 119 MB 的 .bai。
`--reference` `-r`	必填	參考基因組序列。用來找出視窗內所有 CpG 的位置（要先有序列才知道哪裡是 CG），也用來決定染色體長度以裁切視窗邊界。	● `.fai` 索引缺失是 Hard error，程式會直接停。實測：`data/ref/hg38.fa` → GRCh38（3.14 GB）。
`--vcf` `-v`	必填	體細胞突變清單。每一筆記錄 = 一個分析區域。 VCF 裡有幾個突變，就會產生幾個 region 目錄。	● VCF 的檔名會變成輸出目錄的第一層（取 stem）。實測：`filtered_snv_tp.vcf.gz` = 1.15 MB。
`--normal-bam` `-n`	選填	配對正常組織。給了它才會有：germline 甲基化基準線、腫瘤對正常的殘差、跨區域分層。不給的話這些欄位一律是 `NOT_APPLICABLE_TUMOR_ONLY`。	正常 read 是用整個視窗抓（不是點抓），因為它們貢獻的是甲基化基準而非突變判讀。原始碼註解明文：不要把 normal 的 HP 硬改成 0。
`--loh-bed`	選填	雜合性缺失（LOH）區間標註。用來在結果表填三個 LOH 相關欄位，讓後續分析可以把 LOH 區域分層看待。LOH BED 僅供 annotation／stratification； CN/LOH 未整合進 frozen candidate reconstruction，也沒有 CN/LOH-corrected inference。	只硬性要求前 3 欄，第 4 欄之後整行當註解字串。實測檔案 1,094 行，載入時印 `Loaded 1094 LOH regions`。

frozen release baseline ddd8909a 的 29-option census — 以及 兩個文件與程式碼不符的陷阱

這個 29-option 數量是 ddd8909a 的 branch-scoped source census，不是穩定 API。本輪把 --help 的每一行都對照 ArgParser.hpp 與 Config.hpp 逐一比對。大部分一致，但有兩個會讓你算錯資源或拿到非預期結果：

選項	help 說的	實際的	後果
`--threads` / `-j`	`Default: 1`	實際是 16 Config.hpp:43	不給這個旗標時，程式會用 16 條執行緒。資源估算會整整差 16 倍。實跑不帶旗標時 Configuration 印 `Threads: 16`。
`--distance-metric`	`Config.hpp` initializer 是 `BERNOULLI`	未指定時 effective CLI default 是 NHD ArgParser.hpp:86,181-193	parser 以 CLI list 取代 initializer；明示 NHD、BERNOULLI 或重複多值都會保留。多值時只有第一個 metric 驅動 clustering，後續只輸出額外 matrix，因此順序會改變 tree／cluster。

其餘經比對確認一致的預設值： --window-size 1000、--min-mapq 20、--min-read-length 1000、 --min-base-quality 20、--min-common-coverage 3、 --nan-distance-strategy SKIP、--methyl-high 0.8、--methyl-low 0.2、 --output-dir output、--log-level info。

● 另外三個設定是死的（在 Config.hpp 以外零引用）： min_site_coverage、pmd_gating、pmd_bed_path。也就是說「本方法有 PMD gating」「有最小位點覆蓋度過濾」是錯的敘述，不可寫進方法學章節。

● --help 的 exit code 是 1 不是 0 （help 與 parse error 共用同一條 return 路徑）。寫 CI smoke test 時 --help && next 這種寫法會誤判失敗。

02內部流程 — 一條 read 從 BAM 到結果，經過什麼

程式內部共 20 個細部步驟，下圖濃縮成 8 個關鍵階段。每個階段標了實際生效的參數。

圖 1 · InterSubMod 內部處理鏈（每個方框下方是實際生效的參數，非文件宣稱值）

最容易誤解的一步是 ⑥ 二值化：它不是把機率一刀切成 0/1，而是三分 —— 中間帶（0.2～0.8）被視為「不可信」直接丟棄。這代表曖昧的甲基化訊號與「根本沒測到」在下游是同一件事，分析時要意識到這個資訊損失。
另注意：程式裡其實有三套不同的二值化門檻（主線 0.8/0.2、另一模組硬寫死的 0.5、Python 端的 128/255），同一個 CpG 會因為走哪條路而得到不同判定，跨模組的數字不可直接比較。

03輸出 — 每個檔案長什麼樣（真實內容）

輸出分兩層：region 層（每個突變位點一個目錄）與 run 層（整次執行一份總表）。下面所有 header 與資料行都是實際從磁碟 head 出來的，不是照文件抄的。

Region 層 — 每個突變位點一個目錄

reads/reads.tsv 10 欄 — 這個位點有哪些 read

每條 read 一列：內部編號、真正的 read 名字（UUID）、座標、比對品質、單倍型標籤、支持突變型還是正常型、來自腫瘤還是正常、正反股。

read_id	read_name	chr	start	end	mapq	hp	alt_support	is_tumor	strand
0	56c0d051-1557-44af-99cc-46f5b9f48136	chr1	101290545	101347892	60	2-1	ALT	1	+

注意檔案在 reads/ 子目錄下，不在 region 根目錄。

methylation/methylation.csv 核心矩陣 — read × CpG 甲基化機率

每一列一條 read，每一欄一個 CpG，欄名就是基因體座標，格子裡是甲基化機率 0–1，沒觀測到寫 NA。

read_id,101330345,101330431,101330517,101330633,101330755,101330799,...
0,NA,0.0235,0.0275,0.9922,0.0667,0.0078,...

● 這裡有一個會靜默出錯的陷阱第一欄 read_id 是矩陣列號（整數），不是 read 的名字。要知道第 0 列是哪條 read，必須去 reads.tsv 用同一個 read_id 查。也就是說甲基化與 read 身分的綁定完全靠列序約定，程式沒有任何 key 校驗 —— 任何一邊做了過濾或重排而沒同步，就會把甲基化配到別條 read 的單倍型上，而且不會有任何錯誤訊息。

搭配檔 methylation/cpg_sites.tsv（3 欄）把欄號對回實際座標：

cpg_id	chr	position
0	chr1	101330345

distance/<METRIC>/matrix.csv — read 兩兩距離方陣

N×N 方陣，第一欄與 header 都是 read_id，6 位小數，無法計算的填 NA。 <METRIC> 是子目錄名（如 NHD、BERNOULLI），指定幾個 metric 就會有幾個子目錄。

read_id,0,1,2,3,4,5,6,7,8,...
0,0.000000,0.000000,0.142857,0.400000,0.333333,0.250000,...

同目錄另有 stats.txt（距離矩陣的體檢報告：用哪個 metric、多少對 read 有足夠共同 CpG、距離分佈）與分股版本 matrix_forward.csv / matrix_reverse.csv。

clustering/tree.nwk 最容易誤讀 — 階層分群樹

(((((((56c0d051-1557-44af-99cc-46f5b9f48136:0.000001,15a3bad6-8744-45e2-a577-...

● 這棵樹的葉子是 read，不是 clone 葉節點標籤是 read 的名字（BAM 裡的 QNAME，一串 UUID）。這是一棵「read 依甲基化相似度」的階層分群樹，不是細胞亞群的演化譜系樹。兩者是完全不同的東西，非常容易在報告裡被誤引。

同目錄的 leaf_order.txt 是樹上葉子從左到右的順序 —— 畫熱圖時要照這個順序排 read，圖才會跟樹對得上。

● clustering/linkage_matrix.csv 副檔名是 .csv，但實際內容是 TAB 分隔。用 pandas.read_csv() 預設的逗號分隔去讀，整列會變成單一欄位，而且不會報錯。必須寫 sep='\t'。

Run 層 — 整次執行一份

significance_summary.csv 跨版本不相容 — 下游幾乎只吃這一份

一列一個 region，把所有統計欄位攤平。實務上絕大多數下游分析都只讀這個檔。

RegionID,Chr,Pos,Ref,Alt,NumReads,NumCpGs,GlobalP,CramersV,GlobalP_HPFamily,...
0,chr1,877772,G,C,46,17,2.990000e-01,0.0000,3.065000e-01,0.0000,...

● 欄數會隨 binary 版本改變；現有欄位只有 component schema，沒有 whole-file layout version 磁碟上歷史實測到的欄數有 59 / 114 / 117 / 157 / 180 五種，而 frozen release baseline ddd8909a 的 source header 是 199 欄；歷史 73afaeac-dirty audit 不是 release source。意思是不同時期跑出來的結果，欄位位置是錯開的。任何用固定欄號（而非欄名）解析的下游腳本，會靜默讀到錯的欄位。現行欄位含 VerificationSchemaVersion=2 與 RegionStratificationSchemaVersion=1，但兩者都不是整檔 layout version；務必用欄名讀，並記錄 producer commit。

其他 run 層檔案	內容
`significance_statistics.txt`	一頁摘要：總共處理幾個 region、多少個顯著、各染色體分別多少。實測開頭 `=== Significance Analysis Statistics ===`、`Total Regions Processed: 1982`。
`run.log`	這次跑用了什麼輸入與參數，開頭是 `--- Configuration ---` 區塊。它是 provenance 的必要但不充分證據：仍須搭配 producer commit、site profile／tool hash、輸入與輸出 checksum、schema 與 command receipt；不要單靠 log 就宣稱可重現。
`subclone_structure.txt`	Legacy-named deprecation／region-stratification artifact；不是 inferred cellular subclones。只在部分 run 出現，檔名不得當作 biological claim。
`label_first_metrics.tsv`	36 欄，label-first 路線的精簡指標表，用 `chr:pos:ref:alt` 當鍵。只在 canonical 配對輸出見到。

文件與實際不符的一處 .claude/rules/output-structure.md 宣稱 region 目錄會有視覺化 *.png，但本輪在實際輸出目錄中找不到任何 PNG。圖是由 Python 層另外畫的（見分析與呈現層），不是 C++ 產生的。

04怎麼跑 — 三種可攜參數組合

先以 release manifest 的 immutable commit 建立 repo 外 clean build，再用 site profile 定位合法資料。下方結果來自具名內部 HCC1395 receipt；不是任意機器或輸入的固定 runtime／數值承諾。本機絕對路徑只登錄在 machine registry，不寫進公開主命令。

① 最小可跑 — 驗證環境與執行檔正常（runtime 依環境而異）

先從 VCF 抽出單一個突變，只給三個必填參數：

# 先把角括號 placeholder 換成自己的路徑；site profile 不進 Git
REPO_ROOT="<REPO_ROOT>"
DATA_ROOT="<DATA_ROOT>"
BUILD_ROOT="<CLEAN_BUILD_ROOT>"
SITE_PROFILE="<SITE_PROFILE>"
"$REPO_ROOT/scripts/site/doctor" --profile "$SITE_PROFILE" --mode real-preflight

# 準備一個只含第一筆 biallelic SNV 的 smoke-test VCF
SP="${TMPDIR:-/tmp}/ism_demo" && mkdir -p "$SP"
V="$DATA_ROOT/vcf/HCC1395/pileup/filtered_snv_tp.vcf.gz"
bcftools view -h "$V" > "$SP/one_snv.vcf"
bcftools view -m2 -M2 -v snps -H "$V" | head -n 1 >> "$SP/one_snv.vcf"

# 跑
"$BUILD_ROOT/bin/inter_sub_mod" \
  --tumor-bam  "$DATA_ROOT/bam/HCC1395/tumor.bam" \
  --reference  "$DATA_ROOT/ref/hg38.fa" \
  --vcf        "$SP/one_snv.vcf" \
  --output-dir "$SP/out_min"

具名內部 HCC1395 receipt（僅說明曾觀測到的輸出形狀）：

Total regions: 1 / Successful: 1 / Failed: 0
Total reads processed: 85
Forward strand (+): 40 / Reverse strand (-): 45
Total CpG sites found: 11
Metric: NHD / Total valid read pairs: 3443 / Total invalid pairs: 127

② 典型分析 — 指定執行緒、視窗、兩種距離指標

"$BUILD_ROOT/bin/inter_sub_mod" \
  --tumor-bam  "$DATA_ROOT/bam/HCC1395/tumor.bam" \
  --reference  "$DATA_ROOT/ref/hg38.fa" \
  --vcf        "$SP/one_snv.vcf" \
  --output-dir "$SP/out_typical" \
  --window-size 1000 \
  --threads 8 \
  --distance-metric BERNOULLI \
  --distance-metric NHD \
  --min-common-coverage 3 \
  --nan-distance-strategy SKIP \
  --log-level info

Configuration 區塊會印 Threads: 8 與 Distance Metrics: BERNOULLI, NHD，每個 region 目錄下會同時出現 distance/BERNOULLI/ 與 distance/NHD/ 兩個子目錄。 只有第一個 BERNOULLI 驅動 clustering；NHD 在這個順序只是額外 matrix。對調順序可能改變 tree／cluster。

要跑全基因組時把 --vcf 換成完整的 filtered_snv_tp.vcf.gz 即可 —— 但那會產生數萬個 region，屬於長時間計算，請用背景執行並先確認磁碟餘量。（本輪只以單一突變驗證同樣的參數組合可以跑通，未實跑全量。）

③ 完整配對分析 — 加上正常組織與 LOH 標註

"$BUILD_ROOT/bin/inter_sub_mod" \
  --tumor-bam  "$DATA_ROOT/bam/HCC1395/tumor.bam" \
  --normal-bam "$DATA_ROOT/bam/HCC1395/normal.bam" \
  --reference  "$DATA_ROOT/ref/hg38.fa" \
  --vcf        "$SP/one_snv.vcf" \
  --loh-bed    "$DATA_ROOT/loh/HCC1395/tumor_phased_LOH.bed" \
  --output-dir "$SP/out_full" \
  --window-size 1000 --threads 8 \
  --distance-metric BERNOULLI --distance-metric NHD

怎麼確認它真的用到了 normal 與 LOH（實跑結果）：

Loaded 1094 LOH regions from BED file
Total reads processed: 110        <- 85 腫瘤 + 25 正常（不給 normal 時只有 85）

# significance_summary.csv 中：
NTumorReads=85 · NNormalReads=25 · NormalBaseline_Mean=0.0890 · SampleASM_Delta=-0.010062

05陷阱清單 — 會靜默出錯的地方

這節整理「不會報錯、但結果是錯的」那類問題。每一條都有原始碼或實測依據。

嚴重	問題	會怎麼咬你
●	甲基化與 read 靠隱含列序綁定	`methylation.csv` 首欄是列號不是 read 名，唯一的綁定在程式內部。任何過濾或重排不同步 → 甲基化配到別條 read 的單倍型，無任何 assert 會抓到。
●	`significance_summary.csv` 跨 run 欄數不一致	歷史實測有 59/114/117/157/180 欄；frozen release baseline `ddd8909a` 的 source header 為 199 欄。其中第 187 欄 `VerificationSchemaVersion=2`，第 196 欄 `RegionStratificationSchemaVersion=1`；它們是元件 schema，尚非單一的 whole-file layout version。用固定欄號解析仍可能靜默錯位，一律用欄名讀。
中	二值化門檻有三套	主線 0.8/0.2 三分法、另一模組硬寫死 0.5 二分法、Python 端 128/255。同一個 CpG 因為走哪條路而甲基化判定不同，跨模組數字不可比。
中	`linkage_matrix.csv` 其實是 TSV	副檔名與分隔符不符。`pandas.read_csv` 預設逗號 → 整列變單欄，不會報錯。
中	`--threads` 預設是 16 不是 1	help 寫 Default: 1，實際 Config.hpp 是 16。資源估算差 16 倍，在共用機器上平行跑多個 job 時會嚴重超載。
中	`--distance-metric` 的 initializer 不是 CLI default；多值順序有語意	未指定時是 NHD；明示 selection 會保留。多值時只有第一個 metric 驅動 clustering，後續只輸出額外 matrix，因此順序可能改變 tree／cluster。
低	群數最少一定是 2	輪廓係數從 k=2 開始掃，所以 `optimal_k` 永遠不會是 1 —— 它不證明 unsupervised cluster existence/truth。PERMANOVA 只測指定 label 與 read-distance variation／centroid separation 在置換 null 下的 association，且須與 PERMDISP 合讀；Fisher／Cramér's V 只測 label association，均不證明 cellular group 或因果。
低	`--help` 回傳 exit code 1	help 與參數錯誤共用同一條 return 路徑。CI smoke test 會誤判失敗。
低	三個設定是死的	`min_site_coverage`／`pmd_gating`／`pmd_bed_path` 在 Config.hpp 以外零引用。不可在方法學章節宣稱本方法有 PMD gating 或最小位點覆蓋度過濾。

← 回系統全景下一頁：LongLineage 部件 → 深入：ISM 方法本體（既有頁）深入：甲基化讀取與篩選（既有頁）