Reproducibility of a published result - NYCHANES study

library(bugphyzzAnalyses)
library(bugphyzz)
library(mia)
library(dplyr)
library(EnrichmentBrowser)
library(purrr)
library(bugsigdbr)

Expected results

Direct quotes from Tobacco exposure associated with oral microbiota oxygen utilization in the New York City Health and Nutrition Examination Study

These results were obtained using annotated OTUs based on their genus.

ORA

Significant result:

“… differentially abundant OTUs between current smokers and never smokers to be over-represented in aerobic OTUs (Hypergeometric test, p = 0.004).”

GSEA

Significant result:

“…aerobic OTUs were significantly depleted among current smokers relative to never smokers (Enrichment Score test statistic ES = −0.75, p = 0.002, GSEA permutation test).”

Non-significant results:

“Anaerobic OTUs were enriched in smokers relative to never smokers but the difference was not statistically significant (ES = 0.36, p = 0.14, GSEA permutation test).”

“… enrichment of facultative anaerobic OTUs among never smokers compared to current smokers but this result was not statistically significant (ES = −0.29, p = 0.48, GSEA permutation test).”

Data

Data used here is a subset of the original data used for the NYCHANES study. Full data can be found on GitHub. Signatures from this study have been published on BugSigDB.

The data subset presented here only contains abundance values for “Never smoker” vs “Cigarette” samples (“smokingstatus” column in metadata).

NYCHANES data - OTU:

tse <- importNYCHANES()
tse
#> class: TreeSummarizedExperiment 
#> dim: 1291 129 
#> metadata(0):
#> assays(1): counts
#> rownames(1291): New.CleanUp.ReferenceOTU7367
#>   New.CleanUp.ReferenceOTU7623 ... New.CleanUp.ReferenceOTU2608
#>   New.CleanUp.ReferenceOTU5014
#> rowData names(8): Domain Phylum ... Species GenusID
#> colnames(129): NYDH0033 NYDH0026 ... NYDH1460 NYDH0847
#> colData names(714): KEY Burklab_ID ... INC3C DBQ_10_3CAT
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> rowLinks: a LinkDataFrame (1291 rows)
#> rowTree: 1 phylo tree(s) (1291 leaves)
#> colLinks: NULL
#> colTree: NULL

table(tse$smokingstatus)
#> 
#> Never smoker    Cigarette 
#>           43           86

NYCHANES data - Genus:

tse_gn <- splitByRanks(tse)$Genus
## Filter taxa with at least 1 count as abundance in at least 20% of samples
# min_n_samples <- round(ncol(tse_gn) * 0.2)
# tse_gn_subset <- tse_gn[rowSums(assay(tse_gn) >= 1) >= min_n_samples,]
tse_gn
#> class: TreeSummarizedExperiment 
#> dim: 173 129 
#> metadata(1): agglomerated_by_rank
#> assays(1): counts
#> rownames(173): Streptococcus unclassified ... Aerococcus
#>   unclassified_22
#> rowData names(8): Domain Phylum ... Species GenusID
#> colnames(129): NYDH0033 NYDH0026 ... NYDH1460 NYDH0847
#> colData names(714): KEY Burklab_ID ... INC3C DBQ_10_3CAT
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> rowLinks: a LinkDataFrame (173 rows)
#> rowTree: 1 phylo tree(s) (1291 leaves)
#> colLinks: NULL
#> colTree: NULL

bugphyzz data:

discrete_types <- c(
    "binary", "multistate-intersection", "multistate-union"
)
bp <- importBugphyzz(v = 0.5) |> 
    map(~ mutate(.x, NCBI_ID = as.character(NCBI_ID))) |> 
    map(~ {
        attr_type <- unique(.x$Attribute_type)
        if (attr_type %in% discrete_types) {
            df <- .x |> 
                filter(
                    !(Validation < 0.7 & Evidence == "asr")
                )
        } else if (attr_type == "numeric") {
            df <- .x |> 
                filter(
                    !(Validation < 0.5 & Evidence == "asr")
                )
        }
        df
    })

aer <- bp$aerophilicity
dim(aer)
#> [1] 27286    12

