共计 8633 个字符，预计需要花费 22 分钟才能阅读完成。
动动发财的小手，点个赞吧！
简介
ABBA BABA 统计（也称为“D 统计”）为偏离严格的分叉进化历史提供了简略而无力的测试。因而，它们常常用于应用基因组规模的 SNP 数据（例如来自全基因组测序或 RADseq）来测试基因渗入。
在本次实际中，咱们将联合应用可用软件和一些用 R 从头编写的代码来执行 ABBA BABA 剖析。咱们将剖析来自几个 Heliconius 蝴蝶种群的基因组数据。
流程
从多个个体的基因型数据开始，咱们首先推断每个 SNP 的等位基因频率。而后，咱们计算 D 统计量，而后应用 block jackknife 来测试与 D=0 的零冀望的显着偏差。最初咱们预计 f“混合比例”。
数据集
咱们将钻研三个物种的多个种族：Heliconius melpomene、Heliconius timareta 和 Heliconius cydno。这些物种的散布范畴局部重叠，人们认为它们在同源地区产生杂交。咱们的样本集包含来自巴拿马和哥伦比亚安第斯山脉西坡的两对同域种群 H. melpomene 和 H. cydno。在哥伦比亚和秘鲁的安第斯山脉东坡，还有两对同域种群：H. melpomene 和 H. timareta。最初，有两个来自外群物种 Heliconius numata 的样本，这是进行 ABBA BABA 剖析所必须的。
所有样本均应用深度全基因组测序进行测序，并应用规范流程为每个个体的基因组中每个位点获取基因型。数据通过过滤，仅保留双等位基因单核苷酸多态性 (SNP)，并且这些数据已进一步细化，以减小本教程的文件大小。
假如
咱们假如同域物种之间的杂交将导致 H. cydno 和来自东方的 H. melpomene 同域种族之间以及 H. timareta 和来自东部的 H. melpomene 相应的同域种族之间共享遗传变异。安第斯山脉。还有来自法属圭亚那的 H. melpomene 的另一个种族，它是 H. timareta 和 H. cydno 的异源性，它应该没有经验过与这两个物种最近的遗传替换，因而能够作为对照。
除了测试基因渗入的存在之外，咱们还将测试基因组的某些局部比其余局部经验更多基因渗入的假如。具体来说，咱们晓得 Z 性染色体上至多有一个位点会导致这些物种之间的杂交雌性不育，这表明一个物种的常染色体与另一个物种的 Z 染色体之间不相容。因而，与常染色体相比，咱们可能预期 Z 染色体上的渐渗缩小。
全基因组基因渗入测试
在最简略的表述中，ABBA BABA 测试依赖于基因组中与 ABBA 和 BABA 基因型模式相匹配的位点计数。也就是说，给定三个内群体和一个具备关系 (((P1,P2),P3),O) 的外群体，并给定代表每个群体的单个基因组序列（即 H1、H2 和 H3），ABBA 位点是其中 H2 和 H3 共享衍生等位基因（“B”），而 H1 具备先人状态（“A”），如外群样本所定义。同样，BABA 代表 H1 和 H3 共享衍生等位基因的位点。
疏忽反复渐变，只有当基因组的某些局部具备不遵循“物种树”的谱系，而是将 H2 与 H3 或 H1 与 H3 分组时，能力产生两种 SNP 模式。如果种群最近才决裂，因为谱系排序的变动，预计基因组的某些局部会呈现这种“不统一”的谱系。在没有偏离严格分叉拓扑的状况下，咱们预计基因组中大致相同比例的基因组会显示两个不统一的谱系 (((H2,H3),H1),O) 和 (((H1,H3),H2)，O）。因而，通过计算整个基因组（或其中很大一部分）的 ABBA 和 BABA SNP，咱们能够估算出两个不统一的谱系所代表的基因组比例，这意味着咱们预计 ABBA 和 BABA SNP 的比例为 1:1。例如，种群 P3 和 P2 之间的基因流动可能会导致偏差，只管它也可能表明其余突破咱们假如的景象，例如先人种群构造或可变代替率。
为了量化与预期比率的偏差，咱们计算 D，它是整个基因组中 ABBA 和 BABA 模式总和的差除以它们的总和：
D = [sum(ABBA) – sum(BABA)] / [sum(ABBA) + sum(BABA)]
因而，D 的范畴为 -1 到 1，并且在原假如下应等于 0。D>1 示意 ABBA 适量，D<1 示意 BABA 适量。
如果咱们从每个群体中有多个样本，那么计算 ABBA 和 BABA 位点就不那么简略了。一种抉择是仅思考来自同一群体的所有样本共享雷同等位基因的位点，但这会抛弃大量有用数据。更好的抉择是应用每个位点的等位基因频率来量化谱系偏差 ABBA 或 BABA 模式的水平。这实际上相当于应用每个位点的所有可能的四个单倍体基因组组来计算 ABBA 和 BABA SNP。因而，ABBA 和 BABA 不再是二元状态，而是 0 到 1 之间的数字，代表与每个谱系匹配的等位基因组合的频率。它们是依据每个群体中派生等位基因 (p) 和先人等位基因 (1-p) 的频率计算的，如下所示：
ABBA = (1-p1) x p2 x p3 x 1-pO
BABA = p1 x (1-p2) x p3 x 1-pO
实战
筹备
关上终端窗口并导航到将运行练习并存储所有输出和输入数据文件的文件夹。当初创立一个名为“data”的子目录并下载本教程所需的数据文件
mkdir data
cd data
wget https://github.com/simonhmartin/tutorials/raw/master/ABBA_BABA_whole_genome/data/hel92.DP8MP4BIMAC2HET75dist250.geno.gz
wget https://github.com/simonhmartin/tutorials/raw/master/ABBA_BABA_whole_genome/data/hel92.pop.txt
cd ..接下来，在 GitHub 上下载本教程所需的 python 脚本汇合
wget https://github.com/simonhmartin/genomics_general/archive/master.zip
unzip master.zip全基因组等位基因频率
为了依据群体基因组数据计算这些值，咱们须要首先确定每个群体中基因组中每个多态性位点的衍生等位基因的频率。咱们将依据应用 python 脚本提供的 Heliconius 基因型数据来计算这些值。输出文件已被过滤为仅蕴含双等位基因位点。频率脚本要求咱们定义人群。这些在文件 hel92.pop.txt 中定义。
python genomics_general-master/freq.py -g data/hel92.DP8MP4BIMAC2HET75dist250.geno.gz \
-p mel_mel -p mel_ros -p mel_vul -p mel_mal -p mel_ama \
-p cyd_chi -p cyd_zel -p tim_flo -p tim_txn -p num \
--popsFile data/hel92.pop.txt --target derived \
-o data/hel92.DP8MP4BIMAC2HET75dist250.derFreq.tsv.gz通过设置 –target 衍生，咱们取得了每个位点每个群体中衍生等位基因的频率。这是基于应用指定的最终种群（H. numata silvana，或“slv”）作为外群。该种群对于先人状态不固定的地点将被抛弃。
全基因组 ABBA BABA 剖析
为了理解 ABBA BABA 测试的工作原理，咱们将从头开始编写代码来进行测试。启动一个新的 R 脚本。这将使应用不同人群从新运行整个剖析变得容易。
- 用于 ABBA BABA 剖析的 R 函数
首先，咱们定义一个函数来计算每个站点的 ABBA 和 BABA 比例，并应用它们来计算 D atstistic。输出将是群体 P1、P2 和 P3（即 p1、p1 和 p3）中派生等位基因的频率。（外群中先人等位基因的频率在所有位点都为 1，因为咱们应用外群来辨认先人等位基因，因而能够疏忽）。
D.stat <- function(p1, p2, p3) {ABBA <- (1 - p1) * p2 * p3
BABA <- p1 * (1 - p2) * p3
(sum(ABBA, na.rm=T) - sum(BABA, na.rm=T)) / (sum(ABBA, na.rm=T) + sum(BABA, na.rm=T))
}- 读入咱们的等位基因频率数据。
freq_table = read.table("data/hel92.DP8MP4BIMAC2HET75dist250.derFreq.tsv.gz", header=T, as.is=T)这创立了一个名为 freq_table 的对象，其中蕴含每个 SNP 的派生等位基因的频率。咱们能够查看此表中的站点数量，还能够查看前几行以理解数据。
nrow(freq_table)
head(freq_table)请留神，前两列给出了支架的名称（即染色体）和每个位点在染色体上的地位。其余列是不同亚种的等位基因频率，如上图所示。
- D 统计量
当初，为了计算 D，咱们须要定义群体 P1、P2 和 P3。咱们将从一个显著且先前发表的测试案例开始：咱们将询问 H. melpomene rosina (mel_ros) 和 H. cydno chioneus (cyd_chi) 之间是否存在基因渗入的证据。它们别离是 P2 和 P3。P1 将是咱们的异种种群，来自法属圭亚那的 H. melpomene melpomene (mel_mel)。
咱们设置这些群体，而后通过提取这三个群体的所有 SNP 的派生等位基因频率来计算 D。
P1 <- "mel_mel"
P2 <- "mel_ros"
P3 <- "cyd_chi"
D <- D.stat(freq_table[,P1], freq_table[,P2], freq_table[,P3])
print(paste("D =", round(D,4)))咱们失去的 D 统计量（记住 D 从 -1 到 1 变动），表明 ABBA 超过 BABA。这表明来自巴拿马的 H. cydno chioneus (cyd_chi)与来自巴拿马的同域 H. melpomene rosina (mel_ros)比与来自法属圭亚那的异体 H. melpomene melpomene (mel_mel)有更多的遗传变异。这与同源的两个物种之间的杂交和基因流是统一的。
然而，咱们目前不晓得这个后果在统计上是否牢靠。特地是，咱们不晓得适量的 ABBA 是否均匀分布在整个基因组中。如果它是由基因组某一部分的奇怪先人造成的，咱们就不太置信存在重大的入侵。
为了测试统一的全基因组信号，咱们应用了 block-jackknife 程序。
- Block Jackknife
只管站点之间不独立，但 Jackknife 容许咱们计算 D 的方差。更传统的办法，即咱们随机对位点进行从新采样并从新计算 D，是不适合的，因为基因组中左近的位点具备类似的先人，使它们成为非独立的察看后果。
block jackknife 程序预计全基因组均匀 D 的所谓“伪值”的标准偏差，其中每个伪值是通过排除基因组的定义块来计算的，取全基因组均匀 D 和计算的 D 之间的差值当块被省略时。
为了思考链接站点之间的非独立性，块大小须要超过产生自相干的间隔。在咱们的例子中，咱们将应用 1 Mb 的块大小，因为咱们晓得连锁不均衡在远低于 1 Mb 的距离处衰减到背景程度。
运行 jackknife 过程的代码相当简略，但咱们不打算在这里编写它。相同，用于此目标的 R 函数在独自的脚本中提供，咱们当初能够导入该脚本。
source("genomics_general-master/jackknife.R")该过程的第一步是定义在每次迭代中将从基因组中省略的块。jackknife 脚本中的 get_block_indices 函数将执行此操作，并返回与每个块对应的“索引”（即频率表中的行）。它要求咱们指定要剖析的每个位点的块大小以及染色体和地位。
block_indices <- get.block.indices(block_size=1e6,
positions=freq_table$position,
chromosomes=freq_table$scaffold)
n_blocks <- length(block_indices)
print(paste("Genome divided into", n_blocks, "blocks."))当初咱们能够运行 block jackknifing 过程来计算 D 的平均值和标准误差。咱们提供之前创立的 D 统计函数 (D.stat)，该函数将在每次迭代中利用。咱们还提供每个站点的频率以及将用于从给定块中排除所有站点的块索引。
D_jackknife <- block.jackknife(block_indices=block_indices,
FUN=D.stat,
freq_table[,P1], freq_table[,P2], freq_table[,P3])
print(paste("D jackknife mean =", round(D_jackknife$mean,4)))
依据 D 的均值和标准误差的无偏预计，咱们能够计算 Z 分数来测试 D 是否显着偏离零。
D_Z <- D_jackknife$mean / D_jackknife$standard_error
print(paste("D Z score =", round(D_Z,3)))通常，大于 3 或 4 的 Z 分数被视为显着，因而在这种状况下，大量 Z 分数意味着与零的偏差十分显着。
- 预计 admixture proportion
D 统计量为基因渗入提供了强有力的测试，但它并没有量化已共享的基因组的比例。曾经开发出一种相干办法来预计“混合比例”f。
这种办法背地的想法是，咱们将察看到的 ABBA 超过 BABA 位点的适量与齐全混合下的预期进行比拟。为了近似齐全混合下的预期，咱们从新计算 ABBA 和 BABA，但用 P3 物种的第二种群代替 P2。如果您不足第二个群体，您能够简略地将 P3 样本分成两局部。在本例中，咱们有两个种群来代表每个物种，因而如果咱们应用 H. cydno chioneus (cyd_chi) 作为 P3a，咱们能够应用 H. cydno zelinde (cyd_zel) 作为 P3b)。
咱们须要编写本人的函数来计算 f。输出将是每个群体中衍生的等位基因频率，但当初咱们包含 P3a 和 P3b。
f.stat <- function(p1, p2, p3a, p3b) {ABBA_numerator <- (1 - p1) * p2 * p3a
BABA_numerator <- p1 * (1 - p2) * p3a
ABBA_denominator <- (1 - p1) * p3b * p3a
BABA_denominator <- p1 * (1 - p3b) * p3a
(sum(ABBA_numerator, na.rm=TRUE) - sum(BABA_numerator, na.rm=TRUE)) /
(sum(ABBA_denominator, na.rm=TRUE) - sum(BABA_denominator, na.rm=TRUE))
}咱们当初能够抉择 P3a 和 P3b，并预计 f。
P3a <- "cyd_chi"
P3b <- "cyd_zel"
f <- f.stat(freq_table[,P1], freq_table[,P2], freq_table[,P3a], freq_table[,P3b])
print(paste("Admixture proportion =", round(f,4)))这表明，H. melpomene 和 H. cydno 的同源基因组中有超过 25% 的基因组是共享的。混合比例能够解释为任何繁多基因组中外来血统的均匀比例。或者，它能够解释为该群体中基因组中任何给定位点的外来等位基因的预期频率。
咱们能够再次来预计 f 的标准差，并取得置信区间。Jackknife 块索引曾经计算出来，因而咱们能够简略地再次运行 Jackknife 函数，这次将 f 函数指定为运行每次迭代的函数。
f_jackknife <- block.jackknife(block_indices=block_indices,
FUN=f.stat,
freq_table[,P1], freq_table[,P2], freq_table[,P3a], freq_table[,P3b])95% 置信区间是平均值 +/- ~1.96 标准误差。
f_CI_lower <- f_jackknife$mean - 1.96*f_jackknife$standard_error
f_CI_upper <- f_jackknife$mean + 1.96*f_jackknife$standard_error
print(paste("95% confidence interval of f =", round(f_CI_lower,4), round(f_CI_upper,4)))染色体 ABBA BABA 剖析
- 所有染色体都显示基因渗入的证据吗？
下面，咱们钻研了整个基因组的基因渗入水平。假如每条染色体都有足够数量的 SNP，咱们能够在染色体程度上进行相似的剖析，以评估单个染色体上的渗入。
执行此操作的第一步是辨认频率表中与 21 条 Heliconius 染色体中的每一条绝对应的行。
咱们首先应用 unique 函数辨认数据集中存在的所有染色体名称。而后咱们须要辨认表中代表每条染色体的行。为此，咱们应用 lapply 函数，该函数屡次利用一个简略函数以创立 R 列表格局的组合输入。在这种状况下，咱们将应用染色体名称来利用该函数，并且咱们利用的函数将简略地询问表支架列中的哪些值对应于该染色体，利用 R which 函数。
chrom_names <- unique(freq_table$scaffold)
chrom_indices <- lapply(chrom_names, function(chrom) which(freq_table$scaffold == chrom))
names(chrom_indices) <- chrom_names这将创立一个蕴含 21 个元素的列表 – 每个染色体对应一个元素。每个元素都是表中来自该染色体的所有位点的向量。咱们能够通过对咱们刚刚创立的列表利用长度函数来查看每个染色体有多少个 SNP。
sapply(chrom_indices, length)（sapply 与 lapply 相似，只是它会尽可能简化输入，因而这里它返回一个向量，而不是向量列表）。
当初咱们能够应用这些索引来计算每个染色体的 D 值。咱们再次应用 sapply，这次利用 D.stat 函数并仅对每种状况下来自特定染色体的表中的行进行索引。
D_by_chrom <- sapply(chrom_names,
function(chrom) D.stat(freq_table[chrom_indices[[chrom]], P1],
freq_table[chrom_indices[[chrom]], P2],
freq_table[chrom_indices[[chrom]], P3]))咱们还须要来确定每个染色体的 D 是否与零显着不同。首先，咱们将定义用于每个染色体的块。
block_indices_by_chrom <- sapply(chrom_names,
function(chrom) get.block.indices(block_size=1e6,
positions=freq_table$position[freq_table$scaffold==chrom]),
simplify=FALSE)此命令返回列表的列表。这是一个蕴含 21 个元素的列表 – 每条染色体一个。每个元素都是一个列表，给出该染色体内每个块的索引。
咱们能够查看每条染色体的块数，以及每条染色体每个块的 SNP 数量。
sapply(block_indices_by_chrom, length)
lapply(block_indices_by_chrom, sapply, length)当初咱们应用折刀计算每个染色体 D 的 Z 分数。
D_jackknife_by_chrom <- sapply(chrom_names,
function(chrom) block.jackknife(block_indices=block_indices_by_chrom[[chrom]],
FUN=D.stat,
freq_table[chrom_indices[[chrom]], P1],
freq_table[chrom_indices[[chrom]], P2],
freq_table[chrom_indices[[chrom]], P3]))
D_jackknife_by_chrom <- as.data.frame(t(D_jackknife_by_chrom))
D_jackknife_by_chrom$Z <- as.numeric(D_jackknife_by_chrom$mean) / as.numeric(D_jackknife_by_chrom$standard_error)
D_jackknife_by_chrom咱们看到 1-20 号染色体均显示出显着的渐渗证据（Z > 4），而 21 号染色体（Z 性染色体）则没有。事实上，chr21 的 D 为负值，表明异域 H. melpomene 种群与 H. cydno 的变异多于同域 H. melpomene 与 H. cydno 的变异，只管差别并不显着。这表明与基因组的其余部分相比，性染色体上的基因渗入显着缩小，这与针对性染色体上基因渗入的等位基因的强烈抉择统一。如果存在一种或多种不相容性导致波及 Z 染色体上的基因座，这就是咱们所冀望的。
本文由 mdnice 多平台公布
