In RGLab/ImmuneSignatures2: Pre-Processing for HIPC ImmuneSignatures2 Analysis

suppressPackageStartupMessages(library(package = "knitr"))
suppressPackageStartupMessages(library(package = "Biobase"))
suppressPackageStartupMessages(library(package = "tidyverse"))

library (Biobase)
library (limma)
library (affy)
library (stringr)
library (dplyr)
library (tibble)
library (ggplot2)
library (pvca)
library (ExpressionNormalizationWorkflow)
library (sva)
library(ggplot2)
library(ggthemes)
library(reshape2)
library(plyr)
library (RColorBrewer)
library (beeswarm)
library (ggbeeswarm)
library(gridExtra)
library(grid)
library(Matching)
library(pheatmap)
library(janitor)
suppressMessages(library(tidyverse))

Preprocessing Data

library (knitr)
knitr::opts_chunk$set(echo = TRUE)
opts_knit$set (root.dir="~/Dropbox (BCH)/HIPC/IOF/IS2/RDA/CURRENT/CURRENT_2020_05_28")

rm(list=ls())

setwd("~/Dropbox (BCH)/HIPC/IOF/IS2/RDA/CURRENT/CURRENT_2020_05_28")
getwd()
workDir <- dirname(getwd())
IS2_eset_withResponse_norm_young <- readRDS ("2020_05_28_IS2_eset_withResponse_norm_young.rds") #this one
IS2_eset_withResponse_noNorm_young <- readRDS ("2020_05_28_IS2_eset_withResponse_noNorm_young.rds")
IS2_eset_noResponse_norm_young <- readRDS  ("2020_05_28_IS2_eset_noResponse_norm_young.rds") #this one
IS2_eset_noResponse_noNorm_young <- readRDS ("2020_05_28_IS2_eset_noResponse_noNorm_young.rds")
IS2_eset_withResponse_norm_older <- readRDS ("2020_05_28_IS2_eset_withResponse_norm_older.rds") #this one
IS2_eset_noResponse_norm_older <- readRDS ("2020_05_28_IS2_eset_noResponse_norm_older.rds") #this one
IS2_eset_withResponse_noNorm_older <- readRDS ("2020_05_28_IS2_eset_withResponse_noNorm_older.rds")
IS2_eset_noResponse_noNorm_older <- readRDS ("2020_05_28_IS2_eset_noResponse_noNorm_older.rds")
IS2_immdata <- readRDS ("2020_06_04_IS2_immdata.rds") #this one
IS2_immdata_selected <- readRDS ("2020_06_04_IS2_immdata_selected.rds")
immdata <- IS2_immdata

f<-intersect(featureNames(IS2_eset_noResponse_noNorm_older),featureNames(IS2_eset_noResponse_noNorm_young))
s<-intersect(varLabels(IS2_eset_noResponse_noNorm_older),varLabels(IS2_eset_noResponse_noNorm_young))
IS2_eset<-new("ExpressionSet", 
          exprs=cbind(exprs(IS2_eset_noResponse_noNorm_older[f,]), exprs(IS2_eset_noResponse_noNorm_young[f,])),
          phenoData=new('AnnotatedDataFrame',rbind(pData(IS2_eset_noResponse_noNorm_older)[,s],pData(IS2_eset_noResponse_noNorm_young)[,s])))

IS2_eset_noResponse_noNorm_young
IS2_eset_withResponse_noNorm_young
IS2_eset

Figures for Thomas' paper

IS2_all_GE_metaData<- pData (IS2_eset)

#Recode
IS2_all_GE_metaData$vaccine=recode(IS2_all_GE_metaData$vaccine, "Quadrivalent Meningococcal Vaccine plain-polysaccharide conjugate"="MenACWY-PS")
IS2_all_GE_metaData$vaccine=recode(IS2_all_GE_metaData$vaccine, "Quadrivalent Meningococcal Vaccine protein-polysaccharide conjugate"="MenACWY-CRM")
IS2_all_GE_metaData$vaccine =recode(IS2_all_GE_metaData$vaccine, "TIV(2011)"="TIV (2011)")
IS2_all_GE_metaData$vaccine =recode(IS2_all_GE_metaData$vaccine, "TIV(2013)"="TIV (2013)")
IS2_all_GE_metaData$pathogen_vaccinetype=paste0(IS2_all_GE_metaData$pathogen," (",IS2_all_GE_metaData$vaccine_type,")") 
unique (IS2_all_GE_metaData$pathogen_vaccinetype)

#Convert timepoint to numeric
IS2_all_GE_metaData$time_post_last_vax=as.numeric(IS2_all_GE_metaData$time_post_last_vax)
IS2_all_GE_metaData$time_post_last_vax

#Remove saline groups
IS2_all_GE_metaData<-IS2_all_GE_metaData%>% 
  filter(vaccine !="Saline") 
IS2_all_GE_metaData$time_post_last_vax=as.factor(IS2_all_GE_metaData$time_post_last_vax)
table (IS2_all_GE_metaData$time_post_last_vax)
IS2_all_GE_metaData$time_post_last_vax=revalue(IS2_all_GE_metaData$time_post_last_vax,
                                            c("0.166666666666667"="0.167", "21"=">20", "25"=">20", "28"=">20", "37"=">20", "41"=">20", "60"=">20", "70"=">20", "84"=">20"))

