Analysis objectives

  1. Import, recode, and subset data from bugsigdb.org
  2. Create a table of studies
  3. Calculate the frequency of appearance of each taxa in independent signatures and identify the most frequently reported taxa
  4. Estimate the probability of the most frequently identified taxa occuring by chance

Making sure packages are installed

install.packages(c("devtools", "tidyverse", "kableExtra"))
devtools::install_github("waldronlab/bugSigSimple")
devtools::install_github("waldronlab/BugSigDBStats")
devtools::install_github("waldronlab/bugsigdbr")

Load and subset data

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
##  ggplot2 3.3.5      purrr   0.3.4
##  tibble  3.1.6      dplyr   1.0.8
##  tidyr   1.2.0      stringr 1.4.0
##  readr   2.1.2      forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
##  dplyr::filter() masks stats::filter()
##  dplyr::lag()    masks stats::lag()
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
dat <- bugsigdbr::importBugSigDB(cache= TRUE)
## Using cached version from 2022-03-31 12:38:44
dim(dat)
## [1] 2098   48
names(dat)
##  [1] "Study"                      "Study design"              
##  [3] "PMID"                       "DOI"                       
##  [5] "URL"                        "Authors"                   
##  [7] "Title"                      "Journal"                   
##  [9] "Year"                       "Experiment"                
## [11] "Location of subjects"       "Host species"              
## [13] "Body site"                  "UBERON ID"                 
## [15] "Condition"                  "EFO ID"                    
## [17] "Group 0 name"               "Group 1 name"              
## [19] "Group 1 definition"         "Group 0 sample size"       
## [21] "Group 1 sample size"        "Antibiotics exclusion"     
## [23] "Sequencing type"            "16S variable region"       
## [25] "Sequencing platform"        "Statistical test"          
## [27] "Significance threshold"     "MHT correction"            
## [29] "LDA Score above"            "Matched on"                
## [31] "Confounders controlled for" "Pielou"                    
## [33] "Shannon"                    "Chao1"                     
## [35] "Simpson"                    "Inverse Simpson"           
## [37] "Richness"                   "Signature page name"       
## [39] "Source"                     "Curated date"              
## [41] "Curator"                    "Revision editor"           
## [43] "Description"                "Abundance in Group 1"      
## [45] "MetaPhlAn taxon names"      "NCBI Taxonomy IDs"         
## [47] "State"                      "Reviewer"

Subsetting

Subsetting only studies done on humans - 6 of 17 studies were excluded because they were either done on mice or not statistically reliable (Cregger et. al)

subset.dat <-
  dat[which(
    dat$PMID == "30778155" |
      dat$PMID == "32192080" |
      dat$PMID == "31087436" |
      dat$PMID == "26901400" |
      dat$PMID == "33839907" |
      dat$PMID == "32046455" |
      dat$PMID == "33925708" |
      dat$PMID == "32299442" |
      dat$PMID == "33313185" |
      dat$PMID == "34268384" | dat$PMID == "33660232"
  ), ]

All studies

Summary of studies and most frequent taxa increased and decreased in endometriosis patients for all studies

