library(knitr)nlibrary(psych)nlibrary(reshape2)nlibrary(ggplot2)nlibrary(ggbeeswarm)nlibrary(scatterplot3d)nlibrary(useful)nlibrary(ggfortify)nmat_show <- function(matr) {n printmrow <- function(x) {n ret <- paste(paste(x,collapse = " & "),"n")n sprintf(ret)n }n align_str <- paste0('{',paste0(rep('r',ncol(matr)), collapse=""),'}')n format_mat <- apply(matr,1,printmrow)n add_env <- paste("left[begin{array}", align_str, n paste(format_mat, collapse=' '),"end{array}right]")n return(add_env)n}

主成分分析简介

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

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

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

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

主成分分析的意义

简化运算。
在问题研究中，为了全面系统地分析问题，我们通常会收集众多的影响因素也就是众多的变量。这样会使得研究更丰富，通常也会带来较多的冗余数据和复杂的计算量。
比如我们我们测序了100种样品的基因表达谱借以通过分子表达水平的差异对这100种样品进行分类。在这个问题中，研究的变量就是不同的基因。每个基因的表达都可以在一定程度上反应样品之间的差异，但某些基因之间却有着调控、协同或拮抗的关系，表现为它们的表达值存在一些相关性，这就造成了统计数据所反映的信息存在一定程度的冗余。另外假如某些基因如持家基因在所有样本中表达都一样，它们对于解释样本的差异也没有意义。这么多的变量在后续统计分析中会增大运算量和计算复杂度，应用PCA就可以在尽量多的保持变量所包含的信息又能维持尽量少的变量数目，帮助简化运算和结果解释。
去除数据噪音。
比如说我们在样品的制备过程中，由于不完全一致的操作，导致样品的状态有细微的改变，从而造成一些持家基因也发生了相应的变化，但变化幅度远小于核心基因 (一般认为噪音的方差小于信息的方差）。而PCA在降维的过程中滤去了这些变化幅度较小的噪音变化，增大了数据的信噪比。
利用散点图实现多维数据可视化。
在上面的表达谱分析中，假如我们有1个基因，可以在线性层面对样本进行分类；如果我们有2个基因，可以在一个平面对样本进行分类；如果我们有3个基因，可以在一个立体空间对样本进行分类；如果有更多的基因，比如说n个，那么每个样品就是n维空间的一个点，则很难在图形上展示样品的分类关系。利用PCA分析，我们可以选取贡献最大的2个或3个主成分作为数据代表用以可视化。这比直接选取三个表达变化最大的基因更能反映样品之间的差异。（利用Pearson相关系数对样品进行聚类在样品数目比较少时是一个解决办法）
发现隐性相关变量。
我们在合并冗余原始变量得到主成分过程中，会发现某些原始变量对同一主成分有着相似的贡献，也就是说这些变量之间存在着某种相关性，为相关变量。同时也可以获得这些变量对主成分的贡献程度。对基因表达数据可以理解为发现了存在协同或拮抗关系的基因。

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

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

count <- 50nGene1_a <- rnorm(count,5,0.5)nGene1_b <- rnorm(count,20,0.5)ngrp_a <- rep('a', count)ngrp_b <- rep('b', count)ncy_data <- data.frame(Gene1 = c(Gene1_a, Gene1_b), Group=c(grp_a, grp_b))ncy_data <- as.data.frame(cy_data)nlabel <- c(paste0(grp_a, 1:count), paste0(grp_b, 1:count))nrow.names(cy_data) <- labelnlibrary(knitr)nlibrary(psych)n# Add additional column to data only for plottingncy_data$Y <- rep(0,count*2)nkable(headTail(cy_data), booktabs=TRUE, caption="Expression profile for Gene1 in 100 samples")

Expression profile for Gene1 in 100 samples

Gene1

Group

5.99

4.66

4.92

4.79

…

b47

20.78

b48

20.5

b49

20.46

b50

19.94

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

library("ggplot2")nlibrary("ggbeeswarm")n# geom_quasirandom:用于画Jitter Plotn# theme(axis.*.y): 去除Y轴n# xlim, ylim设定坐标轴的区间nggplot(cy_data,aes(Gene1, Y))+geom_quasirandom(aes(color=factor(Group)), groupOnX=FALSE)+n theme(legend.position=c(0.5,0.7)) + theme(legend.title=element_blank()) + n scale_fill_discrete(name="Group") + theme(axis.line.y=element_blank(), n axis.text.y=element_blank(), axis.ticks.y=element_blank(), n axis.title.y=element_blank()) + ylim(-0.5,5) + xlim(0,25)n 一文读懂PCA分析（原理、算法、解释和可视化）

那么如果有2个基因呢？

count <- 50nGene2_a <- rnorm(count,5,0.2)nGene2_b <- rnorm(count,5,0.2)ncy_data2 <- data.frame(Gene1 = c(Gene1_a, Gene1_b), Gene2 = c(Gene2_a, Gene2_b), n Group=c(grp_a, grp_b))ncy_data2 <- as.data.frame(cy_data2)nrow.names(cy_data2) <- labelnkable(headTail(cy_data2), booktabs=T, n caption="Expression profile for Gene1 and Gene2 in 100 samples")

Expression profile for Gene1 and Gene2 in 100 samples

Gene1

Gene2

Group

5.99

5.35

4.66

5.31

4.92

5.12

4.79

5.11

…

b47

20.78

4.95

b48

20.5

4.9

b49

20.46

4.85

b50

19.94

4.87

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

ggplot(cy_data2,aes(Gene1, Gene2))+geom_point(aes(color=factor(Group)))+n theme(legend.position=c(0.5,0.9)) + theme(legend.title=element_blank()) + n ylim(0,10) + xlim(0,25) 一文读懂PCA分析（原理、算法、解释和可视化）

如果有3个基因呢？

count <- 50nGene3_a <- c(rnorm(count/2,5,0.2), rnorm(count/2,15,0.2))nGene3_b <- c(rnorm(count/2,15,0.2), rnorm(count/2,5,0.2))ndata3 <- data.frame(Gene1 = c(Gene1_a, Gene1_b), Gene2 = c(Gene2_a, Gene2_b), n Gene3 = c(Gene3_a, Gene3_b), Group=c(grp_a, grp_b))ndata3 <- as.data.frame(data3)ndata3$Group <- as.factor(data3$Group)nrow.names(data3) <- labelnkable(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

5.99

5.35

5.14

4.66

5.31

5.05

4.92

5.12

4.88

4.79

5.11

4.8

…

b47

20.78

4.95

4.91

b48

20.5

4.9

5.11

b49

20.46

4.85

4.95

b50

19.94

4.87

5.18

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

library(scatterplot3d)ncolorl <- c("#E69F00", "#56B4E9")n# Extract same number of colors as the Group and same Group would have same color.ncolors <- colorl[as.numeric(data3$Group)]nscatterplot3d(data3[,1:3], color=colors, xlim=c(0,25), ylim=c(0,25), zlim=c(0,25), n angle=55, pch=16)nlegend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T) 一文读懂PCA分析（原理、算法、解释和可视化）

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

library(scatterplot3d)ncolorl <- c("#E69F00", "#56B4E9")ncolors <- colorl[as.numeric(data3$Group)]nscatterplot3d(data3[,1:3], color=colors, xlim=c(0,25), ylim=c(0,25), zlim=c(0,25),type='h')nlegend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T) 一文读懂PCA分析（原理、算法、解释和可视化）

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

library(scatterplot3d)ncolorl <- c("#E69F00", "#56B4E9")ncolors <- colorl[as.numeric(data3$Group)]nscatterplot3d(x=data3$Gene1, y= data3$Gene3, z= data3$Gene2, color=colors, n xlim=c(0,25), ylim=c(0,25), zlim=c(0,25),type='h')nlegend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T) 一文读懂PCA分析（原理、算法、解释和可视化）

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

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

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

