通過前兩周的《本地化適應是怎麼發生的?》和《突變是否影響個體的適應性?》了解了群體的核酸多樣性後,我們接下來就開始要著手進行功能基因的定位了。工欲善其事,必先利其器。在我們可以自由選用各類實驗設計前,我們需要了解各種方法的基本原理。讓我們先從連鎖分析開始。
既然群體中產生了多樣性,我們就期望將與性狀相關的基因定位出來。在之前的文章中,我們提到功能基因定位的方法主要包括QTL定位(包含GWAS)和群體遺傳(選擇壓力分析)。這裡的QTL定位是廣義上的QTL定位,包括經典的連鎖分析和關聯分析。在這裡,我們先介紹一下連鎖定位的方法原理。
連鎖分析,之所以被稱為「連鎖分析」,其本質還是利用功能基因與分子標記間的連鎖與重組,實現對功能基因位置的定位。
例如下圖中,Q基因型會導致個體變高,q基因型會導致個體變矮。我們可以看到,鄰近的基因座Bb與Qq基因座連鎖。B總是與Q連鎖,導致B基因型的個體總是更高,對應b基因型的個體更矮。而遠離Qq基因座的Ee基因座則沒有這種現象。由於兩者距離較遠,彼此間沒有必然的連鎖關係(傾向自由重組),因此我們可以看到E基因座對應的既有高的個體,也有矮的個體。
在實際研究中,這些分子標記ABCDE都是位置已知的標記,但我們不知道Qq基因座的位置。如果通過數學的方法,我們發現Bb、Cc基因座與性狀高矮相關,而其他基因座並非如此,我們就可以確定功能基因Qq就位於Bb和Cc之間。
要創造這樣的1個基因型分離的人工群體,我們往往需要使用雜交的策略(兩種不同的親本雜交),在之前小師妹的微信文章《遺傳圖譜各種作圖群體介紹》曾經介紹過各類作圖群體的來源方式,感興趣的讀者可以再看看。
圖1與功能基因的連鎖與自由組合
正如上文所說的,我們需要挖掘確認哪些分子標記與性狀關聯,從而進一步推斷影響性狀的功能基因與這類分子標記連鎖,從而判斷功能基因位於該分子標記附近。在統計學上,我們使用最簡單的方差分析,也可以實現這樣的推斷。
如圖2,我們可以將整個群體按照Bb基因座的基因型分為BB基因型群體和bb基因型群體的兩個子群體。如果我們使用方差分析證明子群體BB的平均身高顯著大於bb,則證明Bb基因座與性狀相關。類似,我們將會發現按照Ee基因座分類的兩個子群體在平均身高上則沒有區別。這樣我們就可以推斷,由於Bb基因座與性狀相關,那麼決定身高的基因座Qq應該位於Bb附近,這樣就實現了QTL初步定位。
圖2 使用方差分析進行單標記分析
再看看圖1的示意圖,我們是否可以將其看成1個線性回歸方程組:
身高 = u+A*GT_A+B*GT_B+C*GT_C+D*GT_D+ E*GT_E
#方程1
其中u為群體均值(也就是方程的截距),係數A是A基因座的遺傳效應,GT_A是Aa基因座的基因型,可能是aa、Aa、AA,當然數學上可以使用0,1,2替代。其中,係數A、B、C、D、E都是待求解的變量。
如果求解這個多元線性方程組,我們將發現A、D、E均為0(效應為0),而B、C則顯著大於0,則一樣推斷Bb和Cc基因座對身高是有貢獻的。那麼,它們為什麼對身高有貢獻呢?因為它們與功能基因連鎖啊,由此我們知道了功能基因的初步位置。這就是QTL定位中的線性回歸模型。
以上的方程組在實際情況中,將可能會面臨自變量的數量(標記數量)大於因變量(樣本數),那麼這個方程是不可準確求得唯一解的。所以,通常會將多元線性回歸方程簡化為一元線性回歸方程組。例如,針對Aa基因座,我們可以構建一個方程組如下:
其中,e是隨機誤差效應。那麼在這裡的案例中,方程1就可以拆解為針對5個不同分子標記的方程2,從而一一求解每個標記/區間的效應。因為,這只是個簡單的一元線性回歸方程,求解起來是非常簡單快速的。
這就是在連鎖分析中常用的區間作圖定位法(interval Mapping)的基本原理。
例如下圖,個體A和B有三個QTL位點的差異。假設紅色的基因型相比褐色的基因型可以將個體的高度提高10釐米。現在我想計算標記Marker1的效應,如果我們只考慮單一標記Marker1的效應(使用方程2),最終我們計算的結果以為A個體 30釐米的身高優勢都來自Marker1的差異,就誤把Marker1的效應計為30釐米(高估)。
但如果我們使用多元線性回歸分析,將Marker2和Marker3併入方程組,在方程組中統一考慮它們的效應,那麼對Marker1效應的估算將會更加準確(三個標記效應都是10釐米)。
圖3 目標QTL與背景QTL
在實際情況下,那些被忽略的背景標記效應會對單標記分析模型帶來各種假陽性和假陰性的情況,所以背景標記的效應是必須被考慮的。
但目前的高密度遺傳圖譜,標記數量成百上千個,如上文提到的,如果每個標記效應都被併入方程,那麼使用標準的方法這個方程組是無法求解的(方程1)。所以,在經典的複合區間作圖定位中,採用了一個折中的方式,大體步驟如下:
a)使用單標記回歸以及逐步回歸的方法,從整個基因組中篩選若干個(例如10個)效應最強的標記。
b)在計算某個標記(區間)效應的時候,將其他區域效應最強的那些標記整合入方程組中,如以下方程式:
身高 = u+A*GT_A+[ B*GT_B+… …+ K*GT_K]+ e
# 方程3
在方差3中,未知變量一共有11個(A~K標記的效應),只要個體足夠多這個方程還是可以求解的。其中目標標記是A(我們期望算出它們的效應)。B~K就是基因組其他區域的效應最強的標記,雖然我們暫時並不關心它們具體的效應,但將它們引入方程會讓我們對A的效應估算更準。我們將B~K標記這種並非我們直接關心,但和自變量(A 標記)一樣,同樣對因變量(身高)有影響的標記稱為協變量(covariant)。
所以,目前某些同行公司使用區間作圖模型進行連鎖分析,實際上是不對的。對於數量性狀,只有使用考慮協變量的模型(例如複合區間作圖)才是合理的方法。
LOD值:這個p value的概念略有不同。P value是這個位點不存在QTL的概率。而LOD=log10(L1/L0),其中L1是這個位點有QTL的概率,L0是這個位點無QTL的概率。如果LOD=3,則意味著這個位點有QLT的概率是無QTL的概率的1000倍。
2-LOD置信區間:QTL定位的結果是1個LOD值在染色體上變化的波形圖(如下圖),QTL區域的LOD值會形成一個信號峰。功能基因理論上就位於信號最強(LOD值最大)的峰尖附近。但功能基因通常只是位於這個區間內,而不是必然位於峰尖。離峰尖距離越遠的位置,LOD值不斷下降,功能基因位於該位置的概率越低。
為了便於後續研究中篩選候選基因,我們通常會設置一個範圍篩選候選基因。一般經驗值會使用2-LOD置信區間。這個名詞的意思就是LOD波動曲線從峰的最大值降低2的時候(Y軸),對應在遺傳圖譜上跨越的區域(X軸)。2-LOD置信區間大概對應99.8%的置信區間,即功能基因有99.8%概率已經落在這個區域內了。
1-LOD置信區間也是類似的概念,對應的置信區間大概是97%。
圖4 QTL定位結果示意圖
以上只是我們對QTL定位數學原理的科普性介紹,我們實際上略微簡化了QTL定位中的模型以便利於初學者理解。要了解更多更專業的介紹,點擊「閱讀原文」可從Omicshare論壇相應的主題帖下載專業附件進行深入閱讀。今天的內容就到這裡啦~
封面丨www.plosone.org