R 语言 ggplot2 PCA 主成分分析（虚拟数据集）

news2024/9/21 14:48:48

生成虚拟数据集

library(ggplot2)

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

head(data.matrix)

以上代码生成了100行基因，10列样本的矩阵
前五列命名wt开头+1-5，表示正常基因
后五列命名ko开头+1-5，表示缺少基因的样本（knock-out）

给每行基因都统一命名gene+1-100

head()函数默认查看前6行

在这里插入图片描述
现在只是定义了矩阵的shape和name，还没填充数值

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)

这段代码的作用是生成一个大小为100x10的数据矩阵data.matrix，其中前5列是"wt"（wild-type）样本的值，后5列是"ko"（knockout）样本的值。

在循环中，对于每个i的取值（从1到100），首先使用sample(x=10:1000, size = 1)从10到1000之间的整数中随机抽取一个数作为泊松分布的参数lambda。然后，使用rpois(5, lambda)函数生成一个具有泊松分布的随机数向量，其中每个元素表示一个基因在"wt"样本中的表达量。同样的过程也用于生成"ko"样本中的表达量。

最后，通过c(wt.values, ko.values)将"wt"和"ko"样本的表达量合并为一个长度为10的向量，并将其赋值给data.matrix的第i行。

用for循环给依次给1-100行的前五列和后五列赋值，填充值介于10-1000之间。

在这里插入图片描述
初始虚拟数据集创建完毕，接下来用prcomp()函数分析各样本之间关系。该函数默认情况下以基因为列，样本为行，和我们创建的矩阵互为转置，因此需要用到转置函数t()

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

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

pca.var <- pca$sdev^2
pca.var.per <- round(pca.var/sum(pca.var)*100, 1)

barplot(pca.var.per, main = "Screen Plot", xlab ="principal component", ylab = "percent variation")

prcomp(, scale = TRUE)表示对数据进行标准化
每一列表示一个基因所对应的10个样本，即一列只有十个数据
plot 生成一个2D的图，前两个主成分的散点图
在这里插入图片描述
pca.var 表示标准差的平方

pca.var.per 表示每个变量所占的百分比，保留小数点后一位
在这里插入图片描述
可以看到前两个成分所占比例最大，尤其是第一个成分

用 barplot 来直观每个成分所占比例

在这里插入图片描述

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("my pac graph")

先按 ggplot2 需要的方式格式化数据，x轴用第一个成分，y轴用第二个成分

在这里插入图片描述
可以发现数据分布在了两侧，我们用prcomp()调用负载得分loading scores 的参数rotation

通过分析pca$rotation，可以了解该主成分与哪些基因相关性较高，哪些基因对主成分的贡献较大。这对于解释主成分分析的结果和理解数据的结构和变化模式非常有帮助。

loading_scores <- pca$rotation[,1]

gene_score <- abs(loading_scores)
gene_score_ranked <- sort(gene_score,decreasing = TRUE)
top_10_genes <- names(gene_score_ranked[1:10])


pca$rotation[top_10_genes, 1]