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.
knitr::opts_chunk$set(echo = TRUE, include = TRUE)
library(ImmuneSignatures2) # vaccine map loaded as `vaccines` library(Biobase) library(dplyr) library(data.table) library(titer) # devtools::install_github("stefanavey/titer") library(limma)
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)
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"
immdata_all <- readRDS(file.path(dataCacheDir, paste0(inputFilenamePrefix, "immdata_all.rds"))) noNormEset <- readRDS(file.path(dataCacheDir, paste0(inputFilenamePrefix, "noNormEset.rds")))
eset.all.norm <- crossStudyNormalize(noNormEset, targetDistributionVendor, targetDistributionExcludedStudies)
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) }
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()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.