## ----message=FALSE, warning=FALSE, include=FALSE------------------------------ knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE) ## ----eval=FALSE--------------------------------------------------------------- # devtools::install_github("DavisLaboratory/standR") ## ----message = FALSE, warning = FALSE----------------------------------------- library(standR) library(SpatialExperiment) library(limma) library(ExperimentHub) ## ----message = FALSE, warning = FALSE----------------------------------------- eh <- ExperimentHub() query(eh, "standR") countFile <- eh[["EH7364"]] sampleAnnoFile <- eh[["EH7365"]] featureAnnoFile <- eh[["EH7366"]] spe <- readGeoMx(countFile, sampleAnnoFile, featureAnnoFile = featureAnnoFile, hasNegProbe = TRUE) ## ----------------------------------------------------------------------------- # check out object spe ## ----------------------------------------------------------------------------- SummarizedExperiment::assayNames(spe) ## ----------------------------------------------------------------------------- dge <- edgeR::SE2DGEList(spe) spe2 <- readGeoMxFromDGE(dge) spe2 ## ----------------------------------------------------------------------------- library(ggalluvial) colData(spe)$regions <- paste0(colData(spe)$region,"_",colData(spe)$SegmentLabel) |> (\(.) gsub("_Geometric Segment","",.))() |> paste0("_",colData(spe)$pathology) |> (\(.) gsub("_NA","_ns",.))() plotSampleInfo(spe, column2plot = c("SlideName","disease_status","regions")) ## ----------------------------------------------------------------------------- spe <- addPerROIQC(spe, rm_genes = TRUE) ## ----------------------------------------------------------------------------- dim(spe) metadata(spe) |> names() ## ----------------------------------------------------------------------------- plotGeneQC(spe, ordannots = "regions", col = regions, point_size = 2) ## ----------------------------------------------------------------------------- plotROIQC(spe, y_threshold = 50000, col = SlideName) ## ----------------------------------------------------------------------------- spe <- spe[,rownames(colData(spe))[colData(spe)$lib_size > 50000]] ## ----eval=FALSE--------------------------------------------------------------- # plotRLExpr(spe) ## ----------------------------------------------------------------------------- plotRLExpr(spe, ordannots = "SlideName", assay = 2, col = SlideName) ## ----------------------------------------------------------------------------- drawPCA(spe, assay = 2, col = SlideName, shape = regions) ## ----------------------------------------------------------------------------- standR::plotMDS(spe, assay = 2, col = SlideName, shape = regions) ## ----------------------------------------------------------------------------- spe <- scater::runUMAP(spe) plotDR(spe, dimred = "UMAP", col = regions) ## ----------------------------------------------------------------------------- plotScreePCA(spe, assay = 2, dims = 10) ## ----------------------------------------------------------------------------- plotPairPCA(spe, col = regions, shape = disease_status, assay = 2) ## ----------------------------------------------------------------------------- plotPCAbiplot(spe, n_loadings = 10, assay = 2, col = regions) ## ----------------------------------------------------------------------------- colData(spe)$biology <- paste0(colData(spe)$disease_status, "_", colData(spe)$regions) spe_tmm <- geomxNorm(spe, method = "TMM") ## ----------------------------------------------------------------------------- plotRLExpr(spe_tmm, assay = 2, color = SlideName) + ggtitle("TMM") ## ----------------------------------------------------------------------------- plotPairPCA(spe_tmm, assay = 2, color = disease_status, shape = regions) ## ----------------------------------------------------------------------------- spe <- findNCGs(spe, batch_name = "SlideName", top_n = 500) metadata(spe) |> names() ## ----------------------------------------------------------------------------- findBestK(spe, maxK = 10, factor_of_int = "biology", NCGs = metadata(spe)$NCGs, factor_batch = "SlideName") ## ----------------------------------------------------------------------------- spe_ruv <- geomxBatchCorrection(spe, factors = "biology", NCGs = metadata(spe)$NCGs, k = 5) ## ----------------------------------------------------------------------------- plotPairPCA(spe_ruv, assay = 2, color = disease_status, shape = regions, title = "RUV4") ## ----------------------------------------------------------------------------- spe_lrb <- geomxBatchCorrection(spe, batch = colData(spe)$SlideName, method = "Limma", design = model.matrix(~ 0 + disease_status + regions, data = colData(spe))) ## ----------------------------------------------------------------------------- plotPairPCA(spe_lrb, assay = 2, color = disease_status, shape = regions, title = "Limma removeBatch") ## ----------------------------------------------------------------------------- spe_list <- list(spe, spe_ruv, spe_lrb) plotClusterEvalStats(spe_list = spe_list, bio_feature_name = "regions", batch_feature_name = "SlideName", data_names = c("Raw","RUV4","Limma")) ## ----------------------------------------------------------------------------- plotRLExpr(spe_ruv, assay = 2, color = SlideName) + ggtitle("RUV4") plotRLExpr(spe_lrb, assay = 2, color = SlideName) + ggtitle("Limma removeBatch") ## ----------------------------------------------------------------------------- sessionInfo()