# 数据集来源于 http://satijalab.org/seurat/old-get-started/n# 原始下载链接 http://www.broadinstitute.org/~rahuls/seurat/seurat_files_nbt.zipn# 为了保证文章的使用，文末附有数据的新下载链接，以防原链接失效ndata4 <- read.table("data/HiSeq301-RSEM-linear_values.txt", header=T, row.names=1,sep="t")ndim(data4)n## [1] 23730 301nlibrary(useful)nkable(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

1.75

0.00

0.40

0.00

0.78

A1BG-AS1

0.00

3.47

0.36

0.00

A1CF

0.00

0.05

0.00

A2LD1

0.00

0.29

0.00

9.19

0.00

A2M

0.00

A2M-AS1

0.00

A2ML1

0.10

0.00

0.14

0.00

A2MP1

0.00

A4GALT

0.57

0.00

0.35

0.00

A4GNT

0.00

AA06

0.00

AAAS

38.95

0.00

4.44

0.00

32.90

0.00

5.58

AACS

0.12

0.00

0.58

1.03

2.16

65.74

AACSP1

0.00

AADAC

0.00

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

#去除表达值全为0的行n#data4_nonzero <- data4[rowSums(data4)!=0,]n#筛选符合要求的表达的行和列n#data4_use <- data4[apply(data4,1,function(row) sum(row>=1)>=5),]n#data4_use <- data4[,apply(data4,2,function(col) sum(col>=1)>=1000),]ndata4_use <- data4[rowSums(data4>=1)>5,colSums(data4>=1)>1000]n# 对于表达谱数据，因为涉及到PCR的指数扩增，一般会取log处理n# 其它数据log处理会降低数据之间的差异，不一定适用ndata4_use_log2 <- log2(data4_use+1)ndim(data4_use_log2)n## [1] 16482 301n# 计算变异系数(标准差除以平均值)度量基因表达变化幅度n#cv <- apply(data4_use_log2,1,sd)/rowMeans(data4_use_log2)n# 根据变异系数排序n#data4_use_log2 <- data4_use_log2[order(cv,decreasing = T),]n# 计算中值绝对偏差 (MAD, median absolute deviation)度量基因表达变化幅度n# 在基因表达中，尽管某些基因很小的变化会导致重要的生物学意义，n# 但是很小的观察值会引入很大的背景噪音，因此也意义不大。nmads <- apply(data4_use_log2, 1, mad)ndata4_use_log2 <- data4_use_log2[rev(order(mads)),]n#筛选前3行ndata_var3 <- data4_use_log2[1:3,]n# 转置矩阵使得每一行为一个样品，每一列为一组变量ndata_var3_forPCA <- t(data_var3)ndim(data_var3_forPCA)n## [1] 301 3nkable(corner(data_var3_forPCA, r=10,c=5), booktabs=TRUE, n 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

