PCA in R

专题介绍：R是一种广泛用于数据分析和统计计算的强大语言，于上世纪90年代开始发展起来。得益于全世界众多爱好者的无尽努力，大家继而开发出了一种基于R但优于R基本文本编辑器的R Studio（用户的界面体验更好）。也正是由于全世界越来越多的数据科学社区和用户对R包的慷慨贡献，让R语言在全球范围内越来越流行。其中一些R包，例如MASS，SparkR， ggplot2，使数据操作,可视化和计算功能越来越强大。R是用于统计分析、绘图的语言和操作环境。R是属于GNU系统的一个自由、免费、源代码开放的软件，它是一个用于统计计算和统计制图的优秀工具。R作为一种统计分析软件，是集统计分析与图形显示于一体的。它可以运行于UNIX、Windows和Macintosh的操作系统上，而且嵌入了一个非常方便实用的帮助系统，相比于其他统计分析软件，R的学术性开发比较早，适合生物学和医学等学术学科的科研人员使用。

PCA in R

PCA in R 介绍如下内容：

1 怎么使用prcomp()函数做PCA

2 使用base graph 和 ggplot2对PCA做可视化

3 每个主成分的方差贡献度

4 通过loading scores决定变量的影响效应

人造数据集生成

人造数据集生成代码

# R包
library(tidyverse)

# 构建人造数据集
data_matrix <- matrix(nrow = 100, ncol = 10)
colnames(data_matrix) <- c(paste0("wt", 1:5, sep=""), paste0("ko", 1:5, sep=""))
rownames(data_matrix) <- paste0("gene", 1:100, sep="")

for(i in 1:100) {
  wt_values = rpois(5, lambda = sample(x=10:1000, size = 1))
  ko_values = rpois(5, lambda = sample(x=10:1000, size = 1))

  data_matrix[i, ] = c(wt_values, ko_values)
}

data_matrix %>% head

利用prcomp函数计算获得主成分

计算主成分的时候，行表示观察，列表式变量。因为人造数据集，行是变量，所以做主成分分析时，需要先对数据集进行转置操作。

参考代码

# 计算PCA
# 行为样本，列为变量
pca <- prcomp(t(data_matrix), scale =TRUE)

主成分可视化分析

绘制第一主成分和第二主成分的可视化。

方式一：利用base graph

代码

plot(pca$x[, 1], pca$x[, 2])

方式二：利用ggplot2

代码

pca_data <- data.frame(Sample=rownames(pca$x), X = pca$x[, 1], Y = pca$x[, 2])
ggplot(data = pca_data, aes(x = X, y = Y, label = Sample)) +
  geom_text() +
  xlab(paste("PC1 -", pca_var_per[1], "%", sep = "")) +
  ylab(paste("PC2 -", pca_var_per[2], "%", sep = "")) +
  theme_bw() +
  ggtitle("第一主成分和第二主成分的可视化")

各个主成分的方差贡献度

计算各个主成分的方差贡献度，代码如下：

pca_var <- pca$sdev^2
pca_var_per <- round(pca_var / sum(pca_var) * 100, 1)

把贡献度做成柱状图

代码

barplot(pca_var_per, main = "Scree plot", xlab = "Principal Component", ylab = "Percent Variation")

通过loading_scores获取关键变量

根据loading_scores的值获得各个变量的效应

loading_scores <- pca$rotation[, 1]
gene_scores <- abs(loading_scores)
gene_score_ranked <- sort(gene_scores, decreasing = TRUE)
top_10_genes <- names(gene_score_ranked[1:10])
top_10_genes
pca$rotation[top_10_genes, 1]