StatQuest生物统计学专题 - MDS

经典MDS（PCoA）的内在思想经典MDS在R中的操作

MDS（多维标度法，Multidimensional scaling）是一种经典的降维方法之一，它可以分为两种：度量MDS和非度量MDS。度量MDS是经典的MDS方法（classical MDS），又称为主坐标分析法（principal coordinate analysis，PCoA）。

本文将只讨论经典MDS，也就是PCoA。PCoA是根据变量间的距离关系进行计算的。

PCoA和PCA（主成分分析）很相似，使用了欧氏距离的PCoA和PCA的结果是一样的。

欧氏距离的计算

两坐标：（x1,y1)与(x2,y2)之间的欧氏距离为sqrt[x1-x2)^2 + (y1-y2)^2]

三坐标：（x1,y1,z1)与(x2,y2,z2)之间的欧氏距离为sqrt[x1-x2)^2 + (y1-y2)^2 + (z1-z2)^2]

其他类似

经典MDS（PCoA）的内在思想

同样从多个细胞的RNA-seq数据着手，每个细胞有9个基因表达数据。需要对基因进行降维，比如将9个基因降维到只需要2个纬度来解释数据。

PCoA会首先求出Cell1与Cell2之间距离，然后利用距离矩阵进行降维计算。

假如是使用的欧氏距离，那么距离就是

同理计算出cell1与cell3，cell1与cell4，cell2与cell3等等全部的两两之间的距离矩阵。

StatQuest-XVI-1

而PCoA的核心思想就是在降维的过程中保持样本点之间的距离不变。

对于上述矩阵，cell1、cell2、cell3等等是不同的数据点，他们目前是被展示在9维空间中（9个基因坐标代表一个点），而PCoA会将其降维到2维或3维（2或3个坐标表示一个点），也就是说cell1、cell2、cell3等等使用2坐标或3坐标来表示。

在此过程中，保持样本点之间的距离不变。比如cell1与cell2之间的欧氏距离距离为5.00，经过PCoA后，cell1与cell2都是2维坐标，则这两个点之间的距离依然保持在5.00（前后两个距离差异尽可能小）。

上述思想数学表达如下：

假定有以下矩阵X，行是不同的细胞，共n个细胞，列是细胞的不同基因，每个细胞d个基因。也就是说共有n个点，每个点由d维数据表示。

StatQuest-XVI-2

那么n个点之间的距离矩阵表示为：

StatQuest-XVI-3

欧式距离的计算公式为：

StatQuest-XVI-4

原始数据矩阵为X，每个细胞由d维数据表示。在经过PCoA降维之后，数据矩阵为Y，每个细胞由p维数据展示。

StatQuest-XVI-5

p<<d，一般而言p为2或3，也就是说将d维的矩阵X降维到2维或3维的矩阵Y。

由于PCoA前后，数据点之间的距离尽可能保持不变，所以PCoA就变成如下求σ(X)最优解的问题：

StatQuest-XVI-6

σ(X)是PCoA降维前后的距离差异，此值应该尽可能小。

具体的求最优解的过程此处不再赘述，详情见参考资料2的“Basic majorization theory”和“柯西-施瓦茨不等式”。

最终结果就是，将n行d列的矩阵X降维为n行2列（或3列）的矩阵Y。

经典MDS在R中的操作

由于计算欧氏距离的PCoA的结果和PCA的结果是一致，本次模拟会首先给出PCA结果与欧式距离下的PCoA的结果对比，然后在计算生物信息学中常见的log fold距离矩阵的结果。

1. 先构造原始数据：

data.matrix <- matrix(nrow=100, ncol=10)
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)
}
data.matrix <- t(data.matrix) #PCA与PCoA都是对列变量进行降维
dim(data.matrix) #10行，100列；行为细胞株，列为基因；

1. 对数据进行PCA分析

pca <- prcomp(data.matrix, scale=TRUE, center=TRUE) 

##pca$sdev代表各个主成分的解释数据波动（标准差）
pca.var <- pca$sdev^2
pca.var.per <- round(pca.var/sum(pca.var)*100, 1)

