pathway.compare: Compare (kegg) pathway abundance

In metamicrobiomeR: Microbiome Data Analysis & Meta-Analysis with GAMLSS-BEZI & Random Effects

Description

This is a slightly modified version of the taxa.compare function. It compares pathway abundance generated from PICRUSt analysis between groups using different methods (apply to pathway summary tables already merged to mapping file and put in a list (e.g.level 1, 2 and 3)). Specifically, it compares relative abundances of bacterial functional pathways at all levels using GAMLSS or LM/LMEM and compares of log(absolute abundances) of bacterial functional pathways at all levels using LM/LMEM.

Usage

pathway.compare(
  pathtab,
  mapfile,
  sampleid = "sampleid",
  pathsum = "rel",
  stat.med = "gamlss",
  transform = "none",
  comvar,
  adjustvar,
  personid = "personid",
  longitudinal = "yes",
  p.adjust.method = "fdr",
  percent.filter = 0.05,
  relabund.filter = 5e-05,
  pooldata = FALSE
)

Arguments

pathtab	list of pathway abundance table of all levels.
mapfile	mapping file or file containing covariates.
sampleid	variable containing sample id to be matched between pathway abundance table and mapping file.
pathsum	type of abundance to be compared. Options are "rel" for relative abundance or "log" for log of absolute abundance. Default is "rel".
stat.med	statistical method for comparison. stat.med can be "lm" for LM/LMEM (usable for both pathsum="rel" or "log") or "gamlss" for GAMLSS with BEZI family (gamlss only make sense if pathsum="rel" ).
transform	transformation of relative abundance data. Options are "none" for no transformation, "asin.sqrt" for arcsine transformation, "logit" for logit transformation. Default is "none".
comvar	main variable for comparison.
adjustvar	variables to be adjusted.
personid	name of variable for person id in mapping file (applicable for longitudinal data)
longitudinal	whether the data is longitudinal. Default is "yes".
p.adjust.method	method for multiple testing adjustment. Available options are those of the p.adjust function. Default is "fdr".
percent.filter	prevalence threshold (the percentage of number of samples the taxa/pathway available). Default is 0.05.
relabund.filter	relative abundance threshold (the minimum of the average relative abundance for a taxa/pathway to be retained). Default is 0.00005.
pooldata	whether the data is pooled from multiple studies. Default is FALSE.

Value

matrice of coefficients, standard errors, p-values and multiple testing adjusted p-values of all variables in the models.

Examples

#Load Bangladesh pathway and metadata
data(kegg.12)
data(covar.rm)
# Comparison of pathway relative abundances for some first pathways of level 1 only
# and assuming crosssectional data (to save running time)
path1<-pathway.compare(pathtab=list(kegg.12[[1]][, 1:2]),
mapfile=covar.rm,sampleid="sampleid",pathsum="rel", stat.med="gamlss",
comvar="gender",adjustvar=c("age.sample","bf"), longitudinal="no",
p.adjust.method="fdr", percent.filter=0.05,relabund.filter=0.00005)
taxcomtab.show(taxcomtab=path1$l1, sumvar="path",tax.lev="l2",
tax.select="none", showvar="genderMale", p.adjust.method="fdr",p.cutoff=1)

Example output

Loading required package: gamlss
Loading required package: splines
Loading required package: gamlss.data

Attaching package: ‘gamlss.data’

The following object is masked from ‘package:datasets’:

    sleep

Loading required package: gamlss.dist
Loading required package: MASS
Loading required package: nlme
Loading required package: parallel
 **********   GAMLSS Version 5.2-0  ********** 
For more on GAMLSS look at https://www.gamlss.com/
Type gamlssNews() to see new features/changes/bug fixes.

[1] 1
******************************************************************
Family:  c("BEZI", "Zero Inflated Beta") 

Call:  gamlss::gamlss(formula = stats::as.formula(paste(pathname[i],  
    paste(c(comvar, adjustvar), collapse = "+"), sep = "~")),  
    family = BEZI, data = testdat, trace = FALSE) 

Fitting method: RS() 

------------------------------------------------------------------
Mu link function:  logit
Mu Coefficients:
                   Estimate Std. Error  t value Pr(>|t|)    
(Intercept)       -2.246380   0.019984 -112.410  < 2e-16 ***
genderMale         0.004272   0.014308    0.299    0.765    
age.sample         0.008404   0.001235    6.804 1.76e-11 ***
bfNon_exclusiveBF -0.022727   0.023196   -0.980    0.327    
bfNo_BF            0.042615   0.038138    1.117    0.264    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
Sigma link function:  log
Sigma Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.36087    0.04449   120.5   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
Nu link function:  logit 
Nu Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   -21.43     866.54  -0.025     0.98

------------------------------------------------------------------
No. of observations in the fit:  995 
Degrees of Freedom for the fit:  7
      Residual Deg. of Freedom:  988 
                      at cycle:  7 
 
Global Deviance:     -4914.618 
            AIC:     -4900.618 
            SBC:     -4866.299 
******************************************************************
******************************************************************
Family:  c("BEZI", "Zero Inflated Beta") 

Call:  gamlss::gamlss(formula = stats::as.formula(paste(pathname[i],  
    paste(c(comvar, adjustvar), collapse = "+"), sep = "~")),  
    family = BEZI, data = testdat, trace = FALSE) 

Fitting method: RS() 

------------------------------------------------------------------
Mu link function:  logit
Mu Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        2.246380   0.019984 112.410  < 2e-16 ***
genderMale        -0.004272   0.014308  -0.299    0.765    
age.sample        -0.008404   0.001235  -6.804 1.76e-11 ***
bfNon_exclusiveBF  0.022727   0.023196   0.980    0.327    
bfNo_BF           -0.042615   0.038138  -1.117    0.264    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
Sigma link function:  log
Sigma Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.36087    0.04449   120.5   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
Nu link function:  logit 
Nu Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   -21.43     866.54  -0.025     0.98

------------------------------------------------------------------
No. of observations in the fit:  995 
Degrees of Freedom for the fit:  7
      Residual Deg. of Freedom:  988 
                      at cycle:  7 
 
Global Deviance:     -4914.618 
            AIC:     -4900.618 
            SBC:     -4866.299 
******************************************************************
Warning messages:
1: In summary.gamlss(gamlss::gamlss(stats::as.formula(paste(pathname[i],  :
  summary: vcov has failed, option qr is used instead

2: In summary.gamlss(gamlss::gamlss(stats::as.formula(paste(pathname[i],  :
  summary: vcov has failed, option qr is used instead

                                    id Estimate.genderMale    ll   ul
1                   Cellular.Processes                   0 -0.02 0.03
2 Environmental.Information.Processing                   0 -0.03 0.02
  Pr(>|t|).genderMale pval.adjust.genderMale
1              0.7654                 0.7654
2              0.7654                 0.7654

metamicrobiomeR documentation built on Nov. 9, 2020, 5:06 p.m.