簡介

ABBA BABA 統(tǒng)計（也稱為 D 統(tǒng)計）為偏離嚴(yán)格的分叉進(jìn)化歷史提供了簡單而有力的檢驗。因此，它們經(jīng)常用于使用基因組規(guī)模的 SNP 數(shù)據(jù)測試基因滲入。

雖然最初開發(fā)用于基因滲入的全基因組測試，但它們也可以應(yīng)用于較小的窗口，從而可以探索基因滲入的基因組景觀。

在本次實踐^[1]中，我們將使用可用的軟件執(zhí)行基于窗口的 ABBA BABA 分析，然后在 R 中編寫代碼來繪制結(jié)果。我們將分析幾個 Heliconius 蝴蝶種群的基因組數(shù)據(jù)。

工作流程

從多個個體的全基因組測序的基因型數(shù)據(jù)開始，我們將運(yùn)行一個腳本來計算每條染色體上各個窗口中的混合比例的度量。然后我們將繪制圖表來測試有關(guān)適應(yīng)性基因滲入的假設(shè)。

數(shù)據(jù)集

我們將研究三個物種的多個種族：Heliconius melpomene、Heliconius timareta 和 Heliconius cydno。這些物種的分布范圍部分重疊，人們認(rèn)為它們在同源地區(qū)發(fā)生雜交。我們的樣本集包括來自巴拿馬和哥倫比亞安第斯山脈西坡的兩對同域種群 H. melpomene 和 H. cydno。在哥倫比亞和秘魯?shù)陌驳谒股矫}東坡，還有兩對同域種群：H. melpomene 和 H. timareta。最后，有兩個來自外群物種 Heliconius numata 的樣本，這是進(jìn)行 ABBA BABA 分析所必需的。

所有樣本均使用深度全基因組測序進(jìn)行測序，并使用標(biāo)準(zhǔn)流程為每個個體的基因組中每個位點(diǎn)獲取基因型。數(shù)據(jù)經(jīng)過過濾，僅保留雙等位基因單核苷酸多態(tài)性 (SNP)。該數(shù)據(jù)集包括來自 18 號染色體的 SNP 數(shù)據(jù)，已知該染色體攜帶特別感興趣的翅膀圖案基因座。

假設(shè)

我們假設(shè)同域物種之間的雜交將導(dǎo)致 H. cydno 和來自西方的 H. melpomene 同域種族之間以及 H. timareta 和來自東部的 H. melpomene 相應(yīng)的同域種族之間共享遺傳變異。

然而，并非基因組的所有部分都會受到同樣的影響。特別是，我們懷疑 18 號染色體上的翅膀圖案基因 optix 受到了強(qiáng)烈的選擇。 optix 的差異調(diào)節(jié)可以導(dǎo)致翅膀上紅色色素沉著的不同分布，如在 H. melpomene 的不同亞種中所見，或者導(dǎo)致紅色色素缺失，如在 H. cydno 中所見。

赫利科尼烏斯翅膀圖案可以警告捕食者它們有毒。有些物種參與繆勒擬態(tài)，有毒物種進(jìn)化得彼此相似，這有助于加強(qiáng)捕食者的學(xué)習(xí)。模仿可以通過相同翅膀圖案上的獨(dú)立收斂來實現(xiàn)，也可以通過適應(yīng)性基因滲入來交換翅膀圖案等位基因來實現(xiàn)。因此，我們預(yù)測不同物種的共模仿種群可能在 optix 附近顯示出過多的基因滲入信號。

我們對具有不同翅膀圖案的種群有著完全不同的期望。如果某個地區(qū)的捕食者認(rèn)為最常見的本地圖案有毒，那么擁有外來翅膀圖案的代價將會很高。同樣，任何具有中間翅膀圖案的混合個體也將面臨更高的捕食風(fēng)險。因此，我們預(yù)測，在具有不同翅型的種群之間，optix附近的基因滲入程度應(yīng)該會減少。

量化整個基因組

簡而言之，該檢驗使用三個群體和一個具有關(guān)系 (((P1,P2),P3),O) 的外群體，并調(diào)查 P2 和 P3 之間是否存在過多的共享變異（與 P1 和 P3 之間共享的變異相比） ）。

這種過量可以用 D 統(tǒng)計量來表示，其范圍從 -1 到 1，并且在無滲入的原假設(shè)下應(yīng)等于 0。 D > 1 表示 P3 和 P2 之間可能存在基因滲入（或其他可能導(dǎo)致偏離嚴(yán)格的分叉物種歷史的因素）。

該測試旨在用于全基因組規(guī)模。 D 統(tǒng)計量不太適合比較整個基因組的混合水平，因為它的絕對值取決于諸如有效種群大小等因素，而有效種群大小可能在整個基因組中有所不同。