bugSigSimple::createStudyTable(subset.dat)
## # A tibble: 11 × 5
##    Study          Condition     Cases Controls `Study Design`                   
##    <chr>          <chr>         <dbl>    <dbl> <chr>                            
##  1 Akiyama 2019   endometriosis    30       39 cross-sectional observational, n…
##  2 Ata 2019       endometriosis    14       14 prospective cohort               
##  3 Chao 2021      endometriosis    37       66 case-control                     
##  4 Chen 2020      endometriosis    14       67 cross-sectional observational, n…
##  5 Hernandes 2020 endometriosis    10       11 case-control                     
##  6 Khan 2016      endometriosis    32       32 case-control                     
##  7 Lee 2021       endometriosis    45       45 case-control                     
##  8 Perrotta 2020  endometriosis    14       21 cross-sectional observational, n…
##  9 Shan 2021      endometriosis    12       12 case-control                     
## 10 Svensson 2021  endometriosis    66      198 case-control                     
## 11 Wei 2020       endometriosis    36       14 cross-sectional observational, n…
getMostFrequentTaxa(subset.dat,n=30)
## 
##                                      k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Atopobiaceae|g__Atopobium 
##                                                                                                                                     7 
##                                          k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Streptococcaceae|g__Streptococcus 
##                                                                                                                                     5 
##                                        k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria|o__Enterobacterales|f__Enterobacteriaceae 
##                                                                                                                                     5 
##                                              k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Coriobacteriaceae 
##                                                                                                                                     4 
##                                                           k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Streptococcaceae 
##                                                                                                                                     4 
##                    k__Bacteria|p__Proteobacteria|c__Epsilonproteobacteria|o__Campylobacterales|f__Campylobacteraceae|g__Campylobacter 
##                                                                                                                                     4 
##            k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria|o__Enterobacterales|f__Enterobacteriaceae|s__Escherichia/Shigella sp. 
##                                                                                                                                     4 
## k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria|o__Pseudomonadales|f__Pseudomonadaceae|g__Pseudomonas|s__Pseudomonas viridiflava 
##                                                                                                                                     4 
##                              k__Bacteria|p__Actinobacteria|c__Actinomycetia|o__Bifidobacteriales|f__Bifidobacteriaceae|g__Gardnerella 
##                                                                                                                                     3 
##                   k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Lactobacillaceae|g__Lactobacillus|s__Lactobacillus iners 
##                                                                                                                                     3 
##                                             k__Bacteria|p__Firmicutes|c__Clostridia|o__Eubacteriales|f__Clostridiaceae|g__Clostridium 
##                                                                                                                                     3 
##                                            k__Bacteria|p__Firmicutes|c__Clostridia|o__Eubacteriales|f__Lachnospiraceae|g__Lachnospira 
##                                                                                                                                     3 
##                                           k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria|o__Pseudomonadales|f__Pseudomonadaceae 
##                                                                                                                                     3 
##                                                                                                         k__Bacteria|p__Actinobacteria 
##                                                                                                                                     2 
##                                  k__Bacteria|p__Actinobacteria|c__Actinomycetia|o__Actinomycetales|f__Actinomycetaceae|g__Actinomyces 
##                                                                                                                                     2 
##                                   k__Bacteria|p__Actinobacteria|c__Actinomycetia|o__Actinomycetales|f__Actinomycetaceae|g__Mobiluncus 
##                                                                                                                                     2 
##                          k__Bacteria|p__Actinobacteria|c__Actinomycetia|o__Bifidobacteriales|f__Bifidobacteriaceae|g__Bifidobacterium 
##                                                                                                                                     2 
##                                                     k__Bacteria|p__Actinobacteria|c__Actinomycetia|o__Micrococcales|f__Micrococcaceae 
##                                                                                                                                     2 
##                                                                          k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales 
##                                                                                                                                     2 
##                                          k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales|f__Prevotellaceae|g__Prevotella 
##                                                                                                                                     2 
##                                                               k__Bacteria|p__Firmicutes|c__Bacilli|o__Bacillales|f__Staphylococcaceae 
##                                                                                                                                     2 
##                                                                         k__Bacteria|p__Firmicutes|c__Bacilli|o__Bacillales|g__Gemella 
##                                                                                                                                     2 
##                                                k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Aerococcaceae|g__Aerococcus 
##                                                                                                                                     2 
##                                            k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Enterococcaceae|g__Enterococcus 
##                                                                                                                                     2 
##                                                           k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Lactobacillaceae 
##                                                                                                                                     2 
##                                          k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Lactobacillaceae|g__Lactobacillus 
##                                                                                                                                     2 
##            k__Bacteria|p__Firmicutes|c__Clostridia|o__Eubacteriales|f__Lachnospiraceae|g__Mediterraneibacter|s__[Ruminococcus] gnavus 
##                                                                                                                                     2 
##                                      k__Bacteria|p__Firmicutes|c__Clostridia|o__Eubacteriales|f__Oscillospiraceae|g__Faecalibacterium 
##                                                                                                                                     2 
##                                          k__Bacteria|p__Firmicutes|c__Clostridia|o__Eubacteriales|f__Oscillospiraceae|g__Oscillospira 
##                                                                                                                                     2 
##                                     k__Bacteria|p__Firmicutes|c__Clostridia|o__Eubacteriales|f__Oscillospiraceae|g__Ruminiclostridium 
##                                                                                                                                     2
getMostFrequentTaxa(subset.dat,sig.type="increased")
## 
##                                          k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Streptococcaceae|g__Streptococcus 
##                                                                                                                                     5 
##                                        k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria|o__Enterobacterales|f__Enterobacteriaceae 
##                                                                                                                                     5 
##                                                           k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Streptococcaceae 
##                                                                                                                                     4 
##            k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria|o__Enterobacterales|f__Enterobacteriaceae|s__Escherichia/Shigella sp. 
##                                                                                                                                     4 
##                                              k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Coriobacteriaceae 
##                                                                                                                                     3 
##                                             k__Bacteria|p__Firmicutes|c__Clostridia|o__Eubacteriales|f__Clostridiaceae|g__Clostridium 
##                                                                                                                                     3 
##                                           k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria|o__Pseudomonadales|f__Pseudomonadaceae 
##                                                                                                                                     3 
## k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria|o__Pseudomonadales|f__Pseudomonadaceae|g__Pseudomonas|s__Pseudomonas viridiflava 
##                                                                                                                                     3 
##                          k__Bacteria|p__Actinobacteria|c__Actinomycetia|o__Bifidobacteriales|f__Bifidobacteriaceae|g__Bifidobacterium 
##                                                                                                                                     2 
##                              k__Bacteria|p__Actinobacteria|c__Actinomycetia|o__Bifidobacteriales|f__Bifidobacteriaceae|g__Gardnerella 
##                                                                                                                                     2
getMostFrequentTaxa(subset.dat,sig.type="decreased")
## 
##                    k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Atopobiaceae|g__Atopobium 
##                                                                                                                   5 
##                 k__Bacteria|p__Actinobacteria|c__Actinomycetia|o__Actinomycetales|f__Actinomycetaceae|g__Mobiluncus 
##                                                                                                                   2 
##                                                        k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales 
##                                                                                                                   2 
##                        k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales|f__Prevotellaceae|g__Prevotella 
##                                                                                                                   2 
## k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Lactobacillaceae|g__Lactobacillus|s__Lactobacillus iners 
##                                                                                                                   2 
##                          k__Bacteria|p__Firmicutes|c__Clostridia|o__Eubacteriales|f__Lachnospiraceae|g__Lachnospira 
##                                                                                                                   2 
##                      k__Bacteria|p__Firmicutes|c__Negativicutes|o__Veillonellales|f__Veillonellaceae|g__Megasphaera 
##                                                                                                                   2 
##                     k__Bacteria|p__Fusobacteria|c__Fusobacteriia|o__Fusobacteriales|f__Leptotrichiaceae|g__Sneathia 
##                                                                                                                   2 
##  k__Bacteria|p__Proteobacteria|c__Epsilonproteobacteria|o__Campylobacterales|f__Campylobacteraceae|g__Campylobacter 
##                                                                                                                   2 
##                                                                                       k__Bacteria|p__Actinobacteria 
##                                                                                                                   1

