3-batch_correction_and_normalization.Rmd

---
title: "PTRC Ex10 Batch correction and normalization"
author: "Michael Nestor"
date: "`r format(Sys.Date(), '%m/%d/%Y')`"
output: html_document
---

```{r setup}
library(vp.misc)
library(sva)
library(dplyr)
library(tibble)
library(tidyr)

library(devtools)
source("../dataProcessing/proteomicsNormalizationMethods.R")
```

# Load metadata

```{r load_metadata}
library(readxl)

phenodata <- read_xlsx("./data/TMTmappingTable.xlsx") %>%
  as.data.frame() %>%
  mutate(Plex = as.factor(Plex)) %>%
  filter(`Sample ID \r\n(abbrev)` != "Ref")

names(phenodata) <- make.names(names(phenodata))

phenodata <- read_xlsx("data/CPTAC_FLT3cohort_UniquePatients_BestSamples_wYield.xlsx") %>%
  select(lab_id, specimen_type, specimen_location, Specimen_access_group_concatenated) %>%
  mutate(InitialAMLDiagnosis = grepl("Initial Acute Leukemia Diagnosis", Specimen_access_group_concatenated),
         PostChemotherapy    = grepl("Post-Chemotherapy", Specimen_access_group_concatenated)) %>%
  left_join(phenodata, ., by = c("Barcode.ID" = "lab_id")) %>%
  dplyr::rename(SampleID.full = Sample.ID..full.,
                SampleID.abbrev = Sample.ID....abbrev.)

phenodata <- phenodata %>%
  select(-New.Top.ID..for.Labeling,
         -Total.peptide....µg.,
         -Volume...uL.,
         -PlexMass)
names(phenodata) <- gsub("_", "\\.", names(phenodata)) 

rownames(phenodata) <- phenodata$`Sample ID \r\n(abbrev)`

write.table(phenodata, file="data/Ex10_metadata.txt",
            quote=F,sep="\t",row.names=F)
```

# Normalize and batch correct global data

```{r normalize_global_data}
normalize_global_data <- function(path_to_crosstab, path_to_phenodata, n.sigfig=3) {

  # Make MSnSet
  crosstab <- read.table(path_to_crosstab, check.names = F)
  m <- MSnSet(as.matrix(crosstab))
  phenodata <- read.table(path_to_phenodata, colClasses="character")
  phenodata$`Loading.Mass` <- as.numeric(phenodata$`Loading.Mass`)
  pData(m) <- phenodata[sampleNames(m),]
  
  # Medpolish and save
  m <- normalizeByMedpolish(m)
  new_path_to_crosstab <- sub("_original", "_medpolish", path_to_crosstab)
  write.table(signif(exprs(m), n.sigfig),
              file = new_path_to_crosstab,
              quote=F, sep="\t")
  
  # Missing value filter
  m <- m %>%
      filterByProportionMissingValues(least_proportion_threshold = 0.5) %>%
      filterByMissingPerBatch("Plex", least_count_threshold = 1L)
      
  # Batch correction
  removed_covariates <- c("Plex", "Loading.Mass")
  m <- correct_batch_effect_empiricalBayesLM(m, removed_covariates)
  
  # Medpolish and save
  m <- normalizeByMedpolish(m)
  new_path_to_crosstab <- sub("_original", "_corrected", path_to_crosstab)
  write.table(signif(exprs(m), n.sigfig),
              file = new_path_to_crosstab,
              quote=F, sep="\t")
}
```

# Normalize phospho data

```{r normalize_phospho_data}
normalize_phospho_data <- function(path_to_crosstab, path_to_phenodata,
                                   path_to_global_crosstab, n.sigfig=3) {
  # Make MSnSet
  crosstab <- read.table(path_to_crosstab, check.names = F)
  m <- MSnSet(as.matrix(crosstab))
  phenodata <- read.table(path_to_phenodata, colClasses="character")
  phenodata$`Loading.Mass` <- as.numeric(phenodata$`Loading.Mass`)
  pData(m) <- phenodata[sampleNames(m),]
  
  # Fetch global sample medians
  global_crosstab <- read.table(path_to_global_crosstab, check.names = F)
  global_coeffs <- apply(global_crosstab,
                         MARGIN = 2, FUN = median, na.rm = T)
  
  # Normalize by global sample medians
  exprs(m) <- sweep(exprs(m), 2, global_coeffs)
  m <- normalizeByMedpolish(m)
  new_path_to_crosstab <- sub("_original", "_medpolish", path_to_crosstab)
  write.table(signif(exprs(m), n.sigfig),
              file = new_path_to_crosstab,
              quote=F, sep="\t")
  
  # Missing value filter
  m <- m %>%
      filterByProportionMissingValues(least_proportion_threshold = 0.5) %>%
      filterByMissingPerBatch("Plex", least_count_threshold = 1L)
      
  # Batch correction
  removed_covariates <- c("Plex", "Loading.Mass")
  m <- correct_batch_effect_empiricalBayesLM(m, removed_covariates)
  
  # Medpolish and save
  m <- normalizeByMedpolish(m)
  new_path_to_crosstab <- sub("_original", "_corrected", path_to_crosstab)
  write.table(signif(exprs(m), n.sigfig),
              file = new_path_to_crosstab,
              quote=F, sep="\t")
  
  
    # Make MSnSet
  crosstab <- read.table(path_to_crosstab, check.names = F)
  m <- MSnSet(as.matrix(crosstab))
  pData(m) <- phenodata[sampleNames(m),]
  
  # Normalize by phospho coefficients
  m <- normalizeByMedpolish(m)
  new_path_to_crosstab <- sub("_original", "_medpolish_phospho_coeffs", path_to_crosstab)
  write.table(signif(exprs(m), n.sigfig),
              file = new_path_to_crosstab,
              quote=F, sep="\t")
  
  # Missing value filter
  m <- m %>%
      filterByProportionMissingValues(least_proportion_threshold = 0.5) %>%
      filterByMissingPerBatch("Plex", least_count_threshold = 1L)
      
  # Batch correction
  removed_covariates <- c("Plex", "Loading.Mass")
  m <- correct_batch_effect_empiricalBayesLM(m, removed_covariates)
  
  # Medpolish and save
  m <- normalizeByMedpolish(m)
  new_path_to_crosstab <- sub("_original", "_corrected_phospho_coeffs", path_to_crosstab)
  write.table(signif(exprs(m), n.sigfig),
              file = new_path_to_crosstab,
              quote=F, sep="\t")
}

```

# Main function calls

```{r main_loop}
t0 <- Sys.time(); print(t0)

lapply(list.files("data/Ex10_global_data/", "_original.txt",
                  full.names=T),
       normalize_global_data,
       path_to_phenodata = "data/Ex10_phenodata.txt")

lapply(list.files("data/Ex10_phospho_data/", "_original.txt",
                  full.names=T),
       normalize_phospho_data,
       path_to_phenodata = "data/Ex10_phenodata.txt",
       path_to_global_crosstab="data/Ex10_global_data/ptrc_ex10_crosstab_global_gene_original.txt")

t0 <- Sys.time(); print(t0)

t1 <- Sys.time(); print(t1); print(t1-t0)
```