create_final_esets.Rmd

---
title: "ImmuneSignatures2: Create Final ExpressionSets from Base ExpressionSet"
author: "Evan Henrich and Helen Miller"
output: html_document
params: 
  outputDir: "/share/files/HIPC/IS2/@files/data/html_outputs"
  dataCacheDir: "/share/files/HIPC/IS2/@files/data"
  timestamp: ""
---

# Overview
The purpose of this vignette is to generate expressionSets for analysis by HIPC collaborators in the ImmuneSignatures2 project.  The base expressionSet includes all age cohorts and is not normalized across study or batch-corrected for issues like different platforms (e.g. Illumina vs Affymetrix).  Therefore, expressionSets for analysis are generated for both young and old cohorts that then are cross-study normalized and batch corrected.  The cross-study normalization is done by taking a target distribution from one platform (Affymetrix) and applying this to all samples. The batch-correction is done by using a linear-model to determining the effects of platform (e.g. Illumina Human HT-V4), platform vendor (e.g. Illumina), cell type (either Whole Blood or PBMC), and study on baseline data, then removing these effects for all timepoints.  Finally each expressionSet comes in a version with or without immune response data.  For subjects with multiple immune response data points (due to multiple assays), the preference is given to HAI then NAb then ELISA data.

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, include = TRUE)
```

```{r load-dependencies}
library(ImmuneSignatures2) # vaccine map loaded as `vaccines`
library(Biobase)
library(dplyr)
library(data.table)
library(titer) # devtools::install_github("stefanavey/titer")
library(limma)
```

```{r output-variables}
outputDir <- params$outputDir
dataCacheDir <- params$dataCacheDir
inputFilenamePrefix <- params$timestamp
outputFilenamePrefix <- params$timestamp
if (!dir.exists(outputDir)) stop("could not find", outputDir)
if (!dir.exists(dataCacheDir)) stop("could not find", dataCacheDir)
```

```{r global-variables}
postVaxDayRanges <- list(
  hai = c(20,46),
  neut_ab_titer = c(28,90),
  elisa = c(21,30)
)

discretizationValues <- list(
    "RBA" = c(0.3, 0.4, 0.5),
    "MFC" = c(0.3, 0.4)
)

young <- c(18,50)
  old <- c(60,91)
  
ageGroups <- list(
  young = young,
  old = old,
  extendedOld = c(max(young), max(old)),
  all = c(min(young), max(old))
)

targetDistributionVendor <- "Affymetrix"
targetDistributionExcludedStudies <- "SDY1293"
```

```{r load-base-data}
immdata_all <- readRDS(file.path(dataCacheDir,  paste0(inputFilenamePrefix, "immdata_all.rds")))
noNormEset <- readRDS(file.path(dataCacheDir, paste0(inputFilenamePrefix, "noNormEset.rds")))
```

```{r create-cross-study-normalized-version}
eset.all.norm <- crossStudyNormalize(noNormEset, 
                                     targetDistributionVendor, 
                                     targetDistributionExcludedStudies)