Excluding feces samples

Summary of studies and most frequent taxa in only samples from female reproductive tract, excluding feces samples

subset.dat2 <-
  dat[which(
    dat$PMID == "30778155" |
      dat$PMID == "32192080" |
      dat$PMID == "31087436" |
      dat$PMID == "26901400" |
      dat$PMID == "32046455" |
      dat$PMID == "33925708" |
      dat$PMID == "32299442" |
      dat$PMID == "33313185" | dat$PMID == "34268384"
  ), ]

reproductive_sigs <-
  subset.dat2[which(subset.dat2$`Body site` != "feces" |
                      is.na(subset.dat2$`Body site`)), ]

bugSigSimple::createStudyTable(reproductive_sigs)
## # A tibble: 9 × 5
##   Study          Condition     Cases Controls `Study Design`                    
##   <chr>          <chr>         <dbl>    <dbl> <chr>                             
## 1 Akiyama 2019   endometriosis    30       39 cross-sectional observational, no…
## 2 Ata 2019       endometriosis    14       14 prospective cohort                
## 3 Chao 2021      endometriosis    37       66 case-control                      
## 4 Chen 2020      endometriosis    14       67 cross-sectional observational, no…
## 5 Hernandes 2020 endometriosis    10       11 case-control                      
## 6 Khan 2016      endometriosis    32       32 case-control                      
## 7 Lee 2021       endometriosis    45       45 case-control                      
## 8 Perrotta 2020  endometriosis    14       21 cross-sectional observational, no…
## 9 Wei 2020       endometriosis    36       14 cross-sectional observational, no…
allfreqs <- bugSigSimple::createTaxonTable(reproductive_sigs, n = 20) %>% #could change number
  arrange(I(decreased_signatures - increased_signatures))
incfreqs <- filter(allfreqs, I(increased_signatures - decreased_signatures) > 0)
decfreqs <- filter(allfreqs, I(increased_signatures - decreased_signatures) < 0)
kableExtra::kbl(allfreqs) %>%
  kable_paper("hover", full_width = FALSE)
