模組 3 · Sequencing principles and the value of long reads

定序原理與長讀的分析價值

比較 Illumina 與 ONT 的讀取方式與錯誤特性,並介紹 homopolymer、coverage、VAF 與 5mC。

約 45 分鐘建議先修:生物學

本模組學習目標

  • 說明長讀在建立長距離分子連結上的優勢與限制
  • 描述 ONT 的常見錯誤特性,並比較 indel 與 SNV 的判讀難度
  • 說明 homopolymer 為何是 indel 假陽性的常見來源
  • 說出 coverage、VAF 與 basecaller 版本各自會如何影響下游結果
  • 辨識 CIGAR 所記錄的 alignment 操作,並說明低支持度過濾如何影響 recall

為什麼重要

M1 已建立以下觀察:若 read 同時涵蓋兩個變異位置,且判讀可靠, 可支持兩者位於同一 DNA 分子。本模組將此觀察連結至定序技術,說明讀長差異、長讀的分析優勢與相應限制。

自學閱讀時,可先比較短讀通常能辨識各位置的變異、但較難由單條 read 建立遠距分子關係, 以及長讀可涵蓋較遠位置的特性;再學習平台名稱與錯誤類型。

瞭解這些限制很重要;後續介紹的 somatic-indel 方法, 即聚焦於降低此類錯誤並改善候選判讀。

概念與互動

提供的長距離資訊

短讀與長讀在跨越兩個變異位點時的差別 上半部是 Illumina 短讀:本示意以約 100 bp 的 read 為例,兩個相距數千鹼基的變異位點無法被同一條 read 跨過, 因此看得到兩個變異各自的存在,卻無法判斷它們是否在同一條分子上。 下半部是 Oxford Nanopore 長讀:單條 read 可達數千至數萬鹼基,可能同時跨過兩個位點, 並可在適當的修飾辨識模型下,由同一條分子推定甲基化訊號。 變異 A 變異 B 相距約 3,000 bp Illumina(二代定序) 本示意讀長約 100 bp · 準確度通常較高 沒有任何一條 read 能同時跨過兩個位點 → 兩個變異是不是在同一條分子上?無法判斷 Oxford Nanopore(三代定序) 讀長可達數千至數萬 bp · 錯誤特性依流程而異,可由模型推定甲基化 品質足夠時,read 同時涵蓋兩變異 → 支持同一分子(cis) 實心/空心小圓 = 甲基化狀態,也在同一條分子上
兩個變異位點相距 3 kb。上:在本示意中,約 100 bp 的 Illumina reads 可分別觀察兩個位點,但沒有單條 read 同時涵蓋兩點。下:一條較長的 ONT read 涵蓋兩點,並可在完成修飾辨識時關聯中間的甲基化訊號。

此圖說明短讀與長讀的主要資訊差異。上半部並非表示短讀無法偵測變異; 在本比較條件下,Illumina 可分別辨識兩個位點,且單鹼基錯誤率通常較低, 但不能由單條 read 建立 3 kb 的分子連結。實際效能依平台、化學版本、basecaller 與變異類型而異。

Illumina
讀長常見約 100–150 bp,依平台與建庫而異可達數千至數萬 bp,分布依流程而異
單鹼基準確度通常較高依化學、basecaller 與變異類型而異; 常較具挑戰
跨越多個變異限於片段長度與建庫設計較能建立長距離相位與分子連結
甲基化bisulfite-seq 是常見方法之一native-DNA 流程可搭配模型由訊號推定修飾
常見優勢高準確度的局部變異偵測長距離 phasing 與分子連結

長讀 indel 判讀的主要限制

ONT 與 Illumina 的錯誤特性不同。 是 ONT indel 錯誤的常見高風險區域, 例如 AAAAAA

原因在判讀原理上:

