library(knitr)library(psych)library(reshape2)library(ggplot2)library(ggbeeswarm)library(scatterplot3d)library(useful)library(ggfortify)mat_show-function(matr){printmrow-function(x){ret-paste(paste(x,collapse=""),"\\\\\n")sprintf(ret)}align_str-paste0('{',paste0(rep('r',ncol(matr)),collapse=""),'}')format_mat-apply(matr,1,printmrow)add_env-paste("\\left[\\begin{array}",align_str,paste(format_mat,collapse=''),"\\{array}\\right]")return(add_env)}主成分分析简介主成分分析(PCA,principalcomponentanalysis)是一种数学降维方法,利用正交变换(orthogonaltransformation)把一系列可能线性相关的变量转换为一组线性不相关的新变量，也称为主成分，从而利用新变量在更小的维度下展示数据的特征。
主成分是原有变量的线性组合，其数目不多于原始变量。组合之后，相当于我们获得了一批新的观测数据，这些数据的含义不同于原有数据，但包含了之前数据的大部分特征，并且有着较低的维度，便于进一步的分析。
在空间上，PCA可以理解为把原始数据投射到一个新的坐标系统，第一主成分为第一坐标轴，它的含义代表了原始数据中多个变量经过某种变换得到的新变量的变化区间；第二成分为第二坐标轴，代表了原始数据中多个变量经过某种变换得到的第二个新变量的变化区间。这样我们把利用原始数据解释样品的差异转变为利用新变量解释样品的差异。
这种投射方式会有很多，为了最大限度保留对原始数据的解释，一般会用最大方差理论或最小损失理论，使得第一主成分有着最大的方差或变异数(就是说其能尽量多的解释原始数据的差异)；随后的每一个主成分都与前面的主成分正交，且有着仅次于前一主成分的最大方差(正交简单的理解就是两个主成分空间夹角为90°，两者之间无线性关联，从而完成去冗余操作)。
主成分分析的意义简化运算。
在问题研究中，为了全面系统地分析问题，我们通常会收集众多的影响因素也就是众多的变量。这样会使得研究更丰富，通常也会带来较多的冗余数据和复杂的计算量。
比如我们我们测序了100种样品的基因表达谱借以通过分子表达水平的差异对这100种样品进行分类。在这个问题中，研究的变量就是不同的基因。每个基因的表达都可以在一定程度上反应样品之间的差异，但某些基因之间却有着调控、协同或拮抗的关系，表现为它们的表达值存在一些相关性，这就造成了统计数据所反映的信息存在一定程度的冗余。另外假如某些基因如持家基因在所有样本中表达都一样，它们对于解释样本的差异也没有意义。这么多的变量在后续统计分析中会增大运算量和计算复杂度，应用PCA就可以在尽量多的保持变量所包含的信息又能维持尽量少的变量数目，帮助简化运算和结果解释。
去除数据噪音。
比如说我们在样品的制备过程中，由于不完全一致的操作，导致样品的状态有细微的改变，从而造成一些持家基因也发生了相应的变化，但变化幅度远小于核心基因(一般认为噪音的方差小于信息的方差）。而PCA在降维的过程中滤去了这些变化幅度较小的噪音变化，增大了数据的信噪比。
利用散点图实现多维数据可视化。
在上面的表达谱分析中，假如我们有1个基因，可以在线性层面对样本进行分类；如果我们有2个基因，可以在一个平面对样本进行分类；如果我们有3个基因，可以在一个立体空间对样本进行分类；如果有更多的基因，比如说n个，那么每个样品就是n维空间的一个点，则很难在图形上展示样品的分类关系。利用PCA分析，我们可以选取贡献最大的2个或3个主成分作为数据代表用以可视化。这比直接选取三个表达变化最大的基因更能反映样品之间的差异。（利用Pearson相关系数对样品进行聚类在样品数目比较少时是一个解决办法）
发现隐性相关变量。
我们在合并冗余原始变量得到主成分过程中，会发现某些原始变量对同一主成分有着相似的贡献，也就是说这些变量之间存在着某种相关性，为相关变量。同时也可以获得这些变量对主成分的贡献程度。对基因表达数据可以理解为发现了存在协同或拮抗关系的基因。
示例展示原始变量对样品的分类假设有一套数据集，包含100个样品中某一基因的表达量。如下所示，每一行为一个样品，每一列为基因的表达值。这也是做PCA分析的基本数据组织方式，每一行代表一个样品，每一列代表一组观察数据即一个变量。
count-50Gene1_a-rnorm(count,5,0.5)Gene1_b-rnorm(count,20,0.5)grp_a-rep('a',count)grp_b-rep('b',count)cy_(Gene1=c(Gene1_a,Gene1_b),Group=c(grp_a,grp_b))cy_(cy_data)label-c(paste0(grp_a,1:count),paste0(grp_b,1:count))(cy_data)-labellibrary(knitr)library(psych)geom_quasirandom:用于画JitterPlotxlim,ylim设定坐标轴的区间ggplot(cy_data,aes(Gene1,Y))+geom_quasirandom(aes(color=factor(Group)),groupOnX=FALSE)+theme(=c(0.5,0.7))+theme(=element_blank())+scale_fill_discrete(name="Group")+theme(=element_blank(),=element_blank(),=element_blank(),=element_blank())+ylim(-0.5,5)+xlim(0,25)那么如果有2个基因呢？
count-50Gene2_a-rnorm(count,5,0.2)Gene2_b-rnorm(count,5,0.2)cy_(Gene1=c(Gene1_a,Gene1_b),Gene2=c(Gene2_a,Gene2_b),Group=c(grp_a,grp_b))cy_(cy_data2)(cy_data2)-labelkable(headTail(cy_data2),booktabs=T,caption="ExpressionprofileforGene1andGene2in100samples")
ExpressionprofileforGene1andGene2in100samples
Gene1
Gene2
Group
a1
5.99
5.35
a
a2
4.66
5.31
a
a3
4.92
5.12
a
a4
4.79
5.11
a
…
…
…
NA
b47
20.78
4.95
b
b48
20.5
4.9
b
b49
20.46
4.85
b
b50
19.94
4.87
b
从下图可以看出，100个样品根据Gene1和Gene2的表达量的不同在坐标轴上被被分为了2类，可以看做是在平面水平的分类。而且在这个例子中，我们可以很容易的看出Gene1对样品分类的贡献要比Gene2大，因为Gene1在样品间的表达差异大。
ggplot(cy_data2,aes(Gene1,Gene2))+geom_point(aes(color=factor(Group)))+theme(=c(0.5,0.9))+theme(=element_blank())+ylim(0,10)+xlim(0,25)
如果有3个基因呢？
count-50Gene3_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))(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)data3$(data3$Group)(data3)-labelkable(headTail(data3),booktabs=T,caption="Expressionprofilefor3genesin100samples")
Expressionprofilefor3genesin100samples
Gene1
Gene2
Gene3
Group
a1
5.99
5.35
5.14
a
a2
4.66
5.31
5.05
a
a3
4.92
5.12
4.88
a
a4
4.79
5.11
4.8
a
…
…
…
…
NA
b47
20.78
4.95
4.91
b
b48
20.5
4.9
5.11
b
b49
20.46
4.85
4.95
b
b50
19.94
4.87
5.18
b
从下图可以看出，100个样品根据Gene1、Gene2和Gene3的表达量的不同在坐标轴上被被分为了4类，可以看做是立体空间的分类。而且在这个例子中，我们可以很容易的看出Gene1和Gene3对样品分类的贡献要比Gene2大。
Geneexpressionmatrix
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的基因视为表达的基因），并把所有基因根据其在不同样品中表达值的变异系数排序。
data4_nonzero-data4[rowSums(data4)!=0,]data4_use-data4[apply(data4,1,function(row)sum(row=1)=5),]对于表达谱数据，因为涉及到PCR的指数扩增，一般会取log处理计算变异系数(标准差除以平均值)度量基因表达变化幅度根据变异系数排序计算中值绝对偏差(MAD,medianabsolutedeviation)度量基因表达变化幅度但是很小的观察值会引入很大的背景噪音，因此也意义不大。mads-apply(data4_use_log2,1,mad)data4_use_log2-data4_use_log2[rev(order(mads)),]转置矩阵使得每一行为一个样品，每一列为一组变量data_var3_forPCA-t(data_var3)dim(data_var3_forPCA)[1]3013kable(corner(data_var3_forPCA,r=10,c=5),booktabs=TRUE,caption="Atableofthe3mostvariablegenes")
Atableofthe3mostvariablegenes
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_split-strsplit(sample,"_")ggplot(data_var3_melt,aes(factor(variable),value))+ylab("Geneexpression")+根据分组数目确定颜色变量colorA-rainbow(length(unique(group)))根据样品分组信息获得leg的颜色colorl-colorA[(unique(group))]产生每个样品的pchsymbolpch-pch_l[(group)]scatterplot3d(data_var3_forPCA[,1:3],color=colors,pch=pch)leg(0,10,leg=levels((group)),col=colorl,pch=pch_l,xpd=T,horiz=F,ncol=6)我们看到图中的样品并没有按照预先设定的标签完全分开。当然我们也可以通过其他方法筛选变异最大的三个基因，最终的分类效果不会相差很大。因为不管怎么筛选，我们都只用到了3个基因的表达量。
假如我们把这个数据用PCA来分类，结果是怎样的呢？
Rowsaresamplesandcolumnsarevariablesdata4_use_log2_t-t(data4_use_log2)Bydefault,,wewouldnormalizedata.`scale`:'sinpcaprint(str(pca))Listof5$sdev:num[1:301]36.630.423.321.619.9$rotation:num[1:16482,1:301]-0.01133-0.01955-0.00199-0.00569-0.0204..-attr(*,"dimnames")=Listof2.$:chr[1:16482]"MT2A""ANXA1""ARHGDIB""RND3".$:chr[1:301]"PC1""PC2""PC3""PC4"$center:Namednum[1:16482]6.55.475.434.234.1..-attr(*,"names")=chr[1:16482]"MT2A""ANXA1""ARHGDIB""RND3"$scale:Namednum[1:16482]5.624.964.784.133.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")=Listof2.$: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(=element_blank(),=element_blank())+theme(="right")采用3个主成分获得的分类效果优于2个主成分，因为这样保留的原始信息更多。
根据每个样品的分组信息获取对应的颜色变量colors-colorA[(group)]获得PCHsymbol列表pch_((unique(group)))mean-centeringdataforcolumnsmean-centeringusingscaleforcolumnsscale(mat,center=T,scale=F)Gene_aGene_bGene_cGene_dGene_eSamp_1-0.55177761.135700-0.7557539-0.002560247-0.58449593Samp_20.2834438-1.0037040.63728480.2578195060.07535752Samp_31.25185321.119149-0.24146650.7441432420.21261872Samp_4-0.9835194-1.2511450.3599356-0.9994025000.29651969attr(,"scaled:center")Gene_aGene_bGene_cGene_dGene_
中位数中心化(mediancentering)：如果数据变换范围很大或有异常值，中位数标准化效果会更好。
CovariancematrixforMeannormalizedmatrixcov(mat_mean_norm)Gene_aGene_bGene_cGene_dGene_eGene_____:matrixmultiplicationcrossprod((mat_mean_norm))/(nrow(mat_mean_norm)-1)Gene_aGene_bGene_cGene_dGene_eGene_____特征值，从大到小排序data3_center_scale_cov_eigen$values[1]1.05800051.00663890.9353605library(scatterplot3d)colorl-c("56B4E9")1row2columnspar(mfrow=c(1,2))scatterplot3d(data3[,1:3],color=colors,angle=55,pch=16,main="Originaldata")leg("top",leg=levels(data3$Group),col=colorl,pch=16,xpd=T,horiz=T)scatterplot3d(data3_new,color=colors,angle=55,pch=16,main="Principlecomponents")leg("top",leg=levels(data3$Group),col=colorl,pch=16,xpd=T,horiz=T)Showwhatsintheresultreturnedbyprcompstr(pca_data3)Listof5$sdev:num[1:3]1.0291.0030.967$rotation:num[1:3,1:3]0.21920.72240.6559-0.9078-0.0953..-attr(*,"dimnames")=Listof2.$:chr[1:3]"Gene1""Gene2""Gene3".$:chr[1:3]"PC1""PC2""PC3"$center:Namednum[1:3]12.464.9810..-attr(*,"names")=chr[1:3]"Gene1""Gene2""Gene3"$scale:Namednum[1:3]7.6020.2025.06..-attr(*,"names")=chr[1:3]"Gene1""Gene2""Gene3"$x:num[1:100,1:3]0.5060.331-0.36-0.409-0.27..-attr(*,"dimnames")=Listof2.$:chr[1:100]"a1""a2""a3""a4".$:chr[1:3]"PC1""PC2""PC3"-attr(*,"class")=chr"prcomp"特征向量，与我们前面计算的一致(特征向量的符号是任意的)pca_data3$
比较手动实现的PCA与prcomp实现的PCA的聚类结果
E69F00","Ex[(data3$Group)]par(mfrow=c(1,1))
自定义PCA计算函数
ct_PCA-function(data,center=TRUE,scale=TRUE){data_norm-scale(data,center=center,scale=scale)data_norm_cov-crossprod((data_norm))/(nrow(data_norm)-1)data_eigen-eigen(data_norm_cov)rotation-data_eigen$vectorslabel-paste0('PC',c(1:ncol(rotation)))colnames(rotation)-labelsdev-sqrt(data_eigen$values)data_new-data_norm%*%rotationcolnames(data_new)-labelct_pca-list('rotation'=rotation,'x'=data_new,'sdev'=sdev)return(ct_pca)}data3_pca_noscale_step=ct_PCA(data3[,1:3],center=TRUE,scale=FALSE)新变量data3_pca_noscale_pc-data3_pca_noscale_step$xE69F00","Ex[(data3$Group)]par(mfrow=c(1,1))PCA结果解释
prcomp函数会返回主成分的标准差、特征向量和主成分构成的新矩阵。接下来，探索下不同主成分对数据差异的贡献和主成分与原始变量的关系。
主成分的平方为为特征值,其含义为每个主成分可以解释的数据差异，计算方式为eigenvalues=(pca$sdev)^2
每个主成分可以解释的数据差异的比例为percent_var=eigenvalues*100/sum(eigenvalues)
可以使用summary(pca)获取以上两条信息。
这两个信息可以判断主成分分析的质量：
成功的降维需要保证在前几个为数不多的主成分对数据差异的解释可以达到80-90%。
指导选择主成分的数目：
选择的主成分足以解释的总方差大于80%(方差比例碎石图)
从前面的协方差矩阵可以看到，自动定标(scale)的变量的方差为1(协方差矩阵对角线的值)。待选择的主成分应该是那些方差大于1的主成分，即其解释的方差大于原始变量（特征值碎石图，方差大于1，特征值也会大于1，反之亦然）。
鉴定核心变量和变量间的隐性关系:
原始变量与主成分的相关性VariablecorrelationwithPCs()=loadings*sdev
原始数据对主成分的贡献度^2/(^2)
在测试数据中，scale后，三个主成分对数据差异的贡献度大都在30%左右，而未scale的数据，三个主成分对数据差异的贡献度相差很大。这是因为三个基因由于自身表达量级所引起的方差的差异导致它们各自对数据的权重差异，从而使主成分偏向于数值大的变量。
PCA应用于测试数据前面用到一组比较大的测试数据集，并做了PCA分析，现在测试不同的处理对结果的影响。
首先回顾下我们用到的数据。
library("gridExtra")RowsaresamplesandcolumnsarevariablesAddgroupcolumnforplottingdata4_use_log2_label$group-groupkable(corner(data4_use_log2_label),digits=3,caption="Singlecellgeneexpressiondata")Singlecellgeneexpressiondata
MT2A
ANXA1
ARHGDIB
RND3
BLVRB
Hi_2338_1
12.211
11.837
9.543
9.762
9.162
Hi_2338_2
11.306
8.099
8.072
10.107
9.423
Hi_2338_3
11.926
10.689
9.721
9.388
9.694
Hi_2338_4
10.974
9.387
9.883
8.792
9.767
Hi_2338_5
14.604
10.375
9.970
8.815
9.350
比较对数运算和scale对样品分类的影响。
ct_pca_2d_plot-function(pca,data_with_label,labelName='group',title='PCA'){SquaretogetvariancepercentVar-pca$sdev^2/sum(pca$sdev^2)data[labelName]-factor(unlist(data[labelName]))level-length(unique(unlist(data_with_label[labelName])))shapes=(1:level)%%30Bydefault,,wewouldnormalizedata.`scale`注意事项一般说来，在PCA之前原始数据需要中心化（centering，数值减去平均值）。中心化的方法很多，除了平均值中心化（mean-centering）外，还包括其它更稳健的方法，比如中位数中心化等。
除了中心化以外，定标(Scale,数值除以标准差)也是数据前处理中需要考虑的一点。如果数据没有定标，则原始数据中方差大的变量对主成分的贡献会很大。数据的方差与其量级成指数关系，比如一组数据(1,2,3,4)的方差是1.67，而(10,20,30,40)的方差就是167,数据变大10倍，方差放大了100倍。
但是定标(scale)可能会有一些负面效果，因为定标后变量之间的权重就是变得相同。如果我们的变量中有噪音的话，我们就在无形中把噪音和信息的权重变得相同，但PCA本身无法区分信号和噪音。在这样的情形下，我们就不必做定标。
一般而言，对于度量单位不同的指标或是取值范围彼此差异非常大的指标不直接由其协方差矩阵出发进行主成分分析，而应该考虑对数据的标准化。比如度量单位不同，有万人、万吨、万元、亿元，而数据间的差异性也非常大，小则几十大则几万，因此在用协方差矩阵求解主成分时存在协方差矩阵中数据的差异性很大。在后面提取主成分时发现，只提取了一个主成分，而此时并不能将所有的变量都解释到，这就没有真正起到降维的作用。此时就需要对数据进行定标(scale)，这样提取的主成分可以覆盖更多的变量，这就实现主成分分析的最终目的。但是对原始数据进行标准化后更倾向于使得各个指标的作用在主成分分析构成中相等。对于数据取值范围不大或是度量单位相同的指标进行标准化处理后，其主成分分析的结果与仍由协方差矩阵出发求得的结果有较大区别。这是因为对数据标准化的过程实际上就是抹杀原有变量离散程度差异的过程，标准化后方差均为1，而实际上方差是对数据信息的重要概括形式，也就是说，对原始数据进行标准化后抹杀了一部分重要信息，因此才使得标准化后各变量在主成分构成中的作用趋于相等。因此，对同度量或是取值范围在同量级的数据还是直接使用非定标数据求解主成分为宜。
中心化和定标都会受数据中离群值（outliers）或者数据不均匀（比如数据被分为若干个小组）的影响，应该用更稳健的中心化和定标方法。
PCA也存在一些限制，例如它可以很好的解除线性相关，但是对于高阶相关性就没有办法了，对于存在高阶相关性的数据，可以考虑KernelPCA，通过Kernel函数将非线性相关转为线性相关，关于这点就不展开讨论了。另外，PCA假设数据各主特征是分布在正交方向上，如果在非正交方向上存在几个方差较大的方向，PCA的效果就大打折扣了。
参考资料PCA教程1
PCA文字化描述
pca1
ggplot2axis
scatterplot3D
稳健PCA
Datacentering
SampleRmarkdown
矩阵特征值，对称矩阵的对角化
Detailusageandvisualizationofprcompresult
ggplot2sidebysideplot
PCA主成分分析实战和可视化|附R代码和测试数据
用了这么多年的PCA可视化竟然是错的！！！
这篇Nature子刊文章的蛋白组学数据PCA分析竟花费了我两天时间来重现|附全过程代码
复现naturecommunicationPCA原图|代码分析（一）
