主编寄语
微生物组发展到现在，相关分析技术和体系已经十分成熟，而微生物组数据的分析方法更是各式各样，我们借助“R语言在微生物组的最佳实践”这篇之前我们发表的论文带领大家从简单到困难，从基础到流程，从单一功能到组合功能，从上游分析到下游分析的全部代码。全套流程全部借助于Rstudio作为平台，在win和linux和mac下都可以完善运行，大家可以自己准备相关硬件资源，跟着我们这批教程进行全面学习。
这篇综述是我们小组在2023年发表在protein cell 期刊，整体内容几乎囊括了微生物组数据分析的全部内容，代码脚本有2w+行，全部属于公开分享阶段，但由于内容很多，为了方便大家理解和使用学习，这里选择分期带大家运行。
微生信生物分享的这批内容经过如下步骤：
同学A进行代码整理运行，基础注释书写，形成分享稿件的基本样貌；
同学B对代码复现一遍，解决其中无法顺利运行的地方，以及再次注释内容无法让自己明白的地方，完善分享稿件；
同学C再次运行一遍，进行全部流程代码学习，继续注释自己不明白的地方，最终形成完善分享稿件。
我们希望通过这个教程带领需要学习的小伙伴开启全民免费学习微生物组联合挖掘时代。给大家尽量做高质量的教程；
写在前面
“R语言在微生物组的最佳实践”的综述流程，已经带大家学习了4节内容，相信大家已经比较熟悉每节的内容，现在我们继续深入学习微生物组数据挖掘
今天带大家学习一下：差异分析（edger，desep2，stamp）
需要注意的是：每次分析之前，需要加载”数据分析准备“ 这一节的代码，才可以开始下游分析
这次的分享内容希望能够帮助大家更加深入的挖掘自己的数据，给大家带来分析数据的便捷，同时也欢迎大家可以引用我们的文章。
下面是文章的具体信息：
Wen, T., et al., The best practice for microbiome analysis using R. Protein & Cell, 2023. 14(10): p. 713-725.
文章链接：https://doi.org/10.1093/procel/pwad024
流程实践数据：https://github.com/taowenmicro/EasyMicrobiomeR
数据代码获取地址
数据地址：https://github.com/taowenmicro/EasyMicrobiomeR/tree/master/data
全套数据代码整套流程下载地址：https://github.com/taowenmicro/EasyMicrobiomeR/tree/master
差异分析edgR、desep2和stamp差异分析它们的共性和差异
DESeq2、edgeR和stemp都是用于生物信息学中的基因表达数据分析的软件包，它们在差异表达分析方面有共性也有差异。以下是它们的共性和差异的概述：
共性
目的：它们都旨在识别在不同条件或处理下表达量显著变化的基因或功能单元（如在STAMP中可能是代谢途径或功能类别）。
适用数据类型：这些工具都可以应用于高通量测序数据，特别是RNA-seq数据（DESeq2和edgeR）以及宏基因组数据（STAMP）。
统计方法：它们都使用统计模型来评估基因表达量的变化，并通过P值和校正方法（如FDR）来控制多重比较中的错误率。
数据可视化：这些软件包都提供了数据可视化的工具，如火山图、热图等，以帮助用户理解分析结果。
差异
专注领域：DESeq2 和 edgeR 主要关注于单细胞或单个生物体的基因表达分析。
STAMP 专注于宏基因组数据，即多个生物体或环境样本中的功能基因分析。
分析方法：DESeq2 使用负二项分布模型和线性模型来分析数据，并提供了数据标准化的方法。edgeR 也使用负二项分布模型，但提供了不同的数据标准化和变异性估计方法。STAMP 使用的方法包括t检验、ANOVA、DESeq2-like方法等，适用于宏基因组数据的功能丰度比较。
应用范围：DESeq2 和 edgeR 主要用于生物医学研究，如疾病机理、药物作用等。STAMP 用于环境微生物学、生态学和宏基因组学研究，如环境适应性、生态系统功能等。
输出结果：DESeq2 和 edgeR 提供基因级别的差异表达结果。STAMP 提供功能类别或代谢途径级别的差异表达结果。
软件环境：DESeq2 和 edgeR 都是R语言的包，需要在R环境中运行。STAMP 是一个独立的软件，可以在命令行环境中运行，也可以通过图形用户界面操作。
总的来说，DESeq2和edgeR更侧重于单个生物体的基因表达差异分析，而STAMP侧重于宏基因组数据中的功能差异分析。尽管它们在分析方法和应用范围上有所不同，但它们都是强大的工具，可以帮助研究人员在各自的领域内识别生物过程中的关键基因或功能。
差异分析edgR、desep2和stemp差异分析它们的特点和功能
edgeR差异分析
edgeR是一种流行的生物信息学软件包，主要用于对RNA测序（RNA-seq）数据进行差异表达基因分析。它是由加拿大不列颠哥伦比亚大学的Michael D. Lawrence等人开发的R语言包。edgeR被广泛用于生物医学研究中，特别是在转录组学和基因表达分析领域。
edgeR的主要特点和功能包括：
负二项分布模型：edgeR使用负二项分布模型来分析计数数据，这是RNA-seq数据中常见的数据类型。该模型可以考虑到RNA-seq数据的离散性和变异性。
对数CPM（Counts Per Million）转换：edgeR提供了一种数据标准化方法，即将原始的读数计数转换为对数CPM值，以便于进行后续的统计分析。
差异表达检验：edgeR提供了一系列的统计方法来确定基因表达量在不同条件或处理组之间的差异是否显著。这些方法包括精确检验和拟合负二项分布模型的一般线性模型（GLM）。
多重比较校正：在进行多个基因的比较时，edgeR可以应用多重比较校正方法，如FDR（False Discovery Rate）控制，以减少假阳性结果的发生。
生物学解释：edgeR的分析结果可以帮助研究人员识别在特定生物学过程、疾病状态或治疗反应中可能发挥作用的基因。
可视化工具：edgeR提供了多种数据可视化工具，如火山图、散点图和热图等，以帮助用户直观地展示和解释分析结果。
DESeq2差异分析
DESeq2是一种用于RNA测序（RNA-seq）数据差异表达分析的生物信息学软件包，它是DESeq的升级版，提供了更精确和强大的统计方法。DESeq2主要在R语言环境中运行，广泛应用于转录组学研究中，特别是在比较不同条件下基因表达水平的变化时。
DESeq2的主要特点和功能包括：
模型基础的分析：DESeq2使用负二项分布和线性模型来分析RNA-seq数据，考虑到数据的离散性、生物学重复和实验设计中的各种因素。
数据标准化：DESeq2在分析之前对数据进行标准化处理，以消除测序深度和技术偏差的影响，使得不同样本之间的比较更加准确。
差异表达检验：DESeq2提供了统计方法来确定基因表达量在不同条件或处理组之间的变化是否具有统计学意义。它通过拟合线性模型并估计参数的方差来实现这一点。
多重比较校正：DESeq2使用多重比较校正方法，如Benjamini-Hochberg程序，来控制假阳性率（FDR），确保结果的可靠性。
生物学解释：通过分析结果，研究人员可以识别在特定生物学过程、疾病状态或治疗反应中可能发挥作用的基因，为后续的实验验证和功能研究提供依据。
可视化工具：DESeq2可以生成多种可视化图表，如火山图、箱线图和MA图等，帮助用户直观地展示和解释分析结果。
stamp差异分析
STAMP（Statistical Analysis of Metagenomic Profiles）是一种生物信息学软件，专门设计用于对宏基因组数据进行统计分析。它提供了一系列工具和方法，使研究人员能够比较不同宏基因组样本之间的功能基因组成，从而识别在特定环境或实验条件下显著变化的微生物功能。
STAMP的主要特点和功能包括：
功能丰度比较：STAMP允许用户比较不同样本中功能基因的相对丰度，以发现与特定条件或处理相关的微生物功能的变化。
多种统计方法：STAMP支持多种统计检验方法，包括t检验、ANOVA、非参数检验等，以及专门为宏基因组数据设计的DESeq2-like方法，这些方法可以帮助评估功能基因丰度差异的统计显著性。
多重比较校正：在进行多重比较时，STAMP提供校正方法（如FDR，False Discovery Rate）来控制错误发现率，确保结果的可靠性。
数据转换和标准化：STAMP可以对输入的数据进行转换和标准化处理，以消除不同样本之间的技术差异，使得比较结果更加准确和可比。
结果可视化：STAMP提供了多种数据可视化选项，如火山图、热图等，帮助用户直观地展示和解释分析结果。
功能注释：STAMP可以结合KEGG（Kyoto Encyclopedia of Genes and Genomes）等数据库，对检测到的功能基因进行注释，从而提供生物学意义上的解释。
下面开始实践操作
edgeR -Volcano
edgeR的火山图（Volcano plot）是一种用于可视化基因表达差异的工具，它适用于多种场景，尤其是在比较不同条件下基因表达变化的研究中。以下是edgeR火山图的一些适用场景：
条件比较：当你想要比较两种或多种不同生物学条件（如健康与疾病、不同时间点、不同处理等）下的基因表达差异时，火山图可以帮助你直观地展示这些差异。
筛选显著变化基因：火山图可以清晰地展示出哪些基因的表达量发生了显著变化。通过设置特定的对数变化倍数（logFC）和P值阈值，你可以快速识别出感兴趣的基因。
发现生物标志物：在疾病研究中，火山图可以帮助研究人员发现潜在的生物标志物或治疗靶点，这些基因在疾病发生或进展中可能起到关键作用。
数据探索：火山图是一种有用的数据探索工具，它可以帮助你在大量数据中发现模式和趋势，为进一步的数据分析提供线索。
方法比较：如果你在研究中使用了不同的统计方法或数据分析流程，火山图可以帮助你比较这些方法的效果，了解它们对结果的影响。
edgeR火山图是一种多功能的工具，适用于各种需要比较和展示基因表达差异的场景。通过火山图，研究人员可以更有效地分析和解释RNA-seq数据，推动生物学研究的进展。
# 数据载入ps = readRDS("./data/dataNEW/ps_16s.rds")diffpath = paste(otupath,"/diff_tax/",sep = "")dir.create(diffpath)diffpath.1 = paste(diffpath,"/DEsep2/",sep = "")dir.create(diffpath.1)diffpath.2 = diffpath.1# 准备脚本group = "Group"pvalue = 0.05lfc =0artGroup = NULLmethod = "TMM"j = 2 # 设置分析的层次为门path = diffpathb = NULL# 读取界门纲目if (j %in% c("OTU","gene","meta")) {ps = ps} else if (j %in% c(1:7)) {ps = ps %>%ggClusterNet::tax_glom_wt(ranks = j)} else if (j %in% c("Kingdom","Phylum","Class","Order","Family","Genus","Species")){} else {ps = psprint("unknown j, checked please")}sub_design <- as.data.frame(phyloseq::sample_data(ps))# 设置sample的分组Desep_group <- as.character(levels(as.factor(sub_design$Group)))Desep_groupif ( is.null(artGroup)) {#--构造两两组合#-----aaa = combn(Desep_group,2)# sub_design <- as.data.frame(sample_data(ps))}if (!is.null(artGroup)) {aaa = as.matrix(b)}otu_table = as.data.frame(ggClusterNet::vegan_otu(ps))count = as.matrix(otu_table)count <- t(count)sub_design <- as.data.frame(phyloseq::sample_data(ps))dim(sub_design)sub_design$SampleType = as.character(sub_design$Group)sub_design$SampleType <- as.factor(sub_design$Group)# create DGE list 进行标准化# DGE表有两个部分，一个是物种在样本内出现的次数，一个是每个样品的文库大小（测序深度）d = edgeR::DGEList(counts=count, group=sub_design$SampleType)d$samplesd = edgeR::calcNormFactors(d,method=method)#默认为TMM标准化# Building experiment matrix建筑实验矩阵design.mat = model.matrix(~ 0 + d$samples$group)colnames(design.mat)=levels(sub_design$SampleType)d2 = edgeR::estimateGLMCommonDisp(d, design.mat)d2 = edgeR::estimateGLMTagwiseDisp(d2, design.mat)fit = edgeR::glmFit(d2, design.mat)#------------根据分组提取需要的差异结果#------------# 开始算groupX与gourpY的差异了# 以火山图的形式出图，并保存至path里for (i in 1:dim(aaa)[2]) {# i = 1Desep_group = aaa[,i]print( Desep_group)# head(design)# 设置比较组写在前面的分组为enrich表明第一个分组含量高# ?limma::makeContrastsgroup = paste(Desep_group[1],Desep_group[2],sep = "-")groupBvsA <- limma::makeContrasts(contrasts = group,levels=c( as.character(levels(as.factor(sub_design$Group)))) )#注意是以GF1为对照做的比较# 组间比较,统计Fold change, Pvaluelrt = edgeR::glmLRT(fit,contrast=BvsA)# FDR检验，控制假阳性率小于5%de_lrt = edgeR::decideTestsDGE(lrt, adjust.method="fdr", p.value=pvalue,lfc=lfc)#lfc=0这个是默认值summary(de_lrt)# 导出计算结果x=lrt$tablex$sig=de_lrthead(x)#------差异结果符合otu表格的顺序row.names(count)[1:6]x <- cbind(x, padj = p.adjust(x$PValue, method = "fdr"))enriched = row.names(subset(x,sig==1))depleted = row.names(subset(x,sig==-1))x$level = as.factor(ifelse(as.vector(x$sig) ==1, "enriched",ifelse(as.vector(x$sig)==-1, "depleted","nosig")))x = data.frame(row.names = row.names(x),logFC = x$logFC,level = x$level,p = x$PValue)head(x)# colnames(x) = paste(group,colnames(x),sep = "")# x = res# head(x)#------差异结果符合otu表格的顺序# x = data.frame(row.names = row.names(x),logFC = x$log2FoldChange,level = x$level,p = x$pvalue)x1 = x %>%dplyr::filter(level %in% c("enriched","depleted","nosig") )head(x1)x1$Genus = row.names(x1)# x$level = factor(x$level,levels = c("enriched","depleted","nosig"))if (nrow(x1)<= 1) {}x2 <- x1 %>%dplyr::mutate(ord = logFC^2) %>%dplyr::filter(level != "nosig") %>%dplyr::arrange(desc(ord)) %>%head(n = 5)file = paste(path,"/",group,j,"_","Edger_Volcano_Top5.csv",sep = "")write.csv(x2,file,quote = F)head(x2)# 想调颜色在这里p <- ggplot(x1,aes(x =logFC ,y = -log2(p), colour=level)) +geom_point() +geom_hline(yintercept=-log10(0.2),linetype=4,color = 'black',size = 0.5) +geom_vline(xintercept=c(-1,1),linetype=3,color = 'black',size = 0.5) +ggrepel::geom_text_repel(data=x2, aes(x =logFC ,y = -log2(p), label=Genus), size=1) +scale_color_manual(values = c('blue2','red2', 'gray30')) +ggtitle(group) + theme_bw()pfile = paste(path,"/",group,j,"_","Edger_Volcano.pdf",sep = "")ggsave(file,p,width = 8,height = 6)file = paste(path,"/",group,j,"_","Edger_Volcano.png",sep = "")ggsave(file,p,width = 8,height = 6)colnames(x) = paste(group,colnames(x),sep = "")if (i ==1) {table =x}if (i != 1) {table = cbind(table,x)}}x = table
火山图是一种直观的工具，用于识别和比较两个条件下显著变化的基因或物种，从而帮助研究人员识别可能在生物学过程中起重要作用的候选基因或物种。在这种情况下，火山图可能用于微生物组学研究，以识别在不同环境或条件下丰度发生变化的微生物群体。
图表标题为 "Group1-Group2"，表示该火山图比较的是两个组别（Group1 和 Group2）之间的数据差异。
水平轴（x轴）表示对数变化倍数（logFC），即两个组别之间相对表达量的比值的对数值。正值表示在Group1中的表达量高于Group2，负值则相反。
垂直轴（y轴）表示显著性水平（level），通常与P值相关，用于衡量表达量变化的统计显著性。"nosig" 表示变化不显著，"depleted" 表示在Group1中的表达量显著低于Group2，而 "enriched" 表示在Group1中的表达量显著高于Group2。
图中的线条和点表示不同的数据点，根据其对数变化倍数和显著性水平进行定位。线条通常用来划分显著性和变化幅度的阈值。
###########添加物种丰度#----------# dim(count)# str(count)count = as.matrix(count)norm = t(t(count)/colSums(count)) #* 100 # normalization to total 100dim(norm)norm1 = norm %>%t() %>% as.data.frame()# head(norm1)#数据分组计算平均值library("tidyverse")head(norm1)iris.split <- split(norm1,as.factor(sub_design$SampleType))iris.apply <- lapply(iris.split,function(x)colMeans(x))# 组合结果iris.combine <- do.call(rbind,iris.apply)norm2= t(iris.combine)#head(norm)str(norm2)norm2 = as.data.frame(norm2)# dim(x)head(norm2)# 得到物种在每个group里的相对数量x = cbind(x,norm2)head(x)#在加入这个文件taxonomy时，去除后面两列不相干的列# 读取taxonomy，并添加各列名称if (!is.null(ps@tax_table)) {taxonomy = as.data.frame(ggClusterNet::vegan_tax(ps))head(taxonomy)# taxonomy <- as.data.frame(tax_table(ps1))#发现这个注释文件并不适用于直接作图。#采用excel将其分列处理，并且删去最后一列，才可以运行if (length(colnames(taxonomy)) == 6) {colnames(taxonomy) = c("kingdom","phylum","class","order","family","genus")}else if (length(colnames(taxonomy)) == 7) {colnames(taxonomy) = c("kingdom","phylum","class","order","family","genus","species")}else if (length(colnames(taxonomy)) == 8) {colnames(taxonomy) = c("kingdom","phylum","class","order","family","genus","species","rep")}# colnames(taxonomy) = c("kingdom","phylum","class","order","family","genus")# Taxonomy排序，并筛选OTU表中存在的library(dplyr)taxonomy$id=rownames(taxonomy)# head(taxonomy)tax = taxonomy[row.names(x),]x = x[rownames(tax), ] # reorder according to taxif (length(colnames(taxonomy)) == 7) {x = x[rownames(tax), ] # reorder according to taxx$phylum = gsub("","",tax$phylum,perl=TRUE)x$class = gsub("","",tax$class,perl=TRUE)x$order = gsub("","",tax$order,perl=TRUE)x$family = gsub("","",tax$family,perl=TRUE)x$genus = gsub("","",tax$genus,perl=TRUE)# x$species = gsub("","",tax$species,perl=TRUE)}else if (length(colnames(taxonomy)) == 8) {x$phylum = gsub("","",tax$phylum,perl=TRUE)x$class = gsub("","",tax$class,perl=TRUE)x$order = gsub("","",tax$order,perl=TRUE)x$family = gsub("","",tax$family,perl=TRUE)x$genus = gsub("","",tax$genus,perl=TRUE)x$species = gsub("","",tax$species,perl=TRUE)}else if (length(colnames(taxonomy)) == 9) {x = x[rownames(tax), ] # reorder according to taxx$phylum = gsub("","",tax$phylum,perl=TRUE)x$class = gsub("","",tax$class,perl=TRUE)x$order = gsub("","",tax$order,perl=TRUE)x$family = gsub("","",tax$family,perl=TRUE)x$genus = gsub("","",tax$genus,perl=TRUE)x$species = gsub("","",tax$species,perl=TRUE)}} else {x = cbind(x,tax)}res = xhead(res)filename = paste(diffpath.2,"/","_",j,"_","edger_all.csv",sep = "")write.csv(res,filename)
edger_all.csv：该文件包含了不同微生物群体在三个不同组（Group1, Group2, Group3）之间的相对丰度比较结果。每个微生物群体的比较结果包括对数变化倍数（logFC）、显著性水平（level）和P值（p）。这些数据通常用于分析微生物组数据，以了解不同环境或条件下微生物群体的相对丰度变化。
对数变化倍数（logFC）：表示Group1与Group2、Group1与Group3、Group2与Group3之间的相对丰度变化。
显著性水平（level）：指示变化是否具有统计学意义，"nosig"表示不显著，"depleted"表示在Group1相对于Group2或Group3中丰度降低，"enriched"表示丰度增加。
P值（p）：用于评估差异的统计显著性，P值越小，表示差异越显著。
DESep2-Volcano
DESeq2-Volcano图是一种基于DESeq2软件包生成的数据可视化工具，用于展示RNA-seq数据分析中基因表达量的差异。与edgeR火山图类似，DESeq2-Volcano图也广泛应用于多种生物学研究场景，特别是在差异表达基因分析中。以下是DESeq2-Volcano图的一些适用场景：
差异表达分析：在比较不同生物学样本（如不同组织、不同疾病状态、不同实验处理等）的基因表达水平时，DESeq2-Volcano图可以直观地展示基因表达量的变化。
统计显著性展示：通过火山图的水平轴（通常表示对数变化倍数logFC）和垂直轴（表示统计显著性，如负对数P值）的分布，研究人员可以快速识别出在统计上显著差异表达的基因。
基因筛选：研究人员可以根据设定的显著性阈值（如调整后的P值FDR）和表达量变化的阈值，从成千上万的基因中筛选出潜在的候选基因进行进一步研究。
生物学假设验证：DESeq2-Volcano图可以帮助验证特定的生物学假设，例如，某个基因是否在特定条件下上调或下调。
DESeq2-Volcano图因其直观性和信息丰富性，在RNA-seq数据分析中被广泛使用。它不仅有助于研究人员识别和解释基因表达的变化，还能够促进科学发现和新知识的产生。
基本流程跟edgR相似
# 设置分析的层次为属j = "Genus"group = "Group"pvalue = 0.05artGroup = NULLpath = diffpath# 读取界门纲目# ps = ps %>%# ggClusterNet::tax_glom_wt(ranks = j)if (j %in% c("OTU","gene","meta")) {ps = ps} else if (j %in% c(1:7)) {ps = ps %>%ggClusterNet::tax_glom_wt(ranks = j)} else if (j %in% c("Kingdom","Phylum","Class","Order","Family","Genus","Species")){} else {ps = psprint("unknown j, checked please")}# 设置分组Desep_group <- ps %>%phyloseq::sample_data() %>%.$Group %>%as.factor() %>%levels() %>%as.character()if ( is.null(artGroup)) {#--构造两两组合#-----aaa = combn(Desep_group,2)# sub_design <- as.data.frame(sample_data(ps))} else if (!is.null(artGroup)) {aaa = as.matrix(b )}count <- ps %>%ggClusterNet::vegan_otu() %>% round(0) %>%t()map = ps %>%phyloseq::sample_data() %>%as.tibble() %>%as.data.frame()dds <- DESeq2::DESeqDataSetFromMatrix(countData = count,colData = map,design = ~ Group)dds2 <- DESeq2::DESeq(dds) ##第二步,标准化# resultsNames(dds2)#------------根据分组提取需要的差异结果#------------for (i in 1:dim(aaa)[2]) {# i = 1Desep_group = aaa[,i]print( Desep_group)# head(design)# 设置比较组写在前面的分组为enrich表明第一个分组含量高group = paste(Desep_group[1],Desep_group[2],sep = "-")group# 将结果用results()函数来获取，赋值给res变量res <- DESeq2::results(dds2, contrast=c("Group",Desep_group ),alpha=0.05)# 导出计算结果x = reshead(x)#------差异结果符合otu表格的顺序x$level = as.factor(ifelse(as.vector(x$padj) < 0.05 & x$log2FoldChange > 0, "enriched",ifelse(as.vector(x$padj) < 0.05 &x$log2FoldChange < 0, "depleted","nosig")))x = data.frame(row.names = row.names(x),logFC = x$log2FoldChange,level = x$level,p = x$pvalue)x1 = x %>%filter(level %in% c("enriched","depleted","nosig") )head(x1)x1$Genus = row.names(x1)# x$level = factor(x$level,levels = c("enriched","depleted","nosig"))x2 <- x1 %>%dplyr::mutate(ord = logFC^2) %>%dplyr::filter(level != "nosig") %>%dplyr::arrange(desc(ord)) %>%head(n = 5)file = paste(path,"/",group,j,"_","DESep2_Volcano_Top5.csv",sep = "")write.csv(x2,file,quote = F)head(x2)#出图，相关颜色在这里改p <- ggplot(x1,aes(x =logFC ,y = -log2(p), colour=level)) +geom_point() +geom_hline(yintercept=-log10(0.2),linetype=4,color = 'black',size = 0.5) +geom_vline(xintercept=c(-1,1),linetype=3,color = 'black',size = 0.5) +ggrepel::geom_text_repel(data=x2, aes(x =logFC ,y = -log2(p), label=Genus), size=1) +scale_color_manual(values = c('blue2','red2', 'gray30')) +ggtitle(group) + theme_bw()pfile = paste(path,"/",group,j,"_","DESep2_Volcano.pdf",sep = "")ggsave(file,p,width = 8,height = 6)file = paste(path,"/",group,j,"_","DESep2_Volcano.png",sep = "")ggsave(file,p,width = 8,height = 6)colnames(x) = paste(group,colnames(x),sep = "")if (i ==1) {table =x}if (i != 1) {table = cbind(table,x)}}
图表标题为 "Group1-Group2"，表示该火山图比较的是两个组别（Group1 和 Group2）之间的数据差异。
水平轴（x轴）表示对数变化倍数（logFC），即两个组别之间相对表达量的比值的对数值。正值表示在Group1中的表达量高于Group2，负值则相反。
垂直轴（y轴）表示显著性水平（level），通常与P值相关，用于衡量表达量变化的统计显著性。"nosig" 表示变化不显著，"depleted" 表示在Group1中的表达量显著低于Group2，而 "enriched" 表示在Group1中的表达量显著高于Group2。
图中特别标注了 "ASV_4375"，这可能意味着这个特定的物种或基因在Group1相对于Group2中表现出显著的变化。根据其在图上的位置，可以判断它是被富集（enriched）还是耗竭（depleted）。
count = as.matrix(count)norm = t(t(count)/colSums(count)) #* 100 # normalization to total 100dim(norm)norm1 = norm %>%t() %>% as.data.frame()# head(norm1)#数据分组计算平均值# library("tidyverse")head(norm1)iris.split <- split(norm1,as.factor(map$Group))iris.apply <- lapply(iris.split,function(x)colMeans(x))# 组合结果iris.combine <- do.call(rbind,iris.apply)norm2= t(iris.combine)#head(norm)str(norm2)norm2 = as.data.frame(norm2)head(norm2)x = cbind(table,norm2)head(x)#在加入这个文件taxonomy时，去除后面两列不相干的列# 读取taxonomy，并添加各列名称if (!is.null(ps@tax_table)) {taxonomy = as.data.frame(ggClusterNet::vegan_tax(ps))head(taxonomy)# taxonomy <- as.data.frame(tax_table(ps1))#发现这个注释文件并不适用于直接作图。#采用excel将其分列处理，并且删去最后一列，才可以运行if (length(colnames(taxonomy)) == 6) {colnames(taxonomy) = c("kingdom","phylum","class","order","family","genus")}else if (length(colnames(taxonomy)) == 7) {colnames(taxonomy) = c("kingdom","phylum","class","order","family","genus","species")}else if (length(colnames(taxonomy)) == 8) {colnames(taxonomy) = c("kingdom","phylum","class","order","family","genus","species","rep")}# colnames(taxonomy) = c("kingdom","phylum","class","order","family","genus")# Taxonomy排序，并筛选OTU表中存在的library(dplyr)taxonomy$id=rownames(taxonomy)# head(taxonomy)tax = taxonomy[row.names(x),]x = x[rownames(tax), ] # reorder according to taxif (length(colnames(taxonomy)) == 7) {x = x[rownames(tax), ] # reorder according to taxx$phylum = gsub("","",tax$phylum,perl=TRUE)x$class = gsub("","",tax$class,perl=TRUE)x$order = gsub("","",tax$order,perl=TRUE)x$family = gsub("","",tax$family,perl=TRUE)x$genus = gsub("","",tax$genus,perl=TRUE)# x$species = gsub("","",tax$species,perl=TRUE)}else if (length(colnames(taxonomy)) == 8) {x$phylum = gsub("","",tax$phylum,perl=TRUE)x$class = gsub("","",tax$class,perl=TRUE)x$order = gsub("","",tax$order,perl=TRUE)x$family = gsub("","",tax$family,perl=TRUE)x$genus = gsub("","",tax$genus,perl=TRUE)x$species = gsub("","",tax$species,perl=TRUE)}else if (length(colnames(taxonomy)) == 9) {x = x[rownames(tax), ] # reorder according to taxx$phylum = gsub("","",tax$phylum,perl=TRUE)x$class = gsub("","",tax$class,perl=TRUE)x$order = gsub("","",tax$order,perl=TRUE)x$family = gsub("","",tax$family,perl=TRUE)x$genus = gsub("","",tax$genus,perl=TRUE)x$species = gsub("","",tax$species,perl=TRUE)}} else {x = x}res = xhead(res)filename = paste(diffpath.1,"/","_",j,"_","DESep2_all.csv",sep = "")write.csv(res,filename,quote = F)
stamp_差异分析
stamp（Statistical Analysis of Metagenomic Profiles）是一种专门用于宏基因组数据差异分析的软件工具。它适用于生物学领域的多种场景，特别是在研究微生物群落结构和功能以及它们在不同环境条件下的变化时。以下是stemp差异分析在生物学领域的一些使用场景：
环境微生物研究：stemp可以用于比较不同环境样本（如土壤、水体、大气等）中的微生物群落组成和功能基因的丰度差异，帮助研究人员理解环境因素如何影响微生物生态系统。
宿主-微生物相互作用：在研究宿主（如人类、动物、植物）与微生物之间的相互作用时，stemp可以用来分析宿主相关微生物群落的变化，以及这些变化如何影响宿主的健康和疾病状态。
疾病与健康比较：stemp可以用来研究疾病状态下与健康状态下微生物群落的差异，识别与疾病相关的微生物标记物或潜在的致病微生物。
药物和治疗效应：在评估药物或其他治疗方法对微生物群落的影响时，stamp可以用来分析治疗前后微生物群落结构和功能的变化，从而揭示治疗的微生物学基础。
生态学和进化研究：stemp可以用于生态学研究，比较不同生态系统中微生物群落的差异，以及这些差异如何反映生态系统的功能和稳定性。
农业和土壤科学：在农业和土壤科学研究中，stemp可以用来分析不同农业实践、施肥或土壤管理措施对土壤微生物群落的影响。
stamp的差异分析功能强大，可以处理大量的宏基因组数据，并提供丰富的统计方法和数据可视化工具，使其成为生物学领域中研究微生物群落结构和功能的有力工具。通过stemp，研究人员可以更好地理解微生物在自然界和人类社会中的多样性、功能和生态作用。
ps = readRDS("./data/dataNEW/ps_16s.rds")diffpath = paste(otupath,"/stemp_diff/",sep = "")dir.create(diffpath)allgroup <- combn(unique(map$Group),2)#只做排列组合里的第一种，在这里改allgroup[,1]ps_sub <- phyloseq::subset_samples(ps,Group %in% allgroup[,1]);ps_subTop = 20#显示名称的突出微生物的数量ranks = 6#做门method = "TMM"test.method = "t.test"#--门水平合并data = ggClusterNet::tax_glom_wt(ps_sub,ranks = phyloseq::rank_names(ps)[ranks]) %>%ggClusterNet::scale_micro(method = method) %>%ggClusterNet::filter_OTU_ps(Top = 200) %>%ggClusterNet::vegan_otu() %>%as.data.frame()tem = colnames(data)data$ID = row.names(data)data <- data %>%dplyr::inner_join(as.tibble(phyloseq::sample_data(ps)),by = "ID")data$Group = as.factor(data$Group)# 做t检验，并检查p值if (test.method == "t.test") {diff <- data[,tem] %>%# dplyr::select_if(is.numeric) %>%purrr::map_df(~ broom::tidy(t.test(. ~ Group,data = data)), .id = 'var')} else if(test.method == "wilcox.test"){diff <- data[,tem] %>%# dplyr::select_if(is.numeric) %>%purrr::map_df(~ broom::tidy(wilcox.test(. ~ Group,data = data)), .id = 'var')}diff$p.value[is.nan(diff$p.value)] = 1diff$p.value <- p.adjust(diff$p.value,"bonferroni")tem = diff$p.value [diff$p.value < 0.05] %>% length()if (tem > 30) {diff <- diff %>%dplyr::filter(p.value < 0.05) %>%head(30)} else {diff <- diff %>%# filter(p.value < 0.05) %>%head(30)}# diff <- diff %>% filter(p.value < 0.05)# diff1$p.value <- p.adjust(diff1$p.value,"bonferroni")# diff1 <- diff1 %>% filter(p.value < 0.05)# 开始出图abun.bar <- data[,c(diff$var,"Group")] %>%tidyr::gather(variable,value,-Group) %>%dplyr::group_by(variable,Group) %>%dplyr::summarise(Mean = mean(value))diff.mean <- diff[,c("var","estimate","conf.low","conf.high","p.value")]diff.mean$Group <- c(ifelse(diff.mean$estimate >0,levels(data$Group)[1],levels(data$Group)[2]))diff.mean <- diff.mean[order(diff.mean$estimate,decreasing = TRUE),]cbbPalette <- c("#E69F00", "#56B4E9")# 设置group的颜色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.4,"lines"),axis.ticks = element_line(color='black'),axis.line = element_line(colour = "black"),axis.title.x=element_text(colour='black', size=12,face = "bold"),axis.text=element_text(colour='black',size=10,face = "bold"),legend.title=element_blank(),# legend.text=element_text(size=12,face = "bold",colour = "black",# margin = margin(r = 20)),legend.position = c(-1,-0.1),legend.direction = "horizontal",legend.key.width = unit(0.8,"cm"),legend.key.height = unit(0.5,"cm"))p1#出的裸图，只有横纵坐标for (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 <- p1 +geom_bar(stat = "identity",position = "dodge",width = 0.7,colour = "black") +scale_fill_manual(values=cbbPalette) + theme(legend.position = "bottom")p1#柱状图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=12,face = "bold"),axis.text=element_text(colour='black',size=10,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 = 15,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.5, size = 0.5) +geom_point(shape = 21,size = 3) +scale_fill_manual(values=cbbPalette) +geom_hline(aes(yintercept = 0), linetype = 'dashed', color = 'black')#p2为cleverland点图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 = 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())#p3为数值标记#合并图片library(patchwork)p <- p1 + p2 + p3 + plot_layout(widths = c(4,6,2))pi = 1# filename = paste(diffpath,"/",paste(allgroup[,i][1],allgroup[,i][2],sep = "_"),"stemp_P_plot.csv",sep = "")# write.csv(diff.mean,filename)filename = paste(diffpath,"/",paste(allgroup[,i][1],allgroup[,i][2],sep = "_"),phyloseq::rank.names(ps)[j],"stemp_plot.pdf",sep = "")ggsave(filename,p,width = 14,height = 6)filename = paste(diffpath,"/",paste(allgroup[,i][1],allgroup[,i][2],sep = "_"),phyloseq::rank.names(ps)[j],"stemp_plot.jpg",sep = "")ggsave(filename,p,width = 14,height = 6)detach("package:patchwork")
这张图像是stemp差异分析结果的一部分，它提供了一种直观的方式来展示和比较不同微生物属在两个组别中的相对丰度和差异显著性。通过这种分析，研究人员可以识别在特定环境条件或实验处理下可能发生变化的关键微生物属，这对于理解微生物群落的动态和功能具有重要意义。
微生物属及其相对丰度：图像列出了一系列微生物属，如Spartobacteria_genera_incertae_sedis、Ralstonia、Dongia、Nitrosopumilus、Dokdonella、Chondromyces 等，并显示了它们在Group1和Group2中的相对丰度百分比。
平均比例：图像的x轴表示两个组别中每个微生物属的平均丰度百分比（Mean proportion (%)），这可以帮助理解每个属在不同组别中的相对重要性。
差异比例：图像的y轴表示两个组别中微生物属丰度的平均差异百分比（Difference in mean proportions (%)），正值表示在Group1中的丰度高于Group2，负值则相反。
置信区间：图像中的线条表示差异比例的95%置信区间，提供了关于丰度差异统计显著性的信息。如果置信区间不包含零，则认为差异在统计上是显著的。
显著性标记：图像中的某些点可能用不同的颜色或标记来表示其显著性水平，例如，显著高于（P < 0.05）或显著低于（P < 0.05）的属可能会有不同的标记。
本次运行内容结束，如果觉得有用的话，欢迎大家来引用我们的文章：Wen, T., et al., The best practice for microbiome analysis using R. Protein & Cell, 2023. 14(10): p. 713-725.
根际互作生物学研究室 简介
根际互作生物学研究室是沈其荣院士土壤微生物与有机肥团队下的一个关注于根际互作的研究小组。本小组由袁军教授带领，主要关注：1.植物和微生物互作在抗病过程中的作用；2 环境微生物大数据整合研究；3 环境代谢组及其与微生物过程研究体系开发和应用。团队在过去三年中在 Nature Communications，ISME J，Microbiome，SCLS，New Phytologist，iMeta，Fundamental Research， PCE，SBB，JAFC（封面），Horticulture Research，SEL（封面），BMC plant biology等期刊上发表了多篇文章。欢迎关注 微生信生物 公众号对本研究小组进行了解。
撰写：文涛 牛国庆 文浩仰 周克成
修改：文涛
排版：杨雯儀
审核：袁军 团队工作及其成果 （点击查看）
了解 交流 合作
小组负责人邮箱 袁军：junyuan@njau.edu.cn；
小组成员文涛：taowen@njau.edu.cn等
团队公众号：微生信生物 添加主编微信，或者后台留言。
1.仅限相关专业或研究方向人员添加，必须实名，不实名则默认忽略。
2.非相关专业的其他人员及推广宣传人员禁止添加。
3.添加主编微信需和简单聊一聊专业相关问题，等待主编判断后，可拉群。
4 微生信生物VIP微信群不受限制，给微生信生物供稿一次即可加入（群里发送推文代码+高效协助解决推文运行等问题+日常问题咨询回复）。
加主编微信 加入群聊
目前营销人员过多，为了维护微生信生物3年来维护的超5500人群聊，目前更新进群要求：
