projects/03_aim3/r_old/v_18.03.21_chmi_aim_03.R
In mvazquezs/chmitools: Package to manage & obtain results of the CHMI Study
##############################################################
### CODE TO ANALYZE CHMI AB DATA FOLLOWING THE PAPER: AIM3 ###
##############################################################
# Start: 2017-06-16
rm(list=ls())
library(Hmisc)
library(plyr)
library(compareGroups)
manual.theme <- theme(panel.background = element_rect(fill = "white",
colour = NA), panel.border = element_rect(fill = NA,
colour = "grey20"), panel.grid.major = element_line(colour = "grey92"),
panel.grid.minor = element_line(colour = "grey92",
size = 0.25), strip.background = element_rect(fill = "grey85",
colour = "grey20"), legend.key = element_rect(fill = "white",
colour = NA), complete = TRUE, plot.title = element_text(hjust = 0.5))
fcn.pvalue.display<-function(p.vec, show.p = FALSE){
p.vec.plot<-as.character(rep(NA,length(p.vec)))
if(any((p.vec >= .01 & p.vec <= .04) | p.vec >= .06)==TRUE) {
p.vec.plot<-as.character(round(p.vec,2))
p.vec.plot[nchar(as.character(p.vec.plot))==3 & p.vec.plot > 0.1]<-paste(as.character(p.vec.plot)[nchar(as.character(p.vec.plot))==3 & p.vec.plot > 0.1],"0",sep="")
}
p.vec.plot<-ifelse ((p.vec >= 0.001 & p.vec < 0.01) | (p.vec > 0.04 & p.vec < 0.06), as.character(round(p.vec,3)),p.vec.plot)
if(show.p){
p.vec.plot[p.vec < 0.001]<-c("P<0.001")
p.vec.plot[p.vec >= 0.001]<-c(paste0("P=", p.vec.plot))
} else {
p.vec.plot[p.vec < 0.001]<-c("<0.001")
}
return(p.vec.plot)
}
localpath <- "/Users/hsanz/Dropbox (M067)/m067_team/working_areas/MALIMMUNO Team Support/3. CHMI study Antibodies"
############################################################################################################################################
## load("/Users/hsanz/Dropbox (M067)/m067_team/working_areas/MALIMMUNO Team Support/3. CHMI study Antibodies/Data/1. All isotypes data 20171107/All_data_20171218.Rdata")
load('~/datasets_ori/v_17.12.18_chmi/v_17.12.18_chmi_dat.RData')
### SOME RECODES
#################
dat <- chmi_ab_l
names(dat) <- tolower(names(dat))
dim(dat)
head(dat)
xxx <- unique(dat[,c("original_id", "dataset", "immune_status", "status")])
with(xxx, table(status, dataset))
with(xxx, table(immune_status, dataset))
xxx <- unique(dat[,c("original_id", "dataset", "hb_status","status")])
with(xxx[xxx$dataset=="L1",], table(hb_status, status))
## Some recodes
dat$timepoint <- as.character(dat$timepoint)
dat$timepoint <- with(dat, ifelse(dataset=="T2" & timepoint=="DM", "D11", timepoint))
dat$hb <- as.character(dat$hb)
dat$hb <- with(dat, ifelse(dataset=="L1", "", hb))
dat$cvac_dose <- with(dat, ifelse(dataset=="T2", cvac_dose, NA))
dat$hb_status_benjamin <- NULL
dat$pre_patent_period <- as.character(dat$pre_patent_period)
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-005" & dataset=="T1", "9.5",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-018" & dataset=="T1", "7",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-039" & dataset=="T1", "11",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-040" & dataset=="T1", "8",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-041" & dataset=="T1", "9",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-043" & dataset=="T1", "7.5",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-044" & dataset=="T1", "10",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-045" & dataset=="T1", "11.5",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-048" & dataset=="T1", "12.5",pre_patent_period))
dat$gender <- ""
dat$gender <- with(dat, ifelse(original_id=="T1-005" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-018" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-039" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-040" & dataset=="T1", "F",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-041" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-043" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-044" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-045" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-048" & dataset=="T1", "F",gender))
datx <- dat
datx$log10_mfi <- NULL
datx$batch <- NULL
datx$hb <- NULL
datx$cvac_dose <- NULL
datx$pre_patent_period <- with(datx, ifelse(is.na(pre_patent_period), "NA", pre_patent_period))
datx$pcr.pos <- with(datx, ifelse(is.na(pcr.pos), 999, pcr.pos))
datx <- datx[, c("original_id", "timepoint","ig","study_number",
"pre_patent_period", "dataset", "antigen","mfi","pcr.pos","hb_status",
"immune_status", "status", "gender", "malaria")]
datx <- aggregate(mfi ~ original_id + timepoint + ig +
+ study_number + malaria
+ pre_patent_period + dataset + antigen + pcr.pos + immune_status
+ status + gender + hb_status, FUN = mean, data=datx)
xxx <- as.data.frame(with(datx[datx$antigen=="CSP",], table(original_id,ig,study_number, timepoint)))
xxx[xxx$Freq>1,]
dat <- datx
dat$log10_mfi <- log10(dat$mfi)
#### PREPARE DATASET FOR ANALYSIS
##################################
tmp.dat <- subset(dat, dataset=="L1")
tmp.dat <- subset(tmp.dat, antigen%nin%c("BSA", "agal"))
with(tmp.dat, table(dataset, useNA = "always"))
with(tmp.dat, table(hb_status, useNA = "always"))
with(tmp.dat, table(hb_status, status, useNA = "always"))
tmp.dat$group <- with(tmp.dat, ifelse(status=="Naive", "Naive", hb_status))
with(tmp.dat, table(hb_status, status, useNA = "always"))
with(tmp.dat, table(group, useNA = "always"))
with(tmp.dat, table(timepoint, useNA = "always"))
tmp.dat$group <- factor(tmp.dat$group, levels=c("Naive", "AA", "AS"))
lantigen <- sort(unique(tmp.dat$antigen))
tmp.dat$timepoint2c <- tmp.dat$timepoint
tmp.dat$timepoint2c <- factor(tmp.dat$timepoint2c, levels=c("C-1", "D7", "D13","D19","D28"))
########################################
########################################
########################################
### AIM 3 Analysis
########################################
########################################
########################################
## Differences in MFIs before challenge
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$antigen==x,]
ff.ans <- try(oneway.test(log10_mfi ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ff.pval(isotype, lantigen)
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x, show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
ymax <- max(dataset[dataset$ig==isotype, "log10_mfi"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(group, log10_mfi))
q <- q + geom_boxplot(outlier.size = NA) + geom_jitter(width = 0.2, height = 0, color="black", alpha=0.25)
q <- q + facet_wrap( ~ antigen) + theme_bw() + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Group") + ylab("log10(MFI)")
q <- q + ggtitle(isotype) + geom_text(data=pvals, aes(x=x.text, y=y.text, label=adj.pvalue.text), size=3)
q <- q + ylim(y.text - (abs(y.text)*0.85), ymax)
ans <- list(q = q, pvals = pvals)
return(ans)
}
igg <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG", lantigen, 2, 2)
ig1 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG1", lantigen, 2, 0.5)
ig2 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG2", lantigen, 2, 1)
ig3 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG3", lantigen, 2, 0.5)
ig4 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG4", lantigen, 2, -2)
igm <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgM", lantigen, 2, 1.5)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "before_challenge_all.pdf"), width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
#############################################################
# After challenge within timepoint: only three timepoints
#############################################################
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen, timepointx){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$timepoint2c==timepointx & dataset$antigen==x,]
ff.ans <- try(oneway.test(log10_mfi ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx,antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ldply(c("D7", "D13", "D28") ,function(x) ff.pval(isotype, lantigen, x))
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x,show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
pvals$log10_mfi <- y.text
pvals$group <- NA
ymax <- max(dataset[dataset$ig==isotype, "log10_mfi"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group))
q <- q + geom_boxplot(outlier.size = NA)
q <- q + geom_point(position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0), alpha=0.25)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
q <- q + geom_text(data=pvals, aes(x=timepoint2c, y=log10_mfi, label=adj.pvalue.text), size=3)
q <- q + ylim(y.text - (abs(y.text)*0.85), ymax)
ans <- list(q = q, pvals = pvals)
return(ans)
}
# D19 removed because of no info of Naive
dataset <- tmp.dat[tmp.dat$timepoint%in%c("D7", "D13", "D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("D7", "D13", "D28"))
igg <- ff.figure(dataset, "IgG", lantigen, NA, 1.5)
ig1 <- ff.figure(dataset, "IgG1", lantigen, 2, 0.5)
ig2 <- ff.figure(dataset, "IgG2", lantigen, NA, 1)
ig3 <- ff.figure(dataset, "IgG3", lantigen, NA, 0.5)
ig4 <- ff.figure(dataset, "IgG4", lantigen, 2, -2)
igm <- ff.figure(dataset, "IgM", lantigen, 2, 1.5)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "after_challenge_within_timepoint_all.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
#############################################################
# After challenge within timepoint: only four timepoints
#############################################################
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen, timepointx){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$timepoint2c==timepointx & dataset$antigen==x,]
ff.ans <- try(oneway.test(log10_mfi ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx,antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ldply(c("C-1","D7", "D13", "D28") ,function(x) ff.pval(isotype, lantigen, x))
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x,show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
pvals$log10_mfi <- y.text
pvals$group <- NA
ymax <- max(dataset[dataset$ig==isotype, "log10_mfi"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group))
q <- q + geom_boxplot(outlier.size = NA)
q <- q + geom_point(position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0), alpha=0.25)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
q <- q + geom_text(data=pvals, aes(x=timepoint2c, y=log10_mfi, label=adj.pvalue.text), size=3)
q <- q + ylim(y.text - (abs(y.text)*0.85), ymax)
ans <- list(q = q, pvals = pvals)
return(ans)
}
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13","D28"))
igg <- ff.figure(dataset, "IgG", lantigen, NA, 1.5)
ig1 <- ff.figure(dataset, "IgG1", lantigen, 2, 0.2)
ig2 <- ff.figure(dataset, "IgG2", lantigen, 2, 1)
ig3 <- ff.figure(dataset, "IgG3", lantigen, 2, 0.5)
ig4 <- ff.figure(dataset, "IgG4", lantigen, 2, -2)
igm <- ff.figure(dataset, "IgM", lantigen, 2, 1.5)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "after_challenge_within_timepoint_all_v02.pdf"),
width = 14, height = 9)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
######################################################
# After challenge within Status: three timepoints
######################################################
ff.figure <- function(dataset, isotype, lantigen){
summ.group <- ddply(dataset[dataset$ig==isotype, ], .(antigen,timepoint2c, group), summarise, med = median(log10_mfi))
summ.group$original_id <- NA
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group, group=original_id))
q <- q + geom_line(alpha=0.2)
# q <- q + geom_point(alpha=0.2)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + geom_point(data=summ.group, aes(timepoint2c, med, color=group), size=2)
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
ans <- list(q = q)
return(ans)
}
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13", "D28"))
igg <- ff.figure(dataset, "IgG", lantigen)
ig1 <- ff.figure(dataset, "IgG1", lantigen)
ig2 <- ff.figure(dataset, "IgG2", lantigen)
ig3 <- ff.figure(dataset, "IgG3", lantigen)
ig4 <- ff.figure(dataset, "IgG4", lantigen)
igm <- ff.figure(dataset, "IgM", lantigen)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "after_challenge_within_status_all.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
######################################################
# After challenge within Status: three timepoints
######################################################
ff.figure <- function(dataset, isotype, lantigen){
summ.group <- ddply(dataset[dataset$ig==isotype, ], .(antigen,timepoint2c, group), summarise, med = median(log10_mfi))
summ.group$original_id <- NA
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group, group=original_id))
q <- q + geom_line(alpha=0.2)
# q <- q + geom_point(alpha=0.2)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + geom_point(data=summ.group, aes(timepoint2c, med, color=group), size=2)
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
ans <- list(q = q)
return(ans)
}
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D19" ,"D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13", "D19" ,"D28"))
igg <- ff.figure(dataset, "IgG", lantigen)
ig1 <- ff.figure(dataset, "IgG1", lantigen)
ig2 <- ff.figure(dataset, "IgG2", lantigen)
ig3 <- ff.figure(dataset, "IgG3", lantigen)
ig4 <- ff.figure(dataset, "IgG4", lantigen)
igm <- ff.figure(dataset, "IgM", lantigen)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "after_challenge_within_status_all_v02.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
########################
## FOLD CHANGE
########################
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D19","D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13", "D19" ,"D28"))
# compute fold change
dataset.wide <- dataset
dataset.wide$timepoint <- NULL
dataset.wide$log10_mfi <- NULL
dataset.wide <- dataset.wide[, c("original_id", "ig", "study_number","antigen", "mfi", "timepoint2c", "group")]
dataset.wide <- reshape(data = dataset.wide, direction="wide",
idvar = c("original_id", "ig", "antigen", "study_number", "group"),
v.names = c("mfi") ,
timevar = "timepoint2c")
dataset.wide$mfi.D7 <- dataset.wide[, "mfi.D7"] / dataset.wide[, "mfi.C-1"]
dataset.wide[, "mfi.D13"] <- dataset.wide[, "mfi.D13"] / dataset.wide[, "mfi.C-1"]
dataset.wide[, "mfi.D19"] <- dataset.wide[, "mfi.D19"] / dataset.wide[, "mfi.C-1"]
dataset.wide[, "mfi.D28"] <- dataset.wide[, "mfi.D28"] / dataset.wide[, "mfi.C-1"]
dataset.ratio <- reshape::melt(data=dataset.wide,
id.vars=c("original_id", "ig", "antigen", "study_number", "group"))
dataset.ratio <- subset(dataset.ratio, variable!="mfi.C-1")
dataset.ratio$timepoint2c <- gsub(pattern = "mfi.", replacement = "", x = dataset.ratio$variable, fixed=TRUE)
dataset.ratio <- subset(dataset.ratio, timepoint2c!="D19")
dataset.ratio$timepoint2c <- factor(dataset.ratio$timepoint2c, levels = c("D7", "D13", "D28"))
# function.to.plot
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen, timepointx){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$timepoint2c==timepointx & dataset$antigen==x,]
ff.ans <- try(oneway.test(value ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx,antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ldply(c("D7", "D13", "D28") ,function(x) ff.pval(isotype, lantigen, x))
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x,show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
pvals$value <- y.text
pvals$group <- NA
ymax <- max(dataset[dataset$ig==isotype, "value"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, value, color=group))
q <- q + geom_boxplot(outlier.size = NA) + scale_y_log10(breaks=c(0.01,1,100,1000), labels=c(0.01,1,100,1000))
q <- q + geom_point(position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0), alpha=0.25)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("MFI ratio (Timepoint/C-1)")
q <- q + geom_text(data=pvals, aes(x=timepoint2c, y=value, label=adj.pvalue.text), size=3)
# q <- q + ylim(log10(y.text) - (abs(log10(y.text))*0.85), log10(ymax))
ans <- list(q = q, pvals = pvals)
return(ans)
}
igg <- ff.figure(dataset.ratio, "IgG", lantigen, 2, 0.1)
ig1 <- ff.figure(dataset.ratio, "IgG1", lantigen, 2, 0.001)
ig2 <- ff.figure(dataset.ratio, "IgG2", lantigen, 2, 0.01)
ig3 <- ff.figure(dataset.ratio, "IgG3", lantigen, 2, 0.01)
ig4 <- ff.figure(dataset.ratio, "IgG4", lantigen, 2, 0.001)
igm <- ff.figure(dataset.ratio, "IgM", lantigen, 2, 0.01)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "after_challenge_within_timepoint_all_fold_change.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
#####################################
#####################################
#####################################
### Repeat figures without NAIVE
#####################################
#####################################
#####################################
## Differences in MFIs before challenge
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$antigen==x,]
ff.ans <- try(oneway.test(log10_mfi ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ff.pval(isotype, lantigen)
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x, show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
ymax <- max(dataset[dataset$ig==isotype, "log10_mfi"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(group, log10_mfi))
q <- q + geom_boxplot(outlier.size = NA) + geom_jitter(width = 0.2, height = 0, color="black", alpha=0.25)
q <- q + facet_wrap( ~ antigen) + theme_bw() + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Group") + ylab("log10(MFI)")
q <- q + ggtitle(isotype) + geom_text(data=pvals, aes(x=x.text, y=y.text, label=adj.pvalue.text), size=3)
q <- q + ylim(y.text - (abs(y.text)*0.85), ymax)
ans <- list(q = q, pvals = pvals)
return(ans)
}
igg <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1" & tmp.dat$group%nin%"Naive",], "IgG", lantigen, 1.5, 2)
ig1 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1" & tmp.dat$group%nin%"Naive",], "IgG1", lantigen, 1.5, 0.5)
ig2 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1" & tmp.dat$group%nin%"Naive",], "IgG2", lantigen, 1.5, 1)
ig3 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1" & tmp.dat$group%nin%"Naive",], "IgG3", lantigen, 1.5, 0.5)
ig4 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1" & tmp.dat$group%nin%"Naive",], "IgG4", lantigen, 1.5, -2)
igm <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1" & tmp.dat$group%nin%"Naive",], "IgM", lantigen, 1.5, 1.5)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_before_challenge_all.pdf"), width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
#############################################################
# After challenge within timepoint: only three timepoints
#############################################################
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen, timepointx){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$timepoint2c==timepointx & dataset$antigen==x,]
ff.ans <- try(oneway.test(log10_mfi ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx,antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ldply(c("D7", "D13", "D28") ,function(x) ff.pval(isotype, lantigen, x))
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x,show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
pvals$log10_mfi <- y.text
pvals$group <- NA
ymax <- max(dataset[dataset$ig==isotype, "log10_mfi"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group))
q <- q + geom_boxplot(outlier.size = NA)
q <- q + geom_point(position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0), alpha=0.25)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
q <- q + geom_text(data=pvals, aes(x=timepoint2c, y=log10_mfi, label=adj.pvalue.text), size=3)
q <- q + ylim(y.text - (abs(y.text)*0.85), ymax)
ans <- list(q = q, pvals = pvals)
return(ans)
}
# D19 removed because of no info of Naive
dataset <- tmp.dat[tmp.dat$timepoint%in%c("D7", "D13", "D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("D7", "D13", "D28"))
igg <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG", lantigen, NA, 1.5)
ig1 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG1", lantigen, 2, 0.5)
ig2 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG2", lantigen, NA, 1)
ig3 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG3", lantigen, NA, 0.5)
ig4 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG4", lantigen, 2, -2)
igm <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgM", lantigen, 2, 1.5)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_after_challenge_within_timepoint_all.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
#############################################################
# After challenge within timepoint: only four timepoints
#############################################################
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen, timepointx){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$timepoint2c==timepointx & dataset$antigen==x,]
ff.ans <- try(oneway.test(log10_mfi ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx,antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ldply(c("C-1","D7", "D13", "D28") ,function(x) ff.pval(isotype, lantigen, x))
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x,show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
pvals$log10_mfi <- y.text
pvals$group <- NA
ymax <- max(dataset[dataset$ig==isotype, "log10_mfi"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group))
q <- q + geom_boxplot(outlier.size = NA)
q <- q + geom_point(position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0), alpha=0.25)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
q <- q + geom_text(data=pvals, aes(x=timepoint2c, y=log10_mfi, label=adj.pvalue.text), size=3)
q <- q + ylim(y.text - (abs(y.text)*0.85), ymax)
ans <- list(q = q, pvals = pvals)
return(ans)
}
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13","D28"))
igg <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG", lantigen, NA, 1.5)
ig1 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG1", lantigen, 2, 0.2)
ig2 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG2", lantigen, 2, 1)
ig3 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG3", lantigen, 2, 0.5)
ig4 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG4", lantigen, 2, -2)
igm <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgM", lantigen, 2, 1.5)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_after_challenge_within_timepoint_all_v02.pdf"),
width = 14, height = 9)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
######################################################
# After challenge within Status: three timepoints
######################################################
ff.figure <- function(dataset, isotype, lantigen){
summ.group <- ddply(dataset[dataset$ig==isotype, ], .(antigen,timepoint2c, group), summarise, med = median(log10_mfi))
summ.group$original_id <- NA
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group, group=original_id))
q <- q + geom_line(alpha=0.2)
# q <- q + geom_point(alpha=0.2)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + geom_point(data=summ.group, aes(timepoint2c, med, color=group), size=2)
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
ans <- list(q = q)
return(ans)
}
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13", "D28"))
igg <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG", lantigen)
ig1 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG1", lantigen)
ig2 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG2", lantigen)
ig3 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG3", lantigen)
ig4 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG4", lantigen)
igm <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgM", lantigen)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_after_challenge_within_status_all.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
######################################################
# After challenge within Status: three timepoints
######################################################
ff.figure <- function(dataset, isotype, lantigen){
summ.group <- ddply(dataset[dataset$ig==isotype, ], .(antigen,timepoint2c, group), summarise, med = median(log10_mfi))
summ.group$original_id <- NA
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group, group=original_id))
q <- q + geom_line(alpha=0.2)
# q <- q + geom_point(alpha=0.2)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + geom_point(data=summ.group, aes(timepoint2c, med, color=group), size=2)
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
ans <- list(q = q)
return(ans)
}
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D19" ,"D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13", "D19" ,"D28"))
igg <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG", lantigen)
ig1 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG1", lantigen)
ig2 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG2", lantigen)
ig3 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG3", lantigen)
ig4 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG4", lantigen)
igm <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgM", lantigen)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_after_challenge_within_status_all_v02.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
########################
## FOLD CHANGE
########################
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D19","D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13", "D19" ,"D28"))
# compute fold change
dataset.wide <- dataset
dataset.wide$timepoint <- NULL
dataset.wide$log10_mfi <- NULL
dataset.wide <- dataset.wide[, c("original_id", "ig", "study_number","antigen", "mfi", "timepoint2c", "group")]
dataset.wide <- reshape(data = dataset.wide, direction="wide",
idvar = c("original_id", "ig", "antigen", "study_number", "group"),
v.names = c("mfi") ,
timevar = "timepoint2c")
dataset.wide$mfi.D7 <- dataset.wide[, "mfi.D7"] / dataset.wide[, "mfi.C-1"]
dataset.wide[, "mfi.D13"] <- dataset.wide[, "mfi.D13"] / dataset.wide[, "mfi.C-1"]
dataset.wide[, "mfi.D19"] <- dataset.wide[, "mfi.D19"] / dataset.wide[, "mfi.C-1"]
dataset.wide[, "mfi.D28"] <- dataset.wide[, "mfi.D28"] / dataset.wide[, "mfi.C-1"]
dataset.ratio <- reshape::melt(data=dataset.wide,
id.vars=c("original_id", "ig", "antigen", "study_number", "group"))
dataset.ratio <- subset(dataset.ratio, variable!="mfi.C-1")
dataset.ratio$timepoint2c <- gsub(pattern = "mfi.", replacement = "", x = dataset.ratio$variable, fixed=TRUE)
dataset.ratio <- subset(dataset.ratio, timepoint2c!="D19")
dataset.ratio$timepoint2c <- factor(dataset.ratio$timepoint2c, levels = c("D7", "D13", "D28"))
# function.to.plot
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen, timepointx){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$timepoint2c==timepointx & dataset$antigen==x,]
ff.ans <- try(oneway.test(value ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx,antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ldply(c("D7", "D13", "D28") ,function(x) ff.pval(isotype, lantigen, x))
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x,show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
pvals$value <- y.text
pvals$group <- NA
ymax <- max(dataset[dataset$ig==isotype, "value"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, value, color=group))
q <- q + geom_boxplot(outlier.size = NA) + scale_y_log10(breaks=c(0.01,1,100,1000), labels=c(0.01,1,100,1000))
q <- q + geom_point(position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0), alpha=0.25)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("MFI ratio (Timepoint/C-1)")
q <- q + geom_text(data=pvals, aes(x=timepoint2c, y=value, label=adj.pvalue.text), size=3)
# q <- q + ylim(log10(y.text) - (abs(log10(y.text))*0.85), log10(ymax))
ans <- list(q = q, pvals = pvals)
return(ans)
}
igg <- ff.figure(dataset.ratio[dataset.ratio$group%nin%"Naive",], "IgG", lantigen, 2, 0.1)
ig1 <- ff.figure(dataset.ratio[dataset.ratio$group%nin%"Naive",], "IgG1", lantigen, 2, 0.001)
ig2 <- ff.figure(dataset.ratio[dataset.ratio$group%nin%"Naive",], "IgG2", lantigen, 2, 0.01)
ig3 <- ff.figure(dataset.ratio[dataset.ratio$group%nin%"Naive",], "IgG3", lantigen, 2, 0.01)
ig4 <- ff.figure(dataset.ratio[dataset.ratio$group%nin%"Naive",], "IgG4", lantigen, 2, 0.001)
igm <- ff.figure(dataset.ratio[dataset.ratio$group%nin%"Naive",], "IgM", lantigen, 2, 0.01)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_after_challenge_within_timepoint_all_fold_change.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
###########################################
###########################################
### Mixed models: NAIVE
###########################################
###########################################
library(nlme)
library(lmerTest)
dataset <- tmp.dat[tmp.dat$timepoint2c%in%c("D7", "D13", "D19", "D28","C-1"),]
dataset$newtimepoint <- with(dataset, ifelse(timepoint2c=="C-1",0,
ifelse(timepoint2c=="D7",7,
ifelse(timepoint2c=="D13",13,
ifelse(timepoint2c=="D19",19,
ifelse(timepoint2c=="D28",28,NA))))))
with(dataset, table(ig))
with(dataset, table(antigen))
with(dataset, table(dataset))
with(dataset, table(newtimepoint))
with(dataset, table(timepoint2c))
with(dataset, table(group))
### All models and isotype
library(multcomp)
ff1 <- function(model){
# Slope Naive
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0,0,0))))
a1.est <- round(a$test$coefficients, 3)
a1.pvalue <- as.numeric(a$test$pvalues)
# Slope AA
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0,1,0))))
a2.est <- round(a$test$coefficients, 3)
a2.pvalue <- as.numeric(a$test$pvalues)
a2.pvalue.int <- summary(model)$coefficients[5,5]
# Slope AS
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0,0,1))))
a3.est <- round(a$test$coefficients, 3)
a3.pvalue <- as.numeric(a$test$pvalues)
a3.pvalue.int <- summary(model)$coefficients[6,5]
# Summary
ans <- data.frame(a1.est, a1.pvalue, a2.est, a2.pvalue, a2.pvalue.int,
a3.est, a3.pvalue, a3.pvalue.int)
}
ans.dat <- list()
zz <- 0
for(i in c("IgG", "IgG1", "IgG2", "IgG3", "IgG4", "IgM") ){
z <- 0
ans1 <- list()
for(j in lantigen){
z <- z + 1
aux.data <- subset(dataset, ig==i & antigen==j) # INCLUDE BASELINE. THINK BASELINE
if(nrow(aux.data)!=0){
# All timepoints
mod <- lmer(log10_mfi ~ newtimepoint*group + (newtimepoint | original_id), data=aux.data)
ans <- try(ff1(mod), silent = TRUE)
if(inherits(ans, "try-error")){
ans <- data.frame(a1.est=NA, a1.pvalue=NA, a2.est=NA, a2.pvalue=NA, a2.pvalue.int=NA,
a3.est=NA, a3.pvalue=NA, a3.pvalue.int=NA)}
ans$isotype <- i
ans$antigen <- j
ans$model <- "three.timepoints"
ans1[[z]] <- ans
}
}
zz <- zz + 1
ans1 <- ldply(ans1)
ans1[, c(2,4,5,7,8)] <- sapply(c(2,4,5,7,8), function(x) fcn.pvalue.display(p.adjust(ans1[,x], method="BH")))
ans1 <- ans1[,c(9,10, 1:8)]
ans.dat[[zz]] <- ans1
}
ans.dat.summ <- ans.dat
#################################
### Export models to word
#################################
library(ReporteRs)
#########################################################
## Create Tables in Word describing
#########################################################
mydoc = docx(file.path(localpath, "Results", "Paper", "aim3", "mfi", "mixed_models.docx"))
xvector <- c("IgG", "IgG1", "IgG2", "IgG3", "IgG4","IgGM")
for(i in 1:6){
my_text = pot(paste0("Table", i, ". Mixed models (log10MFI response models) by antigen fitting the HbS (and Naive) and all timepoints, D7, D13, D19 and D84 (including C-1 information) for ", xvector[i] , " isotype. Coefficient and adjusted p value is shown (reference: Naive)."))
mydoc = addParagraph( mydoc, value = my_text, stylename = "Normal")
mydoc = addParagraph( mydoc, value = "", stylename = "Normal")
ftab = FlexTable(ans.dat.summ[[i]] ,add.rownames = FALSE, header.columns=FALSE,
body.par.props = parCenter(),
header.par.props = parCenter(),
body.text.props = textProperties( font.size = 6))
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","Slope Naive", "Slope AA", "Int. AA",
"Slope AS","Int. AS"),
colspan = c(1,1,2,2,1,2,1), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("Isotype", "Antigen","Coef", "p value", rep(c("Coef", "p value", "p value"),2)),
colspan = rep(1,10), par.properties = parCenter())
mydoc = addFlexTable( mydoc, ftab, par.properties = parCenter())
mydoc <- addPageBreak(mydoc)
}
writeDoc( mydoc, file = file.path(localpath, "Results", "Paper", "aim3", "mfi", "mixed_models.docx"))
###########################################
###########################################
### Mixed models: REMOVING NAIVE
###########################################
###########################################
library(nlme)
library(lmerTest)
dataset <- tmp.dat[tmp.dat$timepoint2c%in%c("D7", "D13", "D19", "D28","C-1"),]
dataset$newtimepoint <- with(dataset, ifelse(timepoint2c=="C-1",0,
ifelse(timepoint2c=="D7",7,
ifelse(timepoint2c=="D13",13,
ifelse(timepoint2c=="D19",19,
ifelse(timepoint2c=="D28",28,NA))))))
dataset <- subset(dataset, group!="Naive")
with(dataset, table(ig))
with(dataset, table(antigen))
with(dataset, table(dataset))
with(dataset, table(newtimepoint))
with(dataset, table(timepoint2c))
with(dataset, table(group))
### All models and isotype
library(multcomp)
ff1 <- function(model){
# Slope AA
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0))))
a1.est <- round(a$test$coefficients, 3)
a1.pvalue <- as.numeric(a$test$pvalues)
# Slope AS
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,1))))
a2.est <- round(a$test$coefficients, 3)
a2.pvalue <- as.numeric(a$test$pvalues)
a2.pvalue.int <- summary(model)$coefficients[4,5]
# Summary
ans <- data.frame(a1.est, a1.pvalue, a2.est, a2.pvalue, a2.pvalue.int)
}
ans.dat <- list()
zz <- 0
for(i in c("IgG", "IgG1", "IgG2", "IgG3", "IgG4", "IgM") ){
z <- 0
ans1 <- list()
for(j in lantigen){
z <- z + 1
aux.data <- subset(dataset, ig==i & antigen==j) # INCLUDE BASELINE. THINK BASELINE
if(nrow(aux.data)!=0){
# All timepoints
mod <- lmer(log10_mfi ~ newtimepoint*group + (newtimepoint | original_id), data=aux.data)
ans <- try(ff1(mod), silent = TRUE)
if(inherits(ans, "try-error")){
ans <- data.frame(a1.est=NA, a1.pvalue=NA, a2.est=NA, a2.pvalue=NA, a2.pvalue.int=NA,
a3.est=NA, a3.pvalue=NA, a3.pvalue.int=NA)}
ans$isotype <- i
ans$antigen <- j
ans$model <- "three.timepoints"
ans1[[z]] <- ans
}
}
zz <- zz + 1
ans1 <- ldply(ans1)
ans1[, c(2,4,5)] <- sapply(c(2,4,5), function(x) fcn.pvalue.display(p.adjust(ans1[,x], method="BH")))
ans1 <- ans1[,c(6,7, 1:5)]
ans.dat[[zz]] <- ans1
}
ans.dat.summ <- ans.dat
#################################
### Export models to word
#################################
library(ReporteRs)
#########################################################
## Create Tables in Word describing
#########################################################
mydoc = docx(file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_mixed_models.docx"))
xvector <- c("IgG", "IgG1", "IgG2", "IgG3", "IgG4","IgGM")
for(i in 1:6){
my_text = pot(paste0("Table", i, ". Mixed models (log10MFI response models) by antigen fitting the HbS (Naive group removed) and all timepoints, D7, D13, D19 and D84 (including C-1 information) for ", xvector[i] , " isotype. Coefficient and adjusted p value is shown (reference: AA)."))
mydoc = addParagraph( mydoc, value = my_text, stylename = "Normal")
mydoc = addParagraph( mydoc, value = "", stylename = "Normal")
ftab = FlexTable(ans.dat.summ[[i]] ,add.rownames = FALSE, header.columns=FALSE,
body.par.props = parCenter(),
header.par.props = parCenter(),
body.text.props = textProperties( font.size = 6))
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","Slope AA", "Slope AS", "Int. AS"),
colspan = c(1,1,2,2,1), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("Isotype", "Antigen","Coef", "p value", rep(c("Coef", "p value", "p value"),1)),
colspan = rep(1,7), par.properties = parCenter())
mydoc = addFlexTable( mydoc, ftab, par.properties = parCenter())
mydoc <- addPageBreak(mydoc)
}
writeDoc( mydoc, file = file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_mixed_models.docx"))
mvazquezs/chmitools documentation built on May 1, 2020, 2:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
##############################################################
### CODE TO ANALYZE CHMI AB DATA FOLLOWING THE PAPER: AIM3 ###
##############################################################
# Start: 2017-06-16
rm(list=ls())
library(Hmisc)
library(plyr)
library(compareGroups)
manual.theme <- theme(panel.background = element_rect(fill = "white",
colour = NA), panel.border = element_rect(fill = NA,
colour = "grey20"), panel.grid.major = element_line(colour = "grey92"),
panel.grid.minor = element_line(colour = "grey92",
size = 0.25), strip.background = element_rect(fill = "grey85",
colour = "grey20"), legend.key = element_rect(fill = "white",
colour = NA), complete = TRUE, plot.title = element_text(hjust = 0.5))
fcn.pvalue.display<-function(p.vec, show.p = FALSE){
p.vec.plot<-as.character(rep(NA,length(p.vec)))
if(any((p.vec >= .01 & p.vec <= .04) | p.vec >= .06)==TRUE) {
p.vec.plot<-as.character(round(p.vec,2))
p.vec.plot[nchar(as.character(p.vec.plot))==3 & p.vec.plot > 0.1]<-paste(as.character(p.vec.plot)[nchar(as.character(p.vec.plot))==3 & p.vec.plot > 0.1],"0",sep="")
}
p.vec.plot<-ifelse ((p.vec >= 0.001 & p.vec < 0.01) | (p.vec > 0.04 & p.vec < 0.06), as.character(round(p.vec,3)),p.vec.plot)
if(show.p){
p.vec.plot[p.vec < 0.001]<-c("P<0.001")
p.vec.plot[p.vec >= 0.001]<-c(paste0("P=", p.vec.plot))
} else {
p.vec.plot[p.vec < 0.001]<-c("<0.001")
}
return(p.vec.plot)
}
localpath <- "/Users/hsanz/Dropbox (M067)/m067_team/working_areas/MALIMMUNO Team Support/3. CHMI study Antibodies"
############################################################################################################################################
## load("/Users/hsanz/Dropbox (M067)/m067_team/working_areas/MALIMMUNO Team Support/3. CHMI study Antibodies/Data/1. All isotypes data 20171107/All_data_20171218.Rdata")
load('~/datasets_ori/v_17.12.18_chmi/v_17.12.18_chmi_dat.RData')
### SOME RECODES
#################
dat <- chmi_ab_l
names(dat) <- tolower(names(dat))
dim(dat)
head(dat)
xxx <- unique(dat[,c("original_id", "dataset", "immune_status", "status")])
with(xxx, table(status, dataset))
with(xxx, table(immune_status, dataset))
xxx <- unique(dat[,c("original_id", "dataset", "hb_status","status")])
with(xxx[xxx$dataset=="L1",], table(hb_status, status))
## Some recodes
dat$timepoint <- as.character(dat$timepoint)
dat$timepoint <- with(dat, ifelse(dataset=="T2" & timepoint=="DM", "D11", timepoint))
dat$hb <- as.character(dat$hb)
dat$hb <- with(dat, ifelse(dataset=="L1", "", hb))
dat$cvac_dose <- with(dat, ifelse(dataset=="T2", cvac_dose, NA))
dat$hb_status_benjamin <- NULL
dat$pre_patent_period <- as.character(dat$pre_patent_period)
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-005" & dataset=="T1", "9.5",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-018" & dataset=="T1", "7",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-039" & dataset=="T1", "11",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-040" & dataset=="T1", "8",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-041" & dataset=="T1", "9",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-043" & dataset=="T1", "7.5",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-044" & dataset=="T1", "10",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-045" & dataset=="T1", "11.5",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-048" & dataset=="T1", "12.5",pre_patent_period))
dat$gender <- ""
dat$gender <- with(dat, ifelse(original_id=="T1-005" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-018" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-039" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-040" & dataset=="T1", "F",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-041" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-043" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-044" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-045" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-048" & dataset=="T1", "F",gender))
datx <- dat
datx$log10_mfi <- NULL
datx$batch <- NULL
datx$hb <- NULL
datx$cvac_dose <- NULL
datx$pre_patent_period <- with(datx, ifelse(is.na(pre_patent_period), "NA", pre_patent_period))
datx$pcr.pos <- with(datx, ifelse(is.na(pcr.pos), 999, pcr.pos))
datx <- datx[, c("original_id", "timepoint","ig","study_number",
"pre_patent_period", "dataset", "antigen","mfi","pcr.pos","hb_status",
"immune_status", "status", "gender", "malaria")]
datx <- aggregate(mfi ~ original_id + timepoint + ig +
+ study_number + malaria
+ pre_patent_period + dataset + antigen + pcr.pos + immune_status
+ status + gender + hb_status, FUN = mean, data=datx)
xxx <- as.data.frame(with(datx[datx$antigen=="CSP",], table(original_id,ig,study_number, timepoint)))
xxx[xxx$Freq>1,]
dat <- datx
dat$log10_mfi <- log10(dat$mfi)
#### PREPARE DATASET FOR ANALYSIS
##################################
tmp.dat <- subset(dat, dataset=="L1")
tmp.dat <- subset(tmp.dat, antigen%nin%c("BSA", "agal"))
with(tmp.dat, table(dataset, useNA = "always"))
with(tmp.dat, table(hb_status, useNA = "always"))
with(tmp.dat, table(hb_status, status, useNA = "always"))
tmp.dat$group <- with(tmp.dat, ifelse(status=="Naive", "Naive", hb_status))
with(tmp.dat, table(hb_status, status, useNA = "always"))
with(tmp.dat, table(group, useNA = "always"))
with(tmp.dat, table(timepoint, useNA = "always"))
tmp.dat$group <- factor(tmp.dat$group, levels=c("Naive", "AA", "AS"))
lantigen <- sort(unique(tmp.dat$antigen))
tmp.dat$timepoint2c <- tmp.dat$timepoint
tmp.dat$timepoint2c <- factor(tmp.dat$timepoint2c, levels=c("C-1", "D7", "D13","D19","D28"))
########################################
########################################
########################################
### AIM 3 Analysis
########################################
########################################
########################################
## Differences in MFIs before challenge
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$antigen==x,]
ff.ans <- try(oneway.test(log10_mfi ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ff.pval(isotype, lantigen)
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x, show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
ymax <- max(dataset[dataset$ig==isotype, "log10_mfi"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(group, log10_mfi))
q <- q + geom_boxplot(outlier.size = NA) + geom_jitter(width = 0.2, height = 0, color="black", alpha=0.25)
q <- q + facet_wrap( ~ antigen) + theme_bw() + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Group") + ylab("log10(MFI)")
q <- q + ggtitle(isotype) + geom_text(data=pvals, aes(x=x.text, y=y.text, label=adj.pvalue.text), size=3)
q <- q + ylim(y.text - (abs(y.text)*0.85), ymax)
ans <- list(q = q, pvals = pvals)
return(ans)
}
igg <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG", lantigen, 2, 2)
ig1 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG1", lantigen, 2, 0.5)
ig2 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG2", lantigen, 2, 1)
ig3 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG3", lantigen, 2, 0.5)
ig4 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG4", lantigen, 2, -2)
igm <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgM", lantigen, 2, 1.5)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "before_challenge_all.pdf"), width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
#############################################################
# After challenge within timepoint: only three timepoints
#############################################################
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen, timepointx){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$timepoint2c==timepointx & dataset$antigen==x,]
ff.ans <- try(oneway.test(log10_mfi ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx,antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ldply(c("D7", "D13", "D28") ,function(x) ff.pval(isotype, lantigen, x))
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x,show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
pvals$log10_mfi <- y.text
pvals$group <- NA
ymax <- max(dataset[dataset$ig==isotype, "log10_mfi"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group))
q <- q + geom_boxplot(outlier.size = NA)
q <- q + geom_point(position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0), alpha=0.25)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
q <- q + geom_text(data=pvals, aes(x=timepoint2c, y=log10_mfi, label=adj.pvalue.text), size=3)
q <- q + ylim(y.text - (abs(y.text)*0.85), ymax)
ans <- list(q = q, pvals = pvals)
return(ans)
}
# D19 removed because of no info of Naive
dataset <- tmp.dat[tmp.dat$timepoint%in%c("D7", "D13", "D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("D7", "D13", "D28"))
igg <- ff.figure(dataset, "IgG", lantigen, NA, 1.5)
ig1 <- ff.figure(dataset, "IgG1", lantigen, 2, 0.5)
ig2 <- ff.figure(dataset, "IgG2", lantigen, NA, 1)
ig3 <- ff.figure(dataset, "IgG3", lantigen, NA, 0.5)
ig4 <- ff.figure(dataset, "IgG4", lantigen, 2, -2)
igm <- ff.figure(dataset, "IgM", lantigen, 2, 1.5)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "after_challenge_within_timepoint_all.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
#############################################################
# After challenge within timepoint: only four timepoints
#############################################################
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen, timepointx){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$timepoint2c==timepointx & dataset$antigen==x,]
ff.ans <- try(oneway.test(log10_mfi ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx,antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ldply(c("C-1","D7", "D13", "D28") ,function(x) ff.pval(isotype, lantigen, x))
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x,show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
pvals$log10_mfi <- y.text
pvals$group <- NA
ymax <- max(dataset[dataset$ig==isotype, "log10_mfi"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group))
q <- q + geom_boxplot(outlier.size = NA)
q <- q + geom_point(position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0), alpha=0.25)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
q <- q + geom_text(data=pvals, aes(x=timepoint2c, y=log10_mfi, label=adj.pvalue.text), size=3)
q <- q + ylim(y.text - (abs(y.text)*0.85), ymax)
ans <- list(q = q, pvals = pvals)
return(ans)
}
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13","D28"))
igg <- ff.figure(dataset, "IgG", lantigen, NA, 1.5)
ig1 <- ff.figure(dataset, "IgG1", lantigen, 2, 0.2)
ig2 <- ff.figure(dataset, "IgG2", lantigen, 2, 1)
ig3 <- ff.figure(dataset, "IgG3", lantigen, 2, 0.5)
ig4 <- ff.figure(dataset, "IgG4", lantigen, 2, -2)
igm <- ff.figure(dataset, "IgM", lantigen, 2, 1.5)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "after_challenge_within_timepoint_all_v02.pdf"),
width = 14, height = 9)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
######################################################
# After challenge within Status: three timepoints
######################################################
ff.figure <- function(dataset, isotype, lantigen){
summ.group <- ddply(dataset[dataset$ig==isotype, ], .(antigen,timepoint2c, group), summarise, med = median(log10_mfi))
summ.group$original_id <- NA
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group, group=original_id))
q <- q + geom_line(alpha=0.2)
# q <- q + geom_point(alpha=0.2)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + geom_point(data=summ.group, aes(timepoint2c, med, color=group), size=2)
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
ans <- list(q = q)
return(ans)
}
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13", "D28"))
igg <- ff.figure(dataset, "IgG", lantigen)
ig1 <- ff.figure(dataset, "IgG1", lantigen)
ig2 <- ff.figure(dataset, "IgG2", lantigen)
ig3 <- ff.figure(dataset, "IgG3", lantigen)
ig4 <- ff.figure(dataset, "IgG4", lantigen)
igm <- ff.figure(dataset, "IgM", lantigen)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "after_challenge_within_status_all.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
######################################################
# After challenge within Status: three timepoints
######################################################
ff.figure <- function(dataset, isotype, lantigen){
summ.group <- ddply(dataset[dataset$ig==isotype, ], .(antigen,timepoint2c, group), summarise, med = median(log10_mfi))
summ.group$original_id <- NA
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group, group=original_id))
q <- q + geom_line(alpha=0.2)
# q <- q + geom_point(alpha=0.2)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + geom_point(data=summ.group, aes(timepoint2c, med, color=group), size=2)
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
ans <- list(q = q)
return(ans)
}
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D19" ,"D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13", "D19" ,"D28"))
igg <- ff.figure(dataset, "IgG", lantigen)
ig1 <- ff.figure(dataset, "IgG1", lantigen)
ig2 <- ff.figure(dataset, "IgG2", lantigen)
ig3 <- ff.figure(dataset, "IgG3", lantigen)
ig4 <- ff.figure(dataset, "IgG4", lantigen)
igm <- ff.figure(dataset, "IgM", lantigen)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "after_challenge_within_status_all_v02.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
########################
## FOLD CHANGE
########################
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D19","D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13", "D19" ,"D28"))
# compute fold change
dataset.wide <- dataset
dataset.wide$timepoint <- NULL
dataset.wide$log10_mfi <- NULL
dataset.wide <- dataset.wide[, c("original_id", "ig", "study_number","antigen", "mfi", "timepoint2c", "group")]
dataset.wide <- reshape(data = dataset.wide, direction="wide",
idvar = c("original_id", "ig", "antigen", "study_number", "group"),
v.names = c("mfi") ,
timevar = "timepoint2c")
dataset.wide$mfi.D7 <- dataset.wide[, "mfi.D7"] / dataset.wide[, "mfi.C-1"]
dataset.wide[, "mfi.D13"] <- dataset.wide[, "mfi.D13"] / dataset.wide[, "mfi.C-1"]
dataset.wide[, "mfi.D19"] <- dataset.wide[, "mfi.D19"] / dataset.wide[, "mfi.C-1"]
dataset.wide[, "mfi.D28"] <- dataset.wide[, "mfi.D28"] / dataset.wide[, "mfi.C-1"]
dataset.ratio <- reshape::melt(data=dataset.wide,
id.vars=c("original_id", "ig", "antigen", "study_number", "group"))
dataset.ratio <- subset(dataset.ratio, variable!="mfi.C-1")
dataset.ratio$timepoint2c <- gsub(pattern = "mfi.", replacement = "", x = dataset.ratio$variable, fixed=TRUE)
dataset.ratio <- subset(dataset.ratio, timepoint2c!="D19")
dataset.ratio$timepoint2c <- factor(dataset.ratio$timepoint2c, levels = c("D7", "D13", "D28"))
# function.to.plot
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen, timepointx){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$timepoint2c==timepointx & dataset$antigen==x,]
ff.ans <- try(oneway.test(value ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx,antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ldply(c("D7", "D13", "D28") ,function(x) ff.pval(isotype, lantigen, x))
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x,show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
pvals$value <- y.text
pvals$group <- NA
ymax <- max(dataset[dataset$ig==isotype, "value"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, value, color=group))
q <- q + geom_boxplot(outlier.size = NA) + scale_y_log10(breaks=c(0.01,1,100,1000), labels=c(0.01,1,100,1000))
q <- q + geom_point(position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0), alpha=0.25)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("MFI ratio (Timepoint/C-1)")
q <- q + geom_text(data=pvals, aes(x=timepoint2c, y=value, label=adj.pvalue.text), size=3)
# q <- q + ylim(log10(y.text) - (abs(log10(y.text))*0.85), log10(ymax))
ans <- list(q = q, pvals = pvals)
return(ans)
}
igg <- ff.figure(dataset.ratio, "IgG", lantigen, 2, 0.1)
ig1 <- ff.figure(dataset.ratio, "IgG1", lantigen, 2, 0.001)
ig2 <- ff.figure(dataset.ratio, "IgG2", lantigen, 2, 0.01)
ig3 <- ff.figure(dataset.ratio, "IgG3", lantigen, 2, 0.01)
ig4 <- ff.figure(dataset.ratio, "IgG4", lantigen, 2, 0.001)
igm <- ff.figure(dataset.ratio, "IgM", lantigen, 2, 0.01)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "after_challenge_within_timepoint_all_fold_change.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
#####################################
#####################################
#####################################
### Repeat figures without NAIVE
#####################################
#####################################
#####################################
## Differences in MFIs before challenge
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$antigen==x,]
ff.ans <- try(oneway.test(log10_mfi ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ff.pval(isotype, lantigen)
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x, show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
ymax <- max(dataset[dataset$ig==isotype, "log10_mfi"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(group, log10_mfi))
q <- q + geom_boxplot(outlier.size = NA) + geom_jitter(width = 0.2, height = 0, color="black", alpha=0.25)
q <- q + facet_wrap( ~ antigen) + theme_bw() + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Group") + ylab("log10(MFI)")
q <- q + ggtitle(isotype) + geom_text(data=pvals, aes(x=x.text, y=y.text, label=adj.pvalue.text), size=3)
q <- q + ylim(y.text - (abs(y.text)*0.85), ymax)
ans <- list(q = q, pvals = pvals)
return(ans)
}
igg <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1" & tmp.dat$group%nin%"Naive",], "IgG", lantigen, 1.5, 2)
ig1 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1" & tmp.dat$group%nin%"Naive",], "IgG1", lantigen, 1.5, 0.5)
ig2 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1" & tmp.dat$group%nin%"Naive",], "IgG2", lantigen, 1.5, 1)
ig3 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1" & tmp.dat$group%nin%"Naive",], "IgG3", lantigen, 1.5, 0.5)
ig4 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1" & tmp.dat$group%nin%"Naive",], "IgG4", lantigen, 1.5, -2)
igm <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1" & tmp.dat$group%nin%"Naive",], "IgM", lantigen, 1.5, 1.5)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_before_challenge_all.pdf"), width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
#############################################################
# After challenge within timepoint: only three timepoints
#############################################################
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen, timepointx){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$timepoint2c==timepointx & dataset$antigen==x,]
ff.ans <- try(oneway.test(log10_mfi ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx,antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ldply(c("D7", "D13", "D28") ,function(x) ff.pval(isotype, lantigen, x))
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x,show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
pvals$log10_mfi <- y.text
pvals$group <- NA
ymax <- max(dataset[dataset$ig==isotype, "log10_mfi"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group))
q <- q + geom_boxplot(outlier.size = NA)
q <- q + geom_point(position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0), alpha=0.25)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
q <- q + geom_text(data=pvals, aes(x=timepoint2c, y=log10_mfi, label=adj.pvalue.text), size=3)
q <- q + ylim(y.text - (abs(y.text)*0.85), ymax)
ans <- list(q = q, pvals = pvals)
return(ans)
}
# D19 removed because of no info of Naive
dataset <- tmp.dat[tmp.dat$timepoint%in%c("D7", "D13", "D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("D7", "D13", "D28"))
igg <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG", lantigen, NA, 1.5)
ig1 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG1", lantigen, 2, 0.5)
ig2 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG2", lantigen, NA, 1)
ig3 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG3", lantigen, NA, 0.5)
ig4 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG4", lantigen, 2, -2)
igm <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgM", lantigen, 2, 1.5)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_after_challenge_within_timepoint_all.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
#############################################################
# After challenge within timepoint: only four timepoints
#############################################################
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen, timepointx){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$timepoint2c==timepointx & dataset$antigen==x,]
ff.ans <- try(oneway.test(log10_mfi ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx,antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ldply(c("C-1","D7", "D13", "D28") ,function(x) ff.pval(isotype, lantigen, x))
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x,show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
pvals$log10_mfi <- y.text
pvals$group <- NA
ymax <- max(dataset[dataset$ig==isotype, "log10_mfi"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group))
q <- q + geom_boxplot(outlier.size = NA)
q <- q + geom_point(position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0), alpha=0.25)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
q <- q + geom_text(data=pvals, aes(x=timepoint2c, y=log10_mfi, label=adj.pvalue.text), size=3)
q <- q + ylim(y.text - (abs(y.text)*0.85), ymax)
ans <- list(q = q, pvals = pvals)
return(ans)
}
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13","D28"))
igg <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG", lantigen, NA, 1.5)
ig1 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG1", lantigen, 2, 0.2)
ig2 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG2", lantigen, 2, 1)
ig3 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG3", lantigen, 2, 0.5)
ig4 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG4", lantigen, 2, -2)
igm <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgM", lantigen, 2, 1.5)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_after_challenge_within_timepoint_all_v02.pdf"),
width = 14, height = 9)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
######################################################
# After challenge within Status: three timepoints
######################################################
ff.figure <- function(dataset, isotype, lantigen){
summ.group <- ddply(dataset[dataset$ig==isotype, ], .(antigen,timepoint2c, group), summarise, med = median(log10_mfi))
summ.group$original_id <- NA
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group, group=original_id))
q <- q + geom_line(alpha=0.2)
# q <- q + geom_point(alpha=0.2)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + geom_point(data=summ.group, aes(timepoint2c, med, color=group), size=2)
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
ans <- list(q = q)
return(ans)
}
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13", "D28"))
igg <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG", lantigen)
ig1 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG1", lantigen)
ig2 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG2", lantigen)
ig3 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG3", lantigen)
ig4 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG4", lantigen)
igm <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgM", lantigen)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_after_challenge_within_status_all.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
######################################################
# After challenge within Status: three timepoints
######################################################
ff.figure <- function(dataset, isotype, lantigen){
summ.group <- ddply(dataset[dataset$ig==isotype, ], .(antigen,timepoint2c, group), summarise, med = median(log10_mfi))
summ.group$original_id <- NA
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=group, group=original_id))
q <- q + geom_line(alpha=0.2)
# q <- q + geom_point(alpha=0.2)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + geom_point(data=summ.group, aes(timepoint2c, med, color=group), size=2)
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("log10(MFI)")
ans <- list(q = q)
return(ans)
}
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D19" ,"D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13", "D19" ,"D28"))
igg <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG", lantigen)
ig1 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG1", lantigen)
ig2 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG2", lantigen)
ig3 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG3", lantigen)
ig4 <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgG4", lantigen)
igm <- ff.figure(dataset[dataset$group%nin%"Naive",], "IgM", lantigen)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_after_challenge_within_status_all_v02.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
########################
## FOLD CHANGE
########################
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D13", "D19","D28"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D13", "D19" ,"D28"))
# compute fold change
dataset.wide <- dataset
dataset.wide$timepoint <- NULL
dataset.wide$log10_mfi <- NULL
dataset.wide <- dataset.wide[, c("original_id", "ig", "study_number","antigen", "mfi", "timepoint2c", "group")]
dataset.wide <- reshape(data = dataset.wide, direction="wide",
idvar = c("original_id", "ig", "antigen", "study_number", "group"),
v.names = c("mfi") ,
timevar = "timepoint2c")
dataset.wide$mfi.D7 <- dataset.wide[, "mfi.D7"] / dataset.wide[, "mfi.C-1"]
dataset.wide[, "mfi.D13"] <- dataset.wide[, "mfi.D13"] / dataset.wide[, "mfi.C-1"]
dataset.wide[, "mfi.D19"] <- dataset.wide[, "mfi.D19"] / dataset.wide[, "mfi.C-1"]
dataset.wide[, "mfi.D28"] <- dataset.wide[, "mfi.D28"] / dataset.wide[, "mfi.C-1"]
dataset.ratio <- reshape::melt(data=dataset.wide,
id.vars=c("original_id", "ig", "antigen", "study_number", "group"))
dataset.ratio <- subset(dataset.ratio, variable!="mfi.C-1")
dataset.ratio$timepoint2c <- gsub(pattern = "mfi.", replacement = "", x = dataset.ratio$variable, fixed=TRUE)
dataset.ratio <- subset(dataset.ratio, timepoint2c!="D19")
dataset.ratio$timepoint2c <- factor(dataset.ratio$timepoint2c, levels = c("D7", "D13", "D28"))
# function.to.plot
ff.figure <- function(dataset, isotype, lantigen, x.text, y.text){
ff.pval <- function(isotype, lantigen, timepointx){
ans <- ldply(lantigen, function(x){
ff.dat <- dataset[dataset$ig==isotype & dataset$timepoint2c==timepointx & dataset$antigen==x,]
ff.ans <- try(oneway.test(value ~ group, data=ff.dat)$p.value, silent=TRUE)
if(class(ff.ans)!="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx, antigen=x, raw.pvalue=ff.ans)
if(class(ff.ans)=="try-error") ans <- data.frame( isotype=isotype, timepoint2c=timepointx,antigen=x, raw.pvalue=NA)
return(ans)})
return(ans)
}
pvals <- ldply(c("D7", "D13", "D28") ,function(x) ff.pval(isotype, lantigen, x))
pvals$adj.pvalue <- p.adjust(pvals$raw.pvalue, method = "BH")
pvals$adj.pvalue <- sapply(pvals$adj.pvalue, function(x) {
if(!is.na(x)) ans <- fcn.pvalue.display(x,show.p=TRUE)
if(is.na(x)) ans <- "--"
return(ans)
})
pvals$adj.pvalue.text <- pvals$adj.pvalue
pvals$value <- y.text
pvals$group <- NA
ymax <- max(dataset[dataset$ig==isotype, "value"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, value, color=group))
q <- q + geom_boxplot(outlier.size = NA) + scale_y_log10(breaks=c(0.01,1,100,1000), labels=c(0.01,1,100,1000))
q <- q + geom_point(position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0), alpha=0.25)
q <- q + facet_wrap(~antigen) + scale_color_discrete(name = "Group",
breaks=c("Naive", "AA", "AS"),
labels=c("Naive", "AA", "AS"))
q <- q + theme_bw() + ggtitle(isotype) + theme(plot.title = element_text(hjust = 0.5))
q <- q + xlab("Timepoint") + ylab("MFI ratio (Timepoint/C-1)")
q <- q + geom_text(data=pvals, aes(x=timepoint2c, y=value, label=adj.pvalue.text), size=3)
# q <- q + ylim(log10(y.text) - (abs(log10(y.text))*0.85), log10(ymax))
ans <- list(q = q, pvals = pvals)
return(ans)
}
igg <- ff.figure(dataset.ratio[dataset.ratio$group%nin%"Naive",], "IgG", lantigen, 2, 0.1)
ig1 <- ff.figure(dataset.ratio[dataset.ratio$group%nin%"Naive",], "IgG1", lantigen, 2, 0.001)
ig2 <- ff.figure(dataset.ratio[dataset.ratio$group%nin%"Naive",], "IgG2", lantigen, 2, 0.01)
ig3 <- ff.figure(dataset.ratio[dataset.ratio$group%nin%"Naive",], "IgG3", lantigen, 2, 0.01)
ig4 <- ff.figure(dataset.ratio[dataset.ratio$group%nin%"Naive",], "IgG4", lantigen, 2, 0.001)
igm <- ff.figure(dataset.ratio[dataset.ratio$group%nin%"Naive",], "IgM", lantigen, 2, 0.01)
pdf(file=file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_after_challenge_within_timepoint_all_fold_change.pdf"),
width = 12, height = 8)
igg$q
ig1$q
ig2$q
ig3$q
ig4$q
igm$q
dev.off()
###########################################
###########################################
### Mixed models: NAIVE
###########################################
###########################################
library(nlme)
library(lmerTest)
dataset <- tmp.dat[tmp.dat$timepoint2c%in%c("D7", "D13", "D19", "D28","C-1"),]
dataset$newtimepoint <- with(dataset, ifelse(timepoint2c=="C-1",0,
ifelse(timepoint2c=="D7",7,
ifelse(timepoint2c=="D13",13,
ifelse(timepoint2c=="D19",19,
ifelse(timepoint2c=="D28",28,NA))))))
with(dataset, table(ig))
with(dataset, table(antigen))
with(dataset, table(dataset))
with(dataset, table(newtimepoint))
with(dataset, table(timepoint2c))
with(dataset, table(group))
### All models and isotype
library(multcomp)
ff1 <- function(model){
# Slope Naive
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0,0,0))))
a1.est <- round(a$test$coefficients, 3)
a1.pvalue <- as.numeric(a$test$pvalues)
# Slope AA
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0,1,0))))
a2.est <- round(a$test$coefficients, 3)
a2.pvalue <- as.numeric(a$test$pvalues)
a2.pvalue.int <- summary(model)$coefficients[5,5]
# Slope AS
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0,0,1))))
a3.est <- round(a$test$coefficients, 3)
a3.pvalue <- as.numeric(a$test$pvalues)
a3.pvalue.int <- summary(model)$coefficients[6,5]
# Summary
ans <- data.frame(a1.est, a1.pvalue, a2.est, a2.pvalue, a2.pvalue.int,
a3.est, a3.pvalue, a3.pvalue.int)
}
ans.dat <- list()
zz <- 0
for(i in c("IgG", "IgG1", "IgG2", "IgG3", "IgG4", "IgM") ){
z <- 0
ans1 <- list()
for(j in lantigen){
z <- z + 1
aux.data <- subset(dataset, ig==i & antigen==j) # INCLUDE BASELINE. THINK BASELINE
if(nrow(aux.data)!=0){
# All timepoints
mod <- lmer(log10_mfi ~ newtimepoint*group + (newtimepoint | original_id), data=aux.data)
ans <- try(ff1(mod), silent = TRUE)
if(inherits(ans, "try-error")){
ans <- data.frame(a1.est=NA, a1.pvalue=NA, a2.est=NA, a2.pvalue=NA, a2.pvalue.int=NA,
a3.est=NA, a3.pvalue=NA, a3.pvalue.int=NA)}
ans$isotype <- i
ans$antigen <- j
ans$model <- "three.timepoints"
ans1[[z]] <- ans
}
}
zz <- zz + 1
ans1 <- ldply(ans1)
ans1[, c(2,4,5,7,8)] <- sapply(c(2,4,5,7,8), function(x) fcn.pvalue.display(p.adjust(ans1[,x], method="BH")))
ans1 <- ans1[,c(9,10, 1:8)]
ans.dat[[zz]] <- ans1
}
ans.dat.summ <- ans.dat
#################################
### Export models to word
#################################
library(ReporteRs)
#########################################################
## Create Tables in Word describing
#########################################################
mydoc = docx(file.path(localpath, "Results", "Paper", "aim3", "mfi", "mixed_models.docx"))
xvector <- c("IgG", "IgG1", "IgG2", "IgG3", "IgG4","IgGM")
for(i in 1:6){
my_text = pot(paste0("Table", i, ". Mixed models (log10MFI response models) by antigen fitting the HbS (and Naive) and all timepoints, D7, D13, D19 and D84 (including C-1 information) for ", xvector[i] , " isotype. Coefficient and adjusted p value is shown (reference: Naive)."))
mydoc = addParagraph( mydoc, value = my_text, stylename = "Normal")
mydoc = addParagraph( mydoc, value = "", stylename = "Normal")
ftab = FlexTable(ans.dat.summ[[i]] ,add.rownames = FALSE, header.columns=FALSE,
body.par.props = parCenter(),
header.par.props = parCenter(),
body.text.props = textProperties( font.size = 6))
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","Slope Naive", "Slope AA", "Int. AA",
"Slope AS","Int. AS"),
colspan = c(1,1,2,2,1,2,1), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("Isotype", "Antigen","Coef", "p value", rep(c("Coef", "p value", "p value"),2)),
colspan = rep(1,10), par.properties = parCenter())
mydoc = addFlexTable( mydoc, ftab, par.properties = parCenter())
mydoc <- addPageBreak(mydoc)
}
writeDoc( mydoc, file = file.path(localpath, "Results", "Paper", "aim3", "mfi", "mixed_models.docx"))
###########################################
###########################################
### Mixed models: REMOVING NAIVE
###########################################
###########################################
library(nlme)
library(lmerTest)
dataset <- tmp.dat[tmp.dat$timepoint2c%in%c("D7", "D13", "D19", "D28","C-1"),]
dataset$newtimepoint <- with(dataset, ifelse(timepoint2c=="C-1",0,
ifelse(timepoint2c=="D7",7,
ifelse(timepoint2c=="D13",13,
ifelse(timepoint2c=="D19",19,
ifelse(timepoint2c=="D28",28,NA))))))
dataset <- subset(dataset, group!="Naive")
with(dataset, table(ig))
with(dataset, table(antigen))
with(dataset, table(dataset))
with(dataset, table(newtimepoint))
with(dataset, table(timepoint2c))
with(dataset, table(group))
### All models and isotype
library(multcomp)
ff1 <- function(model){
# Slope AA
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0))))
a1.est <- round(a$test$coefficients, 3)
a1.pvalue <- as.numeric(a$test$pvalues)
# Slope AS
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,1))))
a2.est <- round(a$test$coefficients, 3)
a2.pvalue <- as.numeric(a$test$pvalues)
a2.pvalue.int <- summary(model)$coefficients[4,5]
# Summary
ans <- data.frame(a1.est, a1.pvalue, a2.est, a2.pvalue, a2.pvalue.int)
}
ans.dat <- list()
zz <- 0
for(i in c("IgG", "IgG1", "IgG2", "IgG3", "IgG4", "IgM") ){
z <- 0
ans1 <- list()
for(j in lantigen){
z <- z + 1
aux.data <- subset(dataset, ig==i & antigen==j) # INCLUDE BASELINE. THINK BASELINE
if(nrow(aux.data)!=0){
# All timepoints
mod <- lmer(log10_mfi ~ newtimepoint*group + (newtimepoint | original_id), data=aux.data)
ans <- try(ff1(mod), silent = TRUE)
if(inherits(ans, "try-error")){
ans <- data.frame(a1.est=NA, a1.pvalue=NA, a2.est=NA, a2.pvalue=NA, a2.pvalue.int=NA,
a3.est=NA, a3.pvalue=NA, a3.pvalue.int=NA)}
ans$isotype <- i
ans$antigen <- j
ans$model <- "three.timepoints"
ans1[[z]] <- ans
}
}
zz <- zz + 1
ans1 <- ldply(ans1)
ans1[, c(2,4,5)] <- sapply(c(2,4,5), function(x) fcn.pvalue.display(p.adjust(ans1[,x], method="BH")))
ans1 <- ans1[,c(6,7, 1:5)]
ans.dat[[zz]] <- ans1
}
ans.dat.summ <- ans.dat
#################################
### Export models to word
#################################
library(ReporteRs)
#########################################################
## Create Tables in Word describing
#########################################################
mydoc = docx(file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_mixed_models.docx"))
xvector <- c("IgG", "IgG1", "IgG2", "IgG3", "IgG4","IgGM")
for(i in 1:6){
my_text = pot(paste0("Table", i, ". Mixed models (log10MFI response models) by antigen fitting the HbS (Naive group removed) and all timepoints, D7, D13, D19 and D84 (including C-1 information) for ", xvector[i] , " isotype. Coefficient and adjusted p value is shown (reference: AA)."))
mydoc = addParagraph( mydoc, value = my_text, stylename = "Normal")
mydoc = addParagraph( mydoc, value = "", stylename = "Normal")
ftab = FlexTable(ans.dat.summ[[i]] ,add.rownames = FALSE, header.columns=FALSE,
body.par.props = parCenter(),
header.par.props = parCenter(),
body.text.props = textProperties( font.size = 6))
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","Slope AA", "Slope AS", "Int. AS"),
colspan = c(1,1,2,2,1), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("Isotype", "Antigen","Coef", "p value", rep(c("Coef", "p value", "p value"),1)),
colspan = rep(1,7), par.properties = parCenter())
mydoc = addFlexTable( mydoc, ftab, par.properties = parCenter())
mydoc <- addPageBreak(mydoc)
}
writeDoc( mydoc, file = file.path(localpath, "Results", "Paper", "aim3", "mfi", "nonaive_mixed_models.docx"))
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.