StatQuest学习笔记15——MDS

前言——主要内容

这篇笔记是StatQuest系列视频教程的第44和第45节，其中，第44讲的是多维尺度变换（multidimensional scaling，MDS）与主坐标分析（principal coordinate analysis,PCoA）。第45节讲的是用R来计算MDS和PCoA。

MDS

MDS的分类

MDS与PCA非常相似，都是一种降维的方法。MDS分为度量MDS（Metric Multidimensional scaling）和非度量MDS（Non-metric Multidimensional scaling），其中非度量MDS也称为经典（Classical）MDS，也叫做主坐标分析（principal coordinate analysis,PCoA），如下所示：

MDS的基本思想

为了更方便地理解MDS，我们还是像介绍PCA那些，先讲一个比较简单的例子，例如，我们在下面的图形中，我们先看一群绿色的点，把它们当成是一群正常的细胞（当然了，如果你不是做生物的，你也可以把它们视为一群人，或一批汽车，或者城市等），如下所示：

虽然这群细胞看起来比较像，但实际上，它们有可能不同，如下所示：

当然了，还有可能像这个样子，如下所示：

但问题就在于，仅从直观的角度来看，我们不太容易区分它们，但我们可以通过类似RNA-Seq这样的手段来发现这些细胞有哪些基因是活跃的，这样就可以理解这些细胞的具体功能，如果这是一群人，我们可以测量他们的血压，身高等数据，如下所示：

我们就当它们是细胞，这是测完RNA-seq的一些数据，如下所示：

其中，第1行表示的是样本，第1列是基因名，其余的部分是不同样本某个基因的表达水平，此时，如果根据我们前面介绍的PCA法对这些数据进行分析，那么我们就能把不同样本之间的相互关系转换到一个二维的坐标图上，如下所示：

那些高度相关的样本就会聚集成一块，如下所示：

此时，我们再回到MDS和PCoA上来，这两种方法非常像PCA，只是它们不是把样本之间的相互关系转换成二维坐标，而是把不同样本之间的距离转换成二维坐标，如下所示：

为了进行MDS和PCoA分析，我们可以计算两两样本之间的距离，就像下面的这个样子，如下所示：

为了方便理解，我们此时就只计算一下Cell1和Cell2的距离，如下所示：

欧氏距离

一种常见的计算两个样本之间距离的方法是欧氏距离(Euclidean Distance) 计算法，如下所示：

如果此时我们只有两个基因，就像下面的这个样子，如下所示：

为了计算这两个样本之间的欧氏距离，我们可以用线段来表示它们基因之间的距离，如下所示：

然后通过勾股定理，计算出欧氏距离，如下所示：

如果我们有更多的基因，那么就按照这种方法，求出两个样本中，相同基因的差值，再求其平方和的平方根，像下面的这个样子：

一旦我们计算好了不同样本之间的欧氏距离，那么我们就可以把这些样本绘制到一个二维坐标图上，如下所示：

利用欧氏距离（Euclidean Distance）来进行的MDS和PCoA分析的一个局限就是，最终生成的二维坐标轴基本上与PCA生成的二维坐标图是一样的，如下所示：

换句话讲，这种基于最小线性距离的聚类（clustering）与最大化线性相关是一样的，如下所示：

不过，计算两两样本之间的距离并非只有欧氏距离（Euclidean Distance）一种方法，如下所示：

其他距离方法

例如，另外一种计算两两样本之间的距离方法就是计算基因间差值的log转换的平均值，例如在下面的案例中，对于Cell1和Cell2的Gene1，它经过log转换后的数值就是log(3/0.25)，同理，Gene2和Gene 8的计算也同样如此，如下所示：

经过计算，再取它们的绝对值，最终我们计算出这些数值的平均值，如下所示：

最终我们看一下分别使用了欧氏距离的MDS图（左侧）和使用了差异倍数log转换的MDS图的坐标轴（右侧），如下所示：

不过，在生物学家眼中，通常会选择表达差异倍数的log转换值来计算两两样本之间的距离，因为他们研究的基因通常是成倍变化的，这样用表达差异倍数的log转换会比较合适，如下所示：

