Skip to contents

Introduction

Robust principal component analysis (RPCA) methods require appropriate filtering of the raw count table before robust centered log-ratio transformation. Low-depth samples and low-abundance or low-prevalence features can otherwise contribute little information and increase computational cost.

The mia::filterRPCAInput() function provides Gemelli-style filters for SummarizedExperiment objects. Samples are retained when their total count is strictly greater than min.sample.count. Features are retained when both their total count and prevalence are strictly greater than the corresponding thresholds.

This vignette demonstrates the filtering step with executable R code. The complete cross-language comparison of mia and Gemelli is maintained as a separate validation report. See the OMAWorkflows package website for additional workflows.

Installation

Install OMAWorkflows and mia from Bioconductor:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
BiocManager::install(c("OMAWorkflows", "mia"))

Load the packages used in the example:

Example data

The GlobalPatterns dataset contains a raw count assay and is small enough for a package vignette.

data(GlobalPatterns, package = "mia")
tse <- GlobalPatterns
count_matrix <- assay(tse, "counts")

data.frame(
    features = nrow(tse),
    samples = ncol(tse),
    total_reads = sum(count_matrix)
)
  features samples total_reads
1    19216      26    28216678

Filter the RPCA input

This example retains samples with more than 10,000 reads and features with more than 100 total reads and prevalence greater than 10%.

min_sample_count <- 10000
min_feature_count <- 100
min_feature_frequency <- 0.10

filtered_tse <- filterRPCAInput(
    tse,
    assay.type = "counts",
    min.sample.count = min_sample_count,
    min.feature.count = min_feature_count,
    min.feature.frequency = min_feature_frequency
)

data.frame(
    object = c("Before filtering", "After filtering"),
    features = c(nrow(tse), nrow(filtered_tse)),
    samples = c(ncol(tse), ncol(filtered_tse))
)
            object features samples
1 Before filtering    19216      26
2  After filtering     6213      26

Verify the filtering rules

The following calculation independently reconstructs the expected identifiers. The checks also demonstrate that all three thresholds use strict greater-than comparisons.

expected_samples <- colSums(count_matrix) > min_sample_count
expected_features <- rowSums(count_matrix) > min_feature_count
feature_prevalence <-
    rowSums(count_matrix > 0) / ncol(count_matrix)
expected_features <- expected_features &
    feature_prevalence > min_feature_frequency

stopifnot(
    identical(
        colnames(filtered_tse),
        colnames(tse)[expected_samples]
    ),
    identical(
        rownames(filtered_tse),
        rownames(tse)[expected_features]
    )
)

data.frame(
    criterion = c(
        "Sample total count",
        "Feature total count",
        "Feature prevalence"
    ),
    threshold = c(
        min_sample_count,
        min_feature_count,
        min_feature_frequency
    ),
    comparison = rep("strictly greater than", 3L)
)
            criterion threshold            comparison
1  Sample total count     1e+04 strictly greater than
2 Feature total count     1e+02 strictly greater than
3  Feature prevalence     1e-01 strictly greater than

Individual filters can be disabled with NULL. This is useful when a threshold has already been applied upstream.

unfiltered_tse <- filterRPCAInput(
    tse,
    assay.type = "counts",
    min.sample.count = NULL,
    min.feature.count = NULL,
    min.feature.frequency = NULL
)

stopifnot(identical(dim(unfiltered_tse), dim(tse)))

Complete mia-Gemelli validation

The complete cross-language validation includes filtering and rCLR agreement, component alignment, sample and feature loadings, explained variation, sample distances, and cross-validation trajectories.

View the complete Joint-RPCA validation report or inspect the reproducibility files.

References

Martino, C., Shenhav, L., Marotz, C. A., et al. (2020). Context-aware dimensionality reduction deconvolutes gut microbial community dynamics. Nature Biotechnology, 39, 165-168. doi:10.1038/s41587-020-0660-7.

Session information

R version 4.6.1 (2026-06-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.4 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=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] mia_1.21.6                      TreeSummarizedExperiment_2.21.0
 [3] Biostrings_2.81.6               XVector_0.53.0
 [5] SingleCellExperiment_1.35.2     MultiAssayExperiment_1.39.0
 [7] SummarizedExperiment_1.43.0     Biobase_2.73.2
 [9] GenomicRanges_1.65.1            Seqinfo_1.3.0
[11] IRanges_2.47.2                  S4Vectors_0.51.6
[13] BiocGenerics_0.59.10            generics_0.1.4
[15] MatrixGenerics_1.25.0           matrixStats_1.5.0

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.1            viridisLite_0.4.3
 [3] dplyr_1.2.1                 vipor_0.4.7
 [5] farver_2.1.2                viridis_0.6.5
 [7] S7_0.2.2                    fastmap_1.2.0
 [9] lazyeval_0.2.3              bluster_1.23.0
[11] digest_0.6.39               rsvd_1.0.5
[13] DirichletMultinomial_1.55.0 lifecycle_1.0.5
[15] cluster_2.1.8.3             tidytree_0.4.8
[17] magrittr_2.0.5              compiler_4.6.1
[19] rlang_1.3.0                 tools_4.6.1
[21] igraph_2.3.3                yaml_2.3.12
[23] knitr_1.51                  S4Arrays_1.13.0
[25] DelayedArray_0.39.3         plyr_1.8.9
[27] RColorBrewer_1.1-3          abind_1.4-8
[29] BiocParallel_1.47.0         purrr_1.2.2
[31] grid_4.6.1                  beachmat_2.29.0
[33] ggplot2_4.0.3               MASS_7.3-66
[35] scales_1.4.0                cli_3.6.6
[37] vegan_2.7-5                 rmarkdown_2.31
[39] crayon_1.5.3                treeio_1.37.0
[41] otel_0.2.0                  reshape2_1.4.5
[43] DelayedMatrixStats_1.35.0   scuttle_1.23.1
[45] DBI_1.3.0                   ggbeeswarm_0.7.3
[47] ape_5.8-1                   stringr_1.6.0
[49] splines_4.6.1               parallel_4.6.1
[51] yulab.utils_0.2.4           vctrs_0.7.3
[53] Matrix_1.7-6                jsonlite_2.0.0
[55] BiocSingular_1.29.0         BiocNeighbors_2.7.2
[57] ggrepel_0.9.8               decontam_1.33.0
[59] irlba_2.3.7                 beeswarm_0.4.0
[61] scater_1.41.2               tidyr_1.3.2
[63] glue_1.8.1                  codetools_0.2-20
[65] stringi_1.8.7               gtable_0.3.6
[67] ScaledMatrix_1.21.0         tibble_3.3.1
[69] pillar_1.11.1               rappdirs_0.3.4
[71] htmltools_0.5.9             R6_2.6.1
[73] sparseMatrixStats_1.25.0    evaluate_1.0.5
[75] lattice_0.22-9              DECIPHER_3.9.2
[77] Rcpp_1.1.2                  permute_0.9-10
[79] gridExtra_2.3.1             SparseArray_1.13.2
[81] nlme_3.1-170                mgcv_1.9-4
[83] xfun_0.60                   fs_2.1.0
[85] pkgconfig_2.0.3