群體結構圖形——structure堆疊圖
1、structure圖的由來
圖1 假設群體亞群數等於3（k=3）的情況下的structure分析結果
“Structure圖”名詞本身來自這種圖形的分析軟件——STRUCTURE。這個軟件是由斯坦福大學Pritchard實驗室開發的一款群體結構分析軟件，最早在2000年發表在《Genetics》上[1]。
圖2 structure驚人的引用次數
Structure軟件分析達到的目的其實也很簡單——分析群體的結構。要獲知群體的結構，其實在前兩篇與大家分享的兩個圖形也是可以解決的——PCA和樹。通過樹形圖或PCA散點圖，我們也可以直觀了解個體間的分類關系（如圖3）。但如果需要解答：這個群體應該分為3個亞群還是4個亞群更合理，群體間是否存在基因交流（直白的說法：是否有雜交），以及每個個體混血程度是多少？對於PCA圖和樹形圖，這解決不了，那么這時候我們就用到了structure圖。
圖3 通過進化樹分析群體的分類關系
2、structure圖形的解讀
Structure圖在開發之初，就使用了與進化樹或PCA截然不同的算法，所以能夠解決后兩者不能解決的問題。具體的算法思路我在下文會解釋，我們先直接解讀structure圖形的意義。當然，我們也不能完全不懂structure的基本原理就直接解讀structure圖。
Structure分析的基本思路過程是：
1、獲得一系列樣本的基因型；
2、在大部分情況下，我並不知道這個群體實際包含幾個亞群。我們把群體的亞群數稱為K值。那么我可以預設群體亞群數等於1~n，即K=1~n。那么軟件就會模擬在K=x的情況下，使用貝葉斯算法推算群體是如何分群的，以及每個個體的血統構成是如何的。
例如圖4，就是在蠶重測序文章中[2]假設群體結構數的先驗值K=2~4的模擬結果。以K=3的模擬結果為例。這個圖形本質上就是1個堆疊圖，每個直方柱子就是代表1個樣本（這里是40個品種的蠶），柱子的顏色以及顏色比例代表這個樣本屬於哪個亞群以及血統構成比例是多少。我們很直觀的通過顏色就可以判斷這40個個體是如何划分為3個亞群體（紅、黃、綠）。
而且，我們可以看到某些樣本的血統是很純正的（1種顏色），而某些樣本卻由2~3個顏色組成。例如放大圖中的樣本D06為例，它對應的直方柱子由兩種顏色構成，大概是40%的黃色和60%的紅色。說明這個個體很有可能是從兩個祖先亞群雜交而來，且兩個祖先亞群各占了40%和60%的血統。
以此類推，我們就可以解讀在k=2、3、4的情況下，這些樣本是如何分群的，以及每個樣本的血統構成是多少。
圖4 家蠶文章中在k=2~4的情況下，structure分析的結果
先驗假設K=2~4，而且每個K值假設都對應1個圖形。到底哪個假設是最正確呢？
因為structure是基於貝葉斯模型的計算方法，對每個K值模擬的結果，都會對應產生最大似然值（likelihood）。這個軟件的最大似然值是取了自然對數后輸出的（ln likelihood）。ln likelihood越大，說明K值越接近於真實情況（記住ln likelihood是小於0的，所以我說的這個值越大越好，就是絕對值越小越好）。當然，一般隨着K值升高，ln likelihood值也會不斷升高，但會慢慢進入平台期。選擇最優K值的目標是要找到那個拐點。
簡單說來，就是要找的一個likelihood最大（越大越可靠）而且K值最小（亞群數最少）的模擬結果，往往這樣的模擬對應的K值是最接近於群體的真實情況的。具體的計算過程可以參考2006年的一篇參考文獻[3]，那篇文章解釋非常清楚。
圖5 通過ln likelihood預估最佳K值的過程
3、structure的計算原理
Structure是與PCA、進化樹相似的方法，就是利用分子標記的基因型信息對一組樣本進行分類，分子標記可以是SNP、indel、SSR… …當然，對於重測序應用的最多的還是SNP。當然，structure本質上使用了與PCA、進化樹完全不同的思路。進化樹和PCA本質上都是計算樣本序列間的差異程度，然后利用兩兩差異度聚類（進化樹）或降維（PCA）來實現對樣本的分類。
但PCA和樹的不足，如上文提到的，無法推算：
1、群體應該划分為幾個亞群；
2、群體間基因交流的程度；
3、某個個體的混血程度。
但structure軟件摒棄的以上的方法，先預設群體由若干亞群（k=x）構成，通過模擬算法找出在k=x的情況下，最合理的樣本分類方法。最后再根據每次模擬的最大似然值，找出最適用這群體的K值。
這個過程的邏輯如下：
1.亞群內符合哈溫平衡
那么，軟件在如何確定樣本的最優分類方法呢？其實基於一個假設：在各個亞群內部個體應該符合哈代-溫伯格平衡（哈溫平衡的概念可以在百度查詢），那么這個亞群內的基因頻率分布應該可通過哈溫平衡檢驗。例如，現在有40個個體的1個SNP位點的基因型，我預設亞群數k=2。我先隨機將40個個體分成兩份，然后檢驗是否符合哈溫平衡。如果不符合，我繼續調整分類策略，直到找到一種最優的分類方法：40個個體被分為了兩份，每個亞群都由若干個體構成，每個亞群內部都最大程度地符合哈溫平衡。
2.每個位點是獨立的
同一個體基因組上的不同SNP可能來源不同亞群體，這是由於雜交混血過程帶來的效應。如下圖的D06個體，就各有部分DNA來自紅色和黃色的亞群體。從另一個維度理解，為了達到哈溫平衡，對不同的位點的分類方法是不同的。例如下圖中，位點1的分類和位點2的分類策略就不同。位點1將D06個體划為紅色亞群，而位點2將D06個體划分為黃色亞群。也就是說，軟件是對每個位點單獨進行分群的。
圖6 structure分析要求位點是獨立的
3.每個樣本的血統構成
既然對每個位點都完成分群了，自然最終就可以計算每個個體的血統構成了。
如果大家明白了這三個步驟了，我再以k=2為例，解釋一下structure是如何找到樣本的最優分類。其實簡單說來，就是利用了計算機超強的運行能力，一開始計算機只是隨機將樣本分為兩份，然后在每個亞群內進行哈溫平衡檢驗。如果不符合哈溫平衡（拍腦袋的分類，一開始當然是慘不忍睹），計算機繼續調整分類，然后繼續檢驗。
如此這般，在計算n次后，計算機再從這一堆結果中找到最佳的分類。這個過程稱為“隱馬科夫-蒙特卡羅鏈”的過程，計算次數n就是這個鏈的長度，這是structure一個重要的參數“Number of MCMC Reps”，需要預先設定。
但因為這個計算的過程是從隨機模擬開始的。如果一開始拍腦袋拍的不好（隨機分類與真實分類差距太大），計算機一黑到底，最后把n次用完了，都沒有找到一個合理的分類。所以，分析軟件往往有個預實驗的過程。
就是在正式進行大規模運算前，計算機先嘗試各種各樣的隨機分類，運行非常短的次數，然后評估哪種隨機分類是最合理的。之后，在根據最優的隨機分類，進行后續的大規模運算。這個過程就稱為burn-in period，預實驗的次數就稱為burin-in的次數。這也是structure分析另外一個重要的參數“length of burn-in period”。
如果你理解了這兩個參數ok，非常好。Structure軟件最重要的兩個參數你已經掌握了。剩下就選擇使用那種模型了。主要涉及兩種模型 no admixture model和admixture model。前者假設亞群間不存在雜交，后者則假設亞群間存在雜交。在絕大部分情況下，當然是選擇admixture 模型更合理了。
4、structure分析的輸入以及軟件
Structure分析，輸入的數據就是樣本的基因型數據，注意：輸入的必須是不存在連鎖不平衡的獨立位點。所以，對重測序的結果來說，輸入所有的SNP是不對的。一來，輸入的SNP數量太大，會大大拖長軟件運行的時間；其次，如果使用大量存在連鎖不平衡的位點，就違背了這個軟件最初的假設——各個位點是獨立的。正確的做法是，根據連鎖不平衡衰減分析的結果，僅從所有SNP中挑選一部分獨立的位點用於structure分析。
Structure分析當然最經典的軟件就是STRUCTURE。但Structure分析還有其他軟件可以選擇[4]。比較經典的軟件包括：ADMIXTURE、FRAPPE。這兩個軟件的運行速度都大大超過STRUCTURE。但FRAPPE的不足沒有提供方法估算最佳K值。ADMIXTURE使用與STRUCTURE相同的模型，而且運行效率也很好，所以是一個比較推薦的軟件。
這些軟件的分析結果都只是一份表格，就是每個樣本的血統構成比例。要把表格變成圖形的話，還有繪圖的步驟。最簡單的畫法，你可以使用excel將這個結果繪制為堆疊圖，或者也可以使用其他專門的圖形化軟件，例如Distruct就是比較推薦的一款將structure分析結果圖形化的軟件。
圖7 structure分析的最初結果只是一張表
Structure軟件是有windows的Java圖形界面版本的。參數就剛剛給大家介紹的lengthof burn-in period、Number of MCMC Reps after burn-in、Admixture model。這個三個參數的設置，可以參考文獻3給出的建議。兩個參數各取10,000。不過隨着樣本數量的增加，需要適當提高參數，會讓結果更可靠。
參考文獻：
[1]Pritchard J K, Stephens M, Donnelly P.Inference of population structure using multilocus genotype data[J]. Genetics,2000, 155(2): 945-959.
[2]Xia Q, Guo Y, Zhang Z, et al. Completeresequencing of 40 genomes reveals domestication events and genes in silkworm(Bombyx)[J]. Science, 2009, 326(5951): 433-436.
[3]Evanno G, Regnaut S, Goudet J. Detecting thenumber of clusters of individuals using the software STRUCTURE: a simulationstudy[J]. Molecular ecology, 2005, 14(8): 2611-2620.
[4] Porras-Hurtado L, Ruiz Y, Santos C, et al. An overview of STRUCTURE:applications, parameter settings, and supporting software[J]. Front Genet,2013, 4: 98.