PCA分析（Principal Component Analysis ）是一重数据降维的分析方式，它可以将复杂的数据通过数学运算变得更加简单。比如，在四维空间中我们要考量时间的影响，正如蚂蚁不懂生活在三维世界的我们，更高维度的世界已经超越了我们的认知。你在实验过程中，一定也遇到过类似情况，面对复杂的数据无从下手，今天我会带大家运用全新的思维解决你在实验过程中遇到的困扰，下面让我们看看吧。
在收集大量数据后，你是不是也面临着找不到头绪的情况。在实验设计愈发多样化的今天，多样化的大数据无疑会为研究和开发提供更加丰富多彩的信息，但是也在一定程度上增加了数据分析的难度。更重要的是在很多条件下，许多变量之间可能存在相关性，这也为数据分析带来不小的挑战。如果对每个指标进行分析，分析往往是孤立的，不能完全利用数据间的相关性，因而盲目减少指标会损失很多有用的信息，从而产生错误的结论。那么如何在不影响结果的情况下，简化分析并做出更加正确的判断，正是我要给大家分享的。
比如我们测了200个单细胞的数据，通过基因表达水平的差异将这200个样品进行分类，我们需要考虑这些样品间每个基因存在的调控、拮抗或者协同等关系，还要考虑一些外界因素的干扰，比如一些housekeeping genes在每个样品中都有一致的表达情况，那么考虑它们对于解释样本差异的作用将毫无意义，所以我们可以通过去掉对结果没有影响的基因来达到简化运算的目的。
医学领域中，我们也经常用PCA图来进行疾病危险因素分析，肠道菌群聚类分析，推断肿瘤亚群之间的进化关系，还用它来观察样本的分组、趋势、剔除异常数据，找到数据中的隐藏模式，找到相关变量。
R语言中可以进行PCA分析的主要有rda()、prcomp()、princomp() 、PCA() 、dudi.pca() 、epPCA() 等包；对于分析结果可视化，factoextra包封装了包括分析结果提取和基于ggplot2的数据可视化的函数。
理论知识我们都知道了，那么在R中如何实现PCA的分析呢，
如果我们考虑 1 个基因，可以用以下代码实现：
count <- 100
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)
# Add additional column to data only for plotting
cy_data$Y <- rep(0,count*2)
kable(headTail(cy_data), booktabs=TRUE, caption="Expression profile for Gene1 in 100 samples")
从上图可以看出，200 个样品根据 Gene1 表达量的不同在横轴上被被分为了 2 类，可以看做是在线性水平的分类。
如果我们考虑 2 个基因，可以用以下代码实现：
count <- 100
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")
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 <- 100
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")
从上图可以看出，200 个样品根据 Gene1、Gene2 和 Gene3 的表达量的不同在坐标轴上被被分为了 4 类，可以看做是
立体空间的分类。而且在这个例子中，我们可以很容易的看出 Gene1 和 Gene3 对样品分类的贡献要比 Gene2 大。
从上图可以看出，200 个样品根据 Gene1 和 Gene2 的表达量的不同在坐标轴上被被分为了 2 类，在这个例子中，我们可以很容易的看出 Gene1 对样品分类的贡献要比 Gene2 大，因为 Gene1 在样品间的表达差异大。
假如我们有 1 个基因，可以在线性层面对样本进行分类；如果我们有 2 个基因，可以在一个平面对样本进行分类；如果我们有 3 个基因，可以在一个立体空间对样本进行分类；如果有更多的基因，比如说 n个，那么每个样品就是 n 维空间的一个点，则很难在图形上展示样品的分类关系。利用 PCA 分析，我们可以选取贡献最大的 2 个或 3 个主成分作为数据代表用以可视化。这比直接选取三个表达变化最大的基因更能反映样品之间的差异。（利用 Pearson 相关系数对样品进行聚类在样品数目比较少时是一个解决办法）。
好了，今天的分享就到这里了，欢迎各位小伙伴留言讨论哦。
推荐阅读：
<u>四维空间到底是什么？通俗又不失严谨的科普 - 知乎 (zhihu.com)</u>
<u>PCA(principal component analysis,主成分分析) - 简书 (jianshu.com)</u>
