Lefser finds features that have greatest differences between classes.
Asya Khleborodova, Ludwig Geistlinger, and Levi
Waldron
School of Public Health, City University of New
York
School of Public Health, City University of New York
2024-03-06
Source:vignettes/lefser.Rmd
lefser.RmdIntroduction
Lefser is metagenomic biomarker discovery tool that is based on LEfSe tool and is
published by Huttenhower
et al. 2011. Lefser is the R implementation of the
LEfSe method.
Using statistical analyses, lefser compares microbial
populations of healthy and diseased subjects to discover differencially
expressed microorganisms. Lefser than computes effect size,
which estimates magnitude of differential expression between the
populations for each differentially expressed microorganism. Subclasses
of classes can also be assigned and used within the analysis.
Installation
To install Bioconductor and the lefser package, run the
following commands.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("lefser")Then load the lefser package.
Overview and example use of lefser
The lefser function can be used with a
SummarizedExperiment.
Load the zeller14 example dataset and exclude ‘adenoma’
conditions.
data(zeller14)
zeller14 <- zeller14[, zeller14$study_condition != "adenoma"]Note. lefser supports only two-group contrasts.
The colData in the SummarizedExperiment
dataset contains the grouping column study_condition which
includes the ‘control’ and ‘CRC’ groups.
table(zeller14$study_condition)##
## control CRC
## 66 91There can be subclasses in each group condition. In the example
dataset we include age_category as a subclass of
study_condition which includes ‘adults’ and ‘seniors’. This
variable will correspond to the blockCol input
argument.
table(zeller14$age_category)##
## adult senior
## 91 66We can create a contingency table for the two categorical variables.
table(zeller14$age_category, zeller14$study_condition)##
## control CRC
## adult 46 45
## senior 20 46We can now use the lefser function. It provides results
as a data.frame with the names of selected microorganisms
and their effect size.
res <- lefser(zeller14, groupCol = "study_condition", blockCol = "age_category")## Warning in lefser(zeller14, groupCol = "study_condition", blockCol =
## "age_category"): Convert counts to relative abundances with 'relativeAb()'## The outcome variable is specified as 'study_condition' and the reference category is 'control'.
## See `?factor` or `?relevel` to change the reference category.
head(res)## Names
## 1 k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Ruminococcaceae|g__Ruminococcus
## 2 k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales
## 3 k__Bacteria|p__Firmicutes|c__Bacilli
## 4 k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Ruminococcaceae|g__Ruminococcus|s__Ruminococcus_sp_5_1_39BFAA
## 5 k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Ruminococcaceae|g__Ruminococcus|s__Ruminococcus_sp_5_1_39BFAA|t__GCF_000159975
## 6 k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Eubacteriaceae|g__Eubacterium|s__Eubacterium_hallii
## scores
## 1 -6.193941
## 2 -5.914844
## 3 -5.903174
## 4 -5.558376
## 5 -5.558376
## 6 -5.507090
Interoperating with phyloseq
When using phyloseq objects, we recommend to extract the
data and create a SummarizedExperiment object as
follows:
##
## Attaching package: 'phyloseq'## The following object is masked from 'package:SummarizedExperiment':
##
## distance## The following object is masked from 'package:Biobase':
##
## sampleNames## The following object is masked from 'package:GenomicRanges':
##
## distance## The following object is masked from 'package:IRanges':
##
## distance
fp <- system.file(
"extdata", "study_1457_split_library_seqs_and_mapping.zip",
package = "phyloseq"
)
kostic <- suppressWarnings({
microbio_me_qiime(fp)
})## Found biom-format file, now parsing it...
## Done parsing biom...
## Importing Sample Metdadata from mapping file...
## Merging the imported objects...
## Successfully merged, phyloseq-class created.
## Returning...
counts <- unclass(otu_table(kostic))
colData <- as(sample_data(kostic), "data.frame")
## create a SummarizedExperiment object
SummarizedExperiment(
assays = list(counts = counts), colData = colData
)## class: SummarizedExperiment
## dim: 2505 190
## metadata(0):
## assays(1): counts
## rownames(2505): 304309 469478 ... 206906 298806
## rowData names(0):
## colnames(190): C0333.N.518126 C0333.T.518046 ... 32I9UNA9.518098
## BFJMKNMP.518102
## colData names(71): X.SampleID BarcodeSequence ... HOST_TAXID
## DescriptionYou may also consider using
makeTreeSummarizedExperimentFromPhyloseq from the
mia package (example not shown).
sessionInfo
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] phyloseq_1.46.0 lefser_1.12.2
## [3] SummarizedExperiment_1.32.0 Biobase_2.62.0
## [5] GenomicRanges_1.54.1 GenomeInfoDb_1.38.6
## [7] IRanges_2.36.0 S4Vectors_0.40.2
## [9] BiocGenerics_0.48.1 MatrixGenerics_1.14.0
## [11] matrixStats_1.2.0 BiocStyle_2.30.0
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-7 permute_0.9-7 sandwich_3.1-0
## [4] rlang_1.1.3 magrittr_2.0.3 multcomp_1.4-25
## [7] ade4_1.7-22 compiler_4.3.2 mgcv_1.9-1
## [10] systemfonts_1.0.5 vctrs_0.6.5 reshape2_1.4.4
## [13] stringr_1.5.1 pkgconfig_2.0.3 crayon_1.5.2
## [16] fastmap_1.1.1 XVector_0.42.0 labeling_0.4.3
## [19] utf8_1.2.4 rmarkdown_2.25 ragg_1.2.7
## [22] purrr_1.0.2 xfun_0.42 modeltools_0.2-23
## [25] zlibbioc_1.48.0 cachem_1.0.8 jsonlite_1.8.8
## [28] biomformat_1.30.0 highr_0.10 rhdf5filters_1.14.1
## [31] DelayedArray_0.28.0 Rhdf5lib_1.24.2 parallel_4.3.2
## [34] cluster_2.1.6 R6_2.5.1 coin_1.4-3
## [37] bslib_0.6.1 stringi_1.8.3 jquerylib_0.1.4
## [40] Rcpp_1.0.12 bookdown_0.38 iterators_1.0.14
## [43] knitr_1.45 zoo_1.8-12 Matrix_1.6-5
## [46] splines_4.3.2 igraph_2.0.2 tidyselect_1.2.0
## [49] abind_1.4-5 yaml_2.3.8 vegan_2.6-4
## [52] codetools_0.2-19 lattice_0.22-5 tibble_3.2.1
## [55] plyr_1.8.9 withr_3.0.0 evaluate_0.23
## [58] desc_1.4.3 survival_3.5-8 Biostrings_2.70.2
## [61] pillar_1.9.0 BiocManager_1.30.22 foreach_1.5.2
## [64] generics_0.1.3 RCurl_1.98-1.14 ggplot2_3.5.0
## [67] munsell_0.5.0 scales_1.3.0 glue_1.7.0
## [70] tools_4.3.2 data.table_1.15.2 fs_1.6.3
## [73] mvtnorm_1.2-4 rhdf5_2.46.1 grid_4.3.2
## [76] ape_5.7-1 libcoin_1.0-10 colorspace_2.1-0
## [79] nlme_3.1-164 GenomeInfoDbData_1.2.11 cli_3.6.2
## [82] textshaping_0.3.7 fansi_1.0.6 S4Arrays_1.2.1
## [85] dplyr_1.1.4 gtable_0.3.4 sass_0.4.8
## [88] digest_0.6.34 SparseArray_1.2.4 TH.data_1.1-2
## [91] farver_2.1.1 memoise_2.0.1 htmltools_0.5.7
## [94] pkgdown_2.0.7 multtest_2.58.0 lifecycle_1.0.4
## [97] MASS_7.3-60.0.1