PCA usage

This is generated by R knitr, please check https://github.com/Tong-Chen/notebook/blob/master/PCA_leanr_knit.Rmd for greatly formatted documents.

library(knitr)
library(psych)
library(reshape2)
library("ggplot2")
library("ggbeeswarm")
library(scatterplot3d)
library(useful)
library(ggfortify)

主成分分析简介

主成分分析 (PCA, principal component analysis)是一种数学降维方法, 利用正交变换 (orthogonal transformation)把一系列可能线性相关的变量转换为一组线性不相关的新变量，也称为主成分，从而利用新变量在更小的维度下展示数据的特征。

主成分是原有变量的线性组合，其数目不多于原始变量。组合之后，相当于我们获得了一批新的观测数据，这些数据的含义不同于原有数据，但包含了之前数据的大部分特征，并且有着较低的维度，便于进一步的分析。

在空间上，PCA可以理解为把原始数据投射到一个新的坐标系统，第一主成分为第一坐标轴，它的含义代表了原始数据中多个变量经过某种变换得到的新变量的变化区间；第二成分为第二坐标轴，代表了原始数据中多个变量经过某种变换得到的第二个新变量的变化区间。这样我们把利用原始数据解释样品的差异转变为利用新变量解释样品的差异。

这种投射方式会有很多，为了最大限度保留对原始数据的解释，一般会用最大方差理论或最小损失理论，使得第一主成分有着最大的方差或变异数 (就是说其能尽量多的解释原始数据的差异)；随后的每一个主成分都与前面的主成分正交，且有着仅次于前一主成分的最大方差 (正交简单的理解就是两个主成分空间夹角为90°，两者之间无线性关联，从而完成去冗余操作)。

主成分分析的意义

简化运算。

在问题研究中，为了全面系统地分析问题，我们通常会收集众多的影响因素也就是众多的变量。这样会使得研究更丰富，通常也会带来较多的冗余数据和复杂的计算量。

比如我们我们测序了100种样品的基因表达谱借以通过分子表达水平的差异对这100种样品进行分类。在这个问题中，研究的变量就是不同的基因。每个基因的表达都可以在一定程度上反应样品之间的差异，但某些基因之间却有着调控、协同或拮抗的关系，表现为它们的表达值存在一些相关性，这就造成了统计数据所反映的信息存在一定程度的冗余。另外假如某些基因如持家基因在所有样本中表达都一样，它们对于解释样本的差异也没有意义。这么多的变量在后续统计分析中会增大运算量和计算复杂度，应用PCA就可以在尽量多的保持变量所包含的信息又能维持尽量少的变量数目，帮助简化运算和结果解释。
去除数据噪音。

