From 3412403e15fafe10b3c59519a23fb288edd8d492 Mon Sep 17 00:00:00 2001 From: Shahin Mohammadi Date: Fri, 27 Jan 2023 14:22:26 -0800 Subject: [PATCH 1/6] Added parametric autocorrelation --- R/RcppExports.R | 8 +++ inst/include/ACTIONet_RcppExports.h | 42 +++++++++++++++ src/ACTIONet | 2 +- src/ACTIONet.cpp | 36 ++++++++++++- src/RcppExports.cpp | 80 +++++++++++++++++++++++++++++ 5 files changed, 166 insertions(+), 2 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 4383a176..1658655e 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1099,6 +1099,14 @@ assess_label_enrichment <- function(G, M, thread_no = 0L) { .Call(`_ACTIONet_assess_label_enrichment`, G, M, thread_no) } +autocorrelation_Moran_parametric_full <- function(G, scores, normalization_method = 4L, thread_no = 0L) { + .Call(`_ACTIONet_autocorrelation_Moran_parametric_full`, G, scores, normalization_method, thread_no) +} + +autocorrelation_Moran_parametric <- function(G, scores, normalization_method = 4L, thread_no = 0L) { + .Call(`_ACTIONet_autocorrelation_Moran_parametric`, G, scores, normalization_method, thread_no) +} + roll_var <- function(X) { .Call(`_ACTIONet_roll_var`, X) } diff --git a/inst/include/ACTIONet_RcppExports.h b/inst/include/ACTIONet_RcppExports.h index 58d81b7a..da021312 100644 --- a/inst/include/ACTIONet_RcppExports.h +++ b/inst/include/ACTIONet_RcppExports.h @@ -2355,6 +2355,48 @@ namespace ACTIONet { return Rcpp::as(rcpp_result_gen); } + inline List autocorrelation_Moran_parametric_full(mat G, mat scores, int normalization_method = 4, int thread_no = 0) { + typedef SEXP(*Ptr_autocorrelation_Moran_parametric_full)(SEXP,SEXP,SEXP,SEXP); + static Ptr_autocorrelation_Moran_parametric_full p_autocorrelation_Moran_parametric_full = NULL; + if (p_autocorrelation_Moran_parametric_full == NULL) { + validateSignature("List(*autocorrelation_Moran_parametric_full)(mat,mat,int,int)"); + p_autocorrelation_Moran_parametric_full = (Ptr_autocorrelation_Moran_parametric_full)R_GetCCallable("ACTIONet", "_ACTIONet_autocorrelation_Moran_parametric_full"); + } + RObject rcpp_result_gen; + { + RNGScope RCPP_rngScope_gen; + rcpp_result_gen = p_autocorrelation_Moran_parametric_full(Shield(Rcpp::wrap(G)), Shield(Rcpp::wrap(scores)), Shield(Rcpp::wrap(normalization_method)), Shield(Rcpp::wrap(thread_no))); + } + if (rcpp_result_gen.inherits("interrupted-error")) + throw Rcpp::internal::InterruptedException(); + if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen)) + throw Rcpp::LongjumpException(rcpp_result_gen); + if (rcpp_result_gen.inherits("try-error")) + throw Rcpp::exception(Rcpp::as(rcpp_result_gen).c_str()); + return Rcpp::as(rcpp_result_gen); + } + + inline List autocorrelation_Moran_parametric(sp_mat G, mat scores, int normalization_method = 4, int thread_no = 0) { + typedef SEXP(*Ptr_autocorrelation_Moran_parametric)(SEXP,SEXP,SEXP,SEXP); + static Ptr_autocorrelation_Moran_parametric p_autocorrelation_Moran_parametric = NULL; + if (p_autocorrelation_Moran_parametric == NULL) { + validateSignature("List(*autocorrelation_Moran_parametric)(sp_mat,mat,int,int)"); + p_autocorrelation_Moran_parametric = (Ptr_autocorrelation_Moran_parametric)R_GetCCallable("ACTIONet", "_ACTIONet_autocorrelation_Moran_parametric"); + } + RObject rcpp_result_gen; + { + RNGScope RCPP_rngScope_gen; + rcpp_result_gen = p_autocorrelation_Moran_parametric(Shield(Rcpp::wrap(G)), Shield(Rcpp::wrap(scores)), Shield(Rcpp::wrap(normalization_method)), Shield(Rcpp::wrap(thread_no))); + } + if (rcpp_result_gen.inherits("interrupted-error")) + throw Rcpp::internal::InterruptedException(); + if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen)) + throw Rcpp::LongjumpException(rcpp_result_gen); + if (rcpp_result_gen.inherits("try-error")) + throw Rcpp::exception(Rcpp::as(rcpp_result_gen).c_str()); + return Rcpp::as(rcpp_result_gen); + } + } #endif // RCPP_ACTIONet_RCPPEXPORTS_H_GEN_ diff --git a/src/ACTIONet b/src/ACTIONet index fb44abb2..f23ae41e 160000 --- a/src/ACTIONet +++ b/src/ACTIONet @@ -1 +1 @@ -Subproject commit fb44abb23d59b04fb21a51c67afa9484e454cd84 +Subproject commit f23ae41efaeaaec113931a0aeeb8a5ecf08d51ea diff --git a/src/ACTIONet.cpp b/src/ACTIONet.cpp index 73caa4ab..594e3252 100644 --- a/src/ACTIONet.cpp +++ b/src/ACTIONet.cpp @@ -2605,4 +2605,38 @@ mat assess_label_enrichment(sp_mat &G, mat &M, int thread_no = 0) mat logPvals = ACTIONet::assess_label_enrichment(G, M, thread_no); return (logPvals); -} \ No newline at end of file +} + +// [[Rcpp::depends(RcppArmadillo)]] +// [[Rcpp::export]] +List autocorrelation_Moran_parametric_full(mat G, mat scores, + int normalization_method = 4, + int thread_no = 0) +{ + field out = ACTIONet::autocorrelation_Moran_parametric(G, scores, normalization_method, thread_no); + + List res; + res["stat"] = out[0]; + res["zscore"] = out[1]; + res["mu"] = out[2]; + res["sigma"] = out[3]; + + return (res); +} + +// [[Rcpp::depends(RcppArmadillo)]] +// [[Rcpp::export]] +List autocorrelation_Moran_parametric(sp_mat G, mat scores, + int normalization_method = 4, + int thread_no = 0) +{ + field out = ACTIONet::autocorrelation_Moran_parametric(G, scores, normalization_method, thread_no); + + List res; + res["stat"] = out[0]; + res["zscore"] = out[1]; + res["mu"] = out[2]; + res["sigma"] = out[3]; + + return (res); +} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 1fff1b1c..bbd0bdbb 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -4192,6 +4192,80 @@ RcppExport SEXP _ACTIONet_assess_label_enrichment(SEXP GSEXP, SEXP MSEXP, SEXP t UNPROTECT(1); return rcpp_result_gen; } +// autocorrelation_Moran_parametric_full +List autocorrelation_Moran_parametric_full(mat G, mat scores, int normalization_method, int thread_no); +static SEXP _ACTIONet_autocorrelation_Moran_parametric_full_try(SEXP GSEXP, SEXP scoresSEXP, SEXP normalization_methodSEXP, SEXP thread_noSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::traits::input_parameter< mat >::type G(GSEXP); + Rcpp::traits::input_parameter< mat >::type scores(scoresSEXP); + Rcpp::traits::input_parameter< int >::type normalization_method(normalization_methodSEXP); + Rcpp::traits::input_parameter< int >::type thread_no(thread_noSEXP); + rcpp_result_gen = Rcpp::wrap(autocorrelation_Moran_parametric_full(G, scores, normalization_method, thread_no)); + return rcpp_result_gen; +END_RCPP_RETURN_ERROR +} +RcppExport SEXP _ACTIONet_autocorrelation_Moran_parametric_full(SEXP GSEXP, SEXP scoresSEXP, SEXP normalization_methodSEXP, SEXP thread_noSEXP) { + SEXP rcpp_result_gen; + { + Rcpp::RNGScope rcpp_rngScope_gen; + rcpp_result_gen = PROTECT(_ACTIONet_autocorrelation_Moran_parametric_full_try(GSEXP, scoresSEXP, normalization_methodSEXP, thread_noSEXP)); + } + Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error"); + if (rcpp_isInterrupt_gen) { + UNPROTECT(1); + Rf_onintr(); + } + bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen); + if (rcpp_isLongjump_gen) { + Rcpp::internal::resumeJump(rcpp_result_gen); + } + Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error"); + if (rcpp_isError_gen) { + SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen); + UNPROTECT(1); + Rf_error(CHAR(rcpp_msgSEXP_gen)); + } + UNPROTECT(1); + return rcpp_result_gen; +} +// autocorrelation_Moran_parametric +List autocorrelation_Moran_parametric(sp_mat G, mat scores, int normalization_method, int thread_no); +static SEXP _ACTIONet_autocorrelation_Moran_parametric_try(SEXP GSEXP, SEXP scoresSEXP, SEXP normalization_methodSEXP, SEXP thread_noSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::traits::input_parameter< sp_mat >::type G(GSEXP); + Rcpp::traits::input_parameter< mat >::type scores(scoresSEXP); + Rcpp::traits::input_parameter< int >::type normalization_method(normalization_methodSEXP); + Rcpp::traits::input_parameter< int >::type thread_no(thread_noSEXP); + rcpp_result_gen = Rcpp::wrap(autocorrelation_Moran_parametric(G, scores, normalization_method, thread_no)); + return rcpp_result_gen; +END_RCPP_RETURN_ERROR +} +RcppExport SEXP _ACTIONet_autocorrelation_Moran_parametric(SEXP GSEXP, SEXP scoresSEXP, SEXP normalization_methodSEXP, SEXP thread_noSEXP) { + SEXP rcpp_result_gen; + { + Rcpp::RNGScope rcpp_rngScope_gen; + rcpp_result_gen = PROTECT(_ACTIONet_autocorrelation_Moran_parametric_try(GSEXP, scoresSEXP, normalization_methodSEXP, thread_noSEXP)); + } + Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error"); + if (rcpp_isInterrupt_gen) { + UNPROTECT(1); + Rf_onintr(); + } + bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen); + if (rcpp_isLongjump_gen) { + Rcpp::internal::resumeJump(rcpp_result_gen); + } + Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error"); + if (rcpp_isError_gen) { + SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen); + UNPROTECT(1); + Rf_error(CHAR(rcpp_msgSEXP_gen)); + } + UNPROTECT(1); + return rcpp_result_gen; +} // roll_var vec roll_var(vec& X); RcppExport SEXP _ACTIONet_roll_var(SEXP XSEXP) { @@ -4333,6 +4407,8 @@ static int _ACTIONet_RcppExport_validate(const char* sig) { signatures.insert("List(*aggregate_genesets)(sp_mat&,sp_mat&,sp_mat&,int,double,int)"); signatures.insert("mat(*run_harmony)(mat&,mat&,vec,int,double,double,int,int,double,double,double,double,double,bool,int)"); signatures.insert("mat(*assess_label_enrichment)(sp_mat&,mat&,int)"); + signatures.insert("List(*autocorrelation_Moran_parametric_full)(mat,mat,int,int)"); + signatures.insert("List(*autocorrelation_Moran_parametric)(sp_mat,mat,int,int)"); } return signatures.find(sig) != signatures.end(); } @@ -4450,6 +4526,8 @@ RcppExport SEXP _ACTIONet_RcppExport_registerCCallable() { R_RegisterCCallable("ACTIONet", "_ACTIONet_aggregate_genesets", (DL_FUNC)_ACTIONet_aggregate_genesets_try); R_RegisterCCallable("ACTIONet", "_ACTIONet_run_harmony", (DL_FUNC)_ACTIONet_run_harmony_try); R_RegisterCCallable("ACTIONet", "_ACTIONet_assess_label_enrichment", (DL_FUNC)_ACTIONet_assess_label_enrichment_try); + R_RegisterCCallable("ACTIONet", "_ACTIONet_autocorrelation_Moran_parametric_full", (DL_FUNC)_ACTIONet_autocorrelation_Moran_parametric_full_try); + R_RegisterCCallable("ACTIONet", "_ACTIONet_autocorrelation_Moran_parametric", (DL_FUNC)_ACTIONet_autocorrelation_Moran_parametric_try); R_RegisterCCallable("ACTIONet", "_ACTIONet_RcppExport_validate", (DL_FUNC)_ACTIONet_RcppExport_validate); return R_NilValue; } @@ -4566,6 +4644,8 @@ static const R_CallMethodDef CallEntries[] = { {"_ACTIONet_aggregate_genesets", (DL_FUNC) &_ACTIONet_aggregate_genesets, 6}, {"_ACTIONet_run_harmony", (DL_FUNC) &_ACTIONet_run_harmony, 15}, {"_ACTIONet_assess_label_enrichment", (DL_FUNC) &_ACTIONet_assess_label_enrichment, 3}, + {"_ACTIONet_autocorrelation_Moran_parametric_full", (DL_FUNC) &_ACTIONet_autocorrelation_Moran_parametric_full, 4}, + {"_ACTIONet_autocorrelation_Moran_parametric", (DL_FUNC) &_ACTIONet_autocorrelation_Moran_parametric, 4}, {"_ACTIONet_roll_var", (DL_FUNC) &_ACTIONet_roll_var, 1}, {"_ACTIONet_computeSparseRowVariances", (DL_FUNC) &_ACTIONet_computeSparseRowVariances, 4}, {"_ACTIONet_RcppExport_registerCCallable", (DL_FUNC) &_ACTIONet_RcppExport_registerCCallable, 0}, From 891a58c3d1db491a0dcc88cbb2af4298203c1ec1 Mon Sep 17 00:00:00 2001 From: Shahin Mohammadi Date: Sun, 29 Jan 2023 22:24:28 -0800 Subject: [PATCH 2/6] Fixed parametric autocorrelation --- src/ACTIONet | 2 +- src/ACTIONet.cpp | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/ACTIONet b/src/ACTIONet index f23ae41e..d2a8e8e5 160000 --- a/src/ACTIONet +++ b/src/ACTIONet @@ -1 +1 @@ -Subproject commit f23ae41efaeaaec113931a0aeeb8a5ecf08d51ea +Subproject commit d2a8e8e52c7af81deb404cec846061e9f41ef099 diff --git a/src/ACTIONet.cpp b/src/ACTIONet.cpp index 594e3252..5aaa79da 100644 --- a/src/ACTIONet.cpp +++ b/src/ACTIONet.cpp @@ -2630,6 +2630,7 @@ List autocorrelation_Moran_parametric(sp_mat G, mat scores, int normalization_method = 4, int thread_no = 0) { + field out = ACTIONet::autocorrelation_Moran_parametric(G, scores, normalization_method, thread_no); List res; From e53a2d241cdc69caeaa44d75a6a2b32a88052053 Mon Sep 17 00:00:00 2001 From: Sebastian Pineda Date: Sat, 18 Feb 2023 23:16:24 -0500 Subject: [PATCH 3/6] Update DESCRIPTION --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index e106e0c5..b51ee87d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: ACTIONet Type: Package Title: Robust multi-resolution analysis of single-cell datasets -Version: 1.0.0 +Version: 3.0.2 Date: 2023-01-01 Authors@R: c( person("Shahin", "Mohammadi", email = "shahin.mohammadi@gmail.com", From 51e9f879a5fee8c2c07419bd14d76f9378119a84 Mon Sep 17 00:00:00 2001 From: Sebastian Pineda Date: Sat, 18 Feb 2023 23:16:54 -0500 Subject: [PATCH 4/6] Update DESCRIPTION --- DESCRIPTION | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index b51ee87d..23300947 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: ACTIONet Type: Package -Title: Robust multi-resolution analysis of single-cell datasets +Title: Multi-resolution analysis of single-cell datasets Version: 3.0.2 -Date: 2023-01-01 +Date: 2022-12-21 Authors@R: c( person("Shahin", "Mohammadi", email = "shahin.mohammadi@gmail.com", role = c("cre", "aut"), comment = c(ORCID = "0000-0001-8734-2326")), From a5bab9cc0663c522a5190e84f823d496d859e3f7 Mon Sep 17 00:00:00 2001 From: Shahin Mohammadi Date: Thu, 20 Apr 2023 19:50:52 -0700 Subject: [PATCH 5/6] Fixed normalization --- R/utils_public.R | 43 ++++++++++++++++++++++++++++++------------- 1 file changed, 30 insertions(+), 13 deletions(-) diff --git a/R/utils_public.R b/R/utils_public.R index 55308809..cf59b136 100644 --- a/R/utils_public.R +++ b/R/utils_public.R @@ -16,20 +16,21 @@ normalize.matrix <- function(S, top_features_frac = 1.0, scale_param = median, transformation = "log", - anchor_features = NULL) { + anchor_features = NULL, + post_rescale = FALSE) { if (!is.matrix(S) && !ACTIONetExperiment:::is.sparseMatrix(S)) { err <- sprintf("`S` must be `matrix` or `sparseMatrix`.\n") stop(err) } if (!is.null(anchor_features)) { - lib_sizes <- Matrix::colMeans(S[anchor_features, ]) + lib_sizes <- Matrix::colSums(S[anchor_features, ]) } else if (top_features_frac < 1.0) { universality <- Matrix::rowMeans(S != 0) selected_features <- which(universality > (1.0 - top_features_frac)) - lib_sizes <- Matrix::colMeans(S[selected_features, ]) + lib_sizes <- Matrix::colSums(S[selected_features, ]) } else { - lib_sizes <- Matrix::colMeans(S) + lib_sizes <- Matrix::colSums(S) } if (is.matrix(S)) { @@ -47,35 +48,51 @@ normalize.matrix <- function(S, S_scaled_norm <- S_scaled * kappa } else if (length(scale_param) == 1) { S_scaled_norm <- S_scaled * scale_param - } else { + } else if (length(scale_param) == ncol(S)) { if (is.matrix(S_scaled)) { - S_scaled_norm <- S_scaled %*% diag(lib_sizes) + S_scaled_norm <- S_scaled %*% diag(x = scale_param) } else { S_scaled_norm <- S_scaled %*% Diagonal(n = length(scale_param), x = scale_param) } + } else { + stop(sprintf("invalid scale_param in the normalize.matrix() function.")) } if (transformation == "log") { - S_scaled_norm <- log1p(S_scaled_norm) + S_scaled_norm_trans <- log1p(S_scaled_norm) } else if (transformation == "tukey") { if (is.matrix(S_scaled_norm)) { - idx <- which(S_scaled_norm > 0) - vv <- S_scaled_norm[idx] + S_scaled_norm_trans <- S_scaled_norm + idx <- which(S_scaled_norm_trans > 0) + vv <- S_scaled_norm_trans[idx] vv_transformed <- sqrt(vv) + sqrt(1 + vv) - S_scaled_norm[idx] <- vv_transformed + S_scaled_norm_trans[idx] <- vv_transformed } else { - S_scaled_norm@x <- sqrt(S_scaled_norm@x) + sqrt(S_scaled_norm@x + 1) + S_scaled_norm_trans@x <- sqrt(S_scaled_norm_trans@x) + sqrt(S_scaled_norm_trans@x + 1) } } else if (transformation == "lsi") { if (is.matrix(S_scaled_norm)) { S_scaled_norm <- as(S_scaled_norm, "dMatrix") } - S_scaled_norm <- ACTIONet::LSI(S_sparse) + S_scaled_norm_trans <- ACTIONet::LSI(S_scaled_norm) } - return(S_scaled_norm) + if (post_rescale == TRUE) { + cs <- Matrix::colSums(S_scaled_norm_trans) + kappa <- median(cs) / cs + + if (is.matrix(S_scaled_norm_trans)) { + S_scaled_norm_trans <- S_scaled_norm_trans %*% diag(x = kappa) + } else { + S_scaled_norm_trans <- S_scaled_norm_trans %*% Diagonal(n = length(kappa), x = kappa) + } + } + + return(S_scaled_norm_trans) } + + .groupedRowSums <- function(S, group_vec) { if (ACTIONetExperiment:::is.sparseMatrix(S)) { mat <- compute_grouped_rowsums(S, sample_assignments = group_vec) From df5b35e5bb20beade528000bf495c7acc7aeb901 Mon Sep 17 00:00:00 2001 From: Shahin Mohammadi <69611001+mohammadi-insitro@users.noreply.github.com> Date: Wed, 3 May 2023 22:17:21 -0700 Subject: [PATCH 6/6] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 696d24e1..e010d58e 100644 --- a/README.md +++ b/README.md @@ -141,7 +141,7 @@ require(ACTIONet) ace = import.ace.from.10X.h5('pbmc_10k_v3.h5', prefilter = T, min_cells_per_feat = 0.01, min_feats_per_cell = 1000) ace = normalize.ace(ace) ace = reduce.ace(ace) -ace = run.ACTIONet(ace) +ace = runACTIONet(ace) # Annotate cell-types data("curatedMarkers_human")