我们常见的计算样本间距离的方法有：曼哈顿距离（Manhattan Distance）、汉明距离（Hamming Distance）、大圆距离（Great Circle Distance）等，如下所示：

总结

PCA构建二维坐标图的方法基础是不同样本之间的相关性，而MDS和PCoA构建二维坐标图的方法基础则是基于不同样本之间的距离，如下所示：

PCA的计算流程为：计算样本间的相关性 -> 复杂数学运算（主要是特征值分解） -> 绘图、主成分变异，Loading Scores。

而MDS和PCoA的计算流程与PCA的计算流程的不同之处就在于，第一步是计算样本间的距离，其余的部分与PCA的计算差异不大。

R与MDS/PCoA

生成模拟数据集

在使用R进行MDS/PCoA分析之前，我们还需要使用R来生成一些随机数据，如下所示：

library(ggplot2)
data.matrix <- matrix(nrow=100, ncol=10)
# generate some fake data
# This dataset consists of a matrix with 10 columns, corresponding to 10 samples
# and 100 rows, corresponding to measurements from 100 genes

colnames(data.matrix) <- c(
  paste("wt", 1:5, sep=""),
  paste("ko", 1:5, sep=""))

rownames(data.matrix) <- paste("gene", 1:100, sep="")

for (i in 1:100){
  wt.values <- rpois(5, lambda = sample(x=10:1000, size = 1))
  ko.values <- rpois(5, lambda = sample(x=10:1000, size = 1))

  data.matrix[i,] <- c(wt.values, ko.values)
}

最终生成的数据如下所示：

> head(data.matrix)
      wt1 wt2 wt3 wt4 wt5 ko1 ko2 ko3 ko4 ko5
gene1 257 263 246 246 233 993 945 961 931 958
gene2 822 888 866 845 855 232 210 217 209 225
gene3 122 127 139 156 140 370 336 366 371 370
gene4 812 806 874 860 896 918 934 963 939 982
gene5 925 870 845 931 921 583 582 547 539 546
gene6 646 611 644 676 644 535 590 530 537 551

PCA分析

为了比较MDS/PCoA分析和PCA分析的区别，在这里，我们先进行一下PCA的分析，代码跟上一篇笔记的PCA分析中的代码类似，如下所示：

pca <- prcomp(t(data.matrix), scale = TRUE, center=TRUE)
pca.var <- pca$sdev^2
pca.var.per <- round(pca.var/sum(pca.var)*100, 1)
pca.var.per
pca.data <- data.frame(Sample=rownames(pca$x),
                       X=pca$x[,1],
                       Y=pca$x[,2])
pca.data

结果如下所示：

> pca.var.per
 [1] 87.8  3.3  2.6  1.5  1.4  1.2  0.8  0.8  0.6  0.0
> pca.data <- data.frame(Sample=rownames(pca$x),
+                        X=pca$x[,1],
+                        Y=pca$x[,2])
> pca.data
    Sample         X          Y
wt1    wt1 -8.698916  1.8065989
wt2    wt2 -8.724745 -2.5129243
wt3    wt3 -8.848156  0.1415564
wt4    wt4 -9.161110  1.2735904
wt5    wt5 -8.993328 -0.7037979
ko1    ko1  9.387947  0.9745217
ko2    ko2  9.093852 -1.2375433
ko3    ko3  8.764754  1.6803429
ko4    ko4  8.584933  1.6705990
ko5    ko5  8.594769 -3.0929438

绘制PCA图，如下所示：

ggplot(data=pca.data, aes(x=X, y=Y, label=Sample))+
  geom_text() +
  xlab(paste("PC1 - ", pca.var.per[1], "%", sep=""))+
  ylab(paste("PC2 - ", pca.var.per[2], "%", sep=""))+
  theme_bw()+
  ggtitle("My PCA Graph")

最终图形如下所示：

其中PC1占据了总变异的87.7%，PC2占据了总变异的3.3%。另外，在图上我们可以看出，wt组都位于左侧，ko组都位于右侧，它们主要是被PC1区分了开来。

MDS/PCoA分析

第一步，构建距离距离，使用dist()函数，如下所示：

