## ----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 <- computeQScore(spe, verbose=TRUE)

## ----compute-QS-safe-run,  message=TRUE---------------------------------------
# safe run using tryCatch directly so errors are caught before evaluation
out <- tryCatch(
  list(result=computeQScore(spe), error=NULL),
  error=function(e) list(result=NULL, error=e)
)

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

## ----compute-QS-different-formula,  message=TRUE------------------------------

set.seed(1713)

model_formula <- "~ log2Ctrl_total_ratio + log2SignalDensity"

spe <- computeQScore(spe, modelFormula=model_formula, verbose=TRUE)

## ----plot-polygons-dbkero-----------------------------------------------------
plotPolygons(spe, colourBy="QScore") + 
    scale_fill_viridis_c(option="plasma")

## ----apply-QS-model, message=TRUE---------------------------------------------
qs_model <- metadata(spe)$QScore_model
spe <- applyQScoreModel(
    spe,
    qsModel=qs_model,
    scoreName="QScore_transferred"
)

## ----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="QScore", mapPointSize=1, 
                fovNumberSize=7, scaleBars=TRUE)

## ----plot-zoom-fov-2----------------------------------------------------------
plotZoomFovsMap(prot, fovs=344, colourBy="Area_um", mapPointSize=1, 
                fovNumberSize=7, scaleBarMap=TRUE, scaleBarPol=FALSE)

## ----plot-centroids-----------------------------------------------------------
plotCentroids(spe, colourBy="QScore", 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()