集团遗传学: 主成分判别分析(DAPC)

通常我们会遇到在没有先前知识的前提下通过确定观察到的群体（集团）数量来推断人口结构。这种情况可以使用几种方法来推断群体，例如K均值聚类、使用STRUCTURE的贝叶斯聚类和多变量方法，如主成分判别分析（DAPC）（Pritchard，Stephens＆Donnelly，2000; Jombart，Devillard＆Balloux，2010; Grünwald和Goss，2011）。类STRUCTURE的方法必须基于标记不相关的假定下，并且人口是全自交的（Pritchard等，2000）。对于部分独立的种群，基于遗传距离的K均值聚类或DAPC是更方便的无模型方法。本文着重于探讨DAPC方法。

流感病毒株H3N2的DAPC分析

DAPC（主成分判别分析）由Jombart等人（Jombart et al., 2010）首创，可以用于推断具有遗传关联的个体的群集数量。在这种多元统计方法中，样本的方差被分为组间和组内两个部分，以最大程度地区分不同群体。在DAPC中，首先使用主成分分析（PCA）对数据进行转换，然后使用判别分析（DA）识别群集。这个教程基于Thibaut Jombart编写的vignette（代码示例。用户还可以通过在R中执行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()函数的参数解释如下：

1. dataset：这个参数指的是正在分析的数据集，本例中是H3N2数据集。它包含了H3N2流感株的遗传数据。
2. var.contrib：这个参数被设置为TRUE，表示我们希望在输出中保留对分析有贡献的变量信息（位点信息）。通过将其设置为TRUE，您可以后续查看哪些位点贡献了群体的分离。
3. center：这个参数被设置为FALSE，表示我们不希望对数据进行重新缩放，使均值等于0。数据居中是一种预处理步骤，涉及将每个变量减去均值，这在某些分析中有帮助。但是，在这种情况下，不需要进行居中处理。
4. n.pca：这个参数确定在主成分分析（PCA）步骤中保留的轴（主成分）的数量。如果设置为NULL（默认值），函数将根据数据自动确定适当的主成分数量。主成分捕捉到了遗传变异的主要方向。
5. n.da：这个参数确定在判别分析（DA）步骤中保留的轴的数量。如果设置为NULL（默认值），函数将根据数据自动确定适当的判别轴的数量。判别轴是从主成分派生出来的，用于识别数据中的群集或群体。

这些参数提供了灵活性和对分析的控制，允许自定义在PCA和DA步骤中使用的成分数量，并指定是否包括变量贡献和居中处理。

可以看出2001到2005的数据根据横坐标一字排开，2006的数据和2001-2005的数据根据在纵坐标上分的很开。

接下来可以看一下，2006样本是不是存在和2001-2005差异很大的等位基因点。

set.seed(4)
contrib <- loadingplot(dapc.H3N2$var.contr, axis = 2, thres = 0.07, lab.jitter = 1)

399和906的SNPs很扎眼。所以接下来可以着重查看这两个点的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)

这个图很好地说明了季节性H3N2流感病毒中突变后随之而来的选择或漂变的影响。