Taxon Name Taxonomic Level total_signatures increased_signatures decreased_signatures Binomial Test pval kingdom phylum class order family genus species n_signatures metaphlan_name
Streptococcus genus 4 4 0 0.12 Bacteria Firmicutes Bacilli Lactobacillales Streptococcaceae Streptococcus NA 4 k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Streptococcaceae|g__Streptococcus
Escherichia/Shigella sp. species 4 4 0 0.12 Bacteria Proteobacteria Gammaproteobacteria Enterobacterales Enterobacteriaceae Escherichia/Shigella sp. NA 4 k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria|o__Enterobacterales|f__Enterobacteriaceae|s__Escherichia/Shigella sp.
Coriobacteriaceae family 3 3 0 0.25 Bacteria Actinobacteria Coriobacteriia Coriobacteriales Coriobacteriaceae NA NA 3 k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Coriobacteriaceae
Streptococcaceae family 3 3 0 0.25 Bacteria Firmicutes Bacilli Lactobacillales Streptococcaceae NA NA 7 k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Streptococcaceae
Enterobacteriaceae family 3 3 0 0.25 Bacteria Proteobacteria Gammaproteobacteria Enterobacterales Enterobacteriaceae NA NA 7 k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria|o__Enterobacterales|f__Enterobacteriaceae
Pseudomonadaceae family 3 3 0 0.25 Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae NA NA 7 k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria|o__Pseudomonadales|f__Pseudomonadaceae
Pseudomonas viridiflava species 4 3 1 0.63 Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas Pseudomonas viridiflava 4 k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria|o__Pseudomonadales|f__Pseudomonadaceae|g__Pseudomonas|s__Pseudomonas viridiflava
Gardnerella genus 2 2 0 0.50 Bacteria Actinobacteria Actinomycetia Bifidobacteriales Bifidobacteriaceae Gardnerella NA 2 k__Bacteria|p__Actinobacteria|c__Actinomycetia|o__Bifidobacteriales|f__Bifidobacteriaceae|g__Gardnerella
Staphylococcaceae family 2 2 0 0.50 Bacteria Firmicutes Bacilli Bacillales Staphylococcaceae NA NA 4 k__Bacteria|p__Firmicutes|c__Bacilli|o__Bacillales|f__Staphylococcaceae
Clostridium genus 2 2 0 0.50 Bacteria Firmicutes Clostridia Eubacteriales Clostridiaceae Clostridium NA 2 k__Bacteria|p__Firmicutes|c__Clostridia|o__Eubacteriales|f__Clostridiaceae|g__Clostridium
Faecalibacterium genus 2 2 0 0.50 Bacteria Firmicutes Clostridia Eubacteriales Oscillospiraceae Faecalibacterium NA 2 k__Bacteria|p__Firmicutes|c__Clostridia|o__Eubacteriales|f__Oscillospiraceae|g__Faecalibacterium
Campylobacter genus 4 2 2 1.00 Bacteria Proteobacteria Epsilonproteobacteria Campylobacterales Campylobacteraceae Campylobacter NA 4 k__Bacteria|p__Proteobacteria|c__Epsilonproteobacteria|o__Campylobacterales|f__Campylobacteraceae|g__Campylobacter
Micrococcaceae family 2 1 1 1.00 Bacteria Actinobacteria Actinomycetia Micrococcales Micrococcaceae NA NA 3 k__Bacteria|p__Actinobacteria|c__Actinomycetia|o__Micrococcales|f__Micrococcaceae
Gemella genus 2 1 1 1.00 Bacteria Firmicutes Bacilli Bacillales Gemella NA NA 2 k__Bacteria|p__Firmicutes|c__Bacilli|o__Bacillales|g__Gemella
Aerococcus genus 2 1 1 1.00 Bacteria Firmicutes Bacilli Lactobacillales Aerococcaceae Aerococcus NA 2 k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Aerococcaceae|g__Aerococcus
Erysipelotrichaceae family 2 1 1 1.00 Bacteria Firmicutes Erysipelotrichia Erysipelotrichales Erysipelotrichaceae NA NA 2 k__Bacteria|p__Firmicutes|c__Erysipelotrichia|o__Erysipelotrichales|f__Erysipelotrichaceae
Lactobacillus iners species 3 1 2 1.00 Bacteria Firmicutes Bacilli Lactobacillales Lactobacillaceae Lactobacillus Lactobacillus iners 3 k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Lactobacillaceae|g__Lactobacillus|s__Lactobacillus iners
Mobiluncus genus 2 0 2 0.50 Bacteria Actinobacteria Actinomycetia Actinomycetales Actinomycetaceae Mobiluncus NA 2 k__Bacteria|p__Actinobacteria|c__Actinomycetia|o__Actinomycetales|f__Actinomycetaceae|g__Mobiluncus
Prevotella genus 2 0 2 0.50 Bacteria Bacteroidetes Bacteroidia Bacteroidales Prevotellaceae Prevotella NA 3 k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales|f__Prevotellaceae|g__Prevotella
Atopobium genus 7 2 5 0.45 Bacteria Actinobacteria Coriobacteriia Coriobacteriales Atopobiaceae Atopobium NA 7 k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Atopobiaceae|g__Atopobium

