projects/02_aim2/r_old/v_18.03.15_chmi_aim_02.R
In mvazquezs/chmitools: Package to manage & obtain results of the CHMI Study
##############################################################
### CODE TO ANALYZE CHMI AB DATA FOLLOWING THE PAPER: AIM2 ###
##############################################################
# 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))
## 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 <- NULL
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",
"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, 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=="T2") # only include Doses and Naive (Placebo)
tmp.dat <- subset(tmp.dat, antigen%nin%c("BSA", "agal"))
tmp.dat$Vaccine <- with(tmp.dat, ifelse(immune_status=="Naive", "Placebo",
ifelse(immune_status=="Vaccinated 12800", "12800",
ifelse(immune_status=="Vaccinated 3200", "3200",
ifelse(immune_status=="Vaccinated 51200", "51200","XXX")))))
with(tmp.dat, table(Vaccine, immune_status))
with(tmp.dat, table(dataset))
tmp.dat$Vaccine <- factor(tmp.dat$Vaccine, levels=c("Placebo", "3200", "12800", "51200"))
lantigen <- sort(unique(tmp.dat$antigen))
tmp.dat$timepoint2c <- tmp.dat$timepoint
tmp.dat$timepoint2c <- factor(tmp.dat$timepoint2c, levels=c("C-1", "D7", "D11", "D84"))
with(tmp.dat, table(timepoint2c, dataset))
########################################
########################################
########################################
### AIM 1 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 ~ Vaccine, 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(Vaccine, 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 + theme(axis.text.x = element_text(angle = 30, hjust = 1)) + xlab("Vaccine") + 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.5, 2)
ig1 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG1", lantigen, 2.5, 0.5)
ig2 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG2", lantigen, 2.5, 1)
ig3 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG3", lantigen, 2.5, 0.5)
ig4 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG4", lantigen, 2.5, -2)
igm <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgM", lantigen, 2.5, 1.5)
pdf(file=file.path(localpath, "Results", "Paper", "aim2", "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 ~ Vaccine, 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", "D11") ,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$Vaccine <- 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=Vaccine))
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 = "Vaccine",
breaks=c("Placebo", "3200", "12800", "51200"),
labels=c("Placebo", "3200", "12800", "51200"))
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("D7", "D11"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("D7", "D11"))
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", "aim2", "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 ~ Vaccine, 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", "D11") ,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$Vaccine <- 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=Vaccine))
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 = "Vaccine",
breaks=c("Placebo", "3200", "12800", "51200"),
labels=c("Placebo", "3200", "12800", "51200"))
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", "D11"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D11"))
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", "aim2", "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, Vaccine), summarise, med = median(log10_mfi))
summ.group$original_id <- NA
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=Vaccine, 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 = "Vaccine",
breaks=c("Placebo", "3200", "12800", "51200"),
labels=c("Placebo", "3200", "12800", "51200"))
q <- q + geom_point(data=summ.group, aes(timepoint2c, med, color=Vaccine), 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", "D11", "D84"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D11", "D84"))
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", "aim2", "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()
########################
## FOLD CHANGE
########################
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D11", "D84"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D11", "D84"))
# 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", "Vaccine")]
dataset.wide <- reshape(data = dataset.wide, direction="wide",
idvar = c("original_id", "ig", "antigen", "study_number", "Vaccine"),
v.names = c("mfi") ,
timevar = "timepoint2c")
dataset.wide$mfi.D7 <- dataset.wide[, "mfi.D7"] / dataset.wide[, "mfi.C-1"]
dataset.wide[, "mfi.D11"] <- dataset.wide[, "mfi.D11"] / dataset.wide[, "mfi.C-1"]
dataset.wide$mfi.D84 <- dataset.wide[, "mfi.D84"] / dataset.wide[, "mfi.C-1"]
dataset.ratio <- reshape::melt(data=dataset.wide,
id.vars=c("original_id", "ig", "antigen", "study_number", "Vaccine"))
dataset.ratio <- subset(dataset.ratio, variable!="mfi.C-1")
dataset.ratio$timepoint2c <- gsub(pattern = "mfi.", replacement = "", x = dataset.ratio$variable, fixed=TRUE)
# 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 ~ Vaccine, 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", "D11", "D84") ,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$Vaccine <- NA
ymax <- max(dataset[dataset$ig==isotype, "value"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, value, color=Vaccine))
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 = "Vaccine",
breaks=c("Placebo", "3200", "12800", "51200"),
labels=c("Placebo", "3200", "12800", "51200"))
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", "aim2", "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()
###########################
###########################
### Mixed models: FOLD CHANGE
###########################
###########################
library(nlme)
library(lmerTest)
dataset <- tmp.dat[tmp.dat$timepoint2c%in%c("D7", "D11", "D84", "C-1"),]
dataset$newtimepoint <- with(dataset, ifelse(timepoint2c=="C-1",0,
ifelse(timepoint2c=="D7",7,
ifelse(timepoint2c=="D11",11,
ifelse(timepoint2c=="D84",84,NA)))))
with(dataset, table(ig))
with(dataset, table(antigen))
with(dataset, table(dataset))
with(dataset, table(newtimepoint))
### All models and isotype
ff1 <- function(model){
# Slope Placebo
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0,0,0,0,0))))
a1.est <- round(a$test$coefficients, 3)
a1.pvalue <- as.numeric(a$test$pvalues)
# Slope 3200
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0,0,1,0,0))))
a2.est <- round(a$test$coefficients, 3)
a2.pvalue <- as.numeric(a$test$pvalues)
a2.pvalue.int <- summary(model)$coefficients[6,5]
# Slope 12800
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0,0,0,1,0))))
a3.est <- round(a$test$coefficients, 3)
a3.pvalue <- as.numeric(a$test$pvalues)
a3.pvalue.int <- summary(model)$coefficients[7,5]
# Slope 51200
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0,0,0,0,1))))
a4.est <- round(a$test$coefficients, 3)
a4.pvalue <- as.numeric(a$test$pvalues)
a4.pvalue.int <- summary(model)$coefficients[8,5]
# Summary
ans <- data.frame(a1.est, a1.pvalue, a2.est, a2.pvalue, a2.pvalue.int,
a3.est, a3.pvalue, a3.pvalue.int,
a4.est, a4.pvalue, a4.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 D7 vs D11 vs D84 (two status)
mod <- lmer(log10_mfi ~ newtimepoint*Vaccine + (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,a4.est=NA, a4.pvalue=NA, a4.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,10,11)] <- sapply(c(2,4,5,7,8,10,11), function(x) fcn.pvalue.display(p.adjust(ans1[,x], method="BH")))
ans1 <- ans1[,c(12,13, 1:11)]
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", "aim2", "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 three doses (and placebo) and three timepoints, D7, D11 and D84 (including C-1 information) for) ", xvector[i] , " isotype. Coefficient and adjusted p value is shown (reference: Placebo)."))
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 Placebo", "Slope 3200", "Int. 3200",
"Slope 12800","Int. 12800","Slope 51200", "Int. 51200"),
colspan = c(1,1,2,2,1,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"),3)),
colspan = rep(1,13), par.properties = parCenter())
mydoc = addFlexTable( mydoc, ftab, par.properties = parCenter())
mydoc <- addPageBreak(mydoc)
}
writeDoc( mydoc, file = file.path(localpath, "Results", "Paper", "aim2", "mfi", "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: AIM2 ###
##############################################################
# 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))
## 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 <- NULL
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",
"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, 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=="T2") # only include Doses and Naive (Placebo)
tmp.dat <- subset(tmp.dat, antigen%nin%c("BSA", "agal"))
tmp.dat$Vaccine <- with(tmp.dat, ifelse(immune_status=="Naive", "Placebo",
ifelse(immune_status=="Vaccinated 12800", "12800",
ifelse(immune_status=="Vaccinated 3200", "3200",
ifelse(immune_status=="Vaccinated 51200", "51200","XXX")))))
with(tmp.dat, table(Vaccine, immune_status))
with(tmp.dat, table(dataset))
tmp.dat$Vaccine <- factor(tmp.dat$Vaccine, levels=c("Placebo", "3200", "12800", "51200"))
lantigen <- sort(unique(tmp.dat$antigen))
tmp.dat$timepoint2c <- tmp.dat$timepoint
tmp.dat$timepoint2c <- factor(tmp.dat$timepoint2c, levels=c("C-1", "D7", "D11", "D84"))
with(tmp.dat, table(timepoint2c, dataset))
########################################
########################################
########################################
### AIM 1 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 ~ Vaccine, 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(Vaccine, 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 + theme(axis.text.x = element_text(angle = 30, hjust = 1)) + xlab("Vaccine") + 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.5, 2)
ig1 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG1", lantigen, 2.5, 0.5)
ig2 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG2", lantigen, 2.5, 1)
ig3 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG3", lantigen, 2.5, 0.5)
ig4 <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgG4", lantigen, 2.5, -2)
igm <- ff.figure(tmp.dat[tmp.dat$timepoint=="C-1",], "IgM", lantigen, 2.5, 1.5)
pdf(file=file.path(localpath, "Results", "Paper", "aim2", "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 ~ Vaccine, 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", "D11") ,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$Vaccine <- 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=Vaccine))
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 = "Vaccine",
breaks=c("Placebo", "3200", "12800", "51200"),
labels=c("Placebo", "3200", "12800", "51200"))
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("D7", "D11"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("D7", "D11"))
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", "aim2", "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 ~ Vaccine, 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", "D11") ,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$Vaccine <- 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=Vaccine))
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 = "Vaccine",
breaks=c("Placebo", "3200", "12800", "51200"),
labels=c("Placebo", "3200", "12800", "51200"))
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", "D11"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D11"))
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", "aim2", "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, Vaccine), summarise, med = median(log10_mfi))
summ.group$original_id <- NA
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, log10_mfi, color=Vaccine, 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 = "Vaccine",
breaks=c("Placebo", "3200", "12800", "51200"),
labels=c("Placebo", "3200", "12800", "51200"))
q <- q + geom_point(data=summ.group, aes(timepoint2c, med, color=Vaccine), 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", "D11", "D84"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D11", "D84"))
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", "aim2", "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()
########################
## FOLD CHANGE
########################
dataset <- tmp.dat[tmp.dat$timepoint%in%c("C-1","D7", "D11", "D84"),]
dataset$timepoint2c <- factor(as.character(dataset$timepoint2c), levels =c("C-1","D7", "D11", "D84"))
# 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", "Vaccine")]
dataset.wide <- reshape(data = dataset.wide, direction="wide",
idvar = c("original_id", "ig", "antigen", "study_number", "Vaccine"),
v.names = c("mfi") ,
timevar = "timepoint2c")
dataset.wide$mfi.D7 <- dataset.wide[, "mfi.D7"] / dataset.wide[, "mfi.C-1"]
dataset.wide[, "mfi.D11"] <- dataset.wide[, "mfi.D11"] / dataset.wide[, "mfi.C-1"]
dataset.wide$mfi.D84 <- dataset.wide[, "mfi.D84"] / dataset.wide[, "mfi.C-1"]
dataset.ratio <- reshape::melt(data=dataset.wide,
id.vars=c("original_id", "ig", "antigen", "study_number", "Vaccine"))
dataset.ratio <- subset(dataset.ratio, variable!="mfi.C-1")
dataset.ratio$timepoint2c <- gsub(pattern = "mfi.", replacement = "", x = dataset.ratio$variable, fixed=TRUE)
# 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 ~ Vaccine, 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", "D11", "D84") ,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$Vaccine <- NA
ymax <- max(dataset[dataset$ig==isotype, "value"], na.rm=TRUE) * 1.2
q <- ggplot(dataset[dataset$ig==isotype, ], aes(timepoint2c, value, color=Vaccine))
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 = "Vaccine",
breaks=c("Placebo", "3200", "12800", "51200"),
labels=c("Placebo", "3200", "12800", "51200"))
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", "aim2", "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()
###########################
###########################
### Mixed models: FOLD CHANGE
###########################
###########################
library(nlme)
library(lmerTest)
dataset <- tmp.dat[tmp.dat$timepoint2c%in%c("D7", "D11", "D84", "C-1"),]
dataset$newtimepoint <- with(dataset, ifelse(timepoint2c=="C-1",0,
ifelse(timepoint2c=="D7",7,
ifelse(timepoint2c=="D11",11,
ifelse(timepoint2c=="D84",84,NA)))))
with(dataset, table(ig))
with(dataset, table(antigen))
with(dataset, table(dataset))
with(dataset, table(newtimepoint))
### All models and isotype
ff1 <- function(model){
# Slope Placebo
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0,0,0,0,0))))
a1.est <- round(a$test$coefficients, 3)
a1.pvalue <- as.numeric(a$test$pvalues)
# Slope 3200
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0,0,1,0,0))))
a2.est <- round(a$test$coefficients, 3)
a2.pvalue <- as.numeric(a$test$pvalues)
a2.pvalue.int <- summary(model)$coefficients[6,5]
# Slope 12800
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0,0,0,1,0))))
a3.est <- round(a$test$coefficients, 3)
a3.pvalue <- as.numeric(a$test$pvalues)
a3.pvalue.int <- summary(model)$coefficients[7,5]
# Slope 51200
a <- summary(glht(model,linfct=rbind(b1=c(0,1,0,0,0,0,0,1))))
a4.est <- round(a$test$coefficients, 3)
a4.pvalue <- as.numeric(a$test$pvalues)
a4.pvalue.int <- summary(model)$coefficients[8,5]
# Summary
ans <- data.frame(a1.est, a1.pvalue, a2.est, a2.pvalue, a2.pvalue.int,
a3.est, a3.pvalue, a3.pvalue.int,
a4.est, a4.pvalue, a4.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 D7 vs D11 vs D84 (two status)
mod <- lmer(log10_mfi ~ newtimepoint*Vaccine + (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,a4.est=NA, a4.pvalue=NA, a4.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,10,11)] <- sapply(c(2,4,5,7,8,10,11), function(x) fcn.pvalue.display(p.adjust(ans1[,x], method="BH")))
ans1 <- ans1[,c(12,13, 1:11)]
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", "aim2", "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 three doses (and placebo) and three timepoints, D7, D11 and D84 (including C-1 information) for) ", xvector[i] , " isotype. Coefficient and adjusted p value is shown (reference: Placebo)."))
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 Placebo", "Slope 3200", "Int. 3200",
"Slope 12800","Int. 12800","Slope 51200", "Int. 51200"),
colspan = c(1,1,2,2,1,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"),3)),
colspan = rep(1,13), par.properties = parCenter())
mydoc = addFlexTable( mydoc, ftab, par.properties = parCenter())
mydoc <- addPageBreak(mydoc)
}
writeDoc( mydoc, file = file.path(localpath, "Results", "Paper", "aim2", "mfi", "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.