比如说我们在样品的制备过程中，由于不完全一致的操作，导致样品的状态有细微的改变，从而造成一些持家基因也发生了相应的变化，但变化幅度远小于核心基因 (一般认为噪音的方差小于信息的方差）。而PCA在降维的过程中滤去了这些变化幅度较小的噪音变化，增大了数据的信噪比。
利用散点图实现多维数据可视化。

在上面的表达谱分析中，假如我们有1个基因，可以在线性层面对样本进行分类；如果我们有2个基因，可以在一个平面对样本进行分类；如果我们有3个基因，可以在一个立体空间对样本进行分类；如果有更多的基因，比如说n个，那么每个样品就是n维空间的一个点，则很难在图形上展示样品的分类关系。利用PCA分析，我们可以选取贡献最大的2个或3个主成分作为数据代表用以可视化。这比直接选取三个表达变化最大的基因更能反映样品之间的差异。（利用Pearson相关系数对样品进行聚类在样品数目比较少时是一个解决办法）
发现隐性相关变量。

我们在合并冗余原始变量得到主成分过程中，会发现某些原始变量对同一主成分有着相似的贡献，也就是说这些变量之间存在着某种相关性，为相关变量。同时也可以获得这些变量对主成分的贡献程度。对基因表达数据可以理解为发现了存在协同或拮抗关系的基因。

示例展示原始变量对样品的分类

假设有一套数据集，包含100个样品中某一基因的表达量。如下所示，每一行为一个样品，每一列为基因的表达值。这也是做PCA分析的基本数据组织方式，每一行代表一个样品，每一列代表一组观察数据即一个变量。

count <- 50
Gene1_a <- rnorm(count,5,0.5)
Gene1_b <- rnorm(count,20,0.5)
grp_a <- rep('a', count)
grp_b <- rep('b', count)
cy_data <- data.frame(Gene1 = c(Gene1_a, Gene1_b), Group=c(grp_a, grp_b))
cy_data <- as.data.frame(cy_data)
label <- c(paste0(grp_a, 1:count), paste0(grp_b, 1:count))
row.names(cy_data) <- label
library(knitr)
library(psych)
kable(headTail(cy_data), booktabs=TRUE, caption="Expression profile for Gene1 in 100 samples")

Expression profile for Gene1 in 100 samples
	Gene1	Group
a1	5.24	a
a2	4.79	a
a3	4.58	a
a4	5.41	a
…	…	NA
b47	20.5	b
b48	19.25	b
b49	19.64	b
b50	20.2	b

# Add additional column to data only for plotting

cy_data$Y <- rep(0,count*2)

从下图可以看出，100个样品根据Gene1表达量的不同在横轴上被被分为了2类，可以看做是在线性水平的分类。

library("ggplot2")
library("ggbeeswarm")

# geom_quasirandom:用于画Jitter Plot
# theme(axis.*.y): 去除Y轴
# xlim, ylim设定坐标轴的区间
ggplot(cy_data,aes(Gene1, Y))+geom_quasirandom(aes(color=factor(Group)))+theme(legend.position=c(0.5,0.7)) + theme(legend.title=element_blank()) + scale_fill_discrete(name="Group") + theme(axis.line.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank(), axis.title.y=element_blank()) + ylim(-0.5,5) + xlim(0,25)

那么如果有2个基因呢？

count <- 50
Gene2_a <- rnorm(count,5,0.2)
Gene2_b <- rnorm(count,5,0.2)

cy_data2 <- data.frame(Gene1 = c(Gene1_a, Gene1_b), Gene2 = c(Gene2_a, Gene2_b), Group=c(grp_a, grp_b))
cy_data2 <- as.data.frame(cy_data2)

row.names(cy_data2) <- label

kable(headTail(cy_data2), booktabs=T, caption="Expression profile for Gene1 and Gene2 in 100 samples")

Expression profile for Gene1 and Gene2 in 100 samples
	Gene1	Gene2	Group
a1	5.24	5.14	a
a2	4.79	4.91	a
a3	4.58	5.08	a
a4	5.41	5.07	a
…	…	…	NA
b47	20.5	4.92	b
b48	19.25	4.99	b
b49	19.64	4.81	b
b50	20.2	4.7	b

从下图可以看出，100个样品根据Gene1和Gene2的表达量的不同在坐标轴上被被分为了2类，可以看做是在平面水平的分类。而且在这个例子中，我们可以很容易的看出Gene1对样品分类的贡献要比Gene2大，因为Gene1在样品间的表达差异大。

ggplot(cy_data2,aes(Gene1, Gene2))+geom_point(aes(color=factor(Group)))+theme(legend.position=c(0.5,0.9)) + theme(legend.title=element_blank()) + ylim(0,10) + xlim(0,25)

如果有3个基因呢？

count <- 50
Gene3_a <- c(rnorm(count/2,5,0.2), rnorm(count/2,15,0.2))
Gene3_b <- c(rnorm(count/2,15,0.2), rnorm(count/2,5,0.2))

data3 <- data.frame(Gene1 = c(Gene1_a, Gene1_b), Gene2 = c(Gene2_a, Gene2_b), Gene3 = c(Gene3_a, Gene3_b), Group=c(grp_a, grp_b))
data3 <- as.data.frame(data3)

row.names(data3) <- label

kable(headTail(data3), booktabs=T, caption="Expression profile for 3 genes in 100 samples")

Expression profile for 3 genes in 100 samples
	Gene1	Gene2	Gene3	Group
a1	5.24	5.14	4.75	a
a2	4.79	4.91	5.13	a
a3	4.58	5.08	4.74	a
a4	5.41	5.07	5.15	a
…	…	…	…	NA
b47	20.5	4.92	4.91	b
b48	19.25	4.99	5.07	b
b49	19.64	4.81	4.63	b
b50	20.2	4.7	4.93	b

从下图可以看出，100个样品根据Gene1、Gene2和Gene3的表达量的不同在坐标轴上被被分为了4类，可以看做是立体空间的分类。而且在这个例子中，我们可以很容易的看出Gene1和Gene3对样品分类的贡献要比Gene2大。

library(scatterplot3d)
colorl <- c("#E69F00", "#56B4E9")
# Extract same number of colors as the Group and same Group would have same color.
colors <- colorl[as.numeric(data3$Group)]
scatterplot3d(data3[,1:3], color=colors, xlim=c(0,25), ylim=c(0,25), zlim=c(0,25), angle=55, pch=16)
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

当我们向由Gene1和Gene2构成的X-Y平面做垂线时，可以很明显的看出，Gene2所在的轴对样品的分类没有贡献。因为投射到X-Y屏幕上的点在Y轴几乎处于同一位置。

library(scatterplot3d)
colorl <- c("#E69F00", "#56B4E9")
colors <- colorl[as.numeric(data3$Group)]
scatterplot3d(data3[,1:3], color=colors, xlim=c(0,25), ylim=c(0,25), zlim=c(0,25),type='h')
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

我们把坐标轴做一个转换，可以看到在由Gene1和Gene3构成的X-Y平面上，样品被分为了4类。Gene2对样品的分类几乎没有贡献，因为几乎所有样品在Gene2维度上的值都一样。

library(scatterplot3d)
colorl <- c("#E69F00", "#56B4E9")
colors <- colorl[as.numeric(data3$Group)]
scatterplot3d(x=data3$Gene1, y= data3$Gene3, z= data3$Gene2, color=colors, xlim=c(0,25), ylim=c(0,25), zlim=c(0,25),type='h')
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

在上述例子中，我们可以很容易的区分出Gene1和Gene3可以作为分类的主成分，而Gene2则对分类没有帮助，可以在计算中去除。

但是如果我们测序了几万个基因的表达时，就很难通过肉眼去看，或者作出一个图供我们筛选哪些基因对样本分类贡献大。这时我们应该怎么做呢？

其中有一个方法是，在这个基因表达矩阵中选出3个变化最大的基因做为3个主成分对样品进行分类。我们试验下效果怎么样。

# 数据集来源于 http://satijalab.org/seurat/old-get-started/
# 原始下载链接 http://www.broadinstitute.org/~rahuls/seurat/seurat_files_nbt.zip
# 为了保证文章的使用，文末附有数据的新下载链接，以防原链接失效
data4 <- read.table("HiSeq301-RSEM-linear_values.txt", header=T, row.names=1,sep="\t")
dim(data4)

## [1] 23730   301

library(useful)
kable(corner(data4,r=15,c=8), booktabs=T, caption="Gene expression matrix")

Gene expression matrix
	Hi_2338_1	Hi_2338_2	Hi_2338_3	Hi_2338_4	Hi_2338_5	Hi_2338_6	Hi_2338_7	Hi_2338_8
A1BG	9.08	0.00	0.00	1.75	0.00	0.40	0.00	0.78
A1BG-AS1	0.00	0.00	3.47	0.36	0.00	0.00	0.00	0.00
A1CF	0.00	0.05	0.00	0.00	0.00	0.00	0.00	0.00
A2LD1	0.00	0.00	0.00	0.29	0.00	9.19	0.00	0.00
A2M	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
A2M-AS1	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
A2ML1	0.10	0.00	0.14	0.00	0.00	0.00	0.00	0.00
A2MP1	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
A4GALT	0.57	0.00	0.00	0.00	0.00	0.00	0.35	0.00
A4GNT	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
AA06	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
AAAS	38.95	0.00	0.00	4.44	0.00	32.90	0.00	5.58
AACS	0.12	0.00	0.00	0.00	0.58	1.03	2.16	65.74
AACSP1	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
AADAC	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00

我们筛选变异系数最大的3个基因。在这之前我们先剔除在少于5个样品中表达的基因和少于1000个表达的基因样品（这里我们把表达值不小于1的基因视为表达的基因），并把所有基因根据其在不同样品中表达值的变异系数排序。

#去除表达值全为0的行
#data4_nonzero <- data4[rowSums(data4)!=0,]

#筛选符合要求的表达的行和列
#data4_use <- data4[apply(data4,1,function(row) sum(row>=1)>=5),]
#data4_use <- data4[,apply(data4,2,function(col) sum(col>=1)>=1000),]
data4_use <- data4[rowSums(data4>=1)>5,colSums(data4>=1)>1000]

# 对于表达谱数据，因为涉及到PCR的指数扩增，一般会取log处理
# 其它数据log处理会降低数据之间的差异，不一定适用
data4_use_log2 <- log2(data4_use+1)

dim(data4_use_log2)

## [1] 16482   301

# 计算变异系数(标准差除以平均值)度量基因表达变化幅度
#cv <- apply(data4_use_log2,1,sd)/rowMeans(data4_use_log2)
# 根据变异系数排序
#data4_use_log2 <- data4_use_log2[order(cv,decreasing = T),]

# 计算中值绝对偏差 (MAD, median absolute deviation)度量基因表达变化幅度
# 在基因表达中，尽管某些基因很小的变化会导致重要的生物学意义，
# 但是很小的观察值会引入很大的背景噪音，因此也意义不大。
mads <- apply(data4_use_log2, 1, mad)
data4_use_log2 <- data4_use_log2[rev(order(mads)),]

#筛选前3列
data_var3 <- data4_use_log2[1:3,]

# 转置矩阵使得每一行为一个样品，每一列为一组变量
data_var3_forPCA <- t(data_var3)

dim(data_var3_forPCA)

## [1] 301   3

kable(corner(data_var3_forPCA, r=10,c=5), booktabs=TRUE, caption="A table of the 3 most variable genes")

A table of the 3 most variable genes
	MT2A	ANXA1	ARHGDIB
Hi_2338_1	12.211493	11.837198	9.543283
Hi_2338_2	11.306216	8.098769	8.071623
Hi_2338_3	11.926226	10.688626	9.720535
Hi_2338_4	10.974207	9.386574	9.883376
Hi_2338_5	14.603994	10.375072	9.970379
Hi_2338_6	6.904604	11.155349	10.093510
Hi_2338_7	12.436719	10.852249	7.742882
Hi_2338_8	9.798375	9.783227	6.270716
Hi_2338_9	11.743673	9.626476	9.250251
Hi_2338_10	11.240016	11.303056	0.000000

# 获得样品分组信息
sample <- rownames(data_var3_forPCA)

# 把样品名字按 <_> 分割，取出其第二部分作为样品的组名
# lapply(X, FUC) 对列表或向量中每个元素执行FUC操作，FUNC为自定义或R自带的函数
## One better way to generate group
group <- unlist(lapply(strsplit(sample, "_"), function(x) x[2]))

##One way to generate group
#sample_split <- strsplit(sample,"_")
#group <- matrix(unlist(sample_split), ncol=3, byrow=T)[,2]
print(sample[1:4])

## [1] "Hi_2338_1" "Hi_2338_2" "Hi_2338_3" "Hi_2338_4"

print(group[1:4])

## [1] "2338" "2338" "2338" "2338"

# 根据分组数目确定颜色变量
colorA <- rainbow(length(unique(group)))

# 根据每个样品的分组信息获取对应的颜色变量
colors <- colorA[as.factor(group)]

# 根据样品分组信息获得legend的颜色
colorl <- colorA[as.factor(unique(group))]

scatterplot3d(data_var3_forPCA[,1:3], color=colors)
legend("top", legend=levels(as.factor(group)), col=colorl, pch=16, xpd=T, horiz=F, ncol=6)

我们看到图中的样品并没有按照预先设定的标签完全分开。当然我们也可以通过其他方法筛选变异最大的三个基因，最终的分类效果不会相差很大。因为不管怎么筛选，我们都只用到了3个基因的表达量。

假如我们把这个数据用PCA来分类，结果是怎样的呢？

# Pay attention to the format of PCA input 
# Rows are samples and columns are variables
data4_use_log2_t <- t(data4_use_log2)

# Add group column for plotting
data4_use_log2_label <- as.data.frame(data4_use_log2_t)
data4_use_log2_label$group <- group

# By default, prcomp will centralized the data using mean.
# Normalize data for PCA by dividing each data by column standard deviation.
# Often, we would normalize data.
# Only when we care about the real number changes other than the trends,
# `scale` can be set to TRUE. 
# We will show the differences of scaling and un-scaling effects.
pca <- prcomp(data4_use_log2_t, scale=T)

# sdev: standard deviation of the principle components.
# Square to get variance
percentVar <- pca$sdev^2 / sum( pca$sdev^2)

# To check what's in pca
print(str(pca))

## List of 5
##  $ sdev    : num [1:301] 36.6 30.4 23.3 21.6 19.9 ...
##  $ rotation: num [1:16482, 1:301] -0.01133 -0.01955 -0.00199 -0.00569 -0.0204 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:16482] "MT2A" "ANXA1" "ARHGDIB" "RND3" ...
##   .. ..$ : chr [1:301] "PC1" "PC2" "PC3" "PC4" ...
##  $ center  : Named num [1:16482] 6.5 5.47 5.43 4.23 4.1 ...
##   ..- attr(*, "names")= chr [1:16482] "MT2A" "ANXA1" "ARHGDIB" "RND3" ...
##  $ scale   : Named num [1:16482] 5.62 4.96 4.78 4.13 3.98 ...
##   ..- attr(*, "names")= chr [1:16482] "MT2A" "ANXA1" "ARHGDIB" "RND3" ...
##  $ x       : num [1:301, 1:301] -37.6 -28.6 -30.6 -43.7 -11.4 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:301] "Hi_2338_1" "Hi_2338_2" "Hi_2338_3" "Hi_2338_4" ...
##   .. ..$ : chr [1:301] "PC1" "PC2" "PC3" "PC4" ...
##  - attr(*, "class")= chr "prcomp"
## NULL

从图中可以看到，数据呈现了一定的分类模式 (当然这个分类结果也不理想，我们随后再进一步优化)。

library(ggfortify)
autoplot(pca, data=data4_use_log2_label, colour="group") + xlab(paste0("PC1 (", round(percentVar[1]*100), "% variance)")) + ylab(paste0("PC2 (", round(percentVar[2]*100), "% variance)")) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + theme(legend.position="right")

采用3个主成分获得的分类效果优于2个主成分，因为这样保留的原始信息更多。

# 根据分组数目确定颜色变量
colorA <- rainbow(length(unique(group)))

# 根据每个样品的分组信息获取对应的颜色变量
colors <- colorA[as.factor(group)]

# 根据样品分组信息获得legend的颜色
colorl <- colorA[as.factor(unique(group))]

# 获得PCH symbol列表
pch_l <- as.numeric(as.factor(unique(group)))
# 产生每个样品的pch symbol
pch <- pch_l[as.factor(group)]

pc <- as.data.frame(pca$x)
scatterplot3d(x=pc$PC1, y=pc$PC2, z=pc$PC3, pch=pch, color=colors, xlab=paste0("PC1 (", round(percentVar[1]*100), "% variance)"), ylab=paste0("PC2 (", round(percentVar[2]*100), "% variance)"), zlab=paste0("PC3 (", round(percentVar[3]*100), "% variance)"))

legend(-3,8, legend=levels(as.factor(group)), col=colorl, pch=pch_l, xpd=T, horiz=F, ncol=6)

PCA的实现原理

在上面的例子中，PCA分析不是简单地选取2个或3个变化最大的基因，而是先把原始的变量做一个评估，计算各个变量各自的变异度(方差)和两两变量的相关度（协方差），得到一个协方差矩阵。在这个协方差矩阵中，对角线的值为每一个变量的方差，其它值为每两个变量的协方差。随后对原变量的协方差矩阵对角化处理，即求解其特征值和特征向量。原变量与特征向量的乘积（对原始变量的线性组合）即为新变量（回顾下线性代数中的矩阵乘法）；新变量的协方差矩阵为对角协方差矩阵且对角线上的方差由大到小排列；然后从新变量中选择信息最丰富也就是方差最大的的前2个或前3个新变量也就是主成分用以可视化。下面我们一步步阐释这是怎么做的。

我们先回忆两个数学概念，方差和协方差。方差用来表示一组一维数据的离散程度。协方差表示2组一维数据的相关性。当协方差为0时，表示两组数据完全独立。当协方差为正时，表示一组数据增加时另外一组也会增加；当协方差为负时表示一组数据增加时另外一组数据会降低（与相关系数类似）。如果我们有很多组一维数据，比如很多基因的表达数据，就会得到很多协方差，这就构成了协方差矩阵。

方差和协方差的计算公式如下：

\[ 方差 Var(X) = \frac{\sum_{i=1}^{n}(X_i-\bar X)^2}{n-1} \]

\[ 协方差 cov(X,Y) = \frac{\sum_{i=1}^{n}(X_i-\bar X)(Y_i-\bar Y)}{n-1} \]

如果数据的均值为0，这个公式可以进一步简化。简化后的公式把计算协方差转变为了矩阵乘法运算。这也是为什么PCA需要中心化数据。

\[ 方差 Var(X) = \frac{\sum_{i=1}^{n}{X_i}^2}{n-1} \]

\[ 协方差 cov(X,Y) = \frac{\sum_{i=1}^{n}X_iY_i}{n-1} \]

\[ 协方差矩阵 cov(X,Y) = \frac{X_{n,m}^T Y_{n,m}}{n-1} \]

假如我们有一个矩阵如下，

mat <- as.data.frame(matrix(rnorm(20,0,1), nrow=4))
colnames(mat) <- paste0("Gene_", letters[1:5])
rownames(mat) <- paste0("Samp_", 1:4)
mat

##             Gene_a     Gene_b    Gene_c     Gene_d     Gene_e
## Samp_1 -1.66213745 -0.5827805 0.5880888  1.0473937 -0.2893133
## Samp_2  1.61285131 -0.1445312 2.4781586 -0.1870570 -0.5177238
## Samp_3 -0.26850170  1.3214449 0.4524715  1.1733962  0.6180341
## Samp_4 -0.07167303 -0.3738574 0.6320398  0.8685486 -0.3860224

平均值中心化 (mean centering)：中心化数据使其平均值为0

# mean-centering data for columns
# Get mean-value matrix first
mat_mean_norm <- mat - rep(colMeans(mat),rep.int(nrow(mat),ncol(mat)))
mat_mean_norm

##             Gene_a     Gene_b     Gene_c     Gene_d     Gene_e
## Samp_1 -1.56477223 -0.6378495 -0.4496009  0.3218233 -0.1455569
## Samp_2  1.71021653 -0.1996002  1.4404689 -0.9126274 -0.3739674
## Samp_3 -0.17113648  1.2663760 -0.5852181  0.4478258  0.7617905
## Samp_4  0.02569219 -0.4289263 -0.4056499  0.1429783 -0.2422661

# mean-centering using scale for columns
scale(mat, center=T, scale=F)

##             Gene_a     Gene_b     Gene_c     Gene_d     Gene_e
## Samp_1 -1.56477223 -0.6378495 -0.4496009  0.3218233 -0.1455569
## Samp_2  1.71021653 -0.1996002  1.4404689 -0.9126274 -0.3739674
## Samp_3 -0.17113648  1.2663760 -0.5852181  0.4478258  0.7617905
## Samp_4  0.02569219 -0.4289263 -0.4056499  0.1429783 -0.2422661
## attr(,"scaled:center")
##      Gene_a      Gene_b      Gene_c      Gene_d      Gene_e 
## -0.09736521  0.05506894  1.03768966  0.72557037 -0.14375634

中位数中心化 (median centering)：如果数据变换范围很大或有异常值，中位数标准化效果会更好。

# median-centering data for columns
mat_median_norm <- mat - rep(apply(mat,2,median),rep.int(nrow(mat),ncol(mat)))
mat_mean_norm

##             Gene_a     Gene_b     Gene_c     Gene_d     Gene_e
## Samp_1 -1.56477223 -0.6378495 -0.4496009  0.3218233 -0.1455569
## Samp_2  1.71021653 -0.1996002  1.4404689 -0.9126274 -0.3739674
## Samp_3 -0.17113648  1.2663760 -0.5852181  0.4478258  0.7617905
## Samp_4  0.02569219 -0.4289263 -0.4056499  0.1429783 -0.2422661

我们可以计算Gene_a的方差为 1.8011002 (var(mat$Gene_a))；Gene_a和Gene_b的协方差为0.1429955。

mat中5组基因的表达值的方差计算如下：

apply(mat,2,var)

##    Gene_a    Gene_b    Gene_c    Gene_d    Gene_e 
## 1.8011002 0.7447927 0.9280412 0.3858166 0.2666853

mat中5组基因表达值的协方差计算如下：

cov(mat)

##            Gene_a     Gene_b     Gene_c     Gene_d     Gene_e
## Gene_a  1.8011002  0.1429955  1.0855890 -0.7124455 -0.1827988
## Gene_b  0.1429955  0.7447927 -0.1892841  0.1608915  0.4120383
## Gene_c  1.0855890 -0.1892841  0.9280412 -0.5931261 -0.2735948
## Gene_d -0.7124455  0.1608915 -0.5931261  0.3858166  0.2003200
## Gene_e -0.1827988  0.4120383 -0.2735948  0.2003200  0.2666853

如果均值为0，数值矩阵的协方差矩阵为矩阵的乘积（实际上是矩阵的转置与其本身的乘积除以变量的维数减1）。

# Covariance matrix for Mean normalized matrix
cov(mat_mean_norm)

##            Gene_a     Gene_b     Gene_c     Gene_d     Gene_e
## Gene_a  1.8011002  0.1429955  1.0855890 -0.7124455 -0.1827988
## Gene_b  0.1429955  0.7447927 -0.1892841  0.1608915  0.4120383
## Gene_c  1.0855890 -0.1892841  0.9280412 -0.5931261 -0.2735948
## Gene_d -0.7124455  0.1608915 -0.5931261  0.3858166  0.2003200
## Gene_e -0.1827988  0.4120383 -0.2735948  0.2003200  0.2666853

# Covariance matrix for Mean normalized matrix 
# crossprod: matrix multiplication
crossprod(as.matrix(mat_mean_norm)) / (nrow(mat_mean_norm)-1)

##            Gene_a     Gene_b     Gene_c     Gene_d     Gene_e
## Gene_a  1.8011002  0.1429955  1.0855890 -0.7124455 -0.1827988
## Gene_b  0.1429955  0.7447927 -0.1892841  0.1608915  0.4120383
## Gene_c  1.0855890 -0.1892841  0.9280412 -0.5931261 -0.2735948
## Gene_d -0.7124455  0.1608915 -0.5931261  0.3858166  0.2003200
## Gene_e -0.1827988  0.4120383 -0.2735948  0.2003200  0.2666853

# Use %*% for matrix multiplication (slower)
t(as.matrix(mat_mean_norm)) %*% as.matrix(mat_mean_norm) / (nrow(mat_mean_norm)-1)

##            Gene_a     Gene_b     Gene_c     Gene_d     Gene_e
## Gene_a  1.8011002  0.1429955  1.0855890 -0.7124455 -0.1827988
## Gene_b  0.1429955  0.7447927 -0.1892841  0.1608915  0.4120383
## Gene_c  1.0855890 -0.1892841  0.9280412 -0.5931261 -0.2735948
## Gene_d -0.7124455  0.1608915 -0.5931261  0.3858166  0.2003200
## Gene_e -0.1827988  0.4120383 -0.2735948  0.2003200  0.2666853

用矩阵形式书写如下，便于理解

\[ \mathbf{cov(mat)} = \frac{1}{3} \left[\begin{array} {rrrr} -1.56 & 1.71 & -0.17 & 0.03 \\ -0.64 & -0.2 & 1.27 & -0.43 \\ -0.45 & 1.44 & -0.59 & -0.41 \\ 0.32 & -0.91 & 0.45 & 0.14 \\ -0.15 & -0.37 & 0.76 & -0.24 \\ \end{array}\right] \left[\begin{array} {rrrrr} -1.56 & -0.64 & -0.45 & 0.32 & -0.15 \\ 1.71 & -0.2 & 1.44 & -0.91 & -0.37 \\ -0.17 & 1.27 & -0.59 & 0.45 & 0.76 \\ 0.03 & -0.43 & -0.41 & 0.14 & -0.24 \\ \end{array}\right] = \left[\begin{array} {rrrrr} 1.8 & 0.14 & 1.09 & -0.71 & -0.18 \\ 0.14 & 0.74 & -0.19 & 0.16 & 0.41 \\ 1.09 & -0.19 & 0.93 & -0.59 & -0.27 \\ -0.71 & 0.16 & -0.59 & 0.39 & 0.2 \\ -0.18 & 0.41 & -0.27 & 0.2 & 0.27 \\ \end{array}\right] \]

根据前面的描述，原始变量的协方差矩阵表示原始变量自身的方差（协方差矩阵的主对角线位置）和原始变量之间的相关程度(非主对角线位置)。如果从这些数据中筛选主成分，则要选择方差大(主对角线值大)，且与其它已选变量之间相关性最小的变量（非主对角线值很小）。如果这些原始变量之间毫不相关，则它们的协方差矩阵在除主对角线处外其它地方的值都为0，这种矩阵成为对角矩阵。

而做PCA分析就是产生一组新的变量，使得新变量的协方差矩阵为对角阵，满足上面的要求。从而达到去冗余的目的。然后再选取方差大的变量，实现降维和去噪。

如果正向推导，这种组合可能会有很多种，一一计算会比较麻烦。那反过来看呢？我们不去寻找这种组合，而是计算如何使原变量的协方差矩阵变为对角阵。

数学推导中谨记的两个概念：

假设: 把未求解到的变量假设出来，用符号代替；这样有助于思考和演算
逆向：如果正向推导求不出，不妨倒着来；尽量多的利用已有信息

前面提到，新变量($Y_{m,k}$)是原始变量($X_{m,n}$)(原始变量的协方差矩阵为($C_{n,n}$))的线性组合，那么假设我们找到了这么一个线性组合(命名为特征矩阵($P_{n,k}$))，得到一组新变量$Y_{m,k}=X_{m,n}P_{n,k}$，并且新变量的协方差矩阵($D_{k,k}$)为对角阵。那么这个特征矩阵($P_{n,k}$)需要符合什么条件呢？

\[ \begin{eqnarray} D_{k,k} &=& \frac{1}{m-1} Y_{m,k}^T Y_{m,k} \\ &=& \frac{1}{m-1} (X_{m,n}P_{n,k})^T (X_{m,n}P_{n,k}) \\ &=& \frac{1}{m-1} P_{n,k}^T X_{m,n}^T X_{m,n} P_{n,k} \\ &=& P_{n,k}^T (\frac{1}{m-1} X_{m,n}^T X_{m,n}) P_{n,k} \\ &=& P_{n,k}^T \frac{1}{m-1} C_{n,n} P_{n,k} \end{eqnarray} \]

从矩阵运算可以看出，最终的特征矩阵($P_{n,k}$)需要把原变量协方差矩阵($C_{n,n}$)转换为对角阵(因为新变量的协方差矩阵($D_{k,k}$)为对角阵)，并且对角元素从大到小排列（保证每个主成分的贡献度依次降低）。

现在就把求解新变量的任务转变为了求解原变量协方差矩阵的对角化问题了。在线性代数中，矩阵对角化的问题就是求解矩阵的特征值和特征向量的问题。

我们举一个例子讲述怎么求解特征值和特征向量。

假设$A_{n,n}$为n阶对称阵，如存在$\lambda$和非零向量$x$，使得$Ax=\lambda x$，则称$\lambda$为矩阵$A_{n,n}$的特征值，非零向量$x$为为矩阵$A_{n,n}$对应于特征值$\lambda$的特征向量。

根据这个定义可以得出$(A_{n,n} - \lambda E)x = 0$，由于$x$为非零向量，所以行列式$|A-\lambda E| = 0$。

\[ \begin{eqnarray} A_{n,n} - \lambda E &=& 0 \\ \left[\begin{array} {cccc} a_{11}-\lambda & a_{12} & ... & a_{1n} \\ a_{21} & a_{22} -\lambda & ... & a_{2n} \\ ... & ... & ... & ... \\ a_{n1} & a_{n2} & ... & a_{nn}-\lambda \end{array}\right] &=& 0 \end{eqnarray} \]

由此求解出n个根$\lambda_{1}, \lambda_{2}, ..., \lambda_{3}$就是矩阵$A$的特征值。

回顾下行列式的计算：

行列式的值为行列式第一列的每一个数乘以它的余子式（余子式是行列式中除去当前元素所在行和列之后剩下的行列式）。
当行列式中存在线性相关的行或列或者有一行或一列元素全为0时，行列式的值为0。
上三角形行列式的值为其主对角线上元素的乘积。
互换行列式的两行或两列，行列式变号。
行列式的某一列（行）乘以同意书加到另一列（列）对应元素上去，行列式不变。

假如我们有一个矩阵 $\mathbf{A} = \left[\begin{array} {cc} 3 & -1 \\ -1 & 3 \end{array}\right]$，如何计算它的特征值和特征向量呢？

\begin{eqnarray} |A_{n,n} - \lambda E| &=& \left|\begin{array} {cc} 3-\lambda & -1 \\ -1 & 3-\lambda \end{array}\right| \\ &=& (3-\lambda)^2-1 \\ &=& 0 \end{eqnarray}

则$\lambda$的值为2或4。

对$\lambda_{1}=2$时，求解$(A-2E)x=\left|\begin{array} {cc} 1 & -1 \\ -1 & 1 \end{array}\right| x= 0$，得$x=k\left|\begin{array}{c} 1 \\ 1 \end{array}\right|$，则对应于$\lambda_{1}=2$时的特征向量$p_{1}=\left|\begin{array}{c} 1 \\ 1 \end{array}\right|$

对$\lambda_{2}=4$时，求解$(A-2E)x=\left|\begin{array} {cc} -1 & -1 \\ -1 & -1 \end{array}\right| x= 0$，得$x=k\left|\begin{array}{c} 1 \\ -1 \end{array}\right|$，则对应于$\lambda_{2}=4$时的特征向量$p_{2}=\left|\begin{array}{c} 1 \\ -1 \end{array}\right|$

以上就完成了PCA的数学推导。

简单的PCA实现

我们使用前面用到的数据data3来演示下如何用R函数实现PCA的计算，并与R中自带的prcomp做个比较。

library(knitr)
kable(headTail(data3), booktabs=T, caption="Expression profile for 3 genes in 100 samples")

Expression profile for 3 genes in 100 samples
	Gene1	Gene2	Gene3	Group
a1	5.24	5.14	4.75	a
a2	4.79	4.91	5.13	a
a3	4.58	5.08	4.74	a
a4	5.41	5.07	5.15	a
…	…	…	…	NA
b47	20.5	4.92	4.91	b
b48	19.25	4.99	5.07	b
b49	19.64	4.81	4.63	b
b50	20.2	4.7	4.93	b

标准化数据

data3_center_scale <- scale(data3[,1:3], center=T, scale=T)
kable(headTail(data3_center_scale), booktabs=T, caption="Normalized expression for 3 genes in 100 samples")

Normalized expression for 3 genes in 100 samples
	Gene1	Gene2	Gene3
a1	-0.95	0.61	-1.04
a2	-1.01	-0.51	-0.96
a3	-1.04	0.28	-1.04
a4	-0.93	0.27	-0.96
…	…	…	…
b47	1.07	-0.48	-1.01
b48	0.9	-0.12	-0.97
b49	0.96	-1.03	-1.06
b50	1.03	-1.57	-1

计算协方差矩阵

data3_center_scale_cov <- cov(data3_center_scale)
kable(data3_center_scale_cov, booktabs=T, caption="Covariance matrix for 3 genes in 100 samples")

Covariance matrix for 3 genes in 100 samples
	Gene1	Gene2	Gene3
Gene1	1.0000000	-0.0053055	0.0030832
Gene2	-0.0053055	1.0000000	-0.0060229
Gene3	0.0030832	-0.0060229	1.0000000

求解特征值和特征向量

data3_center_scale_cov_eigen <- eigen(data3_center_scale_cov)

# 特征值，从大到小排序
data3_center_scale_cov_eigen$values

## [1] 1.0097072 0.9969515 0.9933413

# 特征向量, 每一列为对应特征值的特征向量
data3_center_scale_cov_eigen$vectors

##            [,1]        [,2]      [,3]
## [1,]  0.5268891  0.76546199 0.3693993
## [2,] -0.6370595  0.06797271 0.7678118
## [3,]  0.5626217 -0.63988097 0.5234589

产生新的矩阵

pc_select = 3
label = paste0("PC",c(1:pc_select))
data3_new <- data3_center_scale %*% data3_center_scale_cov_eigen$vectors[,1:pc_select]
colnames(data3_new) <- label
kable(headTail(data3_new), booktabs=T, caption="PCA generated matrix for the expression of 3 genes in 100 samples")

PCA generated matrix for the expression of 3 genes in 100 samples
	PC1	PC2	PC3
a1	-1.48	-0.02	-0.42
a2	-0.75	-0.19	-1.27
a3	-1.31	-0.11	-0.71
a4	-1.2	-0.08	-0.64
…	…	…	…
b47	0.3	1.43	-0.5
b48	0	1.31	-0.27
b49	0.56	1.34	-1
b50	0.98	1.32	-1.35

比较原始数据和新产生的主成分对样品的聚类

#library(scatterplot3d)
colorl <- c("#E69F00", "#56B4E9")
# Extract same number of colors as the Group and same Group would have same color.
colors <- colorl[as.numeric(data3$Group)]

# 1 row 2 columns
par(mfrow=c(1,2))

scatterplot3d(data3[,1:3], color=colors, angle=55, pch=16, main="Original data")
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)