Load HMP data from curatedMetagenomicData

suppressPackageStartupMessages(library(curatedMetagenomicData))
# hmp_vagina_metadata <- filter(sampleMetadata, study_name == "HMP_2012" & body_site == "vagina")
all_healthy_adult_vagina_metadata <- filter(sampleMetadata, body_site == "vagina" & disease == "healthy" & age_category == "adult")
se <-
  curatedMetagenomicData::returnSamples(all_healthy_adult_vagina_metadata, dataType = "relative_abundance", counts = FALSE)

Get matrices of species and genus relative abundance

allranks <- mia::splitByRanks(se)
#species_relab <- t(assay(allranks[["species"]]))
genus_relab <- t(assay(allranks[["genus"]]))
#family_relab <- t(assay(allranks[["family"]]))

Centered log-ratio transformation on genus level relative abundance:

genus_relab_clr <- Hotelling::clr(genus_relab + 1)
df <- dplyr::select(data.frame(genus_relab_clr), c("Streptococcus", "Lactobacillus"))
df$age <- se$age
df$BMI <- se$BMI
fit <- lm(Streptococcus ~ Lactobacillus + age + BMI, data = df)  
summary(fit)
## 
## Call:
## lm(formula = Streptococcus ~ Lactobacillus + age + BMI, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6728 -0.1875 -0.0264  0.0893  3.6729 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.386413   0.436694  -0.885 0.378987    
## Lactobacillus -0.078095   0.047933  -1.629 0.107344    
## age            0.035738   0.010073   3.548 0.000665 ***
## BMI           -0.009094   0.014295  -0.636 0.526581    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4857 on 77 degrees of freedom
##   (14 observations deleted due to missingness)
## Multiple R-squared:  0.1799, Adjusted R-squared:  0.1479 
## F-statistic: 5.629 on 3 and 77 DF,  p-value: 0.00153
summary(lm(Streptococcus ~ Lactobacillus, data = df))
## 
## Call:
## lm(formula = Streptococcus ~ Lactobacillus, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6006 -0.0609 -0.0462 -0.0440  4.1909 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    0.55246    0.20193   2.736  0.00745 **
## Lactobacillus -0.11409    0.04751  -2.401  0.01832 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5642 on 93 degrees of freedom
## Multiple R-squared:  0.05839,    Adjusted R-squared:  0.04826 
## F-statistic: 5.767 on 1 and 93 DF,  p-value: 0.01832
wilcox.test(df$Streptococcus, df$Lactobacillus)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  df$Streptococcus and df$Lactobacillus
## W = 213.5, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
t.test(df$Streptococcus, df$Lactobacillus)
## 
##  Welch Two Sample t-test
## 
## data:  df$Streptococcus and df$Lactobacillus
## t = -28.666, df = 133.93, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -4.258800 -3.709054
## sample estimates:
## mean of x mean of y 
## 0.0879008 4.0718280
secount <-
  curatedMetagenomicData::returnSamples(all_healthy_adult_vagina_metadata, dataType = "relative_abundance", counts = TRUE)