為什麼同一個鹼基連續重複的區段特別容易讀錯長度 奈米孔定序靠電流訊號判讀鹼基。 當序列是各種不同鹼基交替時,每個鹼基造成的電流變化明顯,容易分辨。 但當同一個鹼基連續重複時,電流會維持在同一個平台上, 六個與七個造成的訊號差異可能很小,因此增加長度判讀錯誤。 長度誤判可能形成表觀插入或缺失,是假陽性 indel 的常見來源之一。 一般序列:好判讀 序列 A C G T A C T G 電流訊號 每個鹼基造成明顯的高低變化 → 容易數 連續重複:難判讀 序列 T T T T T T T 電流訊號 一路維持在同一個平台 六個與七個差異可能很小 → 容易數錯 數錯一個的後果:形成表觀插入或缺失 參考序列 A T T T T T T G (六個 T) 讀出來 A T T T T T T T G 多數了一個 → 報成一個「插入」 這類假訊號單靠提高門檻未必能解決;僅憑目前資料可能難以與真實插入缺失區分。 可加入序列以外的證據,例如檢查支持 reads 是否集中於一致的 haplotype。
左上:一般序列的電流訊號隨鹼基組合變化。右上:連續重複時,六個與七個重複單元的訊號差異可能很小,因而增加長度判讀錯誤。下方:一個單元的誤差可形成表觀插入或缺失;僅憑目前證據可能難以與真實 indel 區分。

下表提供效能示例,用於比較 SNV 與 indel 的相對差異。 數值解讀需同時記錄 caller 名稱與版本、資料集、truth set 與指標定義:

CallerSNV F1INDEL F1
ClairS(示例)0.710.32
DeepSomatic(示例)0.810.53

在此示例中,indel F1 低於 SNV F1,表示 somatic-indel 的整體偵測效能仍有改善空間; 後續方法會利用更多 read-level 特徵處理此問題。

補充:重複序列中的等價比對位置

除了 homopolymer,另一個原因是 alignment 位置的模糊性: 同一個 deletion 在重複序列中可有多種等價表示位置,使不同 reads 的訊號分散。 這些差異會反映在 ;CIGAR 是描述 read 與參考序列對齊操作的字串。

三項核心數值

意思本教材示例值
某個位置被多少條 read 覆蓋tumor 50×,normal 25×
支持 alt allele 的 read 比例germline het 在簡化二倍體模型中 0.5;somatic 依 purity、copy number 與 clonality 而異
樣本中腫瘤細胞的比例合成資料用 0.2 / 0.4 / 0.6 / 0.8 / 1.0

以下以簡化數值示例說明低腫瘤比例如何減少可用訊號:

低腫瘤比例為什麼難:把數字算出來就懂了 在本簡化模型中,定序深度相同時,腫瘤 DNA 在樣本中的佔比越低,支持一個真變異的 read 數就越少。 在腫瘤 DNA 佔兩成且符合模型假設的樣本上,一個所有腫瘤細胞都帶有的變異期望約有五條 read 支持。 定序本身也可能在位置產生錯誤 read;背景若接近真訊號,會增加低比例樣本的判讀困難。 同樣 50× 的定序深度,腫瘤 DNA 佔比不同的結果 腫瘤 DNA 佔比 支持真變異的 read 數 好不好判斷 100% 25 較易判讀 60% 15 中等難度 20% 5 較難判讀 定序錯誤 各位置可能有 1–2 模型中的背景假設 在本模型中,腫瘤 DNA 佔 20% 時期望真訊號約 5 條,背景約 1–2 條。 兩者距離較近,過濾門檻會面臨 precision–recall 取捨。 門檻過嚴可能漏掉真變異,過鬆則可能保留較多假陽性;實際結果依資料與模型而異。
每個方塊代表一條支持變異的 read。在本模型中,tumor DNA fraction 為 100% 時期望約 25 條,20% 時約 5 條;下方另以每個位置約 1–2 條錯誤 read 作為背景假設。少量背景錯誤可能與低比例變異訊號重疊。

在二倍體、單拷貝、clonal、無偏取樣且 tumor DNA fraction 為 0.2 的簡化模型中, 期望支持數為 50×0.2÷2=5 條。訊號數量不足是主要限制之一;背景錯誤、品質與演算法也會影響可靠性。

