## ----setup, include=FALSE-----------------------------------------------------
# Set chunk options: suppress echo, messages, and warnings in code output
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

## ----load-cosmx---------------------------------------------------------------
library(SpaceTrooper)
library(ggplot2)

cosmxFolder <- system.file(file.path("extdata", "CosMx_DBKero_Tiny"), 
                        package="SpaceTrooper")

spe <- readCosmxSPE(dirName=cosmxFolder, 
                    sampleName="DBKero_CosMx",
                    coordNames=c("CenterX_global_px", "CenterY_global_px"),
                    countMatFPattern="exprMat_file.csv",
                    metadataFPattern="metadata_file.csv",
                    polygonsFPattern="polygons.csv",
                    fovPosFPattern="fov_positions_file.csv",
                    fovdims=c(xdim = 4256, ydim = 4256)
)

## ----load-cosmx-protein-------------------------------------------------------
protfolder <- system.file( "extdata", "S01_prot", package="SpaceTrooper")

prot <- readCosmxProteinSPE(dirName=protfolder, 
                           sampleName="CosMx_Protein_Tonsil",
                           coordNames=c("CenterX_global_px", "CenterY_global_px"),
                           countMatFPattern="exprMat_file.csv", 
                           metadataFPattern="metadata_file.csv", 
                           polygonsFPattern="polygons.csv", 
                           fovPosFPattern="fov_positions_file.csv", 
                           fovdims=c(xdim = 4256, ydim = 4256)
)

# code line for shift correction
metadata(prot)$fov_positions$y_global_px <- metadata(prot)$fov_positions$y_global_px - 4256

## ----load-xenium, message = TRUE, warning = TRUE------------------------------
xeniumFolder <- system.file( "extdata", "Xenium_small", package="SpaceTrooper")

xen <- readXeniumSPE(dirName=xeniumFolder, 
                     sampleName="Xenium_small",
                     type=c("HDF5", "sparse"),
                     coordNames=c("x_centroid", "y_centroid"),
                     boundariesType=c("parquet", "csv"),
                     computeMissingMetrics=TRUE,
                     keepPolygons=FALSE,
                     countsFilePattern="cell_feature_matrix",
                     metadataFPattern="cells.csv.gz",
                     polygonsFPattern="cell_boundaries",
                     polygonsCol="polygons",
                     txPattern="transcripts",
                     addFOVs=FALSE

)

## ----load-merfish, message = TRUE, warning= TRUE------------------------------
merfishFolder <- system.file("extdata", "Merfish_Tiny", package="SpaceTrooper")

mer <- readMerfishSPE(dirName=merfishFolder, 
                      sampleName="Merfish_Tiny",
                      computeMissingMetrics=TRUE,
                      keepPolygons=FALSE,
                      boundariesType=c("parquet"),
                      countmatFPattern="cell_by_gene.csv",
                      metadataFPattern="cell_metadata.csv",
                      polygonsFPattern="cell_boundaries.parquet",
                      coordNames=c("center_x", "center_y"),
                      polygonsCol="polygons")

## ----load-poly,  message = TRUE-----------------------------------------------
# polygon loading
spe <- readAndAddPolygonsToSPE(spe, 
                               polygonsCol="polygons",
                               keepMultiPol=TRUE,
                               boundariesType=c("csv", "HDF5", "parquet")
)

spe$polygons

## ----load-poly-xen,  message = TRUE-------------------------------------------
# xenium polygon loading
xen <- readAndAddPolygonsToSPE(xen, boundariesType=c("parquet"))

## ----load-poly-mer,  message = TRUE-------------------------------------------
# merfish polygon loading
mer <- readAndAddPolygonsToSPE(mer, boundariesType=c("parquet"))

## ----load-poly-prot,  message = TRUE------------------------------------------
# protein polygon loading
prot <- readAndAddPolygonsToSPE(prot, boundariesType=c("csv"))

## ----analysis-qc-1------------------------------------------------------------
spe <- spatialPerCellQC(spe, rmZeros=TRUE,
        negProbList=c("NegPrb", "Negative", "SystemControl"))

xen <- spatialPerCellQC(xen, rmZeros=TRUE)

mer <- spatialPerCellQC(mer, rmZeros=TRUE)

prot <- spatialPerCellQC(prot, rmZeros=TRUE)

## ----analysis-qc-2------------------------------------------------------------
# check added columns to colData
colnames(colData(spe))

## ----compute-QS,  message = TRUE----------------------------------------------
# MERSCOPE cell_id conversion
# spe$cell_id <- as.character(spe$cell_id)

set.seed(1713)

spe <- computeQCScore(spe, verbose=TRUE)

## ----compute-QS-safe-run,  message = TRUE-------------------------------------
# safe run function
safe_run <- function(expr) {
tryCatch(
  list(result = expr, error = NULL),
  error = function(e) list(result = NULL, error = e)
 )
}

out <- safe_run(computeQCScore(spe))
if (!is.null(out$error)) {
message("Failed: ", out$error$message)
colData(spe)$QC_score <- NA
} else {
spe <- out$result
}

## ----compute-outliers-1,  message = TRUE, warning = TRUE----------------------
# detecting outliers for cell area in Xenium dataset with both Medcouple
# and MAD
xen <- computeSpatialOutlier(xen, computeBy="Area_um", method = "both")

## ----compute-outliers-3,  message = TRUE--------------------------------------
spe <- computeSpatialOutlier(spe, computeBy="Mean.CD68", method = "both")

table(spe$Mean.CD68_outlier_mc)

## ----plot-hist-1--------------------------------------------------------------
spe <- computeSpatialOutlier(spe, computeBy="log2SignalDensity", method = "both")

plotMetricHist(spe, metric="log2SignalDensity", 
               useFences="log2SignalDensity_outlier_mc")

## ----plot-zoom-fov-1----------------------------------------------------------
# plotting FoVs map and zoom in of selected FoV colored by QS
plotZoomFovsMap(spe, fovs=16, colourBy="QC_score", csize=3, scaleBars=TRUE)

## ----plot-zoom-fov-2----------------------------------------------------------
plotZoomFovsMap(prot, fovs = 344, colourBy="Area_um", csize=3, scaleBarMap=TRUE, 
                scaleBarPol=FALSE)

## ----plot-centroids-----------------------------------------------------------
plotCentroids(spe, colourBy="QC_score", size=3, scaleBar=FALSE) +
    scale_fill_viridis_c(option = "plasma") +
    scale_color_viridis_c(option = "plasma")

## ----plot-polygons------------------------------------------------------------
plotPolygons(prot, colourBy="log2SignalDensity", scaleBar=FALSE) + 
    scale_fill_viridis_c(option="mako")

## ----sessionInfo--------------------------------------------------------------
sessionInfo()