LongPhase：核心 phasing 流程

為什麼重要

LongPhaseLongPhase實驗室開發的長 read phasing 工具，是後續所有工具的共同基礎。輸出 germline haplotype。The lab's long-read phasing tool and the common foundation for everything else. Produces germline haplotypes.完整條目 → 是本教材後續工具共用的基礎，處理的是一般（非腫瘤）樣本的 phasing。 開始分析腫瘤資料前，先掌握這個基本流程。

它值得看清楚內部機制，原因是：後面 LongPhase-S 與 LongPhase-TO 的腫瘤專用處理， 都是接在這張圖與這套投票之上。不知道圖長什麼樣，就看不出腫瘤版改了什麼。

概念與互動

輸入與輸出

注意這個順序：先 call variant，再 phase，最後 haplotag。 LongPhase 本身不做 variant calling —— 它需要別人先給它 het 位點當錨點。

四種變異，一張圖

LongPhase 最關鍵的設計選擇，是把四種變異放進同一張圖一起定相， 官方文件稱為 co-phasing。理由來自長讀本身：一條 read 可能同時跨過一個 SNP、 一個甲基化位點、一個小的 indel 與一個 SV —— 這四件事被同一個分子綁在一起， 分開 phasing 就把這個證據丟掉了。

程式內部給每個節點記一個型別，原始碼裡的註解就寫著 0=SNP 1=SV 2=MOD 3=INDEL 4=tandem repeat INDEL。 型別不只是標籤，它會改變後面的判斷 —— 例如串聯重複區的 indel（型別 4）投票權重被降到 0.1， 因為那類位置的定序錯誤率特別高。

一條邊怎麼決定方向，又怎麼被丟掉

邊不是「有」或「沒有」，而是四個加權計數。跨過兩個位置的每一條 read， 依它在兩邊看到的 allele 落進 rr、ra、ar、aa 之一：

兩個細節值得記住，因為它們解釋了很多「為什麼這裡沒接起來」：

• 低品質的鹼基會被打折。邊的任一端 base quality 低於 --baseQuality（預設 12）時， 那條 read 對這條邊只貢獻 --edgeWeight（預設 0.1），而不是 1。
• 甲基化的門檻更嚴。SNP 與甲基化位點之間的邊，門檻從 0.7 收緊到 0.3 —— 甲基化訊號比較吵，要求證據更乾淨才採信。

heuristic phasing：加權投票

圖建好之後，LongPhase 不求全域最佳解。它沿著座標從左往右掃， 一次決定一個變異：前面最多 --connectAdjacent（預設 35）個變異各投一票， 票的權重就是上面那個 20 / 1 / 0.1，兩邊加總大的一邊贏。

這是一個典型的 heuristic：快、在長讀上通常夠準，但代價是它對兩件事敏感 —— 證據剛好平手，以及單一 read 主導。③ 那個特例保護就是為了後者而存在的。 另外，距離超過 --distance（預設 300000，約一個 centromere 的長度）的兩個變異不互相投票。

最後兩道信心門檻

圖定好相之後還有一步，程式把結果反過來檢查一次：先用變異的相位決定每條 read 屬於哪條 haplotype，再用已標好的 read 回頭重算每個變異。兩步各有一道門檻：

三種失敗模式

第一項與第三項在腫瘤樣本上可能更明顯，因為 LOHLOH 異型合子性喪失原本 heterozygous 的區域變成只剩一種 allele。LOH 不等於缺失 —— 也可能是一條 haplotype 遺失後另一條被複製（copy-neutral LOH）。Loss of heterozygosity: a formerly het region retains only one allele. Not necessarily a deletion.完整條目 → 與 copy number 異常會減少可用的 phasing 錨點。 LongPhase-S 與 LongPhase-TO 都針對這些腫瘤資料情境提供額外處理。

真實證據

最小可跑範例

# 1. 準備輸入（見 M4）
samtools faidx reference.fasta
minimap2 --MD -ax map-ont -t 10 reference.fasta reads.fastq -o aln.sam
samtools sort -@ 10 aln.sam -o aln.bam && samtools index aln.bam

# 2. germline calling（Clair3 或 PEPPER-DeepVariant，見 M5）
#    產生 SNP.vcf

