A quick start guide to the standR package
Installation
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("standR")The development version of standR can be installed from
GitHub:
Quick start
Load data for this guide
This is the background for the data:
NanoString GeoMx DSP dataset of diabetic kidney disease (DKD) vs healthy kidney tissue. Seven slides were analyzed, 4 DKD and 3 healthy. Regions of Interest (ROI) were focused two different parts of a kidney’s structure: tubules or glomeruli. Individual glomeruli were identified by a pathologist as either relatively healthy or diseased regardless if the tissue was DKD or healthy. Tubule ROIs were segmented into distal (PanCK) and proximal (neg) tubules. While both distal and proximal tubules are called tubules, they perform very different functions in the kidney.
## ExperimentHub with 3 records
## # snapshotDate(): 2024-10-24
## # $dataprovider: Nanostring
## # $species: NA
## # $rdataclass: data.frame
## # additional mcols(): taxonomyid, genome, description,
## # coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
## # rdatapath, sourceurl, sourcetype
## # retrieve records with, e.g., 'object[["EH7364"]]'
##
## title
## EH7364 | GeomxDKDdata_count
## EH7365 | GeomxDKDdata_sampleAnno
## EH7366 | GeomxDKDdata_featureAnnoQC
metadata visualization
Based on the description of the data, we know that all glomerulus are classified as abnormal and healthy, and tubule are classified as neg and PanCK.
We therefore merge the region-related annotations to avoid collinearity, which can affect the process of batch correction.
colData(spe)$regions <- paste0(colData(spe)$region,"_",colData(spe)$SegmentLabel) |>
(\(.) gsub("_Geometric Segment","",.))() |>
paste0("_",colData(spe)$pathology) |>
(\(.) gsub("_NA","_ns",.))()
library(ggalluvial)
plotSampleInfo(spe, column2plot = c("SlideName","disease_status","regions"))
Gene level QC
Using the plotGeneQC function, we can have a look at
which were the genes removed and the overall distribution of percentage
of non-expressed genes in all ROIs. By default, top 9 genes are plotted
here (arranging by mean expression), user can increase the number of
plotted genes by changing the top_n parameter.
In this case we don’t see any specific biological pattern for the samples expressing this genes (figure above).
ROI level QC
In the ROI level QC, we first aim to identify (if any) ROI(s) that have relatively low library size and low cell count because they are considered as low quality samples due to insufficient sequencing depth or lack of RNA in the chosen region.
In this case, looking at the distribution plots of library size and nuclei count, we don’t see any particular spike in the low ends, rather the distributions are relatively smooth. Looking at the dot plot, library sizes are mostly positively correlate with the nuclei count, with some data have relatively low library size while the nuclei count is reasonable. We therefore can try to draw an filtering threshold at the low end of the library size, in this case 50000. By coloring the dot with their slide names, we find that the ROIs below the threshold are all from slide disease1B, suggesting the reason for this might be some technical issues of slide disease1B.
Since library size of 50000 seems to be a reasonable threshold, here
we subset the spatial experiment object based on the library size in
colData.
Inspection of variations on ROI level
RLE
Here we can see obvious variation from slides to slides, and small variations are also observed within each slide.
PCA
Here we color the PCA with slide information, and shape by regions (tissue). We can see that PC1 is mainly spread out by regions, especially glomerulus and tubule. And grouping based on slide within each tissue are observed. The subtypes in tubule are clearly separated, but different subtypes of glomerulus is still grouping together. Moreover, diseased tissues and control tissues are mixed as well (disease slides and normal slides).
Data normalization
As we observed the technical variations in the data in both RLE and PCA plots. It is necessary to perform normalization in the data.
In the standR package, we offer normalization options
including TMM, RPKM, TPM, CPM, upperquartile and sizefactor. Among them,
RPKM and TPM required gene length information (add
genelength column to the rowData of the
object). For TMM, upperquartile and sizefactor, their normalized factor
will be stored their metadata.
Here we used TMM to normalize the data.
Batch correction
In the Nanostring’s GeoMX DSP protocol, due to the fact that one slide is only big enough for a handful of tissue segments (ROIs), it is common that we see the DSP data being confounded by the batch effect introduced by different slides. In order to establish fair comparison between ROIs later on, it is necessary to remove this batch effect from the data.
To run RUV4 batch correction, we need to provide a list of “negative control genes (NCGs)”.
The function findNCGs allows identifying the NCGs from
the data. In this case, since the batch effect is mostly introduced by
slide, we therefore want to identify NCGs across all slides, so here we
set the batch_name to “SlideName”, and select the top 500
least variable genes across different slides as NCGs.
## [1] "NegProbes" "lcpm_threshold" "genes_rm_rawCount"
## [4] "genes_rm_logCPM" "NCGs"Here we use k of 5 to perform RUV-4 normalization.
We can then inspect the PCA of the corrected data with annotations, to inspect the removal of batch effects, and the retaining of the biological factors.
Moreover, we can also have a look at the RLE plots of the normalized count.
For more detailed analysis pipeline and usage of the standR package, please see https://github.com/DavisLaboratory/GeoMXAnalysisWorkflow
SessionInfo
## R version 4.4.2 (2024-10-31)
## 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] ggalluvial_0.12.5 ggplot2_3.5.1
## [3] ExperimentHub_2.15.0 AnnotationHub_3.15.0
## [5] BiocFileCache_2.15.0 dbplyr_2.5.0
## [7] limma_3.63.2 SpatialExperiment_1.17.0
## [9] SingleCellExperiment_1.29.1 SummarizedExperiment_1.37.0
## [11] Biobase_2.67.0 GenomicRanges_1.59.0
## [13] GenomeInfoDb_1.43.1 IRanges_2.41.1
## [15] S4Vectors_0.45.2 BiocGenerics_0.53.3
## [17] generics_0.1.3 MatrixGenerics_1.19.0
## [19] matrixStats_1.4.1 standR_1.11.0
##
## loaded via a namespace (and not attached):
## [1] sys_3.4.3 jsonlite_1.8.9 magrittr_2.0.3
## [4] ggbeeswarm_0.7.2 magick_2.8.5 farver_2.1.2
## [7] rmarkdown_2.29 zlibbioc_1.52.0 vctrs_0.6.5
## [10] memoise_2.0.1 rstatix_0.7.2 htmltools_0.5.8.1
## [13] S4Arrays_1.7.1 curl_6.0.1 BiocNeighbors_2.1.0
## [16] broom_1.0.7 Formula_1.2-5 SparseArray_1.7.2
## [19] sass_0.4.9 bslib_0.8.0 cachem_1.1.0
## [22] buildtools_1.0.0 mime_0.12 lifecycle_1.0.4
## [25] pkgconfig_2.0.3 rsvd_1.0.5 Matrix_1.7-1
## [28] R6_2.5.1 fastmap_1.2.0 GenomeInfoDbData_1.2.13
## [31] digest_0.6.37 colorspace_2.1-1 patchwork_1.3.0
## [34] AnnotationDbi_1.69.0 scater_1.35.0 irlba_2.3.5.1
## [37] RSQLite_2.3.8 ggpubr_0.6.0 beachmat_2.23.1
## [40] filelock_1.0.3 labeling_0.4.3 fansi_1.0.6
## [43] httr_1.4.7 abind_1.4-8 mgcv_1.9-1
## [46] compiler_4.4.2 bit64_4.5.2 withr_3.0.2
## [49] backports_1.5.0 BiocParallel_1.41.0 carData_3.0-5
## [52] viridis_0.6.5 DBI_1.2.3 ggsignif_0.6.4
## [55] rappdirs_0.3.3 DelayedArray_0.33.2 rjson_0.2.23
## [58] tools_4.4.2 vipor_0.4.7 beeswarm_0.4.0
## [61] glue_1.8.0 nlme_3.1-166 grid_4.4.2
## [64] gtable_0.3.6 tzdb_0.4.0 tidyr_1.3.1
## [67] hms_1.1.3 BiocSingular_1.23.0 ScaledMatrix_1.15.0
## [70] car_3.1-3 utf8_1.2.4 XVector_0.47.0
## [73] ggrepel_0.9.6 BiocVersion_3.21.1 pillar_1.9.0
## [76] vroom_1.6.5 splines_4.4.2 dplyr_1.1.4
## [79] lattice_0.22-6 bit_4.5.0 ruv_0.9.7.1
## [82] tidyselect_1.2.1 locfit_1.5-9.10 maketools_1.3.1
## [85] Biostrings_2.75.1 scuttle_1.17.0 knitr_1.49
## [88] gridExtra_2.3 edgeR_4.5.0 xfun_0.49
## [91] statmod_1.5.0 UCSC.utils_1.3.0 yaml_2.3.10
## [94] evaluate_1.0.1 codetools_0.2-20 tibble_3.2.1
## [97] BiocManager_1.30.25 cli_3.6.3 munsell_0.5.1
## [100] jquerylib_0.1.4 Rcpp_1.0.13-1 png_0.1-8
## [103] parallel_4.4.2 readr_2.1.5 blob_1.2.4
## [106] viridisLite_0.4.2 scales_1.3.0 purrr_1.0.2
## [109] crayon_1.5.3 rlang_1.1.4 cowplot_1.1.3
## [112] KEGGREST_1.47.0