主成份分析

主成份分析是什么

主成份分析 (Principal Component Analysis, PCA) 是一种多元分析技术¹⁾，其它常用的多元分析技术还有回归分析、聚类分析、因子分析等等。主成份分析由Karl Pearson于1901年发明²⁾，其核心思想在于，在尽可能保留数据的差异的前题下，降低数据的维度，也就是抽象出更少的互不相关的变量来描述各数据。可以通过此在线的演示来理解。
比如，一些基因在不同的样本中有个表达量，对于每个基因来说，它的表达情况就是由这些样本中的表达量来定义的。如果在坐标系中把所有的基因画出来，每个样本就是一条坐标轴，这这个表达矩阵数据集就是一群在多维空间中的点，在保持这一群点的相对空间位置不变的情况下，旋转到一个新的坐标系（坐标轴就是各个PC），使得各点在一条新的坐标轴上的坐标（投影）的方差最大，得到PC1，然后在旋转，得到次大，得到PC2，依次类推就得到了一个新的坐标系。
主成份分析可通过求取协方差矩阵的本征向量实现，本征值最大的本征向量即为PC1³⁾, 比如R的princomp函数采样的就是这种方法；也可通过奇异值分解（singular value decomposition，SVD）来实现，比如R的prcomp函数。

怎么做主成份分析

这里采用一组基因表达数据，共九个样本，第一列是基因ID，后边9列是基因表达值，共分为三个组。

dat <- read.delim("ftp://ftp.biolab.wang/data/rpkm.tsv", header=TRUE, sep='\t')  # 读入数据

看看这个数据前边几行，可以看到，第一类是geneID，后边九列是九个样本的表达数据。

nrow(dat)
###  [1] 22195
 
head(dat)
##                geneID           C1          C2           C3         Tx1        Tx2         Tx3          Ty1         Ty2         Ty3
##  1 ENSMUSG00000006169  16.08108065  13.6536716  14.69375959  15.2574267  12.868341  16.1915416  13.92268906  15.0437726  15.0177982
##  2 ENSMUSG00000025190 263.13937306 244.4119889 251.45846978 236.3022661 234.562218 231.2300947 234.18234088 253.3412910 238.5972540
##  3 ENSMUSG00000080725   0.00000000   0.0000000   0.00000000   0.0000000   0.000000   0.0000000   0.00000000   0.0000000   0.0000000
##  4 ENSMUSG00000017499   0.18095450   0.1254891   0.11109292   0.1550804   0.226933   0.1862202   0.14506731   0.2792597   0.2719463
##  5 ENSMUSG00000091383   0.00000000   0.0000000   0.00000000   1.8190357   3.802629   3.3804546   0.00000000   0.0000000   0.0000000
##  6 ENSMUSG00000031710   0.04357382   0.0000000   0.04280186   0.0000000   0.000000   0.0000000   0.05589151   0.0000000   0.0000000
 
# 把一些低表达的去掉，比如我们需要RPKM平均大于0.1的保留，这个是比较武断的，
dat <- dat[rowSums(dat[ , -1]) > 0.9, ]
 
nrow(dat)
##  [1] 15644

描述样本之间的差异

我们这里使用R的prcomp函数，这个函数把每行当做一个因变量，而每列代表一个自变量，也就是每列是一个维度，有N列，那么就是N维的数据，每个因变量的特征由这么N个维度组成。这里我们要看的是样本的情况，也就是说吧样本当做因变量，将基因看成自变量，总的有15644个基因，那么每个样本有15644维，我们很难看出来哪些样本更像，因此希望通过组成分分析，在新的坐标系下来找到一个角度更好的看出样本之间的差异。要把样本作为因变量，那么就要每个样本一行，所以我们在进行PCA之前先要将dat表达矩阵进行转置。

Dat <- t(dat[ , -1])
pca <- prcomp(Dat)  
names(pca)
##  [1] "sdev"     "rotation" "center"   "scale"    "x"

得到的结果为list，含有5个元素：pca$sdev 为各主成份的标准差，其平方为方差,即本征值，可用于碎石图及计算解释方差的百分比；pca$rotation是本征向量(eigenvectors)矩阵；pca$center 为对数据进行centering所减去的值，也就是列的均值； pca$scale是进行scaling所除以的值，也就是标准差，当scale=F的时候则无这项；pca$x PCs的值，也就是原数据在个PC上的投影。平时的作图用到的是x和sdev。

可以通过一些有用的函数来查看

