taxa.compare: Compare taxa relative abundance

Description Usage Arguments Value Examples

View source: R/taxa.compare.R

Description

This function compares taxa relative abundance summary tables at all levels between groups using GAMLSS with BEZI or Linear/Linear Mixed Effect models (LM/LMEM) after filtering (using prevalence and relative abundance thresholds).

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
taxa.compare(
  taxtab,
  propmed.rel = "gamlss",
  transform = "none",
  zeroreplace.method = "none",
  comvar,
  adjustvar,
  personid = "personid",
  longitudinal = "yes",
  percent.filter = 0.05,
  relabund.filter = 5e-05,
  p.adjust.method = "fdr"
)

Arguments

taxtab

taxa relative abundance table (already merged to mapping file) from phylum to species or any preferred highest taxa level.

propmed.rel

statistical method for comparing relative abundance. Options are "lm" for LM/LMEM or "gamlss" for GAMLSS with BEZI family.

transform

transformation of relative abundance data. Options are "none" for no transformation, "gmpr" for Geometric Mean of Pairwise Ratios (GMPR) normalization, "asin.sqrt" for arcsine transformation, "logit" for logit transformation, "clr" for centered log ratio transformation. Default is "none".

zeroreplace.method

Method for zero replacement implemented in R package *zCompositions*. Options are "none" for no replacement, "multKM" for Multiplicative Kaplan-Meier smoothing spline (KMSS) replacement, "multLN" for Multiplicative lognormal replacement, "multRepl" for Multiplicative simple replacement, "lrEM" for Log-ratio EM algorithm, "lrDA" for Log-ratio DA algorithm. Default is "none".

comvar

main variable for comparison

adjustvar

variables to be adjusted.

personid

name of variable for person id (applicable for longitudinal data)

longitudinal

whether data is longitudinal? Options are "yes" or "no". Default is "yes".

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.

p.adjust.method

method for multiple testing adjustment. Options are those of the p.adjust.methods of stats:: p.adjust function. Default for this function is "fdr".

Value

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

Examples

1
2
3
4
5
6
7
#Load summary tables of bacterial taxa relative abundance from Bangladesh data
data(taxtab6)
tab6<-as.data.frame(taxtab6)
tl<-colnames(taxtab6)[grep("k__bacteria.p__fusobacteria",colnames(taxtab6))]
taxacom.ex<-taxa.compare(taxtab=tab6[,c("personid","x.sampleid","bf","age.sample",tl)],
propmed.rel="gamlss",comvar="bf",adjustvar="age.sample",
longitudinal="yes",p.adjust.method="fdr")

Example output

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

Attaching package:gamlss.dataThe following object is masked frompackage: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.

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

Call:  gamlss::gamlss(formula = stats::as.formula(paste(taxname[i],  
    paste(c(comvar, adjustvar, "random(personid)"), 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)       -5.77401    0.30668 -18.827   <2e-16 ***
bfNon_exclusiveBF -0.05981    0.36399  -0.164    0.870    
bfNo_BF           -0.43882    1.06824  -0.411    0.682    
age.sample         0.02862    0.08537   0.335    0.738    
---
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)   4.6285     0.1604   28.85   <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)   1.7177     0.1551   11.07   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
NOTE: Additive smoothing terms exist in the formulas: 
 i) Std. Error for smoothers are for the linear effect only. 
ii) Std. Error for the linear terms may not be reliable. 
------------------------------------------------------------------
No. of observations in the fit:  322 
Degrees of Freedom for the fit:  41.46391
      Residual Deg. of Freedom:  280.5361 
                      at cycle:  20 
 
Global Deviance:     -377.117 
            AIC:     -294.1892 
            SBC:     -137.6816 
******************************************************************
******************************************************************
Family:  c("BEZI", "Zero Inflated Beta") 

Call:  gamlss::gamlss(formula = stats::as.formula(paste(taxname[i],  
    paste(c(comvar, adjustvar, "random(personid)"), 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)       -5.77401    0.30668 -18.827   <2e-16 ***
bfNon_exclusiveBF -0.05981    0.36399  -0.164    0.870    
bfNo_BF           -0.43882    1.06824  -0.411    0.682    
age.sample         0.02862    0.08537   0.335    0.738    
---
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)   4.6285     0.1604   28.85   <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)   1.7177     0.1551   11.07   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
NOTE: Additive smoothing terms exist in the formulas: 
 i) Std. Error for smoothers are for the linear effect only. 
ii) Std. Error for the linear terms may not be reliable. 
------------------------------------------------------------------
No. of observations in the fit:  322 
Degrees of Freedom for the fit:  41.46391
      Residual Deg. of Freedom:  280.5361 
                      at cycle:  20 
 