其他教程中描述的 f 估計更好，因為它根據(jù)定義反映了混合比例，但它對小規(guī)模的隨機(jī)誤差高度敏感。因此，為此目的開發(fā)了一種稱為 fd 的統(tǒng)計數(shù)據(jù)，該統(tǒng)計數(shù)據(jù)對于使用少量 SNP 引入的誤差更加穩(wěn)健。傳統(tǒng)的 f 估計假設(shè) P3 是供體群體，P2 是受者群體，而 fd 則在逐個地點(diǎn)的基礎(chǔ)上推斷供體。

選擇人群

這些統(tǒng)計數(shù)據(jù)的解釋很大程度上取決于所選人群。首先，該測試對從 P3 滲入 P2 最為敏感，而不是相反。

其次，fd 應(yīng)被解釋為 P3 和 P2 之間不與 P1 共享的過度共享變異的量化。如果 P1 和 P2 之間存在持續(xù)的基因流，則任何從 P3 到 P2 的基因滲入都會被低估。最后，我們只能量化比 P1 和 P2 之間的分裂更晚發(fā)生的滲入。

因此，如果我們想要量化整個基因組中發(fā)生的最大可檢測量滲入，我們應(yīng)該選擇異域且與 P2 關(guān)系不太密切的 P1。不過，我們也可以通過測試這個特性來發(fā)揮我們的優(yōu)勢。如果我們選擇與 P2 共享正在進(jìn)行的基因流的 P1，那么測試將顯示 P2 和 P3 共享 P1 不共享的基因組部分。這對于識別翅膀圖案等位基因非常有用，因為這些通常是亞種在基因流中保持獨(dú)特的唯一基因組區(qū)域。

實戰(zhàn)

準(zhǔn)備

打開終端窗口并導(dǎo)航到將運(yùn)行練習(xí)并存儲所有輸入和輸出數(shù)據(jù)文件的文件夾?，F(xiàn)在創(chuàng)建一個名為“data”的子目錄并下載本教程所需的數(shù)據(jù)文件

mkdir?data

cd?data

wget?https://github.com/simonhmartin/tutorials/raw/master/ABBA_BABA_windows/data/hel92.DP8HET75MP9BIminVar2.chr18.geno.gz

wget?https://github.com/simonhmartin/tutorials/raw/master/ABBA_BABA_windows/data/hel92.pop.txt

wget?https://github.com/simonhmartin/tutorials/raw/master/ABBA_BABA_windows/data/chr18.LDhelmet_MLrho.w100.tsv

cd?..

• 接下來，在 GitHub 上下載本教程所需的 python 腳本集合

wget?https://github.com/simonhmartin/genomics_general/archive/master.zip
unzip?master.zip

滑動窗口分析

• 針對兩個不同的情況運(yùn)行分析 python 腳本。在這兩種情況下，P1 都是全域性 H. melpomene melpomene (mel_mel)。 P2 和 P3 是我們期望共享基因的兩個種群。在第一種情況下，我們量化了來自巴拿馬的 H. melpomene rosina (mel_ros) 和 H. cydno chioneus (cyd_chi) 之間的基因滲入。在第二個例子中，我們量化了來自秘魯?shù)?H. melpomene amaryllis (mel_ama) 和 H. timareta thelxinoe (tim_txn) 之間的基因滲入。

python?genomics_general-master/ABBABABAwindows.py?\
-g?data/hel92.DP8HET75MP9BIminVar2.chr18.geno.gz?-f?phased?\
-o?data/hel92.DP8HET75MP9BIminVar2.chr18.ABBABABA_mel_ros_chi_num.w25m250.csv.gz?\
-P1?mel_mel?-P2?mel_ros?-P3?cyd_chi?-O?num?\
--popsFile?data/hel92.pop.txt?-w?25000?-m?250?--T?2

python?genomics_general-master/ABBABABAwindows.py?\
-g?data/hel92.DP8HET75MP9BIminVar2.chr18.geno.gz?-f?phased?\
-o?data/hel92.DP8HET75MP9BIminVar2.chr18.ABBABABA_mel_ama_txn_num.w25m250.csv.gz?\
-P1?mel_mel?-P2?mel_ama?-P3?tim_txn?-O?num?\
--popsFile?data/hel92.pop.txt?-w?25000?-m?250?--T?2

我們?yōu)槟_本提供一個包含基因型數(shù)據(jù) (-g) 的輸入文件、一個輸出文件 (-o)、內(nèi)群體群體和外群體群體（-P1、-P2、-P3 和 -O）以及一個指定每個群體的文件示例位于（--popsFile）中。

我們還給出了窗口的參數(shù)。這些 ae“坐標(biāo)”窗口，這意味著每個窗口相對于參考基因組的長度相同，但每個窗口的 SNP 數(shù)量可能不同。窗口大小 (-w) 將為 25,000 bp。 Windows 需要包含至少 (-m) 250 個 SNP 才能被視為有效。

最后，我們告訴腳本使用兩個線程 (-T)。如果你有一個多核機(jī)器，你可以增加這個值，腳本會運(yùn)行得更快。

• 繪制窗口統(tǒng)計數(shù)據(jù)

我們需要將每個窗口統(tǒng)計文件加載到 R 中。我們將創(chuàng)建一個包含兩個數(shù)據(jù)集的列表。首先輸入輸入文件的名稱