## snapshotDate(): 2021-10-19
## 
## $`2021-03-31.FerrettiP_2018.relative_abundance`
## dropping rows without rowTree matches:
##   k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Coriobacteriaceae|g__Collinsella|s__Collinsella_stercoris
##   k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Coriobacteriaceae|g__Enorma|s__[Collinsella]_massiliensis
##   k__Bacteria|p__Firmicutes|c__Bacilli|o__Bacillales|f__Bacillales_unclassified|g__Gemella|s__Gemella_bergeri
##   k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Carnobacteriaceae|g__Granulicatella|s__Granulicatella_elegans
##   k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Ruminococcaceae|g__Ruminococcus|s__Ruminococcus_champanellensis
##   k__Bacteria|p__Firmicutes|c__Erysipelotrichia|o__Erysipelotrichales|f__Erysipelotrichaceae|g__Bulleidia|s__Bulleidia_extructa
##   k__Bacteria|p__Proteobacteria|c__Betaproteobacteria|o__Burkholderiales|f__Sutterellaceae|g__Sutterella|s__Sutterella_parvirubra
##   k__Bacteria|p__Spirochaetes|c__Spirochaetia|o__Leptospirales|f__Leptospiraceae|g__Leptospira|s__Leptospira_hartskeerlii
##   k__Bacteria|p__Synergistetes|c__Synergistia|o__Synergistales|f__Synergistaceae|g__Cloacibacillus|s__Cloacibacillus_evryensis
## $`2021-03-31.HMP_2012.relative_abundance`
## dropping rows without rowTree matches:
##   k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Atopobiaceae|g__Olsenella|s__Olsenella_profusa
##   k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Coriobacteriaceae|g__Collinsella|s__Collinsella_stercoris
##   k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Coriobacteriaceae|g__Enorma|s__[Collinsella]_massiliensis
##   k__Bacteria|p__Firmicutes|c__Bacilli|o__Bacillales|f__Bacillales_unclassified|g__Gemella|s__Gemella_bergeri
##   k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Carnobacteriaceae|g__Granulicatella|s__Granulicatella_elegans
##   k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Ruminococcaceae|g__Ruminococcus|s__Ruminococcus_champanellensis
##   k__Bacteria|p__Firmicutes|c__Erysipelotrichia|o__Erysipelotrichales|f__Erysipelotrichaceae|g__Bulleidia|s__Bulleidia_extructa
##   k__Bacteria|p__Proteobacteria|c__Betaproteobacteria|o__Burkholderiales|f__Sutterellaceae|g__Sutterella|s__Sutterella_parvirubra
##   k__Bacteria|p__Synergistetes|c__Synergistia|o__Synergistales|f__Synergistaceae|g__Cloacibacillus|s__Cloacibacillus_evryensis
## $`2021-03-31.WampachL_2018.relative_abundance`
## dropping rows without rowTree matches:
##   k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Coriobacteriaceae|g__Collinsella|s__Collinsella_stercoris
##   k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Coriobacteriaceae|g__Enorma|s__[Collinsella]_massiliensis
##   k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Carnobacteriaceae|g__Granulicatella|s__Granulicatella_elegans
##   k__Bacteria|p__Proteobacteria|c__Betaproteobacteria|o__Burkholderiales|f__Sutterellaceae|g__Sutterella|s__Sutterella_parvirubra
##   k__Bacteria|p__Synergistetes|c__Synergistia|o__Synergistales|f__Synergistaceae|g__Cloacibacillus|s__Cloacibacillus_evryensis
allranks <- mia::splitByRanks(secount)
genus_relab_count <- t(assay(allranks[["genus"]]))
df <- dplyr::select(data.frame(genus_relab_count), c("Streptococcus", "Lactobacillus"))
df$age <- secount$age
df$BMI <- secount$BMI
## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileis
fit.ZInegbin <- pscl::zeroinfl(Streptococcus ~ Lactobacillus + age + BMI | 1,
                data = df,
                dist = "negbin")

Create a data frame for regression analysis from healthy controls for log-linear model

library(curatedMetagenomicData)
library(dplyr)
library(DT)
library(purrr)

taxon_strep <- c("1884" = "Streptococcus")
taxon_lact <- c("1578" = "Lactobacillus")

target_sample <- "vagina"
sample_metadata <- sampleMetadata %>% 
    filter(
        grepl(target_sample, body_site) | grepl(target_sample, body_subsite)
    ) %>% 
    discard(~all(is.na(.x))) %>% 
    as_tibble()
