# 基于CRAN安装R包，检测没有则安装p_list = c("tidyverse", "ggsci", "magrittr", "ggh4x", "rstatix", "ggsignif", "ggpubr",           "ggnewscale", "patchwork", "reshape2", "ggplot2")for(p in p_list){if (!requireNamespace(p)){install.packages(p)}    library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)}
# install.packages("devtools")# 基于github安装library(devtools)if(!requireNamespace("ggchicklet", quietly = TRUE))  install_github("hrbrmstr/ggchicklet")
# 加载R包 Load the packagesuppressWarnings(suppressMessages(library(ggchicklet)))suppressWarnings(suppressMessages(library(tidyverse)))suppressWarnings(suppressMessages(library(ggsci)))suppressWarnings(suppressMessages(library(magrittr)))suppressWarnings(suppressMessages(library(ggh4x)))suppressWarnings(suppressMessages(library(rstatix)))suppressWarnings(suppressMessages(library(ggsignif)))suppressWarnings(suppressMessages(library(ggpubr)))suppressWarnings(suppressMessages(library(ggnewscale)))suppressWarnings(suppressMessages(library(patchwork)))suppressWarnings(suppressMessages(library(reshape2)))suppressWarnings(suppressMessages(library(ggplot2)))

# 载入数据# Load datadata <- read.table("data/pathway.txt",header = TRUE,row.names = 1,sep = "\t")group <- read.table("data/group.txt",header = FALSE,sep = "\t")group = group[-1,]# 构建矩阵# Construct matrixrow=as.numeric(length(row.names(data))) col=as.numeric(length(colnames(data))) col_sum=rep(colSums(data), row)col_sum=matrix(col_sum, nrow = col, ncol = row)# 计算相对丰度# Calcualate relative abundancedata2=data/t(col_sum)colSums(data2)#>  Healthy01  Healthy02  Healthy03  Healthy04  Healthy05  Healthy06  Healthy07 #>          1          1          1          1          1          1          1 #>  Healthy08  Healthy09  Healthy10  Healthy11  Healthy12  Healthy13  Healthy14 #>          1          1          1          1          1          1          1 #>  Healthy15  Healthy16  Healthy17  Healthy18  Healthy19  Healthy20 Patients01 #>          1          1          1          1          1          1          1 #> Patients02 Patients03 Patients04 Patients05 Patients06 Patients07 Patients08 #>          1          1          1          1          1          1          1 #> Patients09 Patients10 Patients11 Patients12 Patients13 Patients14 Patients15 #>          1          1          1          1          1          1          1 #> Patients16 Patients17 Patients18 Patients19 Patients20 #>          1          1          1          1          1# 过滤掉平均丰度低于5%的功能分类# Filter out functional categories with an average abundance below 5%data = data2data <- data*100data <- data %>% filter(apply(data,1,mean) > 0.5)data <- t(data)data1 <- data.frame(data,group$V2)colnames(data1) <- c(colnames(data),"Group")data1$Group <- as.factor(data1$Group)# 不符合正态分布，使用非参检验# Does not conform to normal distribution, use non-parametric test# shapiro.test(data1$`Amino acid metabolism`)# shapiro.test(data1$`Biosynthesis of other secondary metabolites`)# shapiro.test(data1$Alistipes_putredinis)# diff <- data1 %>% #     select_if(is.numeric) %>%#     map_df(~ broom::tidy(t.test(. ~ Group,data = data1)), .id = 'var')# # diff$p.value <- p.adjust(diff$p.value,"bonferroni")# diff <- diff %>% filter(p.value < 0.05)# 非参检验# wilcox testlibrary(tidyverse)diff <- data1 %>%     select_if(is.numeric) %>%    map_df(~ broom::tidy(wilcox.test(. ~ Group,data = data1, conf.int = TRUE)), .id = 'var')diff$p.value <- p.adjust(diff$p.value,"BH")diff <- diff %>% filter(p.value < 0.05)# write.csv(diff, "results/Pathway_wixcox_test_0.05_abun_filter_0.05.csv")## 绘图数据构建(Data construction)## 左侧条形图(Left bar plot)abun.bar <- data1[,c(diff$var,"Group")] %>%     rstatix::gather(variable,value,-Group) %>%     group_by(variable,Group) %>%     summarise(Mean = mean(value))## 右侧散点图(Scatter plot on the right)diff.mean <- diff[,c("var","estimate","conf.low","conf.high","p.value")]diff.mean$Group <- c(ifelse(diff.mean$estimate >0,levels(data1$Group)[1],                            levels(data1$Group)[2]))diff.mean <- diff.mean[order(diff.mean$estimate,decreasing = TRUE),]## 左侧条形图(Left bar plot)cbbPalette <- c("#5ebcc2","#d9a285")abun.bar$variable <- factor(abun.bar$variable,levels = rev(diff.mean$var))p1 <- ggplot(abun.bar,aes(variable,Mean,fill = Group)) +    scale_x_discrete(limits = levels(diff.mean$var)) +    coord_flip() +    xlab("") +    ylab("Mean proportion (%)") +    theme(panel.background = element_rect(fill = 'transparent'),          panel.grid = element_blank(),          axis.ticks.length = unit(0.2,"lines"),           axis.ticks = element_line(color='black'),          axis.line = element_line(colour = "black"),          axis.title.x=element_text(colour='black', size=9,face = "bold"),          axis.text=element_text(colour='black',size=7,face = "bold"),          legend.title=element_blank(),          legend.text=element_text(size=12,face = "bold",colour = "black",                                   margin = margin(r = 20)),          #legend.position = c(-0.1,0.1),          legend.position = "top",          legend.direction = "horizontal",          legend.key.width = unit(0.8,"cm"),          legend.key.height = unit(0.5,"cm"))#p1for (i in 1:(nrow(diff.mean) - 1))     p1 <- p1 + annotate('rect', xmin = i+0.5, xmax = i+1.5, ymin = -Inf, ymax = Inf,                         fill = ifelse(i %% 2 == 0, 'white', 'gray95'))p1 <- p1 +     geom_bar(stat = "identity",position = "dodge",width = 0.7,colour = "black",linewidth = 0.1) +    scale_fill_manual(values=cbbPalette)#p1## 右侧散点图(Scatter plot on the right)diff.mean$var <- factor(diff.mean$var,levels = levels(abun.bar$variable))diff.mean$p.value <- signif(diff.mean$p.value,3)diff.mean$p.value <- as.character(diff.mean$p.value)p2 <- ggplot(diff.mean,aes(var,estimate,fill = Group)) +    theme(panel.background = element_rect(fill = 'transparent'),          panel.grid = element_blank(),          axis.ticks.length = unit(0.4,"lines"),           axis.ticks = element_line(color='black'),          axis.line = element_line(colour = "black"),          axis.title.x=element_text(colour='black', size=6,face = "bold"),          axis.text=element_text(colour='black',size=7,face = "bold"),          axis.text.y = element_blank(),          legend.position = "none",          axis.line.y = element_blank(),          axis.ticks.y = element_blank(),          plot.title = element_text(size = 10,face = "bold",colour = "black",hjust = 0.5)) +    scale_x_discrete(limits = levels(diff.mean$var)) +    coord_flip() +    xlab("") +    ylab("Difference in mean proportions (%)") +    labs(title="95% confidence intervals") for (i in 1:(nrow(diff.mean) - 1))     p2 <- p2 + annotate('rect', xmin = i+0.5, xmax = i+1.5, ymin = -Inf, ymax = Inf,                         fill = ifelse(i %% 2 == 0, 'white', 'gray95'))p2 <- p2 +    geom_errorbar(aes(ymin = conf.low, ymax = conf.high),                   position = position_dodge(0.8), width = 0.3, size = 0.40) +    geom_point(shape = 21,size = 1.5) +    scale_fill_manual(values=cbbPalette) +    geom_hline(aes(yintercept = 0), linetype = 'dashed', color = 'black')#p2    p3 <- ggplot(diff.mean,aes(var,estimate,fill = Group)) +    geom_text(aes(y = 0,x = var),label = diff.mean$p.value,              hjust = 0,fontface = "bold",inherit.aes = FALSE,size = 3) +    geom_text(aes(x = nrow(diff.mean)/2 +0.5,y = 0.85),label = "P-value (corrected)",             srt = 90,fontface = "bold",size = 2.5) +    coord_flip() +    ylim(c(0,1)) +    theme(panel.background = element_blank(),          panel.grid = element_blank(),          axis.line = element_blank(),          axis.ticks = element_blank(),          axis.text = element_blank(),          axis.title = element_blank())## 图像拼接(Combined plot)p <- p1 + p2 + p3 + plot_layout(widths = c(2,2,1))#p## 保存图像# Save plotggsave("results/Pathway_stamp_abun_filter_0.05.pdf",p,width = 8,height = 6)