MDS: multidimensional scaling
PCoA: principal coordinate analysis
PCA: principal components analysis
三者最大的区别是PCA计算样本之间的correlation（比如线性关系等等），而MDS和PCoA计算样本之间的距离。

建立模拟数据集

这里用的数据集和PCA中模拟练习的数据集一样，包含10个样本和100个基因，前五个样本是野生型“WT”,后五个样本是敲除“KO”
由于这是我们模拟出来用来练手的数据，暂且就将一百个基因命名为gene1, gene2, ..., gene100
然后生成随机数（这些随机数五个一组符合泊松分布）建立模拟数据集

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)
}
head(data.matrix)
dim(data.matrix)

生成的模拟数据集

MDS or PCoA

euclidean距离

用dist建立距离矩阵

distance.matrix <- dist(scale(t(data.matrix), center=TRUE, scale=TRUE),
  method="euclidean")

就像PCA一样，样本名是行名，基因名是列名，距离矩阵一共有六种方法，在这里我们选择 method="euclidean"。
其实 method="euclidean"的距离计算方式会导致结果跟PCA一模一样，因为这种计算取得是距离的平方相加然后开根号，距离之间的详细介绍可以看这篇博文：https://blog.csdn.net/xxzhangx/article/details/53153821。

用cmdscale()计算MDS

cmdscale, Classical Muti-Dimensional scaling

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

eig=TRUE,eigen value可以帮助我们计算各个变量在距离矩阵每个aixs所占比例。
P.S.也可以用eigen()做MDS分析

计算每个aixs所占比重

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

作图

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

library(ggplot2)
 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")

和PCA一模一样

换一种距离继续计算MDS

遗传学家更常用的the average of the absolute value of the log fold change

计算每个基因log2的值

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

建立基于log2的距离矩阵

先建立一个空的矩阵

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)))

空的矩阵

向矩阵中填写log fold changes 绝对值的平均数

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

因为是对称的，只计算了一半

cmdscale

as.dist()变成机器识别的距离

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.values <- mds.stuff$points
mds.data <- data.frame(Sample=rownames(mds.values),
  X=mds.values[,1],
  Y=mds.values[,2])
mds.data

mds.data

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 fold change)

比较一下两张图，虽然相似，但是不同在于每个axis所占比重

euclidean VS the log fold