本文首发于“生信大碗”公众号，转载请注明出处
1.Principal component analysis(PCA)
Principal component analysis (PCA)，即主成分分析，是最常用的线性降维方法。
主成分分析，通俗的讲，就是将一组数据消除冗余之后得到这组数据中的主要内容，用主要内容来代替原来的数据。例如原来的数据A包含n个变量，是一个k维的数据，那么通过消除冗余，我们将其变为数据B,包含m个变量的e维的数据，这个过程其实也就是一个降维的过程，降维后的维度e<k,且数据B可以被用来代替原始数据A。这里又提到了一个新的概念：降维，降维也就是减少数据中的变量，变量减少了维度也就减少了。这个化繁为简的过程中用到了比较复杂的算法，最终我们用得到的几个变量来代表整个数据，最终的变量也就是所谓的主成分！对于PCA的推导，有两种主流方法，一种是基于最小的投影距离，也就是样本点到超平面的最小距离；另一种是基于最大投影方差。具体的算法说起来就有点复杂了，大家了解即可。
对于PCA,举个简单的例子，如果我们要对授课老师的教学水平进行打分，那么就需要从：授课逻辑性、表达流畅度、课件完整性、课堂活跃度、学生的平时成绩等n多个方面进行综合考虑，那么如果我们现在要找出其中最主要的一些变量，也就是降维得到这些变量中的“主成分”。假如我们从n个变量降维成m个主成分，它们是由原有的变量线性组合而来，而不仅仅只是拿出来几个变量，这些主成分会依据方差的大小进行排序，也就是分别为主成分1（PC1）、主成分2（PC2）, 主成分3（PC3）,主成分4（PC4）等。每个主成分的方差在这一组变量中的总方差中所占的比例，即主成分的贡献度。一般我们仅关注贡献度前2或者前3的主成分，再进行可视化，得到二维或三维的PCA散点图。
在PCA图中，如果数个样本的点非常分散，或者点的连线距离长，就说明样本之间的相似度比较低，也就是差异性很大。如果数个样本的点聚集在一起，或者说点的连线距离短，就说明样本之间的相似度比较高，也就是差异性小。
例如，在RNA-seq分析中，我们获得了各样本中所有基因的表达量后，如果想比较样本之间在基因表达值的整体相似性或者差异程度，但是由于基因数量非常多，不能单独每一个基因去分析。因此要用到"降维"算法，将n维的表达矩阵降维成只有少数几个主成分的形式，以代表所有基因的整体表达情况，进而描述样本间的差异。
如下图所示，我们比较肝癌正常和肿瘤组织间DEG基因的表达差异，进行主成分分析（PCA)并进行可视化。如图，横坐标是PCA1，纵坐标是PCA2,百分比即该主成分的贡献度，两组中的样本在组内相互聚集，说明组内的样本数据非常相似，而组间也有较好的区分度。
2.绘制上图的代码展示
#加载包
install.packages("FactoMineR")
library(FactoMineR)
#加载差异基因
load("差异基因.rdata")
#加载表达矩阵
load("表达矩阵.rdata")
gene=rownames(data)
data1=mRNA_exp[gene,]
data1=t(data1)
#样本分组
normal=rownames(data1)[substr(rownames(data1),14,15)==11]
tumor=rownames(data1[substr(rownames(data1),14,15)!=11]
data1=data1[c(normal,tumor),]
group=data.frame(c(normal,tumor))
group$type=c(rep("normal",59),rep("tumor",526))
colnames(group)=c("sample","type")
#PCA分析
library(ggrepel)
pca <- PCA(data1, ncp = 2, scale.unit = TRUE, graph = TRUE)
pca_sample <- data.frame(pca$ind$coord[ ,1:2])
pca_sample$sample=row.names(pca_sample)
pca_eig1 <- round(pca$eig[1,2], 2)
pca_eig2 <- round(pca$eig[2,2],2 )
pca_sample <- merge(pca_sample, group,by="sample")
#画图
library(ggplot2)
plot <- ggplot(data = pca_sample, aes(x = Dim.1, y = Dim.2)) + geom_point(aes(color = type), size = 2) + scale_color_manual(values = c('red', 'green')) + theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) + labs(x = paste('PCA1:', pca_eig1, '%'), y = paste('PCA2:', pca_eig2, '%'), color = '')+ stat_ellipse(aes(color = type), level = 0.95, show.legend = FALSE)
plot
本文首发于“生信大碗”公众号，转载请注明出处
—END—
网友评论