摘要
- 使用 R 探索哈代-溫伯格平衡的概念
- 初步了解如何使用 R 進(jìn)行編程
導(dǎo)入
本文的大部分章節(jié)將專注于使用 R 代碼探索哈代-溫伯格(HW)模型，并測試與哈代-溫伯格期望(HWE)的偏差蓉冈。HW 模型基本上是關(guān)于等位基因和基因型頻率之間關(guān)系的理想化模型垃杖，在沒有族群內(nèi)繁殖逸吵、遺傳漂變或選擇等進(jìn)化過程的影響下煌恢。因此，該模型作為測試這些進(jìn)程是否發(fā)生的NULL模型(零假設(shè))漩绵。換句話說拴泌，可以將實際數(shù)據(jù)與期望值進(jìn)行比較魏身，并推斷出在真實世界中可能存在的族群演化或進(jìn)化過程。除了 HW 模型和測試 HWE 偏離外蚪腐，還將學(xué)會如何在 R 中模擬基因漂變箭昵，當(dāng)然有在之前的文章簡單的介紹過。本文將利用ggplot2 進(jìn)行可視化回季，另外家制，也會遇到如何編寫 for 循環(huán)、學(xué)習(xí)向量化并了解如何編寫自定義函數(shù)的一些基礎(chǔ)知識泡一。
Hardy-Weinberg模型
HW模型是展示等位基因和基因型頻率之間關(guān)系的一種方式颤殴。數(shù)學(xué)原理很基礎(chǔ)，但非常重要鼻忠，可以幫助我們理解不同的進(jìn)化過程如何改變這些頻率涵但。該模型有五個主要假設(shè)：
- 所有個體在種群中是隨機交配的
- 種群具有無限大小
- 沒有選擇發(fā)生
- 沒有突變發(fā)生
- 沒有基因流(gene flow)發(fā)生
假設(shè)在一個理想種群中，只專注于一個單基因座A，該基因座有兩個等位基因A1和A2矮瘟。我們可以使用p表示A1等位基因的頻率瞳脓，使用q表示A2等位基因的頻率。這意味著如果我們知道p澈侠，我們可以推導(dǎo)出q劫侧，即q=1?p，因為p+q=1哨啃。HW模型是表示等位基因和基因型頻率在世代間的關(guān)系的一種方法烧栋。在我們上述所列的假設(shè)下，決定基因型頻率的唯一因素是它們的概率棘催，這些概率是從等位基因頻率中推導(dǎo)出來的劲弦。
確定這些概率非常簡單。有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")
當(dāng)然也可以用別的方法來實現(xiàn)呼猪，比方說for循環(huán)
# 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循環(huán)并不是一個效率特別高的方式砸琅，追求高效率的話可以考慮向量計算
# 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期望
之前已經(jīng)學(xué)習(xí)了HW模型基本上是一個理想化的情境宋距，其中沒有任何人口統(tǒng)計或進(jìn)化過程在塑造等位基因和基因型頻率之間的關(guān)系。因此症脂，可以將其用作對比我們數(shù)據(jù)的零假設(shè)模型谚赎。可以將觀察到的基因型頻率與在HWE下預(yù)期的頻率進(jìn)行比較诱篷。為了得到預(yù)期頻率壶唤，只需像在本教程的前一部分中那樣從等位基因頻率中生成。
假設(shè)有一個A基因座棕所，有兩種等位基因A1和A2闸盔。這意味著有三種基因型，A1A1琳省，A1A2和A2A2迎吵。從一個群體中抽取150個個體。下面的R代碼塊顯示了我們收集到的每種基因型的數(shù)量针贬，并將它們組合成一個觀察向量击费，以備后用。
# numbers of genotypes
A1A1 <- 80
A1A2 <- 15
A2A2 <- 55
# make into a vector
observed <- c(A1A1, A1A2, A2A2)
接著需要從觀察到的數(shù)據(jù)中計算出等位基因頻率桦他。這很簡單 - 同型合子有兩個等位基因的拷貝蔫巩，而雜合子只有一個。因此，派生每個等位基因的數(shù)量批幌，并將其除以總數(shù)：
# 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)
在上面的代碼中础锐，我們對觀察到的基因型數(shù)量進(jìn)行了求和，然后乘以2荧缘，因為每個個體有2個等位基因（即有150個個體和300個等位基因）〗跃現(xiàn)在我們有了等位基因頻率，因此我們可以使用我們在前一部分中計算出的代碼來計算該群體的預(yù)期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語言的卡方函數(shù)也行
mychi <- chisq.test(observed, p = expected_freq)
