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..7bc3d70
--- /dev/null
+++ b/rebuttal_CellPress/PAWSEY_sccomp_tar_scripts.R
@@ -0,0 +1,926 @@
+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, "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.")
+  #   }
+  # }
+  
+  #-----------------------#
+  # 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, 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 = "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, 
+          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
diff --git a/rebuttal_CellPress/edit_covariates.R b/rebuttal_CellPress/edit_covariates.R
new file mode 100644
index 0000000..0b30213
--- /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) |>
+      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"
+    )
+
+    disease_data_grouped =
+      disease_data_grouped |>
+      left_join(
+        readr::read_csv("/home/zhanchen/From_scratch/Mangiola_ImmuneAtlas/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
+      ) |>
+      as_tibble() |>
+
+      # Center based on adolescence
+      mutate(age_days_scaled = age_days  |> scale(center = 50*365) |> as.numeric())
+
+  }
\ No newline at end of file
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..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
@@ -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({
@@ -11,11 +28,16 @@ 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"
+  ## 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_cell_type"
+  target_cell_type = "cd4 fh em"
+
   tar_option_set(
     
     
@@ -41,10 +63,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),
+          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
         )
       ),
@@ -55,10 +79,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 +95,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
         )
       )
@@ -84,7 +112,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,12 +267,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,6 +297,551 @@ get_adjusted_matrix = function(summary_df, column_adjusted){
     
     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())
+
+  }
+
+  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
@@ -292,30 +865,74 @@ get_adjusted_matrix = function(summary_df, column_adjusted){
     #   
     # ),
     
-    # 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.
+        # reference_sample ------
+    # calculate reference_sample for scalig gene counts
+    # this ensures calculation only done once
     tar_target(
-      # pseudobulk_sample ------
-      pseudobulk_sample,
+      reference_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/")
         
-        # TEMPORARY BECAUSE I FORGOT TO INTEGRATE AGE BINS
-        se = se |> 
+        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(
-            readRDS(metadata_path) |> 
-              distinct(sample_id,  age_days, age_bin) 
-          )
+            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 ]
@@ -345,7 +962,113 @@ get_adjusted_matrix = function(summary_df, column_adjusted){
             which.min()                             # Identify the smallest difference
         ]
         
-        message('TAR: pseudobulk_sample Phase2')
+        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,
+      {
+        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
+        # ]
+        
         se = 
           se |> 
           keep_abundant(design = 
@@ -364,8 +1087,9 @@ get_adjusted_matrix = function(summary_df, column_adjusted){
           ) |> 
           
           # 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),
@@ -377,7 +1101,7 @@ get_adjusted_matrix = function(summary_df, column_adjusted){
           # 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)
+          # library size factor is the reciproque of the multiplier (correction factor)
           mutate(offset = log(1/multiplier)) |> 
           
           # Set intercept
@@ -385,7 +1109,8 @@ get_adjusted_matrix = function(summary_df, column_adjusted){
             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, "Adolescence")
+            age_bin = fct_relevel(age_bin, "Senior_50")
           ) 
         
         # # Add dispersion
@@ -395,16 +1120,17 @@ get_adjusted_matrix = function(summary_df, column_adjusted){
         #   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")),
+      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"
     ),
-
+    
     # pseudobulk_sample_id ------
     # This target extracts unique sample ids from the pseudobulk sample  
     tar_target(
@@ -422,12 +1148,8 @@ get_adjusted_matrix = function(summary_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"
-        # )) |>
-        # head(1) %>% 
+        # filter(.feature %in% readRDS('/hpcfs/groups/phoenix-hpc-mangiola_laboratory/Mangiola_ImmuneAtlas/ning_data/ethnicity_umap_selected_genes.rds')) |>
+        # slice_sample(n=1500) %>% 
         group_by(.feature) |> 
         tar_group(), 
       iteration = "group",
@@ -519,20 +1241,27 @@ get_adjusted_matrix = function(summary_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(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)
+        
+        # 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,189 +1272,296 @@ get_adjusted_matrix = function(summary_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))), 1.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, 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)
+
             # 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,
           family = zero_inflated_negbinomial(),
           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",
+          # sample_prior = TRUE, 
           init = inits,
-          iter = 400  # Increase iterations for better convergence
+          iter = 600  # Increase iterations for better convergence
         )
         
       })) |> 
         
-      # 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")
-      
+      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.
     tar_target(
       summary,
       estimates_chunk |>
-        mutate(summary = 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"
-   #        ),
-
-      # Median instead and mad of mean and sd
-      robust=TRUE,
-      alpha = 0.1
-      )
-        )) |>
-        
-      mutate(
+    #     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"
+    #       #        ),
+          
+    #       # 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(
             brms_fit,  function(x) {
@@ -754,28 +1590,29 @@ get_adjusted_matrix = function(summary_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"))
+      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,11 +1700,29 @@ 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_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") ,