diff --git a/CHANGELOG.md b/CHANGELOG.md index 3054e439..f3033975 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -29,6 +29,8 @@ * Added `method/limma_removebatcheffect` component (PR #79). +* Added `methods/scmerge2` component (PR #63). + ## Minor changes * Un-pin the scPRINT version and update parameters (PR #51) diff --git a/src/methods/scmerge2/config.vsh.yaml b/src/methods/scmerge2/config.vsh.yaml new file mode 100644 index 00000000..3e312f0a --- /dev/null +++ b/src/methods/scmerge2/config.vsh.yaml @@ -0,0 +1,59 @@ +__merge__: /src/api/comp_method.yaml +name: scmerge2 +label: scMerge2 +summary: "scMerge2 integrates single-cell RNA-seq datasets by removing unwanted variation estimated from stably expressed genes and pseudo-replicates." +description: | + scMerge2 corrects batch effects while preserving biological signal. It identifies a set of stably + expressed genes (SEGs), which are assumed to remain consistent across datasets, and uses them as + negative controls in a factor analysis model that estimates and removes unwanted variation. + Pseudo-replicates, constructed from pseudo-bulk profiles of cells grouped within each batch, serve + as anchors for the alignment. + + Cell type labels are optional. When they are supplied the pseudo-replicates are built per cell + type; when they are not, scMerge2 identifies the groupings itself with a mutual nearest cluster + procedure. Both modes are exposed through `--cell_type_aware`. +references: + doi: + - 10.1073/pnas.1820006116 +links: + documentation: https://sydneybiox.github.io/scMerge/articles/scMerge2.html + repository: https://github.com/SydneyBioX/scMerge +info: + method_types: [embedding] + preferred_normalization: log_cp10k + variants: + scmerge2_unsupervised: + scmerge2_semisupervised: + cell_type_aware: true +arguments: + - name: --cell_type_aware + type: boolean + default: false + description: | + Build the pseudo-replicates per cell type, using obs['cell_type']. When false, scMerge2 + identifies the cell groupings itself. + - name: --n_control_genes + type: integer + default: 1000 + description: Number of top-ranked stably expressed genes to use as negative controls. + - name: --n_dim + type: integer + default: 50 + description: Number of principal components in the output embedding. +resources: + - type: r_script + path: script.R +engines: + - type: docker + image: openproblems/base_r:1 + setup: + - type: apt + packages: cmake + - type: r + bioc: + - scMerge +runners: + - type: executable + - type: nextflow + directives: + label: [hightime, highmem, midcpu] diff --git a/src/methods/scmerge2/script.R b/src/methods/scmerge2/script.R new file mode 100644 index 00000000..be8cb0a4 --- /dev/null +++ b/src/methods/scmerge2/script.R @@ -0,0 +1,86 @@ +cat("Loading dependencies\n") +requireNamespace("anndata", quietly = TRUE) +library(Matrix, warn.conflicts = FALSE) +requireNamespace("scMerge", quietly = TRUE) +requireNamespace("BiocParallel", quietly = TRUE) +requireNamespace("BiocSingular", quietly = TRUE) +requireNamespace("methods", quietly = TRUE) + +## VIASH START +par <- list( + input = "resources_test/task_batch_integration/cxg_immune_cell_atlas/dataset.h5ad", + output = "output.h5ad", + cell_type_aware = FALSE, + n_control_genes = 1000L, + n_dim = 50L +) +meta <- list( + name = "scmerge2", + cpus = 1L +) +## VIASH END + +n_cpus <- if (is.null(meta$cpus)) 1L else meta$cpus +bpparam <- BiocParallel::MulticoreParam(workers = n_cpus) + +cat("Read input\n") +adata <- anndata::read_h5ad(par$input) + +# both scSEGIndex and scMerge2 want genes in the rows. scMerge2 only coerces a *dense* matrix to +# CsparseMatrix, so a row-compressed one would slip past that check and break later on +exprs_mat <- methods::as(Matrix::t(adata$layers[["normalized"]]), "CsparseMatrix") +rownames(exprs_mat) <- as.character(adata$var_names) +colnames(exprs_mat) <- as.character(adata$obs_names) + +cat("Select stably expressed genes\n") +seg_df <- scMerge::scSEGIndex(exprs_mat = exprs_mat, BPPARAM = bpparam) +seg_df <- seg_df[order(seg_df$segIdx, decreasing = TRUE), , drop = FALSE] +ctl <- rownames(seg_df)[seq_len(min(par$n_control_genes, nrow(seg_df)))] + +cat("Run scMerge2\n") +out <- scMerge::scMerge2( + exprsMat = exprs_mat, + batch = as.character(adata$obs$batch), + cellTypes = if (par$cell_type_aware) as.character(adata$obs$cell_type) else NULL, + ctl = ctl, + use_bpparam = bpparam, + use_bsparam = BiocSingular::RandomParam(), + verbose = TRUE +) + +cat("Compute embedding\n") +newY <- out$newY +stopifnot(ncol(newY) == adata$n_obs) +if (!is.null(colnames(newY))) { + newY <- newY[, adata$obs_names, drop = FALSE] +} +# subtracting the estimated unwanted variation makes the corrected matrix dense, so let +# BiocSingular stream it rather than materialising it in one go +corrected <- t(newY) +n_dim <- min(par$n_dim, min(dim(corrected)) - 1L) +embedding <- BiocSingular::runPCA( + corrected, + rank = n_dim, + center = TRUE, + scale = FALSE, + BSPARAM = BiocSingular::RandomParam(), + BPPARAM = bpparam +)$x +rownames(embedding) <- adata$obs_names + +cat("Store output\n") +output <- anndata::AnnData( + obs = adata$obs[, c()], + var = adata$var[, c()], + obsm = list( + X_emb = embedding + ), + uns = list( + dataset_id = adata$uns[["dataset_id"]], + normalization_id = adata$uns[["normalization_id"]], + method_id = meta$name + ) +) + +cat("Write output to file\n") +zzz <- output$write_h5ad(par$output, compression = "gzip") diff --git a/src/workflows/run_benchmark/config.vsh.yaml b/src/workflows/run_benchmark/config.vsh.yaml index 93daa059..891abab5 100644 --- a/src/workflows/run_benchmark/config.vsh.yaml +++ b/src/workflows/run_benchmark/config.vsh.yaml @@ -106,6 +106,7 @@ dependencies: - name: methods/scanorama_correct - name: methods/scanorama_integrate - name: methods/scanvi + - name: methods/scmerge2 - name: methods/scgpt_finetuned - name: methods/scgpt_zeroshot - name: methods/scimilarity diff --git a/src/workflows/run_benchmark/main.nf b/src/workflows/run_benchmark/main.nf index 2f3381aa..4f564571 100644 --- a/src/workflows/run_benchmark/main.nf +++ b/src/workflows/run_benchmark/main.nf @@ -36,6 +36,7 @@ methods = [ scanorama_correct, scanorama_integrate, scanvi, + scmerge2, scgpt_finetuned.run( args: [model: file("s3://openproblems-work/cache/scGPT_human.zip")] ),