群落多样性之Beta多样性（三）

2021-02-08 本文已影响0人布莱特杨

0 导语

倘若想要用数值来评价一个美女的身材，需要考察多方面的数据。
我所能想到的大致上需要考察的变量大致上有这些：身高、体重、头围、胸围、腹围、腰围、臀围、腿围、头颈长、上身长、臂长、腿长、肩宽和胯宽等。因此，考察身材的“颜值”至少得需要十几个变量吧。
首先，如果用这些变量直接估计身材的“颜值”，变量有点多。过多变量会导致需要估计的参数过多，造成维度灾难。
另外，这些变量之间也并不是独立的，有些变量存在着高度的相关性，比如身高和腿长就呈现出正相关。女生想拥有一米以上的大长腿，目测起码身高得一米七以上，有着“香港第一美腿”之称的女神万绮雯据说身高就有一米七一。甚至有时候会存在一些冗余变量，比如很明显身高这个指标可以由头颈长、上身长和腿长的加和求出。
因此，在实际的数据分析过程中，我们需要应用降维的手段去除变量中的冗余成分，避免维度灾难。
降维技术有很多种，包括主成分分析（principal component analysis，PCA）、主坐标分析（principal co-ordinates analysis，PCoA）、非度量多维尺度分析（non-metric multidimensional scaling，NMDS）、因子分析（factor analysis）、独立成分分析（independent component analysis，ICA）等^[1]。目前，PCA的应用最为广泛，因此本文主要关注PCA。

1 PCA简介

PCA是一种常用的降维的方法，其基本原理在群落多样性之Beta多样性（一）中用给半盒红塔山香烟拍照的方式描述过，即数据坐标系转换，而转换后的多维新坐标系方向的选择是由数据本身决定的。假设原始数据存在n个特征属性，即n维。第一个新坐标轴，即第一主成分，选择的是原始数据中方差最大的方向；第二主成分选择与第一主成分正交，且具有最大方差的方向，这个过程将一直重复，直到找出第n个主成分，n即原始数据中属性变量的数量。这样最后我们还是会得到一个n维数据，这么看维数并没有降低。然而，这些主成分坐标是降序排列的，它们对数据特征属性的贡献度呈现下降趋势且往往有着很大的差别，在具体的分析过程中我们常常会观察到仅仅前N（ $N << n$ ）个主成分的累计贡献度就可达到较高的百分比，N一般小于5，现实中为了可视化，往往仅用第1~3主成分为坐标绘制二维或三维图像。例如，图1中展示了应用PCA对不同环境下的微生物群落在门和属水平上的比较，仅在第一主成分（PC1）的贡献度就达到了85%以上^[2]。

图1 应用PCA对不同环境下的微生物群落在门和属水平上的比较

2 理解原始数据

为便于理解，我们就从最简单的2维降到1维开始：构造一个仅仅包含两个操作分类单元（OTU）的微生物群落OTU丰度表（见表1）。

表1 最简单的OTU丰度表

Sample	OTU1	OTU2
Samp1	26	34
Samp2	16	22
Samp3	27	23
Samp4	34	47
Samp5	7	0
Samp6	32	45
Samp7	16	22
Samp8	16	44
Samp9	25	12
Samp10	5	11

表1去掉表头和样本名是一个矩阵 $A$ ，后续的计算主要针对的就是矩阵 $A$ 。这里需要强调一下，理解原始数据对于理解PCA分析很重要。

==================如何理解原始数据==================
当年我上大一学的第一门编程语言是FoxPro，一款坑爹的数据库开发工具。当时老师说过，数据库中行表示的是记录（records），列则为变量（variables）。
在具体的宏基因组学项目中的OTU或基因丰度表（如本文表1）的结构与之类似，可以看成是一个矩阵 $A$ 。
其中矩阵 $A$ 的行表示记录，在群落多样性分析中每一行则表示的是就是样本即群落，行的数量用字母 $m$ 表示；区别于矩阵的行，矩阵的列表示变量，又叫属性或者特征或者维，原始变量的数量用 $n$ 表示，降维后的维数则用 $N$ 表示。
由此可见，矩阵 $A$ 是一个 $m\times n$ 的矩阵。
注意：在实际项目中如果遇到列为样本，行是OTU或基因的情况，则需要将原来的表格进行转置，也就是原表的行转换为列，列转换为行，方可进行进一步的计算。

