plink R语言群体间遗传分析 r语言系统聚类

转载

mob6454cc7b19b2 2024-08-22 07:41:38

文章标签 plink R语言群体间遗传分析聚类 python 数据分析 js 文章分类 R语言后端开发

聚类分析就是根据样本变量特征的相似程度将样本分成若干类，每类称为一个簇，一般要求簇内差异最小化，簇间差异最大化。本篇介绍聚类分析方法中的系统聚类法。系统聚类法，也称层次聚类法（Hierarchical Cluster Method）、谱系聚类法。

假设样本数为n，簇数为k。在初始状态下，系统聚类法将每个样本看作一簇，即k = n；然后将最相似的两个样本合并为一簇，得到k = n-1的结果；再在剩余的n-1簇中再次寻找最相似的两簇合并为一簇，得到k = n-2的结果；以此类推，直至k = 1。在这个过程中形成了聚类树。

为了方便演示，本篇从mtcars数据集中随机挑选7个样本进行分析：

library(tidyverse)
set.seed(20210527)
data <- slice_sample(mtcars, n = 7)
rownames(data) <- LETTERS[1:7]

data
##    mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## A 17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
## B 21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
## C 22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
## D 19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
## E 17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
## F 19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
## G 21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1

1 距离函数

系统聚类法在初始状态下依赖于样本间的“距离”矩阵来判断样本间的相似度。生成距离矩阵的的函数是dist：

dist(x, method = "euclidean", 
     diag = FALSE, upper = FALSE, p = 2)

x：包含样本信息的矩阵或数据框；
method：距离计算方法，默认值为euclidean（欧几里得距离），可选项有maximum（切比雪夫距离）、manhattan（曼哈顿距离）、canberra（堪培拉距离）、binary（哑变量距离）、minkowski（闵可夫斯基距离）。

以二维变量空间为例：
闵可夫斯基距离；当时即曼哈顿距离；当时即欧几里得距离；
切比雪夫距离；
堪培拉距离；
binary方法用于计算0-1变量之间的“距离”，两个样本的距离等于取值不相等的变量个数占取值不全为0的变量个数的比例。

diag：取值为TRUE时输出结果为矩阵；
upper：diag取值为FALSE时输出结果为距离矩阵的下三角，若upper参数为TRUE，则输出结果为距离矩阵的上三角。

在进行系统聚类时，的diag和upper参数保持默认值即可。

p：闵可夫斯基距离的p参数。

为了避免变量的量纲和绝对大小对聚类结果的影响，在进行聚类分析前，一般会对原始数据进行标准化，标准化方法可根据需要选取。

# 数据标准化
data0 <- scale(data)
# 计算距离
D = dist(data0, method = "minkowski", p = 1.5)
print(D)

##           A         B         C         D         E         F
## B  6.065828                                                  
## C 10.256391  6.328653                                        
## D  6.155329  5.288102  6.140178                              
## E  6.182098  5.407913  6.526775  1.007835                    
## F  2.238748  6.022026 10.351107  6.991633  7.400449          
## G  8.650247  3.572940  4.448209  5.777199  5.831998  8.897885

2 聚类函数

在k < n后的过程中，每次聚类都会选择簇间差异最小的两簇进行合并。定义簇间“距离”的方法有如下几种：

完全连接法（complete linkage）

簇间距离最远的两个样本的距离。

单连接法（single linkage）

簇间距离最近的两个样本的距离。

平均连接法（average linkage）

簇间所有样本对的距离的平均值。

中间连接法（median linkage）

簇间所有样本对的距离的中位数。

重心法（centroid method）

簇间重心的距离。

离差平方和法（ward method）

该方法由Ward(1963)^[1]提出，也称Ward法，两簇合并后使得簇内总离差平方和的增加量最小化。此方法下的簇间距离等于簇间重心的距离再乘以，其中和

聚类函数为hclust：

hclust(d, method = "complete", members = NULL)

d：dist函数的输出对象；
method：簇间“距离”的计算方式，默认值为complete，可选项有ward.D、ward.D2、single、average、mcquitty、median、centroid。

ward.D为R语言中最初用来实现ward法的参数取值，但后来一篇论文（Murtagh, Fionn and Legendre, Pierre (2014)^[2]）发现该方法实现的Ward法是错误的，所以就又有了ward.D2；
mcquitty法是由McQuitty(1966)^[3]提出的一种相似方法。

tree <- hclust(D, method = "ward.D2")
names(tree)

## [1] "merge"       "height"      "order"       "labels"      "method"     
## [6] "call"        "dist.method"

length(tree) == dim(data)[1]

## [1] TRUE

使用plot函数可以对聚类树进行可视化：

plot(x, labels = NULL, hang = 0.1, check = TRUE,
     axes = TRUE, frame.plot = FALSE, ann = TRUE,
     main = "Cluster Dendrogram",
     sub = NULL, xlab = NULL, ylab = "Height", ...)

hang参数控制样本标签的悬挂的起始位置；为非负数时数值越大标签起始位置越靠下；为任意负数时，标签起始位置都在纵轴为0处。

plot(tree)

plink R语言群体间遗传分析 r语言系统聚类_python

# hang为任意负数
plot(tree, hang = -1)

plink R语言群体间遗传分析 r语言系统聚类_数据分析_02

3 分割树函数

聚类树是一种二叉树，使用水平线可以将其分成若干份，簇数等于水平线与聚类树的交点个数，如图所示：

plink R语言群体间遗传分析 r语言系统聚类_数据分析_03

cutree函数可以对树形结构进行分割：

cutree(tree, k = NULL, h = NULL)

该函数有两种分割方式：
通过k参数指定分割后的簇数；
通过h参数指定分割水平线的高度，即对应的纵轴值。

cutree(tree, k = 3)

## A B C D E F G 
## 1 1 2 2 2 3 2

cluster <- cutree(tree, k = c(2,4))
table(a = cluster[,"2"],  # 交叉统计表
      b = cluster[,"4"])
      
##    b
## a   1 2 3 4
##   1 2 0 0 0
##   2 0 2 1 2

cutree(tree, h = 6)

## A B C D E F G 
## 1 2 2 3 3 1 2

对完整的mtcars数据集进行系统聚类：

library(magrittr)

tree0 <- mtcars %>%
  scale() %>%
  dist(method = "euclidean") %>%
  hclust(method = "ward.D2") %T>%
  plot(hang = -1)

plink R语言群体间遗传分析 r语言系统聚类_聚类_04

cutree(tree0, k = 3)

##           Mazda RX4       Mazda RX4 Wag          Datsun 710      Hornet 4 Drive 
##                   1                   1                   2                   2 
##   Hornet Sportabout             Valiant          Duster 360           Merc 240D 
##                   3                   2                   3                   2 
##            Merc 230            Merc 280           Merc 280C          Merc 450SE 
##                   2                   2                   2                   3 
##          Merc 450SL         Merc 450SLC  Cadillac Fleetwood Lincoln Continental 
##                   3                   3                   3                   3 
##   Chrysler Imperial            Fiat 128         Honda Civic      Toyota Corolla 
##                   3                   2                   2                   2 
##       Toyota Corona    Dodge Challenger         AMC Javelin          Camaro Z28 
##                   2                   3                   3                   3 
##    Pontiac Firebird           Fiat X1-9       Porsche 914-2        Lotus Europa 
##                   3                   2                   2                   2 
##      Ford Pantera L        Ferrari Dino       Maserati Bora          Volvo 142E 
##                   1                   1                   1                   2