diff --git a/.Rbuildignore b/.Rbuildignore index 19be65c..4f53835 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -14,4 +14,6 @@ pipeline_stores ^tests$ ^dev$ _targets -target_framework \ No newline at end of file +target_framework +^\.positai$ +^\.claude$ diff --git a/.github/workflows/rworkflows.yml b/.github/workflows/rworkflows.yml index 86ee2b3..26785ed 100644 --- a/.github/workflows/rworkflows.yml +++ b/.github/workflows/rworkflows.yml @@ -40,6 +40,12 @@ jobs: cont: ~ rspm: ~ steps: + - name: Prefer source installs on macOS as fallback (Bioc 3.23 mac binary gap) + if: runner.os == 'macOS' + run: | + mkdir -p ~/.R + echo 'options(install.packages.check.source = "no", pkgType = "source")' >> ~/.Rprofile + shell: bash - uses: neurogenomics/rworkflows@master with: run_bioccheck: ${{ false }} diff --git a/.gitignore b/.gitignore index 956ee95..0d9decd 100644 --- a/.gitignore +++ b/.gitignore @@ -38,3 +38,4 @@ fibrosis_data !README.md _target* meta +.positai diff --git a/DESCRIPTION b/DESCRIPTION index a930f22..d905d8c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,9 @@ Package: HPCell Title: Massively-Parallel R Native Pipeline for Single-Cell Analysis -Version: 0.5.0 +Version: 0.6.1 Authors@R: c(person("Stefano", "Mangiola", email = "mangiolastefano@gmail.com", + role = c("aut")), + person("Mengyuan", "Shen", email = "shen.m@wehi.edu.au", role = c("aut", "cre")), person("Jiayi", "Si", email = "si.j@wehi.edu.au", role = c("aut")) @@ -40,7 +42,6 @@ Imports: EnsDb.Hsapiens.v86, scater, SingleR, - celldex, scuttle, scDblFinder, magrittr, @@ -74,6 +75,7 @@ Imports: rhdf5 Suggests: testthat (>= 3.0.0), + celldex, qs, Azimuth, CellChat, diff --git a/NAMESPACE b/NAMESPACE index 77d78e2..8889e79 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -156,8 +156,6 @@ importFrom(SummarizedExperiment,rowData) importFrom(biomaRt,getBM) importFrom(biomaRt,useMart) importFrom(callr,r) -importFrom(celldex,BlueprintEncodeData) -importFrom(celldex,MonacoImmuneData) importFrom(crew,crew_controller_local) importFrom(data.table,":=") importFrom(digest,digest) diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 0000000..9b00e02 --- /dev/null +++ b/NEWS.md @@ -0,0 +1,58 @@ +# HPCell NEWS + +## HPCell 0.6.1 + +### Bug fixes + +* Fixed `transform_utility()` to no longer assume the first assay of the input + object is named `"X"`. The canonical output assay name `"X"` is now set + explicitly via a single `assay_name` variable rather than being inferred from + the input, so inputs whose first assay carries any other name are handled + correctly. + +### Warnings + +* `transform_utility()` now emits a warning when the input SCE's assay is not + named `"X"`, reporting the original name before renaming it to `"X"` for + downstream consistency. + +## HPCell 0.6.0 + +### Bug fixes + +* Fixed Azimuth-based cell-type label transfer (`annotation_label_transfer()`) failing + with SeuratObject >= 5.0.0. Azimuth 0.5.0 still calls the defunct + `GetAssayData(slot = ...)` interface; HPCell now temporarily patches the + Azimuth import environment so `slot` is mapped to `layer` during + `RunAzimuth()`. +* Updated `GetAssayData()` calls from deprecated `slot = "counts"` to + `layer = "counts"` in `empty_droplet_id()`, `empty_droplet_threshold()`, and + `alive_identification()` for Seurat 5 compatibility. +* Fixed `transform_utility()` so the `identity` transform also applies + `limit_max_to_scale()` before transformation. Previously only `expm1` was + pre-scaled, which could leave high-count samples unscaled and cause downstream + failures in the assay transformation pipeline. + +### Improvements + +* `non_batch_variation_removal()` now catches `SCTransform()` failures for edge + cases with very few overdispersion genes, emits a warning, and returns + `NULL` instead of stopping the pipeline. +* `initialise_hpc()` gains a new `default_controller` argument to set the + default `crew` controller for targets that do not specify their own + controller via `tar_resources()`. + +### Documentation + +* Expanded and corrected roxygen documentation for `initialise_hpc()`, including + the new `default_controller` argument. + +### Development + +* Added CellNexus 2025 pipeline scripts under `dev/cellnexus-2025-scripts/` + (steps 1–10 for census download, metadata preparation, HPCell execution, + local-cache splitting, pseudobulk preparation, and metadata unification). +* Updated CellNexus 2024 pipeline scripts, including renamed + `step9_unify_and_update_sce_metadata.R` and new steps for pseudobulk + preparation and CellNexus querying. +* Added `.positai` and `.claude` to `.Rbuildignore` and `.gitignore`. diff --git a/R/functions.R b/R/functions.R index eb2cd93..0cb5a83 100644 --- a/R/functions.R +++ b/R/functions.R @@ -48,7 +48,7 @@ empty_droplet_id <- function(input_read_RNA_assay, # Get counts if (inherits(input_read_RNA_assay, "Seurat")) { - counts <- GetAssayData(input_read_RNA_assay, assay, slot = "counts") + counts <- GetAssayData(input_read_RNA_assay, assay, layer = "counts") } else if (inherits(input_read_RNA_assay, "SingleCellExperiment")) { counts <- assay(input_read_RNA_assay, assay) } @@ -269,7 +269,7 @@ empty_droplet_threshold<- function(input_read_RNA_assay, # Get counts if (inherits(input_read_RNA_assay, "Seurat")) { - counts <- GetAssayData(input_read_RNA_assay, assay, slot = "counts") + counts <- GetAssayData(input_read_RNA_assay, assay, layer = "counts") } else if (inherits(input_read_RNA_assay, "SingleCellExperiment")) { counts <- assay(input_read_RNA_assay, assay) } @@ -312,9 +312,6 @@ empty_droplet_threshold<- function(input_read_RNA_assay, #' #' @return A tibble with cell-type annotation data. #' -#' @importFrom celldex BlueprintEncodeData -#' @importFrom celldex MonacoImmuneData -#' #' @importFrom Seurat CreateAssayObject #' @importFrom Seurat SCTransform #' @importFrom Seurat CreateSeuratObject @@ -390,6 +387,12 @@ annotation_label_transfer <- function(input_read_RNA_assay, colnames(input_read_RNA_assay)[2]= "dummy___" } + if (!requireNamespace("celldex", quietly = TRUE)) + stop( + "Package 'celldex' is required for SingleR-based cell-type annotation. ", + "Install it with: BiocManager::install('celldex')" + ) + #snapshotDate(): 2025-10-29 blueprint <- celldex::BlueprintEncodeData( ensembl = feature_nomenclature == "ensembl" @@ -479,7 +482,7 @@ annotation_label_transfer <- function(input_read_RNA_assay, } if(nrow(data_annotated) <= 30 | is.null(reference_azimuth)){ - + # If too little immune cells return(data_annotated) #saveRDS(output_path) @@ -515,11 +518,40 @@ annotation_label_transfer <- function(input_read_RNA_assay, azimuth_annotation = tryCatch({ + # This is necessary because Azimuth relies on Seurat: https://github.com/satijalab/azimuth/issues/195 + if (!"Seurat" %in% .packages()) { + library(Seurat) + } + if(ncol(input_read_RNA_assay)<200) k.weight = 25 else k.weight = 50 - input_read_RNA_assay |> RenameAssays(assay.name = assay, new.assay.name = "RNA") |> - Azimuth::RunAzimuth(reference = reference_azimuth, assay = "RNA", umap.name = "refUMAP") |> + # input_read_RNA_assay |> RenameAssays(assay.name = assay, new.assay.name = "RNA") |> + # Azimuth::RunAzimuth(reference = reference_azimuth, assay = "RNA", umap.name = "refUMAP") + { + tmp <- input_read_RNA_assay |> RenameAssays(assay.name = assay, new.assay.name = "RNA") + + # Below is solved in Azimuth 0.5.1: https://github.com/satijalab/azimuth/issues/294 + # # Azimuth 0.5.0 calls GetAssayData(slot = ...) in ConvertGeneNames(), + # # which is defunct in SeuratObject >= 5.0.0. The binding is in + # # `imports:Azimuth` (parent.env of Azimuth's namespace), confirmed by + # # inspection. Patch it there so all internal Azimuth calls are + # # intercepted, then restore on exit. + # .az_imp <- parent.env(asNamespace("Azimuth")) + # .orig_gad <- get("GetAssayData", envir = .az_imp, inherits = FALSE) + # .shim_gad <- local({ + # orig <- .orig_gad + # function(object, slot = NULL, layer = NULL, ...) { + # if (!is.null(slot) && is.null(layer)) layer <- slot + # orig(object, layer = layer, ...) + # } + # }) + # try(unlockBinding("GetAssayData", .az_imp), silent = TRUE) + # assign("GetAssayData", .shim_gad, envir = .az_imp) + # on.exit(assign("GetAssayData", .orig_gad, envir = .az_imp), add = TRUE) + + Azimuth::RunAzimuth(tmp, reference = reference_azimuth, assay = "RNA", umap.name = "refUMAP") + } %>% as_tibble() |> dplyr::select(.cell, any_of( c( @@ -634,7 +666,7 @@ alive_identification <- function(input_read_RNA_assay, if (inherits(input_read_RNA_assay, "Seurat")) { - counts <- GetAssayData(input_read_RNA_assay, assay = assay, slot = "counts") + counts <- GetAssayData(input_read_RNA_assay, assay = assay, layer = "counts") if (!any(str_which(colnames(input_read_RNA_assay[[]]), nFeature_name)) || !any(str_which(colnames(input_read_RNA_assay[[]]), nCount_name))) { input_read_RNA_assay[[nFeature_name]] <- @@ -661,21 +693,18 @@ alive_identification <- function(input_read_RNA_assay, # Returns a named vector of IDs # Matches the gene id's row by row and inserts NA when it can't find gene names - if (feature_nomenclature == "symbol") { - location <- mapIds( - EnsDb.Hsapiens.v86, - keys=rownames(input_read_RNA_assay), - column="SEQNAME", - keytype="SYMBOL" - ) - } - + location <- mapIds( + EnsDb.Hsapiens.v86, + keys = rownames(input_read_RNA_assay), + column = "SEQNAME", + keytype = if (feature_nomenclature == "symbol") "SYMBOL" else "GENEID" + ) - which_mito = rownames(input_read_RNA_assay) |> str_which("^MT") + which_mito = which(location == "MT") # mitochondrion = # input_read_RNA_assay |> - # GetAssayData( slot = "counts", assay=assay) |> + # GetAssayData( layer = "counts", assay=assay) |> # # # Join mitochondrion statistics # # Compute per-cell quality control metrics for a count matrix or a SingleCellExperiment @@ -726,8 +755,21 @@ alive_identification <- function(input_read_RNA_assay, as_tibble(rownames = ".cell") %>% dplyr::select(-sum, -detected) - # I HAVE TO DROP UNIQUE, AS SOON AS THE BUG IN SEURAT IS RESOLVED. UNIQUE IS BUG PRONE HERE. - percentage_output = PercentageFeatureSet(input_read_RNA_assay, pattern = "^RPS|^RPL", assay = assay) + if (feature_nomenclature == "symbol") { + percentage_output = PercentageFeatureSet(input_read_RNA_assay, pattern = "^RPS|^RPL", assay = assay) + } else { + # Ensembl IDs: resolve ribo gene IDs from biomart reference + data(ensembl_genes_biomart) + ribosome_ensembl_ids <- ensembl_genes_biomart[ + grep("^(RPL|RPS)", ensembl_genes_biomart$external_gene_name), "ensembl_gene_id" + ] + ribosome_features <- intersect(ribosome_ensembl_ids, rownames(input_read_RNA_assay)) + percentage_output = PercentageFeatureSet( + input_read_RNA_assay, + features = if (length(ribosome_features) > 0) ribosome_features else character(0), + assay = assay + ) + } percentage_output = percentage_output[!duplicated(names(percentage_output))] # Compute ribosome statistics ribosome = @@ -765,7 +807,7 @@ alive_identification <- function(input_read_RNA_assay, ribosome |> # Only retrieve metadata so nesting in the next step won't break left_join(input_read_RNA_assay |> as_tibble() |> select(.cell, all_of(cell_type_column)), by = ".cell") |> as_tibble() - + } @@ -808,12 +850,13 @@ alive_identification <- function(input_read_RNA_assay, mitochondrion |> left_join(ribosome) |> mutate(alive = !high_mitochondrion) |> # & !high_ribosome ) |> - # Select informative columns + # Select informative columns select(.cell, {{cell_type_column}}, contains("subsets"), contains("observation"), contains("high"), alive) } + #' Doublet Identification #' #' @description @@ -1105,20 +1148,25 @@ non_batch_variation_removal <- function(input_read_RNA_assay, min_cells <- if (bigm_lgl) 5L else 0L new_min_cells = calculate_num_genes_express_in_cells(m, min_cells) # update min_cells if needed - input_read_RNA_assay <- + input_read_RNA_assay <- tryCatch( input_read_RNA_assay |> - SCTransform( - assay = assay, - return.only.var.genes = FALSE, - residual.features = NULL, - vars.to.regress = factors_to_regress, - vst.flavor = "v2", - scale_factor = 2186, - conserve.memory = TRUE, - min_cells = new_min_cells - ) - # |> - # GetAssayData(assay="SCT") + SCTransform( + assay = assay, + return.only.var.genes = FALSE, + residual.features = NULL, + vars.to.regress = factors_to_regress, + vst.flavor = "v2", + scale_factor = 2186, + conserve.memory = TRUE, + min_cells = new_min_cells + ), + error = function(e) { + warning("HPCell says: post transformed distribution of the sample introduced few overdispersion genes, return NULL because these are very few edge cases.") + NULL + } + ) + + if (is.null(input_read_RNA_assay)) return(NULL) sct_mat <- input_read_RNA_assay |> GetAssayData(assay="SCT") diff --git a/R/modules_grammar_hpc.R b/R/modules_grammar_hpc.R index 4355bbf..5a0a90e 100644 --- a/R/modules_grammar_hpc.R +++ b/R/modules_grammar_hpc.R @@ -1,39 +1,39 @@ #' Initialise an HPCell Targets Pipeline #' #' @description -#' Sets up and writes a `targets` pipeline script for HPCell. Saves input data -#' and configuration to disk, then returns an `HPCell` object that downstream -#' grammar functions (e.g. `remove_empty_DropletUtils`, `evaluate_hpc`) can -#' extend before the pipeline is executed with `evaluate_hpc()`. -#' -#' @param input_hpc Named character vector of paths to input data files, one -#' element per sample. If names are not set, integer indices are used. -#' @param store Directory path where pipeline files and targets store are written. +#' This function sets up and executes a `targets` pipeline for HPCell. It saves input data and configurations, +#' writes a pipeline script, and runs the pipeline using the 'targets' package. +#' +#' @param input_hpc Character vector of input data path for the pipeline. +#' @param store Directory path for storing the pipeline files. #' @param computing_resources A `crew` controller object (or list of controllers) #' specifying the computing back-end. Defaults to a local single-worker controller. +#' @param default_controller Optional character name of the default `crew` +#' controller to use for targets that do not specify their own controller. +#' Passed to `targets::tar_resources(crew = tar_resources_crew(controller = ...))`. +#' `NULL` uses targets defaults. #' @param tier Integer vector (same length as `input_hpc`) assigning each sample #' to a processing tier for tiered execution. Default: all samples in tier 1. -#' @param debug_step Character name of a single target to debug; passed to -#' `targets::tar_option_set(debug = ...)`. `NULL` disables debugging. -#' @param RNA_assay_name Name of the RNA assay in the input Seurat/SCE object. -#' @param gene_nomenclature Character scalar indicating gene identifier type in -#' the input data. One of `"symbol"` or `"ensembl"`. -#' @param data_container_type Character scalar specifying the input data format. -#' Accepted values: `"sce_rds"` (SingleCellExperiment RDS), -#' `"seurat_rds"` (Seurat RDS), `"sce_hdf5"` (HDF5-backed SCE), -#' `"seurat_h5"` (HDF5-backed Seurat). +#' @param debug_step Optional step for debugging. +#' @param RNA_assay_name Name of the RNA assay. +#' @param gene_nomenclature Character vector indicating gene nomenclature in input_data +#' @param data_container_type A character vector of length one specifies the input data type. +#' The accepted input data type are: +#' sce_rds for `SingleCellExperiment` RDS, +#' seurat_rds for `Seurat` RDS, +#' sce_hdf5 for `SingleCellExperiment` HDF5-based object +#' seurat_h5 for `Seurat` HDF5-based object #' @param verbosity Reporter string passed to `targets::tar_make()`. Defaults to #' the current targets configuration value. #' @param error Error-handling strategy passed to `targets::tar_option_set()`. #' `NULL` uses the targets default. #' @param update Cue mode string for `targets::tar_cue()`, controlling when #' targets are re-run. Default: `"thorough"`. -#' @param garbage_collection Numeric interval (in targets) at which R garbage -#' collection is triggered during the pipeline run. Default: `0` (disabled). +#' @param garbage_collection Numeric interval at which R garbage collection is +#' triggered during the pipeline run. Default: `0` (disabled). #' @param workspace_on_error Logical; if `TRUE`, saves a workspace snapshot when #' a target errors. Default: `FALSE`. -#' @return An `HPCell` S3 object containing the initialisation arguments, ready -#' to be extended with pipeline step functions. +#' @return The output of the `targets` pipeline, typically a pre-processed data set. #' #' @importFrom glue glue #' @importFrom targets tar_script @@ -51,6 +51,7 @@ initialise_hpc <- function(input_hpc, store = targets::tar_config_get("store"), computing_resources = crew_controller_local(workers = 1), + default_controller = NULL, tier = rep(1, length(input_hpc)), debug_step = NULL, RNA_assay_name = "RNA", @@ -112,13 +113,14 @@ initialise_hpc <- function(input_hpc, controller = crew_controller_group ( readRDS("temp_computing_resources.rds") ), packages = c("HPCell"), trust_object_timestamps = TRUE, - workspace_on_error = w + workspace_on_error = w, + resources = if (!is.null(dc)) tar_resources(crew = tar_resources_crew(controller = dc)) else tar_resources() ) target_list = list( ) } |> - substitute(env = list(d = debug_step, e = error, u = update, g = garbage_collection, w = workspace_on_error)) |> + substitute(env = list(d = debug_step, e = error, u = update, g = garbage_collection, w = workspace_on_error, dc = default_controller)) |> tar_script_append2(script = glue("{store}.R"), append = FALSE) @@ -621,16 +623,16 @@ calculate_pseudobulk.HPCell = function(input_hpc, group_by = NULL, target_input x = group_by, external_path = glue("{input_hpc$initialisation$store}/external") |> as.character(), container_type = "data_container_type" |> is_target() - ) |> - - # merge - hpc_merge( - target_output = target_output, - user_function = pseudobulk_merge |> quote(), - external_path = glue("{input_hpc$initialisation$store}/external") |> as.character(), - pseudobulk_list = pseudobulk_sample |> is_target(), - packages = c("tidySummarizedExperiment", "HPCell") ) + + # merge: merge step is performed scalably in downstream + # hpc_merge( + # target_output = target_output, + # user_function = pseudobulk_merge |> quote(), + # external_path = glue("{input_hpc$initialisation$store}/external") |> as.character(), + # pseudobulk_list = pseudobulk_sample |> is_target(), + # packages = c("tidySummarizedExperiment", "HPCell") + # ) } diff --git a/R/tranform_assay.R b/R/tranform_assay.R index f10d535..022c4a4 100644 --- a/R/tranform_assay.R +++ b/R/tranform_assay.R @@ -171,14 +171,24 @@ transform_utility = function(input_read_RNA_assay, transform_fx, if(ncol(input_read_RNA_assay) == 0) return(NULL) - # Rename assay names to for consistency - if (length(names(assays(input_read_RNA_assay))) == 1 && - names(assays(input_read_RNA_assay)) != "X") names(assays(input_read_RNA_assay)) <- "X" + # Rename assay to "X" for consistency; warn if the original name differs + if (length(names(assays(input_read_RNA_assay))) == 1 && + names(assays(input_read_RNA_assay)) != "X") { + warning(sprintf( + "Input assay is named '%s', not 'X'. Renaming to 'X' for downstream consistency.", + names(assays(input_read_RNA_assay)) + )) + names(assays(input_read_RNA_assay)) <- "X" + } # strip metadata that we don't need input_read_RNA_assay = input_read_RNA_assay |> - select(.cell, observation_joinid, observation_originalid, donor_id, dataset_id, sample_id, cell_type) + select( + any_of(c(".cell", "observation_joinid", "observation_originalid", + "donor_id", "dataset_id", "sample_id", "cell_type")), + starts_with("cell_type") + ) # Remove reduced dimensions reducedDim(input_read_RNA_assay) = NULL @@ -191,11 +201,11 @@ transform_utility = function(input_read_RNA_assay, transform_fx, dir.create(external_path, showWarnings = FALSE, recursive = TRUE) - # Get the name of the first assay in the data object - assay_name <- names(assays(input_read_RNA_assay))[1] + # Always use "X" as the canonical assay name regardless of the input assay name + assay_name <- "X" - # Extract the counts matrix from the assay - counts <- assay(input_read_RNA_assay, assay_name) + # Extract the counts matrix from the first assay of the input + counts <- assay(input_read_RNA_assay) # Convert transform_method to a function if it is a character string transform_function <- match.fun(transform_fx) @@ -205,11 +215,7 @@ transform_utility = function(input_read_RNA_assay, transform_fx, # Check if the transformation method is not 'identity' and counts exceed counts upper bound # Scale for other transform function is handled internally in `transform_function` if (identical(transform_function, identity) || identical(transform_function, expm1)) { - # For identity: no scaling needed - # For expm1: pre-scale first, then apply - if (identical(transform_function, expm1)) { - counts <- limit_max_to_scale(counts, scale_max) - } + counts <- limit_max_to_scale(counts, scale_max) counts <- transform_function(counts) } else { # identity_with_max_limit and safe_expm1 handle scale_max internally @@ -270,7 +276,7 @@ transform_utility = function(input_read_RNA_assay, transform_fx, # Rebuild the SCE to stay light, and to set the assay with the right name input_read_RNA_assay = SingleCellExperiment( - assays = list(X = input_read_RNA_assay |> assay() ), + assays = setNames(list(input_read_RNA_assay |> assay()), assay_name), colData = colData(input_read_RNA_assay) ) diff --git a/dev/cellnexus-2024-scripts/run_sct_for_cellNexus_test_pipeline_debug_failing_sct.R b/dev/cellnexus-2024-scripts/run_sct_for_cellNexus_test_pipeline_debug_failing_sct.R index f2c6fd3..0a61bc6 100644 --- a/dev/cellnexus-2024-scripts/run_sct_for_cellNexus_test_pipeline_debug_failing_sct.R +++ b/dev/cellnexus-2024-scripts/run_sct_for_cellNexus_test_pipeline_debug_failing_sct.R @@ -7,16 +7,28 @@ library(Matrix) library(purrr) library(glue) -x = sample_names |> head(2) -f = functions |> head(2) -tr = feature_thresh |> head(2) + +samples_to_plot <- sample_target_tbl |> left_join(sample_tbl, by = c("sample_id" = "sample_2")) |> + group_by(dataset_id) |> + slice_head(n=1) |> + ungroup() + +samples_to_plot +samples_to_plot|>saveRDS("~/temp.rds") +samples_to_plot<-readRDS("~/temp.rds") + +x = samples_to_plot |> pull(file_name) |> + set_names(samples_to_plot |> pull(sample_id)) +f = rep("expm1",nrow(samples_to_plot)) +tr = samples_to_plot|>pull(feature_thresh) +ub = rep(10,nrow(samples_to_plot)) job::job({ library(HPCell) x |> initialise_hpc( - store = "~/scratch/test_sct_target_store", + store = "/vast/scratch/users/shen.m//test_sct_target_store", gene_nomenclature = "ensembl", data_container_type = "anndata", # tier = tiers, # WE DON"T NEED AS WE HAVE ELASTIC RESOURCES NOW @@ -29,12 +41,12 @@ job::job({ seconds_idle = 30, crashes_error = 10, options_cluster = crew.cluster::crew_options_slurm( - #memory_gigabytes_required = c(20, 35, 50, 75, 100, 150), + #memory_gigabytes_required = c(20, 35, 50, 75, 100, 150), #memory_gigabytes_required = c(90, 100, 120, 150, 200,240), #memory_gigabytes_required = c(60, 80, 100, 150, 200), - memory_gigabytes_required = c(45, 60, 75, 100, 120, 150), - cpus_per_task = c(2, 2, 5, 10, 20), - time_minutes = c(60*24, 60*24, 60*24, 60*24, 60*24,60*24), + memory_gigabytes_required = c(45, 60, 75, 100, 120, 150), + cpus_per_task = 1, + time_minutes = c(60*4, 60*4, 60*4, 60*4, 60*4,60*4), verbose = T ) ) @@ -42,13 +54,13 @@ job::job({ ), verbosity = "summary", update = "never", - #update = "thorough", + # update = "thorough", error = "continue", garbage_collection = 100, workspace_on_error = TRUE ) |> - transform_assay(fx = f, target_output = "sce_transformed") |> + transform_assay(fx = f, target_output = "sce_transformed", scale_max = ub) |> # # Remove empty outliers based on RNA count threshold per cell remove_empty_threshold(target_input = "sce_transformed", RNA_feature_threshold = tr) |> diff --git a/dev/cellnexus-2024-scripts/samples_reannotate_transformation_for_fail_sct_cellnexus2024.R b/dev/cellnexus-2024-scripts/samples_reannotate_transformation_for_fail_sct_cellnexus2024.R index f31c261..781da59 100644 --- a/dev/cellnexus-2024-scripts/samples_reannotate_transformation_for_fail_sct_cellnexus2024.R +++ b/dev/cellnexus-2024-scripts/samples_reannotate_transformation_for_fail_sct_cellnexus2024.R @@ -152,7 +152,7 @@ job::job({ tar_meta( starts_with("sct_"), store = "~/scratch/test_sct_target_store" ) |> dplyr::count(!is.na(error)) -tar_workspace(sct_matrix_03dd1be4995c69ac,store = "~/scratch/test_sct_target_store" ) +tar_workspace(sct_matrix_bc459084b8534e37,store = "~/scratch/test_sct_target_store" ) sce_transformed_filtered <- sce_transformed |> left_join(empty_tbl, by = ".cell") |> dplyr::filter(!empty_droplet) |> left_join( diff --git a/dev/cellnexus-2024-scripts/step10_query_cellNexus_data_local_and_cloud.R b/dev/cellnexus-2024-scripts/step10_query_cellNexus_data_local_and_cloud.R deleted file mode 100644 index 6829dd4..0000000 --- a/dev/cellnexus-2024-scripts/step10_query_cellNexus_data_local_and_cloud.R +++ /dev/null @@ -1,75 +0,0 @@ -# This scripts test cellNexus API with new data generated in STEP_7_unify_cell_metadata.R -library(dplyr) -library(cellNexus) -library(stringr) -library(zellkonverter) - -cache = "~/scratch/cache_temp" - -x = get_metadata(cache_directory = cache, cloud_metadata = NULL, local_metadata = "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/metadata.2.2.1.parquet") -x = x |> - keep_quality_cells() -x = x |> dplyr::filter( - self_reported_ethnicity == "African" & - assay |> stringr::str_like("%10x%") & - tissue == "lung parenchyma" & - cell_type |> stringr::str_like("%CD4%") - ) - - -# One anndata -anndata = readH5AD("~/scratch/cellNexus/cellxgene/01-07-2024/counts/9e62207287ebeaa020d3e92d17b01f8e___1.h5ad", reader="R",use_hdf5 = T) -anndata - -# Test SCE -sce = x |> - get_single_cell_experiment(cache_directory = "/vast/scratch/users/shen.m/cellNexus", repository = NULL) - -# TEST CPM -cpm = x |> - get_single_cell_experiment(cache_directory = "/vast/scratch/users/shen.m/cellNexus", repository = NULL, assays = "cpm") - -# TEST RANK -rank = x |> - get_single_cell_experiment(cache_directory = "/vast/scratch/users/shen.m/cellNexus", repository = NULL, assays = "rank") - -# TEST SCT -sct = x |> get_single_cell_experiment(cache_directory = "/vast/scratch/users/shen.m/cellNexus", repository = NULL, assays = "sct") - -pseudobulk = x |> - get_pseudobulk(cache_directory = "/vast/scratch/users/shen.m/cellNexus", repository = NULL) - -# TEST metacell_256 -caecum_metacell_256 = get_metadata(cache_directory = cache, - cloud_metadata = NULL, - local_metadata = "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/metadata.2.2.0.parquet") |> - keep_quality_cells() |> - filter(!is.na(metacell_256)) |> - filter(tissue == "caecum epithelium") |> - get_metacell(cache_directory = "/vast/scratch/users/shen.m/cellNexus", cell_aggregation = "metacell_256", repository = NULL) -caecum_metacell_256 - -# Check the number of cells per dataset -x |> dplyr::count(dataset_id) - -# Check the number of cells per sample_id -x |> dplyr::count(sample_id) |> dplyr::count(n>10) - -# Check whether cell_id strategt is implemented -x |> select(cell_id) |> arrange(cell_id) - -# Check QC metrics -x |> dplyr::count(empty_droplet, alive, scDblFinder.class) - -# Check empty droplet ratio -x |> dplyr::count(empty_droplet) |> - collect() |> mutate(n_cells = sum(n), pct = n / n_cells * 100) - -# Check alive ratio -x |> filter(!empty_droplet) |> dplyr::count(alive) |> - collect() |> mutate(n_cells = sum(n), pct = n / n_cells * 100) - -# Check doublet ratio -x |> filter(!empty_droplet, alive) |> dplyr::count(scDblFinder.class) |> - collect() |> mutate(n_cells = sum(n), pct = n / n_cells * 100) - diff --git a/dev/cellnexus-2024-scripts/step1_downloads_census_datasets_from_aws.R b/dev/cellnexus-2024-scripts/step1_downloads_census_datasets_from_aws.R deleted file mode 100644 index 3178e2f..0000000 --- a/dev/cellnexus-2024-scripts/step1_downloads_census_datasets_from_aws.R +++ /dev/null @@ -1,100 +0,0 @@ -# This script identifies CELLxGENE Census datasets stable version from 2024-07-01, and downloads -# their corresponding .h5ad files from AWS S3. -# -# Workflow: -# 1. For every dataset, generate AWS S3 download commands pointing to: -# /vast/scratch/users/shen.m/test_cellnexus_reproducibility/h5ad// -# 2. Write all download commands into a text file. -# 3. Create a GNU Parallel bash script that downloads all new datasets efficiently. -library(cellxgene.census) -library(dplyr) - -# Retrieve previous census stable version dataset_id -census <- open_soma(census_version = "2024-07-01") -# Identify organisms available in the Census -census_data <- census$get("census_data") -org_name <- names(census_data$members)[grepl("homo", names(census_data$members), ignore.case = TRUE)] - -metadata <- census_data$get(org_name)$get("obs") - -selected_columns <- c('assay', 'assay_ontology_term_id','disease', 'disease_ontology_term_id', - 'donor_id', 'sex', 'sex_ontology_term_id', 'self_reported_ethnicity', 'self_reported_ethnicity_ontology_term_id', - 'tissue', 'tissue_ontology_term_id', 'tissue_type', 'development_stage', 'development_stage_ontology_term_id', - 'is_primary_data','dataset_id','observation_joinid', 'suspension_type', - "cell_type", "cell_type_ontology_term_id") - -samples <- metadata$read(column_names = selected_columns, - value_filter = "is_primary_data == 'TRUE'")$concat() - -# Get new datasets -samples <- samples |> as.data.frame() |> distinct() |> - # Add organism - mutate(organism = org_name) - -# saved <- samples |> arrow::write_parquet(glue::glue("/vast/scratch/users/shen.m/test_cellnexus_reproducibility/census_new_datasets_{date}.parquet"), -# compression = "zstd") -# - -# Set the base path where files will be downloaded -date <- "2024-07-01" -path = file.path("/vast/scratch/users/shen.m/test_cellnexus_reproducibility/h5ad/", date) -h5ads.uri = get_census_version_directory() |> tibble::rownames_to_column("version") |> filter(version == date) |> pull(h5ads.uri) - -if (!dir.exists(path)) dir.create(path, recursive = T) - -# Generate the dataset download commands -# Only new datasets will be downloaded to path -dataset_ids_path <- samples |> distinct(dataset_id) |> - #head(2 ) |> - mutate(download_path = paste0("aws s3 cp --no-sign-request ", h5ads.uri, dataset_id, ".h5ad ", path)) |> select(download_path) - -# Path where the script will be saved -output_file_path <- glue::glue( - "/vast/scratch/users/shen.m/test_cellnexus_reproducibility/h5ad/{date}_census_dataset_ids_download_path.txt" -) - -write.table( - dataset_ids_path, - file = output_file_path, - quote = FALSE, - row.names = FALSE, - col.names = FALSE -) - -# Define the content of the bash script -bash_script <- glue::glue(" - -#!/bin/bash - -module load awscli -module load parallel - -# Disable all /dev/tty output from GNU Parallel -export PARALLEL_NO_TTY=1 - -# Config to boost download speed -aws configure set default.s3.max_bandwidth 70MB/s - -# Max 100 parallel HTTPs connections for one file -# Limit AWS CLI internal concurrency PER download (prevents connection storms) -aws configure set default.s3.max_concurrent_requests 8 - -# Path to the file containing AWS S3 copy commands -COMMAND_FILE='{output_file_path}' - -# Number of parallel downloads. Try lower the variable here to avoid TCP limits. -# At most 1-2 simultaneous downloads (prevents firewall congestion) -PARALLEL_DOWNLOADS=2 - -# Execute the download commands in parallel -cat $COMMAND_FILE | parallel -j $PARALLEL_DOWNLOADS --eta --bar --plain --no-notice -") - -# Write the bash script to a file -writeLines(bash_script, "/vast/scratch/users/shen.m/test_cellnexus_reproducibility/parallel_download.sh") - -# Change file permission to make it executable -system("chmod +x /vast/scratch/users/shen.m/test_cellnexus_reproducibility/parallel_download.sh") - -# Execute the bash script -system("/vast/scratch/users/shen.m/test_cellnexus_reproducibility/parallel_download.sh") diff --git a/dev/cellnexus-2024-scripts/step2_cellxgene_to_metadata.R b/dev/cellnexus-2024-scripts/step2_cellxgene_to_metadata.R deleted file mode 100644 index db24929..0000000 --- a/dev/cellnexus-2024-scripts/step2_cellxgene_to_metadata.R +++ /dev/null @@ -1,937 +0,0 @@ -.rs.restartR() -library(tidyverse) -library(targets) -library(glue) - -result_directory = "/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024" - - -tar_script({ - - #-----------------------# - # Input - #-----------------------# - library(tidyverse) - library(targets) - library(tarchetypes) - library(glue) - library(qs) - library(crew) - library(crew.cluster) - - #-----------------------# - # Packages - #-----------------------# - tar_option_set( - packages = c( - "zellkonverter", "cellxgenedp", "CuratedAtlasQueryR", "stringr", "tibble", "tidySingleCellExperiment", "dplyr", "Matrix", - "glue", "qs", "purrr", "tidybulk", "tidySummarizedExperiment", "crew", "magrittr", "digest", "readr", "forcats" - ), - - memory = "transient", - garbage_collection = 100, - storage = "worker", - retrieval = "worker", - error = "continue", - # debug = "dataset_id_sce_b5312463451d7ee3", - cue = tar_cue(mode = "never"), - format = "qs", - - - #-----------------------# - # SLURM - #-----------------------# - - controller = crew_controller_group( - - - crew_controller_slurm( - name = "slurm_1_5", - slurm_memory_gigabytes_per_cpu = 5, - slurm_cpus_per_task = 1, - workers = 200, - tasks_max = 5, - verbose = T, - seconds_idle = 30 - ), - crew_controller_slurm( - name = "slurm_1_10", - slurm_memory_gigabytes_per_cpu = 10, - slurm_cpus_per_task = 1, - workers = 100, - tasks_max = 5, - verbose = T, - seconds_idle = 30 - ), - crew_controller_slurm( - name = "slurm_1_20", - slurm_memory_gigabytes_per_cpu = 20, - slurm_cpus_per_task = 1, - workers = 100, - tasks_max = 5, - verbose = T, , - seconds_idle = 30 - ), - crew_controller_slurm( - name = "slurm_1_40", - slurm_memory_gigabytes_per_cpu = 40, - slurm_cpus_per_task = 1, - workers = 50, - tasks_max = 5, - verbose = T, - seconds_idle = 30 - ), - crew_controller_slurm( - name = "slurm_1_80", - slurm_memory_gigabytes_per_cpu = 80, - slurm_cpus_per_task = 1, - workers = 30, - tasks_max = 5, - verbose = T, - seconds_idle = 30 - ), - crew_controller_slurm( - name = "slurm_1_200", - slurm_memory_gigabytes_per_cpu = 200, - slurm_cpus_per_task = 1, - workers = 5, - tasks_max = 5, - verbose = T, - seconds_idle = 30 - ) - ), - resources = tar_resources(crew = tar_resources_crew("slurm_1_10")) - - ) - - sample_heuristics = function(col_data){ - - col_data |> - - # Sort sample ID - # Fix some sample id missing - when(unique(.$dataset_id) %in% c( - "11b86bc3-6d4d-4e28-903a-0361ea8f6bdf", - "492b0613-ff5b-4fca-a585-503fc4102e4f", - "11b86bc3-6d4d-4e28-903a-0361ea8f6bdf", - "0e8f9ce4-46e5-434e-9ca0-e769d1dd27ea" - ) ~ mutate(., PatientID = glue("{sample} {replicate} {time_point} {target}") |> as.character()) , ~ (.)) %>% - when(unique(.$dataset_id) %in% c( - "0273924c-0387-4f44-98c5-2292dbaab11e", - "a16bec18-5c9f-40ad-8169-12c5199c7506", - "556bb449-bbef-43d3-9487-87031fc0decb" - ) ~ mutate(., PatientID = glue("{Collection.ID} {Genotype} {Location}")|> as.character()) , ~ (.)) %>% - when(unique(.$dataset_id) %in% c( - "b83afdc1-baa1-42c0-bd5b-cb607084757d" - ) ~ mutate(., PatientID = glue("{sex} {development_stage} {disease}")|> as.character()) , ~ (.)) %>% - when(unique(.$dataset_id) %in% c( - "3fe53a40-38ff-4f25-b33b-e4d60f2289ef", - "5c1cc788-2645-45fb-b1d9-2f43d368bba8" - ) ~ mutate(., PatientID = glue("{Batch} {Fetus_id} {Development_day} {sex} {tissue} {disease}")|> as.character()) , ~ (.)) |> - - - mutate_if(is.factor, as.character) |> - type_convert(guess_integer = TRUE) |> - mutate_if(is.integer, as.character) |> - - # Convert types - when("donor_id" %in% colnames(.) ~ mutate(., donor_id = donor_id |> as.character() ), ~(.)) |> - when("Cluster" %in% colnames(.) ~ mutate(., Cluster = Cluster |> as.character() ), ~(.)) |> - when("cluster_id" %in% colnames(.) ~ mutate(., cluster_id = cluster_id |> as.character() ), ~(.)) |> - when("Batch" %in% colnames(.) ~ mutate(., Batch = Batch |> as.character() ), ~(.)) |> - when("batch" %in% colnames(.) ~ mutate(., batch = batch |> as.character() ), ~(.)) |> - when("age" %in% colnames(.) ~ mutate(., age = age |> as.character() ), ~(.)) |> - when("BMI" %in% colnames(.) ~ mutate(., BMI = BMI |> as.character() ), ~(.)) |> - when("donor_BMI" %in% colnames(.) ~ mutate(., donor_BMI = donor_BMI |> as.character() ), ~(.)) |> - when("author_cell_type" %in% colnames(.) ~ mutate(., author_cell_type = author_cell_type |> as.character() ), ~(.)) |> - when("time_point" %in% colnames(.) ~ mutate(., time_point = time_point |> as.character() ), ~(.)) |> - when("cluster" %in% colnames(.) ~ mutate(., cluster = cluster |> as.character() ), ~(.)) |> - when("ClusterID" %in% colnames(.) ~ mutate(., ClusterID = ClusterID |> as.character() ), ~(.)) |> - when("Stage" %in% colnames(.) ~ mutate(., Stage = Stage |> as.character() ), ~(.)) |> - when("individual" %in% colnames(.) ~ mutate(., individual = individual |> as.character() ), ~(.)) |> - when("recurrent_cluster" %in% colnames(.) ~ mutate(., recurrent_cluster = recurrent_cluster |> as.character() ), ~(.)) |> - when("PatientID" %in% colnames(.) ~ mutate(., PatientID = PatientID |> as.character() ), ~(.)) |> - when("PMI" %in% colnames(.) ~ mutate(., PMI = PMI |> as.character() ), ~(.)) |> - when("n_genes" %in% colnames(.) ~ mutate(., n_genes = n_genes |> as.numeric() ), ~(.)) |> - when("n_counts" %in% colnames(.) ~ mutate(., n_counts = n_counts |> as.numeric() ), ~(.)) |> - when("n_genes_by_counts" %in% colnames(.) ~ mutate(., n_genes_by_counts = n_genes_by_counts |> as.numeric() ), ~(.)) |> - when("nUMI" %in% colnames(.) ~ mutate(., nUMI = nUMI |> as.numeric() ), ~(.)) |> - when("percent.cortex" %in% colnames(.) ~ mutate(., percent.cortex = percent.cortex |> as.character() ), ~(.)) |> - when("percent.medulla" %in% colnames(.) ~ mutate(., percent.medulla = percent.medulla |> as.character() ), ~(.)) |> - when("Age" %in% colnames(.) ~ mutate(., Age = Age |> as.numeric() ), ~(.)) |> - when("nCount_RNA" %in% colnames(.) ~ mutate(., nCount_RNA = nCount_RNA |> as.numeric() ), ~(.)) |> - when("is_primary_data" %in% colnames(.) ~ mutate(., is_primary_data = is_primary_data |> as.character() ), ~(.)) |> - - #mutate(across(contains("cluster", ignore.case = TRUE), ~ as.character)) |> - select(-one_of('PCW')) %>% - - # Sort sample ID. It works but not elegant. - # Based on observation of strangely behaving datasets, where sample ID is not clear - when("sampleID" %in% colnames(.) & !"PatientID" %in% colnames(.) ~ - mutate(., PatientID = as.character(sampleID )) |> select(-sampleID), ~(.)) %>% - when("Patient" %in% colnames(.) ~ mutate(., Sample = NA |> as.character()), ~(.)) |> - mutate(sample_placeholder = NA |> as.character()) %>% - when(unique(.$dataset_id)=="e40591e7-0e5a-4bef-9b60-7015abe5b17f" ~ mutate(., sample_placeholder = glue("{batch} {development_stage}") |> as.character()), ~ (.)) %>% - when(unique(.$dataset_id)=="39b6cc45-8c5c-4f7b-944c-58f66da5efb1" ~ mutate(., sample_placeholder =sample_id), ~ (.)) %>% - when(unique(.$dataset_id)=="443d6a0e-dbcb-4002-8af0-628e7d4a18fa" ~ mutate(., sample_placeholder =sample_id), ~ (.)) %>% - when(unique(.$dataset_id)=="a91f075b-52d5-4aa3-8ecc-86c4763a49b3" ~ mutate(., sample_placeholder =sample), ~ (.)) %>% - when(unique(.$dataset_id)=="0af763e1-0e2f-4de6-9563-5abb0ad2b01e" ~ mutate(., sample_placeholder ="only_one_culture"), ~ (.)) %>% - when(unique(.$dataset_id)=="5c64f247-5b7c-4842-b290-65c722a65952" ~ mutate(., sample_placeholder ="only_one_culture"), ~ (.)) %>% - when(unique(.$dataset_id)=="d6f92754-e178-4202-b86f-0f430e965d72" ~ mutate(., sample_placeholder =orig.ident), ~ (.)) %>% - when(unique(.$dataset_id)=="c790ef7a-1523-4627-8603-d6a02f8f4877" ~ mutate(., sample_placeholder =orig.ident), ~ (.)) %>% - when(unique(.$dataset_id)=="1e81a742-e457-4fc6-9c39-c55189ec9dc2" ~ mutate(., sample_placeholder =orig.ident), ~ (.)) %>% - when(unique(.$dataset_id)=="351ef284-b59e-43a5-83ba-0eb907dc282c" ~ mutate(., sample_placeholder =orig.ident), ~ (.)) %>% - when(unique(.$dataset_id)=="f498030e-246c-4376-87e3-90b28c7efb00" ~ mutate(., sample_placeholder =Name), ~ (.)) %>% - - # These are the datasets with too few cells per inferred samples, therefore simplifying - when(unique(.$dataset_id)=="e3a56e00-8417-4d82-9d35-3fab3aac12f2" ~ mutate(., SpecimenID =NA), ~ (.)) %>% - when(unique(.$dataset_id)=="17b34e42-bbd2-494b-bf32-b9229344a3f6" ~ mutate(., Sample =NA), ~ (.)) %>% - - # Fix huge samples for plate experiments - tidyr::extract(.cell, "experiment___", "(^expr?[0-9]+)", remove = F) |> - tidyr::extract(.cell, c("run_from_cell_id"), "(run[a-zA-Z0-9_]+)-.+", remove = FALSE) |> - - mutate(experiment___ = if_else(dataset_id=="3fe53a40-38ff-4f25-b33b-e4d60f2289ef", experiment___, "")) |> - - # If run-based embrio study get sample ID from cell ID - - - # Empirically infer samples from many characteristics - unite("sample_heuristic", one_of( - "sample_placeholder", - "Sample", - "SampleID", - "sample_uuid", - "Sample_ID", - "scRNASeq_sample_ID", - "Sample_Tag", - "Sample.ID", - "sample_names", - "Short_Sample", - "Sample.ID.short", - "Sample.name", - "patient", - "Donor.ID", - "donor_id", - "donor", - "PatientID", - "donor_uuid", - "library_uuid", - "suspension_uuid", - "Patient", - "tissue_section_uuid", - "DonorID", - "specimen", - "SpecimenID", - "Fetus_id", - "individual", - "tissue", - "development_stage", - "assay", - "experiment___", - "disease", - "run_from_cell_id", - "is_primary_data" - ), na.rm = TRUE, sep = "___", remove = F) |> - - - #parquet does not like . prefix - dplyr::rename(cell_ = .cell) |> - - # Add sample hash - mutate(sample_ = getVDigest(algo="md5")(glue("{sample_heuristic}{dataset_id}"))) |> - - # make lighter - mutate_if(is.character, as.factor) |> - - # Some cell ids are dbl - mutate(cell_ = as.factor(cell_)) - - } - - get_metadata = function(.x){ - - cache.path = "/vast/scratch/users/shen.m/cellxgenedp" - dir.create(cache.path, recursive = TRUE, showWarnings = FALSE) - - h5_path = .x |> files_download(dry.run = FALSE, cache.path = cache.path) - - sce = - h5_path |> - readH5AD(use_hdf5 = TRUE, raw = FALSE, skip_assays = TRUE, layers=FALSE, reader = "R" ) - - - if(is.null(sce) || !"donor_id" %in% colnames(colData(sce))) - sce = - h5_path |> - readH5AD(use_hdf5 = TRUE, raw = FALSE, skip_assays = TRUE, layers=FALSE ) - - - - file.remove(h5_path) - - metadata = - sce |> - as_tibble() - - # join the file metadata - column_to_omit_becuse_duplicated = - colnames(.x) |> - intersect(colData(sce) |> colnames()) |> - str_subset("donor_id", negate = TRUE) |> - c("embedding") - - rm(sce) - gc(verbose = FALSE) - - metadata = - metadata |> - left_join( - .x |> - select(!any_of(column_to_omit_becuse_duplicated)) |> - unnest(donor_id) |> - unnest(donor_id) |> - select_if(negate(is.list)) , - by = join_by(donor_id) - ) - - metadata = - metadata |> - sample_heuristics() - - # # delete raw data - # sample_column_to_preserve = - # metadata |> - # slice_sample(n = 500, by = donor_id) |> - # tidybulk::pivot_sample(.sample = sample_) |> - # colnames() - - # # Select only sample_ columns - # metadata = - # metadata |> - # select(sample_, any_of(sample_column_to_preserve)) |> - # distinct() - - - - metadata - } - - select_sample_columns = function(metadata){ - - # delete raw data - sample_column_to_preserve = - metadata |> - slice_sample(n = 500, by = donor_id) |> - distinct(sample_, .keep_all = TRUE) |> - # Drop some clearly cell-wise columns - select(-any_of(c("observation_joinid", "cell_")), -contains("cell_type"), -contains("X_pca"), -contains("X_umap")) |> - #tidybulk::pivot_sample(.sample = sample_) |> - colnames() - - # Select only sample_ columns - metadata |> - select(sample_, any_of(sample_column_to_preserve)) |> - distinct() - - } - #-----------------------# - # Pipeline - #-----------------------# - list( - - # Get rownames - tar_target( - my_db, - db(overwrite=TRUE), - resources = tar_resources(crew = tar_resources_crew("slurm_1_20")) - ), - - # Get SCE SMALL - tar_target( - files_dataset_id, - # cellNexus 2024 uses schema_version 5.1.0 - datasets(my_db) |> - left_join( - files(my_db) |> filter(filetype=="H5AD"), - by = "dataset_id" - ) |> - group_split(dataset_id), - iteration = "list", - resources = tar_resources(crew = tar_resources_crew("slurm_1_80")) - ), - - # Get SCE SMALL - tar_target( - metadata_dataset_id, - get_metadata(files_dataset_id), - pattern = map(files_dataset_id), - iteration = "list", - resources = tar_resources(crew = tar_resources_crew("slurm_1_20")) - ), - - # Get dataset_id cell index dictionary - tar_target( - dataset_cell_dict, - metadata_dataset_id |> - mutate(cell_index = row_number(), - new_cell_id = paste(dataset_id, cell_index, sep = "___"), - cell_id_in_sample = paste(cell_, dataset_id, sep = "___")) |> - select(cell_, dataset_id, new_cell_id, cell_id_in_sample), - pattern = map(metadata_dataset_id), - iteration = "list", - resources = tar_resources(crew = tar_resources_crew("slurm_1_20")) - ), - - # select column that are present in half of the datasets at least, so the common column - tar_target( - common_columns, - metadata_dataset_id |> - map_dfr(~ .x |> colnames() |> as_tibble()) |> - dplyr::count(value) |> - mutate(n_datasets = length(metadata_dataset_id)) |> - filter(n > (n_datasets / 2)) |> - pull(value) , - resources = tar_resources(crew = tar_resources_crew("slurm_1_200")) - ), - - tar_target( - metadata_dataset_id_common_sample_columns, - metadata_dataset_id |> - - # Only get primary data - # filter(is_primary_data=="TRUE") |> - - mutate(cell_ = as.character(cell_)) |> - select(any_of(common_columns)) |> - - # Drop some clearly cell-wise columns - select(-any_of(c("observation_joinid", "cell_")), -contains("cell_type")) |> - - select_sample_columns(), - pattern = map(metadata_dataset_id), - resources = tar_resources(crew = tar_resources_crew("slurm_1_80") ) - ), - - tar_target( - metadata_dataset_id_cell_to_sample_mapping, - metadata_dataset_id |> - - # Only get primary data - # filter(is_primary_data=="TRUE") |> - - mutate( - cell_ = as.character(cell_), - observation_joinid = as.character(observation_joinid) - ) |> - # select(cell_, observation_joinid, sample_, donor_id), - select(observation_joinid, cell_, sample_, donor_id, dataset_id, is_primary_data, sample_heuristic, cell_type, cell_type_ontology_term_id), - pattern = map(metadata_dataset_id), - resources = tar_resources(crew = tar_resources_crew("slurm_1_80")) - ) - - ) - - -}, -ask = FALSE, -script = glue("{result_directory}/_targets.R") -) - -job::job({ - - tar_make( - # callr_function = NULL, - reporter = "summary", - script = glue("{result_directory}/_targets.R"), - store = glue("{result_directory}/_targets") - ) - -}) - - -library(arrow) -library(dplyr) -library(duckdb) - -# Sample metadata -saved <- tar_read(metadata_dataset_id_common_sample_columns, store = glue("{result_directory}/_targets")) |> - write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/sample_metadata.parquet", compression = "zstd") - -# Sample to cell link -saved <- tar_read(metadata_dataset_id_cell_to_sample_mapping, store = glue("{result_directory}/_targets")) |> - write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/cell_ids_for_metadata.parquet", compression = "zstd") - -# Cell id dictionary -tar_read(dataset_cell_dict, store = glue("{result_directory}/_targets"))|>bind_rows() |> - write_parquet("/vast/scratch/users/shen.m/test_cellnexus_reproducibility/dataset_cell_dict.parquet", - compression = "zstd") - - -get_tissue_grouped = function(tissue){ - - list( - - # Respiratory System - "respiratory system" = c( - "lung", "lung parenchyma", "alveolus of lung", "bronchus", - "respiratory airway", "pleura", "pleural effusion", "middle lobe of right lung", - "upper lobe of left lung", "lower lobe of left lung", "upper lobe of right lung", - "lower lobe of right lung", "lingula of left lung", "right lung", "left lung" - ), - - trachea = c( "epithelium of trachea", "trachea"), - - # Cardiovascular System - "cardiovascular system" = c( - "heart", "heart left ventricle", "heart right ventricle", "cardiac ventricle", - "cardiac atrium", "right cardiac atrium", "left cardiac atrium", "apex of heart", - "aorta", "coronary artery", - "venous blood", "anterior wall of left ventricle", "myocardium", "interventricular septum", "ventricular tissue", "basal zone of heart" - ), - - vasculature = c("kidney blood vessel", "artery", "vein", "vasculature", "mesenteric artery"), - # Umbilical Cord Blood - "umbilical cord blood" = "umbilical cord blood", - - # Oesophagus - "oesophagus" = c( - "esophagus", "lower esophagus", "esophagus muscularis mucosa", - "submucosal esophageal gland" - ), - - # Stomach - "stomach" = c( - "stomach", "body of stomach", "cardia of stomach" - ), - - # Small Intestine - "small intestine" = c( - "small intestine", "duodenum", "jejunum", "ileum" - ), - - # Large Intestine - "large intestine" = c( - "large intestine", "colon", "left colon", "right colon", - "sigmoid colon", "descending colon", "transverse colon", - "ascending colon", "hepatic flexure of colon", "caecum", - "rectum", "appendix", "vermiform appendix" - ), - - # Digestive System (General) - "digestive system (general)" = c( - "intestine", "hindgut" - ), - - # Nasal, Oral, and Pharyngeal Regions - "nasal, oral, and pharyngeal regions" = c( - "nasal cavity", "nasopharynx", "oral mucosa", "tongue", "anterior part of tongue", - "posterior part of tongue", "gingiva", "nose", "saliva" - ), - - # Cerebral Lobes and Cortical Areas - "cerebral lobes and cortical areas" = c( - "frontal lobe", "left frontal lobe", "right frontal lobe", "primary motor cortex", - "dorsolateral prefrontal cortex", "superior frontal gyrus", "orbitofrontal cortex", - "medial orbital frontal cortex", "Broca's area", "prefrontal cortex", - "temporal lobe", "left temporal lobe", "right temporal lobe", - "angular gyrus", "entorhinal cortex", - "parietal lobe", "left parietal lobe", "right parietal lobe", "primary somatosensory cortex", - "occipital lobe", "right occipital lobe", "primary visual cortex", - "occipital cortex", "insular cortex", "parietal cortex", "temporal cortex", - "frontal cortex", "Brodmann (1909) area 4", "temporoparietal junction", - "middle temporal gyrus", "cingulate cortex", "brain", "brain white matter", "cerebral cortex", "cerebral nuclei" - ), - - # Limbic and Basal Systems - "limbic and basal systems" = c( - "anterior cingulate cortex", "anterior cingulate gyrus", "hippocampal formation", - "hypothalamus", "thalamic complex", "dentate nucleus", "basal ganglion", - "caudate nucleus", "putamen", "substantia nigra pars compacta", - "lateral ganglionic eminence", "medial ganglionic eminence", - "caudal ganglionic eminence", "ganglionic eminence" - ), - - # Brainstem and Cerebellar Structures - "brainstem and cerebellar structures" = c( - "pons", "midbrain", "myelencephalon", "telencephalon", "forebrain", - "cerebellum", "cerebellum vermis lobule", "cerebellar cortex", - "hemisphere part of cerebellar posterior lobe", "white matter of cerebellum" - ), - - # General Brain and Major Structures - "general brain and major structures" = c( - "spinal cord", "neural tube", "cervical spinal cord white matter" - ), - - # Muscular System (Skeletal Muscles) - "muscular system (skeletal muscles)" = c( - "rectus abdominis muscle", "gastrocnemius", "muscle of abdomen", "muscle organ", - "muscle tissue", "pelvic diaphragm muscle", "skeletal muscle tissue", "muscle of pelvic diaphragm" - ), - - # Connective Tissue - "connective tissue" = c( - "connective tissue", "tendon of semitendinosus", "vault of skull", "bone spine", - "rib" - ), - - # Adipose Tissue - "adipose tissue" = c( - "adipose tissue", "subcutaneous adipose tissue", "visceral abdominal adipose tissue", - "perirenal fat", "omental fat pad", "subcutaneous abdominal adipose tissue", - "abdominal adipose tissue" - ), - - # Endocrine System - "endocrine system" = c( - "thyroid gland", "adrenal tissue", "adrenal gland", "islet of Langerhans", - "endocrine pancreas", "pineal gland" - ), - - # Lymphatic System - "lymphatic system" = c( - "lymph node", "mesenteric lymph node", "thoracic lymph node", - "cervical lymph node", "bronchopulmonary lymph node", "tonsil", "inguinal lymph node" - ), - - # Integumentary System (Skin) - "integumentary system (skin)" = c( - "skin of abdomen", "skin of forearm", "skin of scalp", "skin of face", "skin of leg", - "skin of chest", "skin of back", "skin of hip", "skin of body", "skin of cheek", - "skin of temple", "skin of shoulder", "skin of external ear", "skin of trunk", - "skin of prepuce of penis", "skin epidermis", "arm skin", "lower leg skin", - "hindlimb skin", "zone of skin", "dermis", "skin of nose", "skin of forehead", - "skin of pes", "axilla" - ), - - # Gastrointestinal Accessory Organs - "gallbladder" = "gallbladder", - - # Gastrointestinal Accessory Organs - "pancreas" = c( "pancreas", "exocrine pancreas" ), - - # Gastrointestinal Accessory Organs - "liver" = c( "liver", "caudate lobe of liver", "hepatic cecum" ), - - # Spleen - "spleen" = "spleen", - - # Thymus - "thymus" = "thymus", - - # Blood - "blood" = "blood", - - # Bone Marrow - "bone marrow" = "bone marrow", - - # Female Reproductive System - "female reproductive system" = c( - "uterus", "myometrium", "fallopian tube", "ampulla of uterine tube", - "fimbria of uterine tube", "uterine cervix", "endometrium", - "decidua", "decidua basalis", "placenta", "yolk sac", "isthmus of fallopian tube" - ), - "ovary" = "ovary", - - # Male Reproductive System - "male reproductive system (other)" = c( - "testis", "gonad" - ), - - # Prostate - "prostate" = c( - "prostate gland", "transition zone of prostate", "peripheral zone of prostate" - ), - - # Renal System - "renal system" = c( - "kidney", "cortex of kidney", "renal medulla", "renal papilla", - "renal pelvis", "ureter", "bladder organ" - ), - - # Miscellaneous Glands - "miscellaneous glands" = c( - "parotid gland", "lacrimal gland", "sublingual gland", "mammary gland", - "chorionic villus" - ), - - # Epithelium and Mucosal Tissues - "epithelium and mucosal tissues" = c( - "epithelium of small intestine", "epithelium of esophagus", "caecum epithelium", - "jejunal epithelium", "ileal epithelium", "colonic epithelium", - "submucosa of ascending colon", "submucosa of ileum", "lamina propria", - "lamina propria of large intestine", "lamina propria of small intestine", - "mucosa", "mucosa of colon", "lamina propria of mucosa of colon" - ), - - # Eye and Visual-Related Structures - "sensory-related structures" = c( - "retina", - "retinal neural layer", - "macula lutea", - "macula lutea proper", - "sclera", - "trabecular meshwork", - "conjunctiva", - "pigment epithelium of eye", - "cornea", - "iris", - "ciliary body", - "peripheral region of retina", - "eye trabecular meshwork", - "perifoveal part of retina", - "choroid plexus", - "lens of camera-type eye", - "corneo-scleral junction", - "fovea centralis", - "eye", - "inner ear", - "vestibular system", - "primary auditory cortex" - ), - - # Digestive Tract Junctions and Connections - "digestive tract junctions and connections" = c( - "esophagogastric junction", "duodeno-jejunal junction", "hepatopancreatic ampulla", - "hepatopancreatic duct", "pyloric antrum" - ), - - # Peritoneal and Abdominal Cavity Structures - "peritoneal and abdominal cavity structures" = c( - "peritoneum", "omentum", "retroperitoneum", "mesentery" - ), - - # Breast - "breast" = c( - "breast", "upper outer quadrant of breast" - ) - ) |> - enframe(name ="tissue_groups") |> - distinct() |> - unnest(value) |> - dplyr::rename(tissue = value) |> - mutate() - - # #check - # distinct_tissue = - # tissue |> - # enframe(name = "tissue") |> - # distinct(tissue) - # - # if(nrow(distinct_tissue) != distinct_tissue |> left_join(tissue_grouped_df, copy = TRUE)) - # - # - # tissue |> - # enframe(name = "tissue") |> - # left_join(tissue_grouped_df) -} - -convert_age_labels_to_days <- function(labels) { - # Initialize vector to store age in days - age_days <- rep(NA, length(labels)) - - # Define the mapping for Carnegie stages - carnegie_stages <- c( - '9' = 20, - '10' = 22, - '11' = 24, - '12' = 26, - '13' = 28, - '14' = 32, - '16' = 37, - '17' = 41, - '18' = 44, - '19' = 46, - '20' = 49, - '21' = 51, - '22' = 53, - '23' = 56 - ) - - # Map words to numbers for decades - word_to_num <- c( - 'first' = 0, - 'second' = 10, - 'third' = 20, - 'fourth' = 30, - 'fifth' = 40, - 'sixth' = 50, - 'seventh' = 60, - 'eighth' = 70, - 'ninth' = 80, - 'tenth' = 90 - ) - - # Map ordinal words to numbers - word_ordinal_to_num <- c( - 'first' = 1, - 'second' = 2, - 'third' = 3, - 'fourth' = 4, - 'fifth' = 5, - 'sixth' = 6, - 'seventh' = 7, - 'eighth' = 8, - 'ninth' = 9, - 'tenth' = 10, - 'eleventh' = 11, - 'twelfth' = 12, - 'thirteenth' = 13, - 'fourteenth' = 14, - 'fifteenth' = 15, - 'sixteenth' = 16, - 'seventeenth' = 17, - 'eighteenth' = 18, - 'nineteenth' = 19, - 'twentieth' = 20, - 'twenty-first' = 21, - 'twenty-second' = 22, - 'twenty-third' = 23 - ) - - # Map stages to approximate ages in days - stage_to_age <- list( - 'newborn human' = 0, - 'infant' = 0.5 * 365, - 'child' = 6 * 365, # Midpoint of 2-12 years - 'adolescent' = 15 * 365, # Midpoint of 12-18 years - 'young adult' = 25 * 365, # Approximate age - 'human early adulthood' = 25 * 365, - 'human middle aged' = 50 * 365, - 'human late adulthood' = 70 * 365, - 'human adult' = 40 * 365, - 'human aged' = 75 * 365, - 'mature' = 40 * 365, - 'immature' = 1 * 365, - 'embryonic human' = 28, # Midpoint of embryonic stage - 'organogenesis' = 28, - 'unknown' = NA - ) - - # Loop over labels - for (i in seq_along(labels)) { - label <- labels[i] - - # Initialize age variable - age <- NA - - # Remove leading and trailing whitespaces - label <- trimws(label) - - # 1. Match "unknown" - if (grepl("^unknown$", label, ignore.case = TRUE)) { - age <- NA - } - # 2. Match "[number]-month-old human stage" - else if (grepl("^(\\d+)-month-old human stage$", label)) { - num <- as.numeric(sub("^(\\d+)-month-old human stage$", "\\1", label)) - age <- num * 30 # Average days in a month - } - # 3. Match "[number]-year-old human stage" - else if (grepl("^(\\d+)-year-old human stage$", label)) { - num <- as.numeric(sub("^(\\d+)-year-old human stage$", "\\1", label)) - age <- num * 365 # Average days in a year - } - # 4. Match "[number]th week post-fertilization human stage" - else if (grepl("^(\\d+)(?:st|nd|rd|th) week post-fertilization human stage$", label)) { - num <- as.numeric(sub("^(\\d+)(?:st|nd|rd|th) week post-fertilization human stage$", "\\1", label)) - age <- num * 7 - } - # 5. Match "Carnegie stage [number]" - else if (grepl("^Carnegie stage (\\d+)$", label)) { - num <- sub("^Carnegie stage (\\d+)$", "\\1", label) - if (num %in% names(carnegie_stages)) { - age <- carnegie_stages[[num]] - } else { - age <- NA - } - } - # 6. Match "[number]-[number] year-old human stage" - else if (grepl("^(\\d+)-(\\d+) year-old human stage$", label)) { - num1 <- as.numeric(sub("^(\\d+)-(\\d+) year-old human stage$", "\\1", label)) - num2 <- as.numeric(sub("^(\\d+)-(\\d+) year-old human stage$", "\\2", label)) - avg_years <- (num1 + num2) / 2 - age <- avg_years * 365 - } - # 7. Match "[number]-[number] year-old child stage" - else if (grepl("^(\\d+)-(\\d+) year-old child stage$", label)) { - num1 <- as.numeric(sub("^(\\d+)-(\\d+) year-old child stage$", "\\1", label)) - num2 <- as.numeric(sub("^(\\d+)-(\\d+) year-old child stage$", "\\2", label)) - avg_years <- (num1 + num2) / 2 - age <- avg_years * 365 - } - # 8. Match "under-1-year-old human stage" - else if (grepl("^under-1-year-old human stage$", label)) { - age <- 0.5 * 365 # Assume 0.5 years - } - # 9. Match "[number]-month-old human stage" (again) - else if (grepl("^(\\d+)-month-old human stage$", label)) { - num <- as.numeric(sub("^(\\d+)-month-old human stage$", "\\1", label)) - age <- num * 30 - } - # 10. Match "[ordinal] LMP month human stage" - else if (grepl("^(\\w+) LMP month human stage$", label)) { - ordinal_word <- tolower(sub("^(\\w+) LMP month human stage$", "\\1", label)) - if (ordinal_word %in% names(word_ordinal_to_num)) { - num <- word_ordinal_to_num[ordinal_word] - age <- num * 30 - } else { - age <- NA - } - } - # 11. Match "[ordinal] decade human stage" - else if (grepl("^(\\w+) decade human stage$", label)) { - decade_word <- tolower(sub("^(\\w+) decade human stage$", "\\1", label)) - if (decade_word %in% names(word_to_num)) { - num1 <- word_to_num[decade_word] - num2 <- num1 + 9 - avg_years <- (num1 + num2) / 2 - age <- avg_years * 365 - } else { - age <- NA - } - } - # 12. Match "80 year-old and over human stage" - else if (grepl("^(\\d+).*year-old and over human stage$", label)) { - num <- as.numeric(sub("^(\\d+).*year-old and over human stage$", "\\1", label)) - age <- num * 365 - } - # 13. Match developmental stages - else if (grepl("^(.*) stage$", label)) { - stage <- tolower(sub("^(.*) stage$", "\\1", label)) - if (stage %in% names(stage_to_age)) { - age <- stage_to_age[[stage]] - } else { - age <- NA - } - } - # 14. Default case - else { - age <- NA - } - - # Assign to age_days vector - age_days[i] <- age - } - - return(age_days |> as.integer()) -} - -age_days_tbl = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql("SELECT * FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/sample_metadata.parquet')") - ) |> - distinct(development_stage) |> - as_tibble() |> - mutate(age_days = convert_age_labels_to_days(development_stage)) - -age_days_tbl |> - write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/age_days.parquet") - -tissues_grouped = get_tissue_grouped() - -tissues_grouped |> - write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/tissue_grouped.parquet") - diff --git a/dev/cellnexus-2024-scripts/step3_split-large_samples.R b/dev/cellnexus-2024-scripts/step3_split-large_samples.R deleted file mode 100644 index b934271..0000000 --- a/dev/cellnexus-2024-scripts/step3_split-large_samples.R +++ /dev/null @@ -1,299 +0,0 @@ -# R Script for Processing Sample Metadata in Cellxgene Data -# This script processes sample metadata related to Cellxgene datasets, focusing on Homo sapiens data. -# It filters datasets based on certain criteria like primary data, accepted assays, and large sample size thresholds. -# Additionally, it modifies cell identifiers and merges this information with related datasets to generate final outputs for further analysis. -# The script employs several R packages like arrow, targets, glue, dplyr, and more for data manipulation and storage operations. - - -library(arrow) -library(targets) -library(glue) -library(dplyr) -library(cellxgene.census) -library(stringr) -library(purrr) -library(duckdb) -result_directory = "/vast/scratch/users/shen.m/test_cellnexus_reproducibility" -# # Sample metadata -# sample_meta <- tar_read(metadata_dataset_id_common_sample_columns, store = glue("{result_directory}/_targets")) -# sample_meta |> arrow::write_parquet("~/scratch/Census/sample_meta.parquet", compression = "zstd") - -# Sample to cell link -# sample_to_cell <- tar_read(metadata_dataset_id_cell_to_sample_mapping, store = glue("{result_directory}/_targets")) -# sample_to_cell_primary <- sample_to_cell |> filter(is_primary_data == TRUE) -# sample_to_cell_primary |> arrow::write_parquet("~/scratch/Census/sample_to_cell_primary.parquet", compression = "zstd") - -sample_meta = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql("SELECT * FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/sample_metadata.parquet')") - ) - -sample_to_cell_primary = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql("SELECT * FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/cell_ids_for_metadata.parquet')") - ) - -sample_to_cell_primary_human <- sample_to_cell_primary |> - left_join(sample_meta |> filter(organism == "Homo sapiens"), - by = c("sample_","dataset_id", "donor_id", "is_primary_data", "sample_heuristic"), - copy = T) |> - select(observation_joinid, cell_, sample_, - dataset_id, donor_id, is_primary_data, sample_heuristic, - organism, tissue, development_stage, assay, collection_id, - sex, self_reported_ethnicity, disease) -gc() - -# accepted_assays from census -# accepted_assays <- read.csv("~/git_control/cellNexus/dev/census_accepted_assays_2024-07-01.csv", header=TRUE) this file was used to publish cellNexus, then Census updated assay acceptance -# url <- "https://raw.githubusercontent.com/chanzuckerberg/cellxgene-census/d44bebd3e112ea41d00aa9b2509e2a606402c07d/docs/census_accepted_assays.csv" -# download.file(url, destfile = "./dev/census_accepted_assays_2025-01-30.csv") -accepted_assays <- read.csv("~/git_control/cellNexus/dev/census_accepted_assays_2024-07-01.csv") -colnames(accepted_assays) <- c("id", "assay") - -sample_to_cell_primary_human_accepted_assay <- sample_to_cell_primary_human |> filter(assay %in% accepted_assays$assay) - -large_samples <- sample_to_cell_primary_human_accepted_assay |> - dplyr::count(sample_, assay, collection_id, dataset_id) |> - mutate(above_threshold = n > 15000) - -large_samples_collection_id <- large_samples |> ungroup() |> - dplyr::count(collection_id) |> arrange(desc(n)) - -# function to discard nucleotide in cell_ --------------------------------- -# cell pattern repeated across samples. -# Decision: use modified_cell and sample_ to split data - -# drop cell ID if cell ID is a series of numbers -# ACGT more than 5, drops -# drop cellID if does not have special cahracter : - _ -remove_nucleotides_and_separators <- function(x) { - # convert integer cell ID or contain numerics surrounded by special characters to NA - x[str_detect(x, "^[0-9:_\\-*]+$")] <- NA - - # drop sequence having a consistent stretch of 5 characters from ACGT - modified <- str_replace_all(x, "[ACGT]{5,}", "") - - #remove nucleotides surrounded by optional separators - modified <- str_replace_all(modified, "[:_-]{2,}", "_") -} - -# List of collection IDs -collection_ids <- large_samples_collection_id |> collect() |> pull(collection_id) - -gc() -process_collection <- function(id) { - filtered_data <- sample_to_cell_primary_human_accepted_assay |> - filter(collection_id == id) |> - collect() |> - select(cell_, sample_) - - #filtered_data <- filtered_data |> mutate(cell_modified = remove_nucleotides_and_separators(cell_)) - filtered_data$cell_modified <- remove_nucleotides_and_separators(filtered_data |> pull(cell_)) - filtered_data -} - -final_result <- map(collection_ids, process_collection, .progress = T) -final_result <- reduce(final_result, union_all) - -# conditional generating sample_2 based on whether number of cells > 15K. -sample_to_cell_primary_human_accepted_assay <- sample_to_cell_primary_human_accepted_assay |> - left_join(large_samples, by = c("sample_", "assay","collection_id","dataset_id")) - -sample_to_cell_primary_human_accepted_assay_sample_2 <- - sample_to_cell_primary_human_accepted_assay |> - left_join(final_result, by = c("cell_","sample_"), copy = TRUE) |> - # manual adjust - mutate( - cell_modified = ifelse(dataset_id == "b2dda353-0c96-42df-8dcd-1ea7429a6feb" & sample_ == "5951a81f1d40153bab5d2b808e384f39", - "s14", - cell_modified), - cell_modified = ifelse(dataset_id == "b2dda353-0c96-42df-8dcd-1ea7429a6feb" & sample_ == "7313173de022921da50c34ea2f87c7af", - "s3", - cell_modified) - ) |> - mutate(sample_2 = if_else(above_threshold, - paste(sample_, cell_modified, sep = "___"), - sample_) - ) -# save result -write_parquet_to_parquet = function(data_tbl, output_parquet, compression = "gzip") { - - # Establish connection to DuckDB in-memory database - con_write <- dbConnect(duckdb::duckdb(), dbdir = ":memory:") - - # Register `data_tbl` within the DuckDB connection (this doesn't load it into memory) - duckdb::duckdb_register(con_write, "data_tbl_view", data_tbl) - - # Use DuckDB's COPY command to write `data_tbl` directly to Parquet with compression - copy_query <- paste0(" - COPY data_tbl_view TO '", output_parquet, "' (FORMAT PARQUET, COMPRESSION '", compression, "'); - ") - - # Execute the COPY command - dbExecute(con_write, copy_query) - - # Unregister the temporary view - duckdb::duckdb_unregister(con_write, "data_tbl_view") - - # Disconnect from the database - dbDisconnect(con_write, shutdown = TRUE) -} - -sample_to_cell_primary_human_accepted_assay_sample_2 |> write_parquet_to_parquet("~/scratch/Census_rerun/sample_to_cell_primary_human_accepted_assay_sample_2_modify.parquet") - -gc() - -# Load Census census_version = "2024-07-01" -census <- open_soma(census_version = "2024-07-01") -metadata <- census$get("census_data")$get("homo_sapiens")$get("obs") -selected_columns <- c('assay', 'disease', 'donor_id', 'sex', 'self_reported_ethnicity', 'tissue', 'development_stage','is_primary_data','dataset_id','observation_joinid', - "cell_type", "cell_type_ontology_term_id") -samples <- metadata$read(column_names = selected_columns, - value_filter = "is_primary_data == 'TRUE'")$concat() -samples <- samples |> as.data.frame() |> distinct() -samples |> write_parquet("~/scratch/Census/census_samples_701.parquet") - -######## READ -#sample_to_cell_primary_human_accepted_assay_sample_2 <- arrow::read_parquet("~/scratch/Census_rerun/sample_to_cell_primary_human_accepted_assay_sample_2.parquet") -samples <- - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql("SELECT * FROM read_parquet('~/scratch/Census/census_samples_701.parquet')") - ) - -sample_to_cell_primary_human_accepted_assay_sample_2 <- tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql("SELECT * FROM read_parquet('~/scratch/Census_rerun/sample_to_cell_primary_human_accepted_assay_sample_2_modify.parquet')") -) - -census_samples_to_download <- samples |> - left_join(sample_to_cell_primary_human_accepted_assay_sample_2, - by = c("observation_joinid", "dataset_id"), - relationship = "many-to-many", - copy=TRUE) |> - # Use annotation from census - select(-donor_id.y, - -is_primary_data.y, - -tissue.y, - -development_stage.y, - -assay.y, - -sex.y, - -self_reported_ethnicity.y, - -disease.y) |> - rename(donor_id = donor_id.x, - is_primary_data = is_primary_data.x, - assay = assay.x, - disease = disease.x, - sex = sex.x, - self_reported_ethnicity = self_reported_ethnicity.x, - tissue = tissue.x, - development_stage = development_stage.x - ) |> - #as_tibble() |> - # remove space in the sample_2, as sample_2 will be regarded as filename - mutate(sample_2 = if_else(str_detect(sample_2, " "), str_replace_all(sample_2, " ",""), sample_2)) - -# For query purpose -census_samples_to_download |> write_parquet_to_parquet("/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/census_samples_to_download_MODIFIED.parquet") - - -# This is important: please make sure observation_joinid and cell_ is unique per sample (sample_2) in census_samples_to_download -census_samples_to_download |> dplyr::count(observation_joinid, sample_2) |> dplyr::count(n) -census_samples_to_download |> dplyr::count(cell_, sample_2) |> dplyr::count(n) - -# light version -census_samples_to_download |> group_by(dataset_id, sample_2) |> - summarise(observation_joinid = list(observation_joinid), .groups = "drop") |> as_tibble() |> mutate(list_length = map_dbl(observation_joinid, length)) |> - arrow::write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/census_samples_to_download_groups_MODIFIED.parquet") - -# Establish a connection to DuckDB in memory -job::job({ - - con <- dbConnect(duckdb::duckdb(), dbdir = ":memory:") - - # Create views for each of the datasets in DuckDB - dbExecute(con, " - CREATE VIEW cell_to_refined_sample_from_Mengyuan AS - SELECT cell_, observation_joinid, dataset_id, sample_2 AS sample_id, cell_type, cell_type_ontology_term_id - FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/census_samples_to_download_MODIFIED.parquet') -") - - dbExecute(con, " - CREATE VIEW cell_ids_for_metadata AS - SELECT cell_, observation_joinid, dataset_id, sample_, donor_id - FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/cell_ids_for_metadata.parquet') -") - - dbExecute(con, " - CREATE VIEW sample_metadata AS - SELECT * - FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/sample_metadata.parquet') -") - - dbExecute(con, " - CREATE VIEW age_days_tbl AS - SELECT development_stage, age_days - FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/age_days.parquet') -") - - dbExecute(con, " - CREATE VIEW tissue_grouped AS - SELECT tissue, tissue_groups - FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/tissue_grouped.parquet') -") - - # Perform optimised joins within DuckDB - copy_query <- " -COPY ( - SELECT - cell_to_refined_sample_from_Mengyuan.cell_, - cell_to_refined_sample_from_Mengyuan.observation_joinid, - cell_to_refined_sample_from_Mengyuan.dataset_id, - cell_to_refined_sample_from_Mengyuan.sample_id, - cell_to_refined_sample_from_Mengyuan.cell_type, - cell_to_refined_sample_from_Mengyuan.cell_type_ontology_term_id, - sample_metadata.*, - age_days_tbl.age_days, - tissue_grouped.tissue_groups - - FROM cell_to_refined_sample_from_Mengyuan - - LEFT JOIN cell_ids_for_metadata - ON cell_ids_for_metadata.cell_ = cell_to_refined_sample_from_Mengyuan.cell_ - AND cell_ids_for_metadata.observation_joinid = cell_to_refined_sample_from_Mengyuan.observation_joinid - AND cell_ids_for_metadata.dataset_id = cell_to_refined_sample_from_Mengyuan.dataset_id - - LEFT JOIN sample_metadata - ON cell_ids_for_metadata.sample_ = sample_metadata.sample_ - AND cell_ids_for_metadata.donor_id = sample_metadata.donor_id - AND cell_ids_for_metadata.dataset_id = sample_metadata.dataset_id - - LEFT JOIN age_days_tbl - ON age_days_tbl.development_stage = sample_metadata.development_stage - - LEFT JOIN tissue_grouped - ON tissue_grouped.tissue = sample_metadata.tissue - -) TO '/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/cell_metadata.parquet' -(FORMAT PARQUET, COMPRESSION 'gzip'); -" - - # Execute the final query to write the result to a Parquet file - dbExecute(con, copy_query) - - # Disconnect from the database - dbDisconnect(con, shutdown = TRUE) - -}) - -# system("~/bin/rclone copy /vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/cell_metadata.parquet box_adelaide:/Mangiola_ImmuneAtlas/reannotation_consensus/") - - -cell_metadata = tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql("SELECT * FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/cell_metadata.parquet')") -) - -cell_metadata |> distinct(sample_id) |>dplyr::count() diff --git a/dev/cellnexus-2024-scripts/step4_split_census_anndata_base_on_sample_id.R b/dev/cellnexus-2024-scripts/step4_split_census_anndata_base_on_sample_id.R deleted file mode 100644 index 09baa02..0000000 --- a/dev/cellnexus-2024-scripts/step4_split_census_anndata_base_on_sample_id.R +++ /dev/null @@ -1,140 +0,0 @@ -library(targets) -library(zellkonverter) -library(tibble) -library(dplyr) -library(SummarizedExperiment) -library(tidybulk) -library(tidySingleCellExperiment) -library(stringr) -library(arrow) - -version <- "2024-07-01" -anndata_path_based_on_dataset_id_to_read <- file.path("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/h5ad/", version) -anndata_path_based_on_sample_id_to_save <- file.path("/vast/scratch/users/shen.m/Census/split_h5ad_based_on_sample_id/", version) -dir.create(anndata_path_based_on_sample_id_to_save, recursive = TRUE) - -files <- list.files(anndata_path_based_on_dataset_id_to_read, pattern = "*h5ad", - full.names = TRUE) - -save_data <- function(data, file_name) { - filename <- paste0(file.path(anndata_path_based_on_sample_id_to_save), "/", file_name, ".h5ad") - if(ncol(assay(data)) == 1) { - - # Duplicate the assay to prevent saving errors due to single-column matrices - my_assay = cbind(assay(data), assay(data)) - # Rename the second column to distinguish it - colnames(my_assay)[2] = paste0("DUMMY", "___", colnames(my_assay)[2]) - - cd = colData(data) - cd = cd |> rbind(cd) - rownames(cd)[2] = paste0("DUMMY", "___", rownames(cd)[2]) - - - - data = SingleCellExperiment(assay = list( X = my_assay ), colData = cd) - } - - - zellkonverter::writeH5AD(data, file = filename, compression = "gzip" ) - return(filename) -} - - -#' This function subsets samples from a dataset based on specified join IDs. -#' It reads a single-cell expression dataset stored in the HDF5 AnnData format, -#' removes unneeded metadata to optimize memory usage, and extracts the subset of cells -#' matching the provided observation join IDs. -subset_samples <- function(dataset_id, observation_joinid, sample_id) { - - - # Construct the file path to the dataset - file_path <- paste0(anndata_path_based_on_dataset_id_to_read, "/", dataset_id, ".h5ad") - - # Read the dataset from HDF5 file using zellkonverter package - sce <- zellkonverter::readH5AD(file_path, use_hdf5 = TRUE, reader = "R") - - # Extract the base name of the file and remove the '.h5ad' extension to get dataset_id - sce$dataset_id = file_path |> basename() |> str_remove("\\.h5ad$") - - # Add original cell identifiers from column names to the dataset - sce$observation_originalid = colnames(sce) - - # Create a new identifier by combining the original cell id with the dataset id - cell_modified = paste(sce$observation_originalid, sce$dataset_id, sep = "___") - - # Update column names with the new combined identifier - colnames(sce) <- cell_modified - - # Clear the metadata to reduce memory overhead - S4Vectors::metadata(sce) <- list() - - # Add sample identifier to the dataset - sce <- sce |> mutate(sample_id = !!sample_id) |> - - # Remove the ".h5ad" extension if there is - mutate(sample_id = stringr::str_replace(sample_id, ".h5ad$", "")) - - # Identify cells that match the observation_joinid - cells_to_subset <- which(colData(sce)$observation_joinid %in% unlist(observation_joinid)) - - # Subset and return the SingleCellExperiment object with only the selected cells - sce[, cells_to_subset] -} - -computing_resources = crew.cluster::crew_controller_slurm( - #slurm_memory_gigabytes_per_cpu = 40, - slurm_memory_gigabytes_per_cpu = 25, - slurm_cpus_per_task = 1, - workers = 100, - verbose = TRUE -) - -tar_option_set( - memory = "transient", - garbage_collection = TRUE, - storage = "worker", - retrieval = "worker", - format = "qs", - #cue = tar_cue(mode = "never"), - cue = tar_cue(mode = "thorough"), - error = "continue", - controller = computing_resources -) - -list( - tar_target( - file_paths, - list.files(anndata_path_based_on_dataset_id_to_read, pattern = "*h5ad", - full.names = TRUE), - deployment = "main" - ), - tar_target( - grouped_observation_joinid_per_sample, - # This should be run - read_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/census_samples_to_download_groups_MODIFIED.parquet") |> - - # Note: dataset_id "99950e99-2758-41d2-b2c9-643edcdf6d82" and "9fcb0b73-c734-40a5-be9c-ace7eea401c9" - # from Census does not contain any meaningful data (no observation_joinid in colData), thus produced - # not meaningful samples (0 cells). They need to be deleted. - - filter(!dataset_id %in% c("99950e99-2758-41d2-b2c9-643edcdf6d82", "9fcb0b73-c734-40a5-be9c-ace7eea401c9" )) - ), - tar_target( - sliced_sce, - subset_samples(grouped_observation_joinid_per_sample$dataset_id, - grouped_observation_joinid_per_sample$observation_joinid, - grouped_observation_joinid_per_sample$sample_2) |> - save_data(file_name = grouped_observation_joinid_per_sample$sample_2), - pattern = map(grouped_observation_joinid_per_sample) - ) -) - -# tar_make(store = glue::glue("~/scratch/Census_final_run/{version}_new/split_h5ad_based_on_sample_id_target_store"), -# script = "~/git_control/HPCell/dev/cellnexus-2024-scripts/step4_split_census_anndata_base_on_sample_id.R", -# reporter = "summary") - -# Debug if needed -# tar_errored(store = "~/scratch/Census_final_run/split_h5ad_based_on_sample_id_target_store/") -# tar_meta(store = "~/scratch/Census_final_run/split_h5ad_based_on_sample_id_target_store/") |> -# filter(!is.na(error)) |> pull(error) - diff --git a/dev/cellnexus-2024-scripts/step5_identify_census_sample_counts_distribution.qmd b/dev/cellnexus-2024-scripts/step5_identify_census_sample_counts_distribution.qmd deleted file mode 100644 index cd77cb3..0000000 --- a/dev/cellnexus-2024-scripts/step5_identify_census_sample_counts_distribution.qmd +++ /dev/null @@ -1,434 +0,0 @@ ---- -title: "identify_census_sample_counts_distribution" -author: "Mengyuan Shen" -date: `Sys.Date()` -format: - html: - toc: true - toc-depth: 2 - code-fold: true ---- - - -# Setup -```{r} -# Get sample counts summary ----------------------------------------------- -library(targets) -library(dplyr) -library(stringr) -library(glue) -library(arrow) -library(tidySingleCellExperiment) -set.seed(12345) -Date <- "2024-07-01" -# Identify raw counts range from new census data -summary_store = glue::glue("/vast/scratch/users/shen.m/{Date}_census_sample_raw_counts_summary_target_store") -tar_script({ - library(dplyr) - library(SummarizedExperiment) - library(zellkonverter) - library(crew) - library(crew.cluster) - library(duckdb) - - - # Helper (optional) to avoid repetition - new_elastic <- function(name, mem_gb, time_min, workers, crashes_max, cpus_per_task = 2, backup = NULL) { - crew_controller_slurm( - name = name, - workers = workers, - crashes_max = crashes_max, - seconds_idle = 30, - options_cluster = crew_options_slurm( - memory_gigabytes_required = mem_gb, - cpus_per_task = cpus_per_task, - time_minutes = time_min - ), - backup = backup - ) - } - - # Small → large, with fallbacks to the next size up - elastic_300 <- new_elastic("elastic_300", 300, 60 * 24, workers = 8, crashes_max = 2) - elastic_160 <- new_elastic("elastic_160", 160, 60 * 24, workers = 8, crashes_max = 2, backup = elastic_300) - elastic_120 <- new_elastic("elastic_120", 120, 60 * 4, workers = 16, crashes_max = 1, cpus_per_task = 8, backup = elastic_160) - elastic_80 <- new_elastic("elastic_80", 80, 60 * 4, workers = 24, crashes_max = 1, cpus_per_task = 8, backup = elastic_120) - elastic_40 <- new_elastic("elastic_40", 40, 60 * 4, workers = 32, crashes_max = 1, cpus_per_task = 8, backup = elastic_80) - elastic_20 <- new_elastic("elastic_20", 20, 60 * 4, workers = 48, crashes_max = 1, cpus_per_task = 8, backup = elastic_40) - elastic_10 <- new_elastic("elastic_10", 10, 60 * 4, workers = 150, crashes_max = 6, cpus_per_task = 8, backup = elastic_20) - - elastic_5_minimal <- new_elastic("elastic_5_minimal", 5, 60 * 4, workers = 300, crashes_max = 6, cpus_per_task = 8, backup = elastic_10) - - - # Group for targets (small → large) - controllers <- crew_controller_group( - elastic_10, elastic_20, elastic_40, elastic_80, elastic_120, elastic_160, elastic_300, elastic_5_minimal - ) - tar_option_set( - memory = "transient", - garbage_collection = 100, - storage = "worker", - retrieval = "worker", - error = "continue", - cue = tar_cue(mode = "never"), - format = "qs", - workspace_on_error = TRUE, - controller = controllers, - trust_object_timestamps = TRUE, - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_5_minimal") - ) - ) - - pos_min_med_ratio = function(x){ - mi = x[x > 0] |> min() - me = median(x[x > 0]) - mi / me - } - - get_sample_summary_stats <- function(files) { - sce = readH5AD(files, reader = "R", use_hdf5 = T) - - if (ncol(sce) == 0) return(NULL) - - assay_name = sce@assays |> names() |> magrittr::extract(1) - counts_mat = sce |> assay(assay_name) |> as.matrix() - counts_vec = as.numeric(counts_mat) - sample_id = basename(files) - - # Perform checks - min_val = min(counts_vec, na.rm = TRUE) - max_val = max(counts_vec, na.rm = TRUE) - median_val = median(counts_vec, na.rm = TRUE) - - gap = pos_min_med_ratio(counts_vec) - - has_negative = min_val < 0 - max_gt_10 = max_val > 10 - - # Check if all integers - tol = 1e-4 - all_integer = all(counts_vec == floor(counts_vec), na.rm = TRUE) - - has_floating = !all_integer && all(abs(counts_vec - round(counts_vec)) < tol, na.rm = TRUE) - - tbl = tibble::tibble( - sample_id = sample_id, - min_val = min_val, - median_val = median_val, - max_val = max_val, - counts_gap = gap, - has_negative = has_negative, - max_gt_10 = max_gt_10, - all_integer = all_integer, - has_floating = has_floating, - n_cells = ncol(counts_mat), - n_genes = nrow(counts_mat) - ) - - tbl = tbl |> left_join( - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql("SELECT * FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/census_samples_to_download_groups_MODIFIED.parquet')") - ) |> select(-observation_joinid, list_length) |> - distinct() |> mutate(sample_id = paste0(sample_2, ".h5ad")) |> collect(), - by = "sample_id", - copy = T - ) - - tbl - } - - - list( - tar_target( - files, - list.files("/vast/scratch/users/shen.m/Census/split_h5ad_based_on_sample_id/2024-07-01", full.names = T, pattern = ".h5ad$"), - deployment = "main" - ), - - # Get raw counts matrix with sample_id - tar_target( - sample_summary_df, - get_sample_summary_stats(files), - pattern = map(files), - iteration = "list", - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_20") - ) - ) - ) -}, ask = FALSE, script = glue("{summary_store}/_targets.R")) - - -job::job({ - - tar_make( - # callr_function = NULL, - reporter = "summary", - script = glue("{summary_store}/_targets.R"), - store = glue("{summary_store}/_targets") - ) - -}) - -``` - -```{r} -plot_raw_counts_hist_pdf <- function( - samples_to_plot, - h5ad_dir, - pdf_file, - width = 8, - height = 8, - nrow = 3, - ncol = 3, - max_cells = 5e3, - transform = c("log1p", "none"), - hist_ylim = c(0, 1e5) -) { - - transform <- match.arg(transform) - - pdf(pdf_file, width = width, height = height) - on.exit(dev.off(), add = TRUE) - - par(mfrow = c(nrow, ncol), mar = c(3, 3, 2, 1)) - - samples_to_plot |> - dplyr::pull(sample_id) |> - purrr::map( - ~ { - sce <- zellkonverter::readH5AD( - file = file.path(h5ad_dir, .x), - reader = "R", - use_hdf5 = T - ) - - if (ncol(sce) == 0) return(NULL) - - if (ncol(sce) > max_cells) { - sce <- dplyr::sample_n(sce, max_cells) - } - - assay_name <- names(SummarizedExperiment::assays(sce))[1] - - dataset_id <- sce |> - dplyr::distinct(dataset_id) |> - dplyr::pull() - - x <- SummarizedExperiment::assay(sce, assay_name) |> - as.numeric() - - if (transform == "log1p") { - x <- log1p(x) - } - - hist( - x, - main = dataset_id, - ylim = hist_ylim - ) - - invisible(NULL) - }, - .progress = TRUE - ) - - par(mfrow = c(1, 1)) -} - -``` - -```{r} -sample_summary_df = tar_read(sample_summary_df, store = glue("{summary_store}/_targets")) |> bind_rows() |> - mutate(max_gt_20 = ifelse(max_val > 20, TRUE, FALSE)) - -sample_summary_df |> write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/sample_distribution_summary.parquet") -``` - -```{r echo=FALSE} -sample_summary_df = read_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/sample_distribution_summary.parquet") -``` - -# Sample statstics summary -```{r} -sample_summary_df |> dplyr::count(has_negative, max_gt_20, all_integer, has_floating, counts_gap>0.25) -``` - -# Plot sample in dataset -## Row one -```{r} -sample_summary_df |> dplyr::count(has_negative, max_gt_20, all_integer, has_floating, counts_gap>=0.25) |> dplyr::slice(1) - -samples_to_plot <- sample_summary_df |> filter(!all_integer, !has_negative, !max_gt_20, !has_floating, counts_gap<0.25) |> - arrange(n_cells) |> - group_by(dataset_id) |> - slice_sample(n = 1) |> - ungroup() - -``` - -```{r eval=FALSE} -plot_raw_counts_hist_pdf(samples_to_plot, - h5ad_dir = "/vast/scratch/users/shen.m/Census/split_h5ad_based_on_sample_id/2024-07-01/", - pdf_file = "/home/users/allstaff/shen.m/projects/cellNexus/samples_pos_lt20_decimal_raw_counts_hist.pdf", - transform = "none") - -``` - - -## Row two -```{r} -sample_summary_df |> dplyr::count(has_negative, max_gt_20, all_integer, has_floating, counts_gap>=0.25) |> dplyr::slice(2) -samples_to_plot <- sample_summary_df |> filter( !has_negative, !max_gt_20, !all_integer, !has_floating, counts_gap>=0.25) |> - arrange(n_cells) |> - group_by(dataset_id) |> - slice_sample(n = 1) |> - ungroup() -``` - -```{r eval=FALSE} -plot_raw_counts_hist_pdf(samples_to_plot, - h5ad_dir = "/vast/scratch/users/shen.m/Census/split_h5ad_based_on_sample_id/2024-07-01/", - pdf_file = "/home/users/allstaff/shen.m/projects/cellNexus/plot2.pdf", - transform = "none") -``` - -## Row three -```{r} -sample_summary_df |> dplyr::count(has_negative, max_gt_20, all_integer, has_floating,counts_gap>=0.25) |> dplyr::slice(3) -samples_to_plot <- sample_summary_df |> filter(!has_negative, max_gt_20, !all_integer, !has_floating, counts_gap<0.25) |> - arrange(n_cells) |> - group_by(dataset_id) |> - slice_sample(n = 3) |> - ungroup() -``` - -```{r eval=FALSE} -plot_raw_counts_hist_pdf(samples_to_plot, - h5ad_dir = "/vast/scratch/users/shen.m/Census/split_h5ad_based_on_sample_id/2024-07-01/", - pdf_file = "/home/users/allstaff/shen.m/projects/cellNexus/plot3.pdf", - transform = "none") -``` - -## Row four -```{r} -sample_summary_df |> dplyr::count(has_negative, max_gt_20, all_integer, has_floating,counts_gap>=0.25) |> dplyr::slice(4) -samples_to_plot <- sample_summary_df |> filter(!has_negative, max_gt_20, !all_integer, !has_floating,counts_gap>=0.25) |> - arrange(n_cells) |> - group_by(dataset_id) |> - slice_sample(n = 2) |> - ungroup() -``` - -```{r eval=FALSE} -plot_raw_counts_hist_pdf(samples_to_plot, - h5ad_dir = "/vast/scratch/users/shen.m/Census/split_h5ad_based_on_sample_id/2024-07-01/", - pdf_file = "/home/users/allstaff/shen.m/projects/cellNexus/plot4.pdf", - transform = "none") -``` - -## Row five -```{r} -sample_summary_df |> dplyr::count(has_negative, max_gt_20, all_integer, has_floating, counts_gap>=0.25) |> dplyr::slice(5) -samples_to_plot <- sample_summary_df |> filter(!has_negative,max_gt_20, all_integer, !has_floating, counts_gap<0.25) |> - arrange(n_cells) |> - group_by(dataset_id) |> - slice_sample(n = 2) |> - ungroup() -``` - - -```{r eval=FALSE} -plot_raw_counts_hist_pdf(samples_to_plot, - h5ad_dir = "/vast/scratch/users/shen.m/Census/split_h5ad_based_on_sample_id/2024-07-01/", - pdf_file = "/home/users/allstaff/shen.m/projects/cellNexus/plot5.pdf", - transform = "none") -``` - - -## Row six -```{r} -sample_summary_df |> dplyr::count(has_negative, max_gt_20, all_integer, has_floating, counts_gap>=0.25) |> dplyr::slice(6) -samples_to_plot <- sample_summary_df |> filter(!has_negative, max_gt_20, all_integer, !has_floating, counts_gap>=0.25) |> - arrange(n_cells) |> - group_by(dataset_id) |> - slice_sample(n = 2) |> - ungroup() -``` - - -```{r eval=FALSE} -plot_raw_counts_hist_pdf(samples_to_plot |> head(10), - h5ad_dir = "/vast/scratch/users/shen.m/Census/split_h5ad_based_on_sample_id/2024-07-01/", - pdf_file = "/home/users/allstaff/shen.m/projects/cellNexus/plot6.pdf", - transform = "none") -``` - - -## Row seven -```{r} -sample_summary_df |> dplyr::count(has_negative, max_gt_20, all_integer, has_floating, counts_gap>=0.25) |> dplyr::slice(7) -samples_to_plot <- sample_summary_df |> filter(has_negative, !max_gt_20, !all_integer, !has_floating, counts_gap<0.25) |> - arrange(n_cells) |> - group_by(dataset_id) |> - slice_sample(n = 2) |> - ungroup() -``` - - -```{r eval=FALSE} -plot_raw_counts_hist_pdf(samples_to_plot |> head(10), - h5ad_dir = "/vast/scratch/users/shen.m/Census/split_h5ad_based_on_sample_id/2024-07-01/", - pdf_file = "/home/users/allstaff/shen.m/projects/cellNexus/plot7.pdf", - transform = "none") -``` - -## Row eight -```{r} -sample_summary_df |> dplyr::count(has_negative, max_gt_20, all_integer, has_floating, counts_gap>=0.25) |> dplyr::slice(8) -samples_to_plot <- sample_summary_df |> filter(has_negative, !max_gt_20, !all_integer, !has_floating, counts_gap<0.25) |> - arrange(n_cells) |> - group_by(dataset_id) |> - slice_sample(n = 2) |> - ungroup() -``` - - -```{r eval=FALSE} -plot_raw_counts_hist_pdf(samples_to_plot |> head(10), - h5ad_dir = "/vast/scratch/users/shen.m/Census/split_h5ad_based_on_sample_id/2024-07-01/", - pdf_file = "/home/users/allstaff/shen.m/projects/cellNexus/plot8.pdf", - transform = "none") -``` - -## Row nine -```{r} -sample_summary_df |> dplyr::count(has_negative, max_gt_20, all_integer, has_floating, counts_gap>=0.25) |> dplyr::slice(9) -samples_to_plot <- sample_summary_df |> filter(has_negative, !max_gt_20, !all_integer, !has_floating, counts_gap<0.25) |> - arrange(n_cells) |> - group_by(dataset_id) |> - slice_sample(n = 2) |> - ungroup() -``` - - -```{r eval=FALSE} -plot_raw_counts_hist_pdf(samples_to_plot |> head(10), - h5ad_dir = "/vast/scratch/users/shen.m/Census/split_h5ad_based_on_sample_id/2024-07-01/", - pdf_file = "/home/users/allstaff/shen.m/projects/cellNexus/plot9.pdf", - transform = "none") -``` - -# Final slides and decision - -```{r} -#| label: display-pptx -#| echo: false -#| warning: false -#| message: false - -pptx_file <- "~/projects/cellNexus/census_distribution_decision.pptx" -pptx_file -``` \ No newline at end of file diff --git a/dev/cellnexus-2024-scripts/step6_execute_hpcell_on_census_and_defining_data_tranformation_mengyuan_version.R b/dev/cellnexus-2024-scripts/step6_execute_hpcell_on_census_and_defining_data_tranformation_mengyuan_version.R deleted file mode 100644 index 2a895e6..0000000 --- a/dev/cellnexus-2024-scripts/step6_execute_hpcell_on_census_and_defining_data_tranformation_mengyuan_version.R +++ /dev/null @@ -1,524 +0,0 @@ -# Step 2 -library(dplyr) -library(tibble) -library(glue) -library(purrr) -library(stringr) -library(HPCell) -library(arrow) -library(targets) -library(crew) -library(crew.cluster) -library(duckdb) -directory = "/vast/scratch/users/shen.m/Census/split_h5ad_based_on_sample_id/2024-07-01/" # MODIFY HERE: directory containing per-sample h5ad files -downloaded_samples_tbl <- read_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/census_samples_to_download_groups_MODIFIED.parquet") # MODIFY HERE: input samples metadata parquet -downloaded_samples_tbl <- downloaded_samples_tbl |> - dplyr::rename(cell_number = list_length) |> - mutate(cell_number = cell_number |> as.integer(), - file_name = glue("{directory}{sample_2}.h5ad") |> as.character()) - -# result_directory = "/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024" # MODIFY HERE: directory containing the pre-existing targets store for sample_meta -# -# sample_meta <- tar_read(metadata_dataset_id_common_sample_columns, store = glue("{result_directory}/_targets")) -sample_tbl = downloaded_samples_tbl |> - filter(!dataset_id %in% c("99950e99-2758-41d2-b2c9-643edcdf6d82", "9fcb0b73-c734-40a5-be9c-ace7eea401c9" )) |> - left_join( - cellxgenedp::datasets() |> - select(dataset_id, x_approximate_distribution) |> - distinct(), by = "dataset_id", copy = TRUE) |> - mutate(cell_number = cell_number |> as.integer(), - file_name = glue("{directory}{sample_2}.h5ad") |> as.character()) |> - - left_join( - cellNexus::get_metadata(cache_directory = "/vast/scratch/users/shen.m/cellNexus") |> # MODIFY HERE: cellNexus local cache directory - cellNexus::join_census_table() |> - distinct(sample_id, assay) , - by = c("sample_2" = "sample_id"), - copy = T - ) |> - # Propositional set up expressed genes threshold for panel technologies 500/20K = x/462 - mutate(feature_thresh = ifelse(assay == "BD Rhapsody Targeted mRNA", 11, 200)) - -# Manually updated 300 samples transformation profiles -sample_summary_df = tar_read(sample_summary_df, store = "/vast/scratch/users/shen.m/2024-07-01_census_sample_raw_counts_summary_target_store/_targets") |> # MODIFY HERE: targets store for manually reviewed SCT-failed samples - bind_rows() |> - mutate(max_gt_20 = ifelse(max_val > 20, TRUE, FALSE)) - -impute_x_approximate_distribution <- function(df) { - df |> mutate( - inferred_distribution = case_when( - # 0) When counts gap between 0 and next min value >= 0.25, double log - !has_negative & !max_gt_20 & !all_integer & !has_floating & (counts_gap >= 0.25) ~ "double_log1p" , - - # 1) No negatives, no large values, no integers, no floating - !has_negative & !max_gt_20 & !all_integer & !has_floating & (counts_gap < 0.25) ~ "log1p", - - # 2) No negatives, has large values - !has_negative & max_gt_20 & !all_integer & !has_floating ~ "raw", - - # 3) No negatives, large values, all integer, has floating - !has_negative & max_gt_20 & all_integer & !has_floating ~ "raw", - - # 4) Has negatives, no large values, no integer, no floating. Counts peak at 10 - has_negative & !max_gt_20 & !all_integer & !has_floating ~ "raw_limit_max_to_10", - - # 5) Has negatives and large values - has_negative & max_gt_20 & !all_integer & !has_floating ~ "raw" - ) - ) -} - -sample_summary_df = sample_summary_df |> impute_x_approximate_distribution() |> - mutate(count_upper_bound = case_when( - # 0) When counts gap between 0 and next min value >= 0.25, double log. Max value before exp is 10. - inferred_distribution == "double_log1p" ~ 10, - - # 1) make 10 as max before exp - inferred_distribution == "log1p" ~ 10, - - # 4) Has negatives, no large values, no integer, no floating. Counts peak at 10 - inferred_distribution == "raw_limit_max_to_10" ~ 10, - - # 2,3,5), assign a dummy limit - inferred_distribution == "raw" ~ 9999 - - )) |> - # Inverse distribution - mutate(method_to_apply = case_when(inferred_distribution == "double_log1p" ~ "safe_expm1", - inferred_distribution == "log1p" ~ "expm1", - inferred_distribution == "raw" ~ "identity", - inferred_distribution == "raw_limit_max_to_10" ~ "identity_with_max_limit")) - -sample_tbl = sample_tbl |> left_join(sample_summary_df |> - select(sample_2, - method_to_apply, - dataset_id, - count_upper_bound), - by = c("sample_2", "dataset_id")) - -sample_tbl = sample_tbl |> - - select(file_name, cell_number, dataset_id, sample_2, method_to_apply, assay, count_upper_bound, feature_thresh) - -sample_tbl |> saveRDS("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/updated_transform_sample_tbl_2024_Jul.rds") # MODIFY HERE: output path for sample_tbl RDS - -# # -# # left_join(sample_meta, by = "dataset_id") |> distinct(file_name, tier, cell_number, dataset_id, sample_2, -# # x_normalization, x_approximate_distribution) |> -# mutate(transform_method = case_when(str_like(x_normalization, "C%") ~ "log", -# x_normalization == "none" ~ "log", -# x_normalization == "normalized" ~ "log", -# is.na(x_normalization) & is.na(x_approximate_distribution) ~ "log", -# is.na(x_normalization) & x_approximate_distribution == "NORMAL" ~ "NORMAL", -# is.na(x_normalization) & x_approximate_distribution == "COUNT" ~ "COUNT", -# str_like(x_normalization, "%canpy%") ~ "log1p", -# TRUE ~ x_normalization)) |> -# -# mutate(method_to_apply = case_when(transform_method %in% c("log","LogNormalization","LogNormalize","log-normalization") ~ "exp", -# is.na(x_normalization) & is.na(x_approximate_distribution) ~ "exp", -# str_like(transform_method, "Counts%") ~ "exp", -# str_like(transform_method, "%log2%") ~ "exp", -# transform_method %in% c("log1p", "log1p, base e", "Scanpy", -# "scanpy.api.pp.normalize_per_cell method, scaling factor 10000") ~ "expm1", -# transform_method == "log1p, base 2" ~ "expm1", -# transform_method == "NORMAL" ~ "exp", -# transform_method == "COUNT" ~ "identity", -# is.na(transform_method) ~ "identity" -# ) ) |> -# mutate(comment = case_when(str_like(x_normalization, "Counts%") ~ "a checkpoint for max value of Assay must <= 50", -# is.na(x_normalization) & is.na(x_approximate_distribution) ~ "round negative value to 0", -# x_normalization == "normalized" ~ "round negative value to 0" -# )) -# -# -# # Append assay column -# sample_tbl = sample_tbl |> left_join(cellNexus::get_metadata(cache_directory = "/vast/scratch/users/shen.m/cellNexus") |> # MODIFY HERE: cellNexus local cache directory -# distinct(sample_id, assay) , -# by = c("sample_2" = "sample_id"), -# copy = T) -# -# -# sample_tbl = sample_tbl |> mutate(count_upper_bound = 20, -# # base our filtering on % of expressed genes for panel technologies 500/20K = x/462 -# feature_thresh = ifelse(assay == "BD Rhapsody Targeted mRNA", 11, 200)) -# -# sample_tbl <- saveRDS("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/sample_tbl_2024_Jul.rds") # MODIFY HERE: output path for sample_tbl RDS -# -# sliced_sample_tbl = -# sample_tbl |> -# filter(!dataset_id %in% c("99950e99-2758-41d2-b2c9-643edcdf6d82", "9fcb0b73-c734-40a5-be9c-ace7eea401c9" )) |> -# dplyr::select(file_name, tier, cell_number, dataset_id, sample_2, method_to_apply, assay, count_upper_bound, feature_thresh) - - - - -# #sliced_sample_tbl |> write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/sliced_sample_tbl_2024_Jul.parquet") -# sliced_sample_tbl <- read_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/sliced_sample_tbl_2024_Jul.parquet") # MODIFY HERE: output path for sliced_sample_tbl RDS - -# Enable sample_names.rds to store sample names for the input -sample_names <- - sample_tbl |> - pull(file_name) |> - set_names(sample_tbl |> pull(sample_2)) -functions = sample_tbl |> pull(method_to_apply) -feature_thresh = sample_tbl |> pull(feature_thresh) -count_upper_bound = sample_tbl |> pull(count_upper_bound) - - -my_store = "/vast/scratch/users/shen.m/cellNexus/2024-07-01/process_updated_samples_transform_hpcell_target_store_v1" # MODIFY HERE: HPCell targets store (used throughout this script) -job::job({ - - library(HPCell) - - sample_names |> - initialise_hpc( - store = my_store, - gene_nomenclature = "ensembl", - data_container_type = "anndata", - computing_resources = list( - - crew.cluster::crew_controller_slurm( - name = "elastic", - workers = 300, - tasks_max = 20, - seconds_idle = 30, - crashes_error = 10, - options_cluster = crew.cluster::crew_options_slurm( - #memory_gigabytes_required = c(20, 35, 50, 75, 100, 150), - #memory_gigabytes_required = c(90, 120, 150, 180, 200), - #memory_gigabytes_required = c(70, 80, 100, 150, 200), - memory_gigabytes_required = c(45, 60, 75, 100, 120, 150), - cpus_per_task = c(2, 2, 5, 10, 20), - time_minutes = c(60*4, 60*4, 60*4, 60*4, 60*4), - verbose = T - ) - ) - - ), - verbosity = "summary", - update = "never", - #update = "thorough", - error = "continue", - garbage_collection = 100, - workspace_on_error = TRUE - - ) |> - transform_assay(fx = functions, target_output = "sce_transformed", scale_max = count_upper_bound) |> - # - # # # Remove empty outliers based on RNA count threshold per cell - # remove_empty_threshold(target_input = "sce_transformed", RNA_feature_threshold = feature_thresh) |> - # - # # Annotation - # annotate_cell_type(target_input = "sce_transformed", azimuth_reference = "pbmcref") |> - # - # # Cell type harmonisation - # celltype_consensus_constructor(target_input = "sce_transformed", - # target_output = "cell_type_concensus_tbl") |> - # - # # Alive identification - # remove_dead_scuttle(target_input = "sce_transformed", target_annotation = "cell_type_concensus_tbl", - # group_by = "cell_type_unified_ensemble") |> - # - # # Doublets identification - # remove_doublets_scDblFinder(target_input = "sce_transformed") |> - - # # SCT - # normalise_abundance_seurat_SCT(target_input = "sce_transformed", factors_to_regress = c( - # "subsets_Mito_percent", - # "subsets_Ribo_percent")) |> - - # # Pseudobulk - # calculate_pseudobulk(target_input = "sce_transformed", - # group_by = "cell_type_unified_ensemble") |> - - # # metacell - # cluster_metacell(target_input = "sce_transformed", group_by = "cell_type_unified_ensemble") |> - - # # Cell Chat - # ligand_receptor_cellchat(target_input = "sce_transformed", - # group_by = "cell_type_unified_ensemble") |> - - print() - - -}) - -# # View target metadata if needed -# tar_meta(store = my_store) |> filter(!is.na(error)) |> arrange(desc(time)) |> View() -# tar_meta(store = my_store) |> filter(!is.na(error)) |> distinct(name, error) -# tar_meta(starts_with("annotation_tbl_"), store = "/vast/scratch/users/shen.m/cellNexus_target_store") |> -# filter(!data |> is.na()) |> arrange(desc(time)) |> select(error, name) -# -# # Debug cellchat -# tar_workspace(ligand_receptor_tbl_f1dcd76261dd9c86, store = my_store, -# script = paste0(my_store,".R")) -# debugonce(cell_communication) -# cell_communication(sce_transformed, empty_droplets_tbl = empty_tbl, alive_identification_tbl = alive_tbl, -# doublet_identification_tbl = doublet_tbl, cell_type_tbl = cell_type_concensus_tbl, -# cell_type_column = "cell_type_unified_ensemble", -# feature_nomenclature = gene_nomenclature) - - -#' Pipeline for Lightening Annotations in High-Performance Computing Environment -#' -#' This pipeline is designed to read, process, and "lighten" large annotation tables in an HPC environment. -#' It uses the `targets` package for reproducibility and `crew` for efficient job scheduling on a Slurm cluster. -#' The `lighten_annotation` function selects and processes specific columns from large tables to reduce memory usage. -#' -#' The pipeline consists of: -#' - **Crew Controllers**: Four tiers of Slurm controllers with varying memory allocations to optimize resource usage. -#' - **Targets**: -#' - `my_store`: Defines the path to the target storage directory, ensuring all targets use the correct storage location. -#' - `target_name`: Retrieves metadata to identify branch targets for annotation. -#' - `annotation_tbl_light`: Applies `lighten_annotation` to process each target name, optimally running with `tier_1` resources. -#' -#' @libraries: -#' - `dplyr`, `magrittr`, `tibble`, `targets`, `tarchetypes` for data manipulation and pipeline structure. -#' - `crew`, `crew.cluster` for parallel computation and cluster scheduling in an HPC environment. -#' -#' @options: -#' - Memory settings, garbage collection frequency, and error handling are set to handle large data efficiently. -#' - The `cue` option is set to `never` for forced target updates if needed. -#' - `controller` is a group of Slurm controllers to manage computation across memory tiers. -#' -#' @function `lighten_annotation`: Processes each annotation table target, unnesting and selecting specific columns to reduce data size. -#' -#' @example Usage: -#' The pipeline script is saved as `/vast/scratch/users/shen.m/lighten_annotation_tbl_target.R` by tar_script and can be run using `tar_make()`. -tar_script({ - library(dplyr) - library(magrittr) - library(tibble) - library(targets) - library(tarchetypes) - library(crew) - library(crew.cluster) - tar_option_set( - memory = "transient", - garbage_collection = 100, - storage = "worker", - retrieval = "worker", - error = "continue", - cue = tar_cue(mode = "never"), - controller = crew_controller_group( - list( - crew_controller_slurm( - name = "tier_1", - script_lines = "#SBATCH --mem 8G", - slurm_cpus_per_task = 1, - workers = 200, - tasks_max = 10, - verbose = T, - seconds_idle = 30, - slurm_time_minutes = 480 - ), - - crew_controller_slurm( - name = "tier_2", - script_lines = "#SBATCH --mem 10G", - slurm_cpus_per_task = 1, - workers = 200, - tasks_max = 10, - verbose = T, - seconds_idle = 30, - slurm_time_minutes = 480 - ), - crew_controller_slurm( - name = "tier_3", - script_lines = "#SBATCH --mem 15G", - slurm_cpus_per_task = 1, - workers = 200, - tasks_max = 10, - verbose = T, - seconds_idle = 30, - slurm_time_minutes = 480 - ), - crew_controller_slurm( - name = "tier_4", - script_lines = "#SBATCH --mem 50G", - slurm_cpus_per_task = 1, - workers = 30, - tasks_max = 10, - verbose = T, - seconds_idle = 30, - slurm_time_minutes = 480 - ) - ) - ), - trust_object_timestamps = TRUE - ) - - lighten_annotation = function(target_name, my_store ){ - annotation_tbl = tar_read_raw( target_name, store = my_store ) - if(annotation_tbl |> is.null()) { - warning("this annotation is null -> ", target_name) - return(NULL) - } - - annotation_tbl |> - unnest(blueprint_scores_fine) |> - select(.cell, blueprint_first.labels.fine, monaco_first.labels.fine, any_of("azimuth_predicted.celltype.l2"), monaco_scores_fine, contains("macro"), contains("CD4") ) |> - unnest(monaco_scores_fine) |> - select(.cell, blueprint_first.labels.fine, monaco_first.labels.fine, any_of("azimuth_predicted.celltype.l2"), contains("macro") , contains("CD4"), contains("helper"), contains("Th")) |> - rename(cell_ = .cell) - } - - list( - - # The input DO NOT DELETE - tar_target(my_store, "/vast/scratch/users/shen.m/cellNexus/2024-07-01/process_samples_hpcell_target_store", deployment = "main"), # MODIFY HERE: HPCell targets store (must match my_store above) - - tar_target( - target_name, - tar_meta( - starts_with("annotation_tbl_"), - store = my_store) |> - filter(type=="branch") |> - pull(name), - deployment = "main" - ) , - - tar_target( - annotation_tbl_light, - lighten_annotation(target_name, my_store), - packages = c("dplyr", "tidyr"), - pattern = map(target_name), - resources = tar_resources( - crew = tar_resources_crew(controller = "tier_1") - ) - ) - ) - - -}, script = "/vast/scratch/users/shen.m/lighten_annotation_tbl_target_2024_Jul.R", ask = FALSE) # MODIFY HERE: output path for the tar_script file - -job::job({ - - tar_make( - script = "/vast/scratch/users/shen.m/lighten_annotation_tbl_target_2024_Jul.R", # MODIFY HERE: must match the script path above - store = "/vast/scratch/users/shen.m/lighten_annotation_tbl_target_2024_Jul", # MODIFY HERE: targets store for the lighten-annotation pipeline - reporter = "summary" - ) - -}) - -# Sample metadata -library(arrow) -library(dplyr) -library(duckdb) -library(targets) - -# Write annotation light -# MODIFY HERE: base cell_metadata parquet path inside the SQL string below -cell_metadata <- - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql("SELECT * FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgenedp_Apr_2024/cell_metadata.parquet')") - ) |> - mutate(cell_ = paste0(cell_, "___", dataset_id)) |> - select(cell_, observation_joinid, contains("cell_type"), dataset_id, self_reported_ethnicity, tissue, donor_id, sample_id, is_primary_data, assay) - - -cell_annotation = - tar_read(annotation_tbl_light, store = "/vast/scratch/users/shen.m/lighten_annotation_tbl_target_2024_Jul") |> # MODIFY HERE: lighten-annotation targets store (must match the tar_make store above) - dplyr::rename( - blueprint_first_labels_fine = blueprint_first.labels.fine, - monaco_first_labels_fine = monaco_first.labels.fine, - azimuth_predicted_celltype_l2 = azimuth_predicted.celltype.l2 - ) - -cell_annotation = cell_annotation |> mutate( - blueprint_first_labels_fine = ifelse(is.na(blueprint_first_labels_fine), "Other", blueprint_first_labels_fine), - monaco_first_labels_fine = ifelse(is.na(monaco_first_labels_fine), "Other", monaco_first_labels_fine), - azimuth_predicted_celltype_l2=ifelse(is.na(azimuth_predicted_celltype_l2), "Other", azimuth_predicted_celltype_l2)) - -empty_droplet = - tar_read(empty_tbl, store = "/vast/scratch/users/shen.m/cellNexus/2024-07-01/process_samples_hpcell_target_store") |> # MODIFY HERE: HPCell targets store (must match my_store above) - bind_rows() |> - dplyr::rename(cell_ = .cell) - -alive_cells = - tar_read(alive_tbl, store = "/vast/scratch/users/shen.m/cellNexus/2024-07-01/process_samples_hpcell_target_store") |> # MODIFY HERE: HPCell targets store (must match my_store above) - bind_rows() |> - dplyr::rename(cell_ = .cell) - -doublet_cells = - tar_read(doublet_tbl, store ="/vast/scratch/users/shen.m/cellNexus/2024-07-01/process_samples_hpcell_target_store") |> # MODIFY HERE: HPCell targets store (must match my_store above) - bind_rows() |> - dplyr::rename(cell_ = .cell) - -metacell = - tar_read(metacell_tbl, store = "/vast/scratch/users/shen.m/cellNexus/2024-07-01/process_samples_hpcell_target_store") |> # MODIFY HERE: HPCell targets store (must match my_store above) - bind_rows() |> - dplyr::rename(cell_ = cell) |> - dplyr::rename_with( - ~ stringr::str_replace(.x, "^gamma", "metacell_"), - starts_with("gamma") - ) - -# Save cell type concensus tbl from HPCell output to disk -cell_type_concensus_tbl = tar_read(cell_type_concensus_tbl, store = "/vast/scratch/users/shen.m/cellNexus/2024-07-01/process_samples_hpcell_target_store") |> # MODIFY HERE: HPCell targets store (must match my_store above) - bind_rows() |> - dplyr::rename(cell_ = .cell) - -cell_type_concensus_tbl = cell_type_concensus_tbl |> mutate(cell_type_unified_ensemble = - ifelse(is.na(cell_type_unified_ensemble), - "Unknown", - cell_type_unified_ensemble)) - -# This command needs a big memory machine -cell_metadata_joined = cell_metadata |> - left_join(empty_droplet, copy=TRUE) |> - left_join(cell_type_concensus_tbl, copy=TRUE) |> - left_join(alive_cells, copy=TRUE) |> - left_join(doublet_cells, copy=TRUE) |> - left_join(metacell, copy=TRUE) - -cell_metadata_joined |> filter(is.na(blueprint_first_labels_fine)) - -cell_metadata_joined2 = cell_metadata_joined |> as_tibble() |> - # Match to how pseudobulk annotations get parsed in HPCell/R/functions preprocessing_output() - mutate(cell_type_unified_ensemble = ifelse(cell_type_unified_ensemble |> is.na(), "Unknown", cell_type_unified_ensemble)) |> - mutate(data_driven_ensemble = ifelse(data_driven_ensemble |> is.na(), "Unknown", data_driven_ensemble)) |> - mutate(blueprint_first_labels_fine = ifelse(blueprint_first_labels_fine |> is.na(), "Other", blueprint_first_labels_fine)) |> - mutate(monaco_first_labels_fine = ifelse(monaco_first_labels_fine |> is.na(), "Other", monaco_first_labels_fine)) |> - mutate(azimuth_predicted_celltype_l2 = ifelse(azimuth_predicted_celltype_l2 |> is.na(), "Other", azimuth_predicted_celltype_l2)) |> - mutate(azimuth = ifelse(azimuth |> is.na(), "Other", azimuth)) |> - mutate(blueprint = ifelse(blueprint |> is.na(), "Other", blueprint)) |> - mutate(monaco = ifelse(monaco |> is.na(), "Other", monaco)) - -cell_metadata_joined2 |> - arrow::write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_annotation_2024_Jul.parquet", # MODIFY HERE: output cell annotation parquet (used as input to step6) - compression = "zstd") - -# Cellchat output -ligand_receptor_tbl = tar_read(ligand_receptor_tbl, store = "/vast/scratch/users/shen.m/cellNexus/2024-07-01/process_samples_hpcell_target_store") |> bind_rows() # MODIFY HERE: HPCell targets store (must match my_store above) -# save -con <- dbConnect(duckdb::duckdb(), dbdir = "~/cellxgene_curated/metadata_cellxgene_mengyuan/cellNexus_lr_signaling_pathway_strength.duckdb") # MODIFY HERE: output DuckDB file for ligand-receptor results -duckdb::dbWriteTable(con, "lr_pathway_table", ligand_receptor_tbl, overwrite = TRUE) -dbDisconnect(con) - - -# Helper function to save parquet read by duckdb to parquet on disk -# write_parquet_to_parquet = function(data_tbl, output_parquet, compression = "gzip") { -# -# # Establish connection to DuckDB in-memory database -# con_write <- dbConnect(duckdb::duckdb(), dbdir = ":memory:") -# -# # Register `data_tbl` within the DuckDB connection (this doesn't load it into memory) -# duckdb::duckdb_register(con_write, "data_tbl_view", data_tbl) -# -# # Use DuckDB's COPY command to write `data_tbl` directly to Parquet with compression -# copy_query <- paste0(" -# COPY data_tbl_view TO '", output_parquet, "' (FORMAT PARQUET, COMPRESSION '", compression, "'); -# ") -# -# # Execute the COPY command -# dbExecute(con_write, copy_query) -# -# # Unregister the temporary view -# duckdb::duckdb_unregister(con_write, "data_tbl_view") -# -# # Disconnect from the database -# dbDisconnect(con_write, shutdown = TRUE) -# } diff --git a/dev/cellnexus-2024-scripts/step6_prepare_local_cache_splitting_du_dataset_and_cell_type_mengyuan_version.R b/dev/cellnexus-2024-scripts/step6_prepare_local_cache_splitting_du_dataset_and_cell_type_mengyuan_version.R deleted file mode 100644 index 849490c..0000000 --- a/dev/cellnexus-2024-scripts/step6_prepare_local_cache_splitting_du_dataset_and_cell_type_mengyuan_version.R +++ /dev/null @@ -1,977 +0,0 @@ -# Step6 -# Group samples by dataset_id, cell_type - -# This script sets up a robust and scalable data processing pipeline for single-cell RNA sequencing (scRNA-seq) datasets using the targets package in R, which facilitates reproducible and efficient workflows. Specifically, the code orchestrates the ingestion and preprocessing of multiple SingleCellExperiment objects corresponding to different datasets (dataset_id) and targets (target_name). It leverages high-performance computing resources through the crew package, configuring multiple SLURM-based controllers (tier_1 to tier_4) to handle varying computational loads efficiently. -# -# The pipeline performs several key steps: -# -# 1. Data Retrieval: It reads raw SingleCellExperiment objects for each target, ensuring that only successfully loaded data proceeds further. -# 2. Normalization: Calculates Counts Per Million (CPM) for each cell to normalize gene expression levels across cells and samples. -# 3. Data Aggregation: Groups the data by dataset_id and tar_group, then combines the SingleCellExperiment objects within each group into a single object, effectively consolidating the data for each dataset. -# 4. Metadata Integration: Joins additional metadata, such as cell types, by connecting to a DuckDB database and fetching relevant information from a Parquet file. This enriches the single-cell data with essential annotations. -# 5. Cell Type Segmentation: Splits the combined SingleCellExperiment objects into separate objects based on cell_type, facilitating downstream analyses that are specific to each cell type. -# 6. Data Saving with Error Handling: Generates unique identifiers for each cell type within a dataset and saves both the raw counts and CPM-normalized data to specified directories. It includes special handling for cases where a cell type has only one cell, duplicating the data to prevent errors during the saving process. -# -# By integrating targets, crew, and various data manipulation packages (dplyr, tidyverse, SingleCellExperiment), this script ensures that large-scale scRNA-seq data processing is efficient, reproducible, and capable of leveraging parallel computing resources. It is designed to handle edge cases gracefully and provides a clear framework for preprocessing scRNA-seq data, which is essential for subsequent analyses such as clustering, differential expression, and cell type identification. - - -library(arrow) -library(dplyr) -library(duckdb) - -job::job({ - - get_file_ids = function(cell_annotation ){ - sample_chunk_df = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue::glue("SELECT * FROM read_parquet('{cell_annotation}')")) - ) |> - # Define chunks - dplyr::count(dataset_id, sample_id, name = "cell_count") |> # Ensure unique dataset_id and sample_id combinations - distinct(dataset_id, sample_id, cell_count) |> # Ensure unique dataset_id and sample_id combinations - group_by(dataset_id) |> - dbplyr::window_order(dataset_id, cell_count, sample_id) |> # Ensure order. Note: order cell_count only is not enough because it needs a secondary tie-breaker - mutate(sample_index = row_number()) |> # Create sequential index within each dataset - mutate(sample_chunk = (sample_index - 1) %/% 1000 + 1) |> # Assign chunks (up to 1000 samples per chunk) - mutate(sample_pseudobulk_chunk = (sample_index - 1) %/% 250 + 1) |> # Max combination of dataset_id, sample_pseudobulk_chunk and file_id_pseudobulk up to 10000 - mutate(cell_chunk = cumsum(cell_count) %/% 100000 + 1) |> # max 20K cells per sample - ungroup() - - # Test whether cell_chunk and sample_chunk are unique for this sample - run_chunk_once <- function(column_name, id) { - sample_chunk_df |> filter(sample_id == id) |> pull(!!column_name) - } - - sample_chunk_results <- replicate(20, run_chunk_once("sample_chunk", "d6e942a09a140ee8bb6f0c3da8defea4___exp7-human-150well."), simplify = FALSE) - sample_chunk_identical <- all(sapply(sample_chunk_results[-1], function(x) identical(x, sample_chunk_results[[1]]))) - if (!sample_chunk_identical) { - stop("Inconsistent sample chunk value was generated in multiple runs, this will lead to file id changes") - } - - cell_chunk_results <- replicate(20, run_chunk_once("cell_chunk", "d6e942a09a140ee8bb6f0c3da8defea4___exp7-human-150well."), simplify = FALSE) - cell_chunk_identical <- all(sapply(cell_chunk_results[-1], function(x) identical(x, cell_chunk_results[[1]]))) - if (!cell_chunk_identical) { - stop("Inconsistent cell chunk value was generated in multiple runs, this will lead to file id changes") - } - - - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue::glue("SELECT * FROM read_parquet('{cell_annotation}')")) - ) |> - # Cells in cell_annotation could be more than cells in cell_consensus. In order to avoid NA happens in cell_consensus cell_type column - mutate(cell_type_unified_ensemble = ifelse(cell_type_unified_ensemble |> is.na(), - "Unknown", - cell_type_unified_ensemble)) |> - - left_join(sample_chunk_df |> select(dataset_id, sample_chunk, sample_pseudobulk_chunk, cell_chunk, sample_id), copy=TRUE) |> - - # Define chunks - group_by(dataset_id, sample_chunk, cell_chunk, sample_pseudobulk_chunk, cell_type, sample_id) |> - summarise(cell_count = n(), .groups = "drop") |> - group_by(dataset_id, sample_chunk, cell_chunk, cell_type) |> - dbplyr::window_order(desc(cell_count), sample_id) |> # Important! - mutate(chunk = cumsum(cell_count) %/% 20000 + 1) |> # max 20K cells per sample - ungroup() |> - as_tibble() |> - - # Single cell file ID - mutate(file_id_cellNexus_single_cell = - glue::glue("{dataset_id}___{sample_chunk}___{cell_chunk}___{cell_type}") |> - sapply(digest::digest) |> - paste0("___", chunk, ".h5ad") - ) |> - - # Pseudobulk file id - mutate(file_id_cellNexus_pseudobulk = - glue::glue("{dataset_id}___{sample_pseudobulk_chunk}") |> - sapply(digest::digest) |> - paste0("___", chunk, ".h5ad")) - - } - - get_file_ids( - "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_annotation_2024_Jul.parquet" # MODIFY HERE: input cell annotation parquet - ) |> - write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/file_id_cellNexus_single_cell_2024_Jul.parquet") # MODIFY HERE: output file_id parquet - - gc() - - con <- dbConnect(duckdb::duckdb(), dbdir = ":memory:") - - dir.create("/vast/scratch/users/shen.m/duckdb_tmp", showWarnings = FALSE) # MODIFY HERE: duckdb temp directory - - DBI::dbExecute( - con, - "SET temp_directory='/vast/scratch/users/shen.m/duckdb_tmp';" # MODIFY HERE: duckdb temp directory (must match dir.create above) - ) - - # Create a view for cell_annotation in DuckDB - # MODIFY HERE: cell_metadata parquet path inside the SQL string below - dbExecute(con, " - CREATE VIEW cell_metadata AS - SELECT - CONCAT(cell_, '___', dataset_id) AS cell_, - * EXCLUDE (cell_, dataset_id_1, X_umap1, X_umap2, sample_placeholder, cell_type) -- drop original cell_ and dataset_id_1 - FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata.parquet') -") - - # MODIFY HERE: cell_annotation parquet path inside the SQL string below - dbExecute(con, " - CREATE VIEW hpcell_output_metadata AS - SELECT * EXCLUDE ( - observation_joinid, - cell_type_ontology_term_id, - assay, - donor_id, - is_primary_data, - self_reported_ethnicity, - tissue, - azimuth, - blueprint, - monaco, - subsets_Mito_sum, - subsets_Mito_detected, - ensemble_joinid, - cell_type_unified, - data_driven_ensemble, - observation_originalid -) - - FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_annotation_2024_Jul.parquet') -") - -# MODIFY HERE: file_id parquet path inside the SQL string below (should match the write_parquet output above) -dbExecute(con, " - CREATE VIEW file_id_cellNexus_single_cell AS - SELECT - dataset_id, - sample_chunk, - cell_chunk, - sample_pseudobulk_chunk, - cell_type, - sample_id, - file_id_cellNexus_single_cell, - file_id_cellNexus_pseudobulk - FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/file_id_cellNexus_single_cell_2024_Jul.parquet') -") - -# MODIFY HERE: transformation data frame -dbExecute(con, " - CREATE VIEW sample_distribution_method_tbl AS - SELECT - sample_2 AS sample_id, - count_upper_bound, - feature_thresh AS nfeature_expressed_thresh, - method_to_apply AS inverse_transform - FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/sliced_sample_tbl_2024_Jul.parquet') -") - -# Perform the left join and save to Parquet -copy_query <- " - COPY ( - SELECT - cell_metadata.cell_ AS cell_id, -- Rename cell_ to cell_id - COALESCE(hpcell_output_metadata.alive, FALSE) AS alive, -- Set alive column NULL to FALSE - cell_metadata.* EXCLUDE (cell_), -- drop cell_ since it's already aliased as cell_id - hpcell_output_metadata.* EXCLUDE (cell_, dataset_id, sample_id, alive), -- Deduplicate join keys, and aliased column - file_id_cellNexus_single_cell.* EXCLUDE (sample_id, dataset_id, cell_type), -- Deduplicate join keys - sample_distribution_method_tbl.* EXCLUDE (sample_id) -- Deduplicate join keys - FROM cell_metadata - - LEFT JOIN hpcell_output_metadata - ON hpcell_output_metadata.cell_ = cell_metadata.cell_ - AND hpcell_output_metadata.dataset_id = cell_metadata.dataset_id - - LEFT JOIN file_id_cellNexus_single_cell - ON file_id_cellNexus_single_cell.sample_id = hpcell_output_metadata.sample_id - AND file_id_cellNexus_single_cell.dataset_id = hpcell_output_metadata.dataset_id - AND file_id_cellNexus_single_cell.cell_type = hpcell_output_metadata.cell_type - - LEFT JOIN sample_distribution_method_tbl - ON sample_distribution_method_tbl.sample_id = cell_metadata.sample_id - - WHERE cell_metadata.dataset_id NOT IN ('99950e99-2758-41d2-b2c9-643edcdf6d82', '9fcb0b73-c734-40a5-be9c-ace7eea401c9') -- (THESE TWO DATASETS DOESNT contain meaningful data - no observation_joinid etc), thus was excluded in the final metadata. - - ) TO '/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_0_mengyuan.parquet' -- MODIFY HERE: output merged metadata parquet (v1_2_2) - (FORMAT PARQUET, COMPRESSION 'gzip'); -" - -# Execute the final query to write the result to a Parquet file -dbExecute(con, copy_query) - -# Disconnect from the database -dbDisconnect(con, shutdown = TRUE) - -print("Done.") -}) - -# We decided to make cell_id lighter without re-run everything in HPCell pipeline. Here to swap cell_id in the metadata -# cell_map is processed in a separate target script step6_supp_dataset_cell_map.R -job::job({ - con <- dbConnect(duckdb::duckdb(), dbdir = ":memory:") - - # Create a view for cell_annotation in DuckDB - # MODIFY HERE: v1_2_2 merged metadata parquet path inside the SQL string below (should match the COPY TO output above) - dbExecute(con, " - CREATE VIEW cell_metadata AS - SELECT * - FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_0_mengyuan.parquet') -") - - # MODIFY HERE: cell_id dictionary parquet path inside the SQL string below - dbExecute(con, " - CREATE VIEW cell_map AS - SELECT * - FROM read_parquet('/vast/projects//cellxgene_curated/metadata_cellxgene_mengyuan/file_id_cell_id_dict_v1_1_1_Jul_2024.parquet') -") - - # Perform the left join and save to Parquet - copy_query <- " - COPY ( - SELECT cell_metadata.*, - cell_map.* EXCLUDE (cell_id, file_id_cellNexus_single_cell) - FROM cell_metadata - - LEFT JOIN cell_map - ON cell_metadata.cell_id = cell_map.cell_id - AND cell_metadata.file_id_cellNexus_single_cell = cell_map.file_id_cellNexus_single_cell - - ) TO '/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_1_mengyuan.parquet' -- MODIFY HERE: output final metadata parquet with new cell IDs (v1_3_2) - (FORMAT PARQUET, COMPRESSION 'gzip'); -" - - # Execute the final query to write the result to a Parquet file - dbExecute(con, copy_query) - - # Disconnect from the database - dbDisconnect(con, shutdown = TRUE) - - print("Done.") - - -}) - - - -# MODIFY HERE: final metadata parquet path used for the targets pipeline (should match the COPY TO output above) -cell_metadata = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql("SELECT * FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_1_mengyuan.parquet')") - ) - -library(targets) -library(tidyverse) -store_file_cellNexus = "/vast/scratch/users/shen.m/targets_prepare_database_split_datasets_chunked_1_4_1_single_cell" # MODIFY HERE: targets store directory for this pipeline - -tar_script({ - library(dplyr) - library(magrittr) - library(tibble) - library(targets) - library(tarchetypes) - library(crew) - library(crew.cluster) - - # Helper (optional) to avoid repetition - new_elastic <- function(name, mem_gb, time_min, workers, crashes_max, cpus_per_task = 2, backup = NULL) { - crew_controller_slurm( - name = name, - workers = workers, - crashes_max = crashes_max, - seconds_idle = 30, - options_cluster = crew_options_slurm( - memory_gigabytes_required = mem_gb, - cpus_per_task = cpus_per_task, - time_minutes = time_min - ), - backup = backup - ) - } - - # Small → large, with fallbacks to the next size up - elastic_160 <- new_elastic("elastic_160", 160, 60 * 24, workers = 8, crashes_max = 2) - elastic_120 <- new_elastic("elastic_120", 120, 60 * 4, workers = 16, crashes_max = 1, cpus_per_task = 8, backup = elastic_160) - elastic_80 <- new_elastic("elastic_80", 80, 60 * 4, workers = 24, crashes_max = 1, cpus_per_task = 8, backup = elastic_120) - elastic_40 <- new_elastic("elastic_40", 40, 60 * 4, workers = 32, crashes_max = 1, cpus_per_task = 8, backup = elastic_80) - elastic_20 <- new_elastic("elastic_20", 20, 60 * 4, workers = 48, crashes_max = 1, cpus_per_task = 8, backup = elastic_40) - elastic_10 <- new_elastic("elastic_10", 10, 60 * 4, workers = 150, crashes_max = 6, cpus_per_task = 8, backup = elastic_20) - - elastic_5_minimal <- new_elastic("elastic_5_minimal", 5, 60 * 4, workers = 300, crashes_max = 6, cpus_per_task = 8, backup = elastic_10) - - - # Group for targets (small → large) - controllers <- crew_controller_group( - elastic_10, elastic_20, elastic_40, elastic_80, elastic_120, elastic_160, elastic_5_minimal - ) - tar_option_set( - memory = "transient", - garbage_collection = 100, - storage = "worker", - retrieval = "worker", - error = "continue", - cue = tar_cue(mode = "never"), - format = "qs", - #debug = "dataset_id_sct_ea377f6e2d0ae2b7", - workspace_on_error = TRUE, - controller = controllers, - trust_object_timestamps = TRUE, - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_5_minimal") - ) - ) - - save_anndata = function(dataset_id_sce, cache_directory){ - - dir.create(cache_directory, showWarnings = FALSE, recursive = TRUE) - - .x = dataset_id_sce |> pull(sce) |> _[[1]] - .y = dataset_id_sce |> pull(file_id_cellNexus_single_cell) |> _[[1]] |> str_remove("\\.h5ad") - - .x |> assays() |> names() = "counts" - - # Save the experiment data to the specified counts cache directory - .x |> save_experiment_data(glue("{cache_directory}/{.y}")) - - return(TRUE) # Indicate successful saving - - - } - - # Because they have an inconsistent failure. If I start the pipeline again they might work. Strange. - insistent_save_anndata <- purrr::insistently(save_anndata, rate = purrr::rate_delay(pause = 60, max_times = 3), quiet = FALSE) - - save_anndata_cpm = function(dataset_id_sce, cache_directory){ - - dir.create(cache_directory, showWarnings = FALSE, recursive = TRUE) - - # # Parallelise - dataset_id_sce |> - purrr::transpose() |> - lapply( - FUN = function(x) { - - .x = x[[2]] - .y = x[[1]] |> str_remove("\\.h5ad") - - # Check if the 'sce' has only one cell (column) - if(ncol(assay(.x)) == 1) { - - # Duplicate the assay to prevent saving errors due to single-column matrices - my_assay = cbind(assay(.x), assay(.x)) - # Rename the second column to distinguish it - colnames(my_assay)[2] = paste0("DUMMY", "___", colnames(my_assay)[2]) - - cd = colData(.x) - cd = cd |> rbind(cd) - rownames(cd)[2] = paste0("DUMMY", "___", rownames(cd)[2]) - - - - .x = SingleCellExperiment(assay = list( my_assay ) |> set_names(names(assays(.x))[1]), colData = cd) - } - - - # # TEMPORARY FOR SOME REASON THE MIN COUNTS IS NOT 0 FOR SOME SAMPLES - # .x = HPCell:::check_if_assay_minimum_count_is_zero_and_correct_TEMPORARY(.x, assays(.x) |> names() |> _[1], subset_up_to_number_of_cells = 100) - - # CALCULATE CPM - .x = SingleCellExperiment(assay = list( cpm = calculateCPM(.x, assay.type = names(assays(.x))[1])), colData = colData(.x)) - - # Save the experiment data to the specified counts cache directory - .x |> save_experiment_data(glue("{cache_directory}/{.y}")) - - return(TRUE) # Indicate successful saving - } - - ) - - return("saved") - - } - - # Because they have an inconsistent failure. If I start the pipeline again they might work. Strange. - insistent_save_anndata_cpm <- purrr::insistently(save_anndata_cpm, rate = purrr::rate_delay(pause = 60, max_times = 3), quiet = FALSE) - - - # Function to process matrix in vertical slices - process_matrix_in_slices <- function(h5_matrix, output_filepath, output_filepath_temp, chunk_size = 1000) { - # Load the HDF5 matrix - n_rows <- dim(h5_matrix)[1] - n_cols <- dim(h5_matrix)[2] - - if (file.exists(output_filepath)) { - file.remove(output_filepath) - cat("Existing output file removed.\n") - } - if (file.exists(output_filepath_temp)) { - file.remove(output_filepath_temp) - cat("Existing output file removed.\n") - } - - # Create an empty list to hold the slices - slice_list <- list() - - # Loop through the matrix in chunks - for (start_col in seq(1, n_cols, by = chunk_size)) { - end_col <- min(start_col + chunk_size - 1, n_cols) - cat("Processing columns", start_col, "to", end_col, "\n") - - # Extract a slice of the matrix - matrix_slice <- as.matrix(h5_matrix[, start_col:end_col, drop=FALSE]) - - # Calculate ranks for the slice - ranked_slice <- singscore::rankGenes(matrix_slice) %>% `-` (1) - - # Convert the ranked slice to sparse format - sparse_ranked_slice <- as(ranked_slice, "CsparseMatrix") - - # Write the slice to the output HDF5 file - HDF5Array::writeHDF5Array( - sparse_ranked_slice, - filepath = output_filepath_temp, - name = paste0("rank_", start_col, "_to_", end_col), - as.sparse = TRUE, - H5type = "H5T_STD_I32LE" - ) - - # Store the slice name for later binding - slice_list[[length(slice_list) + 1]] <- paste0("rank_", start_col, "_to_", end_col) - } - - - slice_list |> map(~HDF5Array::HDF5Array(output_filepath_temp, name =.x)) |> do.call(cbind, args=_) - - } - - save_rank_per_cell = function(dataset_id_sce, cache_directory){ - - dir.create(cache_directory, recursive = TRUE, showWarnings = FALSE) - - .x = dataset_id_sce |> pull(sce) |> _[[1]] - .y = dataset_id_sce |> pull(file_id_cellNexus_single_cell) |> _[[1]] |> str_remove("\\.h5ad") - - # Check if the 'sce' has only one cell (column) - if(ncol(assay(.x)) == 1) { - - # Duplicate the assay to prevent saving errors due to single-column matrices - my_assay = cbind(assay(.x), assay(.x)) - # Rename the second column to distinguish it - colnames(my_assay)[2] = paste0("DUMMY", "___", colnames(my_assay)[2]) - - cd = colData(.x) - cd = cd |> rbind(cd) - rownames(cd)[2] = paste0("DUMMY", "___", rownames(cd)[2]) - - - - .x = SingleCellExperiment(assay = list( my_assay ) |> set_names(names(assays(.x))[1]), colData = cd) - } - - - # # TEMPORARY FOR SOME REASON THE MIN COUNTS IS NOT 0 FOR SOME SAMPLES - # .x = HPCell:::check_if_assay_minimum_count_is_zero_and_correct_TEMPORARY(.x, assays(.x) |> names() |> _[1], subset_up_to_number_of_cells = 100) - - print("start ranking") - - # CALCULATE rank - rank_assay = - .x |> - assay() |> - - # This because some datasets are still > 1M cells - process_matrix_in_slices( - paste(c(cache_directory, "/", .y, "_rank_matrix.HDF5Array"), collapse = ""), - paste(c(cache_directory, "/", .y, "_rank_matrix_temp.HDF5Array"), collapse = ""), - chunk_size = 1000 - ) - - print("creating SCE") - - .x = SingleCellExperiment(assay = list( rank = rank_assay), colData = colData(.x)) - - print("saving") - - .x |> save_experiment_data(glue("{cache_directory}/{.y}")) - - # Delete the temp file - file.remove(paste(c(cache_directory, "/", .y, "_rank_matrix_temp.HDF5Array"), collapse = "")) - - return(TRUE) # Indicate successful saving - } - - # Because they have an inconsistent failure. If I start the pipeline again they might work. Strange. - insistent_save_rank_per_cell <- purrr::insistently(save_rank_per_cell, rate = purrr::rate_delay(pause = 60, max_times = 3), quiet = FALSE) - - - save_anndata_sct = function(dataset_id_sce, cache_directory){ - - dir.create(cache_directory, showWarnings = FALSE, recursive = TRUE) - - if (is.null(dataset_id_sce)) return(NULL) - - .x = dataset_id_sce |> pull(sct) |> _[[1]] - - # Fix: check is.null BEFORE ncol() to avoid `argument is of length zero` - if (is.null(.x) || ncol(.x) == 0) return(NULL) - - .y = dataset_id_sce |> pull(file_id_cellNexus_single_cell) |> _[[1]] |> str_remove("\\.h5ad") - - .x |> assays() |> names() = "sct" - - # Wrap save with explicit error logging so the real cause is visible. Strange it shouldnt fail, which passed in debug mode - tryCatch( - .x |> save_experiment_data(glue("{cache_directory}/{.y}")), - error = function(e) { - message(glue::glue("[save_anndata_sct] FAILED for {.y}: {conditionMessage(e)}")) - stop(e) - } - ) - - return(TRUE) - } - - # Because they have an inconsistent failure. If I start the pipeline again they might work. Strange. - insistent_save_anndata_sct <- purrr::insistently(save_anndata_sct, rate = purrr::rate_delay(pause = 60, max_times = 3), quiet = FALSE) - - - cbind_sce_by_dataset_id = function(target_name_grouped_by_dataset_id, file_id_db_file, cell_id_dict, my_store){ - - my_dataset_id = unique(target_name_grouped_by_dataset_id$dataset_id) - my_file_id = unique(target_name_grouped_by_dataset_id$file_id_cellNexus_single_cell) - - file_id_db = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue("SELECT * FROM read_parquet('{file_id_db_file}')")) - ) |> - filter(dataset_id == my_dataset_id) |> - select(cell_id, sample_id, dataset_id, file_id_cellNexus_single_cell) - - file_id_db = - target_name_grouped_by_dataset_id |> - left_join(file_id_db, copy = TRUE) - - - dataset_cell_dict = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue("SELECT * FROM read_parquet('{cell_id_dict}')")) - ) |> - filter(file_id_cellNexus_single_cell == my_file_id) - - file_id_db = - file_id_db |> - left_join(dataset_cell_dict, by = c("file_id_cellNexus_single_cell", "cell_id" ), copy=T ) - - # Parallelise - cores = as.numeric(Sys.getenv("SLURM_CPUS_PER_TASK", unset = 1))-1 - # Respect R CMD CHECK core limit if set - if (nzchar(Sys.getenv("_R_CHECK_LIMIT_CORES_"))) { - cores <- min(cores, 2L) - } - # MulticoreParam need to have enough memory to proceed, otherwise reducer error - bp <- MulticoreParam(workers = cores , progressbar = TRUE) # Adjust the number of workers as needed - - # Begin processing the data pipeline with the initial dataset 'target_name_grouped_by_dataset_id' - sce_df = - file_id_db |> - nest(cells = c(cell_id, new_cell_id)) |> - # Read raw data for each 'target_name' and store it in a new column 'sce' - mutate( - sce = bplapply( - sce_target_name, - FUN = function(x) { - tar_read_raw(x, store = my_store) |> - select(.cell, donor_id, dataset_id, sample_id, cell_type) |> - mutate(sample_id = as.factor(sample_id)) # lighter - }, # Read the raw SingleCellExperiment object - BPPARAM = bp # Use the defined parallel backend - )) |> - # This should not be needed, but there are some data sets with zero cells - filter(!map_lgl(sce, is.null)) |> - mutate(sce = map2(sce, cells, ~ { - - cell_map <- setNames(.y$new_cell_id, .y$cell_id) - - .x |> filter(.cell %in% names(cell_map)) %>% - { - colnames(.) <- cell_map[colnames(.)] - . - } - - }, .progress = TRUE)) - - if(nrow(sce_df) == 0) { - warning("this chunk has no rows for somereason.") - return(NULL) - } - - sce_df |> - - # Group the data by 'dataset_id' and 'tar_group' for further summarization - mutate(sce = map(sce, ~ SingleCellExperiment(assay = assays(.x), colData = colData(.x)) )) |> - - # Combine all 'sce' objects within each group into a single 'sce' object - group_by(file_id_cellNexus_single_cell) |> - summarise( sce = list(do.call(cbind, args = sce) ), - # A step to check missing cells - cells = list(do.call(rbind, args = cells))) - } - - cbind_sct_by_dataset_id = function(target_name_grouped_by_dataset_id, file_id_db_file, cell_id_dict, my_store){ - - my_dataset_id = unique(target_name_grouped_by_dataset_id$dataset_id) - my_file_id = unique(target_name_grouped_by_dataset_id$file_id_cellNexus_single_cell) - - file_id_db = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue("SELECT * FROM read_parquet('{file_id_db_file}')")) - ) |> - filter(dataset_id == my_dataset_id) |> - select(cell_id, sample_id, dataset_id, file_id_cellNexus_single_cell) - - file_id_db = - target_name_grouped_by_dataset_id |> - left_join(file_id_db, copy = TRUE) - - - dataset_cell_dict = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue("SELECT * FROM read_parquet('{cell_id_dict}')")) - ) |> - filter(file_id_cellNexus_single_cell == my_file_id) - - file_id_db = - file_id_db |> - left_join(dataset_cell_dict, by = c("file_id_cellNexus_single_cell", "cell_id"), copy=T ) - - # Parallelise - cores = as.numeric(Sys.getenv("SLURM_CPUS_PER_TASK", unset = 1)) -1 - # Respect R CMD CHECK core limit if set - if (nzchar(Sys.getenv("_R_CHECK_LIMIT_CORES_"))) { - cores <- min(cores, 2L) - } - # MulticoreParam need to have enough memory to proceed, otherwise reducer error - bp <- MulticoreParam(workers = cores , progressbar = TRUE) # Adjust the number of workers as needed - - # Begin processing the data pipeline with the initial dataset 'target_name_grouped_by_dataset_id' - sct_df = file_id_db |> - nest(cells = c(cell_id, new_cell_id)) %>% - # Step 1: Read raw data for each 'target_name' and store it in a new column 'sce' - mutate( - sct = bplapply( - sct_target_name, - FUN = function(x) { - if (is.na(x)) { - return(NULL) # because cant get sample_id and dataset_id from NULL sct_matrix - } - - tar_read_raw(x, store = my_store) |> - select(.cell, donor_id, dataset_id, sample_id, cell_type) |> - mutate(sample_id = as.factor(sample_id)) - - }, # Read the raw SingleCellExperiment object - BPPARAM = bp # Use the defined parallel backend - )) |> - # This should not be needed, but there are some data sets with zero cells - filter(!map_lgl(sct, is.null)) |> - mutate(sct = map2(sct, cells, ~ { - - cell_map <- setNames(.y$new_cell_id, .y$cell_id) - - .x |> filter(.cell %in% names(cell_map)) %>% - { - colnames(.) <- cell_map[colnames(.)] - . - } - - }, .progress = TRUE)) - - if(nrow(sct_df) == 0) { - warning("this chunk has no rows for somereason.") - return(NULL) - } - - sct_df |> - mutate( - sct = map(sct, \(x) { - if (is.null(x)) return(NULL) - SingleCellExperiment(assays = assays(x), colData = colData(x)) - }) - ) |> - group_by(file_id_cellNexus_single_cell) |> - summarise( - sct = { - scts <- compact(sct) # drop NULLs inside each group - - list( - if (length(scts) == 0) { - NULL - } else { - - # A few big samples do not return all features because it reached R limit 2^31-1 in SCTransform - common_genes <- cellNexus:::check_gene_overlap(scts) - - # subset to intersection genes (and keep same order across objects) - scts2 <- map(scts, \(z) z[common_genes, , drop = FALSE]) - - do.call(SummarizedExperiment::cbind, scts2) - } - ) - }, - cells = list(do.call(rbind, cells)), - .groups = "drop" - ) - } - - # Because they have an inconsistent failure. If I start the pipeline again they might work. Strange. - insistent_cbind_sct_by_dataset_id <- purrr::insistently(cbind_sct_by_dataset_id, rate = purrr::rate_delay(pause = 60, max_times = 3), quiet = FALSE) - - - get_dataset_id = function(target_name, my_store){ - # Try reading the target safely (for some failing targets) - sce = tryCatch( - tar_read_raw(target_name, store = my_store), - error = function(e) return(NULL) - ) - - # Still need to catch target_name - if(sce |> is.null()) return(tibble(sample_id = NA_character_, - dataset_id= NA_character_, - target_name= !!target_name)) - - sce |> - - distinct(sample_id, dataset_id) |> mutate(target_name = !!target_name) - } - - create_chunks_for_reading_and_saving = function(dataset_id_sample_id, cell_metadata){ - - # Solve sample_id mismatches because some end with .h5ad suffix while others dont - dataset_id_sample_id |> - - left_join( - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue("SELECT * FROM read_parquet('{cell_metadata}')")) - ) |> - distinct(dataset_id, sample_id, sample_chunk, cell_chunk, file_id_cellNexus_single_cell) |> - as_tibble(), - copy=T - ) - } - - - cbind_sce_by_dataset_id_get_missing_cells = function(dataset_id_sce){ - - dataset_id_sce |> - mutate( - missing_cells = map2( - sce, - cells, - ~{ - cells_in_sce <- .x |> colnames() |> sort() - - cells_in_query <- .y$new_cell_id |> unique() |> sort() - - # Find differences - tibble(cell_id = setdiff(cells_in_query, cells_in_sce)) - } - ) - ) |> - select(file_id_cellNexus_single_cell, missing_cells) - - } - - - list( - - # The input DO NOT DELETE - tar_target(my_store, "/vast/scratch/users/shen.m/cellNexus/2024-07-01/process_samples_hpcell_target_store", deployment = "main"), # MODIFY HERE: HPCell targets store to read SCEs from - tar_target(cache_directory, "/vast/scratch/users/shen.m/cellNexus/cellxgene_2024/0.2.1", deployment = "main"), # MODIFY HERE: output cache directory for saved anndata files - tar_target( - cell_metadata, - "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_1_mengyuan.parquet", # MODIFY HERE: final metadata parquet (should match the COPY TO output above) - packages = c( "arrow","dplyr","duckdb") - - ), - - tar_target( - cell_id_dict, - "/vast/projects//cellxgene_curated/metadata_cellxgene_mengyuan/file_id_cell_id_dict_v1_1_1_Jul_2024.parquet", # MODIFY HERE: cell_id dictionary parquet - packages = c( "arrow","dplyr","duckdb") - ), - - # pre-calculated counts - tar_target( - target_name, - tar_meta( - starts_with("sce_transformed_"), - store = my_store) |> - filter(type=="branch") |> - pull(name), - deployment = "main" - ), - tar_target( - dataset_id_sample_id, - get_dataset_id(target_name, my_store), - packages = "tidySingleCellExperiment", - pattern = map(target_name), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_5_minimal") - ) - ), - - # pre-calculated sct - tar_target( - sct_target_name, - tar_meta( - starts_with("sct_matrix_"), - store = my_store) |> - filter(type=="branch") |> - pull(name), - deployment = "main" - ), - tar_target( - sct_dataset_id_sample_id, - get_dataset_id(sct_target_name, my_store), - packages = "tidySingleCellExperiment", - pattern = map(sct_target_name), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_5_minimal") - ) - ), - - # join - tar_target( - dataset_id_sample_id_target_names, - dataset_id_sample_id |> left_join(sct_dataset_id_sample_id, by = c("sample_id", "dataset_id"), copy=T) |> - dplyr::rename(sce_target_name = target_name.x, - sct_target_name = target_name.y), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_5_minimal") - ) - ), - - tar_target( - target_name_grouped_by_dataset_id, - create_chunks_for_reading_and_saving(dataset_id_sample_id_target_names, cell_metadata) |> - - # # FOR TESTING PURPOSE ONLY - # filter(file_id_cellNexus_single_cell %in% c("e8291e33fc98bc21728255f0a5669015___1.h5ad", - # "4ceb75a970ccdc9aaa9f9e91b931292f___1.h5ad")) |> - - group_by(dataset_id, sample_chunk, cell_chunk, file_id_cellNexus_single_cell) |> - tar_group(), - iteration = "group", - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_5_minimal") - ), - packages = c("arrow", "duckdb", "dplyr", "glue", "targets") - - ), - - tar_target( - dataset_id_sce, - cbind_sce_by_dataset_id(target_name_grouped_by_dataset_id, cell_metadata, cell_id_dict, my_store = my_store), - pattern = map(target_name_grouped_by_dataset_id), - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "digest", "scater", "HDF5Array", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly"), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_40") - ) - ), - - tar_target( - dataset_id_sct, - cbind_sct_by_dataset_id(target_name_grouped_by_dataset_id, cell_metadata, cell_id_dict, my_store = my_store), - pattern = map(target_name_grouped_by_dataset_id), - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "digest", "scater", "HDF5Array", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly"), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_120") - ) - ), - - # This target was run for retrieving missing cells analysis only - tar_target( - missing_cells_tbl, - cbind_sce_by_dataset_id_get_missing_cells(dataset_id_sce), - pattern = map(dataset_id_sce), - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly", "purrr"), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_20") - ) - ), - - - tar_target( - save_anndata, - insistent_save_anndata(dataset_id_sce, paste0(cache_directory, "/counts")), - pattern = map(dataset_id_sce), - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "HPCell", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly"), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_20") - ) - ), - - tar_target( - saved_dataset_cpm, - insistent_save_anndata_cpm(dataset_id_sce, paste0(cache_directory, "/cpm")), - pattern = map(dataset_id_sce), - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "HPCell", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly"), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_20") - ) - ), - - tar_target( - saved_dataset_rank, - insistent_save_rank_per_cell(dataset_id_sce, paste0(cache_directory, "/rank")), - pattern = map(dataset_id_sce), - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "HPCell", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly", "HDF5Array"), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_20") - ) - ), - - tar_target( - saved_sct, - save_anndata_sct(dataset_id_sct, paste0(cache_directory, "/sct")), - pattern = map(dataset_id_sct), - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "HPCell", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly", "HDF5Array"), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_20") - ) - ) - ) - -}, script = paste0(store_file_cellNexus, "_target_script.R"), ask = FALSE) - -job::job({ - - tar_make( - script = paste0(store_file_cellNexus, "_target_script.R"), - store = store_file_cellNexus, - reporter = "summary" #, callr_function = NULL - ) - -}) - -missing_cells_tbl = tar_read(missing_cells_tbl, store = store_file_cellNexus) |> - unnest(missing_cells) - -#missing_cells_tbl |> write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cells_to_remove_in_metadata_Jul_2024.parquet") -missing_cells_tbl <- read_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cells_to_remove_in_metadata_Jul_2024.parquet") - -filtered_cell_metadata = cell_metadata |> anti_join(missing_cells_tbl, by = c("file_id_cellNexus_single_cell", - "new_cell_id" = "cell_id"), copy = T) - -filtered_cell_metadata |> - collect() |> - arrow::write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_1_filtered_missing_cells_mengyuan.parquet", - compression = "zstd") # MODIFY HERE: output parquet after filtering missing cells - diff --git a/dev/cellnexus-2024-scripts/step6_supp_dataset_cell_map.R b/dev/cellnexus-2024-scripts/step6_supp_dataset_cell_map.R deleted file mode 100644 index db2630c..0000000 --- a/dev/cellnexus-2024-scripts/step6_supp_dataset_cell_map.R +++ /dev/null @@ -1,105 +0,0 @@ -library(targets) -store = "/vast/scratch/users/shen.m/cellnexus_dataset_cell_map_Jul_2024_v1_2_1_target_store" -tar_script({ - library(dplyr) - library(magrittr) - library(tibble) - library(targets) - library(tarchetypes) - library(crew) - library(crew.cluster) - - tar_option_set( - memory = "transient", - garbage_collection = 100, - storage = "worker", - retrieval = "worker", - error = "continue", - cue = tar_cue(mode = "never"), - - workspace_on_error = TRUE, - controller = crew_controller_group( - list( - crew_controller_slurm( - name = "elastic", - workers = 300, - tasks_max = 20, - seconds_idle = 30, - crashes_error = 10, - options_cluster = crew_options_slurm( - memory_gigabytes_required = c(25, 35, 40, 80, 160), - cpus_per_task = c(2, 2, 5, 10, 20), - time_minutes = c(30, 30, 30, 60*4, 60*24), - verbose = T - ) - ) - ) - ), - trust_object_timestamps = TRUE - ) - - get_unique_file_ids <- function(cell_metadata){ - tbl(dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue::glue("SELECT * FROM read_parquet('{cell_metadata}')"))) |> - distinct(file_id_cellNexus_single_cell) |> pull() - } - - create_file_id_cell_id_dict <- function(cell_metadata, file_id) { - - - tbl(dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue::glue("SELECT * FROM read_parquet('{cell_metadata}')"))) |> - filter(file_id_cellNexus_single_cell == file_id) |> - dbplyr::window_order(cell_id) |> - mutate(cell_index = row_number()) |> - select(cell_id, file_id_cellNexus_single_cell, - new_cell_id = cell_index) |> - collect() - } - - list( - tar_target(cell_metadata , "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_0_mengyuan.parquet", - deployment = "main"), - tar_target( - unique_file_ids, - # TESTING PURPOSE ONLY - # c("3cef5b6aa0f5772485bb710f71e69456___1.h5ad", - # "cd2caa6de850f73af4ca78a2ea307dd4___1.h5ad") - get_unique_file_ids(cell_metadata) - # |> head(2) - , - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly", "HDF5Array") - # resources = tar_resources( - # crew = tar_resources_crew(controller = "elastic") - # ) - ), - tar_target( - file_id_cell_id_dict, - create_file_id_cell_id_dict(cell_metadata, unique_file_ids), - pattern = map(unique_file_ids), - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly", "HDF5Array") - # resources = tar_resources( - # crew = tar_resources_crew(controller = "elastic") - # ) - ) - ) - -}, script = paste0(store, "_target_script.R"), ask = FALSE) - - -job::job({ - - tar_make( - script = paste0(store, "_target_script.R"), - store = store, - reporter = "summary" - ) - -}) - -file_id_cell_id_dict = tar_read(file_id_cell_id_dict, store = store) -file_id_cell_id_dict |> arrow::write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/file_id_cell_id_dict_v1_1_1_Jul_2024.parquet", - compression = "zstd") -rm(file_id_cell_id_dict) -gc() - diff --git a/dev/cellnexus-2024-scripts/step7_prepare_local_cache_splitting_du_dataset_and_cell_type_mengyuan_version.R b/dev/cellnexus-2024-scripts/step7_prepare_local_cache_splitting_du_dataset_and_cell_type_mengyuan_version.R deleted file mode 100644 index cf5a674..0000000 --- a/dev/cellnexus-2024-scripts/step7_prepare_local_cache_splitting_du_dataset_and_cell_type_mengyuan_version.R +++ /dev/null @@ -1,987 +0,0 @@ -# Step6 -# Group samples by dataset_id, cell_type - -# This script sets up a robust and scalable data processing pipeline for single-cell RNA sequencing (scRNA-seq) datasets using the targets package in R, which facilitates reproducible and efficient workflows. Specifically, the code orchestrates the ingestion and preprocessing of multiple SingleCellExperiment objects corresponding to different datasets (dataset_id) and targets (target_name). It leverages high-performance computing resources through the crew package, configuring multiple SLURM-based controllers (tier_1 to tier_4) to handle varying computational loads efficiently. -# -# The pipeline performs several key steps: -# -# 1. Data Retrieval: It reads raw SingleCellExperiment objects for each target, ensuring that only successfully loaded data proceeds further. -# 2. Normalization: Calculates Counts Per Million (CPM) for each cell to normalize gene expression levels across cells and samples. -# 3. Data Aggregation: Groups the data by dataset_id and tar_group, then combines the SingleCellExperiment objects within each group into a single object, effectively consolidating the data for each dataset. -# 4. Metadata Integration: Joins additional metadata, such as cell types, by connecting to a DuckDB database and fetching relevant information from a Parquet file. This enriches the single-cell data with essential annotations. -# 5. Cell Type Segmentation: Splits the combined SingleCellExperiment objects into separate objects based on cell_type, facilitating downstream analyses that are specific to each cell type. -# 6. Data Saving with Error Handling: Generates unique identifiers for each cell type within a dataset and saves both the raw counts and CPM-normalized data to specified directories. It includes special handling for cases where a cell type has only one cell, duplicating the data to prevent errors during the saving process. -# -# By integrating targets, crew, and various data manipulation packages (dplyr, tidyverse, SingleCellExperiment), this script ensures that large-scale scRNA-seq data processing is efficient, reproducible, and capable of leveraging parallel computing resources. It is designed to handle edge cases gracefully and provides a clear framework for preprocessing scRNA-seq data, which is essential for subsequent analyses such as clustering, differential expression, and cell type identification. - - -library(arrow) -library(dplyr) -library(duckdb) - -job::job({ - - get_file_ids = function(cell_annotation ){ - sample_chunk_df = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue::glue("SELECT * FROM read_parquet('{cell_annotation}')")) - ) |> - # Define chunks - dplyr::count(dataset_id, sample_id, name = "cell_count") |> # Ensure unique dataset_id and sample_id combinations - distinct(dataset_id, sample_id, cell_count) |> # Ensure unique dataset_id and sample_id combinations - group_by(dataset_id) |> - dbplyr::window_order(dataset_id, cell_count, sample_id) |> # Ensure order. Note: order cell_count only is not enough because it needs a secondary tie-breaker - mutate(sample_index = row_number()) |> # Create sequential index within each dataset - mutate(sample_chunk = (sample_index - 1) %/% 1000 + 1) |> # Assign chunks (up to 1000 samples per chunk) - mutate(sample_pseudobulk_chunk = (sample_index - 1) %/% 250 + 1) |> # Max combination of dataset_id, sample_pseudobulk_chunk and file_id_pseudobulk up to 10000 - mutate(cell_chunk = cumsum(cell_count) %/% 100000 + 1) |> # max 20K cells per sample - ungroup() - - # Test whether cell_chunk and sample_chunk are unique for this sample - run_chunk_once <- function(column_name, id) { - sample_chunk_df |> filter(sample_id == id) |> pull(!!column_name) - } - - sample_chunk_results <- replicate(20, run_chunk_once("sample_chunk", "d6e942a09a140ee8bb6f0c3da8defea4___exp7-human-150well."), simplify = FALSE) - sample_chunk_identical <- all(sapply(sample_chunk_results[-1], function(x) identical(x, sample_chunk_results[[1]]))) - if (!sample_chunk_identical) { - stop("Inconsistent sample chunk value was generated in multiple runs, this will lead to file id changes") - } - - cell_chunk_results <- replicate(20, run_chunk_once("cell_chunk", "d6e942a09a140ee8bb6f0c3da8defea4___exp7-human-150well."), simplify = FALSE) - cell_chunk_identical <- all(sapply(cell_chunk_results[-1], function(x) identical(x, cell_chunk_results[[1]]))) - if (!cell_chunk_identical) { - stop("Inconsistent cell chunk value was generated in multiple runs, this will lead to file id changes") - } - - - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue::glue("SELECT * FROM read_parquet('{cell_annotation}')")) - ) |> - # Cells in cell_annotation could be more than cells in cell_consensus. In order to avoid NA happens in cell_consensus cell_type column - mutate(cell_type_unified_ensemble = ifelse(cell_type_unified_ensemble |> is.na(), - "Unknown", - cell_type_unified_ensemble)) |> - - left_join(sample_chunk_df |> select(dataset_id, sample_chunk, sample_pseudobulk_chunk, cell_chunk, sample_id), copy=TRUE) |> - - # Define chunks - group_by(dataset_id, sample_chunk, cell_chunk, sample_pseudobulk_chunk, cell_type, sample_id) |> - summarise(cell_count = n(), .groups = "drop") |> - group_by(dataset_id, sample_chunk, cell_chunk, cell_type) |> - dbplyr::window_order(desc(cell_count), sample_id) |> # Important! - mutate(chunk = cumsum(cell_count) %/% 20000 + 1) |> # max 20K cells per sample - ungroup() |> - as_tibble() |> - - # Single cell file ID - mutate(file_id_cellNexus_single_cell = - glue::glue("{dataset_id}___{sample_chunk}___{cell_chunk}___{cell_type}") |> - sapply(digest::digest) |> - paste0("___", chunk, ".h5ad") - ) |> - - # Pseudobulk file id - mutate(file_id_cellNexus_pseudobulk = - glue::glue("{dataset_id}___{sample_pseudobulk_chunk}") |> - sapply(digest::digest) |> - paste0("___", chunk, ".h5ad")) - - } - - get_file_ids( - "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_annotation_2024_Jul.parquet" # MODIFY HERE: input cell annotation parquet - ) |> - write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/file_id_cellNexus_single_cell_2024_Jul.parquet") # MODIFY HERE: output file_id parquet - - gc() - - con <- dbConnect(duckdb::duckdb(), dbdir = ":memory:") - - dir.create("/vast/scratch/users/shen.m/duckdb_tmp", showWarnings = FALSE) # MODIFY HERE: duckdb temp directory - - DBI::dbExecute( - con, - "SET temp_directory='/vast/scratch/users/shen.m/duckdb_tmp';" # MODIFY HERE: duckdb temp directory (must match dir.create above) - ) - - # Create a view for cell_annotation in DuckDB - # MODIFY HERE: cell_metadata parquet path inside the SQL string below - dbExecute(con, " - CREATE VIEW cell_metadata AS - SELECT - CONCAT(cell_, '___', dataset_id) AS cell_, - * EXCLUDE (cell_, dataset_id_1, X_umap1, X_umap2, sample_placeholder, cell_type) -- drop original cell_ and dataset_id_1 - FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata.parquet') -") - - # MODIFY HERE: cell_annotation parquet path inside the SQL string below - dbExecute(con, " - CREATE VIEW hpcell_output_metadata AS - SELECT * EXCLUDE ( - observation_joinid, - cell_type_ontology_term_id, - assay, - donor_id, - is_primary_data, - self_reported_ethnicity, - tissue, - azimuth, - blueprint, - monaco, - subsets_Mito_sum, - subsets_Mito_detected, - ensemble_joinid, - cell_type_unified, - data_driven_ensemble, - observation_originalid -) - - FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_annotation_2024_Jul.parquet') -") - -# MODIFY HERE: file_id parquet path inside the SQL string below (should match the write_parquet output above) -dbExecute(con, " - CREATE VIEW file_id_cellNexus_single_cell AS - SELECT - dataset_id, - sample_chunk, - cell_chunk, - sample_pseudobulk_chunk, - cell_type, - sample_id, - file_id_cellNexus_single_cell, - file_id_cellNexus_pseudobulk - FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/file_id_cellNexus_single_cell_2024_Jul.parquet') -") - -# MODIFY HERE: transformation data frame -dbExecute(con, " - CREATE VIEW sample_distribution_method_tbl AS - SELECT - sample_2 AS sample_id, - count_upper_bound, - feature_thresh AS nfeature_expressed_thresh, - method_to_apply AS inverse_transform - FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/sliced_sample_tbl_2024_Jul.parquet') -") - -# Perform the left join and save to Parquet -copy_query <- " - COPY ( - SELECT - cell_metadata.cell_ AS cell_id, -- Rename cell_ to cell_id - COALESCE(hpcell_output_metadata.alive, FALSE) AS alive, -- Set alive column NULL to FALSE - cell_metadata.* EXCLUDE (cell_), -- drop cell_ since it's already aliased as cell_id - hpcell_output_metadata.* EXCLUDE (cell_, dataset_id, sample_id, alive), -- Deduplicate join keys, and aliased column - file_id_cellNexus_single_cell.* EXCLUDE (sample_id, dataset_id, cell_type), -- Deduplicate join keys - sample_distribution_method_tbl.* EXCLUDE (sample_id) -- Deduplicate join keys - FROM cell_metadata - - LEFT JOIN hpcell_output_metadata - ON hpcell_output_metadata.cell_ = cell_metadata.cell_ - AND hpcell_output_metadata.dataset_id = cell_metadata.dataset_id - - LEFT JOIN file_id_cellNexus_single_cell - ON file_id_cellNexus_single_cell.sample_id = hpcell_output_metadata.sample_id - AND file_id_cellNexus_single_cell.dataset_id = hpcell_output_metadata.dataset_id - AND file_id_cellNexus_single_cell.cell_type = hpcell_output_metadata.cell_type - - LEFT JOIN sample_distribution_method_tbl - ON sample_distribution_method_tbl.sample_id = cell_metadata.sample_id - - WHERE cell_metadata.dataset_id NOT IN ('99950e99-2758-41d2-b2c9-643edcdf6d82', '9fcb0b73-c734-40a5-be9c-ace7eea401c9') -- (THESE TWO DATASETS DOESNT contain meaningful data - no observation_joinid etc), thus was excluded in the final metadata. - - ) TO '/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_0_mengyuan.parquet' -- MODIFY HERE: output merged metadata parquet (v1_2_2) - (FORMAT PARQUET, COMPRESSION 'gzip'); -" - -# Execute the final query to write the result to a Parquet file -dbExecute(con, copy_query) - -# Disconnect from the database -dbDisconnect(con, shutdown = TRUE) - -print("Done.") -}) - -# We decided to make cell_id lighter without re-run everything in HPCell pipeline. Here to swap cell_id in the metadata -# cell_map is processed in a separate target script step6_supp_dataset_cell_map.R -job::job({ - con <- dbConnect(duckdb::duckdb(), dbdir = ":memory:") - - # Create a view for cell_annotation in DuckDB - # MODIFY HERE: v1_2_2 merged metadata parquet path inside the SQL string below (should match the COPY TO output above) - dbExecute(con, " - CREATE VIEW cell_metadata AS - SELECT * - FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_0_mengyuan.parquet') -") - - # MODIFY HERE: cell_id dictionary parquet path inside the SQL string below - dbExecute(con, " - CREATE VIEW cell_map AS - SELECT * - FROM read_parquet('/vast/projects//cellxgene_curated/metadata_cellxgene_mengyuan/file_id_cell_id_dict_v1_1_1_Jul_2024.parquet') -") - - # Perform the left join and save to Parquet - copy_query <- " - COPY ( - SELECT cell_metadata.*, - cell_map.* EXCLUDE (cell_id, file_id_cellNexus_single_cell) - FROM cell_metadata - - LEFT JOIN cell_map - ON cell_metadata.cell_id = cell_map.cell_id - AND cell_metadata.file_id_cellNexus_single_cell = cell_map.file_id_cellNexus_single_cell - - ) TO '/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_1_mengyuan.parquet' -- MODIFY HERE: output final metadata parquet with new cell IDs (v1_3_2) - (FORMAT PARQUET, COMPRESSION 'gzip'); -" - - # Execute the final query to write the result to a Parquet file - dbExecute(con, copy_query) - - # Disconnect from the database - dbDisconnect(con, shutdown = TRUE) - - print("Done.") - - -}) - - - -# MODIFY HERE: final metadata parquet path used for the targets pipeline (should match the COPY TO output above) -cell_metadata = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql("SELECT * FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_1_mengyuan.parquet')") - ) - -library(targets) -library(tidyverse) -store_file_cellNexus = "/vast/scratch/users/shen.m/targets_prepare_database_split_datasets_chunked_1_5_0_single_cell" # MODIFY HERE: targets store directory for this pipeline - -tar_script({ - library(dplyr) - library(magrittr) - library(tibble) - library(targets) - library(tarchetypes) - library(crew) - library(crew.cluster) - - # Helper (optional) to avoid repetition - new_elastic <- function(name, mem_gb, time_min, workers, crashes_max, cpus_per_task = 2, backup = NULL) { - crew_controller_slurm( - name = name, - workers = workers, - crashes_max = crashes_max, - seconds_idle = 30, - options_cluster = crew_options_slurm( - memory_gigabytes_required = mem_gb, - cpus_per_task = cpus_per_task, - time_minutes = time_min - ), - backup = backup - ) - } - - # Small → large, with fallbacks to the next size up - # elastic_160 <- new_elastic("elastic_160", 160, 60 * 24, workers = 8, crashes_max = 2) - # elastic_120 <- new_elastic("elastic_120", 120, 60 * 4, workers = 16, crashes_max = 1, cpus_per_task = 8, backup = elastic_160) - # elastic_80 <- new_elastic("elastic_80", 80, 60 * 4, workers = 24, crashes_max = 1, cpus_per_task = 8, backup = elastic_120) - # elastic_40 <- new_elastic("elastic_40", 40, 60 * 4, workers = 32, crashes_max = 1, cpus_per_task = 8, backup = elastic_80) - # elastic_20 <- new_elastic("elastic_20", 20, 60 * 4, workers = 48, crashes_max = 1, cpus_per_task = 8, backup = elastic_40) - # elastic_10 <- new_elastic("elastic_10", 10, 60 * 4, workers = 150, crashes_max = 6, cpus_per_task = 8, backup = elastic_20) - # - # elastic_5_minimal <- new_elastic("elastic_5_minimal", 5, 60 * 4, workers = 300, crashes_max = 6, cpus_per_task = 8, backup = elastic_10) - # - elastic_160 <- new_elastic("elastic_160", 160, 60 * 1, workers = 8, crashes_max = 2) - elastic_120 <- new_elastic("elastic_120", 120, 60 * 1, workers = 16, crashes_max = 1, cpus_per_task = 8, backup = elastic_160) - elastic_80 <- new_elastic("elastic_80", 80, 60 * 1, workers = 24, crashes_max = 1, cpus_per_task = 8, backup = elastic_120) - elastic_40 <- new_elastic("elastic_40", 40, 60 * 1, workers = 32, crashes_max = 1, cpus_per_task = 8, backup = elastic_80) - elastic_20 <- new_elastic("elastic_20", 20, 60 * 1, workers = 48, crashes_max = 1, cpus_per_task = 8, backup = elastic_40) - elastic_10 <- new_elastic("elastic_10", 10, 60 * 1, workers = 150, crashes_max = 6, cpus_per_task = 8, backup = elastic_20) - - elastic_5_minimal <- new_elastic("elastic_5_minimal", 5, 60 * 1, workers = 300, crashes_max = 6, cpus_per_task = 8, backup = elastic_10) - - - # Group for targets (small → large) - controllers <- crew_controller_group( - elastic_10, elastic_20, elastic_40, elastic_80, elastic_120, elastic_160, elastic_5_minimal - ) - tar_option_set( - memory = "transient", - garbage_collection = 100, - storage = "worker", - retrieval = "worker", - error = "continue", - cue = tar_cue(mode = "never"), - format = "qs", - #debug = "dataset_id_sct_ea377f6e2d0ae2b7", - workspace_on_error = TRUE, - controller = controllers, - trust_object_timestamps = TRUE, - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_5_minimal") - ) - ) - - save_anndata = function(dataset_id_sce, cache_directory){ - - dir.create(cache_directory, showWarnings = FALSE, recursive = TRUE) - - .x = dataset_id_sce |> pull(sce) |> _[[1]] - .y = dataset_id_sce |> pull(file_id_cellNexus_single_cell) |> _[[1]] |> str_remove("\\.h5ad") - - .x |> assays() |> names() = "counts" - - # Save the experiment data to the specified counts cache directory - .x |> save_experiment_data(glue("{cache_directory}/{.y}")) - - return(TRUE) # Indicate successful saving - - - } - - # Because they have an inconsistent failure. If I start the pipeline again they might work. Strange. - insistent_save_anndata <- purrr::insistently(save_anndata, rate = purrr::rate_delay(pause = 60, max_times = 3), quiet = FALSE) - - save_anndata_cpm = function(dataset_id_sce, cache_directory){ - - dir.create(cache_directory, showWarnings = FALSE, recursive = TRUE) - - # # Parallelise - dataset_id_sce |> - purrr::transpose() |> - lapply( - FUN = function(x) { - - .x = x[[2]] - .y = x[[1]] |> str_remove("\\.h5ad") - - # Check if the 'sce' has only one cell (column) - if(ncol(assay(.x)) == 1) { - - # Duplicate the assay to prevent saving errors due to single-column matrices - my_assay = cbind(assay(.x), assay(.x)) - # Rename the second column to distinguish it - colnames(my_assay)[2] = paste0("DUMMY", "___", colnames(my_assay)[2]) - - cd = colData(.x) - cd = cd |> rbind(cd) - rownames(cd)[2] = paste0("DUMMY", "___", rownames(cd)[2]) - - - - .x = SingleCellExperiment(assay = list( my_assay ) |> set_names(names(assays(.x))[1]), colData = cd) - } - - - # # TEMPORARY FOR SOME REASON THE MIN COUNTS IS NOT 0 FOR SOME SAMPLES - # .x = HPCell:::check_if_assay_minimum_count_is_zero_and_correct_TEMPORARY(.x, assays(.x) |> names() |> _[1], subset_up_to_number_of_cells = 100) - - # CALCULATE CPM - .x = SingleCellExperiment(assay = list( cpm = calculateCPM(.x, assay.type = names(assays(.x))[1])), colData = colData(.x)) - - # Save the experiment data to the specified counts cache directory - .x |> save_experiment_data(glue("{cache_directory}/{.y}")) - - return(TRUE) # Indicate successful saving - } - - ) - - return("saved") - - } - - # Because they have an inconsistent failure. If I start the pipeline again they might work. Strange. - insistent_save_anndata_cpm <- purrr::insistently(save_anndata_cpm, rate = purrr::rate_delay(pause = 60, max_times = 3), quiet = FALSE) - - - # Function to process matrix in vertical slices - process_matrix_in_slices <- function(h5_matrix, output_filepath, output_filepath_temp, chunk_size = 1000) { - # Load the HDF5 matrix - n_rows <- dim(h5_matrix)[1] - n_cols <- dim(h5_matrix)[2] - - if (file.exists(output_filepath)) { - file.remove(output_filepath) - cat("Existing output file removed.\n") - } - if (file.exists(output_filepath_temp)) { - file.remove(output_filepath_temp) - cat("Existing output file removed.\n") - } - - # Create an empty list to hold the slices - slice_list <- list() - - # Loop through the matrix in chunks - for (start_col in seq(1, n_cols, by = chunk_size)) { - end_col <- min(start_col + chunk_size - 1, n_cols) - cat("Processing columns", start_col, "to", end_col, "\n") - - # Extract a slice of the matrix - matrix_slice <- as.matrix(h5_matrix[, start_col:end_col, drop=FALSE]) - - # Calculate ranks for the slice - ranked_slice <- singscore::rankGenes(matrix_slice) %>% `-` (1) - - # Convert the ranked slice to sparse format - sparse_ranked_slice <- as(ranked_slice, "CsparseMatrix") - - # Write the slice to the output HDF5 file - HDF5Array::writeHDF5Array( - sparse_ranked_slice, - filepath = output_filepath_temp, - name = paste0("rank_", start_col, "_to_", end_col), - as.sparse = TRUE, - H5type = "H5T_STD_I32LE" - ) - - # Store the slice name for later binding - slice_list[[length(slice_list) + 1]] <- paste0("rank_", start_col, "_to_", end_col) - } - - - slice_list |> map(~HDF5Array::HDF5Array(output_filepath_temp, name =.x)) |> do.call(cbind, args=_) - - } - - save_rank_per_cell = function(dataset_id_sce, cache_directory){ - - dir.create(cache_directory, recursive = TRUE, showWarnings = FALSE) - - .x = dataset_id_sce |> pull(sce) |> _[[1]] - .y = dataset_id_sce |> pull(file_id_cellNexus_single_cell) |> _[[1]] |> str_remove("\\.h5ad") - - # Check if the 'sce' has only one cell (column) - if(ncol(assay(.x)) == 1) { - - # Duplicate the assay to prevent saving errors due to single-column matrices - my_assay = cbind(assay(.x), assay(.x)) - # Rename the second column to distinguish it - colnames(my_assay)[2] = paste0("DUMMY", "___", colnames(my_assay)[2]) - - cd = colData(.x) - cd = cd |> rbind(cd) - rownames(cd)[2] = paste0("DUMMY", "___", rownames(cd)[2]) - - - - .x = SingleCellExperiment(assay = list( my_assay ) |> set_names(names(assays(.x))[1]), colData = cd) - } - - - # # TEMPORARY FOR SOME REASON THE MIN COUNTS IS NOT 0 FOR SOME SAMPLES - # .x = HPCell:::check_if_assay_minimum_count_is_zero_and_correct_TEMPORARY(.x, assays(.x) |> names() |> _[1], subset_up_to_number_of_cells = 100) - - print("start ranking") - - # CALCULATE rank - rank_assay = - .x |> - assay() |> - - # This because some datasets are still > 1M cells - process_matrix_in_slices( - paste(c(cache_directory, "/", .y, "_rank_matrix.HDF5Array"), collapse = ""), - paste(c(cache_directory, "/", .y, "_rank_matrix_temp.HDF5Array"), collapse = ""), - chunk_size = 1000 - ) - - print("creating SCE") - - .x = SingleCellExperiment(assay = list( rank = rank_assay), colData = colData(.x)) - - print("saving") - - .x |> save_experiment_data(glue("{cache_directory}/{.y}")) - - # Delete the temp file - file.remove(paste(c(cache_directory, "/", .y, "_rank_matrix_temp.HDF5Array"), collapse = "")) - - return(TRUE) # Indicate successful saving - } - - # Because they have an inconsistent failure. If I start the pipeline again they might work. Strange. - insistent_save_rank_per_cell <- purrr::insistently(save_rank_per_cell, rate = purrr::rate_delay(pause = 60, max_times = 3), quiet = FALSE) - - - save_anndata_sct = function(dataset_id_sce, cache_directory){ - - dir.create(cache_directory, showWarnings = FALSE, recursive = TRUE) - - if (is.null(dataset_id_sce)) return(NULL) - - .x = dataset_id_sce |> pull(sct) |> _[[1]] - - # Fix: check is.null BEFORE ncol() to avoid `argument is of length zero` - if (is.null(.x) || ncol(.x) == 0) return(NULL) - - .y = dataset_id_sce |> pull(file_id_cellNexus_single_cell) |> _[[1]] |> str_remove("\\.h5ad") - - .x |> assays() |> names() = "sct" - - # Wrap save with explicit error logging so the real cause is visible. Strange it shouldnt fail, which passed in debug mode - tryCatch( - .x |> save_experiment_data(glue("{cache_directory}/{.y}")), - error = function(e) { - message(glue::glue("[save_anndata_sct] FAILED for {.y}: {conditionMessage(e)}")) - stop(e) - } - ) - - return(TRUE) - } - - # Because they have an inconsistent failure. If I start the pipeline again they might work. Strange. - insistent_save_anndata_sct <- purrr::insistently(save_anndata_sct, rate = purrr::rate_delay(pause = 60, max_times = 3), quiet = FALSE) - - - cbind_sce_by_dataset_id = function(target_name_grouped_by_dataset_id, file_id_db_file, cell_id_dict, my_store){ - - my_dataset_id = unique(target_name_grouped_by_dataset_id$dataset_id) - my_file_id = unique(target_name_grouped_by_dataset_id$file_id_cellNexus_single_cell) - - file_id_db = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue("SELECT * FROM read_parquet('{file_id_db_file}')")) - ) |> - filter(dataset_id == my_dataset_id) |> - select(cell_id, sample_id, dataset_id, file_id_cellNexus_single_cell) - - file_id_db = - target_name_grouped_by_dataset_id |> - left_join(file_id_db, copy = TRUE) - - - dataset_cell_dict = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue("SELECT * FROM read_parquet('{cell_id_dict}')")) - ) |> - filter(file_id_cellNexus_single_cell == my_file_id) - - file_id_db = - file_id_db |> - left_join(dataset_cell_dict, by = c("file_id_cellNexus_single_cell", "cell_id" ), copy=T ) - - # Parallelise - cores = as.numeric(Sys.getenv("SLURM_CPUS_PER_TASK", unset = 1))-1 - # Respect R CMD CHECK core limit if set - if (nzchar(Sys.getenv("_R_CHECK_LIMIT_CORES_"))) { - cores <- min(cores, 2L) - } - # MulticoreParam need to have enough memory to proceed, otherwise reducer error - bp <- MulticoreParam(workers = cores , progressbar = TRUE) # Adjust the number of workers as needed - - # Begin processing the data pipeline with the initial dataset 'target_name_grouped_by_dataset_id' - sce_df = - file_id_db |> - nest(cells = c(cell_id, new_cell_id)) |> - # Read raw data for each 'target_name' and store it in a new column 'sce' - mutate( - sce = bplapply( - sce_target_name, - FUN = function(x) { - tar_read_raw(x, store = my_store) |> - select(.cell, donor_id, dataset_id, sample_id, cell_type) |> - mutate(sample_id = as.factor(sample_id)) # lighter - }, # Read the raw SingleCellExperiment object - BPPARAM = bp # Use the defined parallel backend - )) |> - # This should not be needed, but there are some data sets with zero cells - filter(!map_lgl(sce, is.null)) |> - mutate(sce = map2(sce, cells, ~ { - - cell_map <- setNames(.y$new_cell_id, .y$cell_id) - - .x |> filter(.cell %in% names(cell_map)) %>% - { - colnames(.) <- cell_map[colnames(.)] - . - } - - }, .progress = TRUE)) - - if(nrow(sce_df) == 0) { - warning("this chunk has no rows for somereason.") - return(NULL) - } - - sce_df |> - - # Group the data by 'dataset_id' and 'tar_group' for further summarization - mutate(sce = map(sce, ~ SingleCellExperiment(assay = assays(.x), colData = colData(.x)) )) |> - - # Combine all 'sce' objects within each group into a single 'sce' object - group_by(file_id_cellNexus_single_cell) |> - summarise( sce = list(do.call(cbind, args = sce) ), - # A step to check missing cells - cells = list(do.call(rbind, args = cells))) - } - - cbind_sct_by_dataset_id = function(target_name_grouped_by_dataset_id, file_id_db_file, cell_id_dict, my_store){ - - my_dataset_id = unique(target_name_grouped_by_dataset_id$dataset_id) - my_file_id = unique(target_name_grouped_by_dataset_id$file_id_cellNexus_single_cell) - - file_id_db = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue("SELECT * FROM read_parquet('{file_id_db_file}')")) - ) |> - filter(dataset_id == my_dataset_id) |> - select(cell_id, sample_id, dataset_id, file_id_cellNexus_single_cell) - - file_id_db = - target_name_grouped_by_dataset_id |> - left_join(file_id_db, copy = TRUE) - - - dataset_cell_dict = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue("SELECT * FROM read_parquet('{cell_id_dict}')")) - ) |> - filter(file_id_cellNexus_single_cell == my_file_id) - - file_id_db = - file_id_db |> - left_join(dataset_cell_dict, by = c("file_id_cellNexus_single_cell", "cell_id"), copy=T ) - - # Parallelise - cores = as.numeric(Sys.getenv("SLURM_CPUS_PER_TASK", unset = 1)) -1 - # Respect R CMD CHECK core limit if set - if (nzchar(Sys.getenv("_R_CHECK_LIMIT_CORES_"))) { - cores <- min(cores, 2L) - } - # MulticoreParam need to have enough memory to proceed, otherwise reducer error - bp <- MulticoreParam(workers = cores , progressbar = TRUE) # Adjust the number of workers as needed - - # Begin processing the data pipeline with the initial dataset 'target_name_grouped_by_dataset_id' - sct_df = file_id_db |> - nest(cells = c(cell_id, new_cell_id)) %>% - # Step 1: Read raw data for each 'target_name' and store it in a new column 'sce' - mutate( - sct = bplapply( - sct_target_name, - FUN = function(x) { - if (is.na(x)) { - return(NULL) # because cant get sample_id and dataset_id from NULL sct_matrix - } - - tar_read_raw(x, store = my_store) |> - select(.cell, donor_id, dataset_id, sample_id, cell_type) |> - mutate(sample_id = as.factor(sample_id)) - - }, # Read the raw SingleCellExperiment object - BPPARAM = bp # Use the defined parallel backend - )) |> - # This should not be needed, but there are some data sets with zero cells - filter(!map_lgl(sct, is.null)) |> - mutate(sct = map2(sct, cells, ~ { - - cell_map <- setNames(.y$new_cell_id, .y$cell_id) - - .x |> filter(.cell %in% names(cell_map)) %>% - { - colnames(.) <- cell_map[colnames(.)] - . - } - - }, .progress = TRUE)) - - if(nrow(sct_df) == 0) { - warning("this chunk has no rows for somereason.") - return(NULL) - } - - sct_df |> - mutate( - sct = map(sct, \(x) { - if (is.null(x)) return(NULL) - SingleCellExperiment(assays = assays(x), colData = colData(x)) - }) - ) |> - group_by(file_id_cellNexus_single_cell) |> - summarise( - sct = { - scts <- compact(sct) # drop NULLs inside each group - - list( - if (length(scts) == 0) { - NULL - } else { - - # A few big samples do not return all features because it reached R limit 2^31-1 in SCTransform - common_genes <- cellNexus:::check_gene_overlap(scts) - - # subset to intersection genes (and keep same order across objects) - scts2 <- map(scts, \(z) z[common_genes, , drop = FALSE]) - - do.call(SummarizedExperiment::cbind, scts2) - } - ) - }, - cells = list(do.call(rbind, cells)), - .groups = "drop" - ) - } - - # Because they have an inconsistent failure. If I start the pipeline again they might work. Strange. - insistent_cbind_sct_by_dataset_id <- purrr::insistently(cbind_sct_by_dataset_id, rate = purrr::rate_delay(pause = 60, max_times = 3), quiet = FALSE) - - - get_dataset_id = function(target_name, my_store){ - # Try reading the target safely (for some failing targets) - sce = tryCatch( - tar_read_raw(target_name, store = my_store), - error = function(e) return(NULL) - ) - - # Still need to catch target_name - if(sce |> is.null()) return(tibble(sample_id = NA_character_, - dataset_id= NA_character_, - target_name= !!target_name)) - - sce |> - - distinct(sample_id, dataset_id) |> mutate(target_name = !!target_name) - } - - create_chunks_for_reading_and_saving = function(dataset_id_sample_id, cell_metadata){ - - # Solve sample_id mismatches because some end with .h5ad suffix while others dont - dataset_id_sample_id |> - - left_join( - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue("SELECT * FROM read_parquet('{cell_metadata}')")) - ) |> - distinct(dataset_id, sample_id, sample_chunk, cell_chunk, file_id_cellNexus_single_cell) |> - as_tibble(), - copy=T - ) - } - - - cbind_sce_by_dataset_id_get_missing_cells = function(dataset_id_sce){ - - dataset_id_sce |> - mutate( - missing_cells = map2( - sce, - cells, - ~{ - cells_in_sce <- .x |> colnames() |> sort() - - cells_in_query <- .y$new_cell_id |> unique() |> sort() - - # Find differences - tibble(cell_id = setdiff(cells_in_query, cells_in_sce)) - } - ) - ) |> - select(file_id_cellNexus_single_cell, missing_cells) - - } - - - list( - - # The input DO NOT DELETE - tar_target(my_store, "/vast/scratch/users/shen.m/cellNexus/2024-07-01/process_updated_samples_transform_hpcell_target_store_v1", deployment = "main"), # MODIFY HERE: HPCell targets store to read SCEs from - tar_target(cache_directory, "/vast/scratch/users/shen.m/cellNexus/cellxgene_2024/0.3.0", deployment = "main"), # MODIFY HERE: output cache directory for saved anndata files - tar_target( - cell_metadata, - "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_1_mengyuan.parquet", # MODIFY HERE: final metadata parquet (should match the COPY TO output above) - packages = c( "arrow","dplyr","duckdb") - - ), - - tar_target( - cell_id_dict, - "/vast/projects//cellxgene_curated/metadata_cellxgene_mengyuan/file_id_cell_id_dict_v1_1_1_Jul_2024.parquet", # MODIFY HERE: cell_id dictionary parquet - packages = c( "arrow","dplyr","duckdb") - ), - - # pre-calculated counts - tar_target( - target_name, - tar_meta( - starts_with("sce_transformed_"), - store = my_store) |> - filter(type=="branch") |> - pull(name), - deployment = "main" - ), - tar_target( - dataset_id_sample_id, - get_dataset_id(target_name, my_store), - packages = "tidySingleCellExperiment", - pattern = map(target_name), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_5_minimal") - ) - ), - - # # pre-calculated sct - # tar_target( - # sct_target_name, - # tar_meta( - # starts_with("sct_matrix_"), - # store = my_store) |> - # filter(type=="branch") |> - # pull(name), - # deployment = "main" - # ), - # tar_target( - # sct_dataset_id_sample_id, - # get_dataset_id(sct_target_name, my_store), - # packages = "tidySingleCellExperiment", - # pattern = map(sct_target_name), - # resources = tar_resources( - # crew = tar_resources_crew(controller = "elastic_5_minimal") - # ) - # ), - - # join - tar_target( - dataset_id_sample_id_target_names, - dataset_id_sample_id |> #left_join(sct_dataset_id_sample_id, by = c("sample_id", "dataset_id"), copy=T) |> - dplyr::rename(sce_target_name = target_name), - # dplyr::rename(sce_target_name = target_name.x, - # sct_target_name = target_name.y), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_40") - ) - ), - - tar_target( - target_name_grouped_by_dataset_id, - create_chunks_for_reading_and_saving(dataset_id_sample_id_target_names, cell_metadata) |> - - # # FOR TESTING PURPOSE ONLY - # filter(file_id_cellNexus_single_cell %in% c("e8291e33fc98bc21728255f0a5669015___1.h5ad", - # "4ceb75a970ccdc9aaa9f9e91b931292f___1.h5ad")) |> - - group_by(dataset_id, sample_chunk, cell_chunk, file_id_cellNexus_single_cell) |> - tar_group(), - iteration = "group", - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_5_minimal") - ), - packages = c("arrow", "duckdb", "dplyr", "glue", "targets") - - ), - - tar_target( - dataset_id_sce, - cbind_sce_by_dataset_id(target_name_grouped_by_dataset_id, cell_metadata, cell_id_dict, my_store = my_store), - pattern = map(target_name_grouped_by_dataset_id), - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "digest", "scater", "HDF5Array", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly"), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_40") - ) - ), - - # tar_target( - # dataset_id_sct, - # cbind_sct_by_dataset_id(target_name_grouped_by_dataset_id, cell_metadata, cell_id_dict, my_store = my_store), - # pattern = map(target_name_grouped_by_dataset_id), - # packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "digest", "scater", "HDF5Array", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly"), - # resources = tar_resources( - # crew = tar_resources_crew(controller = "elastic_120") - # ) - # ), - - # This target was run for retrieving missing cells analysis only - # tar_target( - # missing_cells_tbl, - # cbind_sce_by_dataset_id_get_missing_cells(dataset_id_sce), - # pattern = map(dataset_id_sce), - # packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly", "purrr"), - # resources = tar_resources( - # crew = tar_resources_crew(controller = "elastic_20") - # ) - # ), - # - - tar_target( - save_anndata, - insistent_save_anndata(dataset_id_sce, paste0(cache_directory, "/counts")), - pattern = map(dataset_id_sce), - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "HPCell", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly"), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_20") - ) - ), - - tar_target( - saved_dataset_cpm, - insistent_save_anndata_cpm(dataset_id_sce, paste0(cache_directory, "/cpm")), - pattern = map(dataset_id_sce), - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "HPCell", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly"), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_20") - ) - ) - # - # tar_target( - # saved_dataset_rank, - # insistent_save_rank_per_cell(dataset_id_sce, paste0(cache_directory, "/rank")), - # pattern = map(dataset_id_sce), - # packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "HPCell", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly", "HDF5Array"), - # resources = tar_resources( - # crew = tar_resources_crew(controller = "elastic_20") - # ) - # ), - # - # tar_target( - # saved_sct, - # save_anndata_sct(dataset_id_sct, paste0(cache_directory, "/sct")), - # pattern = map(dataset_id_sct), - # packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "HPCell", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly", "HDF5Array"), - # resources = tar_resources( - # crew = tar_resources_crew(controller = "elastic_20") - # ) - # ) - ) - -}, script = paste0(store_file_cellNexus, "_target_script.R"), ask = FALSE) - -job::job({ - - tar_make( - script = paste0(store_file_cellNexus, "_target_script.R"), - store = store_file_cellNexus, - reporter = "summary" #, callr_function = NULL - ) - -}) - -missing_cells_tbl = tar_read(missing_cells_tbl, store = store_file_cellNexus) |> - unnest(missing_cells) - -#missing_cells_tbl |> write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cells_to_remove_in_metadata_Jul_2024.parquet") -missing_cells_tbl <- read_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cells_to_remove_in_metadata_Jul_2024.parquet") - -filtered_cell_metadata = cell_metadata |> anti_join(missing_cells_tbl, by = c("file_id_cellNexus_single_cell", - "new_cell_id" = "cell_id"), copy = T) - -filtered_cell_metadata |> - collect() |> - arrow::write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_1_filtered_missing_cells_mengyuan.parquet", - compression = "zstd") # MODIFY HERE: output parquet after filtering missing cells - diff --git a/dev/cellnexus-2024-scripts/step7_supp_dataset_cell_map.R b/dev/cellnexus-2024-scripts/step7_supp_dataset_cell_map.R deleted file mode 100644 index db2630c..0000000 --- a/dev/cellnexus-2024-scripts/step7_supp_dataset_cell_map.R +++ /dev/null @@ -1,105 +0,0 @@ -library(targets) -store = "/vast/scratch/users/shen.m/cellnexus_dataset_cell_map_Jul_2024_v1_2_1_target_store" -tar_script({ - library(dplyr) - library(magrittr) - library(tibble) - library(targets) - library(tarchetypes) - library(crew) - library(crew.cluster) - - tar_option_set( - memory = "transient", - garbage_collection = 100, - storage = "worker", - retrieval = "worker", - error = "continue", - cue = tar_cue(mode = "never"), - - workspace_on_error = TRUE, - controller = crew_controller_group( - list( - crew_controller_slurm( - name = "elastic", - workers = 300, - tasks_max = 20, - seconds_idle = 30, - crashes_error = 10, - options_cluster = crew_options_slurm( - memory_gigabytes_required = c(25, 35, 40, 80, 160), - cpus_per_task = c(2, 2, 5, 10, 20), - time_minutes = c(30, 30, 30, 60*4, 60*24), - verbose = T - ) - ) - ) - ), - trust_object_timestamps = TRUE - ) - - get_unique_file_ids <- function(cell_metadata){ - tbl(dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue::glue("SELECT * FROM read_parquet('{cell_metadata}')"))) |> - distinct(file_id_cellNexus_single_cell) |> pull() - } - - create_file_id_cell_id_dict <- function(cell_metadata, file_id) { - - - tbl(dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue::glue("SELECT * FROM read_parquet('{cell_metadata}')"))) |> - filter(file_id_cellNexus_single_cell == file_id) |> - dbplyr::window_order(cell_id) |> - mutate(cell_index = row_number()) |> - select(cell_id, file_id_cellNexus_single_cell, - new_cell_id = cell_index) |> - collect() - } - - list( - tar_target(cell_metadata , "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_0_mengyuan.parquet", - deployment = "main"), - tar_target( - unique_file_ids, - # TESTING PURPOSE ONLY - # c("3cef5b6aa0f5772485bb710f71e69456___1.h5ad", - # "cd2caa6de850f73af4ca78a2ea307dd4___1.h5ad") - get_unique_file_ids(cell_metadata) - # |> head(2) - , - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly", "HDF5Array") - # resources = tar_resources( - # crew = tar_resources_crew(controller = "elastic") - # ) - ), - tar_target( - file_id_cell_id_dict, - create_file_id_cell_id_dict(cell_metadata, unique_file_ids), - pattern = map(unique_file_ids), - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly", "HDF5Array") - # resources = tar_resources( - # crew = tar_resources_crew(controller = "elastic") - # ) - ) - ) - -}, script = paste0(store, "_target_script.R"), ask = FALSE) - - -job::job({ - - tar_make( - script = paste0(store, "_target_script.R"), - store = store, - reporter = "summary" - ) - -}) - -file_id_cell_id_dict = tar_read(file_id_cell_id_dict, store = store) -file_id_cell_id_dict |> arrow::write_parquet("/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/file_id_cell_id_dict_v1_1_1_Jul_2024.parquet", - compression = "zstd") -rm(file_id_cell_id_dict) -gc() - diff --git a/dev/cellnexus-2024-scripts/step7_unify_and_update_sce_metadata.R b/dev/cellnexus-2024-scripts/step7_unify_and_update_sce_metadata.R deleted file mode 100644 index 206920d..0000000 --- a/dev/cellnexus-2024-scripts/step7_unify_and_update_sce_metadata.R +++ /dev/null @@ -1,398 +0,0 @@ -# Description: -# This script clean up and generate the ultimate metadata to ship to cellNexus. It reads metadata and dataset-specific information, -# cleans and renames columns, and writes updated data back to disk. The process involves -# connecting to databases in memory, executing SQL queries, and handling data in both -# Parquet and HDF5 formats. - -library(duckdb) -library(dbplyr) -library(dplyr) -library(tidyr) -library(data.table) -library(HDF5Array) -library(SummarizedExperiment) -library(tidySingleCellExperiment) -library(stringr) -library(targets) -library(purrr) -library(arrow) - - -# Add low confidence ethnicity and imputed ethnicity labels to metadata. Both data are from Ning via email -lowConf_ethnicity_df <- zellkonverter::readH5AD("/vast/projects/cellxgene_curated/cellNexus/sce_relabel.h5ad", reader = "R", use_hdf5 = T) |> - colData() |> as_tibble() |> - mutate(low_confidence_ethnicity = ifelse(ethnicity_relabel == "LowConfidenceLabel", TRUE, FALSE) |> as.character()) |> - select(sample_id, ethnicity_flagging_score = score, low_confidence_ethnicity = low_confidence_ethnicity) - -imputed_ethnicity_df <- zellkonverter::readH5AD("/vast/projects/cellxgene_curated/cellNexus/adata_unlabelled_with_predictions.h5ad", reader = "R", use_hdf5 = T)|> - colData() |> as_tibble() |> - select(sample_id, imputed_ethnicity = ethnicity_predictions) |> - mutate(imputed_ethnicity = as.character(imputed_ethnicity)) - -# lowConf_ethnicity_df |> arrow::write_parquet("/vast/projects/cellxgene_curated/cellNexus/lowConf_ethnicity_df.parquet") -# imputed_ethnicity_df |> arrow::write_parquet("/vast/projects/cellxgene_curated/cellNexus/imputed_ethnicity_df.parquet") - -job::job({ - - duckdb_write_parquet <- function(.tbl_sql, path, con) { - - sql_tbl <- - .tbl_sql |> - sql_render() - - # zstd 15 compresses faster than brotli for binary/scientific datasets, whereas brotli reduce could save 100Mb - sql_call <- glue::glue("COPY ({sql_tbl}) TO '{path}' (FORMAT PARQUET, COMPRESSION 'brotli')") - - res <- dbExecute(con, sql_call) - - return(res) - } - - # Single DuckDB connection: do the heavy transforms in SQL (avoid read/write/read on 50M+ rows) - con <- DBI::dbConnect(duckdb::duckdb(), dbdir = ":memory:") - - raw_path <- "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_1_filtered_missing_cells_mengyuan.parquet" # MODIFY HERE: Metadata input parquet path - - DBI::dbExecute(con, glue::glue(" - CREATE VIEW cell_metadata_raw AS - SELECT * - FROM read_parquet({DBI::dbQuoteString(con, raw_path)}, union_by_name=true); - ")) - - raw_cols <- DBI::dbGetQuery(con, "SELECT * FROM cell_metadata_raw LIMIT 0") |> names() - - explicit_drop <- c() - # explicit_drop <- c( - # "cell_", - # "cell__1", - # "dataset_id_1", - # "dataset_id_1_1", - # "cell__2", - # "cell__3", - # "dataset_id_2", - # "dataset_id_3", - # "sample_id_1", - # "sample_id_2", - # "sample_placeholder", - # "cell_type_unified_ensemble_1", - # "cell_type_1", - # "dataset_id_2", - # "observation_joinid_1", - # "self_reported_ethnicity_1", - # "donor_id_1", - # "assay_1", - # "blueprint_first_labels_fine_1", - # "azimuth_predicted_celltype_l2_1", - # "monaco_first_labels_fine_1", - # "dataset_id_3", - # "atlas_id_1", - # "tissue_1", - # "is_primary_data_1", - # "cell_type_ontology_term_id_1", - # "azimuth", - # "blueprint", - # "monaco", - # "alive_1", - # "cell_id_1", - # "dataset_id_4", - # "X_umap1", - # "X_umap2", - # "observation_originalid", - # "subsets_Mito_sum", - # "subsets_Mito_detected", - # "file_id_cellNexus_single_cell_1", - # "ensemble_joinid", - # "cell_type_unified", - # "data_driven_ensemble" - # ) - - pattern_drop <- c( - grep("^scores", raw_cols, value = TRUE), - grep("coarse$", raw_cols, value = TRUE) - ) - - drop_cols <- intersect(unique(c(explicit_drop, pattern_drop)), raw_cols) - - int_cast_cols <- intersect( - unique( - c( - "feature_count", - "nFeature_expressed_in_sample", - "cell_count", - grep("metacell_", raw_cols, value = TRUE), - grep("_chunk", raw_cols, value = TRUE), - grep("subsets_", raw_cols, value = TRUE) - ) - ), - raw_cols - ) - - chr_cast_cols <- intersect(c("published_at", "revised_at"), raw_cols) - - sql_id <- function(x) as.character(DBI::dbQuoteIdentifier(con, x)) - - # Remove originals that we re-add under new names - base_keep <- setdiff( - raw_cols, - c( - drop_cols, - "alive", - "atlas_id", - "blueprint_first_labels_fine", - "monaco_first_labels_fine", - "azimuth_predicted_celltype_l2", - "cell_id", - "new_cell_id" - ) - ) - - select_exprs <- purrr::map_chr(base_keep, function(col) { - col_id <- sql_id(col) - if (col %in% int_cast_cols) { - glue::glue("CAST({col_id} AS INTEGER) AS {col_id}") - } else if (col %in% chr_cast_cols) { - glue::glue("CAST({col_id} AS VARCHAR) AS {col_id}") - } else { - col_id - } - }) - - # alive: NA -> FALSE - if ("alive" %in% raw_cols) { - select_exprs <- c(select_exprs, glue::glue("COALESCE({sql_id('alive')}, FALSE) AS {sql_id('alive')}")) - } else { - select_exprs <- c(select_exprs, glue::glue("FALSE AS {sql_id('alive')}")) - } - - # Rename annotation columns - select_exprs <- c( - select_exprs, - if ("blueprint_first_labels_fine" %in% raw_cols) { - glue::glue("{sql_id('blueprint_first_labels_fine')} AS {sql_id('cell_annotation_blueprint_singler')}") - } else { - glue::glue("NULL::VARCHAR AS {sql_id('cell_annotation_blueprint_singler')}") - }, - if ("monaco_first_labels_fine" %in% raw_cols) { - glue::glue("{sql_id('monaco_first_labels_fine')} AS {sql_id('cell_annotation_monaco_singler')}") - } else { - glue::glue("NULL::VARCHAR AS {sql_id('cell_annotation_monaco_singler')}") - }, - if ("azimuth_predicted_celltype_l2" %in% raw_cols) { - glue::glue("{sql_id('azimuth_predicted_celltype_l2')} AS {sql_id('cell_annotation_azimuth_l2')}") - } else { - glue::glue("NULL::VARCHAR AS {sql_id('cell_annotation_azimuth_l2')}") - } - ) - - # new_cell_id -> cell_id as first column (drop original cell_id entirely) - cell_id_expr <- if ("new_cell_id" %in% raw_cols) { - glue::glue("{sql_id('new_cell_id')} AS {sql_id('cell_id')}") - } else { - glue::glue("NULL::VARCHAR AS {sql_id('cell_id')}") - } - select_exprs <- c(cell_id_expr, select_exprs) - - select_sql <- paste(select_exprs, collapse = ",\n ") - - DBI::dbExecute(con, glue::glue(" - CREATE OR REPLACE VIEW cell_metadata AS - SELECT - {DBI::SQL(select_sql)} - FROM cell_metadata_raw - WHERE dataset_id NOT IN ('99950e99-2758-41d2-b2c9-643edcdf6d82', '9fcb0b73-c734-40a5-be9c-ace7eea401c9'); - ")) - - DBI::dbExecute(con, " - CREATE OR REPLACE VIEW sample_celltype_count AS - SELECT - sample_id, - cell_type_unified_ensemble, - CAST(COUNT(*) AS INTEGER) AS \".aggregated_cells\" - FROM cell_metadata - WHERE empty_droplet = FALSE - AND alive = TRUE - AND \"scDblFinder.class\" != 'doublet' - GROUP BY sample_id, cell_type_unified_ensemble; - ") - - gc() - - dbExecute(con, " - CREATE VIEW lowConf_ethnicity_df AS - SELECT - * - FROM read_parquet('/vast/projects/cellxgene_curated/cellNexus/lowConf_ethnicity_df.parquet') -") - - dbExecute(con, " - CREATE VIEW imputed_ethnicity_df AS - SELECT - * - FROM read_parquet('/vast/projects/cellxgene_curated/cellNexus/imputed_ethnicity_df.parquet') -") - - # Perform left join and save to parquet - # MODIFY HERE: output metadata parquet path and atlas_id - copy_query <- " - COPY ( - SELECT - cell_metadata.*, - lowConf_ethnicity_df.ethnicity_flagging_score, - lowConf_ethnicity_df.low_confidence_ethnicity, - sample_celltype_count.\".aggregated_cells\", - COALESCE(imputed_ethnicity_df.imputed_ethnicity, cell_metadata.self_reported_ethnicity) AS imputed_ethnicity, -- Use imputed_ethnicity if present - 'cellxgene_2024/0.2.1' AS atlas_id - - FROM cell_metadata - - LEFT JOIN lowConf_ethnicity_df - ON cell_metadata.sample_id = lowConf_ethnicity_df.sample_id - - LEFT JOIN imputed_ethnicity_df - ON cell_metadata.sample_id = imputed_ethnicity_df.sample_id - - LEFT JOIN sample_celltype_count - ON cell_metadata.sample_id = sample_celltype_count.sample_id AND cell_metadata.cell_type_unified_ensemble = sample_celltype_count.cell_type_unified_ensemble - - - - - ) TO '/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/metadata.2.2.1.parquet' - (FORMAT PARQUET, COMPRESSION 'zstd'); - " - - # Execute the final query to write the result to a Parquet file - dbExecute(con, copy_query) - - # Disconnect from the database - dbDisconnect(con, shutdown = TRUE) - - print("Done.") - - -}) - -x = tbl(dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql("SELECT * FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/metadata.2.2.1.parquet')") ) # MODIFY HERE: input metadata parquet path - -# Split cell_metadata to cellnexus_metadata, original census_metadata, and metacell_metadata (host Rshiny on smaller file) -# ---- Split: read metadata.x.y.z.parquet once, write smaller derivative Parquets ---- -# This avoids loading the full metadata into R memory and avoids duplicate-column issues in Parquet writes. -job::job({ - - con <- DBI::dbConnect(duckdb::duckdb(), dbdir = ":memory:") - on.exit(DBI::dbDisconnect(con, shutdown = TRUE), add = TRUE) - - input_metadata <- "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/metadata.2.2.1.parquet" # MODIFY HERE: input metadata parquet path - out_dir <- "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan" - - DBI::dbExecute( - con, - glue::glue( - " - CREATE OR REPLACE VIEW metadata AS - SELECT * - FROM read_parquet({DBI::dbQuoteString(con, input_metadata)}, union_by_name=true); - " - ) - ) - - cols <- DBI::dbGetQuery(con, "SELECT * FROM metadata LIMIT 0") |> names() - sql_id <- function(x) as.character(DBI::dbQuoteIdentifier(con, x)) - - # Strip sample annotation from cellnexus annotation doesn't save too much (less than 3Mb), thus keep in one. - remove_cols <- c( - "cell_type", "cell_type_ontology_term_id", "data_driven_ensemble", "ensemble_joinid", - "observation_originalid", "assay", "assay_ontology_term_id", "development_stage", "development_stage_ontology_term_id", - "disease", "disease_ontology_term_id", "donor_id", "is_primary_data", "organism", "organism_ontology_term_id", - "self_reported_ethnicity", "self_reported_ethnicity_ontology_term_id", - "sex", "sex_ontology_term_id", "tissue", "tissue_ontology_term_id", "citation", - "collection_id", "dataset_version_id", "default_embedding", "published_at", "raw_data_location", - "revised_at", "primary_cell_count", "schema_version", "tissue_type", "title", - "tombstone", "x_approximate_distribution", "explorer_url", "cell_count", "feature_count", - "filesize", "filetype", "mean_genes_per_cell", "suspension_type", "url" - ) - - # CellNexus metadata (smaller file for Shiny): drop heavy / internal columns by name patterns - drop_cellnexus <- unique(c( - intersect(remove_cols, cols), - cols[grepl("metacell", cols)] - )) - keep_cellnexus <- setdiff(cols, drop_cellnexus) - select_cellnexus <- paste(sql_id(keep_cellnexus), collapse = ", ") - - # MODIFY HERE: output cellnexus metadata parquet path - DBI::dbExecute( - con, - glue::glue( - " - COPY ( - SELECT {DBI::SQL(select_cellnexus)} - FROM metadata - ) - TO {DBI::dbQuoteString(con, file.path(out_dir, 'cellnexus_metadata.2.2.1.parquet'))} - (FORMAT PARQUET, COMPRESSION 'brotli'); - " - ) - ) - - # Original census-like metadata subset (stable columns) - census_cols <- intersect( - c( - "observation_joinid", "dataset_id", "sample_id", "cell_type", - "cell_type_ontology_term_id", "assay", "assay_ontology_term_id", "development_stage", "development_stage_ontology_term_id", - "disease", "disease_ontology_term_id", "donor_id", "is_primary_data", "organism", "organism_ontology_term_id", - "self_reported_ethnicity", "self_reported_ethnicity_ontology_term_id", - "sex", "sex_ontology_term_id", "tissue", "tissue_ontology_term_id", - "data_driven_ensemble", "ensemble_joinid", "observation_originalid", "citation", - "collection_id", "dataset_version_id", "default_embedding", "published_at", "raw_data_location", - "revised_at", "primary_cell_count", "schema_version", "tissue_type", "title", - "tombstone", "x_approximate_distribution", "explorer_url", "cell_count", "feature_count", - "filesize", "filetype", "mean_genes_per_cell", "suspension_type", "url" - ), - cols - ) - select_census <- paste(sql_id(census_cols), collapse = ", ") - - # MODIFY HERE: output census metadata parquet path - DBI::dbExecute( - con, - glue::glue( - " - COPY ( - SELECT {DBI::SQL(select_census)} - FROM metadata - ) - TO {DBI::dbQuoteString(con, file.path(out_dir, 'census_cell_metadata.2.2.1.parquet'))} - (FORMAT PARQUET, COMPRESSION 'brotli'); - " - ) - ) - - # # Metacell metadata subset - # metacell_cols <- unique(c("cell_id", "sample_id", "dataset_id", cols[grepl("metacell", cols)])) - # metacell_cols <- intersect(metacell_cols, cols) - # select_metacell <- paste(sql_id(metacell_cols), collapse = ", ") - # - # # MODIFY HERE: output metacell metadata parquet path - # DBI::dbExecute( - # con, - # glue::glue( - # " - # COPY ( - # SELECT {DBI::SQL(select_metacell)} - # FROM metadata - # ) - # TO {DBI::dbQuoteString(con, file.path(out_dir, 'metacell_metadata.2.2.1.parquet'))} - # (FORMAT PARQUET, COMPRESSION 'brotli'); - # " - # ) - # ) - - print("Done.") -}) - - -# (Optional) Check whether cellnexus_metadata parquet can be optimised further -# source("~/git_control/cellNexus/dev/data_optimisation_script.R") - diff --git a/dev/cellnexus-2024-scripts/step8_prepare_pseudobulk_local_cache.R b/dev/cellnexus-2024-scripts/step8_prepare_pseudobulk_local_cache.R deleted file mode 100644 index f3c3871..0000000 --- a/dev/cellnexus-2024-scripts/step8_prepare_pseudobulk_local_cache.R +++ /dev/null @@ -1,306 +0,0 @@ -library(targets) -library(tidyverse) -library(cellNexus) -store_file_cellNexus = "/vast/scratch/users/shen.m/targets_prepare_database_split_datasets_chunked_1_4_1_pseudobulk" # MODIFY HERE: targets store directory for this pipeline -my_store = "/vast/scratch/users/shen.m/cellNexus/2024-07-01/process_samples_hpcell_target_store" # MODIFY HERE: HPCell targets store to read pseudobulk SCEs from (used throughout) - -tar_script({ - library(dplyr) - library(magrittr) - library(tibble) - library(targets) - library(tarchetypes) - library(crew) - library(crew.cluster) - - # Helper (optional) to avoid repetition - new_elastic <- function(name, mem_gb, time_min, workers, crashes_max, cpus_per_task = 2, backup = NULL) { - crew_controller_slurm( - name = name, - workers = workers, - crashes_max = crashes_max, - seconds_idle = 30, - options_cluster = crew_options_slurm( - memory_gigabytes_required = mem_gb, - cpus_per_task = cpus_per_task, - time_minutes = time_min - ), - backup = backup - ) - } - - # Small → large, with fallbacks to the next size up - elastic_160 <- new_elastic("elastic_160", 160, 60 * 24, workers = 8, crashes_max = 2) - elastic_120 <- new_elastic("elastic_120", 120, 60 * 4, workers = 16, crashes_max = 1, cpus_per_task = 8, backup = elastic_160) - elastic_80 <- new_elastic("elastic_80", 80, 60 * 4, workers = 24, crashes_max = 1, cpus_per_task = 8, backup = elastic_120) - elastic_40 <- new_elastic("elastic_40", 40, 60 * 4, workers = 32, crashes_max = 1, cpus_per_task = 8, backup = elastic_80) - elastic_20 <- new_elastic("elastic_20", 20, 60 * 4, workers = 48, crashes_max = 1, cpus_per_task = 8, backup = elastic_40) - elastic_10 <- new_elastic("elastic_10", 10, 60 * 4, workers = 150, crashes_max = 6, cpus_per_task = 8, backup = elastic_20) - - elastic_5_minimal <- new_elastic("elastic_5_minimal", 5, 60 * 4, workers = 300, crashes_max = 6, cpus_per_task = 2, backup = elastic_10) - - - # Group for targets (small → large) - controllers <- crew_controller_group( - elastic_10, elastic_20, elastic_40, elastic_80, elastic_120, elastic_160, elastic_5_minimal - ) - tar_option_set( - memory = "transient", - garbage_collection = 100, - storage = "worker", - retrieval = "worker", - error = "continue", - cue = tar_cue(mode = "never"), - - workspace_on_error = TRUE, - controller = controllers, - trust_object_timestamps = TRUE, - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_5_minimal") - ) - ) - - - get_dataset_id = function(target_name, my_store){ - sce = tar_read_raw(target_name, store = my_store) - - if(sce |> is.null()) return(tibble(sample_id = character(), dataset_id= character(), - target_name= target_name)) - - sce |> - - distinct(sample_id, dataset_id) |> mutate(target_name = !!target_name) - } - - create_chunks_for_reading_and_saving = function(dataset_id_sample_id, cell_metadata){ - - # Solve sample_id mismatches because some end with .h5ad suffix while others dont - dataset_id_sample_id |> - - left_join( - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue("SELECT * FROM read_parquet('{cell_metadata}')")) - ) |> - distinct(sample_id, sample_pseudobulk_chunk, cell_chunk, - cell_type_unified_ensemble, - file_id_cellNexus_pseudobulk) |> - as_tibble(), - copy=T - ) - } - - - cbind_sce_by_dataset_id = function(target_name_grouped_by_dataset_id, - file_id_db_file, my_store){ - - #my_dataset_id = unique(target_name_grouped_by_dataset_id$dataset_id) - my_cell_type = unique(target_name_grouped_by_dataset_id$cell_type_unified_ensemble) - - file_id_db = - tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue("SELECT * FROM read_parquet('{file_id_db_file}')")) - ) |> - dplyr::filter(cell_type_unified_ensemble %in% my_cell_type) |> - select(sample_id, dataset_id, cell_type_unified_ensemble, - file_id_cellNexus_pseudobulk) - - - file_id_db = - target_name_grouped_by_dataset_id |> - left_join(file_id_db, copy = TRUE) - - - # Parallelise - cores = as.numeric(Sys.getenv("SLURM_CPUS_PER_TASK", unset = 1)) -1 - # Respect R CMD CHECK core limit if set - if (nzchar(Sys.getenv("_R_CHECK_LIMIT_CORES_"))) { - cores <- min(cores, 2L) - } - bp <- MulticoreParam(workers = cores, progressbar = TRUE) - - # Begin processing the data pipeline with the initial dataset 'target_name_grouped_by_dataset_id' - sce_df = - file_id_db |> - mutate(cell_id = paste(sample_id, cell_type_unified_ensemble, sep = "___")) |> - nest(cells = cell_id) |> - # Step 1: Read raw data for each 'target_name' and store it in a new column 'sce' - mutate( - sce = bplapply( - target_name, - FUN = function(x) tar_read_raw(x, store = my_store), - BPPARAM = bp - ) - ) |> - - # This should not be needed, but there are some data sets with zero cells - filter(!map_lgl(sce, is.null)) |> - - mutate(sce = map2(sce, cells, ~ .x |> - filter(.cell %in% .y$cell_id), - - .progress = TRUE)) - - - - if(nrow(sce_df) == 0) { - warning("this chunk has no rows for somereason.") - return(NULL) - } - - sce_df = sce_df |> - - # THIS SHOULD HAVE BEEN DONE IN THE TRANFORM HPCell - mutate(sce = map(sce, ~ SingleCellExperiment(assay = assays(.x), colData = colData(.x)) )) - - - # Extra Step 1: Harmonize colData columns - Avoid column name mismatch, force cbind - all_col_names <- sce_df$sce %>% - map(~colnames(colData(.x))) %>% - unlist() %>% - unique() - - # Extra Step 2: Standardize colData to have the same columns in each SCE - sce_df$sce <- map(sce_df$sce, function(sce) { - current_cols <- colnames(colData(sce)) - missing_cols <- setdiff(all_col_names, current_cols) - - if (length(missing_cols) > 0) { - - # Fill missing colData columns with NA - for (col in missing_cols) { - # Handle sce with empty cells - if (ncol(sce) == 0) colData(sce)[, col] <- character(0) - else if (ncol(sce) > 0) colData(sce)[, col] <- NA - } - } - - # Ensure the order of columns matches - colData(sce) <- colData(sce)[, all_col_names] - return(sce) - }) - - sce_df |> - - # Step 5: Combine all 'sce' objects within each group into a single 'sce' object - group_by(file_id_cellNexus_pseudobulk) |> - summarise( sce = list(do.call(cbind, args = sce) ) ) - - } - - - - save_anndata = function(dataset_id_sce, cache_directory){ - - dir.create(cache_directory, showWarnings = FALSE, recursive = TRUE) - - .x = dataset_id_sce |> pull(sce) |> _[[1]] - .y = dataset_id_sce |> pull(file_id_cellNexus_pseudobulk) |> _[[1]] |> str_remove("\\.h5ad") - - .x |> assays() |> names() = "counts" - - # Drop list-type columns in colData - cd <- colData(.x) - is_list_col <- vapply(cd, is.list, logical(1)) - colData(.x) <- cd[, !is_list_col, drop = FALSE] - - # Check if there is a memory issue - assays(.x) <- assays(.x) |> map(DelayedArray::realize) - - # Save the experiment data to the specified counts cache directory - .x |> save_experiment_data(glue("{cache_directory}/{.y}")) - - return(TRUE) # Indicate successful saving - - } - - # Because they have an inconsistent failure. If I start the pipeline again they might work. Strange. - insistent_save_anndata <- purrr::insistently(save_anndata, rate = purrr::rate_delay(pause = 60, max_times = 3), quiet = FALSE) - - list( - - # The input DO NOT DELETE - tar_target(my_store, "/vast/scratch/users/shen.m/cellNexus/2024-07-01/process_samples_hpcell_target_store", deployment = "main"), # MODIFY HERE: HPCell targets store (must match my_store above) - tar_target(cache_directory, "/vast/scratch/users/shen.m/cellNexus/cellxgene_2024/0.2.1/pseudobulk", deployment = "main"), # MODIFY HERE: output cache directory for saved pseudobulk anndata files - tar_target( - cell_metadata, - "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_1_mengyuan.parquet", # MODIFY HERE: cell metadata parquet (output of step6/step7) - packages = c( "arrow","dplyr","duckdb") - - ), - tar_target( - target_name, - tar_meta( - starts_with("pseudobulk_se_iterated_"), - store = my_store) |> - filter(type=="branch") |> - pull(name), - deployment = "main" - ), - tar_target( - dataset_id_sample_id, - get_dataset_id(target_name, my_store), - packages = "tidySingleCellExperiment", - pattern = map(target_name), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_10") - ) - ), - - tar_target( - target_name_grouped_by_dataset_id, - create_chunks_for_reading_and_saving(dataset_id_sample_id, cell_metadata) |> - - # # FOR TESTING PURPOSE ONLY - # filter(file_id_cellNexus_pseudobulk %in% c("9722bedfd71d069fe3665b4ae03fbeb9___2.h5ad", - # "2996bb4263f9fb301d8460f4f0450848___2.h5ad")) |> - - group_by(dataset_id, - sample_pseudobulk_chunk, - # When using strategy file_id = dataset_id, dont group by cell_chunk as it will result in returning more than one SCEs for the same dataset_id - #cell_chunk, - file_id_cellNexus_pseudobulk) |> - tar_group(), - iteration = "group", - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_5_minimal") - ), - packages = c("arrow", "duckdb", "dplyr", "glue", "targets") - - ), - - tar_target( - dataset_id_sce, - cbind_sce_by_dataset_id(target_name_grouped_by_dataset_id, cell_metadata, my_store = my_store), - pattern = map(target_name_grouped_by_dataset_id), - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly"), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_20") - ) - ), - tar_target( - get_pseudobulk, - insistent_save_anndata(dataset_id_sce, paste0(cache_directory, "/counts")), - pattern = map(dataset_id_sce), - packages = c("tidySingleCellExperiment", "SingleCellExperiment", "tidyverse", "glue", "HPCell", "digest", "scater", "arrow", "dplyr", "duckdb", "BiocParallel", "parallelly"), - resources = tar_resources( - crew = tar_resources_crew(controller = "elastic_20") - ) - ) - ) - - - -}, script = paste0(store_file_cellNexus, "_target_script.R"), ask = FALSE) - -job::job({ - - tar_make( - script = paste0(store_file_cellNexus, "_target_script.R"), - store = store_file_cellNexus, - reporter = "summary" #, callr_function = NULL - ) - -}) - diff --git a/dev/cellnexus-2024-scripts/step8_unify_and_update_sce_metadata.R b/dev/cellnexus-2024-scripts/step8_unify_and_update_sce_metadata.R deleted file mode 100644 index 206920d..0000000 --- a/dev/cellnexus-2024-scripts/step8_unify_and_update_sce_metadata.R +++ /dev/null @@ -1,398 +0,0 @@ -# Description: -# This script clean up and generate the ultimate metadata to ship to cellNexus. It reads metadata and dataset-specific information, -# cleans and renames columns, and writes updated data back to disk. The process involves -# connecting to databases in memory, executing SQL queries, and handling data in both -# Parquet and HDF5 formats. - -library(duckdb) -library(dbplyr) -library(dplyr) -library(tidyr) -library(data.table) -library(HDF5Array) -library(SummarizedExperiment) -library(tidySingleCellExperiment) -library(stringr) -library(targets) -library(purrr) -library(arrow) - - -# Add low confidence ethnicity and imputed ethnicity labels to metadata. Both data are from Ning via email -lowConf_ethnicity_df <- zellkonverter::readH5AD("/vast/projects/cellxgene_curated/cellNexus/sce_relabel.h5ad", reader = "R", use_hdf5 = T) |> - colData() |> as_tibble() |> - mutate(low_confidence_ethnicity = ifelse(ethnicity_relabel == "LowConfidenceLabel", TRUE, FALSE) |> as.character()) |> - select(sample_id, ethnicity_flagging_score = score, low_confidence_ethnicity = low_confidence_ethnicity) - -imputed_ethnicity_df <- zellkonverter::readH5AD("/vast/projects/cellxgene_curated/cellNexus/adata_unlabelled_with_predictions.h5ad", reader = "R", use_hdf5 = T)|> - colData() |> as_tibble() |> - select(sample_id, imputed_ethnicity = ethnicity_predictions) |> - mutate(imputed_ethnicity = as.character(imputed_ethnicity)) - -# lowConf_ethnicity_df |> arrow::write_parquet("/vast/projects/cellxgene_curated/cellNexus/lowConf_ethnicity_df.parquet") -# imputed_ethnicity_df |> arrow::write_parquet("/vast/projects/cellxgene_curated/cellNexus/imputed_ethnicity_df.parquet") - -job::job({ - - duckdb_write_parquet <- function(.tbl_sql, path, con) { - - sql_tbl <- - .tbl_sql |> - sql_render() - - # zstd 15 compresses faster than brotli for binary/scientific datasets, whereas brotli reduce could save 100Mb - sql_call <- glue::glue("COPY ({sql_tbl}) TO '{path}' (FORMAT PARQUET, COMPRESSION 'brotli')") - - res <- dbExecute(con, sql_call) - - return(res) - } - - # Single DuckDB connection: do the heavy transforms in SQL (avoid read/write/read on 50M+ rows) - con <- DBI::dbConnect(duckdb::duckdb(), dbdir = ":memory:") - - raw_path <- "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/cell_metadata_cell_type_consensus_v1_5_1_filtered_missing_cells_mengyuan.parquet" # MODIFY HERE: Metadata input parquet path - - DBI::dbExecute(con, glue::glue(" - CREATE VIEW cell_metadata_raw AS - SELECT * - FROM read_parquet({DBI::dbQuoteString(con, raw_path)}, union_by_name=true); - ")) - - raw_cols <- DBI::dbGetQuery(con, "SELECT * FROM cell_metadata_raw LIMIT 0") |> names() - - explicit_drop <- c() - # explicit_drop <- c( - # "cell_", - # "cell__1", - # "dataset_id_1", - # "dataset_id_1_1", - # "cell__2", - # "cell__3", - # "dataset_id_2", - # "dataset_id_3", - # "sample_id_1", - # "sample_id_2", - # "sample_placeholder", - # "cell_type_unified_ensemble_1", - # "cell_type_1", - # "dataset_id_2", - # "observation_joinid_1", - # "self_reported_ethnicity_1", - # "donor_id_1", - # "assay_1", - # "blueprint_first_labels_fine_1", - # "azimuth_predicted_celltype_l2_1", - # "monaco_first_labels_fine_1", - # "dataset_id_3", - # "atlas_id_1", - # "tissue_1", - # "is_primary_data_1", - # "cell_type_ontology_term_id_1", - # "azimuth", - # "blueprint", - # "monaco", - # "alive_1", - # "cell_id_1", - # "dataset_id_4", - # "X_umap1", - # "X_umap2", - # "observation_originalid", - # "subsets_Mito_sum", - # "subsets_Mito_detected", - # "file_id_cellNexus_single_cell_1", - # "ensemble_joinid", - # "cell_type_unified", - # "data_driven_ensemble" - # ) - - pattern_drop <- c( - grep("^scores", raw_cols, value = TRUE), - grep("coarse$", raw_cols, value = TRUE) - ) - - drop_cols <- intersect(unique(c(explicit_drop, pattern_drop)), raw_cols) - - int_cast_cols <- intersect( - unique( - c( - "feature_count", - "nFeature_expressed_in_sample", - "cell_count", - grep("metacell_", raw_cols, value = TRUE), - grep("_chunk", raw_cols, value = TRUE), - grep("subsets_", raw_cols, value = TRUE) - ) - ), - raw_cols - ) - - chr_cast_cols <- intersect(c("published_at", "revised_at"), raw_cols) - - sql_id <- function(x) as.character(DBI::dbQuoteIdentifier(con, x)) - - # Remove originals that we re-add under new names - base_keep <- setdiff( - raw_cols, - c( - drop_cols, - "alive", - "atlas_id", - "blueprint_first_labels_fine", - "monaco_first_labels_fine", - "azimuth_predicted_celltype_l2", - "cell_id", - "new_cell_id" - ) - ) - - select_exprs <- purrr::map_chr(base_keep, function(col) { - col_id <- sql_id(col) - if (col %in% int_cast_cols) { - glue::glue("CAST({col_id} AS INTEGER) AS {col_id}") - } else if (col %in% chr_cast_cols) { - glue::glue("CAST({col_id} AS VARCHAR) AS {col_id}") - } else { - col_id - } - }) - - # alive: NA -> FALSE - if ("alive" %in% raw_cols) { - select_exprs <- c(select_exprs, glue::glue("COALESCE({sql_id('alive')}, FALSE) AS {sql_id('alive')}")) - } else { - select_exprs <- c(select_exprs, glue::glue("FALSE AS {sql_id('alive')}")) - } - - # Rename annotation columns - select_exprs <- c( - select_exprs, - if ("blueprint_first_labels_fine" %in% raw_cols) { - glue::glue("{sql_id('blueprint_first_labels_fine')} AS {sql_id('cell_annotation_blueprint_singler')}") - } else { - glue::glue("NULL::VARCHAR AS {sql_id('cell_annotation_blueprint_singler')}") - }, - if ("monaco_first_labels_fine" %in% raw_cols) { - glue::glue("{sql_id('monaco_first_labels_fine')} AS {sql_id('cell_annotation_monaco_singler')}") - } else { - glue::glue("NULL::VARCHAR AS {sql_id('cell_annotation_monaco_singler')}") - }, - if ("azimuth_predicted_celltype_l2" %in% raw_cols) { - glue::glue("{sql_id('azimuth_predicted_celltype_l2')} AS {sql_id('cell_annotation_azimuth_l2')}") - } else { - glue::glue("NULL::VARCHAR AS {sql_id('cell_annotation_azimuth_l2')}") - } - ) - - # new_cell_id -> cell_id as first column (drop original cell_id entirely) - cell_id_expr <- if ("new_cell_id" %in% raw_cols) { - glue::glue("{sql_id('new_cell_id')} AS {sql_id('cell_id')}") - } else { - glue::glue("NULL::VARCHAR AS {sql_id('cell_id')}") - } - select_exprs <- c(cell_id_expr, select_exprs) - - select_sql <- paste(select_exprs, collapse = ",\n ") - - DBI::dbExecute(con, glue::glue(" - CREATE OR REPLACE VIEW cell_metadata AS - SELECT - {DBI::SQL(select_sql)} - FROM cell_metadata_raw - WHERE dataset_id NOT IN ('99950e99-2758-41d2-b2c9-643edcdf6d82', '9fcb0b73-c734-40a5-be9c-ace7eea401c9'); - ")) - - DBI::dbExecute(con, " - CREATE OR REPLACE VIEW sample_celltype_count AS - SELECT - sample_id, - cell_type_unified_ensemble, - CAST(COUNT(*) AS INTEGER) AS \".aggregated_cells\" - FROM cell_metadata - WHERE empty_droplet = FALSE - AND alive = TRUE - AND \"scDblFinder.class\" != 'doublet' - GROUP BY sample_id, cell_type_unified_ensemble; - ") - - gc() - - dbExecute(con, " - CREATE VIEW lowConf_ethnicity_df AS - SELECT - * - FROM read_parquet('/vast/projects/cellxgene_curated/cellNexus/lowConf_ethnicity_df.parquet') -") - - dbExecute(con, " - CREATE VIEW imputed_ethnicity_df AS - SELECT - * - FROM read_parquet('/vast/projects/cellxgene_curated/cellNexus/imputed_ethnicity_df.parquet') -") - - # Perform left join and save to parquet - # MODIFY HERE: output metadata parquet path and atlas_id - copy_query <- " - COPY ( - SELECT - cell_metadata.*, - lowConf_ethnicity_df.ethnicity_flagging_score, - lowConf_ethnicity_df.low_confidence_ethnicity, - sample_celltype_count.\".aggregated_cells\", - COALESCE(imputed_ethnicity_df.imputed_ethnicity, cell_metadata.self_reported_ethnicity) AS imputed_ethnicity, -- Use imputed_ethnicity if present - 'cellxgene_2024/0.2.1' AS atlas_id - - FROM cell_metadata - - LEFT JOIN lowConf_ethnicity_df - ON cell_metadata.sample_id = lowConf_ethnicity_df.sample_id - - LEFT JOIN imputed_ethnicity_df - ON cell_metadata.sample_id = imputed_ethnicity_df.sample_id - - LEFT JOIN sample_celltype_count - ON cell_metadata.sample_id = sample_celltype_count.sample_id AND cell_metadata.cell_type_unified_ensemble = sample_celltype_count.cell_type_unified_ensemble - - - - - ) TO '/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/metadata.2.2.1.parquet' - (FORMAT PARQUET, COMPRESSION 'zstd'); - " - - # Execute the final query to write the result to a Parquet file - dbExecute(con, copy_query) - - # Disconnect from the database - dbDisconnect(con, shutdown = TRUE) - - print("Done.") - - -}) - -x = tbl(dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql("SELECT * FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/metadata.2.2.1.parquet')") ) # MODIFY HERE: input metadata parquet path - -# Split cell_metadata to cellnexus_metadata, original census_metadata, and metacell_metadata (host Rshiny on smaller file) -# ---- Split: read metadata.x.y.z.parquet once, write smaller derivative Parquets ---- -# This avoids loading the full metadata into R memory and avoids duplicate-column issues in Parquet writes. -job::job({ - - con <- DBI::dbConnect(duckdb::duckdb(), dbdir = ":memory:") - on.exit(DBI::dbDisconnect(con, shutdown = TRUE), add = TRUE) - - input_metadata <- "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/metadata.2.2.1.parquet" # MODIFY HERE: input metadata parquet path - out_dir <- "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan" - - DBI::dbExecute( - con, - glue::glue( - " - CREATE OR REPLACE VIEW metadata AS - SELECT * - FROM read_parquet({DBI::dbQuoteString(con, input_metadata)}, union_by_name=true); - " - ) - ) - - cols <- DBI::dbGetQuery(con, "SELECT * FROM metadata LIMIT 0") |> names() - sql_id <- function(x) as.character(DBI::dbQuoteIdentifier(con, x)) - - # Strip sample annotation from cellnexus annotation doesn't save too much (less than 3Mb), thus keep in one. - remove_cols <- c( - "cell_type", "cell_type_ontology_term_id", "data_driven_ensemble", "ensemble_joinid", - "observation_originalid", "assay", "assay_ontology_term_id", "development_stage", "development_stage_ontology_term_id", - "disease", "disease_ontology_term_id", "donor_id", "is_primary_data", "organism", "organism_ontology_term_id", - "self_reported_ethnicity", "self_reported_ethnicity_ontology_term_id", - "sex", "sex_ontology_term_id", "tissue", "tissue_ontology_term_id", "citation", - "collection_id", "dataset_version_id", "default_embedding", "published_at", "raw_data_location", - "revised_at", "primary_cell_count", "schema_version", "tissue_type", "title", - "tombstone", "x_approximate_distribution", "explorer_url", "cell_count", "feature_count", - "filesize", "filetype", "mean_genes_per_cell", "suspension_type", "url" - ) - - # CellNexus metadata (smaller file for Shiny): drop heavy / internal columns by name patterns - drop_cellnexus <- unique(c( - intersect(remove_cols, cols), - cols[grepl("metacell", cols)] - )) - keep_cellnexus <- setdiff(cols, drop_cellnexus) - select_cellnexus <- paste(sql_id(keep_cellnexus), collapse = ", ") - - # MODIFY HERE: output cellnexus metadata parquet path - DBI::dbExecute( - con, - glue::glue( - " - COPY ( - SELECT {DBI::SQL(select_cellnexus)} - FROM metadata - ) - TO {DBI::dbQuoteString(con, file.path(out_dir, 'cellnexus_metadata.2.2.1.parquet'))} - (FORMAT PARQUET, COMPRESSION 'brotli'); - " - ) - ) - - # Original census-like metadata subset (stable columns) - census_cols <- intersect( - c( - "observation_joinid", "dataset_id", "sample_id", "cell_type", - "cell_type_ontology_term_id", "assay", "assay_ontology_term_id", "development_stage", "development_stage_ontology_term_id", - "disease", "disease_ontology_term_id", "donor_id", "is_primary_data", "organism", "organism_ontology_term_id", - "self_reported_ethnicity", "self_reported_ethnicity_ontology_term_id", - "sex", "sex_ontology_term_id", "tissue", "tissue_ontology_term_id", - "data_driven_ensemble", "ensemble_joinid", "observation_originalid", "citation", - "collection_id", "dataset_version_id", "default_embedding", "published_at", "raw_data_location", - "revised_at", "primary_cell_count", "schema_version", "tissue_type", "title", - "tombstone", "x_approximate_distribution", "explorer_url", "cell_count", "feature_count", - "filesize", "filetype", "mean_genes_per_cell", "suspension_type", "url" - ), - cols - ) - select_census <- paste(sql_id(census_cols), collapse = ", ") - - # MODIFY HERE: output census metadata parquet path - DBI::dbExecute( - con, - glue::glue( - " - COPY ( - SELECT {DBI::SQL(select_census)} - FROM metadata - ) - TO {DBI::dbQuoteString(con, file.path(out_dir, 'census_cell_metadata.2.2.1.parquet'))} - (FORMAT PARQUET, COMPRESSION 'brotli'); - " - ) - ) - - # # Metacell metadata subset - # metacell_cols <- unique(c("cell_id", "sample_id", "dataset_id", cols[grepl("metacell", cols)])) - # metacell_cols <- intersect(metacell_cols, cols) - # select_metacell <- paste(sql_id(metacell_cols), collapse = ", ") - # - # # MODIFY HERE: output metacell metadata parquet path - # DBI::dbExecute( - # con, - # glue::glue( - # " - # COPY ( - # SELECT {DBI::SQL(select_metacell)} - # FROM metadata - # ) - # TO {DBI::dbQuoteString(con, file.path(out_dir, 'metacell_metadata.2.2.1.parquet'))} - # (FORMAT PARQUET, COMPRESSION 'brotli'); - # " - # ) - # ) - - print("Done.") -}) - - -# (Optional) Check whether cellnexus_metadata parquet can be optimised further -# source("~/git_control/cellNexus/dev/data_optimisation_script.R") - diff --git a/dev/cellnexus-2024-scripts/step9_aggregate_metacell_from_sce.R b/dev/cellnexus-2024-scripts/step9_aggregate_metacell_from_sce.R deleted file mode 100644 index 09dd4a3..0000000 --- a/dev/cellnexus-2024-scripts/step9_aggregate_metacell_from_sce.R +++ /dev/null @@ -1,301 +0,0 @@ -library(targets) -# aggregate metacell from metadata and save assays -store_file_cellNexus = "/vast/scratch/users/shen.m/targets_prepare_database_split_datasets_chunked_1_4_0_metacell/" # MODIFY HERE: targets store directory for this pipeline -cell_metadata_path = "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/metadata.2.2.0.parquet" # MODIFY HERE: cell metadata parquet (used to dynamically derive metacell column names) - -tar_script({ - library(dplyr) - library(magrittr) - library(tibble) - library(targets) - library(tarchetypes) - library(crew) - library(crew.cluster) - library(tidySingleCellExperiment) - library(SingleCellExperiment) - library(tidyverse) - library(glue) - library(digest) - library(scater) - library(arrow) - library(dplyr) - library(duckdb) - library(cellNexus) - library(BiocParallel) - library(parallelly) - # Helper (optional) to avoid repetition - new_elastic <- function(name, mem_gb, time_min, workers, crashes_max, cpus_per_task = 2, backup = NULL) { - crew_controller_slurm( - name = name, - workers = workers, - crashes_max = crashes_max, - seconds_idle = 30, - options_cluster = crew_options_slurm( - memory_gigabytes_required = mem_gb, - cpus_per_task = cpus_per_task, - time_minutes = time_min - ), - backup = backup - ) - } - - # Small → large, with fallbacks to the next size up - elastic_160 <- new_elastic("elastic_160", 160, 60 * 24, workers = 8, crashes_max = 2) - elastic_120 <- new_elastic("elastic_120", 120, 60 * 4, workers = 16, crashes_max = 1, cpus_per_task = 8, backup = elastic_160) - elastic_80 <- new_elastic("elastic_80", 80, 60 * 4, workers = 24, crashes_max = 1, cpus_per_task = 8, backup = elastic_120) - elastic_40 <- new_elastic("elastic_40", 40, 60 * 4, workers = 32, crashes_max = 1, cpus_per_task = 8, backup = elastic_80) - elastic_20 <- new_elastic("elastic_20", 20, 60 * 4, workers = 48, crashes_max = 1, cpus_per_task = 8, backup = elastic_40) - elastic_10 <- new_elastic("elastic_10", 10, 60 * 4, workers = 150, crashes_max = 6, cpus_per_task = 8, backup = elastic_20) - - elastic_5_minimal <- new_elastic("elastic_5_minimal", 5, 60 * 4, workers = 300, crashes_max = 6, cpus_per_task = 2, backup = elastic_10) - - - # Group for targets (small → large) - controllers <- crew_controller_group( - elastic_10, elastic_20, elastic_40, elastic_80, elastic_120, elastic_160, elastic_5_minimal - ) - tar_option_set( - memory = "transient", - garbage_collection = 100, - storage = "worker", - retrieval = "worker", - error = "continue", - cue = tar_cue(mode = "never"), - format = "qs", - - workspace_on_error = TRUE, - controller = controllers, - trust_object_timestamps = TRUE - #workspaces = "dataset_id_sce_52dbec3c15f98d66" - ) - - get_ids <- function(cell_metadata, metacell_column) { - metacell_column <- as.character(metacell_column) - tbl(dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql(glue("SELECT * FROM read_parquet('{cell_metadata}')"))) |> - filter(!is.na(.data[[metacell_column]])) |> - distinct(file_id_cellNexus_single_cell) |> pull() - } - - get_sce <- function(cell_metadata, id, metacell_column, cache) { - metacell_column <- as.character(metacell_column) - tbl(dbConnect(duckdb::duckdb(), - dbdir = ":memory:"), sql(glue("SELECT * FROM read_parquet('{cell_metadata}')"))) |> - filter(!is.na(.data[[metacell_column]])) |> - filter(empty_droplet == FALSE, alive==TRUE, scDblFinder.class!="doublet") |> # because metacell membership ID was pre-calculated after QC - filter(file_id_cellNexus_single_cell == id) |> - select(cell_id, sample_id, dataset_id, donor_id, file_id_cellNexus_single_cell, cell_type, atlas_id, !!metacell_column) |> - get_single_cell_experiment(cache_directory = cache, repository = NULL) # this assume SCE are not uploaded to cloud - } - - # aggregate_metacell <- function(sce, metacell) { - # cores = as.numeric(Sys.getenv("SLURM_CPUS_PER_TASK", unset = 1)) - 1 - # bp <- MulticoreParam(workers = cores, progressbar = TRUE) - # aggregate_metacell <- aggregateAcrossCells(sce, colData(sce)[, c("sample_id", metacell)], BPPARAM = bp) - # aggregate_metacell = aggregate_metacell |> mutate(cell_id = paste(sample_id, .data[[metacell]], sep = "___")) - # # Assign cell_id to SCE metadata rownames - # rownames(colData(aggregate_metacell)) <- aggregate_metacell$cell_id - # aggregate_metacell = aggregate_metacell |> select(-contains(".1")) - # aggregate_metacell - # - # } - aggregate_metacell <- function(sce, metacell) { - - cores <- max( - 1, - as.numeric(Sys.getenv("SLURM_CPUS_PER_TASK", unset = 1)) - 1 - ) - - bp <- BiocParallel::MulticoreParam( - workers = cores, - progressbar = TRUE - ) - - # aggregate by sample_id + metacell - agg <- scuttle::aggregateAcrossCells( - sce, - ids = colData(sce)[, c("sample_id", metacell), drop = FALSE], - BPPARAM = bp - ) - - # make compact unique cell IDs: _ - cd <- as.data.frame(SummarizedExperiment::colData(agg)) - - cd <- cd |> - dplyr::group_by(.data[[metacell]]) |> - dplyr::mutate( - metacell_id = paste0(.data[[metacell]], "_", dplyr::row_number()) - ) |> - dplyr::ungroup() |> - select(-original_cell_) - - # put back colData - SummarizedExperiment::colData(agg) <- S4Vectors::DataFrame(cd) - - # use compact .cell as column names / colData rownames - colnames(agg) <- cd$metacell_id - rownames(SummarizedExperiment::colData(agg)) <- cd$metacell_id - - # optional: remove duplicated helper columns created by aggregation - keep_cols <- !grepl("\\.1$", colnames(SummarizedExperiment::colData(agg))) - SummarizedExperiment::colData(agg) <- SummarizedExperiment::colData(agg)[, keep_cols, drop = FALSE] - - agg - } - - save_anndata = function(sce, cache_directory) { - dir.create(cache_directory, showWarnings = FALSE, recursive = TRUE) - file_id = pull(distinct(sce, file_id_cellNexus_single_cell)) - cellNexus:::save_sce_as_h5ad(sce, glue("{cache_directory}/{file_id}"), mode = "w") - return(TRUE) - } - - c( - list( - tar_target( - cell_metadata, - "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/metadata.2.2.0.parquet", # MODIFY HERE: cell metadata parquet (must match cell_metadata_path above) - deployment = "main", - packages = c("arrow", "dplyr", "duckdb") - ), - tar_target( - local_cache, - "/vast/scratch/users/shen.m/cellNexus", # MODIFY HERE: local cache directory containing the single-cell h5ad files (input to get_single_cell_experiment) - deployment = "main" - ), - tar_target( - save_cache_directory, - "/vast/scratch/users/shen.m/cellNexus/cellxgene_2024/0.2.0", # MODIFY HERE: output directory where aggregated metacell anndata files are saved - deployment = "main" - ) - ), - tarchetypes::tar_map( - values = tibble( - metacell_column = tbl( - dbConnect(duckdb::duckdb(), dbdir = ":memory:"), - sql("SELECT * FROM read_parquet('/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/metadata.2.2.0.parquet')") # MODIFY HERE: cell metadata parquet path inside SQL (must match cell_metadata_path above) - ) |> select(contains("metacell_")) |> colnames() - ), - names = metacell_column, - unlist = TRUE, - tar_target( - file_ids, - get_ids(cell_metadata, metacell_column) - # |> - # # TEST PURPOSE ONLY - # head(2) - ), - tar_target( - file_id_sce, - get_sce(cell_metadata, file_ids, metacell_column, local_cache), - pattern = map(file_ids), - resources = tar_resources(crew = tar_resources_crew(controller = "elastic_10")) - ), - tar_target( - metacell, - aggregate_metacell(file_id_sce, metacell_column), - pattern = map(file_id_sce), - resources = tar_resources(crew = tar_resources_crew(controller = "elastic_20")) - ), - tar_target( - save_metacell, - save_anndata(metacell, paste0(save_cache_directory, "/", metacell_column, "/counts")), - pattern = map(metacell), - resources = tar_resources(crew = tar_resources_crew(controller = "elastic_20")) - ) - ) - ) - -}, script = paste0(store_file_cellNexus, "_target_script.R"), ask = FALSE) - -job::job({ - - tar_make( - script = paste0(store_file_cellNexus, "_target_script.R"), - store = store_file_cellNexus, - reporter = "summary" #, callr_function = NULL - ) -}) -#tar_invalidate(names = everything(), store = store_file_cellNexus) -tar_meta(store = store_file_cellNexus) |> filter(!is.na(error)) |> distinct(name, error) -tar_errored(store = store_file_cellNexus) -# With tar_map, target names are suffixed by metacell_column, e.g. file_ids_metacell_4, metacell_metacell_4 -# tar_workspace("metacell_metacell_4_", store = store_file_cellNexus, script = paste0(store_file_cellNexus, "_target_script.R")) -# debugonce(aggregate_metacell) -# aggregate_metacell(file_id_sce, "metacell_4") # when debugging a specific branch - -# Check the number of file id should be created for metacell_2 -cache_dir = "/vast/projects/cellxgene_curated/metadata_cellxgene_mengyuan/" # MODIFY HERE: directory used for verification queries below (contains metadata.2.2.0.parquet) -# Define all metacell levels -metacell_levels <- c(2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768, 65536) -metacell_names <- paste0("metacell_", metacell_levels) - -# Function to get file IDs that SHOULD be generated (from metadata) -file_ids <- function(metacell_col_name) { - get_metadata( - cache_directory = cache_dir, - cloud_metadata = NULL, - local_metadata = file.path(cache_dir, "metadata.2.2.0.parquet") - ) |> - filter( - !is.na(.data[[metacell_col_name]]), - empty_droplet == FALSE, - alive == TRUE, - scDblFinder.class != "doublet" - ) |> - distinct(file_id_cellNexus_single_cell) |> - pull(file_id_cellNexus_single_cell) |> - sort() -} - -# Function to get file names that ACTUALLY exist in subfolder -file_names_saved <- function(metacell_name) { - subfolder <- file.path("~/scratch/cellNexus/cellxgene/01-07-2024/", metacell_name, "counts") - if (dir.exists(subfolder)) { - list.files(subfolder, recursive = FALSE) |> sort() - } else { - character(0) - } -} - -# Build the summary tibble -result <- tibble(metacell = metacell_names) |> - mutate( - ids_to_save = map(metacell, file_ids, .progress = "Getting expected file IDs"), - ids_saved = map(metacell, file_names_saved, .progress = "Getting actual file names"), - file_to_save = map_int(ids_to_save, length), - file_saved = map_int(ids_saved, length), - missing = map2(ids_to_save, ids_saved, \(expected, actual) setdiff(expected, actual)), - extra = map2(ids_to_save, ids_saved, \(expected, actual) setdiff(actual, expected)) - ) |> - select(metacell, file_to_save, file_saved, missing, extra) - -result - -# Unit test query lung tissue, metacell 256 -lung_metacell_256 = get_metadata( - cache_directory = cache_dir, - cloud_metadata = NULL, - local_metadata = file.path(cache_dir, "metadata.2.2.0.parquet") -) |> - filter(!is.na(metacell_256), - empty_droplet == FALSE, - alive == TRUE, - scDblFinder.class != "doublet") |> - filter(tissue == "lung") |> - get_metacell(cache_directory = cache_dir, cell_aggregation = "metacell_256") - -# Cell number in lung_metacell_256 should match the count below: -get_metadata( - cache_directory = cache_dir, - cloud_metadata = NULL, - local_metadata = file.path(cache_dir, "metadata.2.2.0.parquet") -) |> - filter(!is.na(metacell_256), - empty_droplet == FALSE, - alive == TRUE, - scDblFinder.class != "doublet") |> - filter(tissue == "lung") |> - distinct(sample_id, metacell_256, cell_type_unified_ensemble) |> dplyr::count() - - diff --git a/man/HPCell-package.Rd b/man/HPCell-package.Rd index 6426254..ffb4967 100644 --- a/man/HPCell-package.Rd +++ b/man/HPCell-package.Rd @@ -9,10 +9,11 @@ Massively-parallel R native pipeline for single-cell analysis. } \author{ -\strong{Maintainer}: Stefano Mangiola \email{mangiolastefano@gmail.com} +\strong{Maintainer}: Mengyuan Shen \email{shen.m@wehi.edu.au} Authors: \itemize{ + \item Stefano Mangiola \email{mangiolastefano@gmail.com} \item Jiayi Si \email{si.j@wehi.edu.au} } diff --git a/man/initialise_hpc.Rd b/man/initialise_hpc.Rd index 67738d7..f1ae12b 100644 --- a/man/initialise_hpc.Rd +++ b/man/initialise_hpc.Rd @@ -8,6 +8,7 @@ initialise_hpc( input_hpc, store = targets::tar_config_get("store"), computing_resources = crew_controller_local(workers = 1), + default_controller = NULL, tier = rep(1, length(input_hpc)), debug_step = NULL, RNA_assay_name = "RNA", @@ -21,29 +22,33 @@ initialise_hpc( ) } \arguments{ -\item{input_hpc}{Named character vector of paths to input data files, one -element per sample. If names are not set, integer indices are used.} +\item{input_hpc}{Character vector of input data path for the pipeline.} -\item{store}{Directory path where pipeline files and targets store are written.} +\item{store}{Directory path for storing the pipeline files.} \item{computing_resources}{A \code{crew} controller object (or list of controllers) specifying the computing back-end. Defaults to a local single-worker controller.} +\item{default_controller}{Optional character name of the default \code{crew} +controller to use for targets that do not specify their own controller. +Passed to \code{targets::tar_resources(crew = tar_resources_crew(controller = ...))}. +\code{NULL} uses targets defaults.} + \item{tier}{Integer vector (same length as \code{input_hpc}) assigning each sample to a processing tier for tiered execution. Default: all samples in tier 1.} -\item{debug_step}{Character name of a single target to debug; passed to -\code{targets::tar_option_set(debug = ...)}. \code{NULL} disables debugging.} +\item{debug_step}{Optional step for debugging.} -\item{RNA_assay_name}{Name of the RNA assay in the input Seurat/SCE object.} +\item{RNA_assay_name}{Name of the RNA assay.} -\item{gene_nomenclature}{Character scalar indicating gene identifier type in -the input data. One of \code{"symbol"} or \code{"ensembl"}.} +\item{gene_nomenclature}{Character vector indicating gene nomenclature in input_data} -\item{data_container_type}{Character scalar specifying the input data format. -Accepted values: \code{"sce_rds"} (SingleCellExperiment RDS), -\code{"seurat_rds"} (Seurat RDS), \code{"sce_hdf5"} (HDF5-backed SCE), -\code{"seurat_h5"} (HDF5-backed Seurat).} +\item{data_container_type}{A character vector of length one specifies the input data type. +The accepted input data type are: +sce_rds for \code{SingleCellExperiment} RDS, +seurat_rds for \code{Seurat} RDS, +sce_hdf5 for \code{SingleCellExperiment} HDF5-based object +seurat_h5 for \code{Seurat} HDF5-based object} \item{verbosity}{Reporter string passed to \code{targets::tar_make()}. Defaults to the current targets configuration value.} @@ -54,19 +59,16 @@ the current targets configuration value.} \item{update}{Cue mode string for \code{targets::tar_cue()}, controlling when targets are re-run. Default: \code{"thorough"}.} -\item{garbage_collection}{Numeric interval (in targets) at which R garbage -collection is triggered during the pipeline run. Default: \code{0} (disabled).} +\item{garbage_collection}{Numeric interval at which R garbage collection is +triggered during the pipeline run. Default: \code{0} (disabled).} \item{workspace_on_error}{Logical; if \code{TRUE}, saves a workspace snapshot when a target errors. Default: \code{FALSE}.} } \value{ -An \code{HPCell} S3 object containing the initialisation arguments, ready -to be extended with pipeline step functions. +The output of the \code{targets} pipeline, typically a pre-processed data set. } \description{ -Sets up and writes a \code{targets} pipeline script for HPCell. Saves input data -and configuration to disk, then returns an \code{HPCell} object that downstream -grammar functions (e.g. \code{remove_empty_DropletUtils}, \code{evaluate_hpc}) can -extend before the pipeline is executed with \code{evaluate_hpc()}. +This function sets up and executes a \code{targets} pipeline for HPCell. It saves input data and configurations, +writes a pipeline script, and runs the pipeline using the 'targets' package. }