表1在R软件中构造成矩阵如下所示：

> MyData <- matrix(c(26,16,27,34,7,32,16,16,25,5,34,22,23,47,0,45,22,44,12,11), ncol=2)
> MyData #原始数据表的矩阵形式
      [,1] [,2]
 [1,]   26   34
 [2,]   16   22
 [3,]   27   23
 [4,]   34   47
 [5,]    7    0
 [6,]   32   45
 [7,]   16   22
 [8,]   16   44
 [9,]   25   12
[10,]    5   11

了解了这些，我们就可以开始进行分析了。

3 PCA分析

3.1 数据偏差

计算数据偏差实际上就是各个数据减去相应属性平均值，从而得到各个属性观察值的偏差。

> mean1 <- mean(MyData[,1]) #计算属性的平均值
> mean2 <- mean(MyData[,2])
> mean1 
[1] 20.4
> OTU1 <- MyData[,1] - mean1 #计算偏差
> OTU2 <- MyData[,2] - mean2
> OTU1
 [1]   5.6  -4.4   6.6  13.6 -13.4  11.6  -4.4  -4.4   4.6 -15.4
> NormMyData <- cbind(OTU1, OTU2) #将两列数据合并
> NormMyData
       OTU1 OTU2
 [1,]   5.6    8
 [2,]  -4.4   -4
 [3,]   6.6   -3
 [4,]  13.6   21
 [5,] -13.4  -26
 [6,]  11.6   19
 [7,]  -4.4   -4
 [8,]  -4.4   18
 [9,]   4.6  -14
[10,] -15.4  -15

3.2 协方差矩阵

协方差主要用于衡量两个属性变量变化的一致性。协方差的计算公式如下：
$COV(X,Y)=E[(X-E(X))(Y-E(Y))]=\frac{1}{n-1}\sum_{i=1}^{n} (x_i-x_{mean})(y_i-y_{mean})$
具体来说协方差反映了下面几个问题：
1）协方差的正负符号分别反映了两个变量在变化过程中是同方向变化还是反方向变化？
2）而协方差的数值，即绝对值则反映了同向或反向协同变化的程度。
在R软件中我们可以用函数cov()计算属性变量之间的协方差矩阵。

> MyCOV <- cov(NormMyData)
> MyCOV
          OTU1     OTU2
OTU1  98.93333 111.3333
OTU2 111.33333 258.6667

在群落多样性分析中，协方差反映了两个物种/OTU丰度的协变方向和程度。多个物种/OTU丰度两两之间的的协方差则形成了协方差矩阵。这里需要强调一点，PCA过程中求得的协方差矩阵是属性变量之间的协方差，而不是样本之间的协方差。另外，协方差矩阵也可由相关

3.3 协方差矩阵的特征向量和特征值

特征值和特征向量的分析是PCA分析过程中的重点和难点，也是线性代数中一项重要的内容，能够通过数据的一般格式来揭示数据的“真实”结构。
设 $A$ 为 $n$ 阶方阵，如果数值 $λ$ 和 $n$ 维非零向量x满足等式
$Aν=λν$
成立，则 $λ$ 称为方阵 $A$ 的特征值，非零向量 $v$ 为方阵 $A$ 对应于特征值 $λ$ 的特征向量。
其中， $v$ 是特征向量， $λ$ 是特征值。特征值都是简单的标量值，因此 $Av = λv$ 代表的是：如果特征向量v被某个矩阵A左乘，那么它就等于某个标量 $λ$ 乘以 $v$ 。
R软件的函数eigen()可用于计算特征值和对应的特征向量，代码如下：

> MyEigen <- eigen(MyCOV,symmetric=T)
> MyEigen
eigen() decomposition
$values
[1] 315.8175  41.7825

$vectors
          [,1]       [,2]
[1,] 0.4566761 -0.8896330
[2,] 0.8896330  0.4566761

关于特征值和特征向量的计算，可参见下面的内容，如果已具备相关知识，可以跳过。

