lefser: a Megatanomic biomarker discovery tool
2024-07-21
Contents
1 Overview
1.1 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 groups of interest, referred as ‘classes’. 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.
1.2 Install and load pacakge
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("lefser")library(lefser)2 Analysis example
2.1 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).
data(zeller14)
zeller14 <- zeller14[, zeller14$study_condition != "adenoma"]The class and subclass information is stored in the colData slot under the
study_condition and age_category columns, respectively.
## Contingency table
table(zeller14$age_category, zeller14$study_condition)#>
#> 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
lefser(zeller14, groupCol = "study_condition", blockCol = "age_category")Warning messages:
1: In lefser(zeller14, groupCol = "study_condition", blockCol = "age_category") :
Convert counts to relative abundances with 'relativeAb()'
2: In lda.default(x, grouping, ...) : variables are collinear2.1.1 Terminal 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.
tn <- get_terminal_nodes(rownames(zeller14))
zeller14tn <- zeller14[tn,]2.1.2 Relative abundance
First warning message informs you that lefser requires relative abundance of
features. You can use relativeAb function to reformat your input.
zeller14tn_ra <- relativeAb(zeller14tn)2.2 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
groupCol = "study_condition",
blockCol = "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__Bacilli|o__Bacillales|f__Bacillales_noname|g__Gemella|s__Gemella_morbillorum|t__GCF_000185645
#> 5 k__Bacteria|p__Fusobacteria|c__Fusobacteriia|o__Fusobacteriales|f__Fusobacteriaceae|g__Fusobacterium|s__Fusobacterium_nucleatum|t__Fusobacterium_nucleatum_unclassified
#> 6 k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Clostridiaceae|g__Clostridium|s__Clostridium_symbiosum|t__Clostridium_symbiosum_unclassified
#> scores
#> 1 -3.336170
#> 2 -3.246547
#> 3 -3.060033
#> 4 -2.893966
#> 5 -2.788997
#> 6 -2.6312982.3 Visualization using lefserPlot
lefserPlot(res)
3 Benchmarking againt other tools
The codes for benchmarking lefser against LEfSe and the other R implementation of LEfSe is available here.
4 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.
mia::makeTreeSummarizedExperimentFromPhyloseq(kostic)#> 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: NULL5 Session Info
sessionInfo()#> R version 4.4.1 (2024-06-14)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.20-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.49.0 lefser_1.15.4
#> [3] SummarizedExperiment_1.35.1 Biobase_2.65.0
#> [5] GenomicRanges_1.57.1 GenomeInfoDb_1.41.1
#> [7] IRanges_2.39.2 S4Vectors_0.43.2
#> [9] BiocGenerics_0.51.0 MatrixGenerics_1.17.0
#> [11] matrixStats_1.3.0 BiocStyle_2.33.1
#>
#> loaded via a namespace (and not attached):
#> [1] splines_4.4.1 tibble_3.2.1
#> [3] rpart_4.1.23 DirichletMultinomial_1.47.0
#> [5] lifecycle_1.0.4 lattice_0.22-6
#> [7] MASS_7.3-61 MultiAssayExperiment_1.31.4
#> [9] backports_1.5.0 magrittr_2.0.3
#> [11] Hmisc_5.1-3 sass_0.4.9
#> [13] rmarkdown_2.27 jquerylib_0.1.4
#> [15] yaml_2.3.9 DBI_1.2.3
#> [17] minqa_1.2.7 ade4_1.7-22
#> [19] multcomp_1.4-26 abind_1.4-5
#> [21] zlibbioc_1.51.1 purrr_1.0.2
#> [23] yulab.utils_0.1.4 nnet_7.3-19
#> [25] TH.data_1.1-2 sandwich_3.1-0
#> [27] GenomeInfoDbData_1.2.12 ggrepel_0.9.5
#> [29] irlba_2.3.5.1 tidytree_0.4.6
#> [31] testthat_3.2.1.1 vegan_2.6-6.1
#> [33] rbiom_1.0.3 permute_0.9-7
#> [35] DelayedMatrixStats_1.27.2 codetools_0.2-20
#> [37] coin_1.4-3 DelayedArray_0.31.9
#> [39] scuttle_1.15.2 tidyselect_1.2.1
#> [41] UCSC.utils_1.1.0 farver_2.1.2
#> [43] lme4_1.1-35.5 ScaledMatrix_1.13.0
#> [45] viridis_0.6.5 base64enc_0.1-3
#> [47] jsonlite_1.8.8 BiocNeighbors_1.23.0
#> [49] multtest_2.61.0 decontam_1.25.0
#> [51] mia_1.13.27 Formula_1.2-5
#> [53] survival_3.7-0 scater_1.33.4
#> [55] iterators_1.0.14 foreach_1.5.2
#> [57] tools_4.4.1 treeio_1.29.0
#> [59] Rcpp_1.0.13 glue_1.7.0
#> [61] gridExtra_2.3 SparseArray_1.5.25
#> [63] xfun_0.46 mgcv_1.9-1
#> [65] TreeSummarizedExperiment_2.13.0 dplyr_1.1.4
#> [67] withr_3.0.0 BiocManager_1.30.23
#> [69] fastmap_1.2.0 boot_1.3-30
#> [71] rhdf5filters_1.17.0 bluster_1.15.0
#> [73] fansi_1.0.6 digest_0.6.36
#> [75] rsvd_1.0.5 R6_2.5.1
#> [77] colorspace_2.1-0 lpSolve_5.6.20
#> [79] utf8_1.2.4 tidyr_1.3.1
#> [81] generics_0.1.3 data.table_1.15.4
#> [83] DECIPHER_3.1.4 httr_1.4.7
#> [85] htmlwidgets_1.6.4 S4Arrays_1.5.5
#> [87] pkgconfig_2.0.3 gtable_0.3.5
#> [89] modeltools_0.2-23 SingleCellExperiment_1.27.2
#> [91] XVector_0.45.0 brio_1.1.5
#> [93] htmltools_0.5.8.1 bookdown_0.40
#> [95] biomformat_1.33.0 scales_1.3.0
#> [97] knitr_1.48 rstudioapi_0.16.0
#> [99] reshape2_1.4.4 checkmate_2.3.1
#> [101] nlme_3.1-165 nloptr_2.1.1
#> [103] cachem_1.1.0 zoo_1.8-12
#> [105] rhdf5_2.49.0 stringr_1.5.1
#> [107] parallel_4.4.1 vipor_0.4.7
#> [109] libcoin_1.0-10 foreign_0.8-87
#> [111] pillar_1.9.0 grid_4.4.1
#> [113] vctrs_0.6.5 slam_0.1-51
#> [115] BiocSingular_1.21.2 beachmat_2.21.4
#> [117] cluster_2.1.6 beeswarm_0.4.0
#> [119] htmlTable_2.4.3 evaluate_0.24.0
#> [121] tinytex_0.52 magick_2.8.4
#> [123] mvtnorm_1.2-5 cli_3.6.3
#> [125] compiler_4.4.1 rlang_1.1.4
#> [127] crayon_1.5.3 labeling_0.4.3
#> [129] mediation_4.5.0 plyr_1.8.9
#> [131] fs_1.6.4 ggbeeswarm_0.7.2
#> [133] stringi_1.8.4 viridisLite_0.4.2
#> [135] BiocParallel_1.39.0 munsell_0.5.1
#> [137] Biostrings_2.73.1 lazyeval_0.2.2
#> [139] Matrix_1.7-0 sparseMatrixStats_1.17.2
#> [141] ggplot2_3.5.1 Rhdf5lib_1.27.0
#> [143] highr_0.11 igraph_2.0.3
#> [145] memoise_2.0.1 RcppParallel_5.1.8
#> [147] bslib_0.7.0 ape_5.8