knitr::opts_chunk$set(echo =TRUE, message = FALSE, warning = FALSE) suppressPackageStartupMessages({ library(Biobase) # library(gridExtra) # library(pheatmap) # library(BiostatsALL) # library(affy) library(ggplot2) # library(preprocessCore) library(limma) library(plyr) # library(sigaR) # library(MergeMaid) }) outputDir <- params$outputDir if (!dir.exists(outputDir)) dir.create(outputDir) dataDir <- params$dataCacheDir prefix <- params$timestamp getPCAs<-function (eset, maxpc = 3) { if (!is.matrix(eset)) { ex <- exprs(eset) out <- pData(eset) } else { out <- NULL ex <- eset } pca.res <- prcomp(t(ex)) x <- pca.res$x[, 1:maxpc] colnames(x) <- paste("PC", 1:maxpc, sep = ".") varex = round(100 * summary(pca.res)$importance[2, ]) db <- cbind.data.frame(x) if (!is.null(out)) db <- cbind(db, out) return(list(db = as.data.frame(db), varex = varex)) } plotPCA<-function(pca.db,title=paste0('PCA'), color.var='study2',var.shape='pathogen',var.size='months2' ){ p1<-ggplot(mutate(as.data.frame(pca.db$db), months=round((study_time_collected/28),1), months2=ifelse(months<0, -1, months)+2), aes_string(x='PC.1',y='PC.2', color=color.var, size=var.size, shape=var.shape)) + geom_jitter() + theme_bw() + scale_color_manual(values=ann.colors[[color.var]]) + labs(x=paste('PC-1 (',pca.db$varex[1],'%)',sep=''), y=paste('PC-2 (',pca.db$varex[2],'%)',sep=''), title=title) return(p1) } my.plot.pvca <- function(pvcaObj, cex.percentage = 1, fname = NULL, ht = 4, wd = 5, title = fname, race.vars=NULL, outputDir = NULL) { require(stringr) require(ggplot2) title <- gsub("/", "", fixed = TRUE, fname) title <- gsub(".", "", fixed = TRUE, fname) title <- gsub("PVCA", "", fixed = TRUE, fname) labels <- gsub(":", " x ", pvcaObj$label) labels <- gsub("_imputed", " ", labels) db <- data.frame(perc = t(pvcaObj$dat), labels = factor(labels, levels = labels)) if (!is.null(race.vars)) { db<-rbind(db, data.frame(perc=sum(subset(db, labels%in%race.vars)$perc), labels='race')) %>% subset(!(labels%in%race.vars)) %>% arrange(1*(labels=='resid'), perc) } p <- ggplot(db, aes(x = reorder(labels, perc), y = perc)) + geom_bar(stat = "identity", fill = "blue") + theme_bw() + scale_y_continuous(limits = c(0, 1.1)) + geom_text(aes(x = labels, y = perc + 0.05, label = paste0(as.character(round(100 * perc, 1)), "%"))) + labs(x = "Effects", y = "Weighted average proportion variance", title = paste("PVCA ", title)) + scale_x_discrete(labels = function(x) str_wrap(x, width = 8)) if (!is.null(fname)) { ggsave(file = file.path(outputDir, paste0(fname, ".pdf")), height = ht, width = wd, plot = p) } return(p) } impute.na<-function(x){x[is.na(x)]<-mean(x,na.rm=T); return(x)} Nfeatures.pca<-5000 race.vars <- c('Hispanic','White','Black','Asian')
eset_noNorm <- readRDS(file.path(dataDir, paste0(prefix, "all_noNorm_noResponse_eset.rds"))) # eset_crossStudyNorm <- readRDS(file.path(dataDir, paste0(prefix, "all_norm_noBatchCorrect_eset.rds"))) # eset_batchCorrected <- readRDS(file.path(dataDir, paste0(prefix, "all_norm_noResponse_eset.rds"))) # eset_noCelltype <- readRDS(file.path(dataDir, paste0(prefix, "all_norm_nocelltype_noResponse_eset.rds"))) # eset_old <- readRDS(file.path(dataDir, paste0(prefix, "old_norm_noResponse_eset.rds"))) # eset_all <- readRDS(file.path(dataDir, paste0(prefix, "all_norm_noResponse_eset.rds"))) # eset_young <- readRDS(file.path(dataDir, paste0(prefix, "young_norm_noResponse_eset.rds"))) # Nfeatures.pca<-5000 # eset_noNorm$pathogen_adjuvant <- interaction(eset_noNorm$pathogen, eset_noNorm$adjuvant) # eset_noNorm$time <- eset_noNorm$study_time_collected # eset_noNorm$age_group <- ifelse(eset_noNorm$age_imputed < 50,'young','old') # eset_noNorm$sex <- eset_noNorm$y_chrom_present # eset_crossStudyNorm$pathogen_adjuvant <- interaction(eset_crossStudyNorm$pathogen, eset_crossStudyNorm$adjuvant) # eset_crossStudyNorm$time <- eset_crossStudyNorm$study_time_collected # eset_crossStudyNorm$sex <- eset_crossStudyNorm$y_chrom_present # eset_crossStudyNorm$age_group <- ifelse(eset_crossStudyNorm$age_imputed < 50,'young','old') # eset_batchCorrected$pathogen_adjuvant <- interaction(eset_batchCorrected$pathogen, eset_batchCorrected$adjuvant) # eset_batchCorrected$time <- eset_batchCorrected$study_time_collected # eset_batchCorrected$sex <- eset_batchCorrected$y_chrom_present # eset_batchCorrected$age_group <- ifelse(eset_batchCorrected$age_imputed < 50,'young','old') # # eset_noCelltype$pathogen_adjuvant <- interaction(eset_noCelltype$pathogen, eset_noCelltype$adjuvant) # eset_noCelltype$time <- eset_noCelltype$study_time_collected # eset_noCelltype$sex <- eset_noCelltype$y_chrom_present # eset_noCelltype$age_group <- ifelse(eset_noCelltype$age_imputed < 50,'young','old') # # # race.vars <- c('Hispanic','White','Black','Asian')
load(here::here("vignettes", "original_code", "adj_path_vt_colors.rda")) # add for batch sc<-c("magenta","purple","orangered3","firebrick2","gray","green","lightblue4","yellow2", "cyan","firebrick4","royalblue3","darkorange2", "brown",'olivedrab','aquamarine','red','gold','dodgerblue','darkgray','lightsalmon1','darkgreen','royalblue3','orchid4', 'tan3','sandybrown','moccasin','pink','magenta','black','green','blue','violetred','snow4') ann.colors<-list(y_chrom_present=c('red','blue'), study_accession=sc, study_accession2=c(sc,'black'), featureSetName2=sc, vaccine=sc[1:length(unique(eset_noNorm$vaccine))]) ann.colors<-c(adj_path_vt_colors,ann.colors) path<-c(Meningitis="#80B1D3","Ebola"='gold', 'Hepatitis B'='purple', 'HIV'='red','Malaria'='cyan') ann.colors$pathogen<-c(path, ann.colors$pathogen[setdiff(names(ann.colors$pathogen), names(path))]) Nfeatures.pca<-5000
eset_noNorm_baseline <- eset_noNorm[, eset_noNorm$study_time_collected == 0] make_the_plots <- function(esetType) { eset <- get(paste0("eset_", esetType)) pcaFileName <- file.path(outputDir, paste0("pca_", esetType, "_baseline.rds")) if (!file.exists(pcaFileName)) { eset_baseline <- eset[, eset$study_time_collected == 0] features.pca.0 <- featureNames(eset_baseline)[ which( (rowMeans(is.na(exprs(eset_baseline))) == 0) & (apply(exprs(eset_baseline),1,sd) > 0) ) ] pca_baseline <- getPCAs(eset_baseline[features.pca.0,]) pca_baseline$db$size=factor(1) saveRDS(pca_baseline, pcaFileName) } else { pca_baseline <- readRDS(pcaFileName) } # PCA for baseline samples, colored by study and shaped by sequencing technology pcaPlot_baseline_v1 <- plotPCA(pca_baseline, title = paste0("PCA baseline ", esetType), color.var = 'study_accession', var.shape = 'featureSetName2', var.size = 'size') + scale_shape_manual(values=c(12,20,21:23,10,11,13:19)) ggsave(file = file.path(outputDir, paste0("pca_", esetType, "_baseline_v1.pdf")), plot = pcaPlot_baseline_v1, height = 6, width = 9) # PCA for baseline samples, colored by technology and shape by sample type pcaPlot_baseline_v2 <- plotPCA(pca_baseline, title = paste0("PCA baseline ", esetType), color.var = 'featureSetName2', var.shape = 'cell_type', var.size = 'size') + scale_shape_manual(values=c(12,20,21:23,10,11,13:19)) ggsave(file = file.path(outputDir, paste0("pca_", esetType, "_baseline_v2.pdf")), plot = pcaPlot_baseline_v2, height = 6, width = 9) pcaPlot_baseline_v3 <- plotPCA(pca_baseline, title = paste0("PCA Baseline by batch ", esetType), color.var = 'featureSetName2', var.shape = 'age_group', var.size = 'size') ggsave(file.path(outputDir, paste0("pca_", esetType, "_baseline_v3.pdf"))) print(pcaPlot_baseline_v1) print(pcaPlot_baseline_v2) print(pcaPlot_baseline_v3) }
library(data.table) eset_young <- readRDS("~/Projects/ImmuneSignatures2/data_cache/2021_02_11/2021_02_11_young_norm_noResponse_eset.rds") eset_young$age_group <- ifelse(eset_young$age_imputed < 50,'young','old') make_the_plots("young") eset_old <- readRDS("~/Projects/ImmuneSignatures2/data_cache/2021_02_11/2021_02_11_old_norm_noResponse_eset.rds") eset_old$age_group <- "old" make_the_plots("old") pd_old <- pData(eset_old) setDT(pd_old) pd_young[, .N, .(featureSetName, featureSetName2)] pd_young <- pData(eset_young) setDT(pd_young) pd_young[featureSetName2 == "HGU133_plus_PM", .N, .(study_accession)] ### Before batchGate ### eset_young_pre <- readRDS("~/Projects/ImmuneSignatures2/data_cache/old/2020_10_27_young_norm_noResponse_eset.rds") eset_old_pre <- readRDS("~/Projects/ImmuneSignatures2/data_cache/old/2020_10_27_old_norm_noResponse_eset.rds") eset_young_pre$age_group <- "young" eset_old_pre$age_group <- "old" make_the_plots("young_pre") make_the_plots("old_pre")
make_the_plots("noNorm")
# make_the_plots("crossStudyNorm")
make_the_plots("noCelltype")
make_the_plots("batchCorrected")
subjects age 18-50
eset_young_batchCorrected <- readRDS(file.path(dataDir, paste0(prefix, "young_norm_noResponse_eset.rds"))) eset_young_batchCorrected$pathogen_adjuvant <- interaction(eset_young_batchCorrected$pathogen, eset_young_batchCorrected$adjuvant) eset_young_batchCorrected$time <- eset_young_batchCorrected$study_time_collected eset_young_batchCorrected$sex <- eset_young_batchCorrected$y_chrom_present eset_young_batchCorrected$age_group <- ifelse(eset_young_batchCorrected$age_imputed < 50,'young','old') race.vars <- c('Hispanic','White','Black','Asian') make_the_plots("young_batchCorrected")
Subjects age 50-91
eset_old2_batchCorrected <- readRDS(file.path(dataDir, paste0(prefix, "old_norm_noResponse_eset.rds"))) eset_old2_batchCorrected$pathogen_adjuvant <- interaction(eset_old2_batchCorrected$pathogen, eset_old2_batchCorrected$adjuvant) eset_old2_batchCorrected$time <- eset_old2_batchCorrected$study_time_collected eset_old2_batchCorrected$sex <- eset_old2_batchCorrected$y_chrom_present eset_old2_batchCorrected$age_group <- ifelse(eset_old2_batchCorrected$age_imputed < 50,'young','old') make_the_plots("old2_batchCorrected")
eset_all_old <- eset_batchCorrected[,eset_batchCorrected$age_impute >= 60 ] make_the_plots("all_old")
eset_all_young <- eset_batchCorrected[,eset_batchCorrected$age_imputed < 60 ] make_the_plots("all_young")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.