scatterplot3d(data3_new, color=colors,angle=55, pch=16, main="Principle components")
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

#par(mfrow=c(1,1))

利用prcomp进行主成分分析

pca_data3 <- prcomp(data3[,1:3], center=TRUE, scale=TRUE)

#Show whats in the result returned by prcomp
str(pca_data3)

## List of 5
##  $ sdev    : num [1:3] 1.005 0.998 0.997
##  $ rotation: num [1:3, 1:3] 0.527 -0.637 0.563 -0.765 -0.068 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:3] "Gene1" "Gene2" "Gene3"
##   .. ..$ : chr [1:3] "PC1" "PC2" "PC3"
##  $ center  : Named num [1:3] 12.42 5.02 9.98
##   ..- attr(*, "names")= chr [1:3] "Gene1" "Gene2" "Gene3"
##  $ scale   : Named num [1:3] 7.556 0.205 5.034
##   ..- attr(*, "names")= chr [1:3] "Gene1" "Gene2" "Gene3"
##  $ x       : num [1:100, 1:3] -1.477 -0.746 -1.311 -1.201 -2.081 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:100] "a1" "a2" "a3" "a4" ...
##   .. ..$ : chr [1:3] "PC1" "PC2" "PC3"
##  - attr(*, "class")= chr "prcomp"

