学习利用prcomp()进行PCA（主成分分析）
利用base graghics和ggplot2来绘制PCA图
确定每个成分的描述的差异
找到影响最大的变量

建立模拟数据集

首先建立一个模拟数据集来进行后续分析
这个数据包干10个样本和100个基因，前五个样本是野生型“WT”,后五个样本是敲除“KO”
由于这是我们模拟出来用来练手的数据，暂且就将一百个基因命名为gene1, gene2, ..., gene100
然后生成随机数（这些随机数五个一组符合泊松分布）建立模拟数据集

data.matrix <- matrix(nrow=100, ncol=10)
colnames(data.matrix) <- c(
  paste("wt", 1:5, sep=""),
  paste("ko", 1:5, sep=""))
rownames(data.matrix) <- paste("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)
}
head(data.matrix)
dim(data.matrix)

生成的模拟数据集

利用prcomp()进行PCA

注意，prcomp默认列名为样本名，行名为基因（变量）名，最终目的计算列之间的关系

pca <- prcomp(t(data.matrix), scale=TRUE)

prcomp()的结果有三个部分x, sdev, rotation
x有10个，因为这里有10个样本，可以用来画图的主成分。在这里我们选x的前两行来画一个二维平面图，。
排名第一的主成分（PC1）能在是个样本中解释大部分变化，一般来说，当我们画二维PC图时，我们会选前两个成分（PC1, PC2)，但有时我们也会选第三个(PC3)或第四个(PC4)
sdev的平方就是SD，可以用来计算每个成分（PC）的解释度
rotation时每个样本在每个PC下基因（变量）的loading score,看影响最大的基因

利用base graghics简单作图

可以做一个柱状图来看每个成分占多大比例

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")

可以看出PC1解释了大多数差异
利用PC1, PC2画图

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

因为PC1最高，我们可以理解在PC1上的差别远大于PC2上的差别

用ggplot2画图

用ggplot2能画出来更好看的图
我们需要有三列分别是样本id，PC1和PC2，即每一行都有样本名，PC1和PC2

画图所需矩阵

library(ggplot2)
#制作输入矩阵id,PC1,PC2 
pca.data <- data.frame(Sample=rownames(pca$x),
  X=pca$x[,1],
  Y=pca$x[,2])
pca.data
 
#aes指定每个点的横纵坐标，并用label命名（在这里命名为id)
#xlab()和ylab()来命名横纵坐标轴
#theme_bw()规定背景为白色
#ggtitle来命名图标
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("My PCA Graph")

通过rotation来理解图片结果

以PC1为例

loading_scores <- pca$rotation[,1]
gene_scores <- abs(loading_scores) ## 我们需要知道的是最远的点，所以在既有正数也有复数的情况下取绝对值abs()
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] ##加上被我们用绝对值去掉的正负号

前十中分数为正的基因造成了五个“KO”样本集中在右侧的结果

前十中分数为负的基因造成了五个“WT”样本集中在右侧的结果

Statquest学习笔记|在R中进行PCA（主成分分析）