526互联

R : PCA 主成分分析

发布时间 2023-12-18 20:46:54作者: 王哲MGG_AI

主成分分析

rm (list = ls ()) 
library(vegan) 
library(tidyverse) 
library(ggalt) 
library(car)
library(ggforce)
library(ggpubr)
library(patchwork) 

# 2. 定义所需的函数。
pairwise.adonis1 <- function(x, factors, p.adjust.m) { #定义了一个名为pairwise.adonis1的函数，该函数接受三个参数：x，factors和p.adjust.m
  x = as.matrix(x)
  co = as.matrix(combn(unique(factors), 2))
  pairs <- F.Model <- R2 <- p.value <- c()
  for (elem in 1:ncol(co)) { 
    ad = adonis(x[factors %in% c(as.character(co[1, elem]), as.character(co[2, elem])), 
                  factors %in% c(as.character(co[1, elem]), as.character(co[2, elem]))] ~ 
                  factors[factors %in% c(as.character(co[1, elem]), as.character(co[2, elem]))], permutations = 999) 
    pairs <- c(pairs, paste(co[1, elem], 'vs', co[2, elem])) 
    F.Model <- c(F.Model, ad$aov.tab[1, 4]) 
    R2 <- c(R2, ad$aov.tab[1, 5]) 
    p.value <- c(p.value, ad$aov.tab[1, 6]) 
  }
  
  # p值调整
  p.adjusted = p.adjust(p.value, method = p.adjust.m) 
  pairw.res = data.frame(pairs, F.Model, R2, p.value, p.adjusted) 
  return(pairw.res)
}

# 3. 读取和处理数据。
setwd("C:\\Users\\Administrator\\Desktop\\PCA") 
otu <- read.table("./otu_table.txt", row.names = 1, sep = "\t", header = TRUE) %>% as.data.frame()
map <- read.table("./metadata3-1.txt", sep = "\t", header = TRUE) 
colnames(map)[1] <- "ID" 
row.names(map) <- map$ID 
idx <- rownames(map) %in% colnames(otu) 
map1 <- map[idx,]
otu <- otu[, rownames(map1)]

# 4. 进行adonis分析并计算统计值。
bray_curtis <- vegan::vegdist(t(otu), method = "bray", na.rm = TRUE)
ado <- adonis(bray_curtis ~ map1$Group, permutations = 999, method = "bray") #进行了Adonis分析。
R2_value <- round(as.data.frame(ado$aov.tab[5])[1, 1], 3)
p_v_value <- as.data.frame(ado$aov.tab[6])[1, 1]
title <- paste("Adonis:R^2 = ", R2_value, " P_value = ", p_v_value, sep = "")

# 5. 绘制PCA图。
#otu_centered <- scale(t(otu), scale = FALSE) # 对数据进行中心化处理
otu_standardized <- scale(t(otu), scale = TRUE) # 对数据进行标准化处理
pca <- prcomp(otu_standardized) # 使用prcomp函数进行了主成分分析（PCA）
#pca <- prcomp(t(otu)) # 使用prcomp函数进行了主成分分析（PCA）
summary_pca <- summary(pca) # 获取PCA的详细结果
points <- as.data.frame(pca$x) %>% dplyr::rename(x = "PC1", y = "PC2") # 将PCA的结果转换为数据框，并重命名了列名.
# pca$x是PCA结果中的坐标，dplyr::rename(x = "PC1", y = "PC2")将列名"PC1"和"PC2"改为"x"和"y"。
eig <- pca$eig
points <- cbind(points, map1[match(rownames(points), map1$ID),]) # 将map1数据框中的元数据添加到了PCA的结果中。

n <- 0.85
colors <- c("d0"="#00BFFF","d0.5"="#00BFFF","d1"="#00BFFF","d3"="#00BFFF","d5"="#00BFFF","d8"="#00BFFF","d10"="#00BFFF","W0"="#FF4500","W0.5"="#FF4500","W1"="#FF4500","W3"="#FF4500","W5"="#FF4500","W8"="#FF4500","W10"="#FF4500")
# 定义了颜色和形状的映射关系，用于后续的可视化。
shapes <- c("d0"=21, "d0.5"=21, "d1"=21, "d3"=21, "d5"=21,"d8"=21,"d10"=21,
            "W0"=21, "W0.5"=21, "W1"=21, "W3"=21, "W5"=21,"W8"=21,"W10"=21)