dim(sample_metadata)
## [1] 96 26
unique(sample_metadata$study_name)
## [1] "FerrettiP_2018" "HMP_2012"       "WampachL_2018"
tse <- returnSamples(
    sampleMetadata = sample_metadata,
    dataType = "relative_abundance",
    counts = FALSE,
    rownames = "short"
)
## snapshotDate(): 2021-10-19
## 
## $`2021-03-31.FerrettiP_2018.relative_abundance`
## dropping rows without rowTree matches:
##   k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Coriobacteriaceae|g__Collinsella|s__Collinsella_stercoris
##   k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Coriobacteriaceae|g__Enorma|s__[Collinsella]_massiliensis
##   k__Bacteria|p__Firmicutes|c__Bacilli|o__Bacillales|f__Bacillales_unclassified|g__Gemella|s__Gemella_bergeri
##   k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Carnobacteriaceae|g__Granulicatella|s__Granulicatella_elegans
##   k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Ruminococcaceae|g__Ruminococcus|s__Ruminococcus_champanellensis
##   k__Bacteria|p__Firmicutes|c__Erysipelotrichia|o__Erysipelotrichales|f__Erysipelotrichaceae|g__Bulleidia|s__Bulleidia_extructa
##   k__Bacteria|p__Proteobacteria|c__Betaproteobacteria|o__Burkholderiales|f__Sutterellaceae|g__Sutterella|s__Sutterella_parvirubra
##   k__Bacteria|p__Spirochaetes|c__Spirochaetia|o__Leptospirales|f__Leptospiraceae|g__Leptospira|s__Leptospira_hartskeerlii
##   k__Bacteria|p__Synergistetes|c__Synergistia|o__Synergistales|f__Synergistaceae|g__Cloacibacillus|s__Cloacibacillus_evryensis
## $`2021-03-31.HMP_2012.relative_abundance`
## dropping rows without rowTree matches:
##   k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Atopobiaceae|g__Olsenella|s__Olsenella_profusa
##   k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Coriobacteriaceae|g__Collinsella|s__Collinsella_stercoris
##   k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Coriobacteriaceae|g__Enorma|s__[Collinsella]_massiliensis
##   k__Bacteria|p__Firmicutes|c__Bacilli|o__Bacillales|f__Bacillales_unclassified|g__Gemella|s__Gemella_bergeri
##   k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Carnobacteriaceae|g__Granulicatella|s__Granulicatella_elegans
##   k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Ruminococcaceae|g__Ruminococcus|s__Ruminococcus_champanellensis
##   k__Bacteria|p__Firmicutes|c__Erysipelotrichia|o__Erysipelotrichales|f__Erysipelotrichaceae|g__Bulleidia|s__Bulleidia_extructa
##   k__Bacteria|p__Proteobacteria|c__Betaproteobacteria|o__Burkholderiales|f__Sutterellaceae|g__Sutterella|s__Sutterella_parvirubra
##   k__Bacteria|p__Synergistetes|c__Synergistia|o__Synergistales|f__Synergistaceae|g__Cloacibacillus|s__Cloacibacillus_evryensis
## $`2021-03-31.WampachL_2018.relative_abundance`
## dropping rows without rowTree matches:
##   k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Coriobacteriaceae|g__Collinsella|s__Collinsella_stercoris
##   k__Bacteria|p__Actinobacteria|c__Coriobacteriia|o__Coriobacteriales|f__Coriobacteriaceae|g__Enorma|s__[Collinsella]_massiliensis
##   k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Carnobacteriaceae|g__Granulicatella|s__Granulicatella_elegans
##   k__Bacteria|p__Proteobacteria|c__Betaproteobacteria|o__Burkholderiales|f__Sutterellaceae|g__Sutterella|s__Sutterella_parvirubra
##   k__Bacteria|p__Synergistetes|c__Synergistia|o__Synergistales|f__Synergistaceae|g__Cloacibacillus|s__Cloacibacillus_evryensis
tse
## class: TreeSummarizedExperiment 
## dim: 902 96 
## metadata(0):
## assays(1): relative_abundance
## rownames(902): Escherichia coli Klebsiella michiganensis ...
##   Pediococcus pentosaceus Rahnella aquatilis
## rowData names(7): superkingdom phylum ... genus species
## colnames(96): CA_C10001MS2003VA_t0M15 CA_C10003MS2031VA_t0M15 ...
##   MV_C118 MV_C120
## colData names(25): study_name subject_id ... family_role body_subsite
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: a LinkDataFrame (902 rows)
## rowTree: 1 phylo tree(s) (10430 leaves)
## colLinks: NULL
## colTree: NULL

Comparing mean to variance of the outcome variable

library(ggplot2)
mean(df$Streptococcus)
## [1] 4049.558
var(df$Streptococcus)
## [1] 509068831
fraczeroes <- nrow(df[df$Streptococcus == 0,])/nrow(df)
fraczeroes
## [1] 0.8736842

#Fit log-linear models

fit.pois <- glm(Streptococcus ~ Lactobacillus + age + BMI,
                data = df,
                family = poisson(link = "log"))

