MARco: Microbiome Analysis RcodeDB
Liu, Po-Yu, Ph.D., National Taiwan University College of Medicine
MARco
MARco is an R code and function set for microbiome analysis with less pain to generate basic plots of microbial ecology, including alpha diversity boxplots, violin plots, beta diversity ordination (PCoA based on Bray-Curtis dissimilarity) with ADONIS test, and several statistical and visualization tools (see below). MARco also contains several discrete color sets for applying to your grouped plots. For advanced analysis, co-occurrence network analysis is also able to conduct by MARco.
All updated codes will be released on https://github.com/poyuliu/MARco.
Please cites:
Liu P-Y. MARco: Microbiome Analysis RcodeDB (2021). 10.5281/zenodo.4589897.
Latest updates: February 10, 2022
Required packages:
vegan with functions for most community ecology analyses
fossil for calculating Chao1 index
vioplot for ploting violin plots
igraph for constructing network
SpiecEasi for calculating SparCC correlation, how to install: https://github.com/zdk123/SpiecEasi
DESeq2 RNA-Seq based differential expression test, see details on Bioconductor: https://bioconductor.org/packages/release/bioc/html/DESeq2.html
pheatmap for plot heatmap
MASS for LDA (linear discriminant analysis) effect size (LEfSe)
gage for GAGE (Generally Applicable Gene-set Enrichment) analysis
cluster and clusterSim for Enterotyping
Load MARco package:
library(MARco)Colorsets
Examples
Regular analyses
Example dataset:
Suzuki T., et al. Host genetic determinants of the gut microbiota of wild mice. Molecular ecology 28 (13), 3197-3207 (2019)
Alpha diversity
alpha.diversity <- alphdiv(data = otu,group = groups,test = T) # calculate alpha diversity indices (Richness, Shannon, Simpson, Chao1)
par(mfrow=c(1,4))
plot.alpha(alpha = alpha.diversity,group = groups,colset = colset.d.4,names = levels(groups),test = T) # output boxplots
vioplot.alpha(alpha = alpha.diversity,group = groups,colset = colset.d.4,names = levels(groups),test = T) # output violin plots
Relative abundance barplots
taxaQ2output <- read.delim("suzuki_taxonomy/taxonomy.tsv")
TaxaLvABar(otu = otu,taxafromQ2tsv = taxaQ2output,level = 2,colsets = colset.d.12)
TaxaLvABar(otu = otu,taxafromQ2tsv = taxaQ2output,level = 3,colsets = colset.d.12)
TaxaLvABar(otu = otu,taxafromQ2tsv = taxaQ2output,level = 5,colsets = colset.d.12)
par(mfrow=c(1,3))
TaxaLvABar(otu = otu,taxafromQ2tsv = taxaQ2output,level = 2,colsets = colset.d.12,grouped = TRUE,groups=groups)
TaxaLvABar(otu = otu,taxafromQ2tsv = taxaQ2output,level = 3,colsets = colset.d.12,grouped = TRUE,groups=groups)
TaxaLvABar(otu = otu,taxafromQ2tsv = taxaQ2output,level = 5,colsets = colset.d.12,grouped = TRUE,groups=groups)
Beta diversity
pcoa <- ordNtest(data = otu,group = groups)
plot.ordi(ordN = pcoa,group = groups,colset = colset.d.4,
ordiellipse = T,ordispider = T,
legend = levels(groups),test = T)
DESeq2 of 2-group comparison
(Statistic/Taxonomy tables &) Volcano plot (Fold change >2, adjust P < 0.05)
DESeq2.FL.MP <- DESq2(data = otu.2,group = group.2)
#DESeq2.sig.table(DESeq.res = head(DESeq2.FL.MP),taxa = taxa ,caption = "MP vs. FL")
volcano.plot(DESeq.res = DESeq2.FL.MP,colset = colset.d.4)
Heatmap with p values
G2.heatmap(otu.sig = otu.2[which(DESeq2.FL.MP$padj<0.0001),],groups = group.2,
p.val = DESeq2.FL.MP$padj[which(DESeq2.FL.MP$padj<0.0001)],
taxa.sig = taxa[match(rownames(otu.2),rownames(taxa)),][which(DESeq2.FL.MP$padj<0.0001),],
otu = otu.2,metadata = as.data.frame(cbind(SampleID=colnames(otu.2))),colset = colset.d.4[1:2])
Co-occurrence network
# prepare data and construct a correlation matrix with SparCC
net.o <- prevalance(data = otu,prevalance = 0.4) # prevalenve ≥ 0.4
ww <- sparccboot(t(net.o),R=9,ncpus = 6) # R: bootsrtapping numbers
w.boot <- pval.sparccboot(ww)
w.mb <- w.boot$cors
p.mb <- w.boot$pvals
w.mb <- boot2matrix(w.mb,net.o)
p.mb <- boot2matrix(p.mb,net.o)
diag(w.mb) <- 0
w.mb[p.mb > 0.05] <- 0 # significance cutoff
rownames(w.mb) <- colnames(w.mb) <- rownames(net.o)or just run a one-step function:
w.mb <- network.pipeline(otu=net.o,prevalence=0.1,alpha=0.05,bootstrap=99,cpu=12)# construct and plot network
g.mb <- make_network(w = w.mb,
data = net.o,
both = F, # only positive correlation
cutoff = quantile(w.mb[w.mb>0],0.95), # correlation coefficients / weights cutoff for removing nodes
degree = 0, # connection degree cutoff for revoming nodes
size = 1.5,
clustering = T)
V(g.mb$graph)$color <- c(colset.d.12,"#cacaca")[g.mb$clusters] # label clusters with color (13 clusters)
plot_network(g.mb,negative = F,pos.col="#5d5d5d55")
Advanced analyses
Example dataset:
Liu P.-Y., et al. Variations in Gut Microbiota of Siberian Flying Squirrels Correspond to Seasonal Phenological Changes in Their Hokkaido Subarctic Forest Ecosystem. Microbial ecology 78, 223–231 (2019)
LEfSe
Linear discriminant analysis (LDA) Effect Size (version 1.0; only for 2-group comparison)
Lefse.temp <- lefse(data = pvo.otu,taxa = pvo.taxa,Groups = meteor$hot,Lefse.factor = 1,FDR = FALSE,alpha = 0.1,LDAscore = 2,plot = TRUE,col = colset.d.2)
print(Lefse.temp$LEfSe)## taxa
## 1 Firmicutes; Clostridia; Clostridiales; Family XIII; [Eubacterium] brachy group
## 13 Firmicutes; Clostridia; Clostridiales; Family XIII; [Eubacterium] brachy group; uncultured organism
## 15 Firmicutes; Clostridia; Clostridiales; Lachnospiraceae; [Ruminococcus] torques group; uncultured bacterium
## 3 Firmicutes; Clostridia; Clostridiales; Lachnospiraceae; [Ruminococcus] torques group
## 10 Firmicutes; Clostridia; Clostridiales; Ruminococcaceae; Ruminococcaceae UCG-010
## 22 Firmicutes; Clostridia; Clostridiales; Lachnospiraceae; Moryella; uncultured bacterium
## 9 Firmicutes; Clostridia; Clostridiales; Lachnospiraceae; Moryella
## 8 Firmicutes; Clostridia; Clostridiales; Lachnospiraceae; Lachnospiraceae UCG-006
## 21 Firmicutes; Clostridia; Clostridiales; Lachnospiraceae; Lachnospiraceae UCG-006; uncultured bacterium
## 11 Firmicutes; Clostridia; Clostridiales; Ruminococcaceae; Ruminococcus 1
## 25 Firmicutes; Clostridia; Clostridiales; Ruminococcaceae; Ruminococcaceae UCG-014; uncultured Clostridiales bacterium
## 7 Firmicutes; Clostridia; Clostridiales; Lachnospiraceae; Lachnospiraceae UCG-004
## 20 Firmicutes; Clostridia; Clostridiales; Lachnospiraceae; Lachnospiraceae UCG-004; uncultured bacterium
## 17 Firmicutes; Clostridia; Clostridiales; Lachnospiraceae; Lachnoclostridium; uncultured organism
## 26 Firmicutes; Clostridia; Clostridiales; Ruminococcaceae; Ruminococcaceae UCG-014; unidentified
## 23 Firmicutes; Clostridia; Clostridiales; Ruminococcaceae; Ruminococcaceae UCG-010; uncultured bacterium
## 12 Firmicutes; Clostridia; Clostridiales; Christensenellaceae; Christensenellaceae R-7 group; uncultured rumen bacterium
## LD1 EffectSize
## 1 -353.1287 -2.547933
## 13 -353.1287 -2.547933
## 15 -1143.1654 -3.058109
## 3 -1143.1654 -3.058109
## 10 -5360.5272 -3.729208
## 22 127.0562 2.103996
## 9 127.0562 2.103996
## 8 324.0076 2.510555
## 21 324.0076 2.510555
## 11 325.1998 2.512150
## 25 411.5078 2.614378
## 7 798.2768 2.902154
## 20 798.2768 2.902154
## 17 1276.2539 3.105937
## 26 1316.7528 3.119504
## 23 4228.3421 3.626170
## 12 4815.3482 3.682628Random Forest
Random Forest regression
Random Forest classification for model prediction and feature selection
rownames(pvo.otu) <- paste0("ASV_",rownames(pvo.otu)) # A number of the first character of feature name is not allowed
rownames(otu) <- paste0("ASV_",rownames(otu))
meteor$hot <- as.factor(ifelse(meteor$Temperature>20,"hot","cool"))
RF <- runRF(otu_table_scaled = log1p(pvo.otu),Y = hot,metadata = meteor,model = "classification")AUROC <- RF2ROC(RF,plot=T)
pvo.otu.2 <- pvo.otu[match(RF$SortFeature$features[1:10],rownames(pvo.otu)),] # select top 10 features
RF2 <- runRF(otu_table_scaled = log1p(pvo.otu.2),Y = hot,metadata = meteor,model = "classification")AUROC2 <- RF2ROC(RF2,plot=T)
Validation cohort cross validation
newd <- log1p(otu[,c(1:6,26:29)])
newdgroup <- as.factor(c(rep("hot",6),rep("cool",4)))
validation.out <- VCCV(RFmodel = RF2,newdata = newd,newdataGroup = newdgroup,RFROC = AUROC2,plot = T)
Longitudinal analysis
Enterotyping (PAM method)
Reference: Gomez, A. et al. Plasticity in the Human Gut Microbiome Defies Evolutionary Constraints. mSphere 4, doi:10.1128/mSphere.00271-19 (2019).
Tutorial: https://enterotype.embl.de/enterotypes.html
source("http://www.lifescipy.net/RcodeDB/FN_enterotyping.R") # Enterotyping
ET <- enterotyping(otu = genus.otu,remove.noise = T,rm.pct = 0.05)
plot.ET(enterotyping = ET,colset = colset.d.4,ordiellipse = T,ordispider = T,test = T)
ET1 <- box.ET(enterotyping = ET,clusterN = 1,out = T,colset = colset.d.4,outline=F)
ET2 <- box.ET(enterotyping = ET,clusterN = 2,out = T,colset = colset.d.4,outline=F)Visualize enterotypes using network
Fitting variable gradients
Distance-based redundancy analysis (db-RDA) principle and examples: https://archetypalecology.wordpress.com/2018/02/21/distance-based-redundancy-analysis-db-rda-in-r/
head(pvo.otu)[,1:7] # view gut microbiota OTU table of Siberian flying squirrels## PVO01 PVO02 PVO03 PVO04 PVO05 PVO06 PVO07
## de5729e94e60bd7f27bef251c0561446 0 54 0 4 0 0 25
## 11b2c89b911fc23c91bac3c019bbd970 0 0 0 0 3 4 0
## 11a1be5c2fab130fa99d87a414c6cc4f 0 0 0 0 41 21 52
## 2ec026b88769ccfbb75ea6ba3ac4c07b 0 0 0 0 0 0 0
## 6705ca1e52fe7706e8916a752a502ac9 0 0 0 0 0 0 0
## 51c5d0318d5fab9825bbb5317da2cb2c 0 0 0 0 0 0 0head(meteor) # view environmental factor metadata## Temperature Precipitation NDVI
## PVO01 20.085714 0.0000000 0.7748209
## PVO02 9.185714 3.4285714 0.3811293
## PVO03 22.400000 0.8571429 0.8465717
## PVO04 8.914286 2.2857143 0.3811293
## PVO05 22.400000 0.8571429 0.8465717
## PVO06 8.914286 2.2857143 0.3811293dbRDA with continuous variable gradients fitting
dbRDA <- dbrdaNtest(data = pvo.otu,ENV = meteor)
print(dbRDA$permucapscale)## Permutation test for capscale under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: vegan::capscale(formula = t(data) ~ Temperature + Precipitation + NDVI, data = ENV, distance = "bray")
## Df SumOfSqs F Pr(>F)
## Temperature 1 0.6118 1.7237 0.003 **
## Precipitation 1 0.3612 1.0177 0.410
## NDVI 1 0.5464 1.5394 0.006 **
## Residual 15 5.3245
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot.dbRDA(ordN = dbRDA, select = meteor$Temperature, colset = WhBuRd.g(20), fit.p = 0.1, legend = TRUE,legend.labels = "Temperature (C)")
terrain20 <- substr(terrain.colors(20,rev = T),1,7) # color for NDVI
plot.dbRDA(ordN = dbRDA, select = meteor$NDVI, colset = terrain20, fit.p = 1, legend = TRUE, legend.labels = "Mean NDVI")
PCoA with continuous variable gradients fitting
pvo <- ordEnvfit(data = pvo.otu,ENV = meteor) # compute PCoA and environmental factor fitting
print(pvo$fit)##
## ***VECTORS
##
## Dim1 Dim2 r2 Pr(>r)
## Temperature 0.92993 -0.36774 0.2875 0.079 .
## Precipitation -0.56959 -0.82193 0.0786 0.521
## NDVI 0.75932 -0.65072 0.4148 0.014 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999plot.ordi.env(ordN = pvo,select = meteor$Temperature,colset = WhBuRd.g(20),fit.p = 0.1,legend = TRUE,legend.labels = "Temperature (C)")
terrain20 <- substr(terrain.colors(20,rev = T),1,7) # color for NDVI
plot.ordi.env(ordN = pvo,select = meteor$NDVI,colset = terrain20,fit.p = 1,legend = TRUE,legend.labels = "Mean NDVI")
OTU correlations
Simple Env correlation
Temperature <- meteor$Temperature
par(mfrow=c(1,2))
SimpleEnvCor <- EnvCorOTU(data = pvo.otu,ENV = Temperature,summary.plot = T)
SimpleEnvCor.loglog <- EnvCorOTU(data = pvo.otu,ENV = Temperature,logdata = TRUE,logENV = TRUE,summary.plot = T)
head(SimpleEnvCor$otu.cor.env)[,1:7]## PVO01 PVO02 PVO03 PVO04 PVO05 PVO06 PVO07
## f968c28ae8e138a530f83f6962e6bd15 0 0 0 17 0 14 0
## 3f6f01010ba6cc374032b85e6fbddf28 0 0 0 7 0 13 0
## 987d5eea139a26d57dec1c94a4f8b55a 0 13 0 0 0 0 0
## 9425c1bd81b35a5e454f3e7104b1ddbd 0 16 0 48 0 0 0
## ef75dc287aa3ae0f09c136e75f433840 0 32 0 10 0 24 17
## 8cb05dce16c998b05ee0849cd75357b3 0 0 0 22 0 17 0taxa[match(rownames(SimpleEnvCor$otu.cor.env),rownames(taxa)),][1:3,]## V1 V2 V3
## f968c28ae8e138a530f83f6962e6bd15 D_0__Bacteria D_1__Firmicutes D_2__Clostridia
## 3f6f01010ba6cc374032b85e6fbddf28 Unassigned Unassigned Unassigned
## 987d5eea139a26d57dec1c94a4f8b55a D_0__Bacteria D_1__Firmicutes D_2__Clostridia
## V4 V5
## f968c28ae8e138a530f83f6962e6bd15 D_3__Clostridiales D_4__Lachnospiraceae
## 3f6f01010ba6cc374032b85e6fbddf28 Unassigned Unassigned
## 987d5eea139a26d57dec1c94a4f8b55a D_3__Clostridiales D_4__Ruminococcaceae
## V6
## f968c28ae8e138a530f83f6962e6bd15 D_5__Lachnospiraceae NK4A136 group
## 3f6f01010ba6cc374032b85e6fbddf28 Unassigned
## 987d5eea139a26d57dec1c94a4f8b55a D_5__Ruminococcaceae UCG-008
## V7
## f968c28ae8e138a530f83f6962e6bd15 D_6__Lachnospiraceae bacterium COE1
## 3f6f01010ba6cc374032b85e6fbddf28 Unassigned
## 987d5eea139a26d57dec1c94a4f8b55a D_6__uncultured bacteriumhead(SimpleEnvCor$sig.corr)## f968c28ae8e138a530f83f6962e6bd15.cor 3f6f01010ba6cc374032b85e6fbddf28.cor
## -0.6280158 -0.4742978
## 987d5eea139a26d57dec1c94a4f8b55a.cor 9425c1bd81b35a5e454f3e7104b1ddbd.cor
## -0.4915133 -0.5582419
## ef75dc287aa3ae0f09c136e75f433840.cor 8cb05dce16c998b05ee0849cd75357b3.cor
## -0.6536772 -0.6656628Coordinates x Env correlation
CoordinateEnvCor <- PCENVcorOTU(data=pvo.otu,Ord=pvo,ENV=Temperature,plot.pc.env = TRUE)
head(CoordinateEnvCor$otu.cor.pc,n = 3)## PVO01 PVO02 PVO03 PVO04 PVO05 PVO06 PVO07
## 4aa13c19e102f5629f07c966da3f2dc8 0 0 0 0 0 0 0
## bc24f4b609982c37de4be846c31e0810 0 0 0 0 0 0 0
## 46a7467c833b41b93b08e2c9d848db27 0 0 0 0 0 0 0
## PVO08 PVO09 PVO10 PVO11 PVO12 PVO13 PVO14
## 4aa13c19e102f5629f07c966da3f2dc8 17 4 6 0 0 0 0
## bc24f4b609982c37de4be846c31e0810 13 7 7 0 0 0 0
## 46a7467c833b41b93b08e2c9d848db27 13 3 5 0 0 0 0
## PVO15 PVO16 PVO17 PVO18 PVO19
## 4aa13c19e102f5629f07c966da3f2dc8 0 0 0 0 0
## bc24f4b609982c37de4be846c31e0810 0 0 0 0 0
## 46a7467c833b41b93b08e2c9d848db27 0 0 0 2 0head(CoordinateEnvCor$sig.corr,n = 5)## 4aa13c19e102f5629f07c966da3f2dc8.cor bc24f4b609982c37de4be846c31e0810.cor
## 0.6321743 0.7181039
## 46a7467c833b41b93b08e2c9d848db27.cor 822c2a4e59a643c0ea64e197f36a49f5.cor
## 0.6619272 0.6913038
## a03a5d4008bc18faf78f4bccaa007cf9.cor
## 0.7356530head(CoordinateEnvCor$sig.p,n = 5)## 4aa13c19e102f5629f07c966da3f2dc8 bc24f4b609982c37de4be846c31e0810
## 0.0036837252 0.0005349764
## 46a7467c833b41b93b08e2c9d848db27 822c2a4e59a643c0ea64e197f36a49f5
## 0.0020206999 0.0010452185
## a03a5d4008bc18faf78f4bccaa007cf9
## 0.0003310058CoordinateEnvCor <- PCENVcorOTU(data=pvo.otu,Ord=pvo,ENV=Temperature,summary.plot = T)
Functional analysis
KO2pathway
Reformat KEGG Orthology (KO) annoated data to pathways or modules
Reformat
(example) original data:
head(picrust2) # functional data input format## Exp.1.1 Exp.1.2 Exp.2.1 Exp.2.2
## K00001 1639.66 1964.01 774.87 2222.67
## K00003 6659.17 3582.34 3016.57 3136.67
## K00004 596.32 1399.94 359.48 542.94
## K00005 1499.32 2188.62 3713.32 2965.53
## K00007 213.25 96.12 288.95 213.02
## K00008 417.50 1226.59 1595.25 1601.50Convert to metabolic pathways
gs.table <- KO2path(predtable = picrust2,GeneSet = "metabolic") # use "metabolic pathway" gene sets
# other gene set options: "pathway" [all KEGG pathways], "module" [KEGG modules]
head(gs.table)## Exp.1.1 Exp.1.2 Exp.2.1
## ko00010 Glycolysis / Gluconeogenesis 166939.18 149334.41 229284.9
## ko00020 Citrate cycle (TCA cycle) 128411.22 95628.92 151968.8
## ko00030 Pentose phosphate pathway 102397.70 91167.36 148110.7
## ko00040 Pentose and glucuronate interconversions 26246.29 41395.86 75169.8
## ko00051 Fructose and mannose metabolism 75283.34 80886.60 139866.0
## ko00052 Galactose metabolism 52751.23 61476.12 85308.7
## Exp.2.2
## ko00010 Glycolysis / Gluconeogenesis 237922.58
## ko00020 Citrate cycle (TCA cycle) 146169.77
## ko00030 Pentose phosphate pathway 147281.08
## ko00040 Pentose and glucuronate interconversions 69958.67
## ko00051 Fructose and mannose metabolism 137782.92
## ko00052 Galactose metabolism 98576.77gs.table.p <- KO2path(predtable = picrust2,GeneSet = "metabolic",prop = TRUE) # present with proportion
head(gs.table.p)## Exp.1.1 Exp.1.2
## ko00010 Glycolysis / Gluconeogenesis 0.014233946 0.015090086
## ko00020 Citrate cycle (TCA cycle) 0.010948888 0.009663202
## ko00030 Pentose phosphate pathway 0.008730865 0.009212366
## ko00040 Pentose and glucuronate interconversions 0.002237871 0.004183008
## ko00051 Fructose and mannose metabolism 0.006418979 0.008173506
## ko00052 Galactose metabolism 0.004497795 0.006212098
## Exp.2.1 Exp.2.2
## ko00010 Glycolysis / Gluconeogenesis 0.014244095 0.014825073
## ko00020 Citrate cycle (TCA cycle) 0.009440909 0.009107910
## ko00030 Pentose phosphate pathway 0.009201231 0.009177156
## ko00040 Pentose and glucuronate interconversions 0.004669849 0.004359159
## ko00051 Fructose and mannose metabolism 0.008689038 0.008585321
## ko00052 Galactose metabolism 0.005299718 0.006142367Convert to KEGG modules
md.table <- KO2path(predtable = picrust2,GeneSet = "module")
head(md.table)## Exp.1.1
## M00001 Glycolysis (Embden-Meyerhof pathway), glucose => pyruvate 64881.43
## M00002 Glycolysis, core module involving three-carbon compounds 45976.57
## M00003 Gluconeogenesis, oxaloacetate => fructose-6P 55561.44
## M00004 Pentose phosphate pathway (Pentose phosphate cycle) 28374.89
## M00005 PRPP biosynthesis, ribose 5P => PRPP 7468.08
## M00006 Pentose phosphate pathway, oxidative phase, glucose 6P => ribulose 5P 6025.82
## Exp.1.2
## M00001 Glycolysis (Embden-Meyerhof pathway), glucose => pyruvate 62881.83
## M00002 Glycolysis, core module involving three-carbon compounds 38598.89
## M00003 Gluconeogenesis, oxaloacetate => fructose-6P 48148.76
## M00004 Pentose phosphate pathway (Pentose phosphate cycle) 36403.99
## M00005 PRPP biosynthesis, ribose 5P => PRPP 4574.15
## M00006 Pentose phosphate pathway, oxidative phase, glucose 6P => ribulose 5P 12192.02
## Exp.2.1
## M00001 Glycolysis (Embden-Meyerhof pathway), glucose => pyruvate 83670.64
## M00002 Glycolysis, core module involving three-carbon compounds 55832.10
## M00003 Gluconeogenesis, oxaloacetate => fructose-6P 66851.20
## M00004 Pentose phosphate pathway (Pentose phosphate cycle) 54434.79
## M00005 PRPP biosynthesis, ribose 5P => PRPP 6409.75
## M00006 Pentose phosphate pathway, oxidative phase, glucose 6P => ribulose 5P 15610.86
## Exp.2.2
## M00001 Glycolysis (Embden-Meyerhof pathway), glucose => pyruvate 90545.76
## M00002 Glycolysis, core module involving three-carbon compounds 57662.32
## M00003 Gluconeogenesis, oxaloacetate => fructose-6P 70362.43
## M00004 Pentose phosphate pathway (Pentose phosphate cycle) 53552.79
## M00005 PRPP biosynthesis, ribose 5P => PRPP 6338.71
## M00006 Pentose phosphate pathway, oxidative phase, glucose 6P => ribulose 5P 15759.05Annotation
Annotate with upper levels
lv.table <- path.lvs(gs.table)
head(lv.table)## Level_1 Level_2
## ko01100 1. Metabolism 1.0 Global and overview maps
## ko01110 1. Metabolism 1.0 Global and overview maps
## ko01120 1. Metabolism 1.0 Global and overview maps
## ko01200 1. Metabolism 1.0 Global and overview maps
## ko01210 1. Metabolism 1.0 Global and overview maps
## ko01212 1. Metabolism 1.0 Global and overview maps
## Level_3 Exp.1.1 Exp.1.2
## ko01100 Metabolic pathways 2821991.84 2373497.91
## ko01110 Biosynthesis of secondary metabolites 1298830.39 1115126.21
## ko01120 Microbial metabolism in diverse environments 830717.90 660121.04
## ko01200 Carbon metabolism 482057.24 357656.68
## ko01210 2-Oxocarboxylic acid metabolism 153679.26 98524.16
## ko01212 Fatty acid metabolism 68005.65 87549.93
## Exp.2.1 Exp.2.2
## ko01100 3963999.0 3890785.3
## ko01110 1734707.4 1717578.4
## ko01120 1094990.9 1117300.4
## ko01200 587878.1 594110.2
## ko01210 164314.5 155855.9
## ko01212 122264.5 129425.9Summarize to level 2
lv2.table <- path.lvs(gs.table,aggregate.lv = 2)
head(lv2.table)## Exp.1.1 Exp.1.2 Exp.2.1
## 1.0 Global and overview maps 6309559.5 5179221.0 8413080.2
## 1.1 Carbohydrate metabolism 1331518.4 1263577.9 2072358.5
## 1.10 Biosynthesis of other secondary metabolites 187888.5 159669.5 228874.4
## 1.11 Xenobiotics biodegradation and metabolism 196861.4 194918.0 308046.2
## 1.2 Energy metabolism 755372.2 529790.4 950335.1
## 1.3 Lipid metabolism 214247.1 262753.0 411775.2
## Exp.2.2
## 1.0 Global and overview maps 8339305.2
## 1.1 Carbohydrate metabolism 2092350.6
## 1.10 Biosynthesis of other secondary metabolites 237803.3
## 1.11 Xenobiotics biodegradation and metabolism 321838.1
## 1.2 Energy metabolism 914836.4
## 1.3 Lipid metabolism 428759.5FSEA
Functional Sets Enrichment Analysis
Methodology reference:
Liu, P-Y, et al. Body-size Scaling is Related to Gut Microbial Diversity, Metabolism and Dietary niche of Arboreal folivorous flying Squirrels. Scientific reports 10, 1-12 (2020)
Within group test
Clutch <- factor(c("clutch.1","clutch.1","clutch.2","clutch.2"))
FSEA.metabolic <- FSEA(Fdata = picrust2,Groups = Clutch,FunctionalSets = "metabolic",SortLv = 2,Test = "within",barcol = colset.d.4[1:2],plot = TRUE)
FSEA.module <- FSEA(Fdata = picrust2,Groups = Clutch,FunctionalSets = "module",SortLv = 3,Test = "within",barcol = colset.d.4[1:2],plot = TRUE)
Between group test
FS4gp <- substr(colnames(tax4fun),1,3)
FS4gp <- gsub("PVO","Small FS",FS4gp)
FS4gp <- gsub("TX0","Medium FS",FS4gp)
FS4gp <- gsub("PAL","Large FS",FS4gp)
FS4gp <- gsub("PPG","Large FS",FS4gp)
FS4gp <- as.factor(FS4gp)
par(mfcol=c(3,1))
FS.FSEA.between.metabolic <- FSEA(Fdata = tax4fun,Groups = FS4gp,FunctionalSets = "metabolic",Test = "between",plot = T,barcol = colset.d.2)
Predict metabolome profile
KO2metabolite(petaurista)