summary(pca)  
##  Importance of components:
##                               PC1        PC2       PC3       PC4       PC5       PC6       PC7     PC8               PC9
##  Standard deviation     4865.6843 1557.80585 864.13727 692.80831 529.65561 441.86436 273.87476 236.647 0.000000000002937
##  Proportion of Variance    0.8475    0.08687   0.02673   0.01718   0.01004   0.00699   0.00269   0.002 0.000000000000000
##  Cumulative Proportion     0.8475    0.93437   0.96110   0.97828   0.98832   0.99531   0.99800   1.000 1.000000000000000

plot(pca)     # 等同于 screeplot(pca)  #碎石图，自己看看得到的图

PCs的值可以看到如下，这个矩阵的列名是PCs，而行名是样本名，每行一个样本。二维图一般用PC1和PC2来做。

pca$x
##            PC1        PC2        PC3        PC4        PC5        PC6        PC7        PC8                    PC9
##  C1  -8217.820  -359.0164  1139.3907 -893.29261   45.03219  325.90944  -75.65703   99.53745  0.0000000000059820948
##  C2   2162.592  1120.4435   346.8057  -29.17881 1203.39744 -371.07497  164.04526   20.12227  0.0000000000014653671
##  C3   1151.664 -1697.3988  1189.6159 1215.84519 -125.97493 -162.65799 -191.36758 -109.08803 -0.0000000000096565178
##  Tx1 -1938.648  -996.6127 -1493.2726  307.35304  339.40794  167.41032 -351.50173  254.73485  0.0000000000015678661
##  Tx2 -1813.871 -1080.2482  -540.2138 -514.89785 -489.98033 -919.74632  162.71762  -19.27196  0.0000000000045162513
##  Tx3 -5793.606  2228.2790  -609.4241  651.34727 -205.47040  170.12816  115.51534 -316.55937 -0.0000000000003566324
##  Ty1  6567.873  1772.4781   162.7143 -663.68684 -347.70749  -20.54642 -427.91542  -96.68661 -0.0000000000021449558
##  Ty2  4665.662 -1976.1140  -415.2146 -446.46695   84.04496  578.11576  271.99849 -260.63401  0.0000000000004959860
##  Ty3  3216.155   988.1894   219.5985  372.97756 -502.74937  232.46202  332.16505  427.84540  0.0000000000024031955

可以通过ggplot2作比较漂亮的图

# 先载入需要用的包，ggplot2用来作图
library(ggplot2)
# 然后把PC做坐标信息转为一个data.frame
PCs <- data.frame(pca$x)
# 添加每个样本的组名，这样方便后边作图的时候根据组来映射到颜色
PCs$Group <- c("C", "C", "C", "Tx", "Tx", "Tx", "Ty", "Ty", "Ty")
# 作图：
ggplot() +
    geom_point(data=PCs, aes(x=PC1, y=PC2, colour = Group)) +
    xlab(sprintf("PC1: %.2f%%", (100*(pca$sdev)^2)[1]/sum((pca$sdev)^2))) +
    ylab(sprintf("PC2: %.2f%%", (100*(pca$sdev)^2)[2]/sum((pca$sdev)^2)))

得到下图：

看起来C和Tx两组间的差异就不是那么明显，不能很好地区分。这可能是基因表达差异巨大，那么可能一些表达很高的基因导致了我们观察的偏差。试试，对表达值进行log转化，看看效果

Dat <- t(dat[ , -1])
Dat <- log10(Dat + 0.01)
pca <- prcomp(Dat)
 
PCs <- data.frame(pca$x)
 
PCs$Group <- c("C", "C", "C", "Tx", "Tx", "Tx", "Ty", "Ty", "Ty")
 
ggplot() +
    geom_point(data=PCs, aes(x=PC1, y=PC2, colour = Group)) +
    xlab(sprintf("PC1: %.2f%%", (100*(pca$sdev)^2)[1]/sum((pca$sdev)^2))) +
    ylab(sprintf("PC2: %.2f%%", (100*(pca$sdev)^2)[2]/sum((pca$sdev)^2)))

很显然，进过log转化之后，效果不一样了。

另外，可以对每个基因在每个样本的表达进行一个标准化，这样看的就是基因的表达趋势在样本之间的差异：

Dat <- t(dat[ , -1])
Dat <- scale(Dat)
pca <- prcomp(Dat)
 
PCs <- data.frame(pca$x)
 