# 获得样品分组信息nsample <- rownames(data_var3_forPCA)n# 把样品名字按 <_> 分割，取出其第二部分作为样品的组名n# lapply(X, FUC) 对列表或向量中每个元素执行FUC操作，FUNC为自定义或R自带的函数n## One better way to generate groupngroup <- unlist(lapply(strsplit(sample, "_"), function(x) x[2]))n##One way to generate groupn#sample_split <- strsplit(sample,"_")n#group <- matrix(unlist(sample_split), ncol=3, byrow=T)[,2]nprint(sample[1:4])n## [1] "Hi_2338_1" "Hi_2338_2" "Hi_2338_3" "Hi_2338_4"nprint(group[1:4])n## [1] "2338" "2338" "2338" "2338"ndata_var3_scatter <- as.data.frame(data_var3_forPCA)ndata_var3_scatter$group <- groupnkable(corner(data_var3_scatter, r=10,c=5), booktabs=TRUE, n caption="A table of the 3 most variable genes")

A table of the 3 most variable genes

MT2A

ANXA1

ARHGDIB

group

Hi_2338_1

12.211493

11.837198

9.543283

2338

Hi_2338_2

11.306216

8.098769

8.071623

2338

Hi_2338_3

11.926226

10.688626

9.720535

2338

Hi_2338_4

10.974207

9.386574

9.883376

2338

Hi_2338_5

14.603994

10.375072

9.970379

2338

Hi_2338_6

6.904604

11.155349

10.093510

2338

Hi_2338_7

12.436719

10.852249

7.742882

2338

Hi_2338_8

9.798375

9.783227

6.270716

2338

Hi_2338_9

11.743673

9.626476

9.250251

2338

Hi_2338_10

11.240016

11.303056

0.000000

2338

