通常我們會(huì)遇到在沒有先前知識(shí)的前提下通過確定觀察到的群體（集團(tuán)）數(shù)量來推斷人口結(jié)構(gòu)抄淑。這種情況可以使用幾種方法來推斷群體，例如K均值聚類挥等、使用STRUCTURE的貝葉斯聚類和多變量方法捕仔，如主成分判別分析（DAPC）（Pritchard键畴，Stephens＆Donnelly，2000; Jombart驻右，Devillard＆Balloux什黑，2010; Grünwald和Goss，2011）堪夭。類STRUCTURE的方法必須基于標(biāo)記不相關(guān)的假定下兑凿，并且人口是全自交的（Pritchard等，2000）茵瘾。對(duì)于部分獨(dú)立的種群礼华，基于遺傳距離的K均值聚類或DAPC是更方便的無模型方法。本文著重于探討DAPC方法拗秘。
流感病毒株H3N2的DAPC分析
DAPC（主成分判別分析）由Jombart等人（Jombart et al., 2010）首創(chuàng)圣絮，可以用于推斷具有遺傳關(guān)聯(lián)的個(gè)體的群集數(shù)量。在這種多元統(tǒng)計(jì)方法中雕旨，樣本的方差被分為組間和組內(nèi)兩個(gè)部分扮匠，以最大程度地區(qū)分不同群體捧请。在DAPC中，首先使用主成分分析（PCA）對(duì)數(shù)據(jù)進(jìn)行轉(zhuǎn)換棒搜，然后使用判別分析（DA）識(shí)別群集疹蛉。這個(gè)教程基于Thibaut Jombart編寫的vignette（代碼示例。用戶還可以通過在R中執(zhí)行adegenetTutorial("dapc")來打開vignette力麸。
# DAPC requires the adegenet package. Let's load this package:
library("adegenet")
data(H3N2) # load the H3N2 influenza data. Type ?H3N2 for more info.
pop(H3N2) <- H3N2$other$epid
dapc.H3N2 <- dapc(H3N2, var.contrib = TRUE, scale = FALSE, n.pca = 30, n.da = nPop(H3N2) - 1)
scatter(dapc.H3N2, cell = 0, pch = 18:23, cstar = 0, mstree = TRUE, lwd = 2, lty = 2)
在上文中使用的dapc()函數(shù)的參數(shù)解釋如下：
- dataset：這個(gè)參數(shù)指的是正在分析的數(shù)據(jù)集可款，本例中是H3N2數(shù)據(jù)集。它包含了H3N2流感株的遺傳數(shù)據(jù)克蚂。
- var.contrib：這個(gè)參數(shù)被設(shè)置為TRUE闺鲸，表示我們希望在輸出中保留對(duì)分析有貢獻(xiàn)的變量信息（位點(diǎn)信息）。通過將其設(shè)置為TRUE埃叭，您可以后續(xù)查看哪些位點(diǎn)貢獻(xiàn)了群體的分離摸恍。
- center：這個(gè)參數(shù)被設(shè)置為FALSE，表示我們不希望對(duì)數(shù)據(jù)進(jìn)行重新縮放赤屋，使均值等于0立镶。數(shù)據(jù)居中是一種預(yù)處理步驟，涉及將每個(gè)變量減去均值类早，這在某些分析中有幫助谜慌。但是，在這種情況下莺奔，不需要進(jìn)行居中處理欣范。
- n.pca：這個(gè)參數(shù)確定在主成分分析（PCA）步驟中保留的軸（主成分）的數(shù)量。如果設(shè)置為NULL（默認(rèn)值）令哟，函數(shù)將根據(jù)數(shù)據(jù)自動(dòng)確定適當(dāng)?shù)闹鞒煞謹(jǐn)?shù)量恼琼。主成分捕捉到了遺傳變異的主要方向。
- n.da：這個(gè)參數(shù)確定在判別分析（DA）步驟中保留的軸的數(shù)量屏富。如果設(shè)置為NULL（默認(rèn)值）晴竞，函數(shù)將根據(jù)數(shù)據(jù)自動(dòng)確定適當(dāng)?shù)呐袆e軸的數(shù)量。判別軸是從主成分派生出來的狠半，用于識(shí)別數(shù)據(jù)中的群集或群體噩死。
這些參數(shù)提供了靈活性和對(duì)分析的控制，允許自定義在PCA和DA步驟中使用的成分?jǐn)?shù)量神年，并指定是否包括變量貢獻(xiàn)和居中處理已维。
可以看出2001到2005的數(shù)據(jù)根據(jù)橫坐標(biāo)一字排開，2006的數(shù)據(jù)和2001-2005的數(shù)據(jù)根據(jù)在縱坐標(biāo)上分的很開已日。
接下來可以看一下垛耳，2006樣本是不是存在和2001-2005差異很大的等位基因點(diǎn)。
set.seed(4)
contrib <- loadingplot(dapc.H3N2$var.contr, axis = 2, thres = 0.07, lab.jitter = 1)
399和906的SNPs很扎眼。所以接下來可以著重查看這兩個(gè)點(diǎn)的SNPs的變化情況堂鲜。
temp <- seploc(H3N2) # seploc {adegenet} creates a list of individual loci.
snp906 <- tab(temp[["906"]]) # tab {adegenet} returns a matrix of genotypes
snp399 <- tab(temp[["399"]])
# The following two commands find the average allele frequencies per population
(freq906 <- apply(snp906, 2, function(e) tapply(e, pop(H3N2), mean, na.rm = TRUE)))
## 906.c 906.t
## 2001 0.000000000 1.0000000
## 2002 0.000000000 1.0000000
## 2003 0.000000000 1.0000000
## 2004 0.000000000 1.0000000
## 2005 0.002155172 0.9978448
## 2006 0.616071429 0.3839286
(freq399 <- apply(snp399, 2, function(e) tapply(e, pop(H3N2), mean, na.rm = TRUE)))
## 399.c 399.t
## 2001 0.000000000 1.0000000
## 2002 0.000000000 1.0000000
## 2003 0.000000000 1.0000000
## 2004 0.001848429 0.9981516
## 2005 0.002079002 0.9979210
## 2006 0.357142857 0.6428571
# First, set the plotting parameters
# mfrow = number of columns, rows
# mar = plot margin size
# las = axis label style (3: always vertical)
par(mfrow = c(1, 2), mar = c(5, 4, 4, 0) + 0.1, las = 3)
matplot(freq906, pch = c("c", "t"), type = "b",
xlab = "year", ylab = "allele frequency", main = "SNP # 906",
xaxt = "n", cex = 1.5)
axis(side = 1, at = 1:6, lab = 2001:2006)
matplot(freq399, pch = c("c", "t"), type = "b",
xlab = "year", ylab = "allele frequency", main = "SNP #399",
xaxt = "n", cex = 1.5)
axis(side = 1, at = 1:6, lab = 2001:2006)
這個(gè)圖很好地說明了季節(jié)性H3N2流感病毒中突變后隨之而來的選擇或漂變的影響栈雳。
