rm(list = ls())
setwd("C:\\Users\\Administrator\\Desktop\\microtable") #设置工作目录
library(microeco)
library(magrittr)
library(tidyverse)
library(aplot)
library(ggsci)
feature_table <- read.table('feature_table_g.txt', header = TRUE, row.names = 1, sep = "\t") #特征表
sample_table <- read.table('sample_table_g.txt', header = TRUE, row.names = 1, sep = "\t") #样本表
tax_table <- read.table('tax_table_g.txt', header = TRUE, row.names = 1, sep = "\t", fill = TRUE) #分类表
dataset <- microtable$new(sample_table = sample_table,
otu_table = feature_table,
tax_table = tax_table)
dataset$tidy_dataset()
dataset$rarefy_samples(sample.size = 1000000)
dataset$cal_betadiv(unifrac = FALSE) #beta 多样性
t1 <- trans_beta$new(dataset = dataset, group = "Subgroup", measure = "bray")
bray_curtis <- dataset$beta_diversity$bray
t1$cal_ordination(ordination = "PCoA")
t1$plot_ordination(plot_color = "Subgroup", plot_shape = "Subgroup", plot_type = c("point", "ellipse"))
#################################
# 执行RMANOVA分析
t1$cal_manova(manova_all = FALSE)
# 获取分析结果
manova_results <- t1$res_manova
# 将结果保存到TXT文件
write.table(manova_results, file = "RMANOVA_results.txt", sep = "\t", quote = FALSE)
#################################
t1$cal_manova(manova_all = TRUE)
manova_results <- t1$res_manova
write.table(manova_results, file = "RMANOVA_results_all.txt", sep = "\t", quote = FALSE)
#################################
t1$cal_manova(manova_all = FALSE, group = "Type", by_group = "Group")
manova_results <- t1$res_manova
write.table(manova_results, file = "RMANOVA_results_Group.txt", sep = "\t", quote = FALSE)
################################
t1$cal_manova(manova_all = FALSE, group = "Group", by_group = "Type")
manova_results <- t1$res_manova
write.table(manova_results, file = "RMANOVA_results_Gene.txt", sep = "\t", quote = FALSE)
#################################
library(vegan)
library(tidyverse)
library(ggalt)
library(car)
library(ggforce)
library(ggpubr)
library(patchwork)
pcoa <- cmdscale(bray_curtis, k = 2, eig = TRUE)
points <- as.data.frame(pcoa$points) %>% dplyr::rename(x = "V1", y = "V2")
eig <- pcoa$eig
rownames(points) <- rownames(sample_table)
points$ID <- rownames(points)
points$Subgroup <- sample_table[match(rownames(points), rownames(sample_table)), "Subgroup"]
n <- 0.85
colors <- c("B73_DAS28"="#8FC9E2","B73_DAS42"="#8FC9E2","B73_DAS56"="#8FC9E2","B73_DAS70"="#8FC9E2","Mo17_DAS28"="#ECC97F","Mo17_DAS42"="#ECC97F","Mo17_DAS56"="#ECC97F","Mo17_DAS70"="#ECC97F")
# 定义了颜色和形状的映射关系，用于后续的可视化。
shapes <- c("B73_DAS28"=24, "B73_DAS42"=22, "B73_DAS56"=21, "B73_DAS70"=23,
"Mo17_DAS28"=24, "Mo17_DAS42"=22, "Mo17_DAS56"=21, "Mo17_DAS70"=23)
levels_order <- c("B73_DAS28", "B73_DAS42", "B73_DAS56", "B73_DAS70", "Mo17_DAS28", "Mo17_DAS42", "Mo17_DAS56", "Mo17_DAS70") #定义顺序
points$Subgroup <- factor(points$Subgroup, levels = levels_order)
# 在ggplot中使用这些形状
p1 <- ggplot(points, aes(x = x, y = y, fill = Subgroup, shape = Subgroup)) + #创建了一个ggplot对象，设置了数据框points作为数据，设置了x和y作为x轴和y轴的变量，设置了Group作为填充色和形状的变量。
geom_point(alpha = .7, size = 5) + #添加了散点图层。alpha = .7设置了点的透明度，size = 5设置了点的大小。
scale_shape_manual(values = shapes) + # 设置了形状和填充色的映射关系。使用了之前定义的shapes和colors变量。
scale_fill_manual(values = colors) +
labs(x = paste("PCo1 (", format(100 * eig[1] / sum(eig), digits = 4), "%)", sep = ""),
y = paste("PCo2 (", format(100 * eig[2] / sum(eig), digits = 4), "%)", sep = "")) +
#设置了x轴、y轴和标题的标签。使用了paste函数将多个字符串连接在一起，使用了format函数将数字格式化为指定的小数位数。
geom_mark_ellipse(aes(fill = Subgroup), 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 = 25),
axis.title.x = element_text(size = 28),
axis.title.y = element_text(size = 28),
plot.title = element_text(size = 16, face = "bold"),
legend.key.size = unit(1, "cm"),
legend.title = element_blank(),
legend.text = element_text(size = 14),
panel.border = element_blank(), # 去掉整个图的边框
axis.ticks = element_line(size = 1.5),
axis.ticks.length = unit(0.3, "cm"),
legend.position = "top",
legend.direction = "horizontal",
axis.line = element_line(color = "black", size = 1)) # 添加此行以显示坐标轴线
p1 + theme(text = element_text(family = "TNR"),
axis.text = element_text(family = "TNR"), # 坐标轴刻度文字
axis.title = element_text(family = "TNR"), # 坐标轴标题
plot.title = element_text(family = "TNR", face = "bold"), # 图表标题
legend.text = element_text(family = "TNR"), # 图例文字
legend.title = element_text(family = "TNR")) # 图例标题
p1
# 使用ggsave保存PCoA图为PNG格式
ggsave(filename = "PCoA_plot.png", plot = p1, width = 10, height = 10, units = "in", dpi = 1200)
