spammR microbiome evaluation
Oct 02, 2026
Package
spammR 0.99.26
This vignette walks you through an example of using spammR for metagenomic data across a geographic area.
1 R and other dependencies
We require a few packages for the mapping requirements that are beyond the basic spammR tools.
The following R packages are required: - spammR - sf relies on other non R-libraries (see Installing for more details on how to install sf including the dependencies)
The following non-R packages may be required: - GEOS - PROJ - GDAL
library(spammR)
library(ggplot2)
library(tibble)
library(magrittr)
library(dplyr)
## these two libraries are only used for geographic data
library(sf)
library(terra)
library(BiocFileCache)2 Metagenomic data formatting requirements
The data that we will be working with is a KEGG ortholog enrichment from the 1000 soils project. The 1000 soils project has the nice feature that each location has two sampling at different depths (‘top’ and ‘bottom’). We will use this circumstance in combination with the KO enrichment to generate a pathway enrichment for each sample location and depth and then visualize the differential between the two depths.
We have downloaded the 1000 soil data and stored it as two separate files that can be loaded directly in the package.
The first file is the actual KO measurements.
We can then create a metadata file that The KO metadata table will contain a mapping from KO identifier - row names in the Omics Measurement table - to a more descriptive annotation.For testing purposes this currently is only the row names & incremented integer values.
oneKSoilsOmicsMeta <- data.frame(KO = character(nrow(oneKSoilsData)))
oneKSoilsOmicsMeta["KO"] <- rownames(oneKSoilsData)
oneKSoilsOmicsMeta["annotation"] <- seq(1, by = 1,
length.out = nrow(oneKSoilsOmicsMeta))
oneKSoilsOmicsMeta["annotation"] <- lapply(oneKSoilsOmicsMeta["annotation"],
as.character)
DT::datatable(head(oneKSoilsOmicsMeta),options = list(scrollX = TRUE))Next we need to import the lat/long information of the locations where the
1000 soil project samples where extracted. The file was retrieved from the
1000 Soils Project Shiny App.
The metadata was retrieved by using the Query>Information tab. On the left side
(the “Information Menu”) we select SAMPLE_ID, latitude & longitude and
download the resulting mapping.
Not all Sample_ID entries have complete lat/long data. To generate a sf
object (see below) which will be needed as input to create a
terra:SpatRaster object the input table can only contain complete cases
(i.e. no NaN values). The sample metadata contains information about the
type of Biome, the layer of the soil sample (‘top’ or ‘bottom’), and the pH,
for example, as well as the coordinates.
data("oneKSoilsMeta")
DT::datatable(head(oneKSoilsMeta),options = list(scrollX = TRUE))ggplot2::ggplot(oneKSoilsMeta, aes(x = Core_Layer, y = pH, fill = biome_name)) +
geom_boxplot()## Warning: Removed 4 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
3 Geographic map formatting requirements
Maps can be downloaded from public repositories as a ‘shape’ file that can then be processed for coordinates and visualization. Here we show how to retrieve a shape file d from the US Census Bureau with the following parameters:
- Year = 2024
- Layer type = “States (and equivalent)”
Or downloaded in a zip file from the [census FTP site] (https://www2.census.gov/geo/tiger/TIGER2024/STATE/).
The code chunk below will not be executed but is displayed as example of how to load a downloaded shape file using the sf package.
us_map <- sf::st_read("../data/tl_2024_us_state/tl_2024_us_state.shp",
quiet = TRUE)
us_mapInstead will load the file from the supplied rda file:
This map contains the entire US, so for visualization purposes we focus only on the continental US.
us_continental <- c(
"WV", "FL", "IL", "MN", "MD", "RI", "ID", "NH", "NC", "VT",
"CT", "DE", "NM", "CA", "NJ", "WI", "OR", "NE", "PA", "WA",
"LA", "GA", "AL", "UT", "OH", "TX", "CO", "SC", "OK", "TN",
"WY", "ND", "KY", "ME", "NY", "NV", "MI", "AR", "MS", "MO",
"MT", "KS", "IN", "SD", "MA", "VA", "DC", "IA", "AZ"
)
s_us_cont <- us_map %>%
.[.$STUSPS %in% us_continental, ] %>%
vect()
s_us_cont## class : SpatVector
## geometry : polygons
## dimensions : 49, 15 (geometries, attributes)
## extent : -124.849, -66.88544, 24.39631, 49.38448 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat NAD83 (EPSG:4269)
## names : REGION DIVISION STATEFP STATENS GEOID GEOIDFQ STUSPS NAME LSAD MTFCC (and 5 more)
## type : <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## values : 3 5 54 01779805 54 0400000US54 WV West Virginia 00 G4000
## 3 5 12 00294478 12 0400000US12 FL Florida 00 G4000
## 2 3 17 01779784 17 0400000US17 IL Illinois 00 G4000
## ...Now that we have the coordinate information we need to convert it to pixels
using the terra package, which requires creating a rasterization template
using the”
terra::SpatRaster. From this object we can generate a rasterized image while
retaining the location information and ability to project the map into a
different map projection.
The terra package can read in the shape file as follows:
template <- terra::rast(
s_us_cont,
res = 0.01
)And then create a rasterized landmass of the continental US:
# us_raster <- rasterize(s_us_cont, template, background=0)
# us_raster <- rasterize(as.lines(s_us_cont), template, field='STATEFP',
# touches=TRUE)
us_raster <- terra::rasterize(s_us_cont, template, field = "STATEFP")3.1 Mapping soil locations to map
Next we generate an sf object from the soils coordinate data so that we can
then map it to the same rasterized image as the US map object:
oneKSoilsMeta <- oneKSoilsMeta[!is.na(oneKSoilsMeta$latitude) &
!is.na(oneKSoilsMeta$longitude), ]
coords <- sf::st_as_sf(oneKSoilsMeta, coords = c("longitude", "latitude"),
remove = TRUE)
sf::st_crs(coords) <- sf::st_crs(s_us_cont)
DT::datatable(coords,options = list(scrollX = TRUE))Using the same rasterization template as we used above, we can make the soils coordinates fit on pixes on the US map.
coords_rast <- terra::rasterize(coords, template)
coords_rast## class : SpatRaster
## size : 2499, 5796, 1 (nrow, ncol, nlyr)
## resolution : 0.01, 0.01 (x, y)
## extent : -124.849, -66.88897, 24.39631, 49.38631 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat NAD83 (EPSG:4269)
## source(s) : memory
## name : last
## min value : 1
## max value : 13.1.1 get the “Pixel” x/y (NOT lat/long x/y) of metagenome samples
Finally to retrieve the “pixel” coordinates instead of the x/y (lat/long)
coordinates that are recorded in the terra::SpatRaster objects the following
procedure is employed:
- use
extract()to get the cells which the coordinates are projected into in the raster image - use
lapply()get the raster coordinates (‘pixel’ x/y) for each row withterra::rowColFromCell() - use
cbind()to combine the cells, ‘pixel’ x/y and original sample coordinatedata.frame - use
dplyr::rename()to clean updata.frame
This process does come with slight differences in lat/long. Most likely due to the rasterization process?
cells <- terra::extract(coords_rast, coords, xy = TRUE, cells = TRUE)
DT::datatable(head(cells),options = list(scrollX = TRUE))cells_tmp <- terra::extract(coords_rast, coords, xy = TRUE, cells = TRUE)
xy_rast <- lapply(cells["cell"], function(x) rowColFromCell(coords_rast, x))
cells <- cells_tmp %>%
cbind(., oneKSoilsMeta, xy_rast) %>%
dplyr::rename(., y_pixels = "cell.1", x_pixels = "cell.2") %>%
subset(., select = -c(x, y, ID, last))
cells$y_pixels <- nrow(coords_rast) - cells$y_pixelsThose need to be finally written out / transferred into the spammR package
for further use.
3.1.2 Exporting the rasterized map
The rasterized map can be exported with the command below.
terra::writeRaster(us_raster, "us_map.png", NAflag = 0, overwrite = TRUE)## Warning in x@pntr$writeRaster(opt): GDAL Message 6: PNG driver doesn't support data type Int32. Only eight bit (Byte) and sixteen bit (UInt16) bands supported. Defaulting to Byte3.1.3 Examplary plotting
coords_test <- sf::st_as_sf(oneKSoilsMeta,
coords = c("longitude", "latitude"),
remove = TRUE)
sf::st_crs(coords_test) <- sf::st_crs(s_us_cont)
plot(us_raster, col = "grey")
# lines(as.polygons(us_raster), col='black')
plot(coords_test, col = "darkred", add = TRUE)## Warning in plot.sf(coords_test, col = "darkred", add = TRUE): ignoring all but
## the first attribute
4 Using spammR on metagenomic geographic data
Now that the data are properly formatted we can put them into a
SpatialExperiment object that is needed by spammR.
The data is comprised of:
- Omics Measurements: species abundances across regions
- Omics Metadata: metadata from these species
- Sample Metadata: information about the sample location
- Image files (generated by the above code): file representing image sample
4.1 Sample Metadata
Here we build the sample metadata table given the previously generated coordinate raster during the rasterization of the map. The SpatialExperiment object requires information about the x and y coordinates and information about each spot in the image for which we have omics data.
dm <- cells
rownames(dm) <- dm$Sample_ID
# dm[,1] <- NULL
dm[, "cell"] <- NULL
dm[, "longitude.1"] <- NULL
dm[, "latitude.1"] <- NULL
dm[, "Image"] <- 0
dm[, "x_max"] <- ncol(coords_rast)
dm[, "y_max"] <- nrow(coords_rast)
dm[, "x_origin"] <- 0
dm[, "y_origin"] <- 0
dm[, "spot_height"] <- 500
dm[, "spot_width"] <- 500
dm[, "above_40_deg_lat"] <- dplyr::if_else(dm[, "latitude"] >= 40,
true = 1, false = 0)
dm[, "psychrophile"] <- dplyr::if_else(dm[, "Soil.Temperature.C"] <= 15,
true = 1, false = 0)
dm[, "Core_Layer_bin"] <- dplyr::if_else(dm[, "Core_Layer"] == "TOP",
true = "TOP", false = "BTM")
dm[, "Desert"] <- dplyr::if_else(dm[, "BIOME"] == "13",
true = "Desert", false = "NonDesert")
data_meta <- dm
DT::datatable(head(data_meta),options = list(scrollX = TRUE))4.2 CLEANING UP OF DATA
Limiting our selfs to only the top layer of the soil column.
data_meta <- subset(data_meta, Sample_ID %in% colnames(oneKSoilsData))
oneKSoilsData <- subset(oneKSoilsData, select = data_meta$Sample_ID)
data_meta <- data_meta[data_meta$Core_Layer == "TOP", ]
oneKSoilsData <- oneKSoilsData[, grep("TOP", x = names(oneKSoilsData))]4.3 Create SpatialExperiment object
Now that our metadata is formatted correctly we can create the spatialExperiment object.
microbiome.spe <- spammR::convert_to_spe(
oneKSoilsData, ## pooled data table
data_meta, ## pooled metadata
oneKSoilsOmicsMeta, ## protein identifiers
feature_meta_colname = "KO", # column name
spatial_coords_colnames = c("x_pixels", "y_pixels"),
image_files = c("us_map.png"),
image_sample_ids = "1000 Soils",
sample_id = "1000 Soils",
image_ids = c("0"),
assay_name = "KO"
)4.4 Metagenomic differential expression and enrichment
The pathway enrichment analysis is done based on the presence and absence
of KO associations found within the individual metagenomes. As such we are
not strictly working with “gene abundances” that would be expected for a
more traditional differential expression and pathway enrichment analysis.
Conversely the “abundance” measurements of the individual KOs also distribute
differently to what would be expected. To this end we are relying on the count
based differential expression analysis that is based on limma::voom.
diff_ex <- spammR::calc_spatial_diff_ex(
microbiome.spe,
assay_name = "KO",
count_based = TRUE,
log_transformed = FALSE,
category_col = "Desert"
)## We found 562 features with a logFC greater than 1 and
## an ajusted p-value less than 0.05sig_ko <- subset(rowData(diff_ex), Desert_vs_NonDesert.adj.P.Val.limma < 0.05)
ups <- subset(sig_ko, Desert_vs_NonDesert.logFC.limma > 0)
downs <- subset(sig_ko, Desert_vs_NonDesert.logFC.limma < 0)
print(paste(
"We found", nrow(sig_ko), "significantly differentally abundant KOs",
nrow(ups), "upregulated KOs and", nrow(downs), "downregulated"
))## [1] "We found 575 significantly differentally abundant KOs 483 upregulated KOs and 92 downregulated"ups %>% as_tibble(.) %>%
arrange(., -Desert_vs_NonDesert.logFC.limma)|>
DT::datatable(options = list(scrollX = TRUE))4.5 Spatial heatmap showing abundance of significant KO
K09758 (aspartate 4-decarboxylase) being one of the “upregulated KOs” we can
visualize the abundances across the different sampling locations with
spammr::spatial_heatmap()
spatial_heatmap(
microbiome.spe,
feature = "K09758",
# feature_type = 'KO',
assay_name = "KO",
sample_id = "1000 Soils",
image_id = "0",
sample_label_size = 2.0,
title_size = 10,
label_column = "Desert",
metric_display = "KO association abundance",
plot_title = "Abundance of aspartate 4-decarboxylase (K09758)",
interactive = FALSE
)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_label()`).
4.6 Pathway enrichment
For the pathway enrichment part of the anaysis we first need to load “genesets” in the required format for leapR that will be called internally. The geneset is effectively a list of KOs that are associated with individual KEGG pathways.
kk <- download.file("https://api.figshare.com/v2/file/download/55158818",
destfile = "kegg.rda", quiet = TRUE)
load("kegg.rda")
file.remove("kegg.rda")## [1] TRUENext we are performing the enrichment analysis with a p-value cut-off of p-value < 0.05
ora.res <- spammR::enrich_ora(diff_ex, geneset = keggPathwayToKo,
geneset_name = "KEGG pathways",
feature_column = "KO", pval_thresh = 0.05)
filtered_res <- ora.res %>%
rownames_to_column(., var = 'KEGG_Pathway') %>%
as_tibble(.) %>%
filter(., oddsratio > 1 & BH_pvalue < 0.05) %>%
arrange(-oddsratio) %>%
select(KEGG_Pathway, oddsratio, BH_pvalue)
pathway_names <- cbind(keggPathwayToKo$name, keggPathwayToKo$desc) %>%
as_tibble() %>%
rename(., KEGG_Pathway = V1, KEGG_pathway_name = V2)## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
## `.name_repair` is omitted as of tibble 2.0.0.
## ℹ Using compatibility `.name_repair`.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.left_join(filtered_res, pathway_names, by = 'KEGG_Pathway')|>
DT::datatable(options = list(scrollX = TRUE))ora.res <- tibble::rownames_to_column(ora.res, "pathways")
pathway_prots <- subset(ora.res, pathways == "map01320") |>
dplyr::select(ingroupnames)
pathway_prots <- unlist(strsplit(pathway_prots[1, ], split = ", "))spatial_heatmap(
microbiome.spe,
feature = pathway_prots,
assay_name = "KO",
sample_id = "1000 Soils",
image_id = "0",
label_column = "Desert",
sample_label_size = 2.0,
title_size = 10,
metric_display = "Enriched KO Groups",
plot_title = "Sulfur cycle",
interactive = FALSE
)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_label()`).
4.7 Session Information
sessionInfo()## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.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] terra_1.9-50 sf_1.1-3
## [3] magrittr_2.0.5 leapR_1.1.2
## [5] BiocFileCache_3.3.0 dbplyr_2.6.0
## [7] ggraph_2.2.2 cowplot_1.2.0
## [9] readxl_1.5.0.1 tidygraph_1.3.1
## [11] lubridate_1.9.5 forcats_1.0.1
## [13] stringr_1.6.0 dplyr_1.2.1
## [15] purrr_1.2.2 readr_2.2.0
## [17] tidyr_1.3.2 tibble_3.3.1
## [19] ggplot2_4.0.3 tidyverse_2.0.0
## [21] spammR_0.99.26 limma_3.99.0
## [23] SpatialExperiment_1.23.0 SingleCellExperiment_1.35.2
## [25] SummarizedExperiment_1.43.0 Biobase_2.73.2
## [27] GenomicRanges_1.65.4 Seqinfo_1.3.2
## [29] IRanges_2.47.5 S4Vectors_0.51.10
## [31] BiocGenerics_0.59.12 generics_0.1.4
## [33] MatrixGenerics_1.25.0 matrixStats_1.5.0
## [35] BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 jsonlite_2.0.0 wk_0.9.5
## [4] magick_2.9.1 farver_2.1.2 rmarkdown_2.32
## [7] vctrs_0.7.3 spdep_1.4-2 memoise_2.0.1
## [10] tinytex_0.61 rstatix_1.1.0 htmltools_0.5.9
## [13] S4Arrays_1.13.2 curl_8.0.0 broom_1.0.13
## [16] cellranger_1.1.0 s2_1.1.13 SparseArray_1.13.4
## [19] Formula_1.2-6 sass_0.4.10 spData_2.3.5
## [22] KernSmooth_2.23-27 bslib_0.12.0 htmlwidgets_1.6.4
## [25] httr2_1.3.0 impute_1.87.0 plotly_4.12.1
## [28] cachem_1.1.0 igraph_2.3.4 lifecycle_1.0.5
## [31] pkgconfig_2.0.3 Matrix_1.7-6 R6_2.6.1
## [34] fastmap_1.2.0 digest_0.6.39 ggnewscale_0.5.2
## [37] crosstalk_1.2.2 RSQLite_3.53.3 ggpubr_1.0.0
## [40] labeling_0.4.3 filelock_1.0.3 timechange_0.4.0
## [43] httr_1.4.9 polyclip_1.10-7 abind_1.4-8
## [46] compiler_4.6.1 proxy_0.4-29 bit64_4.8.6
## [49] withr_3.0.3 S7_0.2.2 backports_1.5.1
## [52] carData_3.0-6 viridis_0.6.5 DBI_1.3.0
## [55] ggforce_0.5.0 ggsignif_0.6.4 MASS_7.3-66
## [58] DelayedArray_0.39.8 rjson_0.2.23 classInt_0.4-11
## [61] tools_4.6.1 units_1.0-1 otel_0.2.0
## [64] glue_1.8.1 grid_4.6.1 gtable_0.3.6
## [67] tzdb_0.5.0 class_7.3-24 data.table_1.18.6.1
## [70] hms_1.1.4 utf8_1.2.6 sp_2.2-3
## [73] car_3.1-5 XVector_0.53.0 ggrepel_0.9.8
## [76] pillar_1.11.1 tweenr_2.0.3 lattice_0.23-1
## [79] bit_4.6.0 deldir_2.0-4 tidyselect_1.2.1
## [82] knitr_1.52 gridExtra_2.3.1 bookdown_0.48
## [85] xfun_0.61 graphlayouts_1.2.5 statmod_1.5.2
## [88] DT_0.34.0 stringi_1.8.9 yaml_2.3.12
## [91] boot_1.3-32 codetools_0.2-20 evaluate_1.0.5
## [94] BiocManager_1.30.27 cli_3.6.6 reticulate_1.47.0
## [97] jquerylib_0.1.4 dichromat_2.0-1 Rcpp_1.1.2
## [100] png_0.1-9 blob_1.3.0 viridisLite_0.4.3
## [103] scales_1.4.0 e1071_1.7-17 rlang_1.3.0