From 06576ea79eeaca3915f5684f562a598396f0a45a Mon Sep 17 00:00:00 2001 From: Chen Zhan Date: Mon, 7 Apr 2025 10:14:52 +0930 Subject: [PATCH 01/12] add additional target to summarise converged param rename function param name set error = null for some targets --- .../de_brms_whole_immune_system.R | 37 ++++++++++++++----- 1 file changed, 28 insertions(+), 9 deletions(-) diff --git a/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R b/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R index cbee7b3..bcf8685 100644 --- a/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R +++ b/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R @@ -239,12 +239,12 @@ tar_script({ bind_cols(fitted_values_ethnicity_tbl) } -get_adjusted_matrix = function(summary_df, column_adjusted){ +get_adjusted_matrix = function(effect_removed_df, column_adjusted){ column_adjusted = enquo(column_adjusted) m = - summary_df |> + effect_removed_df |> unnest(!!column_adjusted) |> # dplyr::filter(analysis == "observed_proportion") |> select(.feature, adjusted___Estimate, sample_id) |> @@ -269,7 +269,7 @@ get_adjusted_matrix = function(summary_df, column_adjusted){ return(m) } - + #-----------------------# # Pipeline #-----------------------# @@ -617,8 +617,8 @@ get_adjusted_matrix = function(summary_df, column_adjusted){ pattern = map(se_df), packages = c( "brms", "glue", "stringr", "dplyr", "purrr", "SummarizedExperiment", "tidySummarizedExperiment"), resources = tar_resources(crew = tar_resources_crew("elastic")), - cue = tar_cue(mode = "never") - + cue = tar_cue(mode = "never"), + error = "null" ), ## summary ----- @@ -627,7 +627,7 @@ get_adjusted_matrix = function(summary_df, column_adjusted){ tar_target( summary, estimates_chunk |> - mutate(summary = map(brms_fit, ~ .x |> hypothesis( + mutate(summary_ethnicity = map(brms_fit, ~ .x |> hypothesis( c( "Europeans" = "(ethnicity_groupsAfrican + ethnicity_groupsEastAsian @@ -772,10 +772,11 @@ get_adjusted_matrix = function(summary_df, column_adjusted){ pattern = map(estimates_chunk), packages = c( "brms", "glue", "dplyr", "purrr", "rstan", "magrittr", "stringr"), - resources = tar_resources(crew = tar_resources_crew("elastic")) + resources = tar_resources(crew = tar_resources_crew("elastic")), + error = "null" ), - ## effect_removed ----- + ## effect_removed ----- # This target generates adjusted model estimates by removing unwanted effects from the fitted Bayesian models, # thereby isolating the effects of interest. Here, nuisance covariates are set to NA and removed from the predictions. # This target produces adjusted estimates from the Bayesian models, removing unwanted effects while retaining @@ -863,9 +864,27 @@ get_adjusted_matrix = function(summary_df, column_adjusted){ pattern = map(estimates_chunk), packages = c( "brms", "glue", "dplyr", "purrr", "rstan"), - resources = tar_resources(crew = tar_resources_crew("elastic")) + resources = tar_resources(crew = tar_resources_crew("elastic")), + error = "null" ), + ## param ----- + tar_target( + param, + estimates_chunk %>% + mutate( + param = map( + brms_fit, + ~ summary(.x$fit) |> as.data.frame() + ) + )%>% + select(-brms_fit), + pattern = map(estimates_chunk), + packages = c( "brms", "dplyr", "purrr", "rstan"), + resources = tar_resources(crew = tar_resources_crew("elastic")), + error = "null" + ), + # adjusted_matrix ----- tar_target( adjusted_assay_ethnicity, From c7b3f9e28908a1df95ec4c484d4fcb4eb823168d Mon Sep 17 00:00:00 2001 From: Chen Zhan Date: Fri, 11 Apr 2025 22:11:57 +0930 Subject: [PATCH 02/12] revise prior to dynamically incorporate intecept --- .../de_brms_whole_immune_system.R | 579 ++++++++++-------- 1 file changed, 328 insertions(+), 251 deletions(-) diff --git a/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R b/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R index bcf8685..79778a9 100644 --- a/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R +++ b/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R @@ -11,7 +11,7 @@ tar_script({ library(qs) library(crew) library(crew.cluster) - + # Set file path ----- hdf5_path = "/hpcfs/groups/phoenix-hpc-mangiola_laboratory/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseudobulk_sample_is_immune" metadata_path = "/hpcfs/groups/phoenix-hpc-mangiola_laboratory/Mangiola_ImmuneAtlas/taskforce_shared_folder/cell_metadata_1_0_6_sccomp_input_counts.rds" @@ -84,7 +84,7 @@ tar_script({ #-----------------------# # Functions #-----------------------# - + #' Remove Unwanted Effects from a brmsfit Model #' #' This function calculates posterior residuals from a \code{brmsfit} model and combines them with @@ -239,7 +239,7 @@ tar_script({ bind_cols(fitted_values_ethnicity_tbl) } -get_adjusted_matrix = function(effect_removed_df, column_adjusted){ + get_adjusted_matrix = function(effect_removed_df, column_adjusted){ column_adjusted = enquo(column_adjusted) @@ -269,7 +269,7 @@ get_adjusted_matrix = function(effect_removed_df, column_adjusted){ return(m) } - + #-----------------------# # Pipeline #-----------------------# @@ -300,111 +300,113 @@ get_adjusted_matrix = function(effect_removed_df, column_adjusted){ # pseudobulk_sample ------ pseudobulk_sample, { - message('TAR: pseudobulk_sample START') - se = - loadHDF5SummarizedExperiment(hdf5_path) |> - filter(is_gene_shared) |> - - #---------------------------------# - # Edit or add more filters here for analyses - #---------------------------------# - filter(is_immune & do_analyse) - - # TEMPORARY BECAUSE I FORGOT TO INTEGRATE AGE BINS - se = se |> - left_join( - readRDS(metadata_path) |> - distinct(sample_id, age_days, age_bin) - ) - - # Filter common genes - se = se[((assay(se, "gene_presence") > 0) |> rowSums() > (ncol(se) * 0.95)),,drop=FALSE ] - - # Filter samples that have enough genes > 0 but not too many - samples_with_right_number_of_detected_genes = - (se |> assay() > 0) |> - colSums() |> - divide_by(nrow(se)) |> - dplyr::between(0.3, 1) - - se = se[,samples_with_right_number_of_detected_genes] - - # Compute mean library size - mean_library_size <- se |> - assay("counts") |> - _[nrow(se) |> seq_len() |> sample(size = 2000), ] |> - colSums() |> - mean() - - # Optional: retrieve the sample name (column name in the SummarizedExperiment) - reference_sample <- colnames(se)[ - se |> - assay("counts") |> - colSums() |> - {\(x) abs(x - mean_library_size)}() |> # Calculate absolute difference from the mean - which.min() # Identify the smallest difference - ] - - message('TAR: pseudobulk_sample Phase2') - se = - se |> - keep_abundant(design = - se |> - - # Discretise the age for the following operation - mutate(is_old_individual = age_days > 50*365) |> - - # This is to resolve some confounders to preserve the genes. - # In this case we care about data variability, not the actual meaning of the variables - resolve_complete_confounders_of_non_interest(tissue_groups, sex, ethnicity_groups, is_old_individual) |> - colData() |> - droplevels() |> - model.matrix(~ tissue_groups + sex___altered + ethnicity_groups___altered + is_old_individual___altered, data = _ ), - minimum_counts = 100 - ) |> - - # Get scaling factor - scale_abundance(method = "TMMwsp", reference_sample = reference_sample) |> - - # Drop sex unknown as causes problem during fit - mutate( - sex = if_else(sex |> is.na(), "unknown", sex), - ethnicity_groups = if_else(ethnicity_groups |> is.na(), "Other/Unknown", ethnicity_groups) - ) |> - filter(sex != "unknown") |> - filter(!age_bin |> is.na()) |> - - # Eliminate complete confounders - tidybulk:::resolve_complete_confounders_of_non_interest(assay_groups, dataset_id, disease_groups) |> - - # sibrary size factor is the reciproque of the multiplier (correction factor) - mutate(offset = log(1/multiplier)) |> - - # Set intercept - mutate( - ethnicity_groups = fct_relevel(ethnicity_groups, "European"), - assay_groups___altered = fct_relevel(assay_groups___altered, "10x Genomics 3"), - disease_groups___altered = fct_relevel(disease_groups___altered, "Normal"), - age_bin = fct_relevel(age_bin, "Adolescence") - ) - - # # Add dispersion - # rowData(se) = - # rowData(se) |> - # as_tibble(rownames = ".feature") |> - # left_join(glmGamPoi_overdispersions |> enframe(name = ".feature", value = "dispersion")) |> - # data.frame(row.names = ".feature") |> DataFrame() - - message('TAR: pseudobulk_sample COMPLETE') - se + # message('TAR: pseudobulk_sample START') + # se = + # loadHDF5SummarizedExperiment(hdf5_path) |> + # filter(is_gene_shared) |> + # + # #---------------------------------# + # # Edit or add more filters here for analyses + # #---------------------------------# + # filter(is_immune & do_analyse) + # + # # TEMPORARY BECAUSE I FORGOT TO INTEGRATE AGE BINS + # se = se |> + # left_join( + # readRDS(metadata_path) |> + # distinct(sample_id, age_days, age_bin) + # ) + # + # # Filter common genes + # se = se[((assay(se, "gene_presence") > 0) |> rowSums() > (ncol(se) * 0.95)),,drop=FALSE ] + # + # # Filter samples that have enough genes > 0 but not too many + # samples_with_right_number_of_detected_genes = + # (se |> assay() > 0) |> + # colSums() |> + # divide_by(nrow(se)) |> + # dplyr::between(0.3, 1) + # + # se = se[,samples_with_right_number_of_detected_genes] + # + # # Compute mean library size + # mean_library_size <- se |> + # assay("counts") |> + # _[nrow(se) |> seq_len() |> sample(size = 2000), ] |> + # colSums() |> + # mean() + # + # # Optional: retrieve the sample name (column name in the SummarizedExperiment) + # reference_sample <- colnames(se)[ + # se |> + # assay("counts") |> + # colSums() |> + # {\(x) abs(x - mean_library_size)}() |> # Calculate absolute difference from the mean + # which.min() # Identify the smallest difference + # ] + # + # message('TAR: pseudobulk_sample Phase2') + # se = + # se |> + # keep_abundant(design = + # se |> + # + # # Discretise the age for the following operation + # mutate(is_old_individual = age_days > 50*365) |> + # + # # This is to resolve some confounders to preserve the genes. + # # In this case we care about data variability, not the actual meaning of the variables + # resolve_complete_confounders_of_non_interest(tissue_groups, sex, ethnicity_groups, is_old_individual) |> + # colData() |> + # droplevels() |> + # model.matrix(~ tissue_groups + sex___altered + ethnicity_groups___altered + is_old_individual___altered, data = _ ), + # minimum_counts = 100 + # ) |> + # + # # Get scaling factor + # scale_abundance(method = "TMMwsp", reference_sample = reference_sample) |> + # + # # Drop sex unknown as causes problem during fit + # mutate( + # sex = if_else(sex |> is.na(), "unknown", sex), + # ethnicity_groups = if_else(ethnicity_groups |> is.na(), "Other/Unknown", ethnicity_groups) + # ) |> + # filter(sex != "unknown") |> + # filter(!age_bin |> is.na()) |> + # + # # Eliminate complete confounders + # tidybulk:::resolve_complete_confounders_of_non_interest(assay_groups, dataset_id, disease_groups) |> + # + # # sibrary size factor is the reciproque of the multiplier (correction factor) + # mutate(offset = log(1/multiplier)) |> + # + # # Set intercept + # mutate( + # ethnicity_groups = fct_relevel(ethnicity_groups, "European"), + # assay_groups___altered = fct_relevel(assay_groups___altered, "10x Genomics 3"), + # disease_groups___altered = fct_relevel(disease_groups___altered, "Normal"), + # age_bin = fct_relevel(age_bin, "Adolescence") + # ) + # + # # # Add dispersion + # # rowData(se) = + # # rowData(se) |> + # # as_tibble(rownames = ".feature") |> + # # left_join(glmGamPoi_overdispersions |> enframe(name = ".feature", value = "dispersion")) |> + # # data.frame(row.names = ".feature") |> DataFrame() + # + # message('TAR: pseudobulk_sample COMPLETE') + # se + # load process data to save time when testing + loadHDF5SummarizedExperiment('/hpcfs/groups/phoenix-hpc-mangiola_laboratory/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseduobulk_sample_tar_load_altered/') }, packages = c("tidybulk", "HDF5Array", "tidySummarizedExperiment", "magrittr", "tibble", "forcats"), - resources = tar_resources(crew = tar_resources_crew("elastic_big")), + resources = tar_resources(crew = tar_resources_crew("elastic_big_30_cores")), memory = "persistent", error = "stop" ), - + # pseudobulk_sample_id ------ # This target extracts unique sample ids from the pseudobulk sample tar_target( @@ -422,11 +424,7 @@ get_adjusted_matrix = function(effect_removed_df, column_adjusted){ pseudobulk_sample |> distinct(.feature)|> # testing genes that ran for long time - # filter(.feature %in% c( - # "ENSG00000175274", "ENSG00000213221", "ENSG00000101405", "ENSG00000003756", - # "ENSG00000236859", "ENSG00000104825", "ENSG00000203497", "ENSG00000143110", - # "ENSG00000077454", "ENSG00000104231" - # )) |> + # filter(.feature %in% readRDS('/hpcfs/groups/phoenix-hpc-mangiola_laboratory/Mangiola_ImmuneAtlas/ning_data/ethnicity_umap_selected_genes.rds')) |> # head(1) %>% group_by(.feature) |> tar_group(), @@ -519,20 +517,27 @@ get_adjusted_matrix = function(effect_removed_df, column_adjusted){ # Using the externally, eBayes inferred overdispersion # shape ~ 1 + offset(log(1/dispersion)) - ) + ) - # prior = c( - # prior(normal(i, 5), class = Intercept), - # prior(normal(0, 2), class = Intercept, dpar = shape), - # prior(normal(0, 5), class = b), - # prior(normal(0, 2), class = b, dpar = shape) - # ) |> - # substitute(env = list(i = mean(log1p(data$counts / exp(data$offset))))) |> - # eval() - + # prior Version 0 + prior = c( + prior(normal(i, 5), class = Intercept), + prior(normal(0, 2), class = Intercept, dpar = shape), + prior(normal(0, 5), class = b), + prior(normal(0, 2), class = b, dpar = shape) + # prior(beta(0.5381488, 10.3577433), class = "zi", lb = 0, ub = 1) # addition zi from V2 + ) |> + substitute(env = list(i = mean(log1p(data$counts / exp(data$offset))))) |> + eval() + # + # chains = 2 + # inits <- list(Intercept = mean(log1p(data$counts / exp(data$offset)))) + # inits <- replicate(chains, inits, simplify = FALSE) + + # test: # prior = prior(normal(-0.0002056948, 0.07690437)) - # HPC pipeline: param V1: + # HPC pipeline: param V1: learned from version 0 # prior = c( # prior(student_t(4.45496, 0.008599254, 1.143344), class = "b"), # prior(student_t(18.16242, 0.07952513, 0.9926044), class = "b", dpar = "shape"), @@ -543,56 +548,127 @@ get_adjusted_matrix = function(effect_removed_df, column_adjusted){ # prior(beta(0.5541155, 9.337894), class = "zi") # ) - # HPC pipeline: param V2: - prior = c( - prior(student_t(6.153327, 0.06161134, 0.9263627), class = "b"), - prior(student_t(40.51669, 0.07603337, 0.8252114), class = "b", dpar = "shape"), - prior(normal(6.057503, 2.438534), class = "Intercept"), - prior(normal(0.4260793, 1.470536), class = "Intercept", dpar = "shape"), - prior(student_t(52.19541 , 0.5703259, 0.4147664), class = "sd", lb = 0), - prior(normal(0.8670409, 0.1779553), class = "sd", dpar = "shape", lb = 0), - prior(beta(0.5381488, 10.3577433), class = "zi", lb = 0, ub = 1) - ) - - chains = 2 + # HPC pipeline: param V2: updated from v1 and set inits + # prior = c( + # prior(student_t(6.153327, 0.06161134, 0.9263627), class = "b"), + # prior(student_t(40.51669, 0.07603337, 0.8252114), class = "b", dpar = "shape"), + # prior(normal(6.057503, 2.438534), class = "Intercept"), + # prior(normal(0.4260793, 1.470536), class = "Intercept", dpar = "shape"), + # prior(student_t(52.19541 , 0.5703259, 0.4147664), class = "sd", lb = 0), + # prior(normal(0.8670409, 0.1779553), class = "sd", dpar = "shape", lb = 0), + # prior(beta(0.5381488, 10.3577433), class = "zi", lb = 0, ub = 1) + # ) + # + # chains = 2 + # + # Kc <- 39 + # Kc_shape <- 28 + # M_1 <- 1; N_1 <- 105 + # M_2 <- 19; N_2 <- 26 + # M_3 <- 1; N_3 <- 26 + # + # inits <- lapply(1:chains, function(i) { + # list( + # # Fixed effects for count part + # b = 0.06161134 + 0.9263627 * rt(Kc, 6.153327), + # Intercept = rnorm(1, 6.057503, 2.438534), + # + # # Fixed effects for shape submodel + # b_shape = 0.07603337 + 0.8252114 * rt(Kc_shape, 40.51669), + # Intercept_shape = rnorm(1, 0.4260793, 1.470536), + # + # # Zero-inflation probability + # zi = rbeta(1, 0.5381488, 10.3577433), + # + # # Group-level standard deviations and effects + # sd_1 = abs(0.5703259 + 0.4147664 * rt(M_1, 52.19541)), # count + # z_1 = replicate(M_1, rnorm(N_1, mean = 0 , sd = 0.08547970), simplify = FALSE), + # + # sd_2 = abs(0.5703259 + 0.4147664 * rt(M_2, 52.19541)), # zi + # z_2 = matrix(rnorm(M_2 * N_2, mean = 0 , sd = 0.08547970), nrow = M_2, ncol = N_2), + # L_2 = diag(M_2), # no correlation (identity) + # + # sd_3 = abs(rnorm(M_3, 0.8670409, 0.1779553)), # shape + # z_3 = replicate(M_3, rnorm(N_3, mean = 0 , sd = 0.08547970), simplify = FALSE) + # ) + # }) - Kc <- 39 - Kc_shape <- 28 - M_1 <- 1; N_1 <- 105 - M_2 <- 19; N_2 <- 26 - M_3 <- 1; N_3 <- 26 + # HPC pipeline: param V3: cap df in v2 and dynamically set mu of intercept + # prior = c( + # # prior(student_t(3, 0.06161134, 0.9263627), class = "b"), + # # prior(student_t(3, 0.07603337, 0.8252114), class = "b", dpar = "shape"), + # prior(normal(i, 2.438534), class = "Intercept"), + # prior(normal(0.4260793, 1.470536), class = "Intercept", dpar = "shape"), + # prior(normal(0, 5), class = b), + # prior(normal(0, 2), class = b, dpar = shape), + # # prior(student_t(3 , 0.5703259, 0.4147664), class = "sd", lb = 0), + # # prior(normal(0.8670409, 0.1779553), class = "sd", dpar = "shape", lb = 0), + # prior(beta(0.5381488, 10.3577433), class = "zi", lb = 0, ub = 1) + # ) |> + # substitute(env = list(i = mean(log1p(data$counts / exp(data$offset))))) |> + # eval() + chains = 2 + + # dynamically extract param from stan data + # code used from brm + bterms <- brmsterms( + formula = brms:::validate_formula( + formula, data = data, family = zero_inflated_negbinomial(), + autocor = NULL, sparse = NULL, cov_ranef = NULL + ) + ) + bframe <- brms:::brmsframe(bterms, data) + sdata <- brms:::.standata( + bframe, data = data, prior = prior, + data2 = NULL, stanvars = NULL, threads = NULL + ) + + Kc <- sdata$Kc + Kc_shape <- sdata$Kc_shape + M_1 <- sdata$M_1; N_1 <- sdata$N_1 + M_2 <- sdata$M_2; N_2 <- sdata$N_2 + M_3 <- sdata$M_3; N_3 <- sdata$N_3 + inits <- lapply(1:chains, function(i) { list( + + #### revert v0 prior # Fixed effects for count part - b = 0.06161134 + 0.9263627 * rt(Kc, 6.153327), - Intercept = rnorm(1, 6.057503, 2.438534), - + b = rnorm(Kc, 0, 5), + # dynamically set mu for intercept + Intercept = rnorm(1, mean(log1p(data$counts / exp(data$offset))), 5), + # Fixed effects for shape submodel - b_shape = 0.07603337 + 0.8252114 * rt(Kc_shape, 40.51669), - Intercept_shape = rnorm(1, 0.4260793, 1.470536), - + b_shape = rnorm(Kc_shape, 0, 2), + Intercept_shape = rnorm(1, 0, 2) + + + # # Fixed effects for count part + # b = 0.06161134 + 0.9263627 * rt(Kc, 3), + # # dynamically set mu for intercept + # Intercept = rnorm(1, mean(log1p(data$counts / exp(data$offset))), 2.438534), + # + # # Fixed effects for shape submodel + # b_shape = 0.07603337 + 0.8252114 * rt(Kc_shape, 3), + # Intercept_shape = rnorm(1, 0.4260793, 1.470536) + # Zero-inflation probability - zi = rbeta(1, 0.5381488, 10.3577433), - - # Group-level standard deviations and effects - sd_1 = abs(0.5703259 + 0.4147664 * rt(M_1, 52.19541)), # count - z_1 = replicate(M_1, rnorm(N_1, mean = 0 , sd = 0.08547970), simplify = FALSE), - - sd_2 = abs(0.5703259 + 0.4147664 * rt(M_2, 52.19541)), # zi - z_2 = matrix(rnorm(M_2 * N_2, mean = 0 , sd = 0.08547970), nrow = M_2, ncol = N_2), - L_2 = diag(M_2), # no correlation (identity) - - sd_3 = abs(rnorm(M_3, 0.8670409, 0.1779553)), # shape - z_3 = replicate(M_3, rnorm(N_3, mean = 0 , sd = 0.08547970), simplify = FALSE) + # zi = rbeta(1, 0.5381488, 10.3577433) + # + # # Group-level standard deviations and effects + # sd_1 = abs(0.5703259 + 0.4147664 * rt(M_1, 3)), # count + # z_1 = replicate(M_1, rnorm(N_1, mean = 0 , sd = 0.08547970), simplify = FALSE), + # + # sd_2 = abs(0.5703259 + 0.4147664 * rt(M_2, 3)), # zi + # z_2 = matrix(rnorm(M_2 * N_2, mean = 0 , sd = 0.08547970), nrow = M_2, ncol = N_2), + # L_2 = diag(M_2), # no correlation (identity) + # + # sd_3 = abs(rnorm(M_3, 0.8670409, 0.1779553)), # shape + # z_3 = replicate(M_3, rnorm(N_3, mean = 0 , sd = 0.08547970), simplify = FALSE) ) }) - # chains = 2 - # inits <- list(Intercept = mean(log1p(data$counts / exp(data$offset)))) - # inits <- replicate(chains, inits, simplify = FALSE) - - brm( formula = formula, data = data, @@ -607,20 +683,21 @@ get_adjusted_matrix = function(effect_removed_df, column_adjusted){ #save_model = glue("{external_directory}~/temp.rds"), #algorithm = "pathfinder", init = inits, - iter = 400 # Increase iterations for better convergence + iter = 400, # Increase iterations for better convergence + sample_prior = TRUE ) })) |> - # Drop data because it is withn the brms object - select(-se), + # Drop data because it is withn the brms object + select(-se), pattern = map(se_df), packages = c( "brms", "glue", "stringr", "dplyr", "purrr", "SummarizedExperiment", "tidySummarizedExperiment"), resources = tar_resources(crew = tar_resources_crew("elastic")), cue = tar_cue(mode = "never"), error = "null" ), - + ## summary ----- # This target summarises the fitted Bayesian models by performing hypothesis tests for ethnicity contrasts # and extracting convergence diagnostics (Rhat) for the ethnicity parameters. @@ -662,70 +739,70 @@ get_adjusted_matrix = function(effect_removed_df, column_adjusted){ + `ethnicity_groupsJapanese` - 5 * ethnicity_groupsHispanicDLatinAmerican ) / 5 = 0", - - "Japanese" = "( + + "Japanese" = "( ethnicity_groupsAfrican + ethnicity_groupsHispanicDLatinAmerican + ethnicity_groupsSouthAsian + ethnicity_groupsEastAsian - 5 * `ethnicity_groupsJapanese` ) / 5 = 0" - ), - # c( - # "African" = "(ethnicity_groupsEuropean - # + ethnicity_groupsEastAsian - # + ethnicity_groupsHispanicDLatinAmerican - # + ethnicity_groupsSouthAsian - # + `ethnicity_groupsJapanese`) / 5 = 0", - # - # "Europeans" = "( - # ethnicity_groupsEastAsian - # + ethnicity_groupsHispanicDLatinAmerican - # + ethnicity_groupsSouthAsian - # + `ethnicity_groupsJapanese` - # - 5 * ethnicity_groupsEuropean - # ) / 5 = 0", - # - # "EastAsian" = "( - # ethnicity_groupsEuropean - # + ethnicity_groupsHispanicDLatinAmerican - # + ethnicity_groupsSouthAsian - # + `ethnicity_groupsJapanese` - # - 5 * ethnicity_groupsEastAsian - # ) / 5 = 0", - # - # "SouthAsian" = "( - # ethnicity_groupsEuropean - # + ethnicity_groupsHispanicDLatinAmerican - # + ethnicity_groupsEastAsian - # + `ethnicity_groupsJapanese` - # - 5 * ethnicity_groupsSouthAsian - # ) / 5 = 0", - # - # "HispanicDLatinAmerican" = "( - # ethnicity_groupsEuropean - # + ethnicity_groupsEastAsian - # + ethnicity_groupsSouthAsian - # + `ethnicity_groupsJapanese` - # - 5 * ethnicity_groupsHispanicDLatinAmerican - # ) / 5 = 0", - # - # "Japanese" = "( - # ethnicity_groupsEuropean - # + ethnicity_groupsHispanicDLatinAmerican - # + ethnicity_groupsSouthAsian - # + ethnicity_groupsEastAsian - # - 5 * `ethnicity_groupsJapanese` - # ) / 5 = 0" - # ), - - # Median instead and mad of mean and sd - robust=TRUE, - alpha = 0.1 - ) + ), + # c( + # "African" = "(ethnicity_groupsEuropean + # + ethnicity_groupsEastAsian + # + ethnicity_groupsHispanicDLatinAmerican + # + ethnicity_groupsSouthAsian + # + `ethnicity_groupsJapanese`) / 5 = 0", + # + # "Europeans" = "( + # ethnicity_groupsEastAsian + # + ethnicity_groupsHispanicDLatinAmerican + # + ethnicity_groupsSouthAsian + # + `ethnicity_groupsJapanese` + # - 5 * ethnicity_groupsEuropean + # ) / 5 = 0", + # + # "EastAsian" = "( + # ethnicity_groupsEuropean + # + ethnicity_groupsHispanicDLatinAmerican + # + ethnicity_groupsSouthAsian + # + `ethnicity_groupsJapanese` + # - 5 * ethnicity_groupsEastAsian + # ) / 5 = 0", + # + # "SouthAsian" = "( + # ethnicity_groupsEuropean + # + ethnicity_groupsHispanicDLatinAmerican + # + ethnicity_groupsEastAsian + # + `ethnicity_groupsJapanese` + # - 5 * ethnicity_groupsSouthAsian + # ) / 5 = 0", + # + # "HispanicDLatinAmerican" = "( + # ethnicity_groupsEuropean + # + ethnicity_groupsEastAsian + # + ethnicity_groupsSouthAsian + # + `ethnicity_groupsJapanese` + # - 5 * ethnicity_groupsHispanicDLatinAmerican + # ) / 5 = 0", + # + # "Japanese" = "( + # ethnicity_groupsEuropean + # + ethnicity_groupsHispanicDLatinAmerican + # + ethnicity_groupsSouthAsian + # + ethnicity_groupsEastAsian + # - 5 * `ethnicity_groupsJapanese` + # ) / 5 = 0" + # ), + + # Median instead and mad of mean and sd + robust=TRUE, + alpha = 0.1 + ) )) |> - mutate( + mutate( summary_tissue = map( brms_fit, function(x) { @@ -754,28 +831,28 @@ get_adjusted_matrix = function(effect_removed_df, column_adjusted){ ) %>% mutate(Rhat_ethnicity = map_dbl(brms_fit, - ~ summary(.x)$fixed |> - as_tibble(rownames = "par") |> - filter(par |> str_detect("ethnicity")) |> - pull(Rhat) |> - max() + ~ summary(.x)$fixed |> + as_tibble(rownames = "par") |> + filter(par |> str_detect("ethnicity")) |> + pull(Rhat) |> + max() )) |> mutate(Rhat_tissue = map_dbl(brms_fit, - ~ summary(.x)$random$tissue_groups |> - as_tibble() |> - pull(Rhat) |> - max() + ~ summary(.x)$random$tissue_groups |> + as_tibble() |> + pull(Rhat) |> + max() )) %>% - + select(-brms_fit), - + pattern = map(estimates_chunk), packages = c( "brms", "glue", "dplyr", "purrr", "rstan", "magrittr", "stringr"), resources = tar_resources(crew = tar_resources_crew("elastic")), error = "null" ), - + ## effect_removed ----- # This target generates adjusted model estimates by removing unwanted effects from the fitted Bayesian models, # thereby isolating the effects of interest. Here, nuisance covariates are set to NA and removed from the predictions. @@ -867,26 +944,26 @@ get_adjusted_matrix = function(effect_removed_df, column_adjusted){ resources = tar_resources(crew = tar_resources_crew("elastic")), error = "null" ), - - ## param ----- - tar_target( - param, - estimates_chunk %>% - mutate( - param = map( - brms_fit, - ~ summary(.x$fit) |> as.data.frame() - ) - )%>% - select(-brms_fit), - pattern = map(estimates_chunk), - packages = c( "brms", "dplyr", "purrr", "rstan"), - resources = tar_resources(crew = tar_resources_crew("elastic")), - error = "null" - ), - - # adjusted_matrix ----- - tar_target( + + ## param ----- + tar_target( + param, + estimates_chunk %>% + mutate( + param = map( + brms_fit, + ~ summary(.x$fit) |> as.data.frame() + ) + )%>% + select(-brms_fit), + pattern = map(estimates_chunk), + packages = c( "brms", "dplyr", "purrr", "rstan"), + resources = tar_resources(crew = tar_resources_crew("elastic")), + error = "null" + ), + + # adjusted_matrix ----- + tar_target( adjusted_assay_ethnicity, get_adjusted_matrix(effect_removed, brms_fit_adjusted_ethnicity_estimate), packages = c( "brms", "glue", "dplyr", "purrr", "rstan", "magrittr", "stringr", "tidySummarizedExperiment") , From e5981d71358e76f77a251a6cc7d9cda84bb47e9e Mon Sep 17 00:00:00 2001 From: Chen Zhan Date: Wed, 16 Apr 2025 14:20:22 +0800 Subject: [PATCH 03/12] create branch for pawsey_pipeline --- .../de_brms_whole_immune_system.R | 192 ++++++++++-------- 1 file changed, 104 insertions(+), 88 deletions(-) diff --git a/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R b/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R index 79778a9..4e367af 100644 --- a/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R +++ b/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R @@ -13,9 +13,15 @@ tar_script({ library(crew.cluster) # Set file path ----- - hdf5_path = "/hpcfs/groups/phoenix-hpc-mangiola_laboratory/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseudobulk_sample_is_immune" - metadata_path = "/hpcfs/groups/phoenix-hpc-mangiola_laboratory/Mangiola_ImmuneAtlas/taskforce_shared_folder/cell_metadata_1_0_6_sccomp_input_counts.rds" + ## Phoenix HPC setting ----- + # hdf5_path = "/hpcfs/groups/phoenix-hpc-mangiola_laboratory/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseudobulk_sample_is_immune" + # metadata_path = "/hpcfs/groups/phoenix-hpc-mangiola_laboratory/Mangiola_ImmuneAtlas/taskforce_shared_folder/cell_metadata_1_0_6_sccomp_input_counts.rds" + ## Pawsey setting ----- + hdf5_path = "/scratch/pawsey1192/zhanchen/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseudobulk_sample_is_immune" + metadata_path = "/scratch/pawsey1192/zhanchen/Mangiola_ImmuneAtlas/taskforce_shared_folder/cell_metadata_1_0_6_sccomp_input_counts.rds" + + tar_option_set( @@ -41,10 +47,12 @@ tar_script({ seconds_idle = 30, crashes_max = 7, options_cluster = crew_options_slurm( - script_lines = '#SBATCH -A saigencir003', + # script_lines = '#SBATCH -A saigencir003', + script_lines = '#SBATCH --account=pawsey1192 \n#SBATCH --time=1-00:00:00 \nsource /software/projects/pawsey1192/zhanchen/miniconda3/bin/activate R_443', memory_gigabytes_required = c(5, 10, 20, 40, 80, 160), cpus_per_task = 2, - time_minutes = c(60*4, 60*4, 60*4, 60*4, 60*24, 60*24), + partition = 'work', # for pawsey + # time_minutes = c(60*4, 60*4, 60*4, 60*4, 60*24, 60*24), verbose = T ) ), @@ -55,10 +63,12 @@ tar_script({ seconds_idle = 30, crashes_max = 5, options_cluster = crew_options_slurm( - script_lines = '#SBATCH -A saigencir003', + # script_lines = '#SBATCH -A saigencir003', + script_lines = '#SBATCH --account=pawsey1192 \n#SBATCH --time=1-00:00:00 \nsource /software/projects/pawsey1192/zhanchen/miniconda3/bin/activate R_443', memory_gigabytes_required = c(80, 160), cpus_per_task = 2, - time_minutes = c(60*24, 60*24), + partition = 'work', # for pawsey + # time_minutes = c(60*24, 60*24), verbose = T ) ), @@ -69,9 +79,11 @@ tar_script({ seconds_idle = 30, crashes_max = 5, options_cluster = crew_options_slurm( - script_lines = '#SBATCH -A saigencir003', + # script_lines = '#SBATCH -A saigencir003', + script_lines = '#SBATCH --account=pawsey1192 \n#SBATCH --time=1-00:00:00 \nsource /software/projects/pawsey1192/zhanchen/miniconda3/bin/activate R_443', memory_gigabytes_required = c(160), cpus_per_task = 30, + partition = 'work', # for pawsey verbose = T ) ) @@ -304,38 +316,38 @@ tar_script({ # se = # loadHDF5SummarizedExperiment(hdf5_path) |> # filter(is_gene_shared) |> - # + # #---------------------------------# # # Edit or add more filters here for analyses # #---------------------------------# # filter(is_immune & do_analyse) - # + # # TEMPORARY BECAUSE I FORGOT TO INTEGRATE AGE BINS # se = se |> # left_join( # readRDS(metadata_path) |> # distinct(sample_id, age_days, age_bin) # ) - # + # # Filter common genes # se = se[((assay(se, "gene_presence") > 0) |> rowSums() > (ncol(se) * 0.95)),,drop=FALSE ] - # + # # Filter samples that have enough genes > 0 but not too many # samples_with_right_number_of_detected_genes = # (se |> assay() > 0) |> # colSums() |> # divide_by(nrow(se)) |> # dplyr::between(0.3, 1) - # + # se = se[,samples_with_right_number_of_detected_genes] - # + # # Compute mean library size # mean_library_size <- se |> # assay("counts") |> # _[nrow(se) |> seq_len() |> sample(size = 2000), ] |> # colSums() |> # mean() - # + # # Optional: retrieve the sample name (column name in the SummarizedExperiment) # reference_sample <- colnames(se)[ # se |> @@ -344,16 +356,16 @@ tar_script({ # {\(x) abs(x - mean_library_size)}() |> # Calculate absolute difference from the mean # which.min() # Identify the smallest difference # ] - # + # message('TAR: pseudobulk_sample Phase2') # se = # se |> # keep_abundant(design = # se |> - # + # # Discretise the age for the following operation # mutate(is_old_individual = age_days > 50*365) |> - # + # # This is to resolve some confounders to preserve the genes. # # In this case we care about data variability, not the actual meaning of the variables # resolve_complete_confounders_of_non_interest(tissue_groups, sex, ethnicity_groups, is_old_individual) |> @@ -362,10 +374,10 @@ tar_script({ # model.matrix(~ tissue_groups + sex___altered + ethnicity_groups___altered + is_old_individual___altered, data = _ ), # minimum_counts = 100 # ) |> - # + # # Get scaling factor # scale_abundance(method = "TMMwsp", reference_sample = reference_sample) |> - # + # # Drop sex unknown as causes problem during fit # mutate( # sex = if_else(sex |> is.na(), "unknown", sex), @@ -373,13 +385,13 @@ tar_script({ # ) |> # filter(sex != "unknown") |> # filter(!age_bin |> is.na()) |> - # + # # Eliminate complete confounders # tidybulk:::resolve_complete_confounders_of_non_interest(assay_groups, dataset_id, disease_groups) |> - # + # # sibrary size factor is the reciproque of the multiplier (correction factor) # mutate(offset = log(1/multiplier)) |> - # + # # Set intercept # mutate( # ethnicity_groups = fct_relevel(ethnicity_groups, "European"), @@ -387,19 +399,20 @@ tar_script({ # disease_groups___altered = fct_relevel(disease_groups___altered, "Normal"), # age_bin = fct_relevel(age_bin, "Adolescence") # ) - # + # # # Add dispersion # # rowData(se) = # # rowData(se) |> # # as_tibble(rownames = ".feature") |> # # left_join(glmGamPoi_overdispersions |> enframe(name = ".feature", value = "dispersion")) |> # # data.frame(row.names = ".feature") |> DataFrame() - # + # message('TAR: pseudobulk_sample COMPLETE') # se # load process data to save time when testing - loadHDF5SummarizedExperiment('/hpcfs/groups/phoenix-hpc-mangiola_laboratory/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseduobulk_sample_tar_load_altered/') + # loadHDF5SummarizedExperiment('/hpcfs/groups/phoenix-hpc-mangiola_laboratory/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseduobulk_sample_tar_load_altered/') + loadHDF5SummarizedExperiment('/scratch/pawsey1192/zhanchen/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseduobulk_sample_tar_load_altered/') }, packages = c("tidybulk", "HDF5Array", "tidySummarizedExperiment", "magrittr", "tibble", "forcats"), resources = tar_resources(crew = tar_resources_crew("elastic_big_30_cores")), @@ -425,7 +438,7 @@ tar_script({ distinct(.feature)|> # testing genes that ran for long time # filter(.feature %in% readRDS('/hpcfs/groups/phoenix-hpc-mangiola_laboratory/Mangiola_ImmuneAtlas/ning_data/ethnicity_umap_selected_genes.rds')) |> - # head(1) %>% + slice_sample(n=1500) %>% group_by(.feature) |> tar_group(), iteration = "group", @@ -520,16 +533,16 @@ tar_script({ ) # prior Version 0 - prior = c( - prior(normal(i, 5), class = Intercept), - prior(normal(0, 2), class = Intercept, dpar = shape), - prior(normal(0, 5), class = b), - prior(normal(0, 2), class = b, dpar = shape) - # prior(beta(0.5381488, 10.3577433), class = "zi", lb = 0, ub = 1) # addition zi from V2 - ) |> - substitute(env = list(i = mean(log1p(data$counts / exp(data$offset))))) |> - eval() - # + # prior = c( + # prior(normal(i, 5), class = Intercept), + # prior(normal(0, 2), class = Intercept, dpar = shape), + # prior(normal(0, 5), class = b), + # prior(normal(0, 2), class = b, dpar = shape) + # # prior(beta(0.5381488, 10.3577433), class = "zi", lb = 0, ub = 1) # addition zi from V2 + # ) |> + # substitute(env = list(i = mean(log1p(data$counts / exp(data$offset))))) |> + # eval() + # chains = 2 # inits <- list(Intercept = mean(log1p(data$counts / exp(data$offset)))) # inits <- replicate(chains, inits, simplify = FALSE) @@ -612,68 +625,71 @@ tar_script({ # dynamically extract param from stan data # code used from brm - bterms <- brmsterms( - formula = brms:::validate_formula( - formula, data = data, family = zero_inflated_negbinomial(), - autocor = NULL, sparse = NULL, cov_ranef = NULL - ) - ) - bframe <- brms:::brmsframe(bterms, data) - sdata <- brms:::.standata( - bframe, data = data, prior = prior, - data2 = NULL, stanvars = NULL, threads = NULL - ) + # bterms <- brmsterms( + # formula = brms:::validate_formula( + # formula, data = data, family = zero_inflated_negbinomial(), + # autocor = NULL, sparse = NULL, cov_ranef = NULL + # ) + # ) + # bframe <- brms:::brmsframe(bterms, data) + # sdata <- brms:::.standata( + # bframe, data = data, prior = prior, + # data2 = NULL, stanvars = NULL, threads = NULL + # ) - Kc <- sdata$Kc - Kc_shape <- sdata$Kc_shape - M_1 <- sdata$M_1; N_1 <- sdata$N_1 - M_2 <- sdata$M_2; N_2 <- sdata$N_2 - M_3 <- sdata$M_3; N_3 <- sdata$N_3 + # Kc <- sdata$Kc + # Kc_shape <- sdata$Kc_shape + # M_1 <- sdata$M_1; N_1 <- sdata$N_1 + # M_2 <- sdata$M_2; N_2 <- sdata$N_2 + # M_3 <- sdata$M_3; N_3 <- sdata$N_3 - inits <- lapply(1:chains, function(i) { - list( + # inits <- lapply(1:chains, function(i) { + # list( - #### revert v0 prior - # Fixed effects for count part - b = rnorm(Kc, 0, 5), - # dynamically set mu for intercept - Intercept = rnorm(1, mean(log1p(data$counts / exp(data$offset))), 5), + # #### revert v0 prior + # # Fixed effects for count part + # b = rnorm(Kc, 0, 5), + # # dynamically set mu for intercept + # Intercept = rnorm(1, mean(log1p(data$counts / exp(data$offset))), 5), - # Fixed effects for shape submodel - b_shape = rnorm(Kc_shape, 0, 2), - Intercept_shape = rnorm(1, 0, 2) + # # Fixed effects for shape submodel + # b_shape = rnorm(Kc_shape, 0, 2), + # Intercept_shape = rnorm(1, 0, 2) - # # Fixed effects for count part - # b = 0.06161134 + 0.9263627 * rt(Kc, 3), - # # dynamically set mu for intercept - # Intercept = rnorm(1, mean(log1p(data$counts / exp(data$offset))), 2.438534), - # - # # Fixed effects for shape submodel - # b_shape = 0.07603337 + 0.8252114 * rt(Kc_shape, 3), - # Intercept_shape = rnorm(1, 0.4260793, 1.470536) + # # # Fixed effects for count part + # # b = 0.06161134 + 0.9263627 * rt(Kc, 3), + # # # dynamically set mu for intercept + # # Intercept = rnorm(1, mean(log1p(data$counts / exp(data$offset))), 2.438534), + # # + # # # Fixed effects for shape submodel + # # b_shape = 0.07603337 + 0.8252114 * rt(Kc_shape, 3), + # # Intercept_shape = rnorm(1, 0.4260793, 1.470536) - # Zero-inflation probability - # zi = rbeta(1, 0.5381488, 10.3577433) - # - # # Group-level standard deviations and effects - # sd_1 = abs(0.5703259 + 0.4147664 * rt(M_1, 3)), # count - # z_1 = replicate(M_1, rnorm(N_1, mean = 0 , sd = 0.08547970), simplify = FALSE), - # - # sd_2 = abs(0.5703259 + 0.4147664 * rt(M_2, 3)), # zi - # z_2 = matrix(rnorm(M_2 * N_2, mean = 0 , sd = 0.08547970), nrow = M_2, ncol = N_2), - # L_2 = diag(M_2), # no correlation (identity) - # - # sd_3 = abs(rnorm(M_3, 0.8670409, 0.1779553)), # shape - # z_3 = replicate(M_3, rnorm(N_3, mean = 0 , sd = 0.08547970), simplify = FALSE) - ) - }) + # # Zero-inflation probability + # # zi = rbeta(1, 0.5381488, 10.3577433) + # # + # # # Group-level standard deviations and effects + # # sd_1 = abs(0.5703259 + 0.4147664 * rt(M_1, 3)), # count + # # z_1 = replicate(M_1, rnorm(N_1, mean = 0 , sd = 0.08547970), simplify = FALSE), + # # + # # sd_2 = abs(0.5703259 + 0.4147664 * rt(M_2, 3)), # zi + # # z_2 = matrix(rnorm(M_2 * N_2, mean = 0 , sd = 0.08547970), nrow = M_2, ncol = N_2), + # # L_2 = diag(M_2), # no correlation (identity) + # # + # # sd_3 = abs(rnorm(M_3, 0.8670409, 0.1779553)), # shape + # # z_3 = replicate(M_3, rnorm(N_3, mean = 0 , sd = 0.08547970), simplify = FALSE) + # ) + # }) + + # script to set path of cmdstan on pawsey + # cmdstanr::set_cmdstan_path("/scratch/pawsey1192/zhanchen/.cmdstan/cmdstan-2.36.0") brm( formula = formula, data = data, family = zero_inflated_negbinomial(), - prior = prior, + # prior = prior, chains = chains, cores = pmax(as.numeric(parallelly::availableCores()), 2), #, threads = 2, warmup = 300, @@ -682,8 +698,8 @@ tar_script({ #sparse = TRUE, #save_model = glue("{external_directory}~/temp.rds"), #algorithm = "pathfinder", - init = inits, - iter = 400, # Increase iterations for better convergence + # init = inits, + iter = 800, # Increase iterations for better convergence sample_prior = TRUE ) From 2f1f14b10e2769e78f5e4aed5f8c1f89663bca0c Mon Sep 17 00:00:00 2001 From: Chen Zhan Date: Wed, 4 Jun 2025 13:59:43 +0800 Subject: [PATCH 04/12] pipeline code updated: new prior, within chain paralle --- .../de_brms_whole_immune_system.R | 1191 +++++++++++++---- 1 file changed, 950 insertions(+), 241 deletions(-) diff --git a/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R b/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R index 4e367af..e5eeee6 100644 --- a/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R +++ b/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R @@ -1,5 +1,22 @@ library(targets) +# cell type: +# [1] "b naive" "cd14 mono" "cd16 mono" "cd4 naive" +# [5] "cd4 tcm" "cd4 th17 em" "cd8 naive" "cd8 tcm" +# [9] "cd8 tem" "t cd4" "t cd8" "tgd" +# [13] "treg" "monocytic" "b memory" "cd4 th1/th17 em" +# [17] "cytotoxic" "cd4 th2 em" "progenitor" "mait" +# [21] "nk" "t" "macrophage" "cd4 fh em" +# [25] "cd4 tem" "cd4 th1 em" "cdc" "b" +# [29] "dc" "plasma" "granulocyte" "erythrocyte" +# [33] "pdc" "ilc" "neuron" "glial" +# [37] "pericyte" "endothelial" "immune" "blood" +# [41] "muscle" "stromal" "mesothelial" "epithelial" +# [45] "liver" "mast" "nkt" "renal" +# [49] "endocrine" "reproductive" "secretory" "fat" +# [53] "pneumocyte" "myoepithelial" "sensory" "lens" +# [57] "epidermal" "cartilage" "bone" + # SET Script ------ tar_script({ @@ -18,9 +35,8 @@ tar_script({ # metadata_path = "/hpcfs/groups/phoenix-hpc-mangiola_laboratory/Mangiola_ImmuneAtlas/taskforce_shared_folder/cell_metadata_1_0_6_sccomp_input_counts.rds" ## Pawsey setting ----- - hdf5_path = "/scratch/pawsey1192/zhanchen/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseudobulk_sample_is_immune" - metadata_path = "/scratch/pawsey1192/zhanchen/Mangiola_ImmuneAtlas/taskforce_shared_folder/cell_metadata_1_0_6_sccomp_input_counts.rds" - + hdf5_path = "/scratch/pawsey1192/zhanchen/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseudobulk_sample_cell_type" + target_cell_type = "cd4 fh em" tar_option_set( @@ -50,7 +66,7 @@ tar_script({ # script_lines = '#SBATCH -A saigencir003', script_lines = '#SBATCH --account=pawsey1192 \n#SBATCH --time=1-00:00:00 \nsource /software/projects/pawsey1192/zhanchen/miniconda3/bin/activate R_443', memory_gigabytes_required = c(5, 10, 20, 40, 80, 160), - cpus_per_task = 2, + cpus_per_task =16, partition = 'work', # for pawsey # time_minutes = c(60*4, 60*4, 60*4, 60*4, 60*24, 60*24), verbose = T @@ -281,6 +297,518 @@ tar_script({ return(m) } + + edit_covariates = function(tbl){ + + tissue_grouped = 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", + + # Epithelium + "epithelium of esophagus" + ), + + # Stomach + "stomach" = c( + "stomach", "body of stomach", "cardia of stomach" + ), + + # Small Intestine + "small intestine" = c( + "small intestine", "duodenum", "jejunum", "ileum", + + # Epithelium + "epithelium of small intestine", "jejunal epithelium", "ileal epithelium", + "submucosa of ileum", "lamina propria of small intestine" + ), + + # 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", + + # epithelium + "colonic epithelium", "submucosa of ascending colon", "lamina propria of large intestine", + "mucosa of colon", "lamina propria of mucosa of colon", "caecum epithelium" + ), + + # Digestive System (General) + "digestive system (general)" = c( + "intestine", "hindgut", "lamina propria", "mucosa" + ), + + # 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" + ), + + # 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) |> + rename(tissue = value) |> + mutate() + + ethnicity_grouped <- tribble( + ~self_reported_ethnicity, ~ethnicity_groups, + "unknown", "Other/Unknown", + "European", "European", + "Korean", "East Asian", + "Asian", "East Asian", + "Japanese", "Japanese", + "African American", "African", + "Hispanic or Latin American", "Hispanic/Latin American", + "Singaporean Chinese", "East Asian", + "Han Chinese", "East Asian", + "Singaporean Indian", "South Asian", + "Singaporean Malay", "Other/Unknown", + "British", "European", + "African", "African", + "South Asian", "South Asian", + "European American", "European", + "East Asian", "East Asian", + "American", "Other/Unknown", + "African American or Afro-Caribbean", "African", + "Oceanian", "Native American & Pacific Islander", + "Jewish Israeli", "Middle Eastern & North African", + "Chinese", "East Asian", + "South East Asian", "Other/Unknown", + "Greater Middle Eastern (Middle Eastern or North African or Persian)", "Middle Eastern & North African", + "Native American", "Native American & Pacific Islander", + "Pacific Islander", "Native American & Pacific Islander", + "Finnish", "European", + "Bangladeshi", "South Asian", + "Native American,Hispanic or Latin American", "Hispanic/Latin American", + "Irish", "European", + "Iraqi", "Middle Eastern & North African", + "European,Asian", "European" + ) + + assay_data_grouped <- tribble( + ~assay, ~assay_groups, + "10x 3' v2", "10x Genomics 3", + "10x 3' v3", "10x Genomics 3", + "10x 5' v2", "10x Genomics 5", + "10x 5' v1", "10x Genomics 5", + "MARS-seq", "Plate based Technologies", + "10x 3' transcription profiling", "10x Genomics 3", + "10x 5' transcription profiling", "10x Genomics 5", + "Smart-seq2", "Smart seq", + "microwell-seq", "Microwell Technologies", + "TruDrop", "TruDrop", + "Drop-seq", "Drop based Technologies", + "Seq-Well S3", "Microwell Technologies", + "GEXSCOPE technology", "Other Technologies", + "Seq-Well", "Microwell Technologies", + "sci-RNA-seq", "Other Technologies", + "10x 3' v1", "10x Genomics 3", + "BD Rhapsody Whole Transcriptome Analysis", "Other Technologies", + "BD Rhapsody Targeted mRNA", "Other Technologies", + "CEL-seq2", "Plate based Technologies", + "SPLiT-seq", "Other Technologies", + "STRT-seq", "Plate based Technologies", + "inDrop", "Drop based Technologies", + "Smart-seq v4", "Smart seq", + "ScaleBio single cell RNA sequencing", "Other Technologies" + ) + + + disease_data_grouped <- tribble( + ~disease, ~disease_groups, + + # Normal control + "normal", "Normal", + + # Isolated Diseases + "COVID-19", "COVID-19 related", + "post-COVID-19 disorder", "COVID-19 related", + "long COVID-19", "COVID-19 related", + "glioblastoma", "Glioblastoma", + "lung adenocarcinoma", "Lung Adenocarcinoma", + "systemic lupus erythematosus", "Systemic Lupus Erythematosus", + + # Infectious and Immune-related Diseases (other than COVID-19) + "Crohn disease", "Infectious and Immune-related Diseases", + "Crohn ileitis", "Infectious and Immune-related Diseases", + "pneumonia", "Infectious and Immune-related Diseases", + "common variable immunodeficiency", "Infectious and Immune-related Diseases", + "toxoplasmosis", "Infectious and Immune-related Diseases", + "Plasmodium malariae malaria", "Infectious and Immune-related Diseases", + "type 1 diabetes mellitus", "Infectious and Immune-related Diseases", + "influenza", "Infectious and Immune-related Diseases", + "chronic rhinitis", "Infectious and Immune-related Diseases", + "periodontitis", "Infectious and Immune-related Diseases", + "localized scleroderma", "Infectious and Immune-related Diseases", + "lymphangioleiomyomatosis", "Infectious and Immune-related Diseases", + "listeriosis", "Infectious and Immune-related Diseases", + + # Cancer (other than isolated cancers) + "squamous cell lung carcinoma", "Cancer", + "small cell lung carcinoma", "Cancer", + "non-small cell lung carcinoma", "Cancer", + "breast carcinoma", "Cancer", + "breast cancer", "Cancer", + "luminal B breast carcinoma", "Cancer", + "luminal A breast carcinoma", "Cancer", + "triple-negative breast carcinoma", "Cancer", + "gastric cancer", "Cancer", + "colorectal cancer", "Cancer", + "colon sessile serrated adenoma/polyp", "Cancer", + "follicular lymphoma", "Cancer", + "B-cell acute lymphoblastic leukemia", "Cancer", + "B-cell non-Hodgkin lymphoma", "Cancer", + "acute myeloid leukemia", "Cancer", + "acute promyelocytic leukemia", "Cancer", + "plasma cell myeloma", "Cancer", + "clear cell renal carcinoma", "Cancer", + "nonpapillary renal cell carcinoma", "Cancer", + "basal cell carcinoma", "Cancer", + "colorectal neoplasm", "Cancer", + "adenocarcinoma", "Cancer", + "chromophobe renal cell carcinoma", "Cancer", + "neuroendocrine carcinoma", "Cancer", + "lung large cell carcinoma", "Cancer", + "tongue cancer", "Cancer", + "Wilms tumor", "Cancer", + "pleomorphic carcinoma", "Cancer", + "blastoma", "Cancer", + + # Neurodegenerative and Neurological Disorders + "dementia", "Neurodegenerative and Neurological Disorders", + "Alzheimer disease", "Neurodegenerative and Neurological Disorders", + "Parkinson disease", "Neurodegenerative and Neurological Disorders", + "amyotrophic lateral sclerosis", "Neurodegenerative and Neurological Disorders", + "multiple sclerosis", "Neurodegenerative and Neurological Disorders", + "Down syndrome", "Neurodegenerative and Neurological Disorders", + "trisomy 18", "Neurodegenerative and Neurological Disorders", + "frontotemporal dementia", "Neurodegenerative and Neurological Disorders", + "temporal lobe epilepsy", "Neurodegenerative and Neurological Disorders", + "Lewy body dementia", "Neurodegenerative and Neurological Disorders", + "amyotrophic lateral sclerosis 26 with or without frontotemporal dementia", "Neurodegenerative and Neurological Disorders", + + # Respiratory Conditions + "pulmonary fibrosis", "Respiratory Conditions", + "respiratory system disorder", "Respiratory Conditions", + "chronic obstructive pulmonary disease", "Respiratory Conditions", + "cystic fibrosis", "Respiratory Conditions", + "interstitial lung disease", "Respiratory Conditions", + "hypersensitivity pneumonitis", "Respiratory Conditions", + "non-specific interstitial pneumonia", "Respiratory Conditions", + "aspiration pneumonia", "Respiratory Conditions", + "pulmonary emphysema", "Respiratory Conditions", + "pulmonary sarcoidosis", "Respiratory Conditions", + + # Cardiovascular Diseases + "myocardial infarction", "Cardiovascular Diseases", + "acute myocardial infarction", "Cardiovascular Diseases", + "dilated cardiomyopathy", "Cardiovascular Diseases", + "heart failure", "Cardiovascular Diseases", + "arrhythmogenic right ventricular cardiomyopathy", "Cardiovascular Diseases", + "congenital heart disease", "Cardiovascular Diseases", + "non-compaction cardiomyopathy", "Cardiovascular Diseases", + "cardiomyopathy", "Cardiovascular Diseases", + "heart disorder", "Cardiovascular Diseases", + + # Metabolic and Other Disorders + "type 2 diabetes mellitus", "Metabolic and Other Disorders", + "chronic kidney disease", "Metabolic and Other Disorders", + "digestive system disorder", "Metabolic and Other Disorders", + "primary sclerosing cholangitis", "Metabolic and Other Disorders", + "gastritis", "Metabolic and Other Disorders", + "acute kidney failure", "Metabolic and Other Disorders", + "tubular adenoma", "Metabolic and Other Disorders", + "benign prostatic hyperplasia", "Metabolic and Other Disorders", + "opiate dependence", "Metabolic and Other Disorders", + "gingivitis", "Metabolic and Other Disorders", + "hyperplastic polyp", "Metabolic and Other Disorders", + "clonal hematopoiesis", "Metabolic and Other Disorders", + "epilepsy", "Metabolic and Other Disorders", + "age related macular degeneration 7", "Metabolic and Other Disorders", + "kidney benign neoplasm", "Metabolic and Other Disorders", + "malignant pancreatic neoplasm", "Metabolic and Other Disorders", + "cataract", "Metabolic and Other Disorders", + "macular degeneration", "Metabolic and Other Disorders", + "hydrosalpinx", "Metabolic and Other Disorders", + "tubulovillous adenoma", "Metabolic and Other Disorders", + "gastric intestinal metaplasia", "Metabolic and Other Disorders", + "Barrett esophagus", "Metabolic and Other Disorders", + + # Other Diseases + "injury", "Other Diseases", + "anencephaly", "Other Diseases", + "primary biliary cholangitis", "Other Diseases", + "keloid", "Other Diseases", + "kidney oncocytoma", "Other Diseases", + "respiratory failure", "Other Diseases", + "pilocytic astrocytoma", "Other Diseases" + ) + + disease_data_grouped = + disease_data_grouped |> + left_join( + readr::read_csv("/home/zhanchen/From_scratch/Mangiola_ImmuneAtlas/disease_data_grouped_further.csv") |> + rename(disease_groups_further = disease_groups) + ) |> + mutate(disease_groups = if_else(!disease_groups_further |> is.na(), disease_groups_further, disease_groups)) |> + select(disease, disease_groups) + + tbl |> + + # TISSUE + select(-any_of("tissue_groups")) |> + left_join(tissue_grouped, copy=TRUE) |> + + # TECH + left_join(assay_data_grouped, copy=TRUE) |> + + # DISEASE + left_join(disease_data_grouped, copy=TRUE) |> + + # make disease tissue specific, omit Normal + mutate(disease_groups = paste(disease_groups, tissue_groups, sep = "_")) |> + mutate(disease_groups = if_else(disease_groups |> str_detect("Normal_.+"), "Normal", disease_groups)) |> + + + # TEMPORARY. de-group pancreas and liver + mutate(tissue_groups = case_when( + + tissue %in% c("gallbladder") ~ "gallbladder", + tissue %in% c("pancreas", "exocrine pancreas") ~ "pancreas", + tissue %in% c("liver", "caudate lobe of liver", "hepatic cecum" ) ~ "liver", + TRUE ~ tissue_groups + )) |> + + # SEX edit + mutate(sex = if_else(sex |> is.na(), "unknown", sex)) |> + + # Age + filter(age_days > 365) |> + mutate(age_years = age_days / 365) |> + mutate(age_bin = dplyr::case_when( + age_years < 3 ~ "Infancy", + age_years < 12 ~ "Childhood", + age_years < 20 ~ "Adolescence", + age_years < 40 ~ "Young Adulthood", + age_years < 50 ~ "Middle Age", + age_years < 60 ~ "Senior_50", + age_years < 70 ~ "Senior_60", + age_years >= 70 ~ "Senior_70", + TRUE ~ NA_character_ + )) |> + + # left_join(age_bin_table, copy=TRUE) |> + + # ETHNICITY + left_join(ethnicity_grouped, copy=TRUE) |> + + dplyr::select( + sample_id, donor_id, dataset_id, title, collection_id, age_days, age_bin, sex, ethnicity_groups, + tissue_groups, tissue, assay_groups, cell_type_unified_ensemble, cell_type, + disease_groups + ) |> + as_tibble() |> + + # Center based on adolescence + mutate(age_days_scaled = age_days |> scale(center = 50*365) |> as.numeric()) + + } #-----------------------# # Pipeline @@ -304,50 +832,201 @@ tar_script({ # # ), + # reference_sample ------ + # calculate reference_sample for scalig gene counts + # this ensures calculation only done once + tar_target( + reference_sample, + { + metadata <- + get_metadata(cache_directory="/home/zhanchen/From_scratch/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseduobulk/") + + metadata <- metadata %>% + select( + sample_id, donor_id, dataset_id, title, collection_id, + age_days, sex, self_reported_ethnicity, + tissue, assay, cell_type_unified_ensemble, cell_type, + disease + ) %>% as_tibble() %>% distinct() %>% + edit_covariates + + #---------------------------------# + # Edit or add more filters here for analyses + #---------------------------------# + # filter(is_gene_shared) |> + # filter(is_immune & do_analyse) + + se <- + loadHDF5SummarizedExperiment(hdf5_path) %>% + filter(cell_type_unified_ensemble == target_cell_type) %>% + filter(do_analyse) %>% + filter(is_gene_shared) |> + select(-c( + age_days, sex, ethnicity_groups, tissue_groups, + assay_groups, disease_groups, age_days_scaled + )) + + se <- se %>% + left_join( + metadata %>% + distinct( + sample_id, + age_days, age_bin, age_days_scaled, + sex, ethnicity_groups, + tissue, tissue_groups, + assay_groups, + disease_groups + ), + by = 'sample_id', + copy = T + ) + + # TEMPORARY BECAUSE I FORGOT TO INTEGRATE AGE BINS + # se = se |> + # left_join( + # readRDS(metadata_path) |> + # filter(age_days > 365) |> + # mutate(age_years = age_days / 365) |> + # mutate(age_bin = dplyr::case_when( + # age_years < 3 ~ "Infancy", + # age_years < 12 ~ "Childhood", + # age_years < 20 ~ "Adolescence", + # age_years < 40 ~ "Young Adulthood", + # age_years < 50 ~ "Middle Age", + # age_years < 60 ~ "Senior_50", + # age_years < 70 ~ "Senior_60", + # age_years >= 70 ~ "Senior_70", + # TRUE ~ NA_character_ + # )) %>% + # distinct(sample_id, age_days, age_bin) + # ) + + # Filter common genes + se = se[((assay(se, "gene_presence") > 0) |> rowSums() > (ncol(se) * 0.95)),,drop=FALSE ] + + # Filter samples that have enough genes > 0 but not too many + samples_with_right_number_of_detected_genes = + (se |> assay() > 0) |> + colSums() |> + divide_by(nrow(se)) |> + dplyr::between(0.3, 1) + + se = se[,samples_with_right_number_of_detected_genes] + + # Compute mean library size + mean_library_size <- se |> + assay("counts") |> + _[nrow(se) |> seq_len() |> sample(size = 2000), ] |> + colSums() |> + mean() + + # Optional: retrieve the sample name (column name in the SummarizedExperiment) + reference_sample <- colnames(se)[ + se |> + assay("counts") |> + colSums() |> + {\(x) abs(x - mean_library_size)}() |> # Calculate absolute difference from the mean + which.min() # Identify the smallest difference + ] + + reference_sample + }, + packages = c("tidybulk", "HDF5Array", "tidySummarizedExperiment", "magrittr", "tibble", "forcats", "readr", 'stringr', "cellNexus"), + resources = tar_resources(crew = tar_resources_crew("elastic_big_30_cores")), + memory = "persistent", + error = "stop" + ), + # This target loads and processes the pseudobulk sample data. It imports a HDF5 SummarizedExperiment, # applies filters to retain shared genes, immune cells, and samples marked for analysis, integrates age metadata, # filters for common genes and samples with an appropriate number of detected genes, computes the mean library size, # selects a reference sample, and performs normalisation and scaling. + # pseudobulk_sample ------ tar_target( - # pseudobulk_sample ------ pseudobulk_sample, { - # message('TAR: pseudobulk_sample START') - # se = - # loadHDF5SummarizedExperiment(hdf5_path) |> - # filter(is_gene_shared) |> - - # #---------------------------------# - # # Edit or add more filters here for analyses - # #---------------------------------# - # filter(is_immune & do_analyse) + metadata <- + get_metadata(cache_directory="/home/zhanchen/From_scratch/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseduobulk/") + + metadata <- metadata %>% + select( + sample_id, donor_id, dataset_id, title, collection_id, + age_days, sex, self_reported_ethnicity, + tissue, assay, cell_type_unified_ensemble, cell_type, + disease + ) %>% as_tibble() %>% distinct() %>% + edit_covariates - # # TEMPORARY BECAUSE I FORGOT TO INTEGRATE AGE BINS - # se = se |> + #---------------------------------# + # Edit or add more filters here for analyses + #---------------------------------# + # filter(is_gene_shared) |> + # filter(is_immune & do_analyse) + + se <- + loadHDF5SummarizedExperiment(hdf5_path) %>% + filter(cell_type_unified_ensemble == target_cell_type) %>% + filter(do_analyse) %>% + filter(is_gene_shared) |> + select(-c( + age_days, sex, ethnicity_groups, tissue_groups, + assay_groups, disease_groups, age_days_scaled + )) + + se <- se %>% + left_join( + metadata %>% + distinct( + sample_id, + age_days, age_bin, age_days_scaled, + sex, ethnicity_groups, + tissue, tissue_groups, + assay_groups, + disease_groups + ), + by = 'sample_id', + copy = T + ) + + # TEMPORARY BECAUSE I FORGOT TO INTEGRATE AGE BINS + # se = se |> # left_join( - # readRDS(metadata_path) |> - # distinct(sample_id, age_days, age_bin) + # readRDS(metadata_path) |> + # filter(age_days > 365) |> + # mutate(age_years = age_days / 365) |> + # mutate(age_bin = dplyr::case_when( + # age_years < 3 ~ "Infancy", + # age_years < 12 ~ "Childhood", + # age_years < 20 ~ "Adolescence", + # age_years < 40 ~ "Young Adulthood", + # age_years < 50 ~ "Middle Age", + # age_years < 60 ~ "Senior_50", + # age_years < 70 ~ "Senior_60", + # age_years >= 70 ~ "Senior_70", + # TRUE ~ NA_character_ + # )) %>% + # distinct(sample_id, age_days, age_bin) # ) - # # Filter common genes - # se = se[((assay(se, "gene_presence") > 0) |> rowSums() > (ncol(se) * 0.95)),,drop=FALSE ] + # Filter common genes + se = se[((assay(se, "gene_presence") > 0) |> rowSums() > (ncol(se) * 0.95)),,drop=FALSE ] - # # Filter samples that have enough genes > 0 but not too many - # samples_with_right_number_of_detected_genes = - # (se |> assay() > 0) |> - # colSums() |> - # divide_by(nrow(se)) |> - # dplyr::between(0.3, 1) + # Filter samples that have enough genes > 0 but not too many + samples_with_right_number_of_detected_genes = + (se |> assay() > 0) |> + colSums() |> + divide_by(nrow(se)) |> + dplyr::between(0.3, 1) - # se = se[,samples_with_right_number_of_detected_genes] + se = se[,samples_with_right_number_of_detected_genes] - # # Compute mean library size + # Compute mean library size # mean_library_size <- se |> # assay("counts") |> # _[nrow(se) |> seq_len() |> sample(size = 2000), ] |> # colSums() |> # mean() - + # # # Optional: retrieve the sample name (column name in the SummarizedExperiment) # reference_sample <- colnames(se)[ # se |> @@ -357,64 +1036,62 @@ tar_script({ # which.min() # Identify the smallest difference # ] - # message('TAR: pseudobulk_sample Phase2') - # se = - # se |> - # keep_abundant(design = - # se |> + se = + se |> + keep_abundant(design = + se |> - # # Discretise the age for the following operation - # mutate(is_old_individual = age_days > 50*365) |> + # Discretise the age for the following operation + mutate(is_old_individual = age_days > 50*365) |> - # # This is to resolve some confounders to preserve the genes. - # # In this case we care about data variability, not the actual meaning of the variables - # resolve_complete_confounders_of_non_interest(tissue_groups, sex, ethnicity_groups, is_old_individual) |> - # colData() |> - # droplevels() |> - # model.matrix(~ tissue_groups + sex___altered + ethnicity_groups___altered + is_old_individual___altered, data = _ ), - # minimum_counts = 100 - # ) |> + # This is to resolve some confounders to preserve the genes. + # In this case we care about data variability, not the actual meaning of the variables + resolve_complete_confounders_of_non_interest(tissue_groups, sex, ethnicity_groups, is_old_individual) |> + colData() |> + droplevels() |> + model.matrix(~ tissue_groups + sex___altered + ethnicity_groups___altered + is_old_individual___altered, data = _ ), + minimum_counts = 100 + ) |> - # # Get scaling factor - # scale_abundance(method = "TMMwsp", reference_sample = reference_sample) |> + # Get scaling factor + scale_abundance(method = "TMMwsp", reference_sample = reference_sample) |> - # # Drop sex unknown as causes problem during fit - # mutate( - # sex = if_else(sex |> is.na(), "unknown", sex), - # ethnicity_groups = if_else(ethnicity_groups |> is.na(), "Other/Unknown", ethnicity_groups) - # ) |> - # filter(sex != "unknown") |> - # filter(!age_bin |> is.na()) |> + # Drop sex unknown as causes problem during fit + mutate( + sex = if_else(sex |> is.na(), "unknown", sex), + ethnicity_groups = if_else(ethnicity_groups |> is.na(), "Other/Unknown", ethnicity_groups) + ) |> + filter(sex != "unknown") |> + filter(!age_bin |> is.na()) |> - # # Eliminate complete confounders - # tidybulk:::resolve_complete_confounders_of_non_interest(assay_groups, dataset_id, disease_groups) |> + # Eliminate complete confounders + tidybulk:::resolve_complete_confounders_of_non_interest(assay_groups, dataset_id, disease_groups) |> - # # sibrary size factor is the reciproque of the multiplier (correction factor) - # mutate(offset = log(1/multiplier)) |> + # library size factor is the reciproque of the multiplier (correction factor) + mutate(offset = log(1/multiplier)) |> - # # Set intercept - # mutate( - # ethnicity_groups = fct_relevel(ethnicity_groups, "European"), - # assay_groups___altered = fct_relevel(assay_groups___altered, "10x Genomics 3"), - # disease_groups___altered = fct_relevel(disease_groups___altered, "Normal"), - # age_bin = fct_relevel(age_bin, "Adolescence") - # ) - - # # # Add dispersion - # # rowData(se) = - # # rowData(se) |> - # # as_tibble(rownames = ".feature") |> - # # left_join(glmGamPoi_overdispersions |> enframe(name = ".feature", value = "dispersion")) |> - # # data.frame(row.names = ".feature") |> DataFrame() - - # message('TAR: pseudobulk_sample COMPLETE') - # se + # Set intercept + mutate( + ethnicity_groups = fct_relevel(ethnicity_groups, "European"), + assay_groups___altered = fct_relevel(assay_groups___altered, "10x Genomics 3"), + disease_groups___altered = fct_relevel(disease_groups___altered, "Normal"), + # age_bin = fct_relevel(age_bin, "Adolescence") + age_bin = fct_relevel(age_bin, "Senior_50") + ) + + # # Add dispersion + # rowData(se) = + # rowData(se) |> + # as_tibble(rownames = ".feature") |> + # left_join(glmGamPoi_overdispersions |> enframe(name = ".feature", value = "dispersion")) |> + # data.frame(row.names = ".feature") |> DataFrame() + + se # load process data to save time when testing # loadHDF5SummarizedExperiment('/hpcfs/groups/phoenix-hpc-mangiola_laboratory/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseduobulk_sample_tar_load_altered/') - loadHDF5SummarizedExperiment('/scratch/pawsey1192/zhanchen/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseduobulk_sample_tar_load_altered/') }, - packages = c("tidybulk", "HDF5Array", "tidySummarizedExperiment", "magrittr", "tibble", "forcats"), + packages = c("tidybulk", "HDF5Array", "tidySummarizedExperiment", "magrittr", "tibble", "forcats", "readr", 'stringr', "cellNexus"), resources = tar_resources(crew = tar_resources_crew("elastic_big_30_cores")), memory = "persistent", error = "stop" @@ -438,7 +1115,7 @@ tar_script({ distinct(.feature)|> # testing genes that ran for long time # filter(.feature %in% readRDS('/hpcfs/groups/phoenix-hpc-mangiola_laboratory/Mangiola_ImmuneAtlas/ning_data/ethnicity_umap_selected_genes.rds')) |> - slice_sample(n=1500) %>% + # slice_sample(n=1500) %>% group_by(.feature) |> tar_group(), iteration = "group", @@ -533,16 +1210,16 @@ tar_script({ ) # prior Version 0 - # prior = c( - # prior(normal(i, 5), class = Intercept), - # prior(normal(0, 2), class = Intercept, dpar = shape), - # prior(normal(0, 5), class = b), - # prior(normal(0, 2), class = b, dpar = shape) - # # prior(beta(0.5381488, 10.3577433), class = "zi", lb = 0, ub = 1) # addition zi from V2 - # ) |> - # substitute(env = list(i = mean(log1p(data$counts / exp(data$offset))))) |> - # eval() - + prior = c( + prior(student_t(3, i, 1.5), class = Intercept), + prior(student_t(3, 0, 1), class = Intercept, dpar = shape), + prior(student_t(3, 0, 5), class = b), + prior(student_t(3, 0, 2), class = b, dpar = shape) + # prior(beta(0.5381488, 10.3577433), class = "zi", lb = 0, ub = 1) # addition zi from V2 + ) |> + substitute(env = list(i = mean(log1p(data$counts / exp(data$offset))))) |> + eval() + # # chains = 2 # inits <- list(Intercept = mean(log1p(data$counts / exp(data$offset)))) # inits <- replicate(chains, inits, simplify = FALSE) @@ -625,82 +1302,80 @@ tar_script({ # dynamically extract param from stan data # code used from brm - # bterms <- brmsterms( - # formula = brms:::validate_formula( - # formula, data = data, family = zero_inflated_negbinomial(), - # autocor = NULL, sparse = NULL, cov_ranef = NULL - # ) - # ) - # bframe <- brms:::brmsframe(bterms, data) - # sdata <- brms:::.standata( - # bframe, data = data, prior = prior, - # data2 = NULL, stanvars = NULL, threads = NULL - # ) + bterms <- brmsterms( + formula = brms:::validate_formula( + formula, data = data, family = zero_inflated_negbinomial(), + autocor = NULL, sparse = NULL, cov_ranef = NULL + ) + ) + bframe <- brms:::brmsframe(bterms, data) + sdata <- brms:::.standata( + bframe, data = data, prior = prior, + data2 = NULL, stanvars = NULL, threads = NULL + ) - # Kc <- sdata$Kc - # Kc_shape <- sdata$Kc_shape - # M_1 <- sdata$M_1; N_1 <- sdata$N_1 - # M_2 <- sdata$M_2; N_2 <- sdata$N_2 - # M_3 <- sdata$M_3; N_3 <- sdata$N_3 + Kc <- sdata$Kc + Kc_shape <- sdata$Kc_shape + M_1 <- sdata$M_1; N_1 <- sdata$N_1 + M_2 <- sdata$M_2; N_2 <- sdata$N_2 + M_3 <- sdata$M_3; N_3 <- sdata$N_3 - # inits <- lapply(1:chains, function(i) { - # list( + inits <- lapply(1:chains, function(i) { + list( - # #### revert v0 prior - # # Fixed effects for count part - # b = rnorm(Kc, 0, 5), - # # dynamically set mu for intercept - # Intercept = rnorm(1, mean(log1p(data$counts / exp(data$offset))), 5), + #### revert v0 prior + # Fixed effects for count part + b = rnorm(Kc, 0, 5), + # dynamically set mu for intercept + Intercept = rnorm(1, mean(log1p(data$counts / exp(data$offset))), 1.5), - # # Fixed effects for shape submodel - # b_shape = rnorm(Kc_shape, 0, 2), - # Intercept_shape = rnorm(1, 0, 2) + # Fixed effects for shape submodel + b_shape = rnorm(Kc_shape, 0, 2), + Intercept_shape = rnorm(1, 0, 1) - # # # Fixed effects for count part - # # b = 0.06161134 + 0.9263627 * rt(Kc, 3), - # # # dynamically set mu for intercept - # # Intercept = rnorm(1, mean(log1p(data$counts / exp(data$offset))), 2.438534), - # # - # # # Fixed effects for shape submodel - # # b_shape = 0.07603337 + 0.8252114 * rt(Kc_shape, 3), - # # Intercept_shape = rnorm(1, 0.4260793, 1.470536) + # # Fixed effects for count part + # b = 0.06161134 + 0.9263627 * rt(Kc, 3), + # # dynamically set mu for intercept + # Intercept = rnorm(1, mean(log1p(data$counts / exp(data$offset))), 2.438534), + # + # # Fixed effects for shape submodel + # b_shape = 0.07603337 + 0.8252114 * rt(Kc_shape, 3), + # Intercept_shape = rnorm(1, 0.4260793, 1.470536) - # # Zero-inflation probability - # # zi = rbeta(1, 0.5381488, 10.3577433) - # # - # # # Group-level standard deviations and effects - # # sd_1 = abs(0.5703259 + 0.4147664 * rt(M_1, 3)), # count - # # z_1 = replicate(M_1, rnorm(N_1, mean = 0 , sd = 0.08547970), simplify = FALSE), - # # - # # sd_2 = abs(0.5703259 + 0.4147664 * rt(M_2, 3)), # zi - # # z_2 = matrix(rnorm(M_2 * N_2, mean = 0 , sd = 0.08547970), nrow = M_2, ncol = N_2), - # # L_2 = diag(M_2), # no correlation (identity) - # # - # # sd_3 = abs(rnorm(M_3, 0.8670409, 0.1779553)), # shape - # # z_3 = replicate(M_3, rnorm(N_3, mean = 0 , sd = 0.08547970), simplify = FALSE) - # ) - # }) - - # script to set path of cmdstan on pawsey - # cmdstanr::set_cmdstan_path("/scratch/pawsey1192/zhanchen/.cmdstan/cmdstan-2.36.0") + # Zero-inflation probability + # zi = rbeta(1, 0.5381488, 10.3577433) + # + # # Group-level standard deviations and effects + # sd_1 = abs(0.5703259 + 0.4147664 * rt(M_1, 3)), # count + # z_1 = replicate(M_1, rnorm(N_1, mean = 0 , sd = 0.08547970), simplify = FALSE), + # + # sd_2 = abs(0.5703259 + 0.4147664 * rt(M_2, 3)), # zi + # z_2 = matrix(rnorm(M_2 * N_2, mean = 0 , sd = 0.08547970), nrow = M_2, ncol = N_2), + # L_2 = diag(M_2), # no correlation (identity) + # + # sd_3 = abs(rnorm(M_3, 0.8670409, 0.1779553)), # shape + # z_3 = replicate(M_3, rnorm(N_3, mean = 0 , sd = 0.08547970), simplify = FALSE) + ) + }) brm( formula = formula, data = data, family = zero_inflated_negbinomial(), - # prior = prior, + prior = prior, chains = chains, - cores = pmax(as.numeric(parallelly::availableCores()), 2), #, threads = 2, - warmup = 300, + cores = pmin(as.numeric(parallelly::availableCores()), chains), + threads = threading(threads = (as.numeric(parallelly::availableCores()) / chains) |> floor()), + warmup = 400, refresh = 10, backend = "cmdstanr", #sparse = TRUE, #save_model = glue("{external_directory}~/temp.rds"), #algorithm = "pathfinder", - # init = inits, - iter = 800, # Increase iterations for better convergence - sample_prior = TRUE + # sample_prior = TRUE, + init = inits, + iter = 600 # Increase iterations for better convergence ) })) |> @@ -720,104 +1395,138 @@ tar_script({ tar_target( summary, estimates_chunk |> - mutate(summary_ethnicity = map(brms_fit, ~ .x |> hypothesis( - c( - "Europeans" = "(ethnicity_groupsAfrican - + ethnicity_groupsEastAsian - + ethnicity_groupsHispanicDLatinAmerican - + ethnicity_groupsSouthAsian - + `ethnicity_groupsJapanese`) / 5 = 0", - "EastAsian" = "( - ethnicity_groupsAfrican - + ethnicity_groupsHispanicDLatinAmerican - + ethnicity_groupsSouthAsian - + `ethnicity_groupsJapanese` - - 5 * ethnicity_groupsEastAsian - ) / 5 = 0", - "SouthAsian" = "( - ethnicity_groupsAfrican - + ethnicity_groupsHispanicDLatinAmerican - + ethnicity_groupsEastAsian - + `ethnicity_groupsJapanese` - - 5 * ethnicity_groupsSouthAsian - ) / 5 = 0", - "African" = "( - ethnicity_groupsEastAsian - + ethnicity_groupsHispanicDLatinAmerican - + ethnicity_groupsSouthAsian - + `ethnicity_groupsJapanese` - - 5 * ethnicity_groupsAfrican - ) / 5 = 0", - "HispanicDLatinAmerican" = "( - ethnicity_groupsAfrican - + ethnicity_groupsEastAsian - + ethnicity_groupsSouthAsian - + `ethnicity_groupsJapanese` - - 5 * ethnicity_groupsHispanicDLatinAmerican - ) / 5 = 0", + # mutate(summary_ethnicity = map(brms_fit, ~ .x |> hypothesis( + # c( + # "Europeans" = "(ethnicity_groupsAfrican + # + ethnicity_groupsEastAsian + # + ethnicity_groupsHispanicDLatinAmerican + # + ethnicity_groupsSouthAsian + # + `ethnicity_groupsJapanese`) / 5 = 0", + # "EastAsian" = "( + # ethnicity_groupsAfrican + # + ethnicity_groupsHispanicDLatinAmerican + # + ethnicity_groupsSouthAsian + # + `ethnicity_groupsJapanese` + # - 5 * ethnicity_groupsEastAsian + # ) / 5 = 0", + # "SouthAsian" = "( + # ethnicity_groupsAfrican + # + ethnicity_groupsHispanicDLatinAmerican + # + ethnicity_groupsEastAsian + # + `ethnicity_groupsJapanese` + # - 5 * ethnicity_groupsSouthAsian + # ) / 5 = 0", + # "African" = "( + # ethnicity_groupsEastAsian + # + ethnicity_groupsHispanicDLatinAmerican + # + ethnicity_groupsSouthAsian + # + `ethnicity_groupsJapanese` + # - 5 * ethnicity_groupsAfrican + # ) / 5 = 0", + # "HispanicDLatinAmerican" = "( + # ethnicity_groupsAfrican + # + ethnicity_groupsEastAsian + # + ethnicity_groupsSouthAsian + # + `ethnicity_groupsJapanese` + # - 5 * ethnicity_groupsHispanicDLatinAmerican + # ) / 5 = 0", - "Japanese" = "( - ethnicity_groupsAfrican - + ethnicity_groupsHispanicDLatinAmerican - + ethnicity_groupsSouthAsian - + ethnicity_groupsEastAsian - - 5 * `ethnicity_groupsJapanese` - ) / 5 = 0" - ), - # c( - # "African" = "(ethnicity_groupsEuropean - # + ethnicity_groupsEastAsian - # + ethnicity_groupsHispanicDLatinAmerican - # + ethnicity_groupsSouthAsian - # + `ethnicity_groupsJapanese`) / 5 = 0", - # - # "Europeans" = "( - # ethnicity_groupsEastAsian - # + ethnicity_groupsHispanicDLatinAmerican - # + ethnicity_groupsSouthAsian - # + `ethnicity_groupsJapanese` - # - 5 * ethnicity_groupsEuropean - # ) / 5 = 0", - # - # "EastAsian" = "( - # ethnicity_groupsEuropean - # + ethnicity_groupsHispanicDLatinAmerican - # + ethnicity_groupsSouthAsian - # + `ethnicity_groupsJapanese` - # - 5 * ethnicity_groupsEastAsian - # ) / 5 = 0", - # - # "SouthAsian" = "( - # ethnicity_groupsEuropean - # + ethnicity_groupsHispanicDLatinAmerican - # + ethnicity_groupsEastAsian - # + `ethnicity_groupsJapanese` - # - 5 * ethnicity_groupsSouthAsian - # ) / 5 = 0", - # - # "HispanicDLatinAmerican" = "( - # ethnicity_groupsEuropean - # + ethnicity_groupsEastAsian - # + ethnicity_groupsSouthAsian - # + `ethnicity_groupsJapanese` - # - 5 * ethnicity_groupsHispanicDLatinAmerican - # ) / 5 = 0", - # - # "Japanese" = "( - # ethnicity_groupsEuropean - # + ethnicity_groupsHispanicDLatinAmerican - # + ethnicity_groupsSouthAsian - # + ethnicity_groupsEastAsian - # - 5 * `ethnicity_groupsJapanese` - # ) / 5 = 0" - # ), + # "Japanese" = "( + # ethnicity_groupsAfrican + # + ethnicity_groupsHispanicDLatinAmerican + # + ethnicity_groupsSouthAsian + # + ethnicity_groupsEastAsian + # - 5 * `ethnicity_groupsJapanese` + # ) / 5 = 0" + # ), + # # c( + # # "African" = "(ethnicity_groupsEuropean + # # + ethnicity_groupsEastAsian + # # + ethnicity_groupsHispanicDLatinAmerican + # # + ethnicity_groupsSouthAsian + # # + `ethnicity_groupsJapanese`) / 5 = 0", + # # + # # "Europeans" = "( + # # ethnicity_groupsEastAsian + # # + ethnicity_groupsHispanicDLatinAmerican + # # + ethnicity_groupsSouthAsian + # # + `ethnicity_groupsJapanese` + # # - 5 * ethnicity_groupsEuropean + # # ) / 5 = 0", + # # + # # "EastAsian" = "( + # # ethnicity_groupsEuropean + # # + ethnicity_groupsHispanicDLatinAmerican + # # + ethnicity_groupsSouthAsian + # # + `ethnicity_groupsJapanese` + # # - 5 * ethnicity_groupsEastAsian + # # ) / 5 = 0", + # # + # # "SouthAsian" = "( + # # ethnicity_groupsEuropean + # # + ethnicity_groupsHispanicDLatinAmerican + # # + ethnicity_groupsEastAsian + # # + `ethnicity_groupsJapanese` + # # - 5 * ethnicity_groupsSouthAsian + # # ) / 5 = 0", + # # + # # "HispanicDLatinAmerican" = "( + # # ethnicity_groupsEuropean + # # + ethnicity_groupsEastAsian + # # + ethnicity_groupsSouthAsian + # # + `ethnicity_groupsJapanese` + # # - 5 * ethnicity_groupsHispanicDLatinAmerican + # # ) / 5 = 0", + # # + # # "Japanese" = "( + # # ethnicity_groupsEuropean + # # + ethnicity_groupsHispanicDLatinAmerican + # # + ethnicity_groupsSouthAsian + # # + ethnicity_groupsEastAsian + # # - 5 * `ethnicity_groupsJapanese` + # # ) / 5 = 0" + # # ), - # Median instead and mad of mean and sd - robust=TRUE, - alpha = 0.1 - ) - )) |> - + # # Median instead and mad of mean and sd + # robust=TRUE, + # alpha = 0.1 + # ) + # )) |> + + mutate( + + summary_ethnicity = map( + + brms_fit, function(x) { + + params = x$fit %>% summary() |> _[[1]] |> rownames() + params = params[grepl("^b_ethnicity_groups", params)] %>% sub("^b_", "", .) %>% setdiff(c("ethnicity_groupsOtherDUnknown", "ethnicity_groupsNativeAmerican&PacificIslander")) %>% paste0("`", . , "`") + ethnicity_groups_names <- sub("`ethnicity_groups(.*)`", "\\1", params) + + equations <- sapply(seq_along(params), function(i) { + this_ethnicity <- ethnicity_groups_names[i] + this_param <- params[i] + other_params <- params[-i] + avg_expr <- paste0("(", paste(other_params, collapse = " + "), ")/", length(other_params) + 1) + eq <- paste0(this_param, " - ", avg_expr, " = 0") + eq + }) + names(equations) <- ethnicity_groups_names + equations = append( + equations, + c('Europeans' = paste0("(", paste(params, collapse = " + "), ")/", length(params), ' = 0')) + ) + + return( + x |> hypothesis(equations, robust=TRUE, alpha = 0.1) + ) + + } + + ) + + ) %>% + mutate( summary_tissue = map( From ea510973ffebfa1a8a494a611d844a5323adccee Mon Sep 17 00:00:00 2001 From: Chen Zhan Date: Tue, 10 Jun 2025 10:49:38 +0800 Subject: [PATCH 05/12] code revision to fit pawsey --- rebuttal_CellPress/PAWSEY_run_sccomp_script.R | 101 ++ .../PAWSEY_sccomp_tar_scripts.R | 927 ++++++++++++++++++ 2 files changed, 1028 insertions(+) create mode 100644 rebuttal_CellPress/PAWSEY_run_sccomp_script.R create mode 100644 rebuttal_CellPress/PAWSEY_sccomp_tar_scripts.R diff --git a/rebuttal_CellPress/PAWSEY_run_sccomp_script.R b/rebuttal_CellPress/PAWSEY_run_sccomp_script.R new file mode 100644 index 0000000..fe5414f --- /dev/null +++ b/rebuttal_CellPress/PAWSEY_run_sccomp_script.R @@ -0,0 +1,101 @@ + +library(targets) +setwd('/home/zhanchen/From_scratch/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2') + +tar_make( + # callr_function = NULL, + script = "/scratch/pawsey1192/zhanchen/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2/_targets.R", + store = "/scratch/pawsey1192/zhanchen/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2/_targets", + reporter = "verbose" #, callr_function = NULL +) + + + + + + + +system("rclone copy /scratch/pawsey1192/zhanchen/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2/estimates_age_bins___L3___disease_TRUE___immune_only_TRUE.rds UofA_Box:/Mangiola_ImmuneAtlas/taskforce_shared_folder/sccomp_on_cellNexus_1_0_10_2/") + + + +tar_meta(store = "/scratch/pawsey1192/zhanchen/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2/_targets") |> + arrange(desc(time)) |> + filter(!error |> is.na()) |> + dplyr::select(name, error) + + + + +library(tidyverse) +library(sccomp) +library(magrittr) +library(glue) +library(forcats) +library(stringr) + +library(arrow) +library(dplyr) +library(duckdb) + +library(targets) + +x = tar_read(input_relative, store = "/vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/_targets") + + +tar_workspace(estimates_ae6e35523d730ab3, + script = "/vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/_targets.R", + store = "/vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/_targets" + ) + +tar_meta(starts_with("estimates_"), store = "/vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/_targets") + +# Age proportion prediciton +estimates_age_bins |> + sccomp_predict( + formula_composition = ~ 1 + age_bin*sex + (1 + age_bin*sex | tissue_groups), + number_of_draws = 100, + summary_instead_of_draws = TRUE + ) |> + mutate(age_bin = factor( + age_bin, + c("Infancy", "Childhood", "Adolescence", "Young Adulthood", "Middle Age", "Senior"), + ordered = TRUE + )) |> + mutate(age_bin_numeric = age_bin |> as.integer()) |> + saveRDS("/vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/prediction_age_bins.rds") + +system("~/bin/rclone copy /vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/prediction_age_bins.rds UofA_Box:/Mangiola_ImmuneAtlas/taskforce_shared_folder/") + +tar_read(formula_df, store = "/vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/_targets") + + +# For Hong +estimate_age_bins = readRDS("/vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/estimates_age_bins___L3.rds") +estimate_age_bins = estimate_age_bins |> dplyr::select(-count_data) +attr(estimate_age_bins, "fit") = NULL +estimate_age_bins |> saveRDS("/vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/estimates_age_bins_effect_tibble_only.rds") +system("~/bin/rclone copy /vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/estimates_age_bins_effect_tibble_only.rds UofA_Box:/immune_map_disease/") + +# Save fit +library(magrittr) +estimate_age_bins = readRDS("/vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/estimates_age_bins.rds") +estimate_age_bins |> attr("fit") %$% save_object(file = "/vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/estimates_age_bins_FIT_FOR_PORTABILITY.rds") +system("~/bin/rclone copy /vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/estimates_age_bins_FIT_FOR_PORTABILITY.rds UofA_Box:/Mangiola_ImmuneAtlas/taskforce_shared_folder/") +estimate_age_bins |> attr("fit") = readRDS("/vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/estimates_age_bins_FIT_FOR_PORTABILITY.rds") +estimate_age_bins |> saveRDS("/vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/estimates_age_bins.rds") +system("~/bin/rclone copy /vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/estimates_age_bins.rds UofA_Box:/Mangiola_ImmuneAtlas/taskforce_shared_folder/") + + +# estimate_age_bins |> saveRDS("/vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/estimates_age_bins.rds") +system("~/bin/rclone copy /vast/projects/mangiola_immune_map/PostDoc/immuneHealthyBodyMap/sccomp_on_cellNexus_1_0_10_2/21_11_2024_sccomp_archive_before_factor_ordering/estimates_age_bins.rds UofA_Box:/Mangiola_ImmuneAtlas/taskforce_shared_folder/") + +# Benchmark +tic() +estimate_age_bins |> sccomp_test(contrasts = c( "respiratory system" = "sexmale + `sexmale___respiratory system`", + "blood" = "sexmale + `sexmale___blood`")) +toc() + + +estimate_age_bins |> + sccomp_test() diff --git a/rebuttal_CellPress/PAWSEY_sccomp_tar_scripts.R b/rebuttal_CellPress/PAWSEY_sccomp_tar_scripts.R new file mode 100644 index 0000000..b15e9a7 --- /dev/null +++ b/rebuttal_CellPress/PAWSEY_sccomp_tar_scripts.R @@ -0,0 +1,927 @@ +library(targets) + +tar_script({ + + #-----------------------# + # Input + #-----------------------# + library(tidyverse) + library(targets) + library(tarchetypes) + library(glue) + library(qs) + library(crew) + library(crew.cluster) + + #-----------------------# + # Packages + #-----------------------# + tar_option_set( + + memory = "transient", + garbage_collection = 100, + storage = "worker", + retrieval = "worker", + workspace_on_error = TRUE, + # workspaces = "estimates", + # format = "qs", + + + #-----------------------# + # SLURM + #-----------------------# + controller = crew_controller_group( + + + + crew_controller_slurm( + name = "slurm_1_80", + tasks_max = 1, + workers = 40, + seconds_idle = 30, + crashes_max = 7, + options_cluster = crew_options_slurm( + script_lines = '#SBATCH --account=pawsey1192 \n#SBATCH --time=1-00:00:00 \nsource /software/projects/pawsey1192/zhanchen/miniconda3/bin/activate R_443', + memory_gigabytes_required = 80, + cpus_per_task = 30, + partition = 'work', # for pawsey + verbose = T + ) + ), + crew_controller_slurm( + name = "slurm_1_200", + tasks_max = 1, + workers = 20, + seconds_idle = 30, + crashes_max = 7, + options_cluster = crew_options_slurm( + script_lines = '#SBATCH --account=pawsey1192 \n#SBATCH --time=1-00:00:00 \nsource /software/projects/pawsey1192/zhanchen/miniconda3/bin/activate R_443', + memory_gigabytes_required = 200, + cpus_per_task = 2, + partition = 'work', # for pawsey + verbose = T + ) + ) + ), + debug = "input_relative", + + resources = tar_resources(crew = tar_resources_crew("slurm_1_80")) + #, # Set the target you want to debug. + # + ) + + #-----------------------# + # FUNCTIONS + #-----------------------# + + edit_covariates = function(tbl){ + + tissue_grouped = 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", + + # Epithelium + "epithelium of esophagus" + ), + + # Stomach + "stomach" = c( + "stomach", "body of stomach", "cardia of stomach" + ), + + # Small Intestine + "small intestine" = c( + "small intestine", "duodenum", "jejunum", "ileum", + + # Epithelium + "epithelium of small intestine", "jejunal epithelium", "ileal epithelium", + "submucosa of ileum", "lamina propria of small intestine" + ), + + # 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", + + # epithelium + "colonic epithelium", "submucosa of ascending colon", "lamina propria of large intestine", + "mucosa of colon", "lamina propria of mucosa of colon", "caecum epithelium" + ), + + # Digestive System (General) + "digestive system (general)" = c( + "intestine", "hindgut", "lamina propria", "mucosa" + ), + + # 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" + ), + + # 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) |> + rename(tissue = value) |> + mutate() + + ethnicity_grouped <- tribble( + ~self_reported_ethnicity, ~ethnicity_groups, + "unknown", "Other/Unknown", + "European", "European", + "Korean", "East Asian", + "Asian", "East Asian", + "Japanese", "Japanese", + "African American", "African", + "Hispanic or Latin American", "Hispanic/Latin American", + "Singaporean Chinese", "East Asian", + "Han Chinese", "East Asian", + "Singaporean Indian", "South Asian", + "Singaporean Malay", "Other/Unknown", + "British", "European", + "African", "African", + "South Asian", "South Asian", + "European American", "European", + "East Asian", "East Asian", + "American", "Other/Unknown", + "African American or Afro-Caribbean", "African", + "Oceanian", "Native American & Pacific Islander", + "Jewish Israeli", "Middle Eastern & North African", + "Chinese", "East Asian", + "South East Asian", "Other/Unknown", + "Greater Middle Eastern (Middle Eastern or North African or Persian)", "Middle Eastern & North African", + "Native American", "Native American & Pacific Islander", + "Pacific Islander", "Native American & Pacific Islander", + "Finnish", "European", + "Bangladeshi", "South Asian", + "Native American,Hispanic or Latin American", "Hispanic/Latin American", + "Irish", "European", + "Iraqi", "Middle Eastern & North African", + "European,Asian", "European" + ) + + assay_data_grouped <- tribble( + ~assay, ~assay_groups, + "10x 3' v2", "10x Genomics 3", + "10x 3' v3", "10x Genomics 3", + "10x 5' v2", "10x Genomics 5", + "10x 5' v1", "10x Genomics 5", + "MARS-seq", "Plate based Technologies", + "10x 3' transcription profiling", "10x Genomics 3", + "10x 5' transcription profiling", "10x Genomics 5", + "Smart-seq2", "Smart seq", + "microwell-seq", "Microwell Technologies", + "TruDrop", "TruDrop", + "Drop-seq", "Drop based Technologies", + "Seq-Well S3", "Microwell Technologies", + "GEXSCOPE technology", "Other Technologies", + "Seq-Well", "Microwell Technologies", + "sci-RNA-seq", "Other Technologies", + "10x 3' v1", "10x Genomics 3", + "BD Rhapsody Whole Transcriptome Analysis", "Other Technologies", + "BD Rhapsody Targeted mRNA", "Other Technologies", + "CEL-seq2", "Plate based Technologies", + "SPLiT-seq", "Other Technologies", + "STRT-seq", "Plate based Technologies", + "inDrop", "Drop based Technologies", + "Smart-seq v4", "Smart seq", + "ScaleBio single cell RNA sequencing", "Other Technologies" + ) + + + disease_data_grouped <- tribble( + ~disease, ~disease_groups, + + # Normal control + "normal", "Normal", + + # Isolated Diseases + "COVID-19", "COVID-19 related", + "post-COVID-19 disorder", "COVID-19 related", + "long COVID-19", "COVID-19 related", + "glioblastoma", "Glioblastoma", + "lung adenocarcinoma", "Lung Adenocarcinoma", + "systemic lupus erythematosus", "Systemic Lupus Erythematosus", + + # Infectious and Immune-related Diseases (other than COVID-19) + "Crohn disease", "Infectious and Immune-related Diseases", + "Crohn ileitis", "Infectious and Immune-related Diseases", + "pneumonia", "Infectious and Immune-related Diseases", + "common variable immunodeficiency", "Infectious and Immune-related Diseases", + "toxoplasmosis", "Infectious and Immune-related Diseases", + "Plasmodium malariae malaria", "Infectious and Immune-related Diseases", + "type 1 diabetes mellitus", "Infectious and Immune-related Diseases", + "influenza", "Infectious and Immune-related Diseases", + "chronic rhinitis", "Infectious and Immune-related Diseases", + "periodontitis", "Infectious and Immune-related Diseases", + "localized scleroderma", "Infectious and Immune-related Diseases", + "lymphangioleiomyomatosis", "Infectious and Immune-related Diseases", + "listeriosis", "Infectious and Immune-related Diseases", + + # Cancer (other than isolated cancers) + "squamous cell lung carcinoma", "Cancer", + "small cell lung carcinoma", "Cancer", + "non-small cell lung carcinoma", "Cancer", + "breast carcinoma", "Cancer", + "breast cancer", "Cancer", + "luminal B breast carcinoma", "Cancer", + "luminal A breast carcinoma", "Cancer", + "triple-negative breast carcinoma", "Cancer", + "gastric cancer", "Cancer", + "colorectal cancer", "Cancer", + "colon sessile serrated adenoma/polyp", "Cancer", + "follicular lymphoma", "Cancer", + "B-cell acute lymphoblastic leukemia", "Cancer", + "B-cell non-Hodgkin lymphoma", "Cancer", + "acute myeloid leukemia", "Cancer", + "acute promyelocytic leukemia", "Cancer", + "plasma cell myeloma", "Cancer", + "clear cell renal carcinoma", "Cancer", + "nonpapillary renal cell carcinoma", "Cancer", + "basal cell carcinoma", "Cancer", + "colorectal neoplasm", "Cancer", + "adenocarcinoma", "Cancer", + "chromophobe renal cell carcinoma", "Cancer", + "neuroendocrine carcinoma", "Cancer", + "lung large cell carcinoma", "Cancer", + "tongue cancer", "Cancer", + "Wilms tumor", "Cancer", + "pleomorphic carcinoma", "Cancer", + "blastoma", "Cancer", + + # Neurodegenerative and Neurological Disorders + "dementia", "Neurodegenerative and Neurological Disorders", + "Alzheimer disease", "Neurodegenerative and Neurological Disorders", + "Parkinson disease", "Neurodegenerative and Neurological Disorders", + "amyotrophic lateral sclerosis", "Neurodegenerative and Neurological Disorders", + "multiple sclerosis", "Neurodegenerative and Neurological Disorders", + "Down syndrome", "Neurodegenerative and Neurological Disorders", + "trisomy 18", "Neurodegenerative and Neurological Disorders", + "frontotemporal dementia", "Neurodegenerative and Neurological Disorders", + "temporal lobe epilepsy", "Neurodegenerative and Neurological Disorders", + "Lewy body dementia", "Neurodegenerative and Neurological Disorders", + "amyotrophic lateral sclerosis 26 with or without frontotemporal dementia", "Neurodegenerative and Neurological Disorders", + + # Respiratory Conditions + "pulmonary fibrosis", "Respiratory Conditions", + "respiratory system disorder", "Respiratory Conditions", + "chronic obstructive pulmonary disease", "Respiratory Conditions", + "cystic fibrosis", "Respiratory Conditions", + "interstitial lung disease", "Respiratory Conditions", + "hypersensitivity pneumonitis", "Respiratory Conditions", + "non-specific interstitial pneumonia", "Respiratory Conditions", + "aspiration pneumonia", "Respiratory Conditions", + "pulmonary emphysema", "Respiratory Conditions", + "pulmonary sarcoidosis", "Respiratory Conditions", + + # Cardiovascular Diseases + "myocardial infarction", "Cardiovascular Diseases", + "acute myocardial infarction", "Cardiovascular Diseases", + "dilated cardiomyopathy", "Cardiovascular Diseases", + "heart failure", "Cardiovascular Diseases", + "arrhythmogenic right ventricular cardiomyopathy", "Cardiovascular Diseases", + "congenital heart disease", "Cardiovascular Diseases", + "non-compaction cardiomyopathy", "Cardiovascular Diseases", + "cardiomyopathy", "Cardiovascular Diseases", + "heart disorder", "Cardiovascular Diseases", + + # Metabolic and Other Disorders + "type 2 diabetes mellitus", "Metabolic and Other Disorders", + "chronic kidney disease", "Metabolic and Other Disorders", + "digestive system disorder", "Metabolic and Other Disorders", + "primary sclerosing cholangitis", "Metabolic and Other Disorders", + "gastritis", "Metabolic and Other Disorders", + "acute kidney failure", "Metabolic and Other Disorders", + "tubular adenoma", "Metabolic and Other Disorders", + "benign prostatic hyperplasia", "Metabolic and Other Disorders", + "opiate dependence", "Metabolic and Other Disorders", + "gingivitis", "Metabolic and Other Disorders", + "hyperplastic polyp", "Metabolic and Other Disorders", + "clonal hematopoiesis", "Metabolic and Other Disorders", + "epilepsy", "Metabolic and Other Disorders", + "age related macular degeneration 7", "Metabolic and Other Disorders", + "kidney benign neoplasm", "Metabolic and Other Disorders", + "malignant pancreatic neoplasm", "Metabolic and Other Disorders", + "cataract", "Metabolic and Other Disorders", + "macular degeneration", "Metabolic and Other Disorders", + "hydrosalpinx", "Metabolic and Other Disorders", + "tubulovillous adenoma", "Metabolic and Other Disorders", + "gastric intestinal metaplasia", "Metabolic and Other Disorders", + "Barrett esophagus", "Metabolic and Other Disorders", + + # Other Diseases + "injury", "Other Diseases", + "anencephaly", "Other Diseases", + "primary biliary cholangitis", "Other Diseases", + "keloid", "Other Diseases", + "kidney oncocytoma", "Other Diseases", + "respiratory failure", "Other Diseases", + "pilocytic astrocytoma", "Other Diseases" + ) + + temp_path = '/scratch/pawsey1192/zhanchen/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2/temp' + system(glue("rclone copy UofA_Box:/minh_immune_map_disease/disease_data_grouped_further.csv {temp_path}/")) + + disease_data_grouped = + disease_data_grouped |> + left_join( + read_csv(glue("{temp_path}/disease_data_grouped_further.csv")) |> + rename(disease_groups_further = disease_groups) + ) |> + mutate(disease_groups = if_else(!disease_groups_further |> is.na(), disease_groups_further, disease_groups)) |> + select(disease, disease_groups) + + tbl |> + + # TISSUE + select(-any_of("tissue_groups")) |> + left_join(tissue_grouped, copy=TRUE) |> + + # TECH + left_join(assay_data_grouped, copy=TRUE) |> + + # DISEASE + left_join(disease_data_grouped, copy=TRUE) |> + + # make disease tissue specific, omit Normal + mutate(disease_groups = paste(disease_groups, tissue_groups, sep = "_")) |> + mutate(disease_groups = if_else(disease_groups |> str_detect("Normal_.+"), "Normal", disease_groups)) |> + + + # TEMPORARY. de-group pancreas and liver + mutate(tissue_groups = case_when( + + tissue %in% c("gallbladder") ~ "gallbladder", + tissue %in% c("pancreas", "exocrine pancreas") ~ "pancreas", + tissue %in% c("liver", "caudate lobe of liver", "hepatic cecum" ) ~ "liver", + TRUE ~ tissue_groups + )) |> + + # SEX edit + mutate(sex = if_else(sex |> is.na(), "unknown", sex)) |> + + # Age + filter(age_days > 365) |> + mutate(age_years = age_days / 365) |> + mutate(age_bin = dplyr::case_when( + age_years < 3 ~ "Infancy", + age_years < 12 ~ "Childhood", + age_years < 20 ~ "Adolescence", + age_years < 40 ~ "Young Adulthood", + age_years < 50 ~ "Middle Age", + age_years < 60 ~ "Senior_50", + age_years < 70 ~ "Senior_60", + age_years >= 70 ~ "Senior_70", + TRUE ~ NA_character_ + )) |> + mutate(age_decade = ceiling(age_years/10) |> as.character()) |> + + # left_join(age_bin_table, copy=TRUE) |> + + # ETHNICITY + left_join(ethnicity_grouped, copy=TRUE) |> + + dplyr::select( + sample_id, donor_id, dataset_id, title, collection_id, age_days, age_bin, age_decade, sex, + ethnicity_groups, tissue_groups, tissue, assay_groups, cell_type_unified_ensemble, + cell_type, disease_groups, is_immune + ) |> + as_tibble() |> + + # Center based on adolescence + mutate(age_days_scaled = age_days |> scale(center = 50*365) |> as.numeric()) + + } + + create_input_cell_counts = function(cellNexus_metadata, drop_sample_df, caq_celltype_level_map, ethnicity_imputed, result_directory){ + + tbl( + dbConnect(duckdb::duckdb(), dbdir = ":memory:"), + sql(glue("SELECT * FROM read_parquet('{cellNexus_metadata}')")) + ) |> + + # Filter empty droplets + filter(!empty_droplet) |> + + # TISSUE + filter(!tissue_groups %in% c( + "muscular system (skeletal muscles)", + "ovary", + "vasculature", + "digestive tract junctions and connections", + "peritoneal and abdominal cavity structures", + "connective tissue", + "miscellaneous glands" + ), + !tissue_groups |> is.na() + ) |> + + # NON immune cells + mutate(cell_type_unified_ensemble = if_else(is_immune, cell_type_unified_ensemble, "non_immune")) |> + + # IMMUNE CELLS + # filter(is_immune) |> + filter(cell_type_unified_ensemble %in% c("non_immune", "cd8 naive", "cd16 mono", "cd4 tcm", "cd4 th17 em", "granulocyte", "cd4 th1/th17 em", "treg", "b memory", "b naive", "nk", "plasma", "cd4 th2 em", "mast", "cd4 th1 em", "cd8 tem", "mait", "tgd", "cdc", "cd4 fh em", "cd4 naive", "nkt", "macrophage", "cytotoxic", "cd8 tcm", "cd14 mono", "pdc", "ilc")) |> + + edit_covariates() |> + + # Here we drop those samples with a low cell type entropy. E.g. one cell type only. + anti_join(drop_sample_df, copy = TRUE) |> + + dplyr::count( + sample_id, donor_id, dataset_id, title, collection_id, age_days, age_bin, age_days_scaled, age_decade, + sex, ethnicity_groups, tissue_groups, tissue, assay_groups, cell_type_unified_ensemble, is_immune, + disease_groups) |> + mutate(n = as.integer(n)) |> + as_tibble() |> + + # Add hierarchy of cell types L1, L2, L3 + left_join( + caq_celltype_level_map, + by = join_by(cell_type_unified_ensemble == cell_type_unified_harmonised) + ) |> + mutate( + across( + matches("^L[0-9]"), # select columns whose names start L0, L1, … L9 + ~ if_else(is_immune, as.character(.), "non_immune") + ) + ) |> + + # Add imputed ethnicities, and assign original if not present (cmposition and DE might have different samples because of filtering) + left_join( + ethnicity_imputed, + by = join_by(sample_id, ethnicity_groups) + ) |> + mutate(ethnicity_groups_imputed = if_else(ethnicity_groups_imputed |> is.na(), ethnicity_groups, ethnicity_groups_imputed)) |> + + # Remove confounders of non interest + tidybulk:::.resolve_complete_confounders_of_non_interest_df(dataset_id, assay_groups, disease_groups) |> + + # Set intercept + mutate( + ethnicity_groups_imputed = fct_relevel(ethnicity_groups_imputed, "European"), + assay_groups___altered = fct_relevel(assay_groups___altered, "10x Genomics 3"), + disease_groups___altered = fct_relevel(disease_groups___altered, "Normal"), + age_bin = fct_relevel(age_bin, "Senior_50"), + age_decade = fct_relevel(age_decade, "50") + ) + + } + + check_rclone_installation = function(){ + rclone_path <- Sys.which("rclone") + if (rclone_path == "") { + stop("rclone is not installed or not found in the system PATH.") + } + } + + #-----------------------# + # Pipeline + #-----------------------# + list( + tar_target( + result_directory, + "/scratch/pawsey1192/zhanchen/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2", + deployment = "main" + ), + tar_target( + drop_samples, + { + # evaluate result_directory for targets + print(result_directory) + + check_rclone_installation() + system(glue("rclone copy UofA_Box:/Mangiola_ImmuneAtlas/dharmesh_shared_mix/drop_samples.csv {result_directory}/")) + + read_csv(glue("{result_directory}/drop_samples.csv")) + + }, packages = c("glue", "readr") + ), + tar_target( + ethnicity_imputed, + { + check_rclone_installation() + temp_path = '/scratch/pawsey1192/zhanchen/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2/temp' + system(glue("rclone copy UofA_Box:/Mangiola_ImmuneAtlas/reports/ning/data/All_pseudobulk_1_0_6_ethnicity_imputed_colData.csv {temp_path}/")) + + read_csv(glue("{temp_path}/All_pseudobulk_1_0_6_ethnicity_imputed_colData.csv")) |> + select(sample_id, ethnicity_groups, ethnicity_groups_imputed = finalEthnicity_groups) |> + mutate(ethnicity_groups_imputed = ethnicity_groups_imputed |> str_replace("_imp$", "_imputed")) + + }, packages = c("glue", "readr", "dplyr", "stringr") + ), + tar_target( + caq_celltype_level_map, + { + check_rclone_installation() + temp_path = '/scratch/pawsey1192/zhanchen/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2/temp' + system(glue("rclone copy UofA_Box:/Mangiola_ImmuneAtlas/reannotation_consensus/caq_celltype_level_map.csv {temp_path}/")) + read_csv(glue("{temp_path}/caq_celltype_level_map.csv")) + + }, packages = c("glue", "readr") + ), + + tar_target( + input_relative, + create_input_cell_counts( + "/scratch/pawsey1192/zhanchen/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseduobulk/metadata.1.0.10.parquet", + drop_samples, + caq_celltype_level_map, + ethnicity_imputed, + result_directory + ), + packages = c( + "dplyr", # Data manipulation (mutate, filter, count, left_join) + "tidyr", + "tibble", # Creating tibbles (tribble function) + "duckdb", # For connecting to DuckDB and reading parquet files + "glue", # For constructing SQL queries with glue syntax + "readr", # For reading CSV files + "forcats", + "tidybulk" + ) + ), + tar_target( + saved_input_relative, + { + # evaluate result_directory for targets + print(result_directory) + + file_name = glue("{result_directory}/cell_metadata_1_0_10_sccomp_input_counts.rds") + input_relative |> + saveRDS(file_name) + + check_rclone_installation() + system(glue("rclone copy {file_name} UofA_Box:/Mangiola_ImmuneAtlas/taskforce_shared_folder/")) + + }, + packages = "glue" + ), + tar_target( + formula_df, + tribble( + ~ formula_composition, ~ formula_variability, ~ name, + + # continuous age + "~ 1 + age_days_scaled*sex + disease_groups___altered + ethnicity_groups_imputed + assay_groups___altered + + (1 | dataset_id___altered) + + (1 + age_days_scaled*sex + ethnicity_groups_imputed | tissue_groups)", + "~ age_days_scaled*sex + disease_groups___altered", + "estimates_continuous_age", + + # discrete + "~ 1 + age_bin + disease_groups___altered + sex + age_bin:sex + ethnicity_groups_imputed + assay_groups___altered + + (1 | dataset_id___altered) + + (1 + age_bin + sex + age_bin:sex + ethnicity_groups_imputed | tissue_groups)", + "~ age_bin + disease_groups___altered", + "estimates_age_bins", + + # discrete decade + "~ 1 + age_decade + disease_groups___altered + sex + age_decade:sex + ethnicity_groups_imputed + assay_groups___altered + + (1 | dataset_id___altered) + + (1 + age_decade + sex + age_decade:sex + ethnicity_groups_imputed | tissue_groups)", + "~ disease_groups___altered", + "estimates_age_decade", + + # discrete + interaction ethnicity sex + "~ 1 + age_bin + disease_groups___altered + sex + age_bin:sex + ethnicity_groups_imputed * sex + assay_groups___altered + + (1 | dataset_id___altered) + + (1 + age_bin + sex + age_bin:sex + ethnicity_groups_imputed * sex | tissue_groups)", + "~ age_bin + disease_groups___altered", + "estimates_age_bins_sex_ethnicity_interaction", + + # continuous + discrete + "~ 1 + age_days_scaled + disease_groups___altered + age_bin*sex + ethnicity_groups_imputed + assay_groups___altered + + (1 | dataset_id___altered) + + (1 + age_days_scaled + age_bin*sex + ethnicity_groups_imputed | tissue_groups)", + "~ age_days_scaled + disease_groups___altered", + "estimates_continuous_age_plus_age_bins", + + # disease tissue specific + "~ 1 + age_bin + disease_groups___altered + sex + age_bin:sex + ethnicity_groups_imputed + assay_groups___altered + + (1 | dataset_id___altered) + + (1 + age_bin + disease_groups___altered + sex + age_bin:sex + ethnicity_groups_imputed | tissue_groups)", + "~ age_bin + disease_groups___altered", + "estimates_age_bins_disease", + ) |> + expand_grid( + cell_type_level = glue("L{0:3}") |> as.character(), + drop_disease = c(TRUE, FALSE), + immune_only = c(TRUE,FALSE) + ) |> + + # Keep NON immune for one model only + filter(immune_only | name == "estimates_age_bins") |> + + mutate(counts = list(input_relative)) |> + + # Filter immune if needed + mutate(counts = map2(counts, immune_only, ~ { + if(.y) .x |> filter(is_immune) + else .x |> + nest(data = -c(sample_id, tissue_groups)) |> + filter(map_lgl(data, ~ .x |> filter(cell_type_unified_ensemble=="non_immune") |> nrow() > 0)) |> + mutate(non_immune_count = map_int(data, ~ .x |> filter(cell_type_unified_ensemble=="non_immune") |> pull(n))) |> + mutate(total_count = map_int(data, ~ .x |> pull(n) |> sum())) |> + mutate(proportion_non_immune = non_immune_count/total_count) |> + filter(proportion_non_immune > 2/3) |> + filter(!tissue_groups %in% c( + "bone marrow", + "lymphatic system", + "spleen", + "thymus", + "blood" + )) |> + unnest(data) + } )) |> + + # Drop Disease if needed + mutate(counts = map2(counts, drop_disease, + ~ { + if(.y) .x |> filter(disease_groups___altered == "Normal") + else .x + })) |> + mutate( + formula_composition = if_else(drop_disease, formula_composition |> str_remove_all("\\+ disease_groups___altered"), formula_composition), + formula_variability = if_else(drop_disease, formula_variability |> str_remove_all("\\+ disease_groups___altered"), formula_variability) + ) |> + mutate(local_file_name = glue("{name}___{cell_type_level}___disease_{!drop_disease}___immune_only_{immune_only}")) |> + group_by(local_file_name) |> + tar_group(), + iteration = "group", + packages = c("tibble", "glue", "targets", "dplyr", "tidyr", "purrr", "stringr") + ), + tar_target( + estimates, + + formula_df$counts[[1]] |> + + # With L0 I have to summarise further because I have counts already + with_groups( + c(sample_id,age_days_scaled,age_bin, sex, disease_groups___altered,ethnicity_groups_imputed, assay_groups___altered, dataset_id___altered, tissue_groups, all_of(formula_df$cell_type_level)), + ~ .x |> summarise(n = sum(n)) + ) |> + # + sccomp_estimate( + formula_composition = formula_df$formula_composition |> as.formula(), + formula_variability = formula_df$formula_variability |> as.formula(), # Differential variability + sample_column = "sample_id", + cell_group_column = formula_df$cell_type_level, # A level of the hierarchy + abundance_column = "n", + cores = as.numeric(Sys.getenv("SLURM_CPUS_PER_TASK", unset = 1)), + mcmc_seed = 42, + verbose = T, + bimodal_mean_variability_association = TRUE, + prior_mean = list(intercept = c(0, 0.8), coefficients = c(0, 3)), + prior_overdispersion_mean_association = list(intercept = c(3.6539176, 0.5), slope = c(-0.5255242, 0.1), standard_deviation = c(20, 40)), + output_directory = "/scratch/pawsey1192/zhanchen/HPC_sccomp/my_draws", + max_sampling_iterations = 5000, + + # # TEMPORARY DEBUG + # max_sampling_iterations = 1, + # warmup_samples =1, + + inference_method = "hmc", + refresh = 1 + ), + pattern = map(formula_df), + resources = tar_resources(crew = tar_resources_crew("slurm_1_80")), + error = "continue", + packages = "sccomp" + #, + + # TEMPORARY + # cue = tar_cue(mode = "never") + ), + tar_target( + saved_and_tranferred, + { + # evaluate result_directory for targets + final_result_directory = "/scratch/pawsey1192/zhanchen/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2/results" + print(result_directory) + + local_file_name = glue("{final_result_directory}/{formula_df$local_file_name}.rds") + local_file_name_FIT_FOR_PORTABILITY = glue("{final_result_directory}/{formula_df$local_file_name}_FIT_FOR_PORTABILITY.rds") + + # Save draws as monolythic + attr(estimates, "fit")$save_object(file = local_file_name_FIT_FOR_PORTABILITY) + + # check_rclone_installation() + # system(glue("rclone copy {local_file_name_FIT_FOR_PORTABILITY} UofA_Box:/Mangiola_ImmuneAtlas/taskforce_shared_folder/sccomp_estimates_1_0_10/")) + estimates |> attr("fit") = readRDS(local_file_name_FIT_FOR_PORTABILITY) + + # Save sccomp estimates + estimates |> sccomp_test() |> saveRDS(local_file_name) + + # check_rclone_installation() + # system(glue("rclone copy {local_file_name} UofA_Box:/Mangiola_ImmuneAtlas/taskforce_shared_folder/sccomp_on_cellNexus_1_0_10_2/")) + + }, + pattern = map(formula_df, estimates), + resources = tar_resources(crew = tar_resources_crew("slurm_1_200")), + error = "continue", + packages = c("glue", "sccomp", "magrittr") + ) + + ) +}, +ask = FALSE, +script = "/scratch/pawsey1192/zhanchen/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2/_targets.R" +) \ No newline at end of file From 5e19ca4dcec4602ab278fcf8e0cb33d4e2181cd7 Mon Sep 17 00:00:00 2001 From: Chen Zhan Date: Fri, 13 Jun 2025 19:36:30 +0800 Subject: [PATCH 06/12] revised normalisation for DE pipeline --- .../PAWSEY_sccomp_tar_scripts.R | 45 +++++++++---------- .../de_brms_whole_immune_system.R | 38 +++++++++++++++- 2 files changed, 58 insertions(+), 25 deletions(-) diff --git a/rebuttal_CellPress/PAWSEY_sccomp_tar_scripts.R b/rebuttal_CellPress/PAWSEY_sccomp_tar_scripts.R index b15e9a7..7bc3d70 100644 --- a/rebuttal_CellPress/PAWSEY_sccomp_tar_scripts.R +++ b/rebuttal_CellPress/PAWSEY_sccomp_tar_scripts.R @@ -517,7 +517,7 @@ tar_script({ ) temp_path = '/scratch/pawsey1192/zhanchen/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2/temp' - system(glue("rclone copy UofA_Box:/minh_immune_map_disease/disease_data_grouped_further.csv {temp_path}/")) + # system(glue("rclone copy UofA_Box:/minh_immune_map_disease/disease_data_grouped_further.csv {temp_path}/")) disease_data_grouped = disease_data_grouped |> @@ -660,17 +660,17 @@ tar_script({ assay_groups___altered = fct_relevel(assay_groups___altered, "10x Genomics 3"), disease_groups___altered = fct_relevel(disease_groups___altered, "Normal"), age_bin = fct_relevel(age_bin, "Senior_50"), - age_decade = fct_relevel(age_decade, "50") + age_decade = fct_relevel(age_decade, "5.0") ) } - check_rclone_installation = function(){ - rclone_path <- Sys.which("rclone") - if (rclone_path == "") { - stop("rclone is not installed or not found in the system PATH.") - } - } + # check_rclone_installation = function(){ + # rclone_path <- Sys.which("rclone") + # if (rclone_path == "") { + # stop("rclone is not installed or not found in the system PATH.") + # } + # } #-----------------------# # Pipeline @@ -678,8 +678,8 @@ tar_script({ list( tar_target( result_directory, - "/scratch/pawsey1192/zhanchen/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2", - deployment = "main" + "/scratch/pawsey1192/zhanchen/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2" + # deployment = "main" ), tar_target( drop_samples, @@ -687,8 +687,8 @@ tar_script({ # evaluate result_directory for targets print(result_directory) - check_rclone_installation() - system(glue("rclone copy UofA_Box:/Mangiola_ImmuneAtlas/dharmesh_shared_mix/drop_samples.csv {result_directory}/")) + # check_rclone_installation() + # system(glue("rclone copy UofA_Box:/Mangiola_ImmuneAtlas/dharmesh_shared_mix/drop_samples.csv {result_directory}/")) read_csv(glue("{result_directory}/drop_samples.csv")) @@ -697,9 +697,9 @@ tar_script({ tar_target( ethnicity_imputed, { - check_rclone_installation() + # check_rclone_installation() temp_path = '/scratch/pawsey1192/zhanchen/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2/temp' - system(glue("rclone copy UofA_Box:/Mangiola_ImmuneAtlas/reports/ning/data/All_pseudobulk_1_0_6_ethnicity_imputed_colData.csv {temp_path}/")) + # system(glue("rclone copy UofA_Box:/Mangiola_ImmuneAtlas/reports/ning/data/All_pseudobulk_1_0_6_ethnicity_imputed_colData.csv {temp_path}/")) read_csv(glue("{temp_path}/All_pseudobulk_1_0_6_ethnicity_imputed_colData.csv")) |> select(sample_id, ethnicity_groups, ethnicity_groups_imputed = finalEthnicity_groups) |> @@ -710,9 +710,9 @@ tar_script({ tar_target( caq_celltype_level_map, { - check_rclone_installation() + # check_rclone_installation() temp_path = '/scratch/pawsey1192/zhanchen/HPC_sccomp/sccomp_on_cellNexus_1_0_10_2/temp' - system(glue("rclone copy UofA_Box:/Mangiola_ImmuneAtlas/reannotation_consensus/caq_celltype_level_map.csv {temp_path}/")) + # system(glue("rclone copy UofA_Box:/Mangiola_ImmuneAtlas/reannotation_consensus/caq_celltype_level_map.csv {temp_path}/")) read_csv(glue("{temp_path}/caq_celltype_level_map.csv")) }, packages = c("glue", "readr") @@ -748,8 +748,8 @@ tar_script({ input_relative |> saveRDS(file_name) - check_rclone_installation() - system(glue("rclone copy {file_name} UofA_Box:/Mangiola_ImmuneAtlas/taskforce_shared_folder/")) + # check_rclone_installation() + # system(glue("rclone copy {file_name} UofA_Box:/Mangiola_ImmuneAtlas/taskforce_shared_folder/")) }, packages = "glue" @@ -855,16 +855,15 @@ tar_script({ # With L0 I have to summarise further because I have counts already with_groups( - c(sample_id,age_days_scaled,age_bin, sex, disease_groups___altered,ethnicity_groups_imputed, assay_groups___altered, dataset_id___altered, tissue_groups, all_of(formula_df$cell_type_level)), - ~ .x |> summarise(n = sum(n)) + c(sample_id,age_days_scaled,age_bin, age_decade, sex, disease_groups___altered,ethnicity_groups_imputed, assay_groups___altered, dataset_id___altered, tissue_groups, all_of(formula_df$cell_type_level)), ~ .x |> summarise(n = sum(n)) ) |> # sccomp_estimate( formula_composition = formula_df$formula_composition |> as.formula(), formula_variability = formula_df$formula_variability |> as.formula(), # Differential variability - sample_column = "sample_id", - cell_group_column = formula_df$cell_type_level, # A level of the hierarchy - abundance_column = "n", + sample = "sample_id", + cell_group = formula_df$cell_type_level, # A level of the hierarchy + abundance = "n", cores = as.numeric(Sys.getenv("SLURM_CPUS_PER_TASK", unset = 1)), mcmc_seed = 42, verbose = T, diff --git a/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R b/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R index e5eeee6..fa402a5 100644 --- a/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R +++ b/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R @@ -809,6 +809,39 @@ tar_script({ mutate(age_days_scaled = age_days |> scale(center = 50*365) |> as.numeric()) } + + offset_calcuation = function(se, method = 'TMMwsp', reference_sample){ + + # Check if package is installed, otherwise install + tidybulk:::check_and_install_packages("edgeR") + + # Drop genes with NAs, as edgeR::calcNormFactors does not accept them + my_counts_filtered = se %>% assays() %>% as.list() %>% .[[1]] %>% na.omit() + # Calcuate library size + library_size_filtered = my_counts_filtered %>% colSums(na.rm = TRUE) + + # Calculate TMM + nf <- + edgeR::calcNormFactors( + my_counts_filtered, + refColumn = reference_sample, + method = method + ) + + # Calculate multiplier + multiplier = library_size_filtered[reference_sample] * nf[reference_sample] %>% divide_by(library_size_filtered * nf) + + # Calcuate offset + offset = log(1/multiplier) + + # Add to sample info + colData(se)$normalisation_factor = nf + colData(se)$multiplier = multiplier + colData(se)$offset = offset + + return(se) + + } #-----------------------# # Pipeline @@ -1054,8 +1087,9 @@ tar_script({ ) |> # Get scaling factor - scale_abundance(method = "TMMwsp", reference_sample = reference_sample) |> - + # scale_abundance(method = "TMMwsp", reference_sample = reference_sample) |> + offset_calcuation(method = "TMMwsp", reference_sample = reference_sample) |> + # Drop sex unknown as causes problem during fit mutate( sex = if_else(sex |> is.na(), "unknown", sex), From c5c213c1212a83eac94222fa9b10cfe167b5d5d0 Mon Sep 17 00:00:00 2001 From: Chen Zhan Date: Wed, 2 Jul 2025 15:27:39 +0800 Subject: [PATCH 07/12] minor change --- .../whole_immune_system_analyses/de_brms_whole_immune_system.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R b/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R index fa402a5..e638b9a 100644 --- a/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R +++ b/rebuttal_CellPress/whole_immune_system_analyses/de_brms_whole_immune_system.R @@ -36,7 +36,7 @@ tar_script({ ## Pawsey setting ----- hdf5_path = "/scratch/pawsey1192/zhanchen/Mangiola_ImmuneAtlas/taskforce_shared_folder/pseudobulk_sample_cell_type" - target_cell_type = "cd4 fh em" + target_cell_type = "cd4 fh em" tar_option_set( From 623dfcaf7d2a3ed8de75fa4a3aa479dae2fef8e8 Mon Sep 17 00:00:00 2001 From: Chen Zhan Date: Wed, 2 Jul 2025 15:30:49 +0800 Subject: [PATCH 08/12] Add edit_covariates.R from origin/master --- rebuttal_CellPress/edit_covariates.R | 511 +++++++++++++++++++++++++++ 1 file changed, 511 insertions(+) create mode 100644 rebuttal_CellPress/edit_covariates.R diff --git a/rebuttal_CellPress/edit_covariates.R b/rebuttal_CellPress/edit_covariates.R new file mode 100644 index 0000000..36e33ef --- /dev/null +++ b/rebuttal_CellPress/edit_covariates.R @@ -0,0 +1,511 @@ +edit_covariates = function(tbl){ + + tissue_grouped = 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", + + # Epithelium + "epithelium of esophagus" + ), + + # Stomach + "stomach" = c( + "stomach", "body of stomach", "cardia of stomach" + ), + + # Small Intestine + "small intestine" = c( + "small intestine", "duodenum", "jejunum", "ileum", + + # Epithelium + "epithelium of small intestine", "jejunal epithelium", "ileal epithelium", + "submucosa of ileum", "lamina propria of small intestine" + ), + + # 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", + + # epithelium + "colonic epithelium", "submucosa of ascending colon", "lamina propria of large intestine", + "mucosa of colon", "lamina propria of mucosa of colon", "caecum epithelium" + ), + + # Digestive System (General) + "digestive system (general)" = c( + "intestine", "hindgut", "lamina propria", "mucosa" + ), + + # 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" + ), + + # 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) |> + rename(tissue = value) |> + mutate() + + ethnicity_grouped <- tribble( + ~self_reported_ethnicity, ~ethnicity_groups, + "unknown", "Other/Unknown", + "European", "European", + "Korean", "East Asian", + "Asian", "East Asian", + "Japanese", "Japanese", + "African American", "African", + "Hispanic or Latin American", "Hispanic/Latin American", + "Singaporean Chinese", "East Asian", + "Han Chinese", "East Asian", + "Singaporean Indian", "South Asian", + "Singaporean Malay", "Other/Unknown", + "British", "European", + "African", "African", + "South Asian", "South Asian", + "European American", "European", + "East Asian", "East Asian", + "American", "Other/Unknown", + "African American or Afro-Caribbean", "African", + "Oceanian", "Native American & Pacific Islander", + "Jewish Israeli", "Middle Eastern & North African", + "Chinese", "East Asian", + "South East Asian", "Other/Unknown", + "Greater Middle Eastern (Middle Eastern or North African or Persian)", "Middle Eastern & North African", + "Native American", "Native American & Pacific Islander", + "Pacific Islander", "Native American & Pacific Islander", + "Finnish", "European", + "Bangladeshi", "South Asian", + "Native American,Hispanic or Latin American", "Hispanic/Latin American", + "Irish", "European", + "Iraqi", "Middle Eastern & North African", + "European,Asian", "European" + ) + + assay_data_grouped <- tribble( + ~assay, ~assay_groups, + "10x 3' v2", "10x Genomics 3", + "10x 3' v3", "10x Genomics 3", + "10x 5' v2", "10x Genomics 5", + "10x 5' v1", "10x Genomics 5", + "MARS-seq", "Plate based Technologies", + "10x 3' transcription profiling", "10x Genomics 3", + "10x 5' transcription profiling", "10x Genomics 5", + "Smart-seq2", "Smart seq", + "microwell-seq", "Microwell Technologies", + "TruDrop", "TruDrop", + "Drop-seq", "Drop based Technologies", + "Seq-Well S3", "Microwell Technologies", + "GEXSCOPE technology", "Other Technologies", + "Seq-Well", "Microwell Technologies", + "sci-RNA-seq", "Other Technologies", + "10x 3' v1", "10x Genomics 3", + "BD Rhapsody Whole Transcriptome Analysis", "Other Technologies", + "BD Rhapsody Targeted mRNA", "Other Technologies", + "CEL-seq2", "Plate based Technologies", + "SPLiT-seq", "Other Technologies", + "STRT-seq", "Plate based Technologies", + "inDrop", "Drop based Technologies", + "Smart-seq v4", "Smart seq", + "ScaleBio single cell RNA sequencing", "Other Technologies" + ) + + + disease_data_grouped <- tribble( + ~disease, ~disease_groups, + + # Normal control + "normal", "Normal", + + # Isolated Diseases + "COVID-19", "COVID-19 related", + "post-COVID-19 disorder", "COVID-19 related", + "long COVID-19", "COVID-19 related", + "glioblastoma", "Glioblastoma", + "lung adenocarcinoma", "Lung Adenocarcinoma", + "systemic lupus erythematosus", "Systemic Lupus Erythematosus", + + # Infectious and Immune-related Diseases (other than COVID-19) + "Crohn disease", "Infectious and Immune-related Diseases", + "Crohn ileitis", "Infectious and Immune-related Diseases", + "pneumonia", "Infectious and Immune-related Diseases", + "common variable immunodeficiency", "Infectious and Immune-related Diseases", + "toxoplasmosis", "Infectious and Immune-related Diseases", + "Plasmodium malariae malaria", "Infectious and Immune-related Diseases", + "type 1 diabetes mellitus", "Infectious and Immune-related Diseases", + "influenza", "Infectious and Immune-related Diseases", + "chronic rhinitis", "Infectious and Immune-related Diseases", + "periodontitis", "Infectious and Immune-related Diseases", + "localized scleroderma", "Infectious and Immune-related Diseases", + "lymphangioleiomyomatosis", "Infectious and Immune-related Diseases", + "listeriosis", "Infectious and Immune-related Diseases", + + # Cancer (other than isolated cancers) + "squamous cell lung carcinoma", "Cancer", + "small cell lung carcinoma", "Cancer", + "non-small cell lung carcinoma", "Cancer", + "breast carcinoma", "Cancer", + "breast cancer", "Cancer", + "luminal B breast carcinoma", "Cancer", + "luminal A breast carcinoma", "Cancer", + "triple-negative breast carcinoma", "Cancer", + "gastric cancer", "Cancer", + "colorectal cancer", "Cancer", + "colon sessile serrated adenoma/polyp", "Cancer", + "follicular lymphoma", "Cancer", + "B-cell acute lymphoblastic leukemia", "Cancer", + "B-cell non-Hodgkin lymphoma", "Cancer", + "acute myeloid leukemia", "Cancer", + "acute promyelocytic leukemia", "Cancer", + "plasma cell myeloma", "Cancer", + "clear cell renal carcinoma", "Cancer", + "nonpapillary renal cell carcinoma", "Cancer", + "basal cell carcinoma", "Cancer", + "colorectal neoplasm", "Cancer", + "adenocarcinoma", "Cancer", + "chromophobe renal cell carcinoma", "Cancer", + "neuroendocrine carcinoma", "Cancer", + "lung large cell carcinoma", "Cancer", + "tongue cancer", "Cancer", + "Wilms tumor", "Cancer", + "pleomorphic carcinoma", "Cancer", + "blastoma", "Cancer", + + # Neurodegenerative and Neurological Disorders + "dementia", "Neurodegenerative and Neurological Disorders", + "Alzheimer disease", "Neurodegenerative and Neurological Disorders", + "Parkinson disease", "Neurodegenerative and Neurological Disorders", + "amyotrophic lateral sclerosis", "Neurodegenerative and Neurological Disorders", + "multiple sclerosis", "Neurodegenerative and Neurological Disorders", + "Down syndrome", "Neurodegenerative and Neurological Disorders", + "trisomy 18", "Neurodegenerative and Neurological Disorders", + "frontotemporal dementia", "Neurodegenerative and Neurological Disorders", + "temporal lobe epilepsy", "Neurodegenerative and Neurological Disorders", + "Lewy body dementia", "Neurodegenerative and Neurological Disorders", + "amyotrophic lateral sclerosis 26 with or without frontotemporal dementia", "Neurodegenerative and Neurological Disorders", + + # Respiratory Conditions + "pulmonary fibrosis", "Respiratory Conditions", + "respiratory system disorder", "Respiratory Conditions", + "chronic obstructive pulmonary disease", "Respiratory Conditions", + "cystic fibrosis", "Respiratory Conditions", + "interstitial lung disease", "Respiratory Conditions", + "hypersensitivity pneumonitis", "Respiratory Conditions", + "non-specific interstitial pneumonia", "Respiratory Conditions", + "aspiration pneumonia", "Respiratory Conditions", + "pulmonary emphysema", "Respiratory Conditions", + "pulmonary sarcoidosis", "Respiratory Conditions", + + # Cardiovascular Diseases + "myocardial infarction", "Cardiovascular Diseases", + "acute myocardial infarction", "Cardiovascular Diseases", + "dilated cardiomyopathy", "Cardiovascular Diseases", + "heart failure", "Cardiovascular Diseases", + "arrhythmogenic right ventricular cardiomyopathy", "Cardiovascular Diseases", + "congenital heart disease", "Cardiovascular Diseases", + "non-compaction cardiomyopathy", "Cardiovascular Diseases", + "cardiomyopathy", "Cardiovascular Diseases", + "heart disorder", "Cardiovascular Diseases", + + # Metabolic and Other Disorders + "type 2 diabetes mellitus", "Metabolic and Other Disorders", + "chronic kidney disease", "Metabolic and Other Disorders", + "digestive system disorder", "Metabolic and Other Disorders", + "primary sclerosing cholangitis", "Metabolic and Other Disorders", + "gastritis", "Metabolic and Other Disorders", + "acute kidney failure", "Metabolic and Other Disorders", + "tubular adenoma", "Metabolic and Other Disorders", + "benign prostatic hyperplasia", "Metabolic and Other Disorders", + "opiate dependence", "Metabolic and Other Disorders", + "gingivitis", "Metabolic and Other Disorders", + "hyperplastic polyp", "Metabolic and Other Disorders", + "clonal hematopoiesis", "Metabolic and Other Disorders", + "epilepsy", "Metabolic and Other Disorders", + "age related macular degeneration 7", "Metabolic and Other Disorders", + "kidney benign neoplasm", "Metabolic and Other Disorders", + "malignant pancreatic neoplasm", "Metabolic and Other Disorders", + "cataract", "Metabolic and Other Disorders", + "macular degeneration", "Metabolic and Other Disorders", + "hydrosalpinx", "Metabolic and Other Disorders", + "tubulovillous adenoma", "Metabolic and Other Disorders", + "gastric intestinal metaplasia", "Metabolic and Other Disorders", + "Barrett esophagus", "Metabolic and Other Disorders", + + # Other Diseases + "injury", "Other Diseases", + "anencephaly", "Other Diseases", + "primary biliary cholangitis", "Other Diseases", + "keloid", "Other Diseases", + "kidney oncocytoma", "Other Diseases", + "respiratory failure", "Other Diseases", + "pilocytic astrocytoma", "Other Diseases" + ) + + disease_data_grouped = + disease_data_grouped |> + left_join( + readr::read_csv("/home/zhanchen/From_scratch/Mangiola_ImmuneAtlas/disease_data_grouped_further.csv") |> + rename(disease_groups_further = disease_groups) + ) |> + mutate(disease_groups = if_else(!disease_groups_further |> is.na(), disease_groups_further, disease_groups)) |> + select(disease, disease_groups) + + tbl |> + + # TISSUE + select(-any_of("tissue_groups")) |> + left_join(tissue_grouped, copy=TRUE) |> + + # TECH + left_join(assay_data_grouped, copy=TRUE) |> + + # DISEASE + left_join(disease_data_grouped, copy=TRUE) |> + + # make disease tissue specific, omit Normal + mutate(disease_groups = paste(disease_groups, tissue_groups, sep = "_")) |> + mutate(disease_groups = if_else(disease_groups |> str_detect("Normal_.+"), "Normal", disease_groups)) |> + + + # TEMPORARY. de-group pancreas and liver + mutate(tissue_groups = case_when( + + tissue %in% c("gallbladder") ~ "gallbladder", + tissue %in% c("pancreas", "exocrine pancreas") ~ "pancreas", + tissue %in% c("liver", "caudate lobe of liver", "hepatic cecum" ) ~ "liver", + TRUE ~ tissue_groups + )) |> + + # SEX edit + mutate(sex = if_else(sex |> is.na(), "unknown", sex)) |> + + # Age + filter(age_days > 365) |> + mutate(age_years = age_days / 365) |> + mutate(age_bin = dplyr::case_when( + age_years < 3 ~ "Infancy", + age_years < 12 ~ "Childhood", + age_years < 20 ~ "Adolescence", + age_years < 40 ~ "Young Adulthood", + age_years < 50 ~ "Middle Age", + age_years < 60 ~ "Senior_50", + age_years < 70 ~ "Senior_60", + age_years >= 70 ~ "Senior_70", + TRUE ~ NA_character_ + )) |> + + # left_join(age_bin_table, copy=TRUE) |> + + # ETHNICITY + left_join(ethnicity_grouped, copy=TRUE) |> + + dplyr::select( + sample_id, donor_id, dataset_id, title, collection_id, age_days, age_bin, sex, ethnicity_groups, + tissue_groups, tissue, assay_groups, cell_type_unified_ensemble, cell_type, + disease_groups + ) |> + as_tibble() |> + + # Center based on adolescence + mutate(age_days_scaled = age_days |> scale(center = 50*365) |> as.numeric()) + + } \ No newline at end of file From 24ff8a6a1aa210f178baf7b33d243b9b29d7c045 Mon Sep 17 00:00:00 2001 From: Chen Zhan Date: Wed, 2 Jul 2025 15:35:27 +0800 Subject: [PATCH 09/12] revert --- rebuttal_CellPress/edit_covariates.R | 511 --------------------------- 1 file changed, 511 deletions(-) delete mode 100644 rebuttal_CellPress/edit_covariates.R diff --git a/rebuttal_CellPress/edit_covariates.R b/rebuttal_CellPress/edit_covariates.R deleted file mode 100644 index 36e33ef..0000000 --- a/rebuttal_CellPress/edit_covariates.R +++ /dev/null @@ -1,511 +0,0 @@ -edit_covariates = function(tbl){ - - tissue_grouped = 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", - - # Epithelium - "epithelium of esophagus" - ), - - # Stomach - "stomach" = c( - "stomach", "body of stomach", "cardia of stomach" - ), - - # Small Intestine - "small intestine" = c( - "small intestine", "duodenum", "jejunum", "ileum", - - # Epithelium - "epithelium of small intestine", "jejunal epithelium", "ileal epithelium", - "submucosa of ileum", "lamina propria of small intestine" - ), - - # 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", - - # epithelium - "colonic epithelium", "submucosa of ascending colon", "lamina propria of large intestine", - "mucosa of colon", "lamina propria of mucosa of colon", "caecum epithelium" - ), - - # Digestive System (General) - "digestive system (general)" = c( - "intestine", "hindgut", "lamina propria", "mucosa" - ), - - # 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" - ), - - # 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) |> - rename(tissue = value) |> - mutate() - - ethnicity_grouped <- tribble( - ~self_reported_ethnicity, ~ethnicity_groups, - "unknown", "Other/Unknown", - "European", "European", - "Korean", "East Asian", - "Asian", "East Asian", - "Japanese", "Japanese", - "African American", "African", - "Hispanic or Latin American", "Hispanic/Latin American", - "Singaporean Chinese", "East Asian", - "Han Chinese", "East Asian", - "Singaporean Indian", "South Asian", - "Singaporean Malay", "Other/Unknown", - "British", "European", - "African", "African", - "South Asian", "South Asian", - "European American", "European", - "East Asian", "East Asian", - "American", "Other/Unknown", - "African American or Afro-Caribbean", "African", - "Oceanian", "Native American & Pacific Islander", - "Jewish Israeli", "Middle Eastern & North African", - "Chinese", "East Asian", - "South East Asian", "Other/Unknown", - "Greater Middle Eastern (Middle Eastern or North African or Persian)", "Middle Eastern & North African", - "Native American", "Native American & Pacific Islander", - "Pacific Islander", "Native American & Pacific Islander", - "Finnish", "European", - "Bangladeshi", "South Asian", - "Native American,Hispanic or Latin American", "Hispanic/Latin American", - "Irish", "European", - "Iraqi", "Middle Eastern & North African", - "European,Asian", "European" - ) - - assay_data_grouped <- tribble( - ~assay, ~assay_groups, - "10x 3' v2", "10x Genomics 3", - "10x 3' v3", "10x Genomics 3", - "10x 5' v2", "10x Genomics 5", - "10x 5' v1", "10x Genomics 5", - "MARS-seq", "Plate based Technologies", - "10x 3' transcription profiling", "10x Genomics 3", - "10x 5' transcription profiling", "10x Genomics 5", - "Smart-seq2", "Smart seq", - "microwell-seq", "Microwell Technologies", - "TruDrop", "TruDrop", - "Drop-seq", "Drop based Technologies", - "Seq-Well S3", "Microwell Technologies", - "GEXSCOPE technology", "Other Technologies", - "Seq-Well", "Microwell Technologies", - "sci-RNA-seq", "Other Technologies", - "10x 3' v1", "10x Genomics 3", - "BD Rhapsody Whole Transcriptome Analysis", "Other Technologies", - "BD Rhapsody Targeted mRNA", "Other Technologies", - "CEL-seq2", "Plate based Technologies", - "SPLiT-seq", "Other Technologies", - "STRT-seq", "Plate based Technologies", - "inDrop", "Drop based Technologies", - "Smart-seq v4", "Smart seq", - "ScaleBio single cell RNA sequencing", "Other Technologies" - ) - - - disease_data_grouped <- tribble( - ~disease, ~disease_groups, - - # Normal control - "normal", "Normal", - - # Isolated Diseases - "COVID-19", "COVID-19 related", - "post-COVID-19 disorder", "COVID-19 related", - "long COVID-19", "COVID-19 related", - "glioblastoma", "Glioblastoma", - "lung adenocarcinoma", "Lung Adenocarcinoma", - "systemic lupus erythematosus", "Systemic Lupus Erythematosus", - - # Infectious and Immune-related Diseases (other than COVID-19) - "Crohn disease", "Infectious and Immune-related Diseases", - "Crohn ileitis", "Infectious and Immune-related Diseases", - "pneumonia", "Infectious and Immune-related Diseases", - "common variable immunodeficiency", "Infectious and Immune-related Diseases", - "toxoplasmosis", "Infectious and Immune-related Diseases", - "Plasmodium malariae malaria", "Infectious and Immune-related Diseases", - "type 1 diabetes mellitus", "Infectious and Immune-related Diseases", - "influenza", "Infectious and Immune-related Diseases", - "chronic rhinitis", "Infectious and Immune-related Diseases", - "periodontitis", "Infectious and Immune-related Diseases", - "localized scleroderma", "Infectious and Immune-related Diseases", - "lymphangioleiomyomatosis", "Infectious and Immune-related Diseases", - "listeriosis", "Infectious and Immune-related Diseases", - - # Cancer (other than isolated cancers) - "squamous cell lung carcinoma", "Cancer", - "small cell lung carcinoma", "Cancer", - "non-small cell lung carcinoma", "Cancer", - "breast carcinoma", "Cancer", - "breast cancer", "Cancer", - "luminal B breast carcinoma", "Cancer", - "luminal A breast carcinoma", "Cancer", - "triple-negative breast carcinoma", "Cancer", - "gastric cancer", "Cancer", - "colorectal cancer", "Cancer", - "colon sessile serrated adenoma/polyp", "Cancer", - "follicular lymphoma", "Cancer", - "B-cell acute lymphoblastic leukemia", "Cancer", - "B-cell non-Hodgkin lymphoma", "Cancer", - "acute myeloid leukemia", "Cancer", - "acute promyelocytic leukemia", "Cancer", - "plasma cell myeloma", "Cancer", - "clear cell renal carcinoma", "Cancer", - "nonpapillary renal cell carcinoma", "Cancer", - "basal cell carcinoma", "Cancer", - "colorectal neoplasm", "Cancer", - "adenocarcinoma", "Cancer", - "chromophobe renal cell carcinoma", "Cancer", - "neuroendocrine carcinoma", "Cancer", - "lung large cell carcinoma", "Cancer", - "tongue cancer", "Cancer", - "Wilms tumor", "Cancer", - "pleomorphic carcinoma", "Cancer", - "blastoma", "Cancer", - - # Neurodegenerative and Neurological Disorders - "dementia", "Neurodegenerative and Neurological Disorders", - "Alzheimer disease", "Neurodegenerative and Neurological Disorders", - "Parkinson disease", "Neurodegenerative and Neurological Disorders", - "amyotrophic lateral sclerosis", "Neurodegenerative and Neurological Disorders", - "multiple sclerosis", "Neurodegenerative and Neurological Disorders", - "Down syndrome", "Neurodegenerative and Neurological Disorders", - "trisomy 18", "Neurodegenerative and Neurological Disorders", - "frontotemporal dementia", "Neurodegenerative and Neurological Disorders", - "temporal lobe epilepsy", "Neurodegenerative and Neurological Disorders", - "Lewy body dementia", "Neurodegenerative and Neurological Disorders", - "amyotrophic lateral sclerosis 26 with or without frontotemporal dementia", "Neurodegenerative and Neurological Disorders", - - # Respiratory Conditions - "pulmonary fibrosis", "Respiratory Conditions", - "respiratory system disorder", "Respiratory Conditions", - "chronic obstructive pulmonary disease", "Respiratory Conditions", - "cystic fibrosis", "Respiratory Conditions", - "interstitial lung disease", "Respiratory Conditions", - "hypersensitivity pneumonitis", "Respiratory Conditions", - "non-specific interstitial pneumonia", "Respiratory Conditions", - "aspiration pneumonia", "Respiratory Conditions", - "pulmonary emphysema", "Respiratory Conditions", - "pulmonary sarcoidosis", "Respiratory Conditions", - - # Cardiovascular Diseases - "myocardial infarction", "Cardiovascular Diseases", - "acute myocardial infarction", "Cardiovascular Diseases", - "dilated cardiomyopathy", "Cardiovascular Diseases", - "heart failure", "Cardiovascular Diseases", - "arrhythmogenic right ventricular cardiomyopathy", "Cardiovascular Diseases", - "congenital heart disease", "Cardiovascular Diseases", - "non-compaction cardiomyopathy", "Cardiovascular Diseases", - "cardiomyopathy", "Cardiovascular Diseases", - "heart disorder", "Cardiovascular Diseases", - - # Metabolic and Other Disorders - "type 2 diabetes mellitus", "Metabolic and Other Disorders", - "chronic kidney disease", "Metabolic and Other Disorders", - "digestive system disorder", "Metabolic and Other Disorders", - "primary sclerosing cholangitis", "Metabolic and Other Disorders", - "gastritis", "Metabolic and Other Disorders", - "acute kidney failure", "Metabolic and Other Disorders", - "tubular adenoma", "Metabolic and Other Disorders", - "benign prostatic hyperplasia", "Metabolic and Other Disorders", - "opiate dependence", "Metabolic and Other Disorders", - "gingivitis", "Metabolic and Other Disorders", - "hyperplastic polyp", "Metabolic and Other Disorders", - "clonal hematopoiesis", "Metabolic and Other Disorders", - "epilepsy", "Metabolic and Other Disorders", - "age related macular degeneration 7", "Metabolic and Other Disorders", - "kidney benign neoplasm", "Metabolic and Other Disorders", - "malignant pancreatic neoplasm", "Metabolic and Other Disorders", - "cataract", "Metabolic and Other Disorders", - "macular degeneration", "Metabolic and Other Disorders", - "hydrosalpinx", "Metabolic and Other Disorders", - "tubulovillous adenoma", "Metabolic and Other Disorders", - "gastric intestinal metaplasia", "Metabolic and Other Disorders", - "Barrett esophagus", "Metabolic and Other Disorders", - - # Other Diseases - "injury", "Other Diseases", - "anencephaly", "Other Diseases", - "primary biliary cholangitis", "Other Diseases", - "keloid", "Other Diseases", - "kidney oncocytoma", "Other Diseases", - "respiratory failure", "Other Diseases", - "pilocytic astrocytoma", "Other Diseases" - ) - - disease_data_grouped = - disease_data_grouped |> - left_join( - readr::read_csv("/home/zhanchen/From_scratch/Mangiola_ImmuneAtlas/disease_data_grouped_further.csv") |> - rename(disease_groups_further = disease_groups) - ) |> - mutate(disease_groups = if_else(!disease_groups_further |> is.na(), disease_groups_further, disease_groups)) |> - select(disease, disease_groups) - - tbl |> - - # TISSUE - select(-any_of("tissue_groups")) |> - left_join(tissue_grouped, copy=TRUE) |> - - # TECH - left_join(assay_data_grouped, copy=TRUE) |> - - # DISEASE - left_join(disease_data_grouped, copy=TRUE) |> - - # make disease tissue specific, omit Normal - mutate(disease_groups = paste(disease_groups, tissue_groups, sep = "_")) |> - mutate(disease_groups = if_else(disease_groups |> str_detect("Normal_.+"), "Normal", disease_groups)) |> - - - # TEMPORARY. de-group pancreas and liver - mutate(tissue_groups = case_when( - - tissue %in% c("gallbladder") ~ "gallbladder", - tissue %in% c("pancreas", "exocrine pancreas") ~ "pancreas", - tissue %in% c("liver", "caudate lobe of liver", "hepatic cecum" ) ~ "liver", - TRUE ~ tissue_groups - )) |> - - # SEX edit - mutate(sex = if_else(sex |> is.na(), "unknown", sex)) |> - - # Age - filter(age_days > 365) |> - mutate(age_years = age_days / 365) |> - mutate(age_bin = dplyr::case_when( - age_years < 3 ~ "Infancy", - age_years < 12 ~ "Childhood", - age_years < 20 ~ "Adolescence", - age_years < 40 ~ "Young Adulthood", - age_years < 50 ~ "Middle Age", - age_years < 60 ~ "Senior_50", - age_years < 70 ~ "Senior_60", - age_years >= 70 ~ "Senior_70", - TRUE ~ NA_character_ - )) |> - - # left_join(age_bin_table, copy=TRUE) |> - - # ETHNICITY - left_join(ethnicity_grouped, copy=TRUE) |> - - dplyr::select( - sample_id, donor_id, dataset_id, title, collection_id, age_days, age_bin, sex, ethnicity_groups, - tissue_groups, tissue, assay_groups, cell_type_unified_ensemble, cell_type, - disease_groups - ) |> - as_tibble() |> - - # Center based on adolescence - mutate(age_days_scaled = age_days |> scale(center = 50*365) |> as.numeric()) - - } \ No newline at end of file From 41ccd1a4488f17a934ec1a26aff8eb10b6e2a987 Mon Sep 17 00:00:00 2001 From: Chen Zhan Date: Wed, 2 Jul 2025 15:36:28 +0800 Subject: [PATCH 10/12] Add edit_covariates.R from origin/master --- rebuttal_CellPress/edit_covariates.R | 514 +++++++++++++++++++++++++++ 1 file changed, 514 insertions(+) create mode 100644 rebuttal_CellPress/edit_covariates.R diff --git a/rebuttal_CellPress/edit_covariates.R b/rebuttal_CellPress/edit_covariates.R new file mode 100644 index 0000000..7a3de89 --- /dev/null +++ b/rebuttal_CellPress/edit_covariates.R @@ -0,0 +1,514 @@ +edit_covariates = function(tbl){ + + tissue_grouped = 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", + + # Epithelium + "epithelium of esophagus" + ), + + # Stomach + "stomach" = c( + "stomach", "body of stomach", "cardia of stomach" + ), + + # Small Intestine + "small intestine" = c( + "small intestine", "duodenum", "jejunum", "ileum", + + # Epithelium + "epithelium of small intestine", "jejunal epithelium", "ileal epithelium", + "submucosa of ileum", "lamina propria of small intestine" + ), + + # 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", + + # epithelium + "colonic epithelium", "submucosa of ascending colon", "lamina propria of large intestine", + "mucosa of colon", "lamina propria of mucosa of colon", "caecum epithelium" + ), + + # Digestive System (General) + "digestive system (general)" = c( + "intestine", "hindgut", "lamina propria", "mucosa" + ), + + # 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" + ), + + # 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) + + ethnicity_grouped <- tribble( + ~self_reported_ethnicity, ~ethnicity_groups, + "unknown", "Other/Unknown", + "European", "European", + "Korean", "East Asian", + "Asian", "East Asian", + "Japanese", "Japanese", + "African American", "African", + "Hispanic or Latin American", "Hispanic/Latin American", + "Singaporean Chinese", "East Asian", + "Han Chinese", "East Asian", + "Singaporean Indian", "South Asian", + "Singaporean Malay", "Other/Unknown", + "British", "European", + "African", "African", + "South Asian", "South Asian", + "European American", "European", + "East Asian", "East Asian", + "American", "Other/Unknown", + "African American or Afro-Caribbean", "African", + "Oceanian", "Native American & Pacific Islander", + "Jewish Israeli", "Middle Eastern & North African", + "Chinese", "East Asian", + "South East Asian", "Other/Unknown", + "Greater Middle Eastern (Middle Eastern or North African or Persian)", "Middle Eastern & North African", + "Native American", "Native American & Pacific Islander", + "Pacific Islander", "Native American & Pacific Islander", + "Finnish", "European", + "Bangladeshi", "South Asian", + "Native American,Hispanic or Latin American", "Hispanic/Latin American", + "Irish", "European", + "Iraqi", "Middle Eastern & North African", + "European,Asian", "European" + ) + + assay_data_grouped <- tribble( + ~assay, ~assay_groups, + "10x 3' v2", "10x Genomics 3", + "10x 3' v3", "10x Genomics 3", + "10x 5' v2", "10x Genomics 5", + "10x 5' v1", "10x Genomics 5", + "MARS-seq", "Plate based Technologies", + "10x 3' transcription profiling", "10x Genomics 3", + "10x 5' transcription profiling", "10x Genomics 5", + "Smart-seq2", "Smart seq", + "microwell-seq", "Microwell Technologies", + "TruDrop", "TruDrop", + "Drop-seq", "Drop based Technologies", + "Seq-Well S3", "Microwell Technologies", + "GEXSCOPE technology", "Other Technologies", + "Seq-Well", "Microwell Technologies", + "sci-RNA-seq", "Other Technologies", + "10x 3' v1", "10x Genomics 3", + "BD Rhapsody Whole Transcriptome Analysis", "Other Technologies", + "BD Rhapsody Targeted mRNA", "Other Technologies", + "CEL-seq2", "Plate based Technologies", + "SPLiT-seq", "Other Technologies", + "STRT-seq", "Plate based Technologies", + "inDrop", "Drop based Technologies", + "Smart-seq v4", "Smart seq", + "ScaleBio single cell RNA sequencing", "Other Technologies" + ) + + + disease_data_grouped <- tribble( + ~disease, ~disease_groups, + + # Normal control + "normal", "Normal", + + # Isolated Diseases + "COVID-19", "COVID-19 related", + "post-COVID-19 disorder", "COVID-19 related", + "long COVID-19", "COVID-19 related", + "glioblastoma", "Glioblastoma", + "lung adenocarcinoma", "Lung Adenocarcinoma", + "systemic lupus erythematosus", "Systemic Lupus Erythematosus", + + # Infectious and Immune-related Diseases (other than COVID-19) + "Crohn disease", "Infectious and Immune-related Diseases", + "Crohn ileitis", "Infectious and Immune-related Diseases", + "pneumonia", "Infectious and Immune-related Diseases", + "common variable immunodeficiency", "Infectious and Immune-related Diseases", + "toxoplasmosis", "Infectious and Immune-related Diseases", + "Plasmodium malariae malaria", "Infectious and Immune-related Diseases", + "type 1 diabetes mellitus", "Infectious and Immune-related Diseases", + "influenza", "Infectious and Immune-related Diseases", + "chronic rhinitis", "Infectious and Immune-related Diseases", + "periodontitis", "Infectious and Immune-related Diseases", + "localized scleroderma", "Infectious and Immune-related Diseases", + "lymphangioleiomyomatosis", "Infectious and Immune-related Diseases", + "listeriosis", "Infectious and Immune-related Diseases", + + # Cancer (other than isolated cancers) + "squamous cell lung carcinoma", "Cancer", + "small cell lung carcinoma", "Cancer", + "non-small cell lung carcinoma", "Cancer", + "breast carcinoma", "Cancer", + "breast cancer", "Cancer", + "luminal B breast carcinoma", "Cancer", + "luminal A breast carcinoma", "Cancer", + "triple-negative breast carcinoma", "Cancer", + "gastric cancer", "Cancer", + "colorectal cancer", "Cancer", + "colon sessile serrated adenoma/polyp", "Cancer", + "follicular lymphoma", "Cancer", + "B-cell acute lymphoblastic leukemia", "Cancer", + "B-cell non-Hodgkin lymphoma", "Cancer", + "acute myeloid leukemia", "Cancer", + "acute promyelocytic leukemia", "Cancer", + "plasma cell myeloma", "Cancer", + "clear cell renal carcinoma", "Cancer", + "nonpapillary renal cell carcinoma", "Cancer", + "basal cell carcinoma", "Cancer", + "colorectal neoplasm", "Cancer", + "adenocarcinoma", "Cancer", + "chromophobe renal cell carcinoma", "Cancer", + "neuroendocrine carcinoma", "Cancer", + "lung large cell carcinoma", "Cancer", + "tongue cancer", "Cancer", + "Wilms tumor", "Cancer", + "pleomorphic carcinoma", "Cancer", + "blastoma", "Cancer", + + # Neurodegenerative and Neurological Disorders + "dementia", "Neurodegenerative and Neurological Disorders", + "Alzheimer disease", "Neurodegenerative and Neurological Disorders", + "Parkinson disease", "Neurodegenerative and Neurological Disorders", + "amyotrophic lateral sclerosis", "Neurodegenerative and Neurological Disorders", + "multiple sclerosis", "Neurodegenerative and Neurological Disorders", + "Down syndrome", "Neurodegenerative and Neurological Disorders", + "trisomy 18", "Neurodegenerative and Neurological Disorders", + "frontotemporal dementia", "Neurodegenerative and Neurological Disorders", + "temporal lobe epilepsy", "Neurodegenerative and Neurological Disorders", + "Lewy body dementia", "Neurodegenerative and Neurological Disorders", + "amyotrophic lateral sclerosis 26 with or without frontotemporal dementia", "Neurodegenerative and Neurological Disorders", + + # Respiratory Conditions + "pulmonary fibrosis", "Respiratory Conditions", + "respiratory system disorder", "Respiratory Conditions", + "chronic obstructive pulmonary disease", "Respiratory Conditions", + "cystic fibrosis", "Respiratory Conditions", + "interstitial lung disease", "Respiratory Conditions", + "hypersensitivity pneumonitis", "Respiratory Conditions", + "non-specific interstitial pneumonia", "Respiratory Conditions", + "aspiration pneumonia", "Respiratory Conditions", + "pulmonary emphysema", "Respiratory Conditions", + "pulmonary sarcoidosis", "Respiratory Conditions", + + # Cardiovascular Diseases + "myocardial infarction", "Cardiovascular Diseases", + "acute myocardial infarction", "Cardiovascular Diseases", + "dilated cardiomyopathy", "Cardiovascular Diseases", + "heart failure", "Cardiovascular Diseases", + "arrhythmogenic right ventricular cardiomyopathy", "Cardiovascular Diseases", + "congenital heart disease", "Cardiovascular Diseases", + "non-compaction cardiomyopathy", "Cardiovascular Diseases", + "cardiomyopathy", "Cardiovascular Diseases", + "heart disorder", "Cardiovascular Diseases", + + # Metabolic and Other Disorders + "type 2 diabetes mellitus", "Metabolic and Other Disorders", + "chronic kidney disease", "Metabolic and Other Disorders", + "digestive system disorder", "Metabolic and Other Disorders", + "primary sclerosing cholangitis", "Metabolic and Other Disorders", + "gastritis", "Metabolic and Other Disorders", + "acute kidney failure", "Metabolic and Other Disorders", + "tubular adenoma", "Metabolic and Other Disorders", + "benign prostatic hyperplasia", "Metabolic and Other Disorders", + "opiate dependence", "Metabolic and Other Disorders", + "gingivitis", "Metabolic and Other Disorders", + "hyperplastic polyp", "Metabolic and Other Disorders", + "clonal hematopoiesis", "Metabolic and Other Disorders", + "epilepsy", "Metabolic and Other Disorders", + "age related macular degeneration 7", "Metabolic and Other Disorders", + "kidney benign neoplasm", "Metabolic and Other Disorders", + "malignant pancreatic neoplasm", "Metabolic and Other Disorders", + "cataract", "Metabolic and Other Disorders", + "macular degeneration", "Metabolic and Other Disorders", + "hydrosalpinx", "Metabolic and Other Disorders", + "tubulovillous adenoma", "Metabolic and Other Disorders", + "gastric intestinal metaplasia", "Metabolic and Other Disorders", + "Barrett esophagus", "Metabolic and Other Disorders", + + # Other Diseases + "injury", "Other Diseases", + "anencephaly", "Other Diseases", + "primary biliary cholangitis", "Other Diseases", + "keloid", "Other Diseases", + "kidney oncocytoma", "Other Diseases", + "respiratory failure", "Other Diseases", + "pilocytic astrocytoma", "Other Diseases" + ) + + temp_path = tempdir() + system(glue("~/bin/rclone copy box_adelaide:/minh_immune_map_disease/disease_data_grouped_further.csv {temp_path}/")) + + disease_data_grouped = + disease_data_grouped |> + left_join( + read_csv(glue("{temp_path}/disease_data_grouped_further.csv")) |> + dplyr::rename(disease_groups_further = disease_groups) + ) |> + mutate(disease_groups = if_else(!disease_groups_further |> is.na(), disease_groups_further, disease_groups)) |> + select(disease, disease_groups) + + tbl |> + + # TISSUE + select(-any_of("tissue_groups")) |> + left_join(tissue_grouped, copy=TRUE) |> + + # TECH + left_join(assay_data_grouped, copy=TRUE) |> + + # DISEASE + left_join(disease_data_grouped, copy=TRUE) |> + + # make disease tissue specific, omit Normal + mutate(disease_groups = paste(disease_groups, tissue_groups, sep = "_")) |> + mutate(disease_groups = if_else(disease_groups |> str_detect("Normal_.+"), "Normal", disease_groups)) |> + + + # TEMPORARY. de-group pancreas and liver + mutate(tissue_groups = case_when( + + tissue %in% c("gallbladder") ~ "gallbladder", + tissue %in% c("pancreas", "exocrine pancreas") ~ "pancreas", + tissue %in% c("liver", "caudate lobe of liver", "hepatic cecum" ) ~ "liver", + TRUE ~ tissue_groups + )) |> + + # SEX edit + mutate(sex = if_else(sex |> is.na(), "unknown", sex)) |> + + # Age + filter(age_days > 365) |> + mutate(age_years = age_days / 365) |> + mutate(age_bin = dplyr::case_when( + age_years < 3 ~ "Infancy", + age_years < 12 ~ "Childhood", + age_years < 20 ~ "Adolescence", + age_years < 40 ~ "Young Adulthood", + age_years < 50 ~ "Middle Age", + age_years < 60 ~ "Senior_50", + age_years < 70 ~ "Senior_60", + age_years >= 70 ~ "Senior_70", + TRUE ~ NA_character_ + )) |> + mutate(age_decade = ceiling(age_years/10) |> as.integer() |> as.character()) |> + + # left_join(age_bin_table, copy=TRUE) |> + + # ETHNICITY + left_join(ethnicity_grouped, copy=TRUE) |> + + dplyr::select( + sample_id, donor_id, dataset_id, title, collection_id, age_days, age_bin, age_decade, sex, + ethnicity_groups, tissue_groups, tissue, assay_groups, cell_type_unified_ensemble, + cell_type, disease_groups, is_immune + ) |> + as_tibble() |> + + # Center based on adolescence + mutate(age_days_scaled = age_days |> scale(center = 50*365) |> as.numeric()) + +} From 961b36b8fd151c5a38c09c5dc6612f1ed9c87011 Mon Sep 17 00:00:00 2001 From: Chen Zhan Date: Wed, 2 Jul 2025 15:41:27 +0800 Subject: [PATCH 11/12] update edit_covariates.R, applicable to Pawsey --- rebuttal_CellPress/edit_covariates.R | 1021 +++++++++++++------------- 1 file changed, 509 insertions(+), 512 deletions(-) diff --git a/rebuttal_CellPress/edit_covariates.R b/rebuttal_CellPress/edit_covariates.R index 7a3de89..8170092 100644 --- a/rebuttal_CellPress/edit_covariates.R +++ b/rebuttal_CellPress/edit_covariates.R @@ -1,514 +1,511 @@ edit_covariates = function(tbl){ - - tissue_grouped = 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", - - # Epithelium - "epithelium of esophagus" - ), - - # Stomach - "stomach" = c( - "stomach", "body of stomach", "cardia of stomach" - ), - - # Small Intestine - "small intestine" = c( - "small intestine", "duodenum", "jejunum", "ileum", - - # Epithelium - "epithelium of small intestine", "jejunal epithelium", "ileal epithelium", - "submucosa of ileum", "lamina propria of small intestine" - ), - - # 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", - - # epithelium - "colonic epithelium", "submucosa of ascending colon", "lamina propria of large intestine", - "mucosa of colon", "lamina propria of mucosa of colon", "caecum epithelium" - ), - - # Digestive System (General) - "digestive system (general)" = c( - "intestine", "hindgut", "lamina propria", "mucosa" - ), - - # 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" - ), - - # 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" + + tissue_grouped = 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", + + # Epithelium + "epithelium of esophagus" + ), + + # Stomach + "stomach" = c( + "stomach", "body of stomach", "cardia of stomach" + ), + + # Small Intestine + "small intestine" = c( + "small intestine", "duodenum", "jejunum", "ileum", + + # Epithelium + "epithelium of small intestine", "jejunal epithelium", "ileal epithelium", + "submucosa of ileum", "lamina propria of small intestine" + ), + + # 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", + + # epithelium + "colonic epithelium", "submucosa of ascending colon", "lamina propria of large intestine", + "mucosa of colon", "lamina propria of mucosa of colon", "caecum epithelium" + ), + + # Digestive System (General) + "digestive system (general)" = c( + "intestine", "hindgut", "lamina propria", "mucosa" + ), + + # 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" + ), + + # 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) |> + rename(tissue = value) + + ethnicity_grouped <- tribble( + ~self_reported_ethnicity, ~ethnicity_groups, + "unknown", "Other/Unknown", + "European", "European", + "Korean", "East Asian", + "Asian", "East Asian", + "Japanese", "Japanese", + "African American", "African", + "Hispanic or Latin American", "Hispanic/Latin American", + "Singaporean Chinese", "East Asian", + "Han Chinese", "East Asian", + "Singaporean Indian", "South Asian", + "Singaporean Malay", "Other/Unknown", + "British", "European", + "African", "African", + "South Asian", "South Asian", + "European American", "European", + "East Asian", "East Asian", + "American", "Other/Unknown", + "African American or Afro-Caribbean", "African", + "Oceanian", "Native American & Pacific Islander", + "Jewish Israeli", "Middle Eastern & North African", + "Chinese", "East Asian", + "South East Asian", "Other/Unknown", + "Greater Middle Eastern (Middle Eastern or North African or Persian)", "Middle Eastern & North African", + "Native American", "Native American & Pacific Islander", + "Pacific Islander", "Native American & Pacific Islander", + "Finnish", "European", + "Bangladeshi", "South Asian", + "Native American,Hispanic or Latin American", "Hispanic/Latin American", + "Irish", "European", + "Iraqi", "Middle Eastern & North African", + "European,Asian", "European" ) - ) |> - enframe(name ="tissue_groups") |> - distinct() |> - unnest(value) |> - dplyr::rename(tissue = value) - - ethnicity_grouped <- tribble( - ~self_reported_ethnicity, ~ethnicity_groups, - "unknown", "Other/Unknown", - "European", "European", - "Korean", "East Asian", - "Asian", "East Asian", - "Japanese", "Japanese", - "African American", "African", - "Hispanic or Latin American", "Hispanic/Latin American", - "Singaporean Chinese", "East Asian", - "Han Chinese", "East Asian", - "Singaporean Indian", "South Asian", - "Singaporean Malay", "Other/Unknown", - "British", "European", - "African", "African", - "South Asian", "South Asian", - "European American", "European", - "East Asian", "East Asian", - "American", "Other/Unknown", - "African American or Afro-Caribbean", "African", - "Oceanian", "Native American & Pacific Islander", - "Jewish Israeli", "Middle Eastern & North African", - "Chinese", "East Asian", - "South East Asian", "Other/Unknown", - "Greater Middle Eastern (Middle Eastern or North African or Persian)", "Middle Eastern & North African", - "Native American", "Native American & Pacific Islander", - "Pacific Islander", "Native American & Pacific Islander", - "Finnish", "European", - "Bangladeshi", "South Asian", - "Native American,Hispanic or Latin American", "Hispanic/Latin American", - "Irish", "European", - "Iraqi", "Middle Eastern & North African", - "European,Asian", "European" - ) - - assay_data_grouped <- tribble( - ~assay, ~assay_groups, - "10x 3' v2", "10x Genomics 3", - "10x 3' v3", "10x Genomics 3", - "10x 5' v2", "10x Genomics 5", - "10x 5' v1", "10x Genomics 5", - "MARS-seq", "Plate based Technologies", - "10x 3' transcription profiling", "10x Genomics 3", - "10x 5' transcription profiling", "10x Genomics 5", - "Smart-seq2", "Smart seq", - "microwell-seq", "Microwell Technologies", - "TruDrop", "TruDrop", - "Drop-seq", "Drop based Technologies", - "Seq-Well S3", "Microwell Technologies", - "GEXSCOPE technology", "Other Technologies", - "Seq-Well", "Microwell Technologies", - "sci-RNA-seq", "Other Technologies", - "10x 3' v1", "10x Genomics 3", - "BD Rhapsody Whole Transcriptome Analysis", "Other Technologies", - "BD Rhapsody Targeted mRNA", "Other Technologies", - "CEL-seq2", "Plate based Technologies", - "SPLiT-seq", "Other Technologies", - "STRT-seq", "Plate based Technologies", - "inDrop", "Drop based Technologies", - "Smart-seq v4", "Smart seq", - "ScaleBio single cell RNA sequencing", "Other Technologies" - ) - - - disease_data_grouped <- tribble( - ~disease, ~disease_groups, - - # Normal control - "normal", "Normal", - - # Isolated Diseases - "COVID-19", "COVID-19 related", - "post-COVID-19 disorder", "COVID-19 related", - "long COVID-19", "COVID-19 related", - "glioblastoma", "Glioblastoma", - "lung adenocarcinoma", "Lung Adenocarcinoma", - "systemic lupus erythematosus", "Systemic Lupus Erythematosus", - - # Infectious and Immune-related Diseases (other than COVID-19) - "Crohn disease", "Infectious and Immune-related Diseases", - "Crohn ileitis", "Infectious and Immune-related Diseases", - "pneumonia", "Infectious and Immune-related Diseases", - "common variable immunodeficiency", "Infectious and Immune-related Diseases", - "toxoplasmosis", "Infectious and Immune-related Diseases", - "Plasmodium malariae malaria", "Infectious and Immune-related Diseases", - "type 1 diabetes mellitus", "Infectious and Immune-related Diseases", - "influenza", "Infectious and Immune-related Diseases", - "chronic rhinitis", "Infectious and Immune-related Diseases", - "periodontitis", "Infectious and Immune-related Diseases", - "localized scleroderma", "Infectious and Immune-related Diseases", - "lymphangioleiomyomatosis", "Infectious and Immune-related Diseases", - "listeriosis", "Infectious and Immune-related Diseases", - - # Cancer (other than isolated cancers) - "squamous cell lung carcinoma", "Cancer", - "small cell lung carcinoma", "Cancer", - "non-small cell lung carcinoma", "Cancer", - "breast carcinoma", "Cancer", - "breast cancer", "Cancer", - "luminal B breast carcinoma", "Cancer", - "luminal A breast carcinoma", "Cancer", - "triple-negative breast carcinoma", "Cancer", - "gastric cancer", "Cancer", - "colorectal cancer", "Cancer", - "colon sessile serrated adenoma/polyp", "Cancer", - "follicular lymphoma", "Cancer", - "B-cell acute lymphoblastic leukemia", "Cancer", - "B-cell non-Hodgkin lymphoma", "Cancer", - "acute myeloid leukemia", "Cancer", - "acute promyelocytic leukemia", "Cancer", - "plasma cell myeloma", "Cancer", - "clear cell renal carcinoma", "Cancer", - "nonpapillary renal cell carcinoma", "Cancer", - "basal cell carcinoma", "Cancer", - "colorectal neoplasm", "Cancer", - "adenocarcinoma", "Cancer", - "chromophobe renal cell carcinoma", "Cancer", - "neuroendocrine carcinoma", "Cancer", - "lung large cell carcinoma", "Cancer", - "tongue cancer", "Cancer", - "Wilms tumor", "Cancer", - "pleomorphic carcinoma", "Cancer", - "blastoma", "Cancer", - - # Neurodegenerative and Neurological Disorders - "dementia", "Neurodegenerative and Neurological Disorders", - "Alzheimer disease", "Neurodegenerative and Neurological Disorders", - "Parkinson disease", "Neurodegenerative and Neurological Disorders", - "amyotrophic lateral sclerosis", "Neurodegenerative and Neurological Disorders", - "multiple sclerosis", "Neurodegenerative and Neurological Disorders", - "Down syndrome", "Neurodegenerative and Neurological Disorders", - "trisomy 18", "Neurodegenerative and Neurological Disorders", - "frontotemporal dementia", "Neurodegenerative and Neurological Disorders", - "temporal lobe epilepsy", "Neurodegenerative and Neurological Disorders", - "Lewy body dementia", "Neurodegenerative and Neurological Disorders", - "amyotrophic lateral sclerosis 26 with or without frontotemporal dementia", "Neurodegenerative and Neurological Disorders", - - # Respiratory Conditions - "pulmonary fibrosis", "Respiratory Conditions", - "respiratory system disorder", "Respiratory Conditions", - "chronic obstructive pulmonary disease", "Respiratory Conditions", - "cystic fibrosis", "Respiratory Conditions", - "interstitial lung disease", "Respiratory Conditions", - "hypersensitivity pneumonitis", "Respiratory Conditions", - "non-specific interstitial pneumonia", "Respiratory Conditions", - "aspiration pneumonia", "Respiratory Conditions", - "pulmonary emphysema", "Respiratory Conditions", - "pulmonary sarcoidosis", "Respiratory Conditions", - - # Cardiovascular Diseases - "myocardial infarction", "Cardiovascular Diseases", - "acute myocardial infarction", "Cardiovascular Diseases", - "dilated cardiomyopathy", "Cardiovascular Diseases", - "heart failure", "Cardiovascular Diseases", - "arrhythmogenic right ventricular cardiomyopathy", "Cardiovascular Diseases", - "congenital heart disease", "Cardiovascular Diseases", - "non-compaction cardiomyopathy", "Cardiovascular Diseases", - "cardiomyopathy", "Cardiovascular Diseases", - "heart disorder", "Cardiovascular Diseases", - - # Metabolic and Other Disorders - "type 2 diabetes mellitus", "Metabolic and Other Disorders", - "chronic kidney disease", "Metabolic and Other Disorders", - "digestive system disorder", "Metabolic and Other Disorders", - "primary sclerosing cholangitis", "Metabolic and Other Disorders", - "gastritis", "Metabolic and Other Disorders", - "acute kidney failure", "Metabolic and Other Disorders", - "tubular adenoma", "Metabolic and Other Disorders", - "benign prostatic hyperplasia", "Metabolic and Other Disorders", - "opiate dependence", "Metabolic and Other Disorders", - "gingivitis", "Metabolic and Other Disorders", - "hyperplastic polyp", "Metabolic and Other Disorders", - "clonal hematopoiesis", "Metabolic and Other Disorders", - "epilepsy", "Metabolic and Other Disorders", - "age related macular degeneration 7", "Metabolic and Other Disorders", - "kidney benign neoplasm", "Metabolic and Other Disorders", - "malignant pancreatic neoplasm", "Metabolic and Other Disorders", - "cataract", "Metabolic and Other Disorders", - "macular degeneration", "Metabolic and Other Disorders", - "hydrosalpinx", "Metabolic and Other Disorders", - "tubulovillous adenoma", "Metabolic and Other Disorders", - "gastric intestinal metaplasia", "Metabolic and Other Disorders", - "Barrett esophagus", "Metabolic and Other Disorders", - - # Other Diseases - "injury", "Other Diseases", - "anencephaly", "Other Diseases", - "primary biliary cholangitis", "Other Diseases", - "keloid", "Other Diseases", - "kidney oncocytoma", "Other Diseases", - "respiratory failure", "Other Diseases", - "pilocytic astrocytoma", "Other Diseases" - ) - - temp_path = tempdir() - system(glue("~/bin/rclone copy box_adelaide:/minh_immune_map_disease/disease_data_grouped_further.csv {temp_path}/")) - - disease_data_grouped = - disease_data_grouped |> - left_join( - read_csv(glue("{temp_path}/disease_data_grouped_further.csv")) |> - dplyr::rename(disease_groups_further = disease_groups) - ) |> - mutate(disease_groups = if_else(!disease_groups_further |> is.na(), disease_groups_further, disease_groups)) |> - select(disease, disease_groups) - - tbl |> - - # TISSUE - select(-any_of("tissue_groups")) |> - left_join(tissue_grouped, copy=TRUE) |> - - # TECH - left_join(assay_data_grouped, copy=TRUE) |> - - # DISEASE - left_join(disease_data_grouped, copy=TRUE) |> - - # make disease tissue specific, omit Normal - mutate(disease_groups = paste(disease_groups, tissue_groups, sep = "_")) |> - mutate(disease_groups = if_else(disease_groups |> str_detect("Normal_.+"), "Normal", disease_groups)) |> - - - # TEMPORARY. de-group pancreas and liver - mutate(tissue_groups = case_when( - - tissue %in% c("gallbladder") ~ "gallbladder", - tissue %in% c("pancreas", "exocrine pancreas") ~ "pancreas", - tissue %in% c("liver", "caudate lobe of liver", "hepatic cecum" ) ~ "liver", - TRUE ~ tissue_groups - )) |> - - # SEX edit - mutate(sex = if_else(sex |> is.na(), "unknown", sex)) |> - - # Age - filter(age_days > 365) |> - mutate(age_years = age_days / 365) |> - mutate(age_bin = dplyr::case_when( - age_years < 3 ~ "Infancy", - age_years < 12 ~ "Childhood", - age_years < 20 ~ "Adolescence", - age_years < 40 ~ "Young Adulthood", - age_years < 50 ~ "Middle Age", - age_years < 60 ~ "Senior_50", - age_years < 70 ~ "Senior_60", - age_years >= 70 ~ "Senior_70", - TRUE ~ NA_character_ - )) |> - mutate(age_decade = ceiling(age_years/10) |> as.integer() |> as.character()) |> - - # left_join(age_bin_table, copy=TRUE) |> - - # ETHNICITY - left_join(ethnicity_grouped, copy=TRUE) |> - - dplyr::select( - sample_id, donor_id, dataset_id, title, collection_id, age_days, age_bin, age_decade, sex, - ethnicity_groups, tissue_groups, tissue, assay_groups, cell_type_unified_ensemble, - cell_type, disease_groups, is_immune - ) |> - as_tibble() |> - - # Center based on adolescence - mutate(age_days_scaled = age_days |> scale(center = 50*365) |> as.numeric()) - -} + + assay_data_grouped <- tribble( + ~assay, ~assay_groups, + "10x 3' v2", "10x Genomics 3", + "10x 3' v3", "10x Genomics 3", + "10x 5' v2", "10x Genomics 5", + "10x 5' v1", "10x Genomics 5", + "MARS-seq", "Plate based Technologies", + "10x 3' transcription profiling", "10x Genomics 3", + "10x 5' transcription profiling", "10x Genomics 5", + "Smart-seq2", "Smart seq", + "microwell-seq", "Microwell Technologies", + "TruDrop", "TruDrop", + "Drop-seq", "Drop based Technologies", + "Seq-Well S3", "Microwell Technologies", + "GEXSCOPE technology", "Other Technologies", + "Seq-Well", "Microwell Technologies", + "sci-RNA-seq", "Other Technologies", + "10x 3' v1", "10x Genomics 3", + "BD Rhapsody Whole Transcriptome Analysis", "Other Technologies", + "BD Rhapsody Targeted mRNA", "Other Technologies", + "CEL-seq2", "Plate based Technologies", + "SPLiT-seq", "Other Technologies", + "STRT-seq", "Plate based Technologies", + "inDrop", "Drop based Technologies", + "Smart-seq v4", "Smart seq", + "ScaleBio single cell RNA sequencing", "Other Technologies" + ) + + + disease_data_grouped <- tribble( + ~disease, ~disease_groups, + + # Normal control + "normal", "Normal", + + # Isolated Diseases + "COVID-19", "COVID-19 related", + "post-COVID-19 disorder", "COVID-19 related", + "long COVID-19", "COVID-19 related", + "glioblastoma", "Glioblastoma", + "lung adenocarcinoma", "Lung Adenocarcinoma", + "systemic lupus erythematosus", "Systemic Lupus Erythematosus", + + # Infectious and Immune-related Diseases (other than COVID-19) + "Crohn disease", "Infectious and Immune-related Diseases", + "Crohn ileitis", "Infectious and Immune-related Diseases", + "pneumonia", "Infectious and Immune-related Diseases", + "common variable immunodeficiency", "Infectious and Immune-related Diseases", + "toxoplasmosis", "Infectious and Immune-related Diseases", + "Plasmodium malariae malaria", "Infectious and Immune-related Diseases", + "type 1 diabetes mellitus", "Infectious and Immune-related Diseases", + "influenza", "Infectious and Immune-related Diseases", + "chronic rhinitis", "Infectious and Immune-related Diseases", + "periodontitis", "Infectious and Immune-related Diseases", + "localized scleroderma", "Infectious and Immune-related Diseases", + "lymphangioleiomyomatosis", "Infectious and Immune-related Diseases", + "listeriosis", "Infectious and Immune-related Diseases", + + # Cancer (other than isolated cancers) + "squamous cell lung carcinoma", "Cancer", + "small cell lung carcinoma", "Cancer", + "non-small cell lung carcinoma", "Cancer", + "breast carcinoma", "Cancer", + "breast cancer", "Cancer", + "luminal B breast carcinoma", "Cancer", + "luminal A breast carcinoma", "Cancer", + "triple-negative breast carcinoma", "Cancer", + "gastric cancer", "Cancer", + "colorectal cancer", "Cancer", + "colon sessile serrated adenoma/polyp", "Cancer", + "follicular lymphoma", "Cancer", + "B-cell acute lymphoblastic leukemia", "Cancer", + "B-cell non-Hodgkin lymphoma", "Cancer", + "acute myeloid leukemia", "Cancer", + "acute promyelocytic leukemia", "Cancer", + "plasma cell myeloma", "Cancer", + "clear cell renal carcinoma", "Cancer", + "nonpapillary renal cell carcinoma", "Cancer", + "basal cell carcinoma", "Cancer", + "colorectal neoplasm", "Cancer", + "adenocarcinoma", "Cancer", + "chromophobe renal cell carcinoma", "Cancer", + "neuroendocrine carcinoma", "Cancer", + "lung large cell carcinoma", "Cancer", + "tongue cancer", "Cancer", + "Wilms tumor", "Cancer", + "pleomorphic carcinoma", "Cancer", + "blastoma", "Cancer", + + # Neurodegenerative and Neurological Disorders + "dementia", "Neurodegenerative and Neurological Disorders", + "Alzheimer disease", "Neurodegenerative and Neurological Disorders", + "Parkinson disease", "Neurodegenerative and Neurological Disorders", + "amyotrophic lateral sclerosis", "Neurodegenerative and Neurological Disorders", + "multiple sclerosis", "Neurodegenerative and Neurological Disorders", + "Down syndrome", "Neurodegenerative and Neurological Disorders", + "trisomy 18", "Neurodegenerative and Neurological Disorders", + "frontotemporal dementia", "Neurodegenerative and Neurological Disorders", + "temporal lobe epilepsy", "Neurodegenerative and Neurological Disorders", + "Lewy body dementia", "Neurodegenerative and Neurological Disorders", + "amyotrophic lateral sclerosis 26 with or without frontotemporal dementia", "Neurodegenerative and Neurological Disorders", + + # Respiratory Conditions + "pulmonary fibrosis", "Respiratory Conditions", + "respiratory system disorder", "Respiratory Conditions", + "chronic obstructive pulmonary disease", "Respiratory Conditions", + "cystic fibrosis", "Respiratory Conditions", + "interstitial lung disease", "Respiratory Conditions", + "hypersensitivity pneumonitis", "Respiratory Conditions", + "non-specific interstitial pneumonia", "Respiratory Conditions", + "aspiration pneumonia", "Respiratory Conditions", + "pulmonary emphysema", "Respiratory Conditions", + "pulmonary sarcoidosis", "Respiratory Conditions", + + # Cardiovascular Diseases + "myocardial infarction", "Cardiovascular Diseases", + "acute myocardial infarction", "Cardiovascular Diseases", + "dilated cardiomyopathy", "Cardiovascular Diseases", + "heart failure", "Cardiovascular Diseases", + "arrhythmogenic right ventricular cardiomyopathy", "Cardiovascular Diseases", + "congenital heart disease", "Cardiovascular Diseases", + "non-compaction cardiomyopathy", "Cardiovascular Diseases", + "cardiomyopathy", "Cardiovascular Diseases", + "heart disorder", "Cardiovascular Diseases", + + # Metabolic and Other Disorders + "type 2 diabetes mellitus", "Metabolic and Other Disorders", + "chronic kidney disease", "Metabolic and Other Disorders", + "digestive system disorder", "Metabolic and Other Disorders", + "primary sclerosing cholangitis", "Metabolic and Other Disorders", + "gastritis", "Metabolic and Other Disorders", + "acute kidney failure", "Metabolic and Other Disorders", + "tubular adenoma", "Metabolic and Other Disorders", + "benign prostatic hyperplasia", "Metabolic and Other Disorders", + "opiate dependence", "Metabolic and Other Disorders", + "gingivitis", "Metabolic and Other Disorders", + "hyperplastic polyp", "Metabolic and Other Disorders", + "clonal hematopoiesis", "Metabolic and Other Disorders", + "epilepsy", "Metabolic and Other Disorders", + "age related macular degeneration 7", "Metabolic and Other Disorders", + "kidney benign neoplasm", "Metabolic and Other Disorders", + "malignant pancreatic neoplasm", "Metabolic and Other Disorders", + "cataract", "Metabolic and Other Disorders", + "macular degeneration", "Metabolic and Other Disorders", + "hydrosalpinx", "Metabolic and Other Disorders", + "tubulovillous adenoma", "Metabolic and Other Disorders", + "gastric intestinal metaplasia", "Metabolic and Other Disorders", + "Barrett esophagus", "Metabolic and Other Disorders", + + # Other Diseases + "injury", "Other Diseases", + "anencephaly", "Other Diseases", + "primary biliary cholangitis", "Other Diseases", + "keloid", "Other Diseases", + "kidney oncocytoma", "Other Diseases", + "respiratory failure", "Other Diseases", + "pilocytic astrocytoma", "Other Diseases" + ) + + disease_data_grouped = + disease_data_grouped |> + left_join( + readr::read_csv("/home/zhanchen/From_scratch/Mangiola_ImmuneAtlas/disease_data_grouped_further.csv") |> + rename(disease_groups_further = disease_groups) + ) |> + mutate(disease_groups = if_else(!disease_groups_further |> is.na(), disease_groups_further, disease_groups)) |> + select(disease, disease_groups) + + tbl |> + + # TISSUE + select(-any_of("tissue_groups")) |> + left_join(tissue_grouped, copy=TRUE) |> + + # TECH + left_join(assay_data_grouped, copy=TRUE) |> + + # DISEASE + left_join(disease_data_grouped, copy=TRUE) |> + + # make disease tissue specific, omit Normal + mutate(disease_groups = paste(disease_groups, tissue_groups, sep = "_")) |> + mutate(disease_groups = if_else(disease_groups |> str_detect("Normal_.+"), "Normal", disease_groups)) |> + + + # TEMPORARY. de-group pancreas and liver + mutate(tissue_groups = case_when( + + tissue %in% c("gallbladder") ~ "gallbladder", + tissue %in% c("pancreas", "exocrine pancreas") ~ "pancreas", + tissue %in% c("liver", "caudate lobe of liver", "hepatic cecum" ) ~ "liver", + TRUE ~ tissue_groups + )) |> + + # SEX edit + mutate(sex = if_else(sex |> is.na(), "unknown", sex)) |> + + # Age + filter(age_days > 365) |> + mutate(age_years = age_days / 365) |> + mutate(age_bin = dplyr::case_when( + age_years < 3 ~ "Infancy", + age_years < 12 ~ "Childhood", + age_years < 20 ~ "Adolescence", + age_years < 40 ~ "Young Adulthood", + age_years < 50 ~ "Middle Age", + age_years < 60 ~ "Senior_50", + age_years < 70 ~ "Senior_60", + age_years >= 70 ~ "Senior_70", + TRUE ~ NA_character_ + )) |> + mutate(age_decade = ceiling(age_years/10) |> as.integer() |> as.character()) |> + + # left_join(age_bin_table, copy=TRUE) |> + + # ETHNICITY + left_join(ethnicity_grouped, copy=TRUE) |> + + dplyr::select( + sample_id, donor_id, dataset_id, title, collection_id, age_days, age_bin, age_decade, sex, + ethnicity_groups, tissue_groups, tissue, assay_groups, cell_type_unified_ensemble, + cell_type, disease_groups + ) |> + as_tibble() |> + + # Center based on adolescence + mutate(age_days_scaled = age_days |> scale(center = 50*365) |> as.numeric()) + + } \ No newline at end of file From 1288796bdc5bda07f3193b6331b7963b27dd50ad Mon Sep 17 00:00:00 2001 From: Chen Zhan Date: Fri, 18 Jul 2025 20:54:37 +0800 Subject: [PATCH 12/12] dplyr::rename --- rebuttal_CellPress/edit_covariates.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rebuttal_CellPress/edit_covariates.R b/rebuttal_CellPress/edit_covariates.R index 8170092..0b30213 100644 --- a/rebuttal_CellPress/edit_covariates.R +++ b/rebuttal_CellPress/edit_covariates.R @@ -243,7 +243,7 @@ edit_covariates = function(tbl){ enframe(name ="tissue_groups") |> distinct() |> unnest(value) |> - rename(tissue = value) + dplyr::rename(tissue = value) ethnicity_grouped <- tribble( ~self_reported_ethnicity, ~ethnicity_groups, @@ -443,7 +443,7 @@ edit_covariates = function(tbl){ disease_data_grouped |> left_join( readr::read_csv("/home/zhanchen/From_scratch/Mangiola_ImmuneAtlas/disease_data_grouped_further.csv") |> - rename(disease_groups_further = disease_groups) + dplyr::rename(disease_groups_further = disease_groups) ) |> mutate(disease_groups = if_else(!disease_groups_further |> is.na(), disease_groups_further, disease_groups)) |> select(disease, disease_groups)