step2_batch_effect_correction_withChallengeData.Rmd

---
title: "CMI-PB Team's Data Processing Approach for 3rd (Public) Challenge Dataset: All data batch-effect correction"
#date: "Last update: `r format(Sys.time(), '%d %B, %Y')`" 
output:
  BiocStyle::html_document:
    toc: true
    toc_float:
        collapsed: true
        smooth_scroll: true
    toc_depth: 3
    fig_caption: yes
    code_folding: show
    number_sections: false
date: "2024-08-20"
---

The training dataset for the Public challenge comprises three multi-omics datasets (designated as 2020, 2021 and 2022) and challenge dataset (2023) that require processing and normalization to generate computable matrices suitable for subsequent model development. While the procedures for data processing and normalization are inherently user-specific, the CMI-PB team has devised a data processing method, drawing inspiration from the approach used in the 2nd CMI-PB challenge. The codebase is also available on GitHub at [<https://github.com/CMI-PB/cmi-pb-3rd-public-challenge-data-prep>]. If you have specific questions, please contact us via Solutions center.

## Download and read Public challenge data from CMI-PB website

The data files for the Public CMI-PB challenge can be accessed at [<https://www.cmi-pb.org/downloads/cmipb_challenge_datasets/current/3rd_challenge/raw_datasets/>]. They are available for direct file download or as R data objects. For our notebook, we chose to download the data as R data objects. These objects contain both demographical metadata of subjects and experimental data, including:

-   **Plasma antibody titers**: Plasma antibodies against Tdap were measured at all time points using Luminex assay,
-   **Plasma cytokine concentrations by Olink**: Plasma cytokine concentration analysis using OLINK assay,
-   **Plasma cytokine concentrations by Legendplex**: Plasma cytokine concentration analysis using Legendplex assay,
-   **PBMC gene expression:** Gene expression analysis (RNAseq) of bulk peripheral blood mononuclear cells (PBMCs),
-   **PBMC cell frequency:** Cell frequency analysis of PBMC subsets,
-   **t cell activation:** T cell polarization using FluoroSpot assay,
-   **t cell polarization:** T cell activation using AIM assay

were performed before and after booster vaccination.

```{r message=FALSE, warning=FALSE, paged.print=TRUE}

base_dir = "/home/pshinde/repos/cmi-pb/public_challenge/training_preprocess/"
#base_dir = "../"

dir_raw_training <- paste0(base_dir, "data/raw_training_dataset/")
dir_raw_prediction <- paste0(base_dir, "data/raw_challenge_dataset/")
dir_RDS_objects <- paste0(base_dir, "data/output/")

## `codebase.R` installs required packages and all house keeping functions
source(paste0(base_dir, "scripts/codebase.R"))

master_database_data <- readRDS(paste0(dir_RDS_objects, "master_harmonized_data_v20240825.RDS"))

#training_dataset <- subset_dataset(master_database_data$training, c("2020_dataset", "2021_dataset", "2022_dataset"))
training_dataset <- master_database_data$training
challenge_dataset <- master_database_data$challenge

#subject_specimen <- master_database_data$subject_specimen  %>% 
#  mutate(timepoint = planned_day_relative_to_boost)

training_subject_specimen <- training_dataset$subject_specimen %>%
  select(c("specimen_id","subject_id","dataset","timepoint","infancy_vac","biological_sex","date_of_boost"))

challenge_subject_specimen <- challenge_dataset$subject_specimen %>%
  select(c("specimen_id","subject_id","dataset","timepoint","infancy_vac","biological_sex","date_of_boost"))

subject_specimen = training_subject_specimen %>%
  rbind(challenge_subject_specimen)

gene_90_38_export <- read_tsv(paste0(base_dir, "data/gene_90_38_export.tsv"))

mito_genes <- gene_90_38_export %>%
  filter(substr(display_label, 1,3) == "MT-")

gene_90_38_shortlist <- gene_90_38_export %>%
  filter(biotype == "protein_coding") %>%
  filter(!versioned_ensembl_gene_id %in% mito_genes$versioned_ensembl_gene_id)

batch.factors = c("timepoint","infancy_vac","biological_sex","dataset")
data_obj = training_dataset

challenge_subject_specimen_baseline <- subject_specimen %>% 
  filter(dataset %in% c("2023_dataset")) %>% 
  filter(timepoint <= 40)

subject_specimen_baseline <- subject_specimen %>% 
  filter(timepoint <= 40)

```
## t_cell_polarization (Flurospot) Analysis
```{r}
batch.factors = c("timepoint","dataset")

tcellpol_wide_before_wide <- training_dataset$t_cell_polarization$wide %>%
  rbind(challenge_dataset$t_cell_polarization$wide) %>%
  select(-c("DMSO_P01579", "DMSO_Q16552", "DMSO_P05113")) 

# Calculate the percentage of missing values for each column
missing_percent <- colMeans(is.na(tcellpol_wide_before_wide)) * 100
# Identify columns with more than 80% missing values
columns_with_high_na <- which(missing_percent > 80)

tcellpol_wide_before_wide = tcellpol_wide_before_wide %>%
  filter(!specimen_id %in% c(1159)) %>%
  column_to_rownames("specimen_id")%>%
  t() 

tcellpol_wide_before_wide_imputed = tcellpol_wide_before_wide[rowMeans(is.na(tcellpol_wide_before_wide)) < 1, ] %>%
    as.matrix() %>%
    impute.knn() %>%
    .$data

pvca_analysis(tcellpol_wide_before_wide_imputed, subject_specimen, batch.factors, plot_title = "T cell polarization:  Raw data")


tcell_polarization_processed_data = list(
      
      raw_data = tcellpol_wide_before_wide,
      normalized_data =  "Note: normalized_data data was not included, raw_data shows not batch effects. We suggest using raw_data data for t_cell_polarization assay",
      batchCorrected_data =  "Note: Batch-corrected data was not included, raw_data shows not batch effects. We suggest using raw_data data for t_cell_polarization assay"
    )

```
## t_cell_activation (AIM) Analysis
```{r}
batch.factors = c("timepoint","dataset")

tcell_activation_wide_before_wide <- training_dataset$t_cell_activation$wide %>%
  rbind(challenge_dataset$t_cell_activation$wide) %>%
  filter(specimen_id %in% subject_specimen$specimen_id) %>%
  select(-c("DMSO")) %>%
  column_to_rownames("specimen_id")%>%
  t() 

# Calculate the percentage of missing values for each column
missing_percent <- colMeans(is.na(tcell_activation_wide_before_wide)) * 100
# Identify columns with more than 80% missing values
columns_with_high_na <- which(missing_percent > 80)

pvca_analysis(tcell_activation_wide_before_wide, subject_specimen, batch.factors, plot_title = "T cell Activation:  Raw data")


tcell_activation_processed_data = list(
      
      raw_data = tcell_activation_wide_before_wide,
      normalized_data =  "Note: normalized_data data was not included, raw_data show no dataset/batch effects. We suggest using raw_data data for t_cell_activation assay",
      batchCorrected_data =  "Note: Batch-corrected data was not included, raw_data show no dataset/batch effects. We suggest using raw_data data for t_cell_activation assay"
    )

```
## Antibody titers
```{r}
abtiter_wide_before <- training_dataset$plasma_antibody_levels$wide %>%
  rbind(challenge_dataset$plasma_antibody_levels$wide) %>%
  filter(specimen_id %in% subject_specimen$specimen_id)  %>%
  column_to_rownames("specimen_id")%>%
  t() 

abtiter_wide_before_long <- training_dataset$plasma_antibody_levels$long %>%
  rbind(challenge_dataset$plasma_antibody_levels$long) 

pvca_analysis(abtiter_wide_before, subject_specimen, batch.factors, plot_title = "Plasma Antibody titer:  Raw data")

ab_data_obj = list(
  plasma_antibody_levels = list(
    long = abtiter_wide_before_long,
    wide = abtiter_wide_before
  ),
  subject_specimen = subject_specimen
)


## Apply data normalization and batch correction
abtiter_data_processed = processAbtiter(ab_data_obj, BatchCorrection = TRUE)

pvca_analysis(abtiter_data_processed$normalized_data, subject_specimen, batch.factors, plot_title = "Plasma Antibody titer: Normalization")
pvca_analysis(abtiter_data_processed$batchCorrected_data, subject_specimen, batch.factors, plot_title = "Plasma Antibody titer:  Normalization and batch effect correction")

abtiter_processed_data = list(
      
      #metadata = rnaseq_metaData,
      raw_data = abtiter_data_processed$raw_data,
      normalized_data = abtiter_data_processed$normalized_data,
      batchCorrected_data = abtiter_data_processed$batchCorrected_data
    )
```
## Cell frequency Analysis
```{r}
## Before normalization
cell_wide_before_wide <-  training_dataset$pbmc_cell_frequency$wide %>%
  filter(specimen_id %in% subject_specimen$specimen_id) %>%
  rbind(challenge_dataset$pbmc_cell_frequency$wide) %>%
  distinct() %>%
  column_to_rownames("specimen_id")%>%
  t() 

pvca_analysis(cell_wide_before_wide, subject_specimen, batch.factors, plot_title = "PBMC Cell frequency:  Raw data")

## Normalized
cell_wide_before_long <- cell_wide_before_wide %>%
  as.data.frame()  %>%
  rownames_to_column("cell_type_name") %>%
  pivot_longer(!cell_type_name, names_to = "specimen_id", values_to = "percent_live_cell") %>%
  mutate(specimen_id = as.numeric(specimen_id))
  
cell_data_obj = list(
  
  pbmc_cell_frequency = list(
    
    long = cell_wide_before_long,
    wide = cell_wide_before_wide
  ),
  subject_specimen = subject_specimen %>%
    mutate(planned_day_relative_to_boost = timepoint)
)

## Assemble data
count_data_long = cell_wide_before_long
df_subject_specimen = subject_specimen %>%
    mutate(planned_day_relative_to_boost = timepoint)
  
## Perform median normalization
cytof_median_D0 <- count_data_long %>%
    left_join(df_subject_specimen[c("specimen_id", "dataset")]) %>%
    filter(specimen_id %in% unique(df_subject_specimen[df_subject_specimen$planned_day_relative_to_boost == 0,]$specimen_id)) %>%
    group_by(dataset, cell_type_name)  %>%
    summarise(median = median(percent_live_cell, na.rm = T))
  
cell_long_normalized_pre <-  count_data_long  %>%
    left_join(df_subject_specimen[c("specimen_id", "dataset")]) %>%
    left_join(cytof_median_D0) %>%
    mutate(percent_live_cell_normalized = if_else(is.na(percent_live_cell) == T, NA, percent_live_cell/median))
  
## Reshape dataframe in wide format
cell_wide_normalized_pre <- cell_long_normalized_pre  %>%
    dplyr::select(cell_type_name, specimen_id, percent_live_cell_normalized) %>%
    pivot_wider(names_from = "cell_type_name", values_from = percent_live_cell_normalized) %>%
    column_to_rownames("specimen_id")%>%
    t() 
  
cellFreq_normalized_imputed = cell_wide_normalized_pre[rowMeans(is.na(cell_wide_normalized_pre)) < 1, ] %>%
    as.matrix() %>%
    impute.knn() %>%
    .$data

pvca_analysis(cellFreq_normalized_imputed, subject_specimen, batch.factors, plot_title = "PBMC Cell frequency:  Normalized data")
  
## Batch correction
batch_lebels = as.data.frame(colnames(cellFreq_normalized_imputed)) %>%
      rename(specimen_id = starts_with("colnames")) %>%
      mutate(specimen_id = as.integer(specimen_id)) %>%
      left_join(subject_specimen) %>%
      dplyr::select(dataset)
  
 cellFreq_batchCorrected = ComBat(cellFreq_normalized_imputed, batch = batch_lebels$dataset)
 cellFreq_batchCorrected_imputed = cellFreq_batchCorrected[rowMeans(is.na(cellFreq_batchCorrected)) < 1, ] %>%
      as.matrix() %>%
      impute.knn() %>%
      .$data

pvca_analysis(cellFreq_batchCorrected_imputed, subject_specimen, batch.factors, plot_title = "PBMC Cell frequency: Normalization and batch effect correction")


cell_frequency_processed_data = list(
      
      raw_data = cell_wide_before_wide,
      normalized_data = cellFreq_normalized_imputed,
      batchCorrected_data = cellFreq_batchCorrected_imputed
    )

```
## Cytokine concetrations By Olink Analysis
```{r}
## Before batch correction
olink_wide_before <- training_dataset$plasma_cytokine_concentrations_by_olink$wide  %>%
  rbind(challenge_dataset$plasma_cytokine_concentrations_by_olink$wide) %>%
  filter(specimen_id %in% subject_specimen$specimen_id) %>%
  column_to_rownames("specimen_id")%>%
  t() 

pvca_analysis(olink_wide_before, subject_specimen, batch.factors, plot_title = "Cytokine concetrations By Olink assay:  Raw data")

olink_wide_before_long <- olink_wide_before %>%
  as.data.frame()  %>%
  rownames_to_column("protein_id") %>%
  pivot_longer(!protein_id, names_to = "specimen_id", values_to = "concentration") %>%
  mutate(specimen_id = as.numeric(specimen_id))

## Assemble data
count_data_long_olink = olink_wide_before_long
df_subject_specimen = subject_specimen  %>%
    mutate(planned_day_relative_to_boost = timepoint)
  
## Perform median normalization
cytokine_median_D0_olink <- count_data_long_olink %>%
    left_join(df_subject_specimen[c("specimen_id", "dataset")]) %>%
    filter(specimen_id %in% unique(df_subject_specimen[df_subject_specimen$planned_day_relative_to_boost == 0,]$specimen_id)) %>%
    group_by(dataset, protein_id)  %>%
    summarise(median = median(concentration, na.rm = T))
  
cytokine_long_normalized_pre_olink <-  count_data_long_olink  %>%
    left_join(df_subject_specimen[c("specimen_id", "dataset")]) %>%
    left_join(cytokine_median_D0_olink) %>%
    mutate(concentration_normalized = if_else(is.na(concentration) == T, NA, concentration/median))
  
## Reshape dataframe in wide format
cytokine_wide_normalized_pre_olink <- cytokine_long_normalized_pre_olink  %>%
    dplyr::select(protein_id, specimen_id, concentration_normalized) %>%
    pivot_wider(names_from = "protein_id", values_from = concentration_normalized) %>%
    column_to_rownames("specimen_id")%>%
    t() 
  
cytokineFreq_normalized_imputed_olink = cytokine_wide_normalized_pre_olink[rowMeans(is.na(cytokine_wide_normalized_pre_olink)) < 1, ] %>%
    as.matrix() %>%
    impute.knn() %>%
    .$data
  
pvca_analysis(cytokineFreq_normalized_imputed_olink, subject_specimen, batch.factors, plot_title = "Cytokine concetrations By OLINK: Normalization")

batch_lebels = as.data.frame(colnames(cytokineFreq_normalized_imputed_olink)) %>%
      rename(specimen_id = starts_with("colnames")) %>%
      mutate(specimen_id = as.integer(specimen_id)) %>%
      left_join(df_subject_specimen) %>%
      dplyr::select(dataset)
    
cytokineFreq_batchCorrected_olink = ComBat(cytokineFreq_normalized_imputed_olink, batch = batch_lebels$dataset)
    
cytokineFreq_batchCorrected_imputed_olink = cytokineFreq_batchCorrected_olink[rowMeans(is.na(cytokineFreq_batchCorrected_olink)) < 1, ] %>%
      as.matrix() %>%
      impute.knn() %>%
      .$data
    

pvca_analysis(cytokineFreq_batchCorrected_imputed_olink, subject_specimen, batch.factors, plot_title = "Cytokine concetrations By OLINK: Normalization and batch effect correction")


olink_processed_data = list(
      
      raw_data = olink_wide_before,
      normalized_data = cytokineFreq_normalized_imputed_olink,
      batchCorrected_data = cytokineFreq_batchCorrected_imputed_olink
    )

```


## plasma_cytokine_concentrations_by_legendplex Analysis
```{r}
#subject_specimen_legendplex = data_obj$subject_specimen %>%
#  filter(!dataset %in% c("2020_dataset"))

## Before batch correction
legendplex_wide_before <- training_dataset$plasma_cytokine_concentrations_by_legendplex$wide  %>%
  rbind(challenge_dataset$plasma_cytokine_concentrations_by_legendplex$wide)  %>%
  column_to_rownames("specimen_id")%>%
  t() 

pvca_analysis(legendplex_wide_before, subject_specimen, batch.factors, plot_title = "Cytokine concentrations By legendplex assay:  Raw data")

legendplex_wide_before_long <- legendplex_wide_before %>%
  as.data.frame()  %>%
  rownames_to_column("protein_id") %>%
  pivot_longer(!protein_id, names_to = "specimen_id", values_to = "concentration") %>%
  mutate(specimen_id = as.numeric(specimen_id))

## Assemble data
count_data_long = legendplex_wide_before_long
  df_subject_specimen = subject_specimen  %>%
    mutate(planned_day_relative_to_boost = timepoint)
  
## Perform median normalization
cytokine_median_D0 <- count_data_long %>%
    left_join(df_subject_specimen[c("specimen_id", "dataset")]) %>%
    filter(specimen_id %in% unique(df_subject_specimen[df_subject_specimen$planned_day_relative_to_boost == 0,]$specimen_id)) %>%
    group_by(dataset, protein_id)  %>%
    summarise(median = median(concentration, na.rm = T))
  
cytokine_long_normalized_pre <-  count_data_long  %>%
    left_join(df_subject_specimen[c("specimen_id", "dataset")]) %>%
    left_join(cytokine_median_D0) %>%
    mutate(concentration_normalized = if_else(is.na(concentration) == T, NA, concentration/median))
  
## Reshape dataframe in wide format
cytokine_wide_normalized_pre <- cytokine_long_normalized_pre  %>%
    dplyr::select(protein_id, specimen_id, concentration_normalized) %>%
    pivot_wider(names_from = "protein_id", values_from = concentration_normalized) %>%
    column_to_rownames("specimen_id")%>%
    t() 
  
cytokineFreq_normalized_imputed = cytokine_wide_normalized_pre[rowMeans(is.na(cytokine_wide_normalized_pre)) < 1, ] %>%
    as.matrix() %>%
    impute.knn() %>%
    .$data
  
pvca_analysis(cytokineFreq_normalized_imputed, subject_specimen, batch.factors, plot_title = "Cytokine concentrations By legendplex assay:  Raw data")

## Apply data normalization and batch correction
batch_lebels = as.data.frame(colnames(cytokineFreq_normalized_imputed)) %>%
      rename(specimen_id = starts_with("colnames")) %>%
      mutate(specimen_id = as.integer(specimen_id)) %>%
      left_join(df_subject_specimen) %>%
      dplyr::select(dataset)
    
cytokineFreq_batchCorrected = ComBat(cytokineFreq_normalized_imputed, batch = batch_lebels$dataset)
    
cytokineFreq_batchCorrected_imputed = cytokineFreq_batchCorrected[rowMeans(is.na(cytokineFreq_batchCorrected)) < 1, ] %>%
      as.matrix() %>%
      impute.knn() %>%
      .$data

pvca_analysis(cytokineFreq_batchCorrected_imputed, subject_specimen, batch.factors, plot_title = "Cytokine concetrations By legendplex: Normalization and batch effect correction")
## Note: Batch-corrected data was not included, as it resulted in over-correction and introduced unexpected variations.

legendplex_processed_data = list(
      
      raw_data = legendplex_wide_before,
      normalized_data = cytokineFreq_normalized_imputed,
      batchCorrected_data = "Note: Batch-corrected data was not included, as it resulted in over-correction and introduced unexpected variations. We suggest using normalized_data data for legendplex assay"
    )

```

## Gene expression data Analysis:: Raw Count
```{r}
rnaseq_countData <- training_dataset$pbmc_gene_expression$wide_raw_count %>%
  rbind(challenge_dataset$pbmc_gene_expression$wide_raw_count)  %>%
  column_to_rownames("specimen_id") %>%
  t()  %>%
  as.data.frame()  

colnames(rnaseq_countData) = as.integer(colnames(rnaseq_countData))

rnaseq_metaData <- subject_specimen %>%
  filter(specimen_id %in% colnames(rnaseq_countData)) %>%
  mutate(specimen_id1 = specimen_id) %>%
  column_to_rownames("specimen_id1") 

#' Identify genes with rawcount >=1 that are present in at least 80% of either the aP(wP) cohort and absent in at least 80% of wP(aP) cohort.
threshold_proportion_greater_than_1 = 0.8

rawcount_sum_infancy_subgroup <- rnaseq_countData %>%
  rownames_to_column("versioned_ensembl_gene_id") %>%
  filter(versioned_ensembl_gene_id %in% gene_90_38_shortlist$versioned_ensembl_gene_id) %>%
  pivot_longer(!versioned_ensembl_gene_id, values_to = "rawcount", names_to = "specimen_id") %>%
  mutate(specimen_id = as.integer(specimen_id)) %>%
  left_join(subject_specimen) %>%
  group_by(dataset, versioned_ensembl_gene_id, infancy_vac) %>%
  #group_by(versioned_ensembl_gene_id, infancy_vac) %>%
  summarise(proportion_greater_than_1 = mean(rawcount >= 1)) %>%
  pivot_wider(names_from = infancy_vac, values_from = proportion_greater_than_1)  %>%
  mutate(gene_meets_criterion_aP = aP >= threshold_proportion_greater_than_1 & wP <= (1 - threshold_proportion_greater_than_1),
         gene_meets_criterion_wP = wP >= threshold_proportion_greater_than_1 & aP <= (1 - threshold_proportion_greater_than_1)
         )  %>%
  filter((gene_meets_criterion_aP == TRUE & gene_meets_criterion_wP == FALSE) || (gene_meets_criterion_aP == FALSE & gene_meets_criterion_wP == TRUE))
  
#' Create a shortlist of genes (rawcount >= 1) in at least 30% of the specimens.
rawcount_shortlist <- rnaseq_countData %>%
  rownames_to_column("versioned_ensembl_gene_id") %>%
  filter(versioned_ensembl_gene_id %in% gene_90_38_shortlist$versioned_ensembl_gene_id) %>%
  pivot_longer(!versioned_ensembl_gene_id, values_to = "rawcount", names_to = "specimen_id") %>%
  mutate(specimen_id = as.integer(specimen_id)) %>%
  left_join(subject_specimen) %>%
  group_by(versioned_ensembl_gene_id) %>%
  #group_by(versioned_ensembl_gene_id, infancy_vac) %>%
  summarise(proportion = mean(rawcount >= 50))  %>%
  filter(proportion >= 0.3)


## Before batch correction
rnaseq_countData_v2 <- rnaseq_countData %>%
  rownames_to_column("versioned_ensembl_gene_id") %>%
  filter(versioned_ensembl_gene_id %in% gene_90_38_shortlist$versioned_ensembl_gene_id) %>%
  filter(!versioned_ensembl_gene_id %in% rawcount_sum_infancy_subgroup$versioned_ensembl_gene_id) %>%
  filter(versioned_ensembl_gene_id %in% rawcount_shortlist$versioned_ensembl_gene_id) %>%
  column_to_rownames("versioned_ensembl_gene_id")

mad_2020 <- mad_calculations(rnaseq_countData_v2, subject_specimen, c("2020_dataset"))
mad_2021 <- mad_calculations(rnaseq_countData_v2, subject_specimen, c("2021_dataset"))
mad_2022 <- mad_calculations(rnaseq_countData_v2, subject_specimen, c("2022_dataset"))
mad_2023 <- mad_calculations(rnaseq_countData_v2, subject_specimen, c("2023_dataset"))

mad_shotlisted_genes <- intersect(intersect(mad_2020$gene_id, mad_2021$gene_id), mad_2022$gene_id)
## mad_2023$gene_id genes were not included

rnaseq_countData_v3 <- rnaseq_countData_v2 %>%
  rownames_to_column("versioned_ensembl_gene_id") %>%
  filter(versioned_ensembl_gene_id %in% mad_shotlisted_genes) %>%
  column_to_rownames("versioned_ensembl_gene_id")
```


## RNASeq plot: Raw Count
```{r}
pvca_analysis_rnaseq(rnaseq_countData_v3, subject_specimen, batch.factors, plot_title = "RNASeq Raw count: Raw data")

batch_labels = as.data.frame(colnames(rnaseq_countData_v3)) %>%
      rename(specimen_id = starts_with("colnames")) %>%
      mutate(specimen_id = as.integer(specimen_id)) %>%
      left_join(rnaseq_metaData) %>%
      dplyr::select(dataset)
    
rnaseq_batchCorrected_rawcount = sva::ComBat_seq(as.matrix(rnaseq_countData_v3), batch = batch_labels$dataset)
    
pvca_analysis_rnaseq(rnaseq_batchCorrected_rawcount, subject_specimen, batch.factors, plot_title = "RNASeq Raw count: Batch correction")

rnaseq_normalised_data_rawcount = list(
      
      raw_data = as.matrix(rnaseq_countData_v3),
      batchCorrected_data = rnaseq_batchCorrected_rawcount
    )


```
## Gene expression data Analysis:: TPM Count
```{r}
rnaseq_tpmcountData <- training_dataset$pbmc_gene_expression$wide_tpm %>%
  rbind(challenge_dataset$pbmc_gene_expression$wide_tpm)  %>%
  column_to_rownames("specimen_id") %>%
  t()  %>%
  as.data.frame()  

colnames(rnaseq_tpmcountData) = as.integer(colnames(rnaseq_tpmcountData))

rnaseq_metaData <- subject_specimen %>%
  filter(specimen_id %in% colnames(rnaseq_tpmcountData)) %>%
  mutate(specimen_id1 = specimen_id) %>%
  column_to_rownames("specimen_id1") 

#' Identify genes with tpm >=1 that are present in at least 80% of either the aP(wP) cohort and absent in at least 80% of wP(aP) cohort.
threshold_proportion_greater_than_1 = 0.8

tpm_sum_infancy_subgroup <- rnaseq_tpmcountData %>%
  rownames_to_column("versioned_ensembl_gene_id") %>%
  filter(versioned_ensembl_gene_id %in% gene_90_38_shortlist$versioned_ensembl_gene_id) %>%
  pivot_longer(!versioned_ensembl_gene_id, values_to = "tpm", names_to = "specimen_id") %>%
  mutate(specimen_id = as.integer(specimen_id)) %>%
  left_join(subject_specimen) %>%
  group_by(dataset, versioned_ensembl_gene_id, infancy_vac) %>%
  #group_by(versioned_ensembl_gene_id, infancy_vac) %>%
  summarise(proportion_greater_than_1 = mean(tpm >= 1)) %>%
  pivot_wider(names_from = infancy_vac, values_from = proportion_greater_than_1)  %>%
  mutate(gene_meets_criterion_aP = aP >= threshold_proportion_greater_than_1 & wP <= (1 - threshold_proportion_greater_than_1),
         gene_meets_criterion_wP = wP >= threshold_proportion_greater_than_1 & aP <= (1 - threshold_proportion_greater_than_1)
         )  %>%
  filter((gene_meets_criterion_aP == TRUE & gene_meets_criterion_wP == FALSE) || (gene_meets_criterion_aP == FALSE & gene_meets_criterion_wP == TRUE))
  
#' Create a shortlist of genes (tpm >= 1) in at least 30% of the specimens.
tpm_shortlist <- rnaseq_tpmcountData %>%
  rownames_to_column("versioned_ensembl_gene_id") %>%
  filter(versioned_ensembl_gene_id %in% gene_90_38_shortlist$versioned_ensembl_gene_id) %>%
  pivot_longer(!versioned_ensembl_gene_id, values_to = "tpm", names_to = "specimen_id") %>%
  mutate(specimen_id = as.integer(specimen_id)) %>%
  left_join(subject_specimen) %>%
  group_by(versioned_ensembl_gene_id) %>%
  #group_by(versioned_ensembl_gene_id, infancy_vac) %>%
  summarise(proportion = mean(tpm >= 5))  %>%
  filter(proportion >= 0.3)


## Before batch correction
rnaseq_tpmcountData_v2 <- rnaseq_tpmcountData %>%
  rownames_to_column("versioned_ensembl_gene_id") %>%
  filter(versioned_ensembl_gene_id %in% gene_90_38_shortlist$versioned_ensembl_gene_id) %>%
  filter(!versioned_ensembl_gene_id %in% tpm_sum_infancy_subgroup$versioned_ensembl_gene_id) %>%
  filter(versioned_ensembl_gene_id %in% tpm_shortlist$versioned_ensembl_gene_id) %>%
  column_to_rownames("versioned_ensembl_gene_id")

mad_2020 <- mad_calculations(rnaseq_countData_v2, subject_specimen, c("2020_dataset"))
mad_2021 <- mad_calculations(rnaseq_countData_v2, subject_specimen, c("2021_dataset"))
mad_2022 <- mad_calculations(rnaseq_countData_v2, subject_specimen, c("2022_dataset"))
mad_2023 <- mad_calculations(rnaseq_countData_v2, subject_specimen, c("2023_dataset"))

mad_shotlisted_genes <- intersect(intersect(mad_2020$gene_id, mad_2021$gene_id), mad_2022$gene_id)
## mad_2023$gene_id genes were not included

rnaseq_tpmcountData_v3 <- rnaseq_tpmcountData_v2 %>%
  rownames_to_column("versioned_ensembl_gene_id") %>%
  filter(versioned_ensembl_gene_id %in% mad_shotlisted_genes) %>%
  column_to_rownames("versioned_ensembl_gene_id")
```
## RNASeq plot: tpm Count
```{r}
pvca_analysis_rnaseq(rnaseq_tpmcountData_v3, subject_specimen, batch.factors, plot_title = "RNASeq tpm: Raw data")

batch_labels = as.data.frame(colnames(rnaseq_tpmcountData_v3)) %>%
      rename(specimen_id = starts_with("colnames")) %>%
      mutate(specimen_id = as.integer(specimen_id)) %>%
      left_join(rnaseq_metaData) %>%
      dplyr::select(dataset)
    
rnaseq_tpmbatchCorrected = sva::ComBat_seq(as.matrix(rnaseq_tpmcountData_v3), batch = batch_labels$dataset)
    
pvca_analysis_rnaseq(rnaseq_tpmbatchCorrected, subject_specimen, batch.factors, plot_title = "RNASeq tpm: Batch correction")

rnaseq_normalised_data_tpm = list(
      
      raw_data = as.matrix(rnaseq_tpmcountData_v3),
      batchCorrected_data = rnaseq_tpmbatchCorrected
    )

rnaseq_normalised_data = list(
  
  raw_count = rnaseq_normalised_data_rawcount,
  tpm = rnaseq_normalised_data_tpm
)

```


## Save normalized data

```{r}
master_processed_data <- list(
  
  subject_specimen = subject_specimen,
  plasma_ab_titer = abtiter_processed_data,
  plasma_cytokine_concentrations_by_olink = olink_processed_data,
  plasma_cytokine_concentrations_by_legendplex = legendplex_processed_data,
  pbmc_cell_frequency = cell_frequency_processed_data,
  pbmc_gene_expression = rnaseq_normalised_data,
  t_cell_polarization= tcell_polarization_processed_data,
  t_cell_activation = tcell_activation_processed_data
  
)

#sapply(master_normalized_data$subject_specimen, dim)
sapply(master_processed_data$abtiter_wide, dim)
sapply(master_processed_data$plasma_cytokine_concentrations_by_olink, dim)
sapply(master_processed_data$plasma_cytokine_concentrations_by_legendplex, dim)
sapply(master_processed_data$pbmc_cell_frequency, dim)
sapply(master_processed_data$pbmc_gene_expression, dim)
sapply(master_processed_data$t_cell_polarization, dim)
sapply(master_processed_data$t_cell_activation, dim)


```

## Save data as RDS and individual TSV File

```{r}
saveRDS(master_processed_data, file = paste0(dir_RDS_objects, "master_processed_data_v20240825.RDS"))
#master_processed_data = readRDS(file = paste0(dir_RDS_objects, "master_allData_batchCorrected_v20240825.RDS"))

# Recursive function to save data frames (and lists): master_allData_batchCorrected_TSV
dir_rds_objects = dir_RDS_objects
save_dataframes_to_tsv(master_processed_data)


```

## session_info()

```{r}
sessioninfo::session_info()
```