# 新的数据，与前面计算的抑制
data3_pca_new <- pca_data3$x
kable(headTail(data3_pca_new), booktabs=T, caption="PCA generated matrix usig princomp for the expression of 3 genes in 100 samples")

PCA generated matrix usig princomp for the expression of 3 genes in 100 samples
	PC1	PC2	PC3
a1	-1.48	0.02	-0.42
a2	-0.75	0.19	-1.27
a3	-1.31	0.11	-0.71
a4	-1.2	0.08	-0.64
…	…	…	…
b47	0.3	-1.43	-0.5
b48	0	-1.31	-0.27
b49	0.56	-1.34	-1
b50	0.98	-1.32	-1.35

# 特征向量，与我们前面计算的一致(特征向量的符号是任意的)
pca_data3$rotation

##              PC1         PC2       PC3
## Gene1  0.5268891 -0.76546199 0.3693993
## Gene2 -0.6370595 -0.06797271 0.7678118
## Gene3  0.5626217  0.63988097 0.5234589

比较手动实现的PCA与prcomp实现的PCA的聚类结果

#library(scatterplot3d)
colorl <- c("#E69F00", "#56B4E9")
# Extract same number of colors as the Group and same Group would have same color.
colors <- colorl[as.numeric(data3$Group)]

# 1 row 2 columns
par(mfrow=c(1,2))

scatterplot3d(data3_new, color=colors,angle=55, pch=16, main="PCA by steps")
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

scatterplot3d(data3_pca_new, color=colors,angle=55, pch=16, main="PCA by prcomp")
legend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)

招聘

Tips

Lists

Cat.

Tags

About

Home

PCA usage

October 30, 2016

PCA usage

主成分分析简介

主成分分析的意义

示例展示原始变量对样品的分类

PCA的实现原理

简单的PCA实现