==============特征值和特征向量的计算================
计算下面矩阵A的特征值和特征向量
$A= \begin{bmatrix} -2 & 1 & 1 \\ 0 & 2 & 0 \\ -4 & 1 & 3 \end{bmatrix}$
根据特征值和特征向量的定义公式（ $Aν=λν$ ）有：
$(A−λI)ν=0$
则 $A−λI$ 的行列式为0，即
$|A−λ∗I| = \begin{vmatrix} \begin{bmatrix} -2 & 1 & 1 \\ 0 & 2 & 0 \\ -4 & 1 & 3 \end{bmatrix}-\begin{bmatrix} λ & 0 & 0 \\ 0 & λ & 0 \\ 0 & 0 & λ \end{bmatrix} \end{vmatrix}=0$
进一步计算得：
$\quad \ \begin{vmatrix} -2-λ & 1 & 1 \\ 0 & 2-λ & 0 \\ -4 & 1 & 3-λ \end{vmatrix}$
$=-(λ+1)(λ-2)^2=0$
解这个方程得方阵A的特征值 $λ$ 为： $λ_1=-1, λ_2=λ_3=2$ 。
然后分别将各个λ代入方程 $(A−λI)ν=0$ 求各个特征值对应的特征向量，以 $λ_1=-1$ 为例，将其代入 $(A−λI)ν=0$ 得
$\left ( \begin{bmatrix} -2 & 1 & 1 \\ 0 & 2 & 0 \\ -4 & 1 & 3 \end{bmatrix}-\begin{bmatrix} λ & 0 & 0 \\ 0 & λ & 0 \\ 0 & 0 & λ \end{bmatrix}\right )v=\begin{bmatrix} -2-λ & 1 & 1 \\ 0 & 2-λ & 0 \\ -4 & 1 & 3-λ \end{bmatrix}\begin{bmatrix} v_{1,1} \\ v_{1,2} \\ v_{1,3} \end{bmatrix}$
$\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \ =\begin{bmatrix} -1 & 1 & 1 \\ 0 & 3 & 0 \\ -4 & 1 & 4 \end{bmatrix}\begin{bmatrix} v_{1,1} \\ v_{1,2} \\ v_{1,3} \end{bmatrix}=0$
因此得到方程组：
$-1v_{1,1}+1v_{1,2}+1v_{1,3}=0$
$0v_{1,1}+3v_{1,2}+0v_{1,3}=0$
$-4v_{1,1}+1v_{1,2}+4v_{1,3}=0$
可以得出该方程组的一组基础解：
$v_{1,1}=1$
$v_{1,2}=0$
$v_{1,3}=1$
即对应 $λ=-1$ 的一个基础的特征向量为
$v_1=\begin{bmatrix} 1 \\ 0 \\ 1 \end{bmatrix}$
故对应 $λ=-1$ 的所有特征向量，即向量空间为 $kv_1 \quad (k≠0)$ 。
PCA分析需要求单位特征向量，也就是满足欧式范数（ $l_2$ ），即 $||x||_2=\sqrt{x_1^2+x_2^2+...+x_n^2}=1$ 的特征向量，则有
$\begin{Vmatrix} \begin{bmatrix} 1 \\ 0 \\ 1 \end{bmatrix} \end{Vmatrix}=\sqrt{1^2+0^2+1^2}=\sqrt{2}$
其单位特征向量为
$\begin{bmatrix} 1 \\ 0 \\ 1 \end{bmatrix}/\sqrt{2}=\begin{bmatrix} 0.70710678 \\ 0 \\ 0.70710678 \end{bmatrix}$
依此方法继续求出对应另外两个特征值（ $λ_2=λ_3=2$ ）的特征空间和单位特征向量。

3.4 特征值排序和特征向量提取

现在我们知道，每个特征向量都对应一个特征值，最终通过计算得到的特征值的数量和原始数据的属性变量的数量是相等的，都是 $n$ ，也就是说尽管经过了坐标变换但是维数并没有减少。此时，我们需要对特征值的从大到小进行排序，保留特征值最大的N个特征向量（N维），也就是第1~N个主成分，这就是降维，那么最终的维数降了多少呢？即 $n-N$ 。
由于本文的案例是2维降为1维，因此我们只需要提取出对应最大特征值的那个1维的特征向量，即第一主成分。

> MyEigen$vectors[,1]
[1] 0.4566761 0.8896330