```

```{r create-variants-of-final-eset}
for(ageGroupName in names(ageGroups)){
  ages <- ageGroups[[ageGroupName]]
  
  # No Cross-Study Normalization, No Immune Response Calls
  eset.noNorm.noResponse <- filterEsetByAge(noNormEset, ages)
  eset.noNorm.noResponse <- removeAllNArows(eset.noNorm.noResponse)
  res <- testFinalEset(eset.noNorm.noResponse, 
                       expectResponse = FALSE, 
                       expectNormalization = FALSE, 
                       ages)
  checkFinalTestResults(res)
  
  fullSuffix <- paste0(ageGroupName, "_noNorm_eset.rds")
  filename <- file.path(dataCacheDir, paste0(outputFilenamePrefix, fullSuffix))
  saveRDS(eset.noNorm.noResponse, filename)
  write_data_metadata(file.path(dataCacheDir, "dataset_metadata.csv"),
                      dataset_name = fullSuffix,
                      data_path = filename, 
                      data = eset.noNorm.noResponse, 
                      include_counts = TRUE)
  
  # With Cross-Study Normalization, No Immune Response Calls
  
  if(grepl("old", ageGroupName, ignore.case = TRUE)){
      eset.young.norm <- filterEsetByAge(eset.all.norm, ages = ageGroups[["young"]])
      eset.young.norm <- removeAllNArows(eset.young.norm)
      
      eset.old.norm <- filterEsetByAge(eset.all.norm, ages = ageGroups[[ageGroupName]])
      eset.old.norm <- removeAllNArows(eset.old.norm)
      
      # Remove studies that do not have young cohorts and cannot be modeled
      eset.old.norm <- eset.old.norm[, !eset.old.norm$study_accession %in% c("SDY1368","SDY67") ]
      
      eset.corr.noResponse <- batchCorrect.importedModel(
        modelEset = eset.young.norm,
        targetEset = eset.old.norm,
        batch.vars=c('cell_type',
                     'featureSetVendor',
                     'featureSetName2',
                     'geBatchName'),
        covariates = c('y_chrom_present')) 
  }else{
    eset.norm.noResponse <- filterEsetByAge(eset.all.norm, ages)
    eset.norm.noResponse <- removeAllNArows(eset.norm.noResponse)
    eset.norm.noResponse$age_group <- 
      factor(ifelse(eset.norm.noResponse$age_imputed < 50, "young", "extendedOld"), 
             levels = c("young", "extendedOld"))
    eset.corr.noResponse <- batchCorrect(
      eset.norm.noResponse, 
      batch.vars = c('cell_type',
                     'featureSetVendor',
                     'featureSetName2',
                     'geBatchName'),
      covariates = c('y_chrom_present', 
                     'age_group'))
    
  }
  
  eset.corr.noResponse <- removeAllNArows(eset.corr.noResponse)
  eset.corr.noResponse <- eset.corr.noResponse[complete.cases(exprs(eset.corr.noResponse)), ]
  
  res <- testFinalEset(eset.corr.noResponse, 
                       expectResponse = FALSE, 
                       expectNormalization = TRUE, 
                       ages)
  checkFinalTestResults(res)
  
  if (grepl("old", ageGroupName, ignore.case = TRUE)) {
    fullSuffix <- paste0(ageGroupName, "_norm_batchCorrectedFromYoung_eset.rds")
  } else {
    fullSuffix <- paste0(ageGroupName, "_norm_eset.rds")
  }
  filename <- file.path(dataCacheDir, paste0(outputFilenamePrefix, fullSuffix))
  saveRDS(eset.corr.noResponse, filename)  
  write_data_metadata(file.path(dataCacheDir, "dataset_metadata.csv"),
                      dataset_name = fullSuffix,
                      data_path = filename, 
                      data = eset.corr.noResponse, 
                      include_counts = TRUE)
  
  # Generate Immune Response Calls
  filteredImmdata <- filterImmdataByAge(immdata_all, ages)
  immdataWithResponses <- lapply(names(filteredImmdata), function(assay){
    dataWithResponses <- generateResponseCall(
      assay = assay,
      data = filteredImmdata[[assay]],
      postVaxDayRange = postVaxDayRanges[[assay]],
      discretizationValues = discretizationValues
    )
  })
  selectedImmdata <- selectResponsesToUse(immdataWithResponses)
  if (ageGroupName == "all") {
    filename <- file.path(dataCacheDir, "all_immdata_with_response.rds")
    saveRDS(selectedImmdata, filename)
    write_data_metadata(file.path(dataCacheDir, "dataset_metadata.csv"),
                        dataset_name = "all_immdata_with_response.rds",
                        data_path = filename, 
                        data = selectedImmdata, 
                        include_counts = TRUE)
  }
    
  # No Cross-Study Normalization, with Immune Response Calls
  eset.noNorm.withResponse <- addResponseData(eset.noNorm.noResponse,
                                              selectedImmdata)
  res <- testFinalEset(eset.noNorm.withResponse, 
                       expectResponse = TRUE, 
                       expectNormalization = FALSE, 
                       ages)
  checkFinalTestResults(res)
  
  fullSuffix <- paste0(ageGroupName, "_noNorm_withResponse_eset.rds")
  filename <- file.path(dataCacheDir, paste0(outputFilenamePrefix, fullSuffix))
  saveRDS(eset.noNorm.withResponse, filename)  
  write_data_metadata(file.path(dataCacheDir, "dataset_metadata.csv"),
                      dataset_name = fullSuffix,
                      data_path = filename, 
                      data = eset.noNorm.withResponse, 
                      include_counts = TRUE)
  
  
  # With Cross-Study Normalization, with Immune Response Calls
  eset.corr.withResponse <- addResponseData(eset.corr.noResponse, 
                                            selectedImmdata)
  res <- testFinalEset(eset.corr.withResponse, 
                       expectResponse = TRUE, 
                       expectNormalization = TRUE, 
                       ages)
  checkFinalTestResults(res)
  
  if (grepl("old", ageGroupName, ignore.case = TRUE)) {
    fullSuffix <- paste0(ageGroupName, "_norm_batchCorrectedFromYoung_withResponse_eset.rds")
  } else {
    fullSuffix <- paste0(ageGroupName, "_norm_withResponse_eset.rds")
  }
  filename <- file.path(dataCacheDir, paste0(outputFilenamePrefix, fullSuffix))
  saveRDS(eset.corr.withResponse, filename)
  write_data_metadata(file.path(dataCacheDir, "dataset_metadata.csv"),
                      dataset_name = fullSuffix,
                      data_path = filename, 
                      data = eset.corr.withResponse, 
                      include_counts = TRUE)
}
```

```{r create-diagnostic-plots, eval = FALSE}
allEset <- readRDS(file.path(dataCacheDir, paste0(outputFilenamePrefix, "all_noNorm_noResponse_eset.rds")))
ychromQCPlot <- qualityControl.failedYchromQC(allEset)
pdf(file = file.path(outputDir, paste0(outputFilenamePrefix, "ychromQC.pdf")),
    width = 8.5,
    height = 11)