##pca$x是主成分矩阵，列为各个主成分
#构建数据框用于ggplot2绘图
pca.data <- data.frame(Sample=rownames(pca$x),
  X=pca$x[,1],
  Y=pca$x[,2])

library（ggplot2）
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("PCA Graph")+
  theme(plot.title = element_text(hjust = 0.5))

此时的PCA图如下图所示：

StatQuest-XVI-7

1. 对数据进行欧氏距离的PCoA

## 求解距离矩阵，计算距离前先使用scale函数标准化
distance.matrix <- dist(scale(data.matrix, center=TRUE, scale=TRUE),
                        method="euclidean")
## cmdscale进行MDS（PCoA），eig代表返回特征矩阵和特征向量，x.ret带包返回双中心对称距离矩阵
mds.stuff <- cmdscale(distance.matrix, eig=TRUE, x.ret=TRUE)

## 计算各个主坐标成分所解释的数据波动百分比，
#eig为特征值，其值的平方等于各主坐标成分解释的数据波动
mds.var.per <- round(mds.stuff$eig/sum(mds.stuff$eig)*100, 1)

## 构建数据框用于ggplot2绘图
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() +
  theme_bw() +
  xlab(paste("MDS1 - ", mds.var.per[1], "%", sep="")) +
  ylab(paste("MDS2 - ", mds.var.per[2], "%", sep="")) +
  ggtitle("MDS plot using Euclidean distance")+
  theme(plot.title = element_text(hjust = 0.5))

此时的成分图如下，可以看出不论是坐标轴解释的数据波动比例，还是图形都是和PCA完全一致的。也就是说使用欧氏距离的MDA（或PCoA）和PCA的结果是一样的。

StatQuest-XVI-8

1. log(fold change)距离矩阵下的MDS降维

dist函数可以计算包括欧氏距离在内共6种距离，不过这6种距离并不包括log fold change距离矩阵，需要手动构建：log(fold change) = mean(abs(log2(FC)))=mean(abs(log2 X1-log2 X2))。其中X1与X2代表不同的细胞株，其他类同。

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

## 创建空的距离矩阵，10行10列的矩阵，元素为0
log2.distance.matrix <- matrix(0,
                               nrow=nrow(log2.data.matrix),
                               ncol=nrow(log2.data.matrix),
                               dimnames=list(rownames(log2.data.matrix),
                                             rownames(log2.data.matrix)))

## 计算log fold矩阵，不同细胞株之间的logFC距离为mean(abs(log2(FC)))
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,]))
  }
}
## 将距离矩阵转换为R中的标准距离矩阵格式（不显示对角线的下三角矩阵）
log2.distance.matrix <- as.dist(log2.distance.matrix)

## mds求解
mds.stuff.logfc <- cmdscale(log2.distance.matrix,
                      eig=TRUE,
                      x.ret=TRUE)

## 计算各个坐标成分所解释的数据波动
mds.var.per <- round(mds.stuff.logfc$eig/sum(mds.stuff.logfc$eig)*100, 1)

## 构建数据框用于ggplot2绘图
mds.values <- mds.stuff.logfc$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() +
  theme_bw() +
  xlab(paste("MDS1 - ", mds.var.per[1], "%", sep="")) +
  ylab(paste("MDS2 - ", mds.var.per[2], "%", sep="")) +
  ggtitle("MDS plot using logFC as the distance")+
  theme(plot.title = element_text(hjust = 0.5))

此时的成分图如下图所示，使用logFC距离矩阵，不论是坐标轴解释的数据波动比例，还是图形都和上述的PCA和使用欧氏距离的PCoA不再一样。此时MDS1几乎可以解释全部数据变异（99.3%）。

StatQuest-XVI-9

专题以往文章

参考资料

StatQuest课程：https://statquest.org/video-index/

MDS（multidimensional scaling）多维尺度分析：https://blog.csdn.net/u010705209/article/details/53518772

---

猜你喜欢

还有更多文章，请移步公众号阅读

▼ 如果你生信基本技能已经入门，需要提高自己，请关注下面的生信技能树，看我们是如何完善生信技能，成为一个生信全栈工程师。

▼ 如果你是初学者，请关注下面的生信菜鸟团，了解生信基础名词，概念，扎实的打好基础，争取早日入门。