library(reshape2)nlibrary(ggplot2)ndata_var3_melt <- melt(data_var3_scatter, id.vars=c("group"))nkable(corner(data_var3_melt, r=10,c=5), booktabs=TRUE, n caption="A table of the 3 most variable genes in melted format")

A table of the 3 most variable genes in melted format

group

variable

value

2338

MT2A

12.211493

2338

MT2A

11.306216

2338

MT2A

11.926226

2338

MT2A

10.974207

2338

MT2A

14.603994

2338

MT2A

6.904604

2338

MT2A

12.436719

2338

MT2A

9.798375

2338

MT2A

11.743673

2338

MT2A

11.240016

ggplot(data_var3_melt, aes(factor(variable),value))+ylab("Gene expression")+n geom_violin(aes(fill=factor(group)), stat="ydensity", position="dodge",n scale="width", trim=TRUE) +xlab(NULL) 一文读懂PCA分析（原理、算法、解释和可视化）

高颜值免费在线绘图工具升级版来了~~~

#ggplot(data_var3_melt, aes(factor(variable),value))+ylab("Gene expression")+ n #geom_quasirandom(aes(color=factor(group))) +xlab(NULL)n# 根据分组数目确定颜色变量ncolorA <- rainbow(length(unique(group)))n# 根据每个样品的分组信息获取对应的颜色变量ncolors <- colorA[as.factor(group)]n# 根据样品分组信息获得legend的颜色ncolorl <- colorA[as.factor(unique(group))]n# 获得PCH symbol列表npch_l <- as.numeric(as.factor(unique(group)))n# 产生每个样品的pch symbolnpch <- pch_l[as.factor(group)]nscatterplot3d(data_var3_forPCA[,1:3], color=colors, pch=pch)nlegend(0,10, legend=levels(as.factor(group)), col=colorl, pch=pch_l, xpd=T, horiz=F, ncol=6) 一文读懂PCA分析（原理、算法、解释和可视化）

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

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

# Pay attention to the format of PCA input n# Rows are samples and columns are variablesndata4_use_log2_t <- t(data4_use_log2)n# Add group column for plottingndata4_use_log2_label <- as.data.frame(data4_use_log2_t)ndata4_use_log2_label$group <- groupn# By default, prcomp will centralized the data using mean.n# Normalize data for PCA by dividing each data by column standard deviation.n# Often, we would normalize data.n# Only when we care about the real number changes other than the trends,n# `scale` can be set to TRUE. n# We will show the differences of scaling and un-scaling effects.npca <- prcomp(data4_use_log2_t, scale=T)n# sdev: standard deviation of the principle components.n# Square to get variancenpercentVar <- pca$sdev^2 / sum( pca$sdev^2)n# To check what's in pcanprint(str(pca))n## List of 5n## $ sdev : num [1:301] 36.6 30.4 23.3 21.6 19.9 ...n## $ rotation: num [1:16482, 1:301] -0.01133 -0.01955 -0.00199 -0.00569 -0.0204 ...n## ..- attr(*, "dimnames")=List of 2n## .. ..$ : chr [1:16482] "MT2A" "ANXA1" "ARHGDIB" "RND3" ...n## .. ..$ : chr [1:301] "PC1" "PC2" "PC3" "PC4" ...n## $ center : Named num [1:16482] 6.5 5.47 5.43 4.23 4.1 ...n## ..- attr(*, "names")= chr [1:16482] "MT2A" "ANXA1" "ARHGDIB" "RND3" ...n## $ scale : Named num [1:16482] 5.62 4.96 4.78 4.13 3.98 ...n## ..- attr(*, "names")= chr [1:16482] "MT2A" "ANXA1" "ARHGDIB" "RND3" ...n## $ x : num [1:301, 1:301] -37.6 -28.6 -30.6 -43.7 -11.4 ...n## ..- attr(*, "dimnames")=List of 2n## .. ..$ : chr [1:301] "Hi_2338_1" "Hi_2338_2" "Hi_2338_3" "Hi_2338_4" ...n## .. ..$ : chr [1:301] "PC1" "PC2" "PC3" "PC4" ...n## - attr(*, "class")= chr "prcomp"n## NULL

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

