如何使用 hap.py stratification？

發表於

基因體結構會影響變異分析的效度。舉例來說，在 homopolymer 或是 segmental duplicate 區域，依賴 PCR 的 Indel 分析準確性往往低於其他區域 1。因此，若能針對基因體不同區域個別評估效度，將有助於了解技術或流程的限制，進而提升效度分析的鑑別度。

為此，GIAB (Genome in a Bottle consortium) 維護了一系列 BED 檔，記錄基因體上的功能性區域、重複性區域以及高度多樣性區域等。用戶可以使用 hap.py 等工具，配合這些 BED 檔來評估其分析流程在各個區域的表現。

使用方法

要分區呈現效度分析結果，首先得準備一份 BED 檔。這個檔案每一列都記錄一個基因體區域，四個欄位分別表示染色體名稱 (chrom)、起始位置 (start)3、結束位置 (end) 以及區域名稱 (name)。

1
2
chr1    90 110 intron
chr1    290 310 exon

在執行 hap.py 時，需要添加 --stratification-region <tag>:<bed> 選項來讀取 BED 檔。其中 <bed> 是指 BED 檔的存放路徑，<tag> 則用於標示這些區域的分類，例如我們可以用 functional_regions 來代表 intron 和 exon 兩類區域。範例的指令如下：

1
2
3
4
hap.py truth.vcf query.vcf\
    -o output_prefix\
    -r reference.fa\
    --stratification-region functional_regions:input.bed

分析結果會記錄在 *extended.csv，用戶可以依據 “Subset” 欄位來篩選想關注的區域。”Subset” 的格式為 <tag>_<name>，結合了分類標籤與區域名稱。

Type Subset
SNP functional_regions
SNP functional_regions_intron
SNP functional_regions_exon

如果提供的 BED 檔沒有標記區域名稱，那麼整份檔案會被視為相同區域，”Subset” 將只會顯示 <tag>。

Type Subset
SNP functional_regions

Confident regions 和 stratification regions 有什麼差異？

Confident regions 和 stratification regions 的用途與定位不同。

Confident regions 用於界定可靠的分析範圍。由於定序與分析技術限制，即使經過嚴格品管與專家審查，仍無法確保某些區域發現的變異是否可靠。因此，GIAB 在建立 benchmark VCF 時，基於先驗知識排除了高度重複或含有結構變異的片段，定義出所謂的 high confident regions。

用戶可以透過 -f/--false-positives <bed> 指定這些區域，來限縮效度分析的範圍，唯有位在這些區域的變異會納入比較。

Stratification regions 則反映了比較的對象。用戶可以透過 --stratification-region <tag>:<bed> 呈現不同區域的效度分析結果。

換句話說，-f/--false-positives 影響資料篩選，--stratification-regions 則影響結果呈現。當同時指定這兩個選項時，hap.py 會先挑出 confident regions 內的變異來計算效度指標，再根據 stratification regions 呈現各區域的分析結果。

值得留意的是，即使 stratification regions 和 cofident regions 沒有交集，最終結果仍會完整呈現 stratification regions 所有區域的效度指標值。

舉例來說，假設我們以兩個內容相同的 VCF 檔作為 truth.vcf 和 query.vcf，

1
2
chr1	100	var1	A	C	1000	PASS	.	GT	0/1
chr1	300	var2	T	G	1000	PASS	.	GT	0/1

同時也定義了涵蓋這些變異的 Stratification regions。

1
2
chr1    90 110 intron
chr1    290 310 exon

如果 confident regions 介於一號染色體第 50 到 250 個鹼基之間，那麼 var2（以及 exon）不在 confident region 範圍內，所以不會納入效度計算。

1
chr1    50 250

然而，hap.py 仍會完整呈現所有 stratification regions 的效度統計量。

Type Subset TRUTH.TOTAL TRUTH.TP
SNP functional_regions 1 1
SNP functional_regions_intron 1 1
SNP functional_regions_exon 0 0

批次執行

hap.py 也支援批次分析，用戶需要準備一個 TSV 檔，在第一欄標註 <tag>，並在第二欄提供 <bed> 對應路徑。

1
2
functional_regions  functional_regions_regions.bed
mappability mappability_regions.bed

接著執行：

1
2
3
4
hap.py truth.vcf query.vcf\
    -o output_prefix\
    -r reference.fa\
    --stratification stratification.tsv

便能獲得以下分析結果。<tag> 的用途在批次分析時更為明顯，它能幫助用戶區別基因區域的分類。

Type Subset
SNP functional_regions
SNP functional_regions_intron
SNP functional_regions_exon
SNP mappability
SNP mappability_easy
SNP mappability_hard

如同單一 BED 檔的分析，如果沒有標註區域名稱，那麼整個 BED 檔案都會視為同樣區域。

GIAB genome stratifications

GIAB 提供多種 BED 檔，涵蓋了低複雜度、XY染色體、GC比例較高等性質不同的區域，這些都是現有技術較難準確分析的地方。

GIAB 也合併前述區域為 Union BED 檔，只分為 difficult 和 non-difficult 兩類，能較直觀地評估分析流程在具挑戰性區域的表現。

用戶能從 GIAB 的資料庫 下載對應版本與參考基因體的 BED 檔案，再從 Github 取得批次分析的 TSV 檔。

由於這些 BED 檔較大，建議挑選必要的內容下載，自己編輯批次檔來執行較為便利。

  1. 1.Olson et al. (2022). PrecisionFDA Truth Challenge V2: Calling variants from short and long reads in difficult-to-map regions. Cell genomics, 2(5).↩
  2. 2.Dwarshuis et al. (2024). The GIAB genomic stratifications resource for human reference genomes. Nature Communications, 15(1), 9029.↩
  3. 3.BED 檔的座標是從 0 開始計算。↩