summary (fit.pois)
## 
## Call:
## glm(formula = Streptococcus ~ Lactobacillus + age + BMI, family = poisson(link = "log"), 
##     data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -206.14   -36.50   -19.90   -10.99   541.19  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -5.313e-01  1.873e-02  -28.36   <2e-16 ***
## Lactobacillus  5.141e-08  9.828e-10   52.32   <2e-16 ***
## age            3.599e-01  4.136e-04  870.25   <2e-16 ***
## BMI           -1.573e-01  6.453e-04 -243.77   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1929049  on 80  degrees of freedom
## Residual deviance:  632598  on 77  degrees of freedom
##   (14 observations deleted due to missingness)
## AIC: 632704
## 
## Number of Fisher Scoring iterations: 20
library(pscl)
fit.ZIpois <- zeroinfl(Streptococcus ~ Lactobacillus + age + BMI | 1,
                data = df,
                dist = "poisson")

summary (fit.ZIpois)
## 
## Call:
## zeroinfl(formula = Streptococcus ~ Lactobacillus + age + BMI | 1, data = df, 
##     dist = "poisson")
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3753 -0.3753 -0.3752 -0.3750 27.8168 
## 
## Count model coefficients (poisson with log link):
##                 Estimate Std. Error z value Pr(>|z|)
## (Intercept)    6.169e+00        NaN     NaN      NaN
## Lactobacillus  1.321e-07        NaN     NaN      NaN
## age            2.420e-01        NaN     NaN      NaN
## BMI           -2.116e-01        NaN     NaN      NaN
## 
## Zero-inflation model coefficients (binomial with logit link):
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)     1.96        NaN     NaN      NaN
## 
## Number of iterations in BFGS optimization: 1 
## Log-likelihood: -1.02e+05 on 5 Df
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
fit.negbin <-
  glm.nb(
    Streptococcus ~ Lactobacillus + age + BMI,
    data = df,
    control = glm.control(
      epsilon = 1e-8,
      maxit = 25,
      trace = FALSE
    )
  )
## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: algorithm did not converge
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning: glm.fit: algorithm did not converge
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: algorithm did not converge
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning: glm.fit: algorithm did not converge
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: algorithm did not converge
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning: glm.fit: algorithm did not converge
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: algorithm did not converge
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning: glm.fit: algorithm did not converge
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: algorithm did not converge
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning: glm.fit: algorithm did not converge
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in sqrt(1/i): NaNs produced
## Warning in glm.nb(Streptococcus ~ Lactobacillus + age + BMI, data = df, :
## alternation limit reached
summary(fit.negbin)
## 
## Call:
## glm.nb(formula = Streptococcus ~ Lactobacillus + age + BMI, data = df, 
##     control = glm.control(epsilon = 1e-08, maxit = 25, trace = FALSE), 
##     init.theta = 0.01199558153, link = log)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -206.1   -36.5   -19.9   -11.0   541.2  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -5.313e-01  1.873e-02  -28.36   <2e-16 ***
## Lactobacillus  5.141e-08  9.828e-10   52.32   <2e-16 ***
## age            3.599e-01  4.136e-04  870.25   <2e-16 ***
## BMI           -1.573e-01  6.453e-04 -243.77   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(5.380153e+14) family taken to be 1)
## 
##     Null deviance: 1929048  on 80  degrees of freedom
## Residual deviance:  632598  on 77  degrees of freedom
##   (14 observations deleted due to missingness)
## AIC: 280.39
## 
## Number of Fisher Scoring iterations: 5
## 
## 
##               Theta:  0.01200 
##           Std. Err.:  0.00395 
## Warning while fitting theta: alternation limit reached 
## 
##  2 x log-likelihood:  -270.38700
fit.ZInegbin <- zeroinfl(Streptococcus ~ Lactobacillus + age + BMI | 1,
                data = df,
                dist = "negbin")

summary(fit.ZInegbin)
## 
## Call:
## zeroinfl(formula = Streptococcus ~ Lactobacillus + age + BMI | 1, data = df, 
##     dist = "negbin")
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -0.1898 -0.1898 -0.1898 -0.1898  7.0566 
## 
## Count model coefficients (negbin with log link):
##                 Estimate Std. Error z value Pr(>|z|)
## (Intercept)    1.945e+01        NaN     NaN      NaN
## Lactobacillus -2.787e-07        NaN     NaN      NaN
## age            2.837e-02        NaN     NaN      NaN
## BMI           -4.085e-01        NaN     NaN      NaN
## Log(theta)    -9.837e-01        NaN     NaN      NaN
## 
## Zero-inflation model coefficients (binomial with logit link):
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)    1.921        NaN     NaN      NaN
## 
## Theta = 0.3739 
## Number of iterations in BFGS optimization: 14 
## Log-likelihood: -131.8 on 6 Df