Signatures

Create signatures at the genus level:

gn_sigs <- makeSignatures(
    dat = aer, taxIdType = 'Taxon_name', taxLevel = 'genus'
) |> 
    list_flatten(name_spec = "{inner") |> 
    discard(is.null)
gn_sigs_bk <- c(gn_sigs, list(background = rownames(tse_gn)))
map(gn_sigs_bk, head)
#> $`bugphyz:aerophilicity|aerobic`
#> [1] "Flavihumibacter"  "Parasegetibacter" "Ferruginibacter"  "Fimbriimonas"    
#> [5] "Camelimonas"      "Thiothrix"       
#> 
#> $`bugphyz:aerophilicity|anaerobic`
#> [1] "Microaerobacter"      "Desulfitispora"       "Desulfurispira"      
#> [4] "Pseudoflavonifractor" "Chromatium"           "Ectothiorhodospira"  
#> 
#> $`bugphyz:aerophilicity|facultatively anaerobic`
#> [1] "Capnocytophaga" "Kistimonas"     "Trueperella"    "Telmatobacter" 
#> [5] "Alishewanella"  "Muricauda"     
#> 
#> $background
#> [1] "Streptococcus"  "unclassified"   "Blautia"        "Gemella"       
#> [5] "unclassified_1" "Abiotrophia"

Create signatures at the OTU level based on genus annotations:

genera <- rowData(tse)$Genus 
names(genera) <- rownames(tse)
otu_sigs <- lapply(gn_sigs, function(x) names(genera)[which(genera %in% x)])
otu_sigs_bk <- c(otu_sigs, list(background = rownames(tse)))
map(otu_sigs_bk, head)
#> $`bugphyz:aerophilicity|aerobic`
#> [1] "New.CleanUp.ReferenceOTU9251"  "New.CleanUp.ReferenceOTU11022"
#> [3] "New.CleanUp.ReferenceOTU11549" "New.CleanUp.ReferenceOTU11122"
#> [5] "New.CleanUp.ReferenceOTU8077"  "New.CleanUp.ReferenceOTU11146"
#> 
#> $`bugphyz:aerophilicity|anaerobic`
#> [1] "New.CleanUp.ReferenceOTU8003"  "New.CleanUp.ReferenceOTU10697"
#> [3] "New.CleanUp.ReferenceOTU10070" "New.CleanUp.ReferenceOTU8703" 
#> [5] "New.CleanUp.ReferenceOTU3192"  "New.CleanUp.ReferenceOTU3925" 
#> 
#> $`bugphyz:aerophilicity|facultatively anaerobic`
#> [1] "New.CleanUp.ReferenceOTU7367" "New.CleanUp.ReferenceOTU7623"
#> [3] "New.ReferenceOTU195"          "New.CleanUp.ReferenceOTU2376"
#> [5] "New.CleanUp.ReferenceOTU3539" "New.CleanUp.ReferenceOTU7335"
#> 
#> $background
#> [1] "New.CleanUp.ReferenceOTU7367" "New.CleanUp.ReferenceOTU7623"
#> [3] "New.ReferenceOTU195"          "New.CleanUp.ReferenceOTU2376"
#> [5] "New.CleanUp.ReferenceOTU3539" "New.CleanUp.ReferenceOTU7335"

max_sig_size <- max(map_int(c(otu_sigs_bk, gn_sigs_bk), length))
max_sig_size
#> [1] 2422

Current value:

configEBrowser(key="GS.MAX.SIZE")
#> [1] 500

New max.size value:

configEBrowser(key="GS.MAX.SIZE", value = max_sig_size)
configEBrowser(key="GS.MAX.SIZE")
#> [1] 2422

Enrichment analysis at the OTU level

Differential abundance:

tse$GROUP <- ifelse(
  tse$smokingstatus == 'Never smoker', 0, 1
)
edger <- deAna(
  expr = tse, de.method = 'edgeR', padj.method = 'fdr', 
  filter.by.expr = FALSE
)
edger_norm <- limmaVoom(edger)

Enrichment (GSEA)