timepoint <- factor(IS2_all_GE_metaData$time_post_last_vax, levels=sort(as.numeric(unique(IS2_all_GE_metaData$time_post_last_vax))))
#Created variable for the counts to be used for plotting
count=rep(1,length(timepoint))
#Combined variables into one dataset

IS2_met=as.data.frame(cbind(IS2_all_GE_metaData, timepoint, count))

Manuscript Figure 1B (Thomas)

colourCount= length (unique (IS2_met$pathogen_vaccinetype))
#getPalette = colorRampPalette(brewer.pal(13, "Set3"))
getPalette= colorRampPalette (c(brewer.pal(name="Set3", n = 10), brewer.pal(name="Dark2", n = 4)))

#ggsave("Fig1BThomasIOF.png", width = 16, height =12, dpi=900)
ggplot(IS2_met, aes(x=time_post_last_vax, y=count, fill=pathogen_vaccinetype))+
  geom_bar(stat="identity")+
  #scale_x_continuous(breaks=seq(0,14))+
  #scale_y_continuous(breaks=c(0,250,500,750,1000))+
  scale_fill_manual(values=getPalette (colourCount)) +
  xlab("Day Post-Vaccination")+
  #ylab("# of Datasets")+
  ylab("Number of Samples")+
  labs(fill='Pathogen (Vaccine Type)') + theme_bw() +
  theme(text = element_text(size=20),
        axis.text.x = element_text(hjust=1, size=22, angle=45),
        axis.text.y = element_text(hjust=1, size=22)) +
  theme (legend.title = element_text(size=14, face="bold"), 
         legend.text = element_text(size = 12)) + 
  theme(legend.key=element_blank(), legend.key.size=unit(10,"point")) +
  theme (legend.position = "bottom") +
#  theme(plot.margin=unit(c(1,1,1,1),"cm")) +
  theme(legend.margin=margin(c(1,2,1,2), unit='cm')) +
  guides(fill=guide_legend(ncol=2,byrow=TRUE)) + ggtitle ("Summary of Pathogen (Vaccine Type) for Transcriptomics Dataset")+theme(plot.title = element_text(size=16))

Figure 1C

####Focus on Day 0 subjects####
IS2_met$age_reported[is.na(IS2_met$age_reported)] <- 0
IS2_met$age_reported <-as.integer (as.character(IS2_met$age_reported))
table (IS2_met$age_reported)
#IS2_all_GE_metaData$age_reported <- cut(IS2_all_GE_metaData$age_reported, breaks=c(-Inf,0, 18, 45, 65, Inf), 
#labels=c("Not Specified", "Age below 18","Age 18-45","Age 45-65", "Age over 65"))
IS2_met$age_imputed[is.na(IS2_met$age_imputed)] <- 0
IS2_met$age_imputed <-as.integer (as.character(IS2_met$age_imputed))
#table (IS2_met$age_imputed)
#colnames (IS2_eset)
D0<- grepl(pattern = '_0_Days', colnames(IS2_eset))
D0eset <- IS2_eset[, D0]
#dim (D0eset)
D0pdata <-IS2_met[match(colnames(D0eset), table=IS2_met$uid),]
#dim (D0pdata) #[1] 1222 36
#colnames (D0pdata)
#table (D0pdata$study_accession)

#dim (D0eset) 
#Features  Samples 
#25860    1222
#options(max.print=1000000)
#colnames(D0eset)
D0eset <- as.data.frame(D0eset)

#Filter those with 0 age for now
D0pdata <- D0pdata %>%
  filter (age_imputed !=0) 

D0pdata<-D0pdata %>% 
  filter(vaccine != "Saline" ) 
vaccine <- D0pdata %>%
  group_by(vaccine) %>%
  dplyr::summarise(counts = n())

age_reported <-D0pdata$age_reported
age_imputed <- D0pdata$age_imputed
race <- D0pdata$race
gender_imputed <-D0pdata$gender_imputed
matrix <- D0pdata$matrix
#matrix
#head (D0pdata)[1:20]

ggplot(D0pdata, aes(x=pathogen_vaccinetype, y=age_imputed)) + geom_boxplot(outlier.shape=NA) + 
  theme_light() + #scale_shape_manual(values=1:nlevels(vaccine)) +
  theme(axis.text.x = element_text(size=8, angle=60, hjust = 1, face = "bold")) +
  theme(axis.text.y = element_text(size = 11, hjust = 1, face = "bold")) +
  scale_y_continuous(breaks = pretty(age_imputed, n = 10)) +
  ylab("Age") + xlab ("Study")+
  #ggtitle("HIPC Signatures 2 Summary: Comparative Meta-analysis Study Metadata") +
  theme(plot.title = element_text(face="bold", size=11)) +
  theme(axis.title = element_text(face="bold", size=9)) +
  scale_fill_manual(values=getPalette (colourCount)) +
  #scale_color_viridis_d(option="plasma", direction=-1) +
  geom_jitter(aes(color=vaccine, shape=gender_imputed), position=position_jitter(width=.15, height=.1))+
  theme (legend.title = element_text(size=10, face="bold"), 
         legend.text = element_text(size = 10)) + 
  theme(legend.key=element_blank(), legend.key.size=unit(10,"point")) +
  theme (legend.position = "bottom") #+
#  theme(plot.margin=unit(c(1,1,1,1),"cm")) +
  #theme(legend.margin=margin(c(1,2,1,2), unit='cm'))