# 3. phasing
longphase-s phase \
  -s SNP.vcf \
  -b aln.bam \
  -r reference.fasta \
  -t 8 \
  -o phased \
  --ont

# 4. haplotagging
longphase-s haplotag \
  -r reference.fasta \
  -s phased.vcf \
  -b aln.bam \
  -t 8 \
  -o tagged

co-phasing：把其他三種變異一起帶進去

# 甲基化要先用 modcall 產生 VCF，才能餵給 phase
longphase-s modcall \
  -b aln.bam \
  -r reference.fasta \
  -t 8 \
  -o modcall

# 四種一起 co-phase；每種變異會各自輸出一份 phased VCF
longphase-s phase \
  -s SNP.vcf \
  -b aln.bam \
  -r reference.fasta \
  --indels \
  --sv-file SV.vcf \
  --mod-file modcall.vcf \
  -t 8 \
  -o phased \
  --ont

注意：modcall 的 help 顯示 -b, --bam-file， 但選項表註冊的長選項其實是 --methylbamfile；短選項 -b 兩者都通。 同理 phase 的 help 寫 --sv-window 與 --sv-threshold， 選項表註冊的是 --svWindow 與 --svThreshold。

把圖直接印出來看

前面那三張圖不必只當示意圖 —— --dot 會為每一條染色體輸出一個 dot 檔，裡面就是實際建出來的邊：

longphase-s phase -s SNP.vcf -b aln.bam -r reference.fasta --ont --dot -o phased

# 每一行是一條邊：<位置>.<allele> -> <位置>.<allele>
# allele 相同 = 平行連接，相反 = 交叉連接
head -20 phased.chr20.dot

如何驗證輸出

# 在已定相的 heterozygous records 檢查 | 與 PS
grep -v '^##' phased.vcf | head -5

# BAM 裡應該有 HP tag
samtools view tagged.bam | grep -o 'HP:i:[12]' | sort | uniq -c

# 有多少 read 沒被標到？
samtools view -c tagged.bam
samtools view tagged.bam | grep -c 'HP:i:'

最後兩行的差就是沒過 --readConfidence（0.65）那道門檻的 read。 請將差值除以總 read 數，再依資料品質與分析目標判讀；比例偏高時，回到 M6 的排查清單。

參數對應到哪一步

參數很多，但每一個都落在前面某一張圖上。按階段記比按字母記容易得多：

判讀練習

實作練習

執行最小流程

# 僅執行小區域測試，以縮短重複測試時間
REGION="chr20:1000000-2000000"

longphase-s phase \
  -s SNP.vcf \
  -b aln.bam \
  -r reference.fasta \
  -t 8 \
  -o phased \
  --ont

longphase-s haplotag \
  -r reference.fasta \
  -s phased.vcf \
  -b aln.bam \
  --region "${REGION}" \
  -t 8 \
  -o tagged

把兩道門檻的效果量出來

這個練習把前面那張「兩道信心門檻」的圖變成實際數字。 調鬆 --readConfidence，看有多少 read 從沒有 HP 變成有：

# 預設 0.65
longphase-s phase -s SNP.vcf -b aln.bam -r reference.fasta --ont -o strict
# 調鬆到 0.55
longphase-s phase -s SNP.vcf -b aln.bam -r reference.fasta --ont -m 0.55 -o loose

# 各自 haplotag 之後比較 HP 標記率
for p in strict loose; do
  total=$(samtools view -c ${p}.tagged.bam)
  tagged=$(samtools view ${p}.tagged.bam | grep -c 'HP:i:')
  echo "${p}: ${tagged} / ${total}"
done

標記率會上升 —— 但那不代表結果更好。放寬門檻是把不確定的 read 也貼上標籤， 標記率與正確率是兩件事，要用有 truth 的資料才分得出來（見 M10）。

確認參數名稱

# 看目前版本接受的長選項
longphase-s phase --help

# help 與實際行為不符時，以選項表為準（檔案在 repo 根目錄）
grep -A40 'static const struct option longopts' Phasing.cpp

原始文獻與程式碼

LongPhase：Lin JH、Chen LC、Yu SQ、Huang YT， LongPhase: an ultra-fast chromosome-scale phasing algorithm for small and large variants，Bioinformatics，2022。程式碼在 github.com/twolinin/LongPhase。