A quick start guide to the standR package
Ning Liu, Dharmesh Bhuva, Ahmed Mohamed, Chin Wee Tan, Melissa Davis
2023-06-13
1 Installation
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("standR")The development version of standR can be installed from GitHub:
2 Quick start
2.1 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(): 2023-04-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_featureAnno2.2 QC
2.2.1 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"))
2.2.2 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).
2.2.3 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.
2.3 Inspection of variations on ROI level
2.3.1 RLE
Here we can see obvious variation from slides to slides, and small variations are also observed within each slide.
2.3.2 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).
3 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.
4 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
5 SessionInfo
## R version 4.3.0 RC (2023-04-13 r84269)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.2 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.17-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] ggalluvial_0.12.5 ggplot2_3.4.2
## [3] ExperimentHub_2.8.0 AnnotationHub_3.8.0
## [5] BiocFileCache_2.8.0 dbplyr_2.3.2
## [7] limma_3.56.2 SpatialExperiment_1.10.0
## [9] SingleCellExperiment_1.22.0 SummarizedExperiment_1.30.2
## [11] Biobase_2.60.0 GenomicRanges_1.52.0
## [13] GenomeInfoDb_1.36.0 IRanges_2.34.0
## [15] S4Vectors_0.38.1 BiocGenerics_0.46.0
## [17] MatrixGenerics_1.12.2 matrixStats_1.0.0
## [19] standR_1.4.2
##
## loaded via a namespace (and not attached):
## [1] jsonlite_1.8.5 magrittr_2.0.3
## [3] ggbeeswarm_0.7.2 magick_2.7.4
## [5] farver_2.1.1 rmarkdown_2.22
## [7] zlibbioc_1.46.0 vctrs_0.6.2
## [9] memoise_2.0.1 DelayedMatrixStats_1.22.1
## [11] RCurl_1.98-1.12 rstatix_0.7.2
## [13] htmltools_0.5.5 S4Arrays_1.0.4
## [15] curl_5.0.1 broom_1.0.5
## [17] BiocNeighbors_1.18.0 Rhdf5lib_1.22.0
## [19] rhdf5_2.44.0 sass_0.4.6
## [21] bslib_0.5.0 cachem_1.0.8
## [23] mime_0.12 lifecycle_1.0.3
## [25] pkgconfig_2.0.3 rsvd_1.0.5
## [27] Matrix_1.5-4.1 R6_2.5.1
## [29] fastmap_1.1.1 GenomeInfoDbData_1.2.10
## [31] shiny_1.7.4 digest_0.6.31
## [33] colorspace_2.1-0 patchwork_1.1.2
## [35] AnnotationDbi_1.62.1 scater_1.28.0
## [37] irlba_2.3.5.1 dqrng_0.3.0
## [39] RSQLite_2.3.1 ggpubr_0.6.0
## [41] beachmat_2.16.0 filelock_1.0.2
## [43] labeling_0.4.2 fansi_1.0.4
## [45] abind_1.4-5 httr_1.4.6
## [47] mgcv_1.8-42 compiler_4.3.0
## [49] bit64_4.0.5 withr_2.5.0
## [51] backports_1.4.1 BiocParallel_1.34.2
## [53] carData_3.0-5 viridis_0.6.3
## [55] DBI_1.1.3 highr_0.10
## [57] HDF5Array_1.28.1 R.utils_2.12.2
## [59] ggsignif_0.6.4 rappdirs_0.3.3
## [61] DelayedArray_0.26.3 rjson_0.2.21
## [63] tools_4.3.0 vipor_0.4.5
## [65] beeswarm_0.4.0 interactiveDisplayBase_1.38.0
## [67] httpuv_1.6.11 R.oo_1.25.0
## [69] glue_1.6.2 nlme_3.1-162
## [71] rhdf5filters_1.12.1 promises_1.2.0.1
## [73] grid_4.3.0 generics_0.1.3
## [75] gtable_0.3.3 tzdb_0.4.0
## [77] R.methodsS3_1.8.2 tidyr_1.3.0
## [79] hms_1.1.3 car_3.1-2
## [81] ScaledMatrix_1.8.1 BiocSingular_1.16.0
## [83] utf8_1.2.3 XVector_0.40.0
## [85] ggrepel_0.9.3 BiocVersion_3.17.1
## [87] pillar_1.9.0 vroom_1.6.3
## [89] later_1.3.1 splines_4.3.0
## [91] dplyr_1.1.2 lattice_0.21-8
## [93] bit_4.0.5 ruv_0.9.7.1
## [95] tidyselect_1.2.0 locfit_1.5-9.8
## [97] Biostrings_2.68.1 scuttle_1.10.1
## [99] knitr_1.43 gridExtra_2.3
## [101] edgeR_3.42.4 xfun_0.39
## [103] DropletUtils_1.20.0 yaml_2.3.7
## [105] evaluate_0.21 codetools_0.2-19
## [107] archive_1.1.5 tibble_3.2.1
## [109] BiocManager_1.30.21 cli_3.6.1
## [111] xtable_1.8-4 munsell_0.5.0
## [113] jquerylib_0.1.4 Rcpp_1.0.10
## [115] png_0.1-8 parallel_4.3.0
## [117] ellipsis_0.3.2 readr_2.1.4
## [119] blob_1.2.4 sparseMatrixStats_1.12.0
## [121] bitops_1.0-7 viridisLite_0.4.2
## [123] scales_1.2.1 purrr_1.0.1
## [125] crayon_1.5.2 rlang_1.1.1
## [127] cowplot_1.1.1 KEGGREST_1.40.0