摘要
- 使用 R 探索哈代-温伯格平衡的概念
- 初步了解如何使用 R 进行编程
导入
本文的大部分章节将专注于使用 R 代码探索哈代-温伯格(HW)模型，并测试与哈代-温伯格期望(HWE)的偏差。HW 模型基本上是关于等位基因和基因型频率之间关系的理想化模型，在没有族群内繁殖、遗传漂变或选择等进化过程的影响下。因此，该模型作为测试这些进程是否发生的NULL模型(零假设)。换句话说，可以将实际数据与期望值进行比较，并推断出在真实世界中可能存在的族群演化或进化过程。除了 HW 模型和测试 HWE 偏离外，还将学会如何在 R 中模拟基因漂变，当然有在之前的文章简单的介绍过。本文将利用ggplot2 进行可视化，另外，也会遇到如何编写 for 循环、学习向量化并了解如何编写自定义函数的一些基础知识。
Hardy-Weinberg模型
HW模型是展示等位基因和基因型频率之间关系的一种方式。数学原理很基础，但非常重要，可以帮助我们理解不同的进化过程如何改变这些频率。该模型有五个主要假设：
- 所有个体在种群中是随机交配的
- 种群具有无限大小
- 没有选择发生
- 没有突变发生
- 没有基因流(gene flow)发生
假设在一个理想种群中，只专注于一个单基因座A，该基因座有两个等位基因A1和A2。我们可以使用p表示A1等位基因的频率，使用q表示A2等位基因的频率。这意味着如果我们知道p，我们可以推导出q，即q=1−p，因为p+q=1。HW模型是表示等位基因和基因型频率在世代间的关系的一种方法。在我们上述所列的假设下，决定基因型频率的唯一因素是它们的概率，这些概率是从等位基因频率中推导出来的。
确定这些概率非常简单。有2个等位基因A1和A2，我们有三种可能的基因型，可以通过以下三种方式确定它们的概率：
等位基因频率的可视化
可以用R些个程序并且简单的可视化一下，
library(tidyverse)
p <- seq(0, 1, 0.01)
q <- 1 - p
p
A1A1_e <- p^2
A1A2_e <- 2 * (p * q)
A2A2_e <- q^2
geno_freq <- as.tibble(cbind(p, q, A1A1_e, A1A2_e, A2A2_e))
geno_freq <- gather(geno_freq, key = "genotype", value = "freq", -p, -q)
a <- ggplot(geno_freq, aes(p, freq, colour = genotype)) + geom_line()
a <-a + ylab("Genotype frequency") + xlab("p frequency")
a + theme_light() + theme(legend.position = "bottom")
当然也可以用别的方法来实现，比方说for循环
# set up p
p_vec <- seq(0, 1, 0.01)
# initiate null vectors for the expected frequencies
A1A1_e <- NULL
A1A2_e <- NULL
A2A2_e <- NULL
# run the for loop
for(p in p_vec){
# define q
q <- 1 - p
# calculate genotype frequencies
A1A1_e <- append(A1A1_e, p^2)
A1A2_e <- append(A1A2_e, 2 * (p * q))
A2A2_e <- append(A2A2_e, q^2)
}
# combine them into a matrix (as an example here)
cbind(A1A1_e, A1A2_e, A2A2_e)
然鹅在R语言里，for循环并不是一个效率特别高的方式，追求高效率的话可以考虑向量计算
# set up p
p_vec <- seq(0, 1, 0.01)
geno_freq <- sapply(p_vec, function(p){
# define q
q <- 1 - p
# calculate probabilities
A1A1_e <- p^2
A1A2_e <- 2 * (p * q)
A2A2_e <- q^2
# write out
c(A1A1_e, A1A2_e, A2A2_e)
})
# transpose the matrix to get it the right way up!
geno_freq <- t(geno_freq)
检验是否偏离Hardy-Weinberg期望
之前已经学习了HW模型基本上是一个理想化的情境，其中没有任何人口统计或进化过程在塑造等位基因和基因型频率之间的关系。因此，可以将其用作对比我们数据的零假设模型。可以将观察到的基因型频率与在HWE下预期的频率进行比较。为了得到预期频率，只需像在本教程的前一部分中那样从等位基因频率中生成。
假设有一个A基因座，有两种等位基因A1和A2。这意味着有三种基因型，A1A1，A1A2和A2A2。从一个群体中抽取150个个体。下面的R代码块显示了我们收集到的每种基因型的数量，并将它们组合成一个观察向量，以备后用。
# numbers of genotypes
A1A1 <- 80
A1A2 <- 15
A2A2 <- 55
# make into a vector
observed <- c(A1A1, A1A2, A2A2)
接着需要从观察到的数据中计算出等位基因频率。这很简单 - 同型合子有两个等位基因的拷贝，而杂合子只有一个。因此，派生每个等位基因的数量，并将其除以总数：
# calculate total number of alleles
n <- 2*sum(observed)
# calculate frequency of A1 or p
p <- (2*(A1A1) + A1A2)/n
# calculate frequency of A2 or q
q <- (2*(A2A2) + A1A2)/n
# print frequencies to screen
c(p, q)
在上面的代码中，我们对观察到的基因型数量进行了求和，然后乘以2，因为每个个体有2个等位基因（即有150个个体和300个等位基因）。现在我们有了等位基因频率，因此我们可以使用我们在前一部分中计算出的代码来计算该群体的预期Hardy-Weinberg频率。
# generate the expected genotype frequencies
A1A1_e <- p^2
A1A2_e <- 2 * (p * q)
A2A2_e <- q^2
# combine into a vector
expected_freq <- c(A1A1_e, A1A2_e, A2A2_e)
计算出理论期待值，
expected <- expected_freq * 150
# combine observed and expected frequencies into a matrix
mydata <- cbind(observed, expected)
# add rownames
rownames(mydata) <- c("A1A1", "A1A2", "A2A2")
于是我们就得到了期待值和实际观测值，接下来要做的就很简单了，用卡方检验来判断是否正确。
chi <- sum(((expected - observed)^2)/expected)
1 - pchisq(chi, df = 2)
直接用R语言的卡方函数也行
mychi <- chisq.test(observed, p = expected_freq)