levels_order <- c("d0", "d0.5", "d1", "d3", "d5", "d8", "d10","W0","W0.5","W1","W3","W5","W8","W10") #定义顺序
points$Group <- factor(points$Group, levels = levels_order)
# 在ggplot中使用这些形状
p1 <- ggplot(points, aes(x = x, y = y, fill = Group, shape = Group)) + 
  geom_point(alpha = .7, size = 6) + 
  scale_shape_manual(values = shapes) + 
  scale_fill_manual(values = colors) +
  labs(x = paste("PC1 (", format(summary_pca$importance[2, 1] * 100, digits = 4), "%)", sep = ""),
       y = paste("PC2 (", format(summary_pca$importance[2, 2] * 100, digits = 4), "%)", sep = ""), title = title) +
  geom_mark_ellipse(aes(fill = Group, label = Group), alpha = 0.1, color = "grey", linetype = 3) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.text = element_text(color = "black", size = 9),
        axis.title = element_text(size = 20), # 改变坐标轴标题的大小
        legend.text = element_text(size = 14), # 改变图例文字的大小
        legend.key.size = unit(1, "cm")) +# 改变图例图形的大小
  coord_cartesian(xlim = c(-max(abs(points$x)) * 1.1, max(abs(points$x)) * 1.1), ylim = c(-max(abs(points$y)) * 1.1, max(abs(points$y)) * 1.1)) + # 改变x轴和y轴的范围
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") + # 添加垂直虚线
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") # 添加水平虚线

# 显示绘制的图
p1
# 6. 输出pairwise adonis结果。
pair_bray_adonis <- pairwise.adonis1(bray_curtis, map1$Group, p.adjust.m = "bonferroni") #进行了成对的Adonis分析
# 存储为文本文件
write.table(as.data.frame(pair_bray_adonis), "table.txt", sep = "\t", quote = FALSE, row.names = FALSE) #将Adonis分析的结果保存为了一个文本文件。
tab2 <- ggtexttable(pair_bray_adonis, rows = NULL)
p2 <- tab2
p2
ggsave(filename = "PCA_plot3.png", plot = p1, width = 12, height = 10, units = "in", dpi = 600)

1. 初始化和加载库

  • rm(list = ls())：清除R环境中的所有变量，以便开始新的分析。
  • library(...)：加载所需的R包。其中，vegan用于生态多样性分析，tidyverse提供了数据处理和可视化的强大工具，ggalt、car、ggforce、ggpubr和patchwork用于数据可视化和图形展示。

2. 定义函数 pairwise.adonis1

这个函数用于执行成对的Adonis分析（一种非参数多元方差分析方法）。它接收三个参数：数据矩阵x、因子factors和p值调整方法p.adjust.m。函数内部，它首先将数据转换为矩阵格式，然后对每一对唯一的因子组合进行Adonis分析，最后返回一个包含分析结果的数据框。

3. 数据读取和处理

  • setwd()：设置工作目录。
  • read.table()：读取OTU表和元数据文件。OTU表包含样本和操作分类单元（OTU）的计数，元数据文件包含样本信息。
  • dplyr::rename()、dplyr::filter()等：对数据进行处理，使其适合后续分析。

4. Adonis分析和PCA

  • 使用vegan::vegdist()计算基于Bray-Curtis距离的距离矩阵。
  • 进行Adonis分析，以测试群落组成是否因Group而异。
  • 进行PCA分析，这是一种降维技术，用于探索样本间的关系。

5. 数据可视化

  • 使用ggplot2和相关包创建PCA图，展示样本在主要成分上的分布。
  • 对PCA结果中的点进行美化和标记，以提供更多的可视化信息。

6. 输出结果

  • 执行pairwise.adonis1函数，以获得成对Adonis分析的结果。
  • 将结果保存为文本文件，并使用ggpubr创建一个表格，显示成对Adonis分析的结果。
  • 使用ggsave保存PCA图像。