AB_files?<-?c("data/hel92.DP8HET75MP9BIminVar2.chr18.ABBABABA_mel_ros_chi_num.w25m250.csv.gz",
????????????????"data/hel92.DP8HET75MP9BIminVar2.chr18.ABBABABA_mel_ama_txn_num.w25m250.csv.gz")

AB_tables?=?lapply(AB_files,?read.csv)

head(AB_tables[[1]])

當(dāng) D 為負(fù)數(shù)時，fd 毫無意義，因為它旨在僅在存在過量時量化 ABBA 相對于 BABA 的過量。因此，我們在 D 為負(fù)數(shù)的位置將所有 fd 值轉(zhuǎn)換為 0。

for?(x?in?1:length(AB_tables)){
AB_tables[[x]]$fd?=?ifelse(AB_tables[[x]]$D?<?0,?0,?AB_tables[[x]]$fd)
????}

然后，我們可以為我們分析的兩種情況繪制染色體上的 fd 圖。

par(mfrow=c(length(AB_tables),?1),?mar?=?c(4,4,1,1))

for?(x?in?1:length(AB_tables)){
????plot(AB_tables[[x]]$mid,?AB_tables[[x]]$fd,
????type?=?"l",?xlim=c(0,17e6),ylim=c(0,1),ylab="Admixture?Proportion",xlab="Position")
????rect(1000000,0,1250000,1,?col?=?rgb(0,0,0,0.2),?border=NA)
????}

這表明染色體滲入程度存在相當(dāng)大的異質(zhì)性。如果我們考慮 optix 周圍的區(qū)域，我們就會看到 H. melpomene rosina 和 H. cydno chioneus 之間基因滲入減少的證據(jù)，正如我們預(yù)測的那樣。相比之下，我們看到了 H. melpomene amaryllis 和 H. timareta thelxinoe 之間基因滲入程度升高的證據(jù)，這表明它們共享的翅膀模式可能是適應(yīng)性基因滲入的結(jié)果。鑒于這一證據(jù)，建議對 optix 周圍區(qū)域進(jìn)行系統(tǒng)發(fā)育，以測試 H. timareta 等位基因是否似乎“嵌套”在 H. melpomene 進(jìn)化枝內(nèi)。

基因滲入和重組率之間的關(guān)聯(lián)

理論預(yù)測，如果存在許多選擇基因滲入的“屏障位點(diǎn)”，我們應(yīng)該會看到低重組區(qū)域基因滲入減少的趨勢，因為中性基因滲入等位基因和有害基因之間的聯(lián)系更強(qiáng)。我們可以通過檢查重組率和整個染色體的 fd 之間的關(guān)系來檢驗這個假設(shè)。我們有一個先前生成的數(shù)據(jù)文件（已提供），給出了該染色體上 100 kb 窗口中估計的群體重組率。

現(xiàn)在我們將制作一個 100 kb 窗口的 fd 匹配數(shù)據(jù)集，這里僅使用顯示滲入率最高的物種對：來自巴拿馬的 H. melpomene rosina 和 H. cydno chioneus。

python?~/Research/genomics_general/ABBABABAwindows.py?\
-g?data/hel92.DP8HET75MP9BIminVar2.chr18.geno.gz?-f?phased?\
-o?data/hel92.DP8HET75MP9BIminVar2.chr18.ABBABABA_mel_ros_chi_num.w100m1.csv.gz?\
-P1?mel_mel?-P2?mel_ros?-P3?cyd_chi?-O?num?\
--popsFile?data/hel92.pop.txt?-w?100000?-m?1000?--T?2

現(xiàn)在，回到 R 中，讀入這個新數(shù)據(jù)文件。

AB_table_w100?<-?read.csv("data/hel92.DP8HET75MP9BIminVar2.chr18.ABBABABA_mel_ros_chi_num.w100m1.csv.gz")

和之前一樣，我們將 D 為負(fù)的窗口的所有 fd 值轉(zhuǎn)換為 0。

AB_table_w100$fd?=?ifelse(AB_table_w100$D?<?0,?0,?AB_table_w100$fd)

現(xiàn)在我們讀取 100 kb 窗口的重組率表。這里，ML_rho 列給出了每個窗口的群體重組率的最大似然估計。

rec_table?<-?read.table("data/chr18.LDhelmet_MLrho.w100.tsv",?header=T)
head(rec_table)

由于數(shù)據(jù)的噪聲性質(zhì)，我們想要比較不同重組率的 bin 中的 fd 值。我們將使用剪切函數(shù)將窗口分成三個具有低、中和高重組率的箱。

rec_bin?<-?cut(rec_table$ML_rho,?3)

最后，我們可以制作箱線圖來比較這些箱之間推斷的混合水平 (fd)。

boxplot(AB_table_w100$fd?~?rec_bin)

這證實了基因滲入水平確實隨著重組率的增加而增加，這與大量屏障位點(diǎn)在全基因組范圍內(nèi)選擇基因滲入的模型一致。

Reference