Wright-Fisher模型解釋：為什麼瀕危動物容易滅絕——R語言：如何收集for循環中的每個數據並進行分析

01-26

最近師弟問了我一個關於群體遺傳學的問題，想了一天多，看樣子解決了。並且從這個題目中明白了為什數量大群體基因資源更加豐富，更不容易滅絕；而數量小的群體，基因資源更加匱乏，更加的脆弱這個自然現象。那麼，現在我就和大家分享一下我的學習心得。

這個作業的題目是這樣的：

2:In this problem,you are to use Monte Carlo simulations to check some of our theory concerning a Wright-Fisher model without mutation.Consider a model with N=40 diploids,p0=0.4 and t=55 generations. Run a simulation of this process,at least 1,000 times and calculate the mean and variance of the frequency at the 55th generation and compare to our theoretical expectations.

[Under the Wright-Fisher neutral model: Var(p_t) is approximately p0(1-p0)*(1-exp(-t/2N))]Also check the mean heterozygosity at generation 55 to see if that agrees with our theory and make a histogram of the distribution of p values in the 55th generation.

意思是一個使用 Monte Carlo模擬檢驗Wright-Fisher model，在種群中不發生突變的情況下，一個群體數量為40，雜合率0.4，並繁殖了55代的群體，通過將這個過程運行至少1000次，將得到的第55代的數據取平均值和算方差，從而驗證Wright-Fisher中的一些理論預期。

關於Wright-Fisher模型，維基百科的解釋是這樣的：

Wright–Fisher model

Think of a gene with two alleles, A or a. In diploid populations consisting of N individuals there are 2N copies of each gene. An individual can have two copies of the same allele or two different alleles. We can call the frequency of one allele p and the frequency of the other q. The Wright–Fisher model (named after Sewall Wright and Ronald Fisher) assumes that generations do not overlap (for example, annual plants have exactly one generation per year) and that each copy of the gene found in the new generation is drawn independently at random from all copies of the gene in the old generation.

認為一個基因里存在兩種等位基因，A或a。在二倍體的群體中包含N個個體，因此有每種基因都有2N個拷貝。一個個體在一個基因位點上可能有兩種類型的拷貝AA，aa（純和）或（雜合）。我們把純和的基因頻率計為p，雜合的基因頻率計為q。Wright–Fisher模型( 是用 Sewall Wright 和 Ronald Fisher這兩個名字命名的)，它假設每一代繁殖不存在重疊（比如說，一年生植物一生只繁殖一次),同時新生群體中的每個基因拷貝是基於上一代的基因拷貝情況隨機的產生的（也就是說，上一代的等位基因會隨機組合，它可能產生更多的純合體也可能產生更多的雜合體，這是隨機的）

The formula to calculate the probability of obtaining k copies of an allele that had frequency p in the last generation is then. Where the symbol "!" signifies the factorial function.

這裡有個公式用於計算一個某種基因型（雜合／純和）為k個拷貝且基因頻率為p的群體，在最後一代中，這個群體的基因頻率。這裡的「！」表示階乘功能。

（公式是不是看起來很複雜？但想靠自己去將上一代基因隨機分配次給下一代，並且有不斷重複此操作55次拿到最終結果，並且為了最終結果的客觀性，這個過程還要重複1000次，並且再求它們的均值和方差，靠自己手算，就算有計算器，也不知道要算到哪年，而且計算過程中能保證不出錯嗎？所以，是時候體驗計算機強大計算能力的魅力了）

本人是一個剛涉足生物信息學的大四／研究生，目前僅會用R語言來解決問題、得出結果並分析，如有不足或者更好的方法，請各位看官不吝賜教。我的代碼如下：

setwd("~/R/homework把循環得到的數據放進數據框里") #這是設置工作路徑，它指向的文

#件夾是用來存放你的R腳本，環境變數或得到的數據，方便在次使用和複習