basecaller 是實驗變數

是把電流訊號轉成 A/C/G/T 的那一步。不同版本的 basecaller 會產生不同的錯誤特性,所以:

同一個細胞株,不同 basecaller,可視為兩份不同的資料集

此記錄並非次要細節。跨 basecaller 的同株資料可作為穩健性評估的一部分;例如 HCC1395_HKUHCC1395_NYGC 來自同一細胞株,但流程不同。 因此分析時應將 (細胞株、平台、basecaller、、 caller 版本、benchmark 來源與混合方式)與結果一併記錄。

真實證據

合成不同 purity 的樣本

若要評估方法在不同 tumor DNA fraction 下的效能,可依計算後的比例混合 tumor 與 normal BAM, 建立預先設定的合成梯度:

# 1. 先量出兩個 BAM 的實際 coverage(本例量到 tumor 50×、normal 25×)
mosdepth --by 500 tumor_cov  tumor.bam
mosdepth --by 500 normal_cov normal.bam

# 2. 抽樣率不是目標 fraction,要從實測深度回推:
#      抽樣率 = (目標 fraction × 合成後總深度)÷ 來源深度
#    總深度取 25×(兩邊都拿得出來),目標 DNA fraction 0.6:
#      tumor  需要 0.6×25 = 15× → 15 ÷ 50 = 0.30
#      normal 需要 0.4×25 = 10× → 10 ÷ 25 = 0.40
samtools view -s 1.30 -b tumor.bam  > tumor_sub.bam
samtools view -s 1.40 -b normal.bam > normal_sub.bam

# 3. 合併成指定 DNA fraction 的樣本
samtools merge dna_frac_0.6.bam tumor_sub.bam normal_sub.bam
samtools index dna_frac_0.6.bam

-s 的小數部分是抽樣率,不是你要的 fraction。 兩個來源深度不同時直接填 0.60/0.40 會得到 30×:10× = 0.75,不是 0.6。 以 tumor 50×、normal 25×、合成後總深度 25× 為例,五個梯度的抽樣率是:

目標 tumor DNA fractiontumor 抽樣率normal 抽樣率
0.2-s 1.10-s 1.80
0.4-s 1.20-s 1.60
0.6-s 1.30-s 1.40
0.8-s 1.40-s 1.20
1.0-s 1.50不混入 normal

這樣控制的是 tumor DNA/read fraction,五個樣本的總深度也才一致 —— 否則深度與 fraction 兩個變數會同時在動,之後看到效能下降就分不清是哪一個造成的。

判讀練習

同樣是 50× coverage 的樣本,你在某個位點看到 3 條 read 支持一個 1 bp 的 deletion, 而且那個位置位於一段 TTTTTTT 裡面。

你會怎麼看待這個候選點?

展開答案

此候選應列為低可信並以其他證據複核,因為三項風險同時存在:

  • 只有 3/50 條 read 支持(VAF 0.06),支持數偏少、可信度較低
  • 它是 indel,在 ONT 資料中通常較具挑戰
  • 它落在 homopolymer 裡,是 ONT 較容易發生長度判讀誤差的區域

但請注意:列為低可信候選不等於可以直接刪除。 低 purity 樣本中的 genuine somatic 變異也可能只有少數 reads 支持。 僅依低支持度刪除候選可能降低

因此需加入更細緻的證據,例如支持 reads 是否具有一致的 haplotype 來源,而非只依數量判定。

此問題也可能受資料可辨識性限制:某些 false positives 在現有特徵下與真變異非常相近, 使模型無法可靠區分。 若 truth set 在這種困難區域的標註不一致,特徵相近的位點可能得到相反標籤, 模型因而難以學得穩定規則。評估時應將 label uncertainty 納入解讀。

實作練習

檢查資料讀長

讀長會影響可連結的變異距離,仍需同時考慮 coverage、品質與演算法。取得新資料時可先檢查:

# 讀長的統計摘要
samtools stats aln.bam | grep -E '^SN' | grep -E 'length|reads mapped'

