Lefser finds features that have greatest differences between classes.
2023-11-22
Contents
1 Introduction
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.
2 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.
library(lefser)3 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)##
## CRC control
## 91 66There 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)##
## CRC control
## adult 45 46
## senior 46 20We 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 'CRC'.
## See `?factor` or `?relevel` to change the reference category.head(res)## Names
## 1 k__Bacteria|p__Fusobacteria|c__Fusobacteriia|o__Fusobacteriales|f__Fusobacteriaceae
## 2 k__Bacteria|p__Fusobacteria|c__Fusobacteriia|o__Fusobacteriales|f__Fusobacteriaceae|g__Fusobacterium
## 3 k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales|f__Porphyromonadaceae|g__Porphyromonas
## 4 k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales|f__Porphyromonadaceae|g__Porphyromonas|s__Porphyromonas_asaccharolytica|t__Porphyromonas_asaccharolytica_unclassified
## 5 k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales|f__Porphyromonadaceae|g__Porphyromonas|s__Porphyromonas_asaccharolytica
## 6 k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Peptostreptococcaceae|g__Peptostreptococcus|s__Peptostreptococcus_stomatis|t__GCF_000147675
## scores
## 1 -5.220469
## 2 -5.184039
## 3 -5.053910
## 4 -4.972276
## 5 -4.947255
## 6 -4.8557564 Visualizing results with lefserPlot
lefserPlot(res)
5 Interoperating with phyloseq
When using phyloseq objects, we recommend to extract the data and create a
SummarizedExperiment object as follows:
library(phyloseq)##
## 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':
##
## distancefp <- 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).
5.1 sessionInfo
sessionInfo()## R version 4.3.2 Patched (2023-11-13 r85521)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [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: America/New_York
## 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.1
## [3] SummarizedExperiment_1.32.0 Biobase_2.62.0
## [5] GenomicRanges_1.54.1 GenomeInfoDb_1.38.1
## [7] IRanges_2.36.0 S4Vectors_0.40.1
## [9] BiocGenerics_0.48.1 MatrixGenerics_1.14.0
## [11] matrixStats_1.1.0 BiocStyle_2.30.0
##
## loaded via a namespace (and not attached):
## [1] ade4_1.7-22 tidyselect_1.2.0 libcoin_1.0-10
## [4] dplyr_1.1.4 farver_2.1.1 Biostrings_2.70.1
## [7] bitops_1.0-7 fastmap_1.1.1 RCurl_1.98-1.13
## [10] TH.data_1.1-2 digest_0.6.33 lifecycle_1.0.4
## [13] cluster_2.1.4 survival_3.5-7 magrittr_2.0.3
## [16] compiler_4.3.2 rlang_1.1.2 sass_0.4.7
## [19] tools_4.3.2 igraph_1.5.1 utf8_1.2.4
## [22] yaml_2.3.7 data.table_1.14.8 knitr_1.45
## [25] S4Arrays_1.2.0 labeling_0.4.3 DelayedArray_0.28.0
## [28] plyr_1.8.9 multcomp_1.4-25 abind_1.4-5
## [31] withr_2.5.2 grid_4.3.2 fansi_1.0.5
## [34] multtest_2.58.0 biomformat_1.30.0 colorspace_2.1-0
## [37] Rhdf5lib_1.24.0 ggplot2_3.4.4 scales_1.2.1
## [40] iterators_1.0.14 MASS_7.3-60 cli_3.6.1
## [43] mvtnorm_1.2-3 vegan_2.6-4 rmarkdown_2.25
## [46] crayon_1.5.2 generics_0.1.3 reshape2_1.4.4
## [49] rhdf5_2.46.0 ape_5.7-1 cachem_1.0.8
## [52] stringr_1.5.1 zlibbioc_1.48.0 modeltools_0.2-23
## [55] splines_4.3.2 parallel_4.3.2 BiocManager_1.30.22
## [58] XVector_0.42.0 vctrs_0.6.4 Matrix_1.6-3
## [61] sandwich_3.0-2 jsonlite_1.8.7 bookdown_0.36
## [64] magick_2.8.1 foreach_1.5.2 jquerylib_0.1.4
## [67] glue_1.6.2 codetools_0.2-19 stringi_1.8.1
## [70] gtable_0.3.4 munsell_0.5.0 tibble_3.2.1
## [73] pillar_1.9.0 rhdf5filters_1.14.1 htmltools_0.5.7
## [76] GenomeInfoDbData_1.2.11 R6_2.5.1 evaluate_0.23
## [79] lattice_0.22-5 highr_0.10 bslib_0.6.0
## [82] Rcpp_1.0.11 permute_0.9-7 nlme_3.1-163
## [85] SparseArray_1.2.2 mgcv_1.9-0 coin_1.4-3
## [88] xfun_0.41 zoo_1.8-12 pkgconfig_2.0.3