lefser: a Megatanomic biomarker discovery tool
Overview
Background
lefser is the R implementation of the Linear discriminant analysis (LDA) Effect Size (LEfSe), a Python package for metagenomic biomarker discovery and explanation. (Huttenhower et al. 2011).
The original software utilizes standard statistical significance tests along with supplementary tests that incorporate biological consistency and the relevance of effects to identity the features (e.g., organisms, clades, OTU, genes, or functions) that are most likely to account for differences between the two sample classes of interest. While LEfSe is widely used and available in different platform such as Galaxy UI and Conda, there is no convenient way to incorporate it in R-based workflows. Thus, we re-implement LEfSe as an R/Bioconductor package, lefser. Following the LEfSe‘s algorithm including Kruskal-Wallis test, Wilcoxon-Rank Sum test, and Linear Discriminant Analysis, with some modifications, lefser successfully reproduces and improves the original statistical method and the associated plotting functionality.
Analysis example
Prepare input
lefser package include the demo dataset,
zeller14, which is the microbiome data from colorectal
cancer (CRC) patients and controls (Zeller
et al. 2014).
In this vignette, we excluded the ‘adenoma’ condition and used control/CRC as the main classes and age category as sub-classes (adult vs. senior) with different numbers of samples: control-adult (n = 46), control-senior (n = 20), CRC-adult (n = 45), and CRC-senior (n = 46).
The class and subclass information is stored in the
colData slot under the study_condition and
age_category columns, respectively.
#>
#> CRC control
#> adult 45 46
#> senior 46 20If you try to run lefser directly on the ‘zeller14’
data, you will get the following warning messages
Warning messages:
1: In lefser(zeller14, classCol = "study_condition", subclassCol = "age_category") :
Convert counts to relative abundances with 'relativeAb()'
2: In lda.default(x, classing, ...) : variables are collinearTerminal node
When working with taxonomic data, including both terminal and non-terminal nodes in the analysis can lead to collinearity problems. Non-terminal nodes (e.g., genus) are often linearly dependent on their corresponding terminal nodes (e.g., species) since the species-level information is essentially a subset or more specific representation of the genus-level information. This collinearity can violate the assumptions of certain statistical methods, such as linear discriminant analysis (LDA), and can lead to unstable or unreliable results. By using only terminal nodes, you can effectively eliminate this collinearity issue, ensuring that your analysis is not affected by linearly dependent or highly correlated variables. Additionally, you can benefit of avoiding redundancy, increasing specificity, simplifying data, and reducing ambiguity, using only terminal nodes.
You can select only the terminal node using
get_terminal_nodes function.
Run lefser
The lefser function returns a data.frame
with two columns - the names of selected features (the
features column) and their effect size (the
scores column).
There is a random number generation step in the lefser
algorithm to ensure that more than half of the values for each features
are unique. In most cases, inputs are sparse, so in practice, this step
is handling 0s. So to reproduce the identical result, you should set the
seed before running lefser.
set.seed(1234)
res <- lefser(zeller14tn_ra, # relative abundance only with terminal nodes
classCol = "study_condition",
subclassCol = "age_category")
head(res)#> features
#> 1 k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Oscillospiraceae|g__Oscillibacter|s__Oscillibacter_unclassified
#> 2 k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Peptostreptococcaceae|g__Peptostreptococcus|s__Peptostreptococcus_stomatis|t__GCF_000147675
#> 3 k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales|f__Porphyromonadaceae|g__Porphyromonas|s__Porphyromonas_asaccharolytica|t__Porphyromonas_asaccharolytica_unclassified
#> 4 k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Clostridiaceae|g__Clostridium|s__Clostridium_symbiosum|t__Clostridium_symbiosum_unclassified
#> 5 k__Bacteria|p__Firmicutes|c__Bacilli|o__Bacillales|f__Bacillales_noname|g__Gemella|s__Gemella_morbillorum|t__GCF_000185645
#> 6 k__Bacteria|p__Fusobacteria|c__Fusobacteriia|o__Fusobacteriales|f__Fusobacteriaceae|g__Fusobacterium|s__Fusobacterium_nucleatum|t__Fusobacterium_nucleatum_unclassified
#> scores
#> 1 -3.336170
#> 2 -2.941230
#> 3 -2.834329
#> 4 -2.706471
#> 5 -2.579108
#> 6 -2.431915
Benchmarking againt other tools
The codes for benchmarking lefser against LEfSe and the other R implementation of LEfSe is available here.
Interoperating with phyloseq
When using phyloseq objects, we recommend to extract the
data and create a SummarizedExperiment object as
follows:
library(phyloseq)
library(SummarizedExperiment)
## Load phyloseq object
fp <- system.file("extdata",
"study_1457_split_library_seqs_and_mapping.zip",
package = "phyloseq")
kostic <- 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...## Split data tables
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.
#> class: TreeSummarizedExperiment
#> dim: 2505 190
#> metadata(0):
#> assays(1): counts
#> rownames(2505): 304309 469478 ... 206906 298806
#> rowData names(7): Kingdom Phylum ... Genus Species
#> colnames(190): C0333.N.518126 C0333.T.518046 ... 32I9UNA9.518098
#> BFJMKNMP.518102
#> colData names(71): X.SampleID BarcodeSequence ... HOST_TAXID
#> Description
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> rowLinks: NULL
#> rowTree: NULL
#> colLinks: NULL
#> colTree: NULLSession Info
#> R version 4.4.1 (2024-06-14)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 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.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 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: 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.50.0 lefser_1.17.1
#> [3] SummarizedExperiment_1.36.0 Biobase_2.67.0
#> [5] GenomicRanges_1.59.0 GenomeInfoDb_1.43.0
#> [7] IRanges_2.41.0 S4Vectors_0.44.0
#> [9] BiocGenerics_0.53.0 MatrixGenerics_1.19.0
#> [11] matrixStats_1.4.1 BiocStyle_2.35.0
#>
#> loaded via a namespace (and not attached):
#> [1] splines_4.4.1 ggplotify_0.1.2
#> [3] tibble_3.2.1 rpart_4.1.23
#> [5] DirichletMultinomial_1.49.0 lifecycle_1.0.4
#> [7] lattice_0.22-6 MASS_7.3-61
#> [9] MultiAssayExperiment_1.33.0 backports_1.5.0
#> [11] magrittr_2.0.3 Hmisc_5.2-0
#> [13] sass_0.4.9 rmarkdown_2.28
#> [15] jquerylib_0.1.4 yaml_2.3.10
#> [17] DBI_1.2.3 buildtools_1.0.0
#> [19] minqa_1.2.8 ade4_1.7-22
#> [21] multcomp_1.4-26 abind_1.4-8
#> [23] zlibbioc_1.52.0 purrr_1.0.2
#> [25] yulab.utils_0.1.7 nnet_7.3-19
#> [27] TH.data_1.1-2 sandwich_3.1-1
#> [29] GenomeInfoDbData_1.2.13 ggrepel_0.9.6
#> [31] irlba_2.3.5.1 tidytree_0.4.6
#> [33] maketools_1.3.1 testthat_3.2.1.1
#> [35] vegan_2.6-8 rbiom_1.0.3
#> [37] permute_0.9-7 DelayedMatrixStats_1.29.0
#> [39] codetools_0.2-20 coin_1.4-3
#> [41] DelayedArray_0.33.1 scuttle_1.16.0
#> [43] tidyselect_1.2.1 aplot_0.2.3
#> [45] UCSC.utils_1.2.0 farver_2.1.2
#> [47] lme4_1.1-35.5 ScaledMatrix_1.14.0
#> [49] viridis_0.6.5 base64enc_0.1-3
#> [51] jsonlite_1.8.9 BiocNeighbors_2.1.0
#> [53] multtest_2.63.0 decontam_1.27.0
#> [55] mia_1.15.0 Formula_1.2-5
#> [57] survival_3.7-0 scater_1.34.0
#> [59] iterators_1.0.14 foreach_1.5.2
#> [61] tools_4.4.1 treeio_1.30.0
#> [63] Rcpp_1.0.13 glue_1.8.0
#> [65] gridExtra_2.3 SparseArray_1.6.0
#> [67] xfun_0.48 mgcv_1.9-1
#> [69] TreeSummarizedExperiment_2.14.0 dplyr_1.1.4
#> [71] withr_3.0.2 BiocManager_1.30.25
#> [73] fastmap_1.2.0 boot_1.3-31
#> [75] rhdf5filters_1.18.0 bluster_1.17.0
#> [77] fansi_1.0.6 digest_0.6.37
#> [79] rsvd_1.0.5 R6_2.5.1
#> [81] gridGraphics_0.5-1 colorspace_2.1-1
#> [83] lpSolve_5.6.21 utf8_1.2.4
#> [85] tidyr_1.3.1 generics_0.1.3
#> [87] data.table_1.16.2 DECIPHER_3.3.0
#> [89] httr_1.4.7 htmlwidgets_1.6.4
#> [91] S4Arrays_1.6.0 pkgconfig_2.0.3
#> [93] gtable_0.3.6 modeltools_0.2-23
#> [95] SingleCellExperiment_1.28.0 XVector_0.46.0
#> [97] sys_3.4.3 brio_1.1.5
#> [99] htmltools_0.5.8.1 biomformat_1.35.0
#> [101] scales_1.3.0 ggfun_0.1.7
#> [103] knitr_1.48 rstudioapi_0.17.1
#> [105] reshape2_1.4.4 checkmate_2.3.2
#> [107] nlme_3.1-166 nloptr_2.1.1
#> [109] cachem_1.1.0 zoo_1.8-12
#> [111] rhdf5_2.50.0 stringr_1.5.1
#> [113] parallel_4.4.1 vipor_0.4.7
#> [115] libcoin_1.0-10 foreign_0.8-87
#> [117] pillar_1.9.0 grid_4.4.1
#> [119] vctrs_0.6.5 slam_0.1-54
#> [121] BiocSingular_1.23.0 beachmat_2.23.0
#> [123] cluster_2.1.6 beeswarm_0.4.0
#> [125] htmlTable_2.4.3 evaluate_1.0.1
#> [127] mvtnorm_1.3-1 cli_3.6.3
#> [129] compiler_4.4.1 rlang_1.1.4
#> [131] crayon_1.5.3 labeling_0.4.3
#> [133] mediation_4.5.0 plyr_1.8.9
#> [135] fs_1.6.5 ggbeeswarm_0.7.2
#> [137] stringi_1.8.4 viridisLite_0.4.2
#> [139] BiocParallel_1.41.0 munsell_0.5.1
#> [141] Biostrings_2.75.0 lazyeval_0.2.2
#> [143] Matrix_1.7-1 patchwork_1.3.0
#> [145] sparseMatrixStats_1.18.0 ggplot2_3.5.1
#> [147] Rhdf5lib_1.28.0 highr_0.11
#> [149] igraph_2.1.1 RcppParallel_5.1.9
#> [151] bslib_0.8.0 ggtree_3.15.0
#> [153] ape_5.8