library(ggfortify)nautoplot(pca, data=data4_use_log2_label, colour="group") + n xlab(paste0("PC1 (", round(percentVar[1]*100), "% variance)")) + n ylab(paste0("PC2 (", round(percentVar[2]*100), "% variance)")) + theme_bw() + n theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + n theme(legend.position="right") 一文读懂PCA分析（原理、算法、解释和可视化）

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

# 根据分组数目确定颜色变量ncolorA <- rainbow(length(unique(group)))n# 根据每个样品的分组信息获取对应的颜色变量ncolors <- colorA[as.factor(group)]n# 根据样品分组信息获得legend的颜色ncolorl <- colorA[as.factor(unique(group))]n# 获得PCH symbol列表npch_l <- as.numeric(as.factor(unique(group)))n# 产生每个样品的pch symbolnpch <- pch_l[as.factor(group)]npc <- as.data.frame(pca$x)nscatterplot3d(x=pc$PC1, y=pc$PC2, z=pc$PC3, pch=pch, color=colors, n xlab=paste0("PC1 (", round(percentVar[1]*100), "% variance)"), n ylab=paste0("PC2 (", round(percentVar[2]*100), "% variance)"), n zlab=paste0("PC3 (", round(percentVar[3]*100), "% variance)"))nlegend(-3,8, legend=levels(as.factor(group)), col=colorl, pch=pch_l, xpd=T, horiz=F, ncol=6) 一文读懂PCA分析（原理、算法、解释和可视化）

PCA的实现原理

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

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

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

mat <- as.data.frame(matrix(rnorm(20,0,1), nrow=4))ncolnames(mat) <- paste0("Gene_", letters[1:5])nrownames(mat) <- paste0("Samp_", 1:4)nmatn## Gene_a Gene_b Gene_c Gene_d Gene_en## Samp_1 -1.2207554 0.7608654 -0.2921542 0.2781680 -0.5515104n## Samp_2 -0.3855339 -1.3785382 1.1008845 0.5385477 0.1083430n## Samp_3 0.5828754 0.7443152 0.2221331 1.0248715 0.2456042n## Samp_4 -1.6524972 -1.6259796 0.8235352 -0.7186743 0.3295052

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

# mean-centering data for columnsn# Get mean-value matrix firstnmat_mean_norm <- mat - rep(colMeans(mat),rep.int(nrow(mat),ncol(mat)))nmat_mean_normn## Gene_a Gene_b Gene_c Gene_d Gene_en## Samp_1 -0.5517776 1.135700 -0.7557539 -0.002560247 -0.58449593n## Samp_2 0.2834438 -1.003704 0.6372848 0.257819506 0.07535752n## Samp_3 1.2518532 1.119149 -0.2414665 0.744143242 0.21261872n## Samp_4 -0.9835194 -1.251145 0.3599356 -0.999402500 0.29651969n# mean-centering using scale for columnsnscale(mat, center=T, scale=F)n## Gene_a Gene_b Gene_c Gene_d Gene_en## Samp_1 -0.5517776 1.135700 -0.7557539 -0.002560247 -0.58449593n## Samp_2 0.2834438 -1.003704 0.6372848 0.257819506 0.07535752n## Samp_3 1.2518532 1.119149 -0.2414665 0.744143242 0.21261872n## Samp_4 -0.9835194 -1.251145 0.3599356 -0.999402500 0.29651969n## attr(,"scaled:center")n## Gene_a Gene_b Gene_c Gene_d Gene_e n## -0.66897778 -0.37483430 0.46359965 0.28072824 0.03298553

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

# median-centering data for columnsnmat_median_norm <- mat - rep(apply(mat,2,median),rep.int(nrow(mat),ncol(mat)))nmat_mean_normn## Gene_a Gene_b Gene_c Gene_d Gene_en## Samp_1 -0.5517776 1.135700 -0.7557539 -0.002560247 -0.58449593n## Samp_2 0.2834438 -1.003704 0.6372848 0.257819506 0.07535752n## Samp_3 1.2518532 1.119149 -0.2414665 0.744143242 0.21261872n## Samp_4 -0.9835194 -1.251145 0.3599356 -0.999402500 0.29651969

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

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

apply(mat,2,var)n## Gene_a Gene_b Gene_c Gene_d Gene_e n## 0.9730820 1.7050318 0.3883852 0.5396773 0.1601483

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