PCs$Group <- c("C", "C", "C", "Tx", "Tx", "Tx", "Ty", "Ty", "Ty")
 
ggplot() +
    geom_point(data=PCs, aes(x=PC1, y=PC2, colour = Group)) +
    xlab(sprintf("PC1: %.2f%%", (100*(pca$sdev)^2)[1]/sum((pca$sdev)^2))) +
    ylab(sprintf("PC2: %.2f%%", (100*(pca$sdev)^2)[2]/sum((pca$sdev)^2)))

描述基因之间的差异

描述基因之间的差异的时候，把基因看做因变量，样本看做自变量，原来的表达矩阵就是每行一个基因，做主成份分析的时候不需要再做转置了。

# 加上行名，这样最后的结果是会带上基因名称的
row.names(dat) <- dat$geneID
 
pca <- prcomp(dat[ , -1])
 
PCs <- data.frame(pca$x)
 
ggplot() +
    geom_point(data=PCs, aes(x=PC1, y=PC2)) +
    xlab(sprintf("PC1: %.2f%%", (100*(pca$sdev)^2)[1]/sum((pca$sdev)^2))) +
    ylab(sprintf("PC2: %.2f%%", (100*(pca$sdev)^2)[2]/sum((pca$sdev)^2)))

对每个基因在每个样本的表达进行一个标准化，然后看看效果

Dat <- t(scale(t(dat[ , -1])))
 
pca <- prcomp(Dat)
 
PCs <- data.frame(pca$x)
 
ggplot() +
    geom_point(data=PCs, aes(x=PC1, y=PC2)) +
    xlab(sprintf("PC1: %.2f%%", (100*(pca$sdev)^2)[1]/sum((pca$sdev)^2))) +
    ylab(sprintf("PC2: %.2f%%", (100*(pca$sdev)^2)[2]/sum((pca$sdev)^2)))

基因之间的差异（示例数据）

我们这里做一个示例的数据集，更好地看出展示基因之间的差异效果，这里虚拟一个有20个样本，250个基因的表达矩阵。

Dat <- matrix(data = c(rnorm(10^3, mean = 1, sd = 1),    # 一千个数字，对应50行 
                       rnorm(10^3, mean = 3, sd = 1),   
                       rnorm(10^3, mean = 5, sd = 1), 
                       rnorm(10^3, mean = 6, sd = 2), 
                       rnorm(10^3, mean = 9, sd = 4)),   # 方差越大，对应后边点越散
                nrow = 250, 
                ncol = 20, 
                byrow = TRUE)                             
# 每50行对应一簇，名称中给带同样的字母，方便后边提取簇的名称
rownames(Dat) <- sprintf("%s%d", rep(letters[1:5], each=50), rep(1:50, times=5))
# 样本名称
colnames(Dat) <- LETTERS[1:20]

pca <- prcomp(Dat)  
PCs <- data.frame(pca$x)
PCs$cluster <- substr(rownames(pca$x),1,1)
ggplot() +
    geom_point(data=PCs, aes(x=PC1, y=PC2, colour=cluster)) +
    xlab(sprintf("PC1: %.2f%%", (100*(pca$sdev)^2)[1]/sum((pca$sdev)^2))) +
    ylab(sprintf("PC2: %.2f%%", (100*(pca$sdev)^2)[2]/sum((pca$sdev)^2))) +
    theme_bw() +
    theme(legend.position="none")

得到下图：

我们能够很清楚看到，表达量相近的基因会聚到一起。而方差越大，点之间越散，相互之间的相似性就越低。这里我们差不多只需要PC1就能够把这些基因分成五个簇。
很显然，如果按照20个样本作为坐标系来描述这些点，那么我们很难通过任意两个坐标轴来把这些基因之间的差异展示出来。比如我们看看A和B两个样本展示情况：

DAT <- as.data.frame(Dat) 
DAT$cluster <- substr(rownames(DAT),1,1)
ggplot() +
    geom_point(data=DAT, aes(x=A, y=B, colour=cluster)) +
    theme_bw() +
    theme(legend.position="none")

显然不容易把这些基因区分成五个簇。

本网页已被访问9154次。

¹⁾

http://www.intechopen.com/books/principal-component-analysis/

²⁾

http://stat.smmu.edu.cn/history/pearson1901.pdf

³⁾

https://www.semanticscholar.org/paper/A-tutorial-on-Principal-Components-Analysis-Smith

侧边栏

生物阁

本站推荐

最多访问

本页目录

目录

主成份分析