distance.matrix <- dist(scale(t(data.matrix), center=TRUE, scale=TRUE),
                        method="euclidean")
# parameter of method has 6 types different distnace to choose

第二步：在距离矩阵的基础上进行多维尺度变换，这里使用cmdscale()函数，这个函数的名称比较奇怪，其实它是Classical Muti-Dimensional Scaling的缩写，代码如下所示：

mds.stuff <- cmdscale(distance.matrix, eig=TRUE, x.ret=TRUE)
# eig=TRUE willl return eigen values

其中使用eig=TRUE会返回特征值，用这些特征来计算最终MDS图中的每个坐标占距离矩阵总变异的比例。x.ret=TRUE用于表示是否返回双中心对称的距离矩阵。如果使用eigen()函数替代cmdscale()函数来演示做MDS分析时，这就比较有用（不太理解这一段的意思），如下所示：

第三步：计算MDS图中的每个坐标占总变异的比例，如下所示：

mds.var.per <- round(mds.stuff$eig/sum(mds.stuff$eig)*100,1)
mds.var.per

结果如下所示：

> mds.var.per
 [1] 87.8  3.3  2.6  1.5  1.4  1.2  0.8  0.8  0.6  0.0

第四步：绘图，代码如下所示：

mds.values <- mds.stuff$points
mds.data <- data.frame(Sample=rownames(mds.values),
                       X=mds.values[,1],
                       Y=mds.values[,2])

ggplot(data=mds.data, aes(x=X, y=Y, label=Sample))+
  geom_text() +
  xlab(paste("MDS1 - ", mds.var.per[1], "%", sep=""))+
  ylab(paste("MDS2 - ", mds.var.per[2], "%", sep=""))+
  theme_bw()+
  ggtitle("MDS plot using Euclidean distance")

最终图形如下所示：

从MDS的这个图形中我们可以看出来，它的结果与PCA比较类似，都是wt组位于左侧，ko组位于右侧，其中MDS1所占总变异的87.8%，MDS2占总变异的3.3%。

实地上，PCA的图与MDS的图不仅看上去相似，实际上它们是一样的，我们把PCA的图与MDS的图放一起来比较看一下：

我生成的图形与教程中不太一样，教程中的这两张图基本上是一样的，并且发现，貌似MDS中的点的位置正好与PCA中的相应点的位置正好相反，此时，还不清楚怎么解决，先放这里。

使用其它的距离

前面使用的欧氏距离，此时我们使用一下其它的距离，例如表达差异的log值转换，其实对于基因表达差异的常用分析，通常是用edgeR来完成的，此处我们使用plotMDS()函数进行表达差异的log值计算，如下所示：

第一步：计算每个基因表达是不是的log2值，如下所示：

log2.data.matrix <- log2(data.matrix)

由于差异倍数的log转换值的平均值并不是一种距离度量，无法直接被dist()函数调用，因此我们需要手动来计算距离，下面的代码是创建一个空矩阵，如下所示：

log2.data.matrix <- log2(data.matrix)

log2.distance.matrix<- matrix(0,
                              nrow=ncol(log2.data.matrix),
                              ncol=ncol(log2.data.matrix),
                              dimnames = list(colnames(log2.data.matrix),
                                              colnames(log2.data.matrix)))
log2.distance.matrix

log2.distance.matrix的这个矩阵如下所示：

> log2.distance.matrix
    wt1 wt2 wt3 wt4 wt5 ko1 ko2 ko3 ko4 ko5
wt1   0   0   0   0   0   0   0   0   0   0
wt2   0   0   0   0   0   0   0   0   0   0
wt3   0   0   0   0   0   0   0   0   0   0
wt4   0   0   0   0   0   0   0   0   0   0
wt5   0   0   0   0   0   0   0   0   0   0
ko1   0   0   0   0   0   0   0   0   0   0
ko2   0   0   0   0   0   0   0   0   0   0
ko3   0   0   0   0   0   0   0   0   0   0
ko4   0   0   0   0   0   0   0   0   0   0
ko5   0   0   0   0   0   0   0   0   0   0

然后，向这个矩阵填表达差异倍数log转换值绝对值平均值，如下所示：