3.5 空间变换

我们知道，第1~N主成分分别是从大到小排列为1~N特征值对应的特征向量，它们可以构成一个转换矩阵。得到了主成分特征向量，接下来需要用去平均值后的原始数据右乘转换矩阵进行转换。
$\begin{bmatrix} d_{1,1} & d_{1,2} & ... & d_{1,n} \\ d_{2,1} & d_{2,2} & ... & d_{2,n} \\ ... & ... & ... & ... \\ d_{m,1} & d_{n,2} & ... & d_{m,n}\end{bmatrix} \begin{bmatrix} v_{1,1} & v_{1,2} & ... & v_{1,N} \\ v_{2,1} & v_{2,2} & ... & v_{2,N} \\ ... & ... & ... & ... \\ v_{n,1} & v_{n,2} & ... & v_{n,N}\end{bmatrix}=\begin{bmatrix} p_{1,1} & p_{1,2} & ... & p_{1,N} \\ p_{2,1} & p_{2,2} & ... & p_{2,N} \\ ... & ... & ... & ... \\ p_{n,1} & p_{n,2} & ... & p_{m,N}\end{bmatrix}$
最终所得到转换后的各个样本在第1~N个主成分空间中的坐标。一般的分析选择 $N=2$ ，即将多维数据降维成2维平面。
在本文的案例中，转换矩阵是一个1维的列向量，故此应用去平均值之后的原始数据右乘转换矩阵，得到各个样本在第1主成分1维空间中的坐标（PC1）。

> PC1 <- NormMyData %*% MyEigen$vectors[,1] #原始数据矩阵右乘转换矩阵（这里为特征向量）
> PC1
            [,1]
 [1,]   9.674450
 [2,]  -5.567907
 [3,]   0.345163
 [4,]  24.893089
 [5,] -29.249919
 [6,]  22.200470
 [7,]  -5.567907
 [8,]  14.004020
 [9,] -10.354153
[10,] -20.377307

4 写在后面

前面我们在群落多样性之Beta多样性（二）中介绍了群落间相似性/距离的计算，如果群落只有两个，那么它们之间的距离可以很直观的看出来，就是一个数值。然而，现实的情况是我们经常需要对多个群落进行距离的计算，从而得到一个相似性/距离矩阵，这很不直观，难以观察各个群落的聚集和分散情况。应用PCA分析方法对数据降维可以解决群落距离的展示问题。PCA主要针对属性的协方差进行降维，然而，（二）中涉及到近80种（或许不止）距离的计算方式。因此，应用这些计算方式可计算属性距离矩阵替代协方差矩阵，而各类计算方式有着不同的应用场景，如何选择就要依靠诸位大哥的经验了，当然也可以根据具体情况测试以寻找最合适的距离算法。

附 PCA分析测试的R代码

#PCA R code
MyData <- matrix(c(26,16,27,34,7,32,16,16,25,5,34,22,23,47,0,45,22,44,12,11),ncol=2) #构造测试数据矩阵
#MyData <- matrix(round(runif(20,0,50)),ncol=2) #随机构造测试数据矩阵
OTU1 <- MyData[,1] - mean(MyData[,1])
OTU2 <- MyData[,2] - mean(MyData[,2])
NormMyData <- cbind(OTU1,OTU2)
MyCOV <- cov(NormMyData)
MyEigen <- eigen(MyCOV,symmetric=T)
PC1 <- NormMyData %*% MyEigen$vectors[,1]

参考文献
[1] Ju F, Zhang T. 16S rRNA gene high-throughput sequencing data mining of microbial diversity and interactions [J]. Applied Microbiology and Biotechnology, 2015, 99(10): 4119-4129.
[2] Cui J, et al. Metagenomic insights into a cellulose-rich niche reveal microbial cooperation in cellulose degradation[J]. Frontiers in Microbiology, 2019, 10(1): 618.
[3] https://blog.csdn.net/HLBoy_happy/article/details/77146012
[4] https://blog.csdn.net/m0_37788308/article/details/78115209
[5] https://wenku.baidu.com/view/f14c18215901020207409c97.html
[6] Gilbert strang. Linear Algebra and Its Applications 4th Edition[M]. May 9, 2011.