rm(list = ls())
setwd("C:\\Users\\Administrator\\Desktop\\New_microtable") #设置工作目录
library(microeco)
library(magrittr)
library(dplyr)
library(tibble)
feature_table <- read.table('Bac_genus.txt', header = TRUE, row.names = 1, sep = "\t", fill = TRUE) #特征表
# 检查并处理缺失值
if (any(is.na(feature_table))) {
feature_table[is.na(feature_table)] <- 0 # 将缺失值填充为零
}
sample_table <- read.table('sample_table.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$sample_sums() %>% range #计算并查看样本总数的范围
dataset$rarefy_samples(sample.size = 1000000) #执行重采样，标准化样本中的测序深度
dataset$sample_sums() %>% range #计算并查看标准化后样本总数的范围
dataset$cal_abund() #计算每个分类等级的分类群丰度
#class(dataset$taxa_abund)
dataset$taxa_abund$Genus[1:5, 1:5]
t1 <- trans_abund$new(dataset = dataset, taxrank = "Genus", ntaxa = 10, groupmean = "Subgroup")
#t1 <- trans_abund$new(dataset = dataset, taxrank = "Genus", ntaxa = 10, groupmean = "Group")
#t1 <- trans_abund$new(dataset = dataset, taxrank = "Genus", ntaxa = 10)
write.table(t1$data_abund, file = "trans_abund_data_DAS.txt", sep = "\t", row.names = TRUE, col.names = NA, quote = FALSE)
########################################################################################
# 清空当前环境中的所有对象
rm (list = ls ())
setwd("C:\\Users\\Administrator\\Desktop\\New_microtable") #设置工作目录
library(ggplot2) # 用于绘图
library(ggalluvial) # 用于绘制柱状图背后的条带
library(grid) # 用于访问unit函数
ra <- as.matrix(read.table("Genus_top15.txt", row.names =1,
header = F, sep = "\t")) # 读入相对丰度数据并转换为矩阵方便后续数据整理
group <- c("B73", "Mo17") # 品种变量
code <- c("B73_DAS28","B73_DAS42","B73_DAS56","B73_DAS70","Mo17_DAS28","Mo17_DAS42","Mo17_DAS56","Mo17_DAS70") # 8组处理变量
code_rep <- rep(code, each = 15)
taxa_rep <- rep(rownames(ra), 8)
cultivar_rep <- rep(group, each = 60)
abundance_vec <- as.vector(ra)
dat <- data.frame(code = code_rep,
Phylum = taxa_rep,
cultivar = cultivar_rep,
abundance = abundance_vec)
dat$Phylum <- factor(dat$Phylum, levels = c("Streptomyces","Sphingobium","Sphingomonas","Bradyrhizobium",
"Variovorax","Lysobacter","Pseudomonas","Acidovorax","Nocardioides",
"Burkholderia","Stenotrophomonas","Mesorhizobium","Microbacterium","Rhizobium","Roseateles"))
ggplot(dat, aes(x = code, y = abundance, fill = Phylum))+
geom_bar(stat = "identity", width = 0.6, alpha = 1)+
geom_flow(aes(alluvium = Phylum), alpha = 0.4) +
scale_fill_manual(values = c(
"#F2A2C7", "#E9C46A", "#88C695", "#F4E3A1", "#76B4BD", "#C9A8C5", "#6D6DB5", "#5BA8A0",
"#9CAD60", "#C0C0C0", "#BC8F8F", "#76B4BD", "#6AB47B", "#89CFF0", "#8A2BE2"
))+
theme_bw()+
facet_grid(.~cultivar, scales = "free_x", space = "free_x")+
xlab("")+
ylab("Relative abundance (%)")+
theme(panel.grid.major.x = element_blank(),
axis.text.x = element_text(size = rel(1.3), angle = 0, face = "plain", color = "black", family = "Times New Roman"),
axis.text.y = element_text(size=rel(1.8), face = "plain", color = "black", family = "Times New Roman"),
legend.text = element_text(size = rel(1), family = "Times New Roman"),
strip.text.x = element_text(size = rel(2), face = "bold", family = "Times New Roman"),
strip.background = element_rect(fill = "lightblue", size = 1),
panel.border = element_rect(colour = "black", fill = NA, size = 1.2),
axis.title.y = element_text(size = 18, family = "Times New Roman"),
legend.position = "bottom",
legend.title = element_blank(),
legend.key.size = unit(0.6, "cm"),
legend.spacing.x = unit(0.1, "cm"),
legend.box.margin = margin(t = -20, unit = "pt")) +
guides(fill = guide_legend(ncol = 7)) +
scale_y_continuous(expand = c(0, 0.01))
ggsave("Abundance2.png", width = 10, height = 7, dpi = 1200) # 图片导出，导出为pdf文件，设置图片长和宽标签：abund,microtable,library,dataset,堆叠,柱状图,丰度,计算 From： https://www.cnblogs.com/wzbzk/p/18234372