for(i in 1:ncol(log2.distance.matrix)){
  for(j in 1:i){
    log2.distance.matrix[i,j]<-
      mean(abs(log2.data.matrix[,i] - log2.data.matrix[,j]))
  }
}

填充后的矩阵如下所示：

> log2.distance.matrix
           wt1       wt2        wt3        wt4      wt5        ko1        ko2
wt1 0.00000000 0.0000000 0.00000000 0.00000000 0.000000 0.00000000 0.00000000
wt2 0.11484382 0.0000000 0.00000000 0.00000000 0.000000 0.00000000 0.00000000
wt3 0.09904503 0.1019653 0.00000000 0.00000000 0.000000 0.00000000 0.00000000
wt4 0.09221168 0.1037076 0.09493001 0.00000000 0.000000 0.00000000 0.00000000
wt5 0.09885778 0.1108717 0.09948529 0.09023129 0.000000 0.00000000 0.00000000
ko1 1.31334287 1.2943370 1.29838285 1.29338723 1.282153 0.00000000 0.00000000
ko2 1.29918744 1.2802797 1.28559603 1.28138419 1.268754 0.09949390 0.00000000
ko3 1.29038890 1.2721498 1.27904237 1.27301615 1.263254 0.09166573 0.09141794
ko4 1.30741667 1.2850593 1.29054531 1.28755635 1.273766 0.09707109 0.09566890
ko5 1.29516621 1.2726247 1.27697549 1.27384192 1.260347 0.10435400 0.10674809
           ko3        ko4 ko5
wt1 0.00000000 0.00000000   0
wt2 0.00000000 0.00000000   0
wt3 0.00000000 0.00000000   0
wt4 0.00000000 0.00000000   0
wt5 0.00000000 0.00000000   0
ko1 0.00000000 0.00000000   0
ko2 0.00000000 0.00000000   0
ko3 0.00000000 0.00000000   0
ko4 0.08960051 0.00000000   0
ko5 0.10106993 0.09935345   0

此时，使用新的矩阵进行MDS分析，如下所示：

mds.stuff <- cmdscale(as.dist(log2.distance.matrix),
                      eig=TRUE,
                      x.ret=TRUE)

接着计算MDS图中每个坐标所占总变异的比例，如下所示：

mds.var.per <- round(mds.stuff$eig/sum(mds.stuff$eig)*100,1)
mds.var.per

结果如下所示：

> mds.var.per
 [1] 99.1  0.3  0.2  0.1  0.1  0.1  0.1  0.1  0.0  0.0

绘图，绘图前先把数据转换成ggplot能识别的格式，代码如下所示：

mds.values <- mds.stuff$points
mds.data <- data.frame(Sample=rownames(mds.values),
                       X=mds.values[,1],
                       Y=mds.values[,2])
mds.data

结果如下：

> mds.data
    Sample          X           Y
wt1    wt1 -0.6577291  0.05852077
wt2    wt2 -0.6369172 -0.01660477
wt3    wt3 -0.6427329 -0.02665479
wt4    wt4 -0.6386594  0.01136691
wt5    wt5 -0.6262231 -0.02766676
ko1    ko1  0.6529016  0.01159350
ko2    ko2  0.6394952  0.02054915
ko3    ko3  0.6323059  0.04457294
ko4    ko4  0.6455593 -0.02576489
ko5    ko5  0.6319998 -0.04991206

此时，使用ggplot绘图，代码如下：

ggplot(data=mds.data, aes(x=X, y=Y, label=Sample)) +
  geom_text()+
  theme_bw()+
  xlab(paste("MDS1 - ", mds.var.per[1],"%", sep=""))+
  ylab(paste("MDS2 - ", mds.var.per[2],"%", sep=""))+
  ggtitle("MDS plot using avg(logFC) as the distance")

图形如下所示：

此时，我们把利用欧氏距离绘制的MDS图和表达差异log转换后的MDS图比较一下，如下所示：

从上图中可以发现，这两个图比较相似，但是不一样。我们看一下两张图的坐标轴，右侧MDS-1的是99.1%，左图是87.8%，也就是说使用了log转换后的表达倍数差异值所占总变异的比例更高。