# 讀長分布(看有多少條真的夠長)
# -F 0x900 排除 secondary 與 supplementary,否則同一條 read 會被數好幾次
samtools view -F 0x900 aln.bam | awk '{print length($10)}' | sort -n | \
  awk '{a[NR]=$1} END {print "中位數", a[int(NR*0.5)]; print "第 90 百分位", a[int(NR*0.9)]; print "最長", a[NR]}'

# 覆蓋度
mosdepth --by 500 out aln.bam && head out.mosdepth.summary.txt

估算簡化模型中的預期支持數

可使用資料的 coverage,估算在二倍體、單拷貝、clonal 與無偏取樣假設下的預期支持 read 數; 這不是經驗證的 detection limit:

預期支持 read 數=覆蓋度×tumor DNA fraction÷2

例:覆蓋度 50×、tumor DNA fraction 為 20% 時 50×0.2÷2=5 條。

若期望值僅為兩至三條,可靠偵測低頻變異可能較困難。資料量是限制之一, 仍需比較增加定序深度、改善品質與調整模型的相對效益。

學習檢核

本模組術語

CIGAR
描述一條 read 如何對上參考基因體的緊湊字串。例如 5M1I5M 表示兩端各有 5 個對齊欄位,中間有 1 個相對於參考序列的插入;M 可能代表吻合,也可能代表錯配。
ONT(奈米孔定序)
Oxford Nanopore Technologies。以電流訊號讀取 DNA,可產生長 read;在 native-DNA 流程與適當模型下,可由訊號推定 5mC,通常不需 bisulfite 轉換。
VAF(變異等位基因頻率)
在某個位點上,支持 alt allele 的 read 佔全部 read 的比例。VAF 不等於帶有這個突變的細胞比例
aneuploidy(非整倍體)
染色體數目異常。腫瘤裡非常普遍,也可能使 cellular purity 與 DNA fraction 不一致。
basecalling(鹼基判讀)
把定序儀的原始訊號(ONT 是電流)轉成 A/C/G/T 字母的步驟。不同 basecaller 版本會產生不同的錯誤特性 —— 所以「同一個細胞株、不同 basecaller」是兩份不同的資料集。
coverage(覆蓋度/深度)
某個位置被多少條 read 覆蓋。50× 表示平均每個位置約有 50 條 read 覆蓋,不代表它們均支持同一 allele。
homopolymer(同聚物)
同一個鹼基連續重複的區段,例如 AAAAAA。nanopore 在此類區段較容易發生長度判讀錯誤,是 indel 假陽性的常見來源之一。
indel(插入/缺失)
短的插入(insertion)或缺失(deletion)。在 nanopore 資料上通常比 SNV 更難可靠判讀,因為 homopolymer 與 alignment 位置漂移都可能製造假訊號。
long read(長讀)
單條可達數千至數萬鹼基的定序片段(例如 ONT、PacBio)。若可靠地同時覆蓋多個 variant,可提供它們位於同一 DNA 分子上的直接觀測證據。
provenance(資料來源履歷)
記錄資料來源、細胞株、定序平台、basecaller、reference build、caller、benchmark 與混樣方式。在本實驗室中,這些資訊是重要實驗變數,而非附註。
recall(召回率)
在指定評估範圍內,所有真陽性中被找出的比例:TP / (TP + FN)。若只對固定的 caller 候選集做後處理,recall 只能維持或下降,不能恢復 caller 從未輸出的真陽性;報告時應說清楚這個候選集範圍。
reference genome(參考基因體)
用於比對的標準序列(人類資料常用 GRCh38)。所有座標均相對於指定版本;更換 build 會改變座標系,座標不可直接比較。
tumour DNA fraction(腫瘤 DNA 比例)
樣本 DNA 中源自腫瘤的比例。與 tumor purity(細胞比例)在 aneuploid 或 WGD 的情況下會不一樣。
tumour purity(腫瘤純度)
樣本中腫瘤細胞所佔的比例。purity 越低,somatic 訊號被正常細胞稀釋得越嚴重,偵測越困難。