gsea <- sbea(
  method = 'gsea', se = edger_norm, gs = otu_sigs, perm = 1000, alpha = 0.1 # This is for DA taxa
)
gsea_res <- as.data.frame(gsea$res.tbl) |> 
    mutate(GENE.SET = ifelse(PVAL < 0.05, paste0(GENE.SET, "*"), GENE.SET)) |> 
    mutate(PVAL = round(PVAL, 3))
myDataTable(gsea_res)

Enrichment (ORA)

Significant (mixed) vs Non-Significant

ora <- sbea(
  method = 'ora', se = edger, gs = otu_sigs_bk, perm = 0, 
  alpha = 0.1 # This is for DA taxa
)
ora_res <- as.data.frame(ora$res.tbl) |> 
    mutate(GENE.SET = ifelse(PVAL < 0.05, paste0(GENE.SET, "*"), GENE.SET)) |> 
    filter(GENE.SET != "background") |> 
    mutate(PVAL = round(PVAL, 3))
myDataTable(ora_res)

Significant over-abundant vs Reference

edger_up <- edger
rowData(edger_up)$ADJ.PVAL <- ifelse(
  rowData(edger_up)$FC < 0, 1, rowData(edger_up)$ADJ.PVAL
)
ora_up <- sbea(
  method = 'ora', se = edger_up, gs = otu_sigs_bk, perm = 0, alpha = 0.1 # For DA taxa
)
ora_up_res <- as.data.frame(ora_up$res.tbl) |> 
    mutate(GENE.SET = ifelse(PVAL < 0.05, paste0(GENE.SET, "*"), GENE.SET)) |> 
    filter(GENE.SET != "background") |> 
    mutate(PVAL = round(PVAL, 3))
myDataTable(ora_up_res)

Significant under-abundant vs Reference

edger_down <- edger
rowData(edger_down)$ADJ.PVAL <- ifelse(
  rowData(edger_down)$FC < 1, 0, rowData(edger_down)$ADJ.PVAL
)
ora_down <- sbea(
  method = 'ora', se = edger_down, gs = otu_sigs_bk, perm = 0, 
  alpha = 0.1 # For DA taxa
)
ora_down_res <- as.data.frame(ora_down$res.tbl) |> 
    mutate(GENE.SET = ifelse(PVAL < 0.05, paste0(GENE.SET, "*"), GENE.SET)) |> 
    filter(GENE.SET != "background") |> 
    mutate(PVAL = round(PVAL, 3))
myDataTable(ora_down_res)

Analysis at the genus level with bugphyzz signatures

Differential abundance analysis

tse_gn$GROUP <- ifelse(
  tse_gn$smokingstatus == 'Never smoker', 0, 1
)
edger_gn <- deAna(
    expr = tse_gn, de.method = 'edgeR', padj.method = 'BH', 
    filter.by.expr = FALSE, 
)
## This step is to avoid filtering out several taxa from the
## next analysis (GSEA)
edger_gn <- limmaVoom(edger_gn)

Enrichment GSEA

gsea_gn <- sbea(
  method = 'gsea', se = edger_gn, gs = gn_sigs, perm = 1000,
  alpha = 0.1 # This FDR is for DA taxa
)
gsea_gn_res <- as.data.frame(gsea_gn$res.tbl) |> 
    mutate(GENE.SET = ifelse(PVAL < 0.05, paste0(GENE.SET, "*"), GENE.SET)) |> 
    mutate(PVAL = round(PVAL, 3))
myDataTable(gsea_gn_res)

Enrichment ORA

Significant (mixing up and down) vs Non-Significant

gn_sigs_bk <- c(gn_sigs, list(background = rownames(edger_gn)))
ora_gn <- sbea(
  method = 'ora', se = edger_gn, gs = gn_sigs_bk, perm = 0, 
  alpha = 0.1 # This is for significant genes in the SE
)
ora_gn_res <- as.data.frame(ora_gn$res.tbl) |> 
    mutate(GENE.SET = ifelse(PVAL < 0.05, paste0(GENE.SET, "*"), GENE.SET)) |> 
    filter(GENE.SET != "background") |> 
    mutate(PVAL = round(PVAL, 3))
myDataTable(ora_gn_res)