drift <- function(N,p0,generation,pronumber){ #這裡是構建一個函數，N代表群體

#數量，p0是雜合體（或純合體，為了方便說明，後面都將p0認為是雜合體的基因型頻率）的#基因頻率，generation是繁殖的代數，pronumber是重複群體繁殖55代這個過程的次數#（我們也可以看成是有「pronumber」個群體進行了55代繁殖y）

M <- data.frame() #建立一個空的數據框用於後面收集後面for循環中產生的數據

for (m in 1:pronumber){ #從第 1 次開始，到第pronumber次

A <- 2*N*p0 #A為雜合體的等位基因數量

for(i in 1:generation){ #進行55代繁殖循

p <- A/(2*N) #p為雜合體基因型在每一代中的比例

A <- rbinom(1,2*N,p) #rbinom函數的作用是在2N個等位基因里，根據上一代

#的雜合基因型頻率，隨機的抽取出p所對應基因型的下一代的等位基因數量

M[m,i] <- p #把本次運行得到的數據放進數據框對應的位置上

}

return(M) #返回數據框M

}

Compara.value=drift(40,0.4,55,1000) #調用函數，並把結果放進一個向量里

#現在我們已經得到數據了，下面就可以「為所欲為」的分析了

V55=Compara.value[,55] #一般認為，第55代時，群體的變異趨於穩定，所以只使用每次運#行後得到的第55代雜合基因型頻率進行計算，即表格的第55列，將它們賦給一個向量

hist(V55,main="Wright-Fisher model",xlab="P-Value") #柱形圖觀察結果的分布

#能看到既有趨近於零的，也有趨近於一的，還有很多分布在0.1到0.9之間，但不是很突出。

#單憑這個圖並不能看出什麼直觀的結果，因此需要再作其它圖看看

boxplot(V55) #畫箱線圖

#可以看到中位數在0.4附近，上下四分衛數與中位數的距離相近，說明數據的分布時比較均勻#的，同時由圖預測第55代的均值應該在0.4附近

mean(Numeric_comapra)#計算均值

[1] 0.39315

var(V55)#計算方差

[1] 0.1146824

#可以看到均值接近0.4，就是我們開始用來計算的群體雜合率，方差大於0.1，結果中變異比較大。

#從這裡我們就可以看出小的群體容易發生遺傳漂變，等位基因中可能發生三種走向：A丟失a富集，a丟失A富集，A和a的基因頻率維持在原有水平附近。前兩種情況的產生可能是自然選擇對優勢等位基因選擇的結果，後者可能自然選擇對這一等位基因不明顯；然而，大自然是變幻莫測的，比如突然的地震、海嘯、瘟疫等自然災害，可能會讓這個群體的等位基因進行一次嚴峻的考驗，進而造成等位基因頻率發生改變，引發等位基因的丟失，即使這個小群體能在這次災害中倖免下來，那麼在後面的自然選擇中可能也會很艱難；因為先前丟失的等位基因可能有利於後面的生存，而先前保留下來的等位基因在之後的生存競爭中並不一定具有優勢，於是這個小群體除非發生基因突變，產生能夠適應新的自然選擇的基因，否則就只能無法生存，進而滅絕，但基因突變的頻率是很低的，而且大多的突變都是有害的，所以這個小群體滅絕的可能性很大，不努力提高群體數量，放在自然環境中，滅絕只是早晚的事。
#相對的，一個同一物種的大群體，比如說由上面所說的1000個小群體組成的群體，他們以小群體的形式獨立的存在，並且每個小群體見不發生基因交流。這樣的大群體在自然環緊中也會受到選擇，它們中有的小群體的A基因頻率提高，有的a基因頻率提高，有的A，a基因頻率基本維持不變，但因為它們屬於同一物種，所以整體來看，A，a的基因頻率在這個物種中沒有發生改變，僅在某個特定值附近波動。並且，自然情況下，大群體中的小群體間一般是會發生基因交流，這就使得等位基因能夠很好的保存下來，不容易發生丟失，即使遇到自然災害，群體中的一部分受到嚴重的自然選擇，某個基因頻率明顯提高，但當災害過去，這個大群體中未受到選擇的那部分又會將因為受到選擇造成而改變的基因頻率稀釋回原來的狀態，當新的自然選擇放生時，上面的過程又會重演，這是群體的自我調節。大的群體同樣可能發生整體性的基因頻率改變，但要要想實現這個結果，那麼一定是發生了大面積乃至全球性的自然災害，使大群體分布的主要區域都受到影響，但發生這種災害的可能性就如大型隕石撞擊地球一樣，發生的可能性極低。所以，大的群體是更加穩定的，不容易滅絕的。
#自然界中，自然選擇的往往不只是一個等位基因，同樣，某個適應環境的性壯也不是由一個基因就能決定的，這往往是多基因控制。所以小的群體，基因更易丟失，基因資源液更加的匱乏；而大的群體基因不易丟失，基因資源更加豐富，並且大的群體也更容易產生突變基因，雖然突變的結果大都不好，但龐大的群體能夠承受突變帶來的負面影響，並篩選出優勢突變基因，並予以保留。所以，我在此呼籲大家，無論如何都不要捕殺瀕危動物，你可能獲得一時的口福，能得到一絲的優越感，能滿足一時的新奇感。但你也加速了這個物種的滅絕，讓這個自然面前脆弱的群體雪上加霜，滅絕一個物種，一個種族，這都巨大罪孽。通過這個分析，我也希望大家能夠明白，為什麼我們要將一些瀕危動物（比如：華南虎、藏羚羊等）圈養或者劃定保護區，因為我們不進行人工干預的話，這些物種在自然環境中很難生存，並且因為無知的偷獵等認為破壞也會加把它們推向滅絕的邊緣。

討論：

其實，大群體不只是在第55代基因頻率穩定，無論第幾代，它們的基因頻率都是相對穩定，通過下面簡單的分析及可以看出來。

#畫曲線圖分析結果

plot(c(1:55),Matrix_compara[1,],type="l",ylim=0:1,xlab="Generation",ylab="Allele frequency") #作一個小群體的55代等位基因隨機變化圖

for(i in 2:1000){ #把剩下小群體的55代等位基因隨機變化圖

points(c(1:55),Matrix_compara[i,],type="l")

}

#我可以改變以上的群體數字來觀察這1000個群體的55代繁殖過程中等位基因變化情況：

#一個群體：

#十個群體：

#一百個群體：

#後面因為數據重疊導致圖像模糊，辨別不出線條，所以1000個數據的圖沒必要放上來。從圖上我們看一看出每根線條的走向都是不一定的，但它們在0.4兩側的分布似乎是均等，所以，我猜想在一個由1000個小群體的大群體中，等位基因的頻率應該是不變的，為了驗證我的想法我做了一下的計算：

Matrix_compara=as.matrix(Compara.value) #數據框轉換成矩陣

Numeric_comapra=as.numeric(Matrix_compara) #轉化數據類型

mean(Numeric_comapra)

[1] 0.3913907

var(Numeric_comapra)

[1] 0.06184234

#可以看到，全部1000個群體中的55代數據的平均數仍然在0.4附近，並且群體數據的方差雖然大於0.05，但小於0.1，可以說數據間的變異不大，因此我認為我的猜想合理。另外，我們還可以把這個1000看成自然界中影響一個正常群體的1000中因素，當這1000中因素相互作用時，群體中的等位基因不變，但如果某種因素丟失後，群體中的這個等位基因將受到嚴重影響。而自然因素遭到破壞的很多原因往往是由於人類的干擾。

（我懷疑沒人會把整篇看完，如果你看完了或者掌握了如何收集for循環中的每個數據以及這個群體遺傳學模型的意思，那點個贊讓我認識你）