ychromQCPlot
dev.off()

young <- readRDS(file.path(dataCacheDir, paste0(outputFilenamePrefix, "young_norm_noResponse_eset.rds")))

pdf(file = file.path(outputDir, paste0(outputFilenamePrefix, "young_mds_byStudy.pdf")),
    width = 8.5,
    height = 11)
qualityControl.samplePlot(young)
dev.off()

pdf(file = file.path(outputDir, paste0(outputFilenamePrefix, "young_mds_byYchrom.pdf")),
    width = 8.5,
    height = 11)
qualityControl.samplePlot(young, colorCol = "y_chrom_present")
dev.off()

pdf(file = file.path(outputDir, paste0(outputFilenamePrefix, "young_pca_byYchrom.pdf")),
    width = 8.5,
    height = 11)
qualityControl.samplePlot(young, method = "PCA", colorCol = "y_chrom_present")
dev.off()

old <- readRDS(file.path(dataCacheDir, paste0(outputFilenamePrefix, "old_norm_noResponse_eset.rds")))

pdf(file = file.path(outputDir, paste0(outputFilenamePrefix, "old_mds_byStudy.pdf")),
    width = 8.5,
    height = 11)
qualityControl.samplePlot(old)
dev.off()

pdf(file = file.path(outputDir, paste0(outputFilenamePrefix, "old_mds_byYchrom.pdf")),
    width = 8.5,
    height = 11)
qualityControl.samplePlot(old, colorCol = "y_chrom_present")
dev.off()

pdf(file = file.path(outputDir, paste0(outputFilenamePrefix, "old_pca_byYchrom.pdf")),
    width = 8.5,
    height = 11)
qualityControl.samplePlot(old, method = "PCA", colorCol = "y_chrom_present")
dev.off()

extendedOld <- readRDS(file.path(dataCacheDir, paste0(outputFilenamePrefix, "extendedOld_norm_noResponse_eset.rds")))

pdf(file = file.path(outputDir, paste0(outputFilenamePrefix, "extendedOld_mds_byStudy.pdf")),
    width = 8.5,
    height = 11)
qualityControl.samplePlot(extendedOld)
dev.off()
```