diff --git a/src/metrics/ks_statistic_sc_features/script.R b/src/metrics/ks_statistic_sc_features/script.R index 0ea1e16..7cf142b 100644 --- a/src/metrics/ks_statistic_sc_features/script.R +++ b/src/metrics/ks_statistic_sc_features/script.R @@ -9,11 +9,13 @@ requireNamespace("SingleCellExperiment", quietly = TRUE) requireNamespace("SummarizedExperiment", quietly = TRUE) -packages <- c("SingleCellExperiment", "SummarizedExperiment", "concaveman", "sp", - "Matrix", "methods", "ggplot2", "ggcorrplot", "MuSiC", "fields", - "MCMCpack", "dplyr", "sf", "RANN", "stats", "reshape2", "RColorBrewer", - "scatterpie", "grDevices", "nnls", "pbmcapply", "spatstat", "gtools", - "RcppML", "NMF") +packages <- c("SingleCellExperiment", "SummarizedExperiment", + "concaveman", "sp", "Matrix", "methods", + "ggplot2", "ggcorrplot", "MuSiC", "fields", + "MCMCpack", "dplyr", "sf", "RANN", "stats", + "reshape2", "RColorBrewer", "scatterpie", + "grDevices", "nnls", "pbmcapply", "spatstat", + "gtools", "RcppML", "NMF") # Load each package lapply(packages, library, character.only = TRUE) @@ -22,7 +24,8 @@ lapply(packages, library, character.only = TRUE) par <- list( input_spatial_dataset = "resources_test/spatialsimbench_mobnew/dataset_sp.h5ad", input_singlecell_dataset = "resources_test/spatialsimbench_mobnew/dataset_sc.h5ad", - input_simulated_dataset = "resources_test/spatialsimbench_mobnew/simulated_dataset.h5ad", output = "output.h5ad" + input_simulated_dataset = "resources_test/spatialsimbench_mobnew/simulated_dataset.h5ad", + output = "output.h5ad" ) meta <- list( name = "ks_statistic", @@ -41,7 +44,7 @@ real_prob_matrix <- input_real_sp$obsm[["celltype_proportions"]] colnames(real_prob_matrix) <- paste0("ct", seq_len(ncol(real_prob_matrix))) sim_log_count <- scuttle::normalizeCounts(t(input_simulated_sp$layers[["counts"]])) -real_log_count <- real_log_count[ , colnames(real_log_count) %in% rownames(real_prob_matrix) ] +real_log_count <- real_log_count[ , colnames(real_log_count) %in% rownames(real_prob_matrix)] sim_sce <- scater::logNormCounts(SingleCellExperiment::SingleCellExperiment( @@ -55,11 +58,12 @@ sim_log_count <- SummarizedExperiment::assay(sim_sce, "logcounts") # build cell type deconvolution in simulated data ## error sim_prob_matrix <- CARD_processing(input_real_sp, input_sc) -sim_log_count <- sim_log_count[ , colnames(sim_log_count) %in% rownames(sim_prob_matrix) ] +sim_log_count <- sim_log_count[ , colnames(sim_log_count) %in% rownames(sim_prob_matrix)] feat_types <- c("L_stats", "celltype_interaction", "nn_correlation", "morans_I") -real_scfeatures_result <- scFeatures::scFeatures(real_log_count, +real_scfeatures_result <- scFeatures::scFeatures( + real_log_count, sample = rep("sample1", ncol(real_log_count)), spatialCoords = input_real_sp$obs[c("row", "col")], feature_types = feat_types, @@ -68,15 +72,18 @@ real_scfeatures_result <- scFeatures::scFeatures(real_log_count, spotProbability = t(real_prob_matrix) ) -sim_scfeatures_result <- scFeatures::scFeatures(sim_log_count, - sample = rep("sample1", ncol(sim_log_count)), - spatialCoords = list( as.numeric( unlist( lapply( strsplit( colnames(sim_log_count) , "x"), `[`, 1) ) ) , - as.numeric( unlist( lapply( strsplit( colnames(sim_log_count) , "x"), `[`, 2) ) ) ) , - feature_types = feat_types, - type = "spatial_t", - species = sc_species, - spotProbability = t(sim_prob_matrix) ) - +sim_scfeatures_result <- scFeatures::scFeatures( + sim_log_count, + sample = rep("sample1", ncol(sim_log_count)), + spatialCoords = list( + as.numeric(unlist(lapply(strsplit(colnames(sim_log_count) , "x"), `[`, 1))), + as.numeric(unlist(lapply(strsplit(colnames(sim_log_count) , "x"), `[`, 2))) + ), + feature_types = feat_types, + type = "spatial_t", + species = sc_species, + spotProbability = t(sim_prob_matrix) +) ks_statistic_L_stats <- ks::kde.test(x1 = as.numeric(real_scfeatures_result$L_stats), x2 = as.numeric(sim_scfeatures_result$L_stats)) ks_statistic_celltype_interaction <- ks::kde.test(x1 = as.numeric(real_scfeatures_result$celltype_interaction), x2 = as.numeric(sim_scfeatures_result$celltype_interaction)) @@ -92,10 +99,10 @@ uns_metric_ids <- c( ) uns_metric_values <- c( - ks_statistic_L_stats, - ks_statistic_celltype_interaction, - ks_statistic_nn_correlation, - ks_statistic_morans_I + ks_statistic_L_stats$zstat, + ks_statistic_celltype_interaction$zstat, + ks_statistic_nn_correlation$zstat, + ks_statistic_morans_I$zstat ) cat("Writing output AnnData to file\n")