Global Deviance:     -377.117 
            AIC:     -294.1892 
            SBC:     -137.6816 
******************************************************************
******************************************************************
Family:  c("BEZI", "Zero Inflated Beta") 

Call:  gamlss::gamlss(formula = stats::as.formula(paste(taxname[i],  
    paste(c(comvar, adjustvar, "random(personid)"), 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)       -5.77401    0.30668 -18.827   <2e-16 ***
bfNon_exclusiveBF -0.05981    0.36399  -0.164    0.870    
bfNo_BF           -0.43882    1.06824  -0.411    0.682    
age.sample         0.02862    0.08537   0.335    0.738    
---
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)   4.6285     0.1604   28.85   <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)   1.7177     0.1551   11.07   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
NOTE: Additive smoothing terms exist in the formulas: 
 i) Std. Error for smoothers are for the linear effect only. 
ii) Std. Error for the linear terms may not be reliable. 
------------------------------------------------------------------
No. of observations in the fit:  322 
Degrees of Freedom for the fit:  41.46391
      Residual Deg. of Freedom:  280.5361 
                      at cycle:  20 
 
Global Deviance:     -377.117 
            AIC:     -294.1892 
            SBC:     -137.6816 
******************************************************************
******************************************************************
Family:  c("BEZI", "Zero Inflated Beta") 

Call:  gamlss::gamlss(formula = stats::as.formula(paste(taxname[i],  
    paste(c(comvar, adjustvar, "random(personid)"), 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)       -5.78117    0.32065 -18.029   <2e-16 ***
bfNon_exclusiveBF -0.06259    0.36574  -0.171    0.864    
bfNo_BF           -0.44786    1.06816  -0.419    0.675    
age.sample         0.02855    0.08629   0.331    0.741    
---
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)   4.6510     0.1622   28.67   <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)   1.7419     0.1565   11.13   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
NOTE: Additive smoothing terms exist in the formulas: 
 i) Std. Error for smoothers are for the linear effect only. 
ii) Std. Error for the linear terms may not be reliable. 
------------------------------------------------------------------
No. of observations in the fit:  322 
Degrees of Freedom for the fit:  40.49761
      Residual Deg. of Freedom:  281.5024 
                      at cycle:  20 
 
Global Deviance:     -366.568 
            AIC:     -285.5728 
            SBC:     -132.7125 
******************************************************************
******************************************************************
Family:  c("BEZI", "Zero Inflated Beta") 

Call:  gamlss::gamlss(formula = stats::as.formula(paste(taxname[i],  
    paste(c(comvar, adjustvar, "random(personid)"), 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)       -5.76889    0.37948 -15.202   <2e-16 ***
bfNon_exclusiveBF -0.11675    0.42666  -0.274    0.785    
bfNo_BF           -0.51402    1.09430  -0.470    0.639    
age.sample         0.03977    0.09919   0.401    0.689    
---
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)   4.6293     0.1898   24.39   <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)    2.104      0.179   11.75   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
NOTE: Additive smoothing terms exist in the formulas: 
 i) Std. Error for smoothers are for the linear effect only. 
ii) Std. Error for the linear terms may not be reliable. 
------------------------------------------------------------------
No. of observations in the fit:  322 
Degrees of Freedom for the fit:  32.7031
      Residual Deg. of Freedom:  289.2969 
                      at cycle:  20 
 
Global Deviance:     -244.1 
            AIC:     -178.6938 
            SBC:     -55.2543 
******************************************************************
Warning messages:
1: In RS() : Algorithm RS has not yet converged
2: In summary.gamlss(gamlss::gamlss(stats::as.formula(paste(taxname[i],  :
  summary: vcov has failed, option qr is used instead

3: In RS() : Algorithm RS has not yet converged
4: In summary.gamlss(gamlss::gamlss(stats::as.formula(paste(taxname[i],  :
  summary: vcov has failed, option qr is used instead

5: In RS() : Algorithm RS has not yet converged
6: In summary.gamlss(gamlss::gamlss(stats::as.formula(paste(taxname[i],  :
  summary: vcov has failed, option qr is used instead

7: In RS() : Algorithm RS has not yet converged
8: In summary.gamlss(gamlss::gamlss(stats::as.formula(paste(taxname[i],  :
  summary: vcov has failed, option qr is used instead

9: In RS() : Algorithm RS has not yet converged
10: In summary.gamlss(gamlss::gamlss(stats::as.formula(paste(taxname[i],  :
  summary: vcov has failed, option qr is used instead

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