Hybrid Indexを算出・可視化する

系統間または種間の交雑は進化、種間相互作用、保全など生物学の様々な側面において重要な現象です。交雑の実態を明らかにするためには、主に遺伝解析が用いられます。現在では、全ゲノム解析やゲノム縮約法（MIG-seqやRAD-seq、GRAS-Diなど）によって得られた一塩基多型（SNP）に基づいて、交雑の歴史や方向性、強さなどを推定する様々な方法が提案されています。ここでは、最もシンプルでよく使われるHybrid Index（hi）の算出およびトライアングルプロットによる可視化までの手順を紹介します。初心者向けの解説です！

環境・使用したプログラム

macOS 14.5

Stacks 2.64
vcftools 0.1.16
bcftools 1.21
R version 4.1.1

手順

1. データファイルの出力
2. 各種・系統の親個体（純系）を選定
3. フィルタリング1回目：親種・系統間でFstの高い座位を選別
4. フィルタリング2回目：ミッシングの多い座位、個体を除外
5. 残ったSNPに基づいてhiを計算

ひとつずつ確認していきましょう。
なお、フィルタリングの順番はどちらが先でもよいです。

データファイルの出力

ここでは、StacksからVCFファイルを出力するところからスタートします。最初から厳しめのフィルタリングをかけておきます。

populations -P ./stacks -M ./popmap.txt -O ./populations --write-single-snp --vcf -R 0.95 --min-mac 2 --max-obs-het 0.75 -t 12 

Hybrid Indexの算出の際には、ミッシング（欠損）の多い座位およびミッシングの多い個体をなるべく少なくすることが重要です。出力結果（Stacksの場合は*.log.distribsファイルや*.snps.vcfファイル）を見ながら、フィルタリングの強さやサンプル選定を検討します。当然、フィルタリングを強くすると利用可能な遺伝子座数が減るので、妥協ラインを探ります。

ミッシングの多い個体は除去は、popmap.txtから除く → populationsの再解析、という手順で行います。座位のミッシング率と個体のミッシング率は共に<5%を目指したいところ。

親個体の選定

次に、種間または系統間でFstの高い座位を選定するために親個体を選定します。一般的には、種間・系統間交雑が生じ得ない地域から得られた個体を使ったり、PCA解析やADMIXTURE解析の結果に基づいて交雑の痕跡がない個体を使うことが多いと思います。今回利用するtriangulaRでは、それぞれの親について、複数個体選定できます。

ここで言う「親」とは父親母親ではなく、今回比較する2つの種や系統の純系個体（交雑を経験していない個体）を指します。

次に選定した各親について、選定された個体をリストしたテキストファイルを作成します。ここでは、それぞれ5個体を選定し、parentA.txtとparentB.txtというファイルを作成しました。

parentA.txt

Ind_A01
Ind_A02
Ind_A04
Ind_A06
Ind_A11

parentファイルが準備できたら、すべての遺伝子座について、親種間のFstを選定した個体に基づいて計算します。ここでは、VCFtoolsを使いました。

vcftools --vcf ./testdata.snps.vcf --weir-fst-pop ./parentA.txt --weir-fst-pop ./parentB.txt --out FST

すると、FST.weir.fstのようなファイルが出力されます。これに遺伝子IDとFstの値が書かれていますので、エクセルやRで読み込みます。Fstの値でソートし、Fst = 1.0の遺伝子座数が十分（>50）あるか確認します。座数があまりに少ない場合、Fstの閾値を少し下げます。

遺伝子座数が十分に多い場合、Fstの値は1に近ければ近いほど良い。得られる座数が少ない場合でも、最低でもFst > 0.75 とする。

解析に使えそうな遺伝子座が選定できたら、ホワイトリスト（選定された遺伝子座の番号をリストしたテキストファイル）を作成します。

whitelist.txt

544871
1116302
1116356
2781341
3955052
5109923
5109949
5300952
5300964
（略）...

次にホワイトリストに挙げられた遺伝子座だけをVCFファイルから抜き出します。

bcftools view -T whitelist.txt -Ov -o testdata.filtered.vcf testdata.snps.vcf

これで得られたtestdata.filtered.vcfがHybrid Index解析に用いるインプットデータとなります。
このファイルには以下のフィルタリングを通過したSNPのみが含まれているはずです。

• ミッシングの多い遺伝子座を除去（populationsの-Rの値を高くする）
• ミッシングの多い個体を除去（popmap.txtから除いて再解析）
• 選定された親個体に基づいて計算されたFstの低い遺伝子座を除去（ホワイトリストに含めないことで除外）

triangulaRでの解析

ここからはRで解析を行います。
使うライブラリは最近論文が出たtriangulaRです。

非常にユーザーフレンドリーです。

# ライブラリの読み込み
library(vcfR)
library(triangulaR)

# データファイルの読み込み
data <- read.vcfR("testdata.filtered.vcf", verbose = F)
popmap <- read.table("popmap.txt")
colnames(popmap) <- c("id", "pop")

# Hybrid Indexの計算
hi.het <- hybridIndex(vcfR = data, pm = popmap, p1 = "sp1", p2 = "sp2")

# プロットの色の設定（サンプル、親種1、親種2）
cols <- c("#878787", "#1b7837", "#762a83")

# トライアングルプロット
fig_hi <- triangle.plot(hi.het, colors = cols, cex =2, alpha =0.75)

# トライアングルプロット（ミッシング率確認用）
fig_hi_missing <- missing.plot(hi.het, cex =2, alpha =0.8)

# 図の保存
ggsave("hi_R=0.95_FST>0.9.pdf", fig_hi, 
       width = 10, height = 7,
       units = "cm", dpi = 300)

Hybrid Index（hi）やヘテロ接合度（he）の値はhi.hetに格納されています。

結果の解釈

公式のこの画像を見たらすべて分かると思います（手抜き）。

• 純系の個体は右下と左下の頂点にきます（図のP1とP2）
  • ヘテロ接合度は0。なぜなら、親種の間で分化している（Fstが1に近い）SNPのみのデータセットを用いているからです
• F1は中央上の頂点、F2は三角形の中央付近にきます
  • 理論上、F1はほとんどすべてのSNPがヘテロになり、F2ではその割合は平均的には半分になります
• バッククロス（戻し交配）は交配した純系個体に寄ったところに位置します（図のBC1とBC2）
  • どちらの親種と交配したかによって、プロットの位置が大きく変わります。

データセットにミッシングが含まれている場合、プロットの位置が三角形の中央に寄る傾向があるようです。また、データ量（遺伝子座数）が少ない場合、データのばらつきが大きくなります。詳しくはWiens et al. (2025)のFig. 4を参照。

参考文献

• Wiens, B.J., DeCicco, L.H. & Colella, J.P. (2025) triangulaR: an R package for identifying AIMs and building triangle plots using SNP data from hybrid zones. Heredity. https://doi.org/10.1038/s41437-025-00760-2