cov(mat)n## Gene_a Gene_b Gene_c Gene_d Gene_en## Gene_a 0.97308195 0.5734631 -0.01954727 0.66299331 0.10613531n## Gene_b 0.57346308 1.7050318 -0.73950788 0.60717437 -0.29082853n## Gene_c -0.01954727 -0.7395079 0.38838520 -0.12438895 0.18171565n## Gene_d 0.66299331 0.6071744 -0.12438895 0.53967732 -0.03906622n## Gene_e 0.10613531 -0.2908285 0.18171565 -0.03906622 0.16014830

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

# Covariance matrix for Mean normalized matrixncov(mat_mean_norm)n## Gene_a Gene_b Gene_c Gene_d Gene_en## Gene_a 0.97308195 0.5734631 -0.01954727 0.66299331 0.10613531n## Gene_b 0.57346308 1.7050318 -0.73950788 0.60717437 -0.29082853n## Gene_c -0.01954727 -0.7395079 0.38838520 -0.12438895 0.18171565n## Gene_d 0.66299331 0.6071744 -0.12438895 0.53967732 -0.03906622n## Gene_e 0.10613531 -0.2908285 0.18171565 -0.03906622 0.16014830n# Covariance matrix for Mean normalized matrix n# crossprod: matrix multiplicationncrossprod(as.matrix(mat_mean_norm)) / (nrow(mat_mean_norm)-1)n## Gene_a Gene_b Gene_c Gene_d Gene_en## Gene_a 0.97308195 0.5734631 -0.01954727 0.66299331 0.10613531n## Gene_b 0.57346308 1.7050318 -0.73950788 0.60717437 -0.29082853n## Gene_c -0.01954727 -0.7395079 0.38838520 -0.12438895 0.18171565n## Gene_d 0.66299331 0.6071744 -0.12438895 0.53967732 -0.03906622n## Gene_e 0.10613531 -0.2908285 0.18171565 -0.03906622 0.16014830n# Use %*% for matrix multiplication (slower)nt(as.matrix(mat_mean_norm)) %*% as.matrix(mat_mean_norm) / (nrow(mat_mean_norm)-1)n## Gene_a Gene_b Gene_c Gene_d Gene_en## Gene_a 0.97308195 0.5734631 -0.01954727 0.66299331 0.10613531n## Gene_b 0.57346308 1.7050318 -0.73950788 0.60717437 -0.29082853n## Gene_c -0.01954727 -0.7395079 0.38838520 -0.12438895 0.18171565n## Gene_d 0.66299331 0.6071744 -0.12438895 0.53967732 -0.03906622n## Gene_e 0.10613531 -0.2908285 0.18171565 -0.03906622 0.16014830

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

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

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

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

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

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

前面提到，新变量(Ym, k)是原始变量(Xm, n)(原始变量的协方差矩阵为(Cn, n))的线性组合，那么假设我们找到了这么一个线性组合(命名为特征矩阵(Pn, k))，得到一组新变量Ym, k = Xm, nPn, k，并且新变量的协方差矩阵(Dk, k)为对角阵。那么这个特征矩阵(Pn, k)需要符合什么条件呢？

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

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

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

假设An, n为n阶对称阵，如存在λ和非零向量x，使得A**x = λ**x，则称λ为矩阵An, n的特征值，非零向量x为为矩阵An, n对应于特征值λ的特征向量。

根据这个定义可以得出(An, n − λ**E)x = 0，由于x为非零向量，所以行列式|A − λ**E| = 0。

由此求解出n个根λ1, λ2, …, λ3就是矩阵A的特征值。

回顾下行列式的计算：

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

简单的PCA实现

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

