StatQuest生物統(tǒng)計學專題 - MDS

微笑如酒 2018-08-01

展開全文

經(jīng)典MDS（PCoA）的內(nèi)在思想經(jīng)典MDS在R中的操作

MDS（多維標度法，Multidimensional scaling）是一種經(jīng)典的降維方法之一，它可以分為兩種：度量MDS和非度量MDS。度量MDS是經(jīng)典的MDS方法（classical MDS），又稱為主坐標分析法（principal coordinate analysis，PCoA）。

本文將只討論經(jīng)典MDS，也就是PCoA。PCoA是根據(jù)變量間的距離關(guān)系進行計算的。

PCoA和PCA（主成分分析）很相似，使用了歐氏距離的PCoA和PCA的結(jié)果是一樣的。

歐氏距離的計算
兩坐標：（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]
其他類似

經(jīng)典MDS（PCoA）的內(nèi)在思想

同樣從多個細胞的RNA-seq數(shù)據(jù)著手，每個細胞有9個基因表達數(shù)據(jù)。需要對基因進行降維，比如將9個基因降維到只需要2個緯度來解釋數(shù)據(jù)。

PCoA會首先求出Cell1與Cell2之間距離，然后利用距離矩陣進行降維計算。

假如是使用的歐氏距離，那么距離就是

同理計算出cell1與cell3，cell1與cell4，cell2與cell3等等全部的兩兩之間的距離矩陣。

而PCoA的核心思想就是在降維的過程中保持樣本點之間的距離不變。

對于上述矩陣，cell1、cell2、cell3等等是不同的數(shù)據(jù)點，他們目前是被展示在9維空間中（9個基因坐標代表一個點），而PCoA會將其降維到2維或3維（2或3個坐標表示一個點），也就是說cell1、cell2、cell3等等使用2坐標或3坐標來表示。

在此過程中，保持樣本點之間的距離不變。比如cell1與cell2之間的歐氏距離距離為5.00，經(jīng)過PCoA后，cell1與cell2都是2維坐標，則這兩個點之間的距離依然保持在5.00（前后兩個距離差異盡可能小）。

上述思想數(shù)學表達如下：

假定有以下矩陣X，行是不同的細胞，共n個細胞，列是細胞的不同基因，每個細胞d個基因。也就是說共有n個點，每個點由d維數(shù)據(jù)表示。

那么n個點之間的距離矩陣表示為：

歐式距離的計算公式為：

原始數(shù)據(jù)矩陣為X，每個細胞由d維數(shù)據(jù)表示。在經(jīng)過PCoA降維之后，數(shù)據(jù)矩陣為Y，每個細胞由p維數(shù)據(jù)展示。

p<>

由于PCoA前后，數(shù)據(jù)點之間的距離盡可能保持不變，所以PCoA就變成如下求σ(X)最優(yōu)解的問題：

σ(X)是PCoA降維前后的距離差異，此值應該盡可能小。

具體的求最優(yōu)解的過程此處不再贅述，詳情見參考資料2的“Basic majorization theory”和“柯西-施瓦茨不等式”。

最終結(jié)果就是，將n行d列的矩陣X降維為n行2列（或3列）的矩陣Y。

經(jīng)典MDS在R中的操作

由于計算歐氏距離的PCoA的結(jié)果和PCA的結(jié)果是一致，本次模擬會首先給出PCA結(jié)果與歐式距離下的PCoA的結(jié)果對比，然后在計算生物信息學中常見的log fold距離矩陣的結(jié)果。

先構(gòu)造原始數(shù)據(jù)：

data.matrix <>100, ncol=10)
colnames(data.matrix) <>
  paste('wt', 1:5, sep=''),
  paste('ko', 1:5, sep=''))
rownames(data.matrix) <>'gene', 1:100, sep='')
for (i in 1:100) {
  wt.values <>5, lambda=sample(x=10:1000, size=1))
  ko.values <>5, lambda=sample(x=10:1000, size=1))

  data.matrix[i,] <>
}
data.matrix <>#PCA與PCoA都是對列變量進行降維
dim(data.matrix) #10行，100列；行為細胞株，列為基因；

對數(shù)據(jù)進行PCA分析

pca <>TRUE, center=TRUE) 

##pca$sdev代表各個主成分的解釋數(shù)據(jù)波動（標準差）
pca.var <>2
pca.var.per <- round(pca.var>100, 1)

##pca$x是主成分矩陣，列為各個主成分
#構(gòu)建數(shù)據(jù)框用于ggplot2繪圖
pca.data <>
  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圖如下圖所示：

對數(shù)據(jù)進行歐氏距離的PCoA

## 求解距離矩陣，計算距離前先使用scale函數(shù)標準化
distance.matrix <>TRUE, scale=TRUE),
                        method='euclidean')
## cmdscale進行MDS（PCoA），eig代表返回特征矩陣和特征向量，x.ret帶包返回雙中心對稱距離矩陣
mds.stuff <>TRUE, x.ret=TRUE)

## 計算各個主坐標成分所解釋的數(shù)據(jù)波動百分比，
#eig為特征值，其值的平方等于各主坐標成分解釋的數(shù)據(jù)波動
mds.var.per <- round(mds.stuff$eig>100, 1)

## 構(gòu)建數(shù)據(jù)框用于ggplot2繪圖
mds.values <>
mds.data <>
                       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))

此時的成分圖如下，可以看出不論是坐標軸解釋的數(shù)據(jù)波動比例，還是圖形都是和PCA完全一致的。也就是說使用歐氏距離的MDA（或PCoA）和PCA的結(jié)果是一樣的。

log(fold change)距離矩陣下的MDS降維

dist函數(shù)可以計算包括歐氏距離在內(nèi)共6種距離，不過這6種距離并不包括log fold change距離矩陣，需要手動構(gòu)建：log(fold change) = mean(abs(log2(FC)))=mean(abs(log2 X1-log2 X2))。其中X1與X2代表不同的細胞株，其他類同。

log2.data.matrix <>

## 創(chuàng)建空的距離矩陣，10行10列的矩陣，元素為0
log2.distance.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,]))
  }
}
## 將距離矩陣轉(zhuǎn)換為R中的標準距離矩陣格式（不顯示對角線的下三角矩陣）
log2.distance.matrix <>

## mds求解
mds.stuff.logfc <>
                      eig=TRUE,
                      x.ret=TRUE)

## 計算各個坐標成分所解釋的數(shù)據(jù)波動
mds.var.per <- round(mds.stuff.logfc$eig>100, 1)

## 構(gòu)建數(shù)據(jù)框用于ggplot2繪圖
mds.values <>
mds.data <>
                       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))