indicator简介

指示物种分析（Indicator Species Analysis）是扩增子数据中使用非常广泛的一种经典分析方式，可在许多高分文章主图中见到。通过indicator value 寻找处理的指示微生物(类似于Lefse分析寻找Biomarker)，即对处理存在显著正响应的微生物，从而进行处理效应的研究。同时可以基于indicator值来比较微生物和处理关联性大小。一下列举indicator分析的例子和代码，以及可视化过程，以供大家学习！！！

Fig. 3 Bipartite networks displaycropping system specific OTUs in the soil and root bacterial and fungal communities as determined using indicator species analysis.(Hartmann et al., 2018, Microbiome)

Table 4 Indicator species analysis for control soils at each site along the ambient N deposition gradient.

(Jessica et al., 2020， GCB)

出图结果

分析流程

########### 寻找处理的indicator和bipartite网络可视化

library(ggplot2)

library(vegan)

library(indicspecies)

library(edgeR)

#加载函数（可后台获取）

setwd("D:/R_input&#​34;)

source("vennDia.R&#​34;)

### 导入样品

otu_16s<-read.csv("D://test_otu.csv&#​34;,row.names = 1)

design<-read.csv("D://test_design.csv&#​34;)

#选择otu和物种

otu_filter_16s<-otu_16s[,1:27]

tax_filter_16s<-otu_16s[,28:34]

##edgeR鉴定差异OTU

model_matsoil_16s <- model.matrix(~CBBpare, data=design)

edgeR_16s_soil_fsyst <- DGEList(counts=otu_filter_16s, group=design$CBBpare, genes=tax_filter_16s)

edgeR_16s_soil_fsyst <- calcNormFactors(edgeR_16s_soil_fsyst)

dge_soilfsyst_16s <- estimateGLMRobustDisp(edgeR_16s_soil_fsyst, design=model_matsoil_16s)

fit_soilfsyst_16s <- glmFit(dge_soilfsyst_16s, design = model_matsoil_16s)

lrt_soilfsyst_16s <- glmLRT(fit_soilfsyst_16s, coef=2:3)

fsyst_soil_16s <- topTags(lrt_soilfsyst_16s, n=Inf, p.value=0.05)

fsyst_soil_16s <- fsyst_soil_16s$table

## 鉴定指示物种indicator

indic_soil_16s <- as.data.frame(t(otu_filter_16s))

indic_soil_groups_16s <- design$CBBpare

length(unique(indic_soil_groups_16s))

set.seed(8046)

indicatorsp_soil_16s <- multipatt(indic_soil_16s,indic_soil_groups_16s,func = "r.g&#​34;,control=how(nperm=99))

summary(indicatorsp_soil_16s,alpha=1,indvalcBBp=T)

indic_soil_df_16s <- indicatorsp_soil_16s$sign

## indicator阈值

net_16s <- as.matrix(indic_soil_df_16s[which(indic_soil_df_16s$p.value < 0.05),])

# 获得indicator物种表

indicator_taxa<-subset(otu_16s,rownames(otu_16s) %in% rownames(net_16s))

indicator<-cbind(net_16s ,indicator_taxa)

#导出

#write.csv(indicator,file=&#​34;D://a_rice/indic_aggregate.csv" )

## 交叉筛选获得存在显著丰度差异的indicator

indic_edge_16s_soil <- intersect(rownames(net_16s),rownames(fsyst_soil_16s))

##韦恩图可视化交叉物种

venndiagram(rownames(net_16s),rownames(fsyst_soil_16s),type=2,

            printsub=F,lcol="black&#​34;,tcol="black&#​34;,diacol="black&#​34;,lines="black&#​34;,

            labels=c("Indicator Species Analysis&#​34;,"edgeR Analysis&#​34;),title="Sensitive \nSoil Bacteria OTUs")

###获得显著丰度差异的indicator表

indicator_16s<-subset(net_16s,rownames(net_16s) %in%indic_edge_16s_soil)

indicator_taxa<-subset(otu_16s,rownames(otu_16s) %in% rownames(indicator_16s))

indicator_taxa<-cbind(indicator_16s,indicator_taxa)

#导出

write.csv(indicator_taxa,file="D://root/bb/indic_fungi_bulk.csv" )

#转化数据框

indicator_16s<-as.data.frame(indicator_16s)

#构建网络文件

soil_bipartite_16s <- data.frame(frBB= c(rep("AA_bulk&#​34;,length(which(net_16s[,"s.bulk_AA&#​34;]==1))),

                                         rep("AA_root&#​34;,length(which(net_16s[,"s.root_AA&#​34;]==1))),

                                         rep("AA_rhizo&#​34;,length(which(net_16s[,"s.rhizo_AA&#​34;]==1))),

                                         rep("BB_bulk&#​34;,length(which(net_16s[,"s.bulk_BB&#​34;]==1))),

                                         rep("BB_root&#​34;,length(which(net_16s[,"s.root_BB&#​34;]==1))),

                                         rep("BB_rhizo&#​34;,length(which(net_16s[,"s.rhizo_BB&#​34;]==1))),

                                         rep("CC_bulk&#​34;,length(which(net_16s[,"s.bulk_CC&#​34;]==1))),

                                         rep("CC_root&#​34;,length(which(net_16s[,"s.root_CC&#​34;]==1))),

                                         rep("CC_rhizo&#​34;,length(which(net_16s[,"s.rhizo_CC&#​34;]==1)))),

                                 to= c(rownames(net_16s)[which(net_16s[,"s.bulk_AA&#​34;]==1)],

                                       rownames(net_16s)[which(net_16s[,"s.rhizo_AA&#​34;]==1)],

                                       rownames(net_16s)[which(net_16s[,"s.root_AA&#​34;]==1)],

                                       rownames(net_16s)[which(net_16s[,"s.bulk_BB&#​34;]==1)],

                                       rownames(net_16s)[which(net_16s[,"s.rhizo_BB&#​34;]==1)],

                                       rownames(net_16s)[which(net_16s[,"s.root_BB&#​34;]==1)],

                                       rownames(net_16s)[which(net_16s[,"s.bulk_CC&#​34;]==1)],

                                       rownames(net_16s)[which(net_16s[,"s.rhizo_CC&#​34;]==1)],

                                       rownames(net_16s)[which(net_16s[,"s.root_CC&#​34;]==1)]),

                                 r= c(net_16s[which(net_16s[,"s.bulk_AA&#​34;]==1),"stat&#​34;],

                                      net_16s[which(net_16s[,"s.rhizo_AA&#​34;]==1),"stat&#​34;],

                                      net_16s[which(net_16s[,"s.root_AA&#​34;]==1),"stat&#​34;],

                                      net_16s[which(net_16s[,"s.bulk_BB&#​34;]==1),"stat&#​34;],

                                      net_16s[which(net_16s[,"s.rhizo_BB&#​34;]==1),"stat&#​34;],

                                      net_16s[which(net_16s[,"s.root_BB&#​34;]==1),"stat&#​34;],

                                      net_16s[which(net_16s[,"s.bulk_CC&#​34;]==1),"stat&#​34;],

                                      net_16s[which(net_16s[,"s.rhizo_CC&#​34;]==1),"stat&#​34;],

                                      net_16s[which(net_16s[,"s.root_CC&#​34;]==1),"stat&#​34;]))

#导出网络文件，gephi可视化，from和to修改为Source和Target

write.csv(soil_bipartite_16s ,file="D://aa_soil_bipartite_16s.csv")