Last updated: 2025-05-23
Checks: 7 0
Knit directory: CosMx_pipeline_LGA/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20250517) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 51b3104. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Ignored files:
Ignored: .Rproj.user/
Ignored: NBClust-Plots/
Ignored: analysis/figure/
Ignored: data/flatFiles/CoronalHemisphere/Run1000_S1_Half_exprMat_file.csv
Ignored: data/flatFiles/CoronalHemisphere/Run1000_S1_Half_fov_positions_file.csv
Ignored: data/flatFiles/CoronalHemisphere/Run1000_S1_Half_metadata_file.csv
Ignored: output/processed_data/Log/
Ignored: output/processed_data/RC/
Ignored: output/processed_data/SCT/
Ignored: output/processed_data/exprMat_unfiltered.RDS
Ignored: output/processed_data/fov_positions_unfiltered.RDS
Ignored: output/processed_data/metadata_unfiltered.RDS
Ignored: output/processed_data/negMat_unfiltered.RDS
Ignored: output/processed_data/seu_filtered.RDS
Ignored: output/processed_data/seu_semifiltered.RDS
Untracked files:
Untracked: analysis/images/
Unstaged changes:
Modified: analysis/_site.yml
Modified: output/performance_reports/0.0_data_loading_PR.csv
Modified: output/performance_reports/1.0_qc_and_filtering_PR.csv
Modified: output/performance_reports/2.0_normalization_PR.csv
Modified: output/performance_reports/3.0_dimensional_reduction_PR.csv
Modified: output/performance_reports/4.0_insitutype_cell_typing_PR.csv
Modified: output/performance_reports/4.1_insitutype_unsup_clustering_PR.csv
Modified: output/performance_reports/4.2_seurat_unsup_clustering_PR.csv
Modified: output/performance_reports/5.0_RC_normalization_PR.csv
Modified: output/performance_reports/5.1_RC_dimensional_reduction_PR.csv
Modified: output/performance_reports/6.0_Log_normalization_PR.csv
Modified: output/performance_reports/6.1_Log_dimensional_reduction_PR.csv
Modified: output/performance_reports/pipeline_PR.csv
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown (analysis/0.0_data_loading.Rmd) and
HTML (docs/0.0_data_loading.html) files. If you’ve
configured a remote Git repository (see ?wflow_git_remote),
click on the hyperlinks in the table below to view the files as they
were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| html | 51b3104 | lidiaga | 2025-05-23 | Build site. |
| html | 8a9f079 | lidiaga | 2025-05-20 | Build site. |
| Rmd | 8410a10 | lidiaga | 2025-05-17 | Add Rmds files in analysis |
library(data.table) # Efficient data management
library(Matrix) # Sparse matrices
library(here) # Enhanced file referencing in project-oriented workflows
library(vecsets) # Enhanced version of the Set Operations from "base"
library(dplyr) # For the use of pipes %>%
library(kableExtra) # For table formatting
library(Seurat) # Seurat objectFirst of all, flat files directory needs to be located in the folder structure of our project.
## Code adapted from Scratch Space vignette
# Locate flat files
flatFiles_dir <- here("data", "flatFiles")
slideNames <- dir(flatFiles_dir)In this case, there is only one slide, this will simplify early pipeline construction. Future improvements could include: adapt the code to work with multiple slides and perform integration.
## Code adapted from Scratch Space vignette
# Locate slide files
slide <- slideNames[1]
slide_dir <- here(flatFiles_dir, slide)
slide_files <- dir(slide_dir, full.names = TRUE)
# Load slide files
### Fov_positions
fov_positions <- fread(slide_files[grep("fov_positions_file.csv", slide_files)], header = TRUE)
## fread() is an efficient version of read.csv, that is able to read gzip files
## and creates a data.table object for efficient management
### Metadata
metadata <- fread(slide_files[grep("metadata_file.csv", slide_files)], header = TRUE)
### Counts matrix (efficient loading)
exprMat_dir <- slide_files[grep("exprMat_file.csv", slide_files)]
exprMat <- loadExprMat(exprMat_dir, slide_number = 1)
## loadExprMat() is a self created function that reads the expression matrix
## in chunks and converts it to sparse (adapted from Scratch Space vignette)
## available code is in "code/aux_functions.R" fov x_global_px y_global_px
<int> <num> <num>
1: 1 -498383.3 6475.000
2: 2 -498383.3 2208.333
3: 3 -498383.3 -2058.333
4: 4 -498383.3 -6325.000
5: 5 -498383.3 -10591.667
6: 6 -498383.3 -14858.333At first sight, there are no signs of previous pre-processing or missing FOVs.
fov cell_ID Area AspectRatio CenterX_local_px CenterY_local_px
<int> <int> <int> <num> <int> <int>
1: 1 1 6073 0.47 4222 4066
2: 1 2 5675 0.72 4182 3938
3: 1 3 12896 1.26 2156 3864
4: 1 4 8234 0.51 4108 3608
5: 1 5 9852 0.88 4162 3506
6: 1 6 13372 0.90 4167 3301In sight of the file columns, there are no signs of previous pre-processing or missing FOVs.
# Check presence of rows with cell_id of 0
cells_0 <- grep("_0$", rownames(exprMat), value = TRUE)
# Check presence of Negative and SystemControl probes columns
negPrb_col <- grep("^Neg", colnames(exprMat), value = TRUE) # Neg columns
sysCon_col <- grep("^SystemControl", colnames(exprMat), value = TRUE) # SystemControl columns
# Show head expression matrix
as.matrix(exprMat[1:6, 1:10]) Chrna4 Slc6a1 Cd109 Ldha Aldoc Drd1 Tank Rit2 Prkag2 Lpar3
c_1_1_0 9 58 11 30 938 6 28 20 15 2
c_1_1_1 0 0 0 0 0 0 0 0 0 0
c_1_1_2 0 0 0 0 0 0 0 0 1 0
c_1_1_3 5 1 0 1 2 0 0 0 0 0
c_1_1_4 1 1 0 1 0 0 0 2 1 0
c_1_1_5 0 0 0 2 2 0 0 0 0 0Although the number of rows/cells in the metadata and the expression matrix should coincide, the expression matrix presents 130 extra rows. This difference is due to the presence of cells of id “0”, which, in CosMx SMI data, are used to store the transcripts counts that can not be assigned to any cell (see an example ReadMe SMI Data File). Therefore, those extra rows represent the sum of unassigned transcripts from the correspondent FOV and can be used to evaluate segmentation.
Apart from that, some manipulation of the data is suspected, since the negative probes, which are named “Negative” on the raw AtoMx output, are found as “NegPrb”, and the System Control codes columns are not present, suggesting its previous removal.
Before creating the Seurat object, some further formatting of the data might be helpful:
For further analysis, the expression matrix must only have the interest features in columns and real cells in rows, with no additional information.
As previously stated, rows with a “cell_ID” of 0 do not represent actual cells, and false positive controls (Negative probes and System Control codes) do not represent features of interest, therefore, they must be removed from the expression matrix. However, they posses interesting information that should be kept in the metadata for future exploration.
# Calculate number of transcripts not assigned to any cell (by FOV)
cells_0 <- grep("_0$", rownames(exprMat)) # Rows with cell_id 0
negPrb_col <- grep("^NegPrb", colnames(exprMat)) # NegPrb columns
unassigned <- data.table(cell_id = character(0), # Data.table to save the info
fov = numeric(0),
unassigned_tx = numeric(0))
for (i in cells_0) {
cell_id <- rownames(exprMat)[i]
fov <- as.numeric(strsplit(cell_id, "_")[[1]][3])
unassigned_tx <- as.numeric(rowSums(exprMat[i, -negPrb_col, drop = FALSE])) # Only real Tx
new_row <- list(cell_id, fov, unassigned_tx)
unassigned <- rbindlist(list(unassigned, new_row))
}
# Remove rows of unassigned transcripts from exprMat
exprMat <- exprMat[-cells_0,]
# Calculate nCount of Negative and System Control probes
sysCon_col <- grep("^SystemControl", colnames(exprMat)) # SystemControl columns
falsePrb <- data.table(cell_id = character(0), # Data.table to save the info
nCount_negprobes = numeric(0),
nCount_falsecode = numeric(0))
nCount_negprobes <- rowSums(exprMat[, negPrb_col])
nCount_falsecode <- rowSums(exprMat[, sysCon_col])
if (identical(names(nCount_falsecode), names(nCount_negprobes))){
cell_id <- names(nCount_negprobes)
new_row <- list(cell_id, nCount_negprobes, nCount_falsecode)
falsePrb <- rbindlist(list(falsePrb, new_row))
} else {
stop("nCount_negprobes and nCount_falsecode don't share cell ids!")
}
# Save Negative probes as a separate matrix (for future Insitutype)
negMat <- exprMat[, negPrb_col]
# Remove Negative and SystemControl probes from exprMat
exprMat <- exprMat[, -c(negPrb_col, sysCon_col)]After removing the cells with id “0”, the number of rows/cells match that of the metadata. Also, the number of columns/features has been diminish after removing the negative probes.
Depending on the AtoMx version from which the data was exported, this file might contain the FOV coordinates in px (older versions) or in mm (newer versions) (see CosMxLite).
The following code will check if the px coordinates are present or not, and calculate them if needed.
## Code adapted from CosMxLite vignette
# Check and/or calculate px coordinates
if (all(c("x_global_px", "y_global_px") %in% colnames(fov_positions))) {
result <- "FOV coordinates in pixels are already present."
} else {
fov_positions$x_global_px <- fov_positions$x_global_mm*1000/0.1202809 # converting mm to px
fov_positions$y_global_px <- fov_positions$y_global_mm*1000/0.1202809 # converting mm to px
result <- "FOV coordinates in pixels have been calculated."
}FOV coordinates in pixels are already present.
In some cases, the exported metadata file from AtoMx might contain some extra information on analysis performed in said software. This information increases the file size and might not be used when performing the analysis from the beginning. Therefore, it can me removed by selecting only the basic columns, or those of interest for your analysis.
The following code adapted from CosMxLite will subset the metadata by selecting the basic columns only. Additionally, several new columns will be created to facilitate further subsetting and visualizations:
Note, the “cell_ID” original column won’t be kept as it only provides FOV-specific identification of the cells and might be misleading.
## Code inspired by CosMxLite vignette
# Add new columns: slide_ID and cell_id
metadata[, slide_ID_numeric := 1]
slide_fov_cell <- paste0("c_", metadata$slide_ID_numeric, "_", metadata$fov, "_", metadata$cell_ID)
metadata[, cell_id := slide_fov_cell]
# Merge X and Y FOV coordinates by "fov"
metadata[fov_positions, on = "fov", `:=` (x_FOV_px = i.x_global_px,
y_FOV_px = i.y_global_px)]
# Merge unassigned_tx by "fov"
metadata[unassigned, on = "fov", `:=` (unassignedTranscripts = i.unassigned_tx)]
# Merge nCount_negprobes and nCount_falsecodes by "cell_id"
metadata[falsePrb, on = "cell_id", `:=` (nCount_negprobes = i.nCount_negprobes,
nCount_falsecode = i.nCount_falsecode)]
# Aisle basic columns (from CosMxLite) + recently added columns
columns_to_keep <- c("slide_ID_numeric", "fov", "cell_id", "Area", "AspectRatio",
"x_FOV_px", "y_FOV_px", "Width", "Height", "CenterX_local_px",
"CenterY_local_px", "CenterX_global_px", "CenterY_global_px",
colnames(metadata)[grep("Mean|Max", colnames(metadata))],
"unassignedTranscripts", "nCount_negprobes", "nCount_falsecode")
metadata <- metadata[, ..columns_to_keep]
# Re-name rows as the slide_fov_cell
rownames(metadata) <- metadata$cell_id
# Ensure that cell content and order in the ExprMat and metadata matches
if (!vsetequal(rownames(metadata), rownames(exprMat))) {
stop("Rows do not match between metadata and exprMat.")
} else if (!identical(rownames(metadata), rownames(exprMat))) {
exprMat <- exprMat[match(rownames(metadata), rownames(exprMat)), ]
}## Code adapted from Scratch Space vignette
if(!dir.exists(here("output","processed_data"))){
dir.create(here("output","processed_data"))
}
saveRDS(exprMat, here("output","processed_data","exprMat_unfiltered.RDS"))
saveRDS(metadata, here("output","processed_data","metadata_unfiltered.RDS"))
saveRDS(fov_positions, here("output","processed_data","fov_positions_unfiltered.RDS"))
saveRDS(negMat, here("output","processed_data","negMat_unfiltered.RDS"))The Seurat package has two built-in functions to read the flat files from CosMx:
In this pipeline, the Seurat object will be created with the “CreateSeuratObject” function, as seen in the CosMxLite vignette, using the data previously loaded and processed.
## Code adapted from CosMxLite vignette
# Create a basic Seurat object from the processed exprMat and metadata
seu <- CreateSeuratObject(counts = t(exprMat), project = "CosMx",
meta.data = metadata, min.cells = 3, min.features = 1)When creating a Seurat object, “nCount” and “nFeature” columns are automatically calculated and added to the meta.data slot. Now, the “unassignedTranscripts” column can be recalculated to represent the fraction of total transcripts.
# Calculate total transcripts per FOV and add to meta.data
nCount_FOV <- seu@meta.data %>%
group_by(fov) %>%
summarise(nCount_FOV = sum(nCount_RNA))
seu@meta.data <- seu@meta.data %>%
left_join(nCount_FOV, by = "fov")
# Re-assign rownames to metadata, lost by dplyr functions
rownames(seu@meta.data) <- seu$cell_id
# Re-calculate unassignedTranscripts to fraction through dividing by
# total detected transcripts: nCount_FOV + unassignedTranscripts
seu$unassignedTranscripts <- seu$unassignedTranscripts / (seu$nCount_FOV + seu$unassignedTranscripts)Finally, it is important to include cell coordinates for visualizations. In the CosMxLite vignette two alternatives on adding the spatial coordinates are shown:
In this pipeline, the “global FOV” methods will be used to maintain interoperability with the AtoMx Seurat output.
## Code adapted from CosMxLite vignette: in this vignette two alternatives
## on adding the spatial coordinates are shown.
# Add cell coordinates as a "global FOV"
# Note: x and y coordinates need to be swapped to orient the tissue correctly
cell_centroids <- data.frame(y = seu$CenterX_global_px,
x = seu$CenterY_global_px,
cell = colnames(seu))
cents <- list("centroids" = CreateCentroids(cell_centroids))
coords <- CreateFOV(coords = cents, type = c("centroids"), assay = "RNA")
seu[["globalFOV"]] <- coords
# Add cell coordinates as a "reduction" (left as an example)
# coordinates <- as.matrix(data.frame(CosMx_1 = seu$CenterX_global_px,
# CosMx_2 = seu$CenterY_global_px))
# seu[['Centroids']] <- CreateDimReducObject(embeddings = coordinates, key = "CosMx_",
# global = TRUE, assay = "RNA")| Chunk | Time_sec | Memory_Mb |
|---|---|---|
| Libraries | 2.48 | 147.0 |
| LocateFiles | 0.20 | 0.0 |
| LoadSlide_1 | 4.80 | 219.4 |
| ExploringFOVpositions | 0.21 | 0.7 |
| ExploringMetadata | 0.18 | 0.1 |
| ExploringExprMat | 0.23 | 0.0 |
| PreparingExprMat | 4.25 | 2.8 |
| PreparingFOVpositions | 0.18 | 0.0 |
| PreparingMetadata | 0.68 | 2.7 |
| SavingData | 11.37 | 0.0 |
| CreateSeuObj | 2.24 | 215.5 |
| UnassignedTxFrac | 0.24 | 3.3 |
| Coords | 0.42 | 4.9 |
| SavingSeuObj | 15.29 | 0.0 |
| Total | 42.77 | 596.4 |
R version 4.4.3 (2025-02-28 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 22000)
Matrix products: default
locale:
[1] LC_COLLATE=Spanish_Spain.utf8 LC_CTYPE=Spanish_Spain.utf8
[3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C
[5] LC_TIME=Spanish_Spain.utf8
time zone: Europe/Madrid
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] SeuratObject_4.1.4 Seurat_4.4.0 kableExtra_1.4.0 dplyr_1.1.4
[5] vecsets_1.4 here_1.0.1 Matrix_1.7-2 data.table_1.17.0
[9] workflowr_1.7.1
loaded via a namespace (and not attached):
[1] deldir_2.0-4 pbapply_1.7-2 gridExtra_2.3
[4] rlang_1.1.6 magrittr_2.0.3 git2r_0.35.0
[7] RcppAnnoy_0.0.22 spatstat.geom_3.3-5 matrixStats_1.5.0
[10] ggridges_0.5.6 compiler_4.4.3 getPass_0.2-4
[13] reshape2_1.4.4 png_0.1-8 systemfonts_1.2.1
[16] callr_3.7.6 vctrs_0.6.5 stringr_1.5.1
[19] pkgconfig_2.0.3 fastmap_1.2.0 promises_1.3.2
[22] rmarkdown_2.29 pracma_2.4.4 ps_1.9.0
[25] purrr_1.0.4 xfun_0.51 cachem_1.1.0
[28] jsonlite_1.9.1 goftest_1.2-3 later_1.4.1
[31] spatstat.utils_3.1-2 irlba_2.3.5.1 parallel_4.4.3
[34] cluster_2.1.8 R6_2.6.1 ica_1.0-3
[37] spatstat.data_3.1-6 bslib_0.9.0 stringi_1.8.4
[40] RColorBrewer_1.1-3 reticulate_1.41.0.1 spatstat.univar_3.1-2
[43] parallelly_1.43.0 scattermore_1.2 lmtest_0.9-40
[46] jquerylib_0.1.4 Rcpp_1.0.14 knitr_1.50
[49] tensor_1.5 future.apply_1.11.3 zoo_1.8-13
[52] sctransform_0.4.1 httpuv_1.6.15 splines_4.4.3
[55] igraph_2.1.4 tidyselect_1.2.1 abind_1.4-8
[58] rstudioapi_0.17.1 yaml_2.3.10 spatstat.random_3.3-2
[61] spatstat.explore_3.3-4 codetools_0.2-20 miniUI_0.1.1.1
[64] processx_3.8.6 listenv_0.9.1 plyr_1.8.9
[67] lattice_0.22-6 tibble_3.2.1 shiny_1.10.0
[70] ROCR_1.0-11 evaluate_1.0.3 Rtsne_0.17
[73] future_1.34.0 survival_3.8-3 polyclip_1.10-7
[76] xml2_1.3.7 fitdistrplus_1.2-2 pillar_1.10.1
[79] whisker_0.4.1 KernSmooth_2.23-26 plotly_4.10.4
[82] generics_0.1.3 rprojroot_2.0.4 sp_2.2-0
[85] ggplot2_3.5.1 munsell_0.5.1 scales_1.3.0
[88] globals_0.16.3 xtable_1.8-4 glue_1.8.0
[91] lazyeval_0.2.2 tools_4.4.3 RANN_2.6.2
[94] fs_1.6.5 leiden_0.4.3.1 cowplot_1.1.3
[97] grid_4.4.3 tidyr_1.3.1 colorspace_2.1-1
[100] nlme_3.1-167 patchwork_1.3.0 cli_3.6.4
[103] spatstat.sparse_3.1-0 viridisLite_0.4.2 svglite_2.1.3
[106] uwot_0.2.3 gtable_0.3.6 sass_0.4.9
[109] digest_0.6.37 progressr_0.15.1 ggrepel_0.9.6
[112] htmlwidgets_1.6.4 farver_2.1.2 htmltools_0.5.8.1
[115] lifecycle_1.0.4 httr_1.4.7 mime_0.12
[118] MASS_7.3-65