library(knitr)nkable(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

5.99

5.35

5.14

4.66

5.31

5.05

4.92

5.12

4.88

4.79

5.11

4.8

…

b47

20.78

4.95

4.91

b48

20.5

4.9

5.11

b49

20.46

4.85

4.95

b50

19.94

4.87

5.18

标准化数据

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

Normalized expression for 3 genes in 100 samples

Gene1

Gene2

Gene3

-0.85

1.83

-0.96

-1.03

1.66

-0.98

-0.99

0.72

-1.01

0.67

-1.03

…

b47

1.09

-0.11

-1.01

b48

1.06

-0.38

-0.97

b49

1.05

-0.61

-1

b50

0.98

-0.52

-0.95

计算协方差矩阵

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

Covariance matrix for 3 genes in 100 samples

Gene1

Gene2

Gene3

Gene1

1.0000000

0.0255834

-0.0087919

Gene2

0.0255834

1.0000000

0.0553298

Gene3

-0.0087919

0.0553298

1.0000000

求解特征值和特征向量

data3_center_scale_cov_eigen <- eigen(data3_center_scale_cov)n# 特征值，从大到小排序ndata3_center_scale_cov_eigen$valuesn## [1] 1.0580005 1.0066389 0.9353605n# 特征向量, 每一列为对应特征值的特征向量ndata3_center_scale_cov_eigen$vectorsn## [,1] [,2] [,3]n## [1,] 0.2192050 0.90784568 0.3574429n## [2,] 0.7223522 0.09525968 -0.6849327n## [3,] 0.6558631 -0.40834033 0.6349030

产生新的矩阵

pc_select = 3nlabel = paste0("PC",c(1:pc_select))ndata3_new <- data3_center_scale %*% data3_center_scale_cov_eigen$vectors[,1:pc_select]ncolnames(data3_new) <- labelnkable(headTail(data3_new), booktabs=T, n 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

0.51

-0.21

-2.17

0.33

-0.37

-2.12

-0.36

-0.42

-1.49

-0.41

-0.43

-1.47

…

b47

-0.5

1.39

-0.17

b48

-0.68

1.32

0.03

b49

-0.86

1.31

0.16

b50

-0.79

1.23

0.1

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

#library(scatterplot3d)ncolorl <- c("#E69F00", "#56B4E9")n# Extract same number of colors as the Group and same Group would have same color.ncolors <- colorl[as.numeric(data3$Group)]n# 1 row 2 columnsnpar(mfrow=c(1,2))nscatterplot3d(data3[,1:3], color=colors, angle=55, pch=16, main="Original data")nlegend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)nscatterplot3d(data3_new, color=colors,angle=55, pch=16, main="Principle components")nlegend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T) 一文读懂PCA分析（原理、算法、解释和可视化）

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

利用prcomp进行主成分分析

pca_data3 <- prcomp(data3[,1:3], center=TRUE, scale=TRUE)n#Show whats in the result returned by prcompnstr(pca_data3)n## List of 5n## $ sdev : num [1:3] 1.029 1.003 0.967n## $ rotation: num [1:3, 1:3] 0.2192 0.7224 0.6559 -0.9078 -0.0953 ...n## ..- attr(*, "dimnames")=List of 2n## .. ..$ : chr [1:3] "Gene1" "Gene2" "Gene3"n## .. ..$ : chr [1:3] "PC1" "PC2" "PC3"n## $ center : Named num [1:3] 12.46 4.98 10n## ..- attr(*, "names")= chr [1:3] "Gene1" "Gene2" "Gene3"n## $ scale : Named num [1:3] 7.602 0.202 5.06n## ..- attr(*, "names")= chr [1:3] "Gene1" "Gene2" "Gene3"n## $ x : num [1:100, 1:3] 0.506 0.331 -0.36 -0.409 -0.27 ...n## ..- attr(*, "dimnames")=List of 2n## .. ..$ : chr [1:100] "a1" "a2" "a3" "a4" ...n## .. ..$ : chr [1:3] "PC1" "PC2" "PC3"n## - attr(*, "class")= chr "prcomp"n# 新的数据，与前面计算的抑制ndata3_pca_new <- pca_data3$xnkable(headTail(data3_pca_new), booktabs=T, n 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

0.51

0.21

-2.17

0.33

0.37

-2.12

-0.36

0.42

-1.49

-0.41

0.43

-1.47

…

b47

-0.5

-1.39

-0.17

b48

-0.68

-1.32

0.03

b49

-0.86

-1.31

0.16

b50

-0.79

-1.23

0.1

# 特征向量，与我们前面计算的一致(特征向量的符号是任意的)npca_data3$rotationn## PC1 PC2 PC3n## Gene1 0.2192050 -0.90784568 0.3574429n## Gene2 0.7223522 -0.09525968 -0.6849327n## Gene3 0.6558631 0.40834033 0.6349030

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

#library(scatterplot3d)ncolorl <- c("#E69F00", "#56B4E9")n# Extract same number of colors as the Group and same Group would have same color.ncolors <- colorl[as.numeric(data3$Group)]n# 1 row 2 columnsnpar(mfrow=c(1,2))nscatterplot3d(data3_new, color=colors,angle=55, pch=16, main="PCA by steps")nlegend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)nscatterplot3d(data3_pca_new, color=colors,angle=55, pch=16, main="PCA by prcomp")nlegend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T) 一文读懂PCA分析（原理、算法、解释和可视化）

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

自定义PCA计算函数

ct_PCA <- function(data, center=TRUE, scale=TRUE){n data_norm <- scale(data, center=center, scale=scale)n data_norm_cov <- crossprod(as.matrix(data_norm)) / (nrow(data_norm)-1)n data_eigen <- eigen(data_norm_cov)n rotation <- data_eigen$vectorsn label <- paste0('PC', c(1:ncol(rotation)))n colnames(rotation) <- labeln sdev <- sqrt(data_eigen$values)n data_new <- data_norm %*% rotationn colnames(data_new) <- labeln ct_pca <- list('rotation'=rotation, 'x'=data_new, 'sdev'=sdev)n return(ct_pca)n}

比较有无scale对聚类的影响，从图中可以看到，如果不对数据进行scale处理，样品的聚类结果更像原始数据，本身数值大的基因对主成分的贡献会大。如果关注的是每个变量自身的实际方差对样品分类的贡献，则不应该SCALE；如果关注的是变量的相对大小对样品分类的贡献，则应该SCALE，以防数值高的变量导入的大方差引入的偏见。

data3_pca_noscale_step = ct_PCA(data3[,1:3], center=TRUE, scale=FALSE)n# 特征向量ndata3_pca_noscale_step$rotationn## PC1 PC2 PC3n## [1,] 0.9999446062 0.010502677 -0.0006915646n## [2,] 0.0006682513 0.002223103 0.9999973056n## [3,] -0.0105041862 0.999942374 -0.0022159613n# 新变量ndata3_pca_noscale_pc <- data3_pca_noscale_step$xn#library(scatterplot3d)ncolorl <- c("#E69F00", "#56B4E9")n# Extract same number of colors as the Group and same Group would have same color.ncolors <- colorl[as.numeric(data3$Group)]n# 1 row 2 columnsnpar(mfrow=c(2,2))nscatterplot3d(data3[,c(1,3,2)], color=colors, angle=55, pch=16, main="Original data")nlegend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)nscatterplot3d(data3_pca_noscale_pc, color=colors,angle=55, pch=16, main="PCA (no scale)")nlegend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)nscatterplot3d(data3_center_scale[,c(1,3,2)], color=colors, angle=55, pch=16, n main="Original data (scale)")nlegend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T)nscatterplot3d(data3_new, color=colors,angle=55, pch=16, main="PCA (scale)")nlegend("top", legend=levels(data3$Group), col=colorl, pch=16, xpd=T, horiz=T) 一文读懂PCA分析（原理、算法、解释和可视化）

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

PCA结果解释

prcomp函数会返回主成分的标准差、特征向量和主成分构成的新矩阵。接下来，探索下不同主成分对数据差异的贡献和主成分与原始变量的关系。

主成分的平方为为特征值,其含义为每个主成分可以解释的数据差异，计算方式为eigenvalues = (pca$sdev)^2
每个主成分可以解释的数据差异的比例为percent_var = eigenvalues*100/sum(eigenvalues)
可以使用summary(pca)获取以上两条信息。

这两个信息可以判断主成分分析的质量：

主成分分析简介

主成分分析的意义

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

PCA的实现原理

简单的PCA实现

PCA结果解释

相关文章

央媒谈“文盲演员”引发饭圈甩锅大战

野生大熊猫深夜遛达 镜头前呆萌打卡

张萌谈演员到底需不需要有文化

情人节多地迎领证高峰

婚房烂尾10年业主住进毛坯房

多地首套房贷利率降至4%以下 为什么急于提前还贷

野生大熊猫深夜遛达镜头前呆萌打卡

多地首套房贷利率降至4%以下为什么急于提前还贷