## ----echo=FALSE------------------------------------------------------------
.standalone=FALSE
## ---- message=FALSE, warning=FALSE, include=!.standalone-------------------
library("AnnotationDbi")
library("abind")
library("beeswarm")
library("Biobase")
library("biomaRt")
library("broom")
library("colorspace")
library("cowplot")
library("dendsort")
library("DESeq2")
library("doParallel")
library("dplyr")
library("foreach")
library("forestplot")
library("genefilter")
library("ggbeeswarm")
library("ggdendro")
library("ggplot2")
library("ggtern")
library("glmnet")
library("grid")
library("gridExtra")
library("gtable")
library("hexbin")
library("IHW")
library("ipflasso")
library("knitr")
library("limma")
library("magrittr")
library("maxstat")
library("nat")
library("org.Hs.eg.db")
library("BloodCancerMultiOmics2017")
library("pheatmap")
library("piano")
library("readxl")
library("RColorBrewer")
library("reshape2")
library("Rtsne")
library("scales")
library("SummarizedExperiment")
library("survival")
library("tibble")
library("tidyr")
library("tidyverse")
library("xtable")
## ---- message=FALSE, include=!exists(".standalone"), eval=!exists(".standalone")----
# library("BloodCancerMultiOmics2017")
# library("Biobase")
# library("ggplot2")
# library("gtable")
# library("grid")
# library("tidyverse")
# library("gridExtra")
## ----echo=FALSE------------------------------------------------------------
plotDir = ifelse(exists(".standalone"), "", "part01/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
## --------------------------------------------------------------------------
options(stringsAsFactors=FALSE)
## --------------------------------------------------------------------------
data("drpar", "drugs", "patmeta", "mutCOM")
## --------------------------------------------------------------------------
# PATIENTS
patM = colnames(drpar)
# DRUGS
drM = rownames(drpar)
drM = drM[!drM %in% "D_CHK"] # remove combintation of 2 drugs: D_CHK
## --------------------------------------------------------------------------
bwScale = c("0"="white","1"="black","N.A."="grey90")
lfsize = 16 # legend font size
## --------------------------------------------------------------------------
drugs$target_category = as.character(drugs$target_category)
drugs$group = NA
drugs$group[which(drugs$approved_042016==1)] = "FDA approved"
drugs$group[which(drugs$devel_042016==1)] = "clinical development/\ntool compound"
## ----echo=FALSE------------------------------------------------------------
res = table(drugs[,c("target_category","group")])
knitr::kable(res[order(res[,1], decreasing=TRUE),])
## ----echo=FALSE------------------------------------------------------------
goM = BloodCancerMultiOmics2017:::plotPathways(dat=drugs[drM,])
## ----dr_desc_M, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=goM[["figure"]][["width"]], fig.height=goM[["figure"]][["height"]]----
#FIG# 1C
grid.draw(goM[["figure"]][["plot"]])
## ---- echo=FALSE-----------------------------------------------------------
knitr::knit_hooks$set(crop=NULL)
## ----dr_desc_M_leg, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=goM[["legend"]][["width"]], fig.height=goM[["legend"]][["height"]], out.height=120, out.width=240----
#FIG# 1C
grid.draw(goM[["legend"]][["plot"]])
## ---- echo=FALSE-----------------------------------------------------------
knitr::knit_hooks$set(crop=knitr:::hook_pdfcrop)
## ----echo=FALSE------------------------------------------------------------
knitr::kable(data.frame(sort(table(patmeta[patM, "Diagnosis"]), decreasing=TRUE)))
## ----echo=FALSE------------------------------------------------------------
goM = BloodCancerMultiOmics2017:::plotPatientStat(pats=patM, gap=c(30,160))
## ----pt_desc_M, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=goM[["figure"]][["width"]], fig.height=goM[["figure"]][["height"]]----
#FIG# 1B
grid.draw(goM[["figure"]][["plot"]])
## ---- echo=FALSE-----------------------------------------------------------
knitr::knit_hooks$set(crop=NULL)
## ----pt_desc_M_leg, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=goM[["legend"]][["width"]], fig.height=goM[["legend"]][["height"]], out.height=120, out.width=240----
#FIG# 1B
grid.draw(goM[["legend"]][["plot"]])
## ---- echo=FALSE-----------------------------------------------------------
knitr::knit_hooks$set(crop=knitr:::hook_pdfcrop)
## --------------------------------------------------------------------------
# select CLL samples
patM = patM[patmeta[patM,"Diagnosis"]=="CLL"]
ighv = factor(setNames(patmeta[patM,"IGHV"], nm=patM), levels=c("U","M"))
mut1 = c("del17p13", "del11q22.3", "trisomy12", "del13q14_any")
mut2 = c("TP53", "ATM", "SF3B1", "NOTCH1", "MYD88")
mc = assayData(mutCOM)$binary[patM,]
## SELECTION OF MUTATIONS
# # include mutations with at least incidence of 4
mut2plot = names(which(sort(colSums(mc, na.rm=TRUE), decreasing=TRUE)>3))
# remove chromothrypsis
mut2plot = mut2plot[-grep("Chromothripsis", mut2plot)]
# divide mutations into gene mut and cnv
mut2plotSV = mut2plot[grep("[[:lower:]]", mut2plot)]
mut2plotSP = mut2plot[grep("[[:upper:]]", mut2plot)]
# remove some other things (it is quite manual thing, so be careful)
# IF YOU WANT TO REMOVE SOME MORE MUTATIONS JUST ADD THE LINES HERE!
mut2plotSV = mut2plotSV[-grep("del13q14_mono", mut2plotSV)]
mut2plotSV = mut2plotSV[-grep("del13q14_bi", mut2plotSV)]
mut2plotSV = mut2plotSV[-grep("del14q24.3", mut2plotSV)]
# rearrange the top ones to match the order in mut1 and mut2
mut2plotSV = c(mut1, mut2plotSV[!mut2plotSV %in% mut1])
mut2plotSP = c(mut2, mut2plotSP[!mut2plotSP %in% mut2])
factors = data.frame(assayData(mutCOM)$binary[patM, c(mut2plotSV, mut2plotSP)],
check.names=FALSE)
# change del13q14_any to del13q14
colnames(factors)[which(colnames(factors)=="del13q14_any")] = "del13q14"
mut2plotSV = gsub("del13q14_any", "del13q14", mut2plotSV)
# change it to factors
for(i in 1:ncol(factors)) {
factors[,i] = factor(factors[,i], levels=c(1,0))
}
ord = order(factors[,1], factors[,2], factors[,3], factors[,4], factors[,5],
factors[,6], factors[,7], factors[,8], factors[,9], factors[,10],
factors[,11], factors[,12], factors[,13], factors[,14],
factors[,15], factors[,16], factors[,17], factors[,18],
factors[,19], factors[,20], factors[,21], factors[,22],
factors[,23], factors[,24], factors[,25], factors[,26],
factors[,27], factors[,28], factors[,29], factors[,30],
factors[,31], factors[,32])
factorsord = factors[ord,]
patM = patM[ord]
(c(mut2plotSV, mut2plotSP))
## --------------------------------------------------------------------------
plotDF = meltWholeDF(factorsord)
plotDF$Mut =
ifelse(sapply(plotDF$X,
function(x) grep(x, list(mut2plotSV, mut2plotSP)))==1,"SV","SP")
plotDF$Status = "N.A."
plotDF$Status[plotDF$Measure==1 & plotDF$Mut=="SV"] = "1a"
plotDF$Status[plotDF$Measure==1 & plotDF$Mut=="SP"] = "1b"
plotDF$Status[plotDF$Measure==0] = "0"
plotDF$Status = factor(plotDF$Status, levels=c("1a","1b","0","N.A."))
plotDF$Y = factor(plotDF$Y, levels=patM)
plotDF$X = factor(plotDF$X, levels=rev(colnames(factorsord)))
mutPL = ggplotGrob(
ggplot(data=plotDF, aes(x=Y, y=X, fill=Status)) + geom_tile() +
scale_fill_manual(
values=c("0"="white","1a"="forestgreen","1b"="navy","N.A."="grey90"),
name="Mutation", labels=c("CNV","Gene mutation","WT","NA")) +
ylab("") + xlab("") +
geom_vline(xintercept=seq(0.5,length(patM)+1,5), colour="grey60") +
geom_hline(yintercept=seq(0.5,ncol(factorsord)+1,1), colour="grey60") +
scale_y_discrete(expand=c(0,0)) + scale_x_discrete(expand=c(0,0)) +
theme(axis.ticks=element_blank(), axis.text.x=element_blank(),
axis.text.y=element_text(
size=60, face=ifelse(levels(plotDF$X) %in% mut2plotSV,
"plain","italic")),
axis.text=element_text(margin=unit(0.5,"cm"), colour="black"),
legend.key = element_rect(colour = "black"),
legend.text=element_text(size=lfsize),
legend.title=element_text(size=lfsize)))
res = table(plotDF[,c("X","Measure")])
knitr::kable(res[order(res[,2], decreasing=TRUE),])
## --------------------------------------------------------------------------
ageDF = data.frame(Factor="Age",
PatientID=factor(patM, levels=patM),
Value=patmeta[patM,c("Age4Main")])
agePL = ggplotGrob(
ggplot(ageDF, aes(x=PatientID, y=Factor, fill=Value)) + geom_tile() +
scale_fill_gradient(low = "gold", high = "#3D1F00", na.value="grey92",
name="Age", breaks=c(40,60,80)) +
theme(axis.ticks=element_blank(),
axis.text=element_text(size=60, colour="black",
margin=unit(0.5,"cm")),
legend.text=element_text(size=lfsize),
legend.title=element_text(size=lfsize)))
hist(ageDF$Value, col="slategrey", xlab="Age", main="")
## --------------------------------------------------------------------------
sexDF = data.frame(Factor="Sex", PatientID=factor(patM, levels=patM),
Value=patmeta[patM, "Gender"])
sexPL = ggplotGrob(
ggplot(sexDF, aes(x=PatientID, y=Factor, fill=Value)) + geom_tile() +
scale_fill_manual(values=c("f"="maroon","m"="royalblue4","N.A."="grey90"),
name="Sex", labels=c("Female","Male","NA")) +
theme(axis.ticks=element_blank(),
axis.text=element_text(size=60, colour="black",
margin=unit(0.5,"cm")),
legend.key = element_rect(colour = "black"),
legend.text=element_text(size=lfsize),
legend.title=element_text(size=lfsize)))
table(sexDF$Value)
## --------------------------------------------------------------------------
treatDF = data.frame(Factor="Treated", PatientID=factor(patM, levels=patM),
Value=ifelse(patmeta[patM, "IC50beforeTreatment"], 0, 1))
treatDF$Value[is.na(treatDF$Value)] = "N.A."
treatDF$Value = factor(treatDF$Value, levels=c("0","1","N.A."))
treatPL = ggplotGrob(
ggplot(treatDF, aes(x=PatientID, y=Factor, fill=Value)) +geom_tile() +
scale_fill_manual(values=bwScale, name="Treated",
labels=c("0"="No","1"="Yes","N.A."="NA")) +
theme(axis.ticks=element_blank(),
axis.text=element_text(size=60, colour="black",
margin=unit(0.5,"cm")),
legend.key = element_rect(colour = "black"),
legend.text=element_text(size=lfsize),
legend.title=element_text(size=lfsize)))
table(treatDF$Value)
## --------------------------------------------------------------------------
ighvDF = data.frame(Factor="IGHV", PatientID=factor(patM, levels=patM),
Value=patmeta[patM, "IGHV"])
ighvDF$Value = ifelse(ighvDF$Value=="M", 1, 0)
ighvDF$Value[is.na(ighvDF$Value)] = "N.A."
ighvDF$Value = factor(ighvDF$Value, levels=c("0","1","N.A."))
ighvPL = ggplotGrob(
ggplot(ighvDF, aes(x=PatientID, y=Factor, fill=Value)) + geom_tile() +
scale_fill_manual(values=bwScale, name="IGHV",
labels=c("0"="Unmutated","1"="Mutated","N.A."="NA")) +
theme(axis.ticks=element_blank(),
axis.text=element_text(size=60, colour="black", margin=unit(0.5,"cm")),
legend.key=element_rect(colour = "black"),
legend.text=element_text(size=lfsize),
legend.title=element_text(size=lfsize)))
table(ighvDF$Value)
## ----echo=FALSE------------------------------------------------------------
nX = length(patM)
nY = ncol(factorsord)
unY1 = 0.6*1.6
unY2 = 0.6*1.8
unX = 0.2
sp = 0.001
wdths = c(6, unX*nX, sp)
hghts = c(sp, unY1,unY1,unY1,unY1, 0.8, sp, sp ,unY2*nY, sp)
gt = gtable(widths=unit(wdths, "in"), heights=unit(hghts, "in"))
# add the plots
gt = gtable_add_grob(gt, sexPL$grobs[[whichInGrob(sexPL, "panel")]], 2, 2)
gt = gtable_add_grob(gt, treatPL$grobs[[whichInGrob(treatPL, "panel")]], 3, 2)
gt = gtable_add_grob(gt, agePL$grobs[[whichInGrob(agePL, "panel")]], 4, 2)
gt = gtable_add_grob(gt, ighvPL$grobs[[whichInGrob(ighvPL, "panel")]], 5, 2)
gt = gtable_add_grob(gt, mutPL$grobs[[whichInGrob(mutPL, "panel")]], 9, 2)
# add x axis
gt = gtable_add_grob(gt, mutPL$grobs[[whichInGrob(mutPL, "axis-b")]], 10, 2)
# add y axis
gt = gtable_add_grob(gt, sexPL$grobs[[whichInGrob(sexPL, "axis-l")]], 2, 1)
gt = gtable_add_grob(gt, treatPL$grobs[[whichInGrob(treatPL, "axis-l")]], 3, 1)
gt = gtable_add_grob(gt, agePL$grobs[[whichInGrob(agePL, "axis-l")]], 4, 1)
gt = gtable_add_grob(gt, ighvPL$grobs[[whichInGrob(ighvPL, "axis-l")]], 5, 1)
gt = gtable_add_grob(gt, mutPL$grobs[[whichInGrob(mutPL, "axis-l")]], 9, 1)
## ----part1, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=sum(wdths), fig.height=sum(hghts)----
#FIG# 1D
grid.draw(gt)
## ----charLegend, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=10, fig.height=2----
#FIG# 1D
BloodCancerMultiOmics2017:::drawLegends(plobj=list(agePL,sexPL,treatPL,ighvPL,mutPL))
## --------------------------------------------------------------------------
data("lpdAll")
## --------------------------------------------------------------------------
plotTab = pData(lpdAll) %>%
transmute(x=log10(ATPday0), y=log10(ATP48h), diff=ATP48h/ATPday0) %>%
filter(!is.na(x))
## ---- fig.width=8, fig.height=7--------------------------------------------
lm_eqn <- function(df){
m <- lm(y ~ 1, df, offset = x)
ypred <- predict(m, newdata = df)
r2 = sum((ypred - df$y)^2)/sum((df$y - mean(df$y)) ^ 2)
eq <- substitute(italic(y) == italic(x) + a*","~~italic(r)^2~"="~r2,
list(a = format(coef(m)[1], digits = 2),
r2 = format(r2, digits = 2)))
as.character(as.expression(eq))
}
plotTab$ypred <- predict(lm(y~1,plotTab, offset = x), newdata = plotTab)
sca <- ggplot(plotTab, aes(x= x, y = y)) + geom_point(size=3) +
geom_smooth(data = plotTab, mapping = aes(x=x, y = ypred), method = lm, se = FALSE, formula = y ~ x) +
geom_text(x = 5.2, y = 6.2, label = lm_eqn(plotTab), parse = TRUE, size =8) +
xlab("log10(day0 ATP luminescence)") + ylab("log10(48h ATP luminescence)") +
theme_bw() + theme(axis.title = element_text(size = 15, face = "bold"),
axis.text = element_text(size=15), legend.position = "none") +
coord_cartesian(xlim = c(4.6,6.3), ylim = c(4.6,6.3))
## --------------------------------------------------------------------------
histo <- ggplot(plotTab, aes(x = diff)) + geom_histogram(col = "red", fill = "red", bins=30, alpha = 0.5) + theme_bw() +
theme(axis.title = element_text(size = 15, face = "bold"), axis.text = element_text(size=15), legend.position = "none") +
xlab("(48h ATP luminescence) / (day0 ATP luminescence)")
## ----ATPcount_combine, fig.path=plotDir, dev=c("png", "pdf"), fig.width=13, fig.height=6----
grid.arrange(sca, histo, ncol=2)
## --------------------------------------------------------------------------
data("day23rep")
## --------------------------------------------------------------------------
maxXY = 125
plottingDF = do.call(rbind, lapply(c("day2","day3"), function(day) {
tmp = merge(
meltWholeDF(assayData(day23rep)[[paste0(day,"rep1")]]),
meltWholeDF(assayData(day23rep)[[paste0(day,"rep2")]]),
by=c("X","Y"))
colnames(tmp) = c("PatientID", "DrugID", "ViabX", "ViabY")
tmp[,c("ViabX", "ViabY")] = tmp[,c("ViabX", "ViabY")] * 100
tmp$Day = ifelse(day=="day2", "48h", "72h")
tmp
}))
plottingDF$Shape =
ifelse(plottingDF$ViabX > maxXY | plottingDF$ViabY > maxXY, "B", "A")
## --------------------------------------------------------------------------
annotation =
do.call(rbind,
tapply(1:nrow(plottingDF),
paste(plottingDF$PatientID,
plottingDF$Day, sep="_"),
function(idx) {
data.frame(X=110, Y=10,
Shape="A",
PatientID=plottingDF$PatientID[idx[1]],
Day=plottingDF$Day[idx[1]],
Cor=cor(plottingDF$ViabX[idx],
plottingDF$ViabY[idx],
method="pearson"))
}))
## ----pilotRep, fig.path=plotDir, dev=c("png", "pdf"), fig.width=7, fig.height=5----
#FIG# S31
ggplot(data=plottingDF,
aes(x=ifelse(ViabX>maxXY,maxXY,ViabX), y=ifelse(ViabY>maxXY,maxXY,ViabY),
shape=Shape)) +
facet_grid(Day ~ PatientID) + theme_bw() +
geom_hline(yintercept=100, linetype="dashed",color="darkgrey") +
geom_vline(xintercept=100, linetype="dashed",color="darkgrey") +
geom_abline(intercept=0, slope=1, colour="grey") +
geom_point(size=1.5, alpha=0.6) +
scale_x_continuous(limits=c(0,maxXY), breaks=seq(0,maxXY,25)) +
scale_y_continuous(limits=c(0,maxXY), breaks=seq(0,maxXY,25)) +
xlab("% viability - replicate 1") + ylab("% viability - replicate 2") +
coord_fixed() + expand_limits(x = 0, y = 0) +
theme(axis.title.x=element_text(size = rel(1), vjust=-1),
axis.title.y=element_text(size = rel(1), vjust=1),
strip.background=element_rect(fill="gainsboro")) +
guides(shape=FALSE, size=FALSE) +
geom_text(data=annotation,
aes(x=X, y=Y, label=format(Cor, digits=2), size=1.2),
colour="maroon", hjust=0.2)
## ---- include=!exists(".standalone"), eval=!exists(".standalone")----------
# sessionInfo()
## ---- message=FALSE, warning=FALSE, include=FALSE--------------------------
rm(list=ls())
## ---- message=FALSE, include=!exists(".standalone"), eval=!exists(".standalone")----
# library("BloodCancerMultiOmics2017")
# library("reshape2") # melt
# library("Biobase")
# library("dplyr")
# library("RColorBrewer")
# library("ggplot2")
# library("ggdendro")
# library("gtable")
# library("grid")
# library("Rtsne")
# library("ggbeeswarm")
## ----echo=FALSE------------------------------------------------------------
plotDir = ifelse(exists(".standalone"), "", "part02/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
## --------------------------------------------------------------------------
options(stringsAsFactors=FALSE)
## --------------------------------------------------------------------------
data("lpdAll")
## --------------------------------------------------------------------------
#select drug screening data on patient samples
lpd <- lpdAll[fData(lpdAll)$type == "viab", pData(lpdAll)$Diagnosis != "hMNC"]
viabTab <- Biobase::exprs(lpd)
viabTab <- viabTab[,complete.cases(t(viabTab))]
viabTab <- reshape2::melt(viabTab)
viabTab$Concentration <- fData(lpd)[viabTab$Var1,"subtype"]
viabTab <- viabTab[viabTab$Concentration %in% c("1","2","3","4","5"),]
viabTab$drugName <- fData(lpd)[viabTab$Var1,"name"]
viabTab <- viabTab[order(viabTab$Concentration),]
#order drug by mean viablitity
drugOrder <- group_by(viabTab, drugName) %>%
summarise(meanViab = mean(value)) %>%
arrange(meanViab)
viabTab$drugName <- factor(viabTab$drugName, levels = drugOrder$drugName)
## ----ViabilityScatter_main, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=10, fig.height=5, warning=FALSE, out.width=560, out.height=280----
#FIG# S1
#Color for each concentration
colorCode <- rev(brewer.pal(7,"Blues")[3:7])
names(colorCode) <- unique(viabTab$Concentration)
ggplot(viabTab, aes(x=drugName,y=value, color=Concentration)) +
geom_jitter(size=1, na.rm = TRUE, alpha=0.8, shape =16) +
scale_color_manual(values = colorCode) +
ylab("Viability") + ylim(c(0,1.2)) + xlab("") +
guides(color = guide_legend(override.aes = list(size=3,alpha=1),
title = "concentration index")) +
theme_bw() +
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 90, hjust = 1, vjust=0.5),
legend.key = element_blank())
## --------------------------------------------------------------------------
data("drpar", "patmeta", "drugs")
## --------------------------------------------------------------------------
givePatientID4Diag = function(pts, diag=NA) {
pts = if(is.na(diag)) {
names(pts)
} else {
names(pts[patmeta[pts,"Diagnosis"]==diag])
}
pts
}
## --------------------------------------------------------------------------
giveViabMatrix = function(diag, screen, chnnl) {
data = if(screen=="main") drpar
else print("Incorrect screen name.")
pid = colnames(data)
if(!is.na(diag))
pid = pid[patmeta[pid,"Diagnosis"]==diag]
return(assayData(data)[[chnnl]][,pid])
}
## --------------------------------------------------------------------------
palette.cor1 = c(rev(brewer.pal(9, "Blues"))[1:8],
"white","white","white","white",brewer.pal(7, "Reds"))
palette.cor2 = c(rev(brewer.pal(9, "Blues"))[1:8],
"white","white","white","white",brewer.pal(7, "YlOrRd"))
## --------------------------------------------------------------------------
main.cll.tpll = BloodCancerMultiOmics2017:::makeCorrHeatmap(
mt=giveViabMatrix(diag="CLL", screen="main", chnnl="viaraw.4_5"),
mt2=giveViabMatrix(diag="T-PLL", screen="main", chnnl="viaraw.4_5"),
colsc=palette.cor2, concNo="one")
## ----main.CLL.T-PLL, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=main.cll.tpll[["figure"]][["width"]], fig.height=main.cll.tpll[["figure"]][["height"]]----
#FIG# 2A
grid.draw(main.cll.tpll[["figure"]][["plot"]])
## ----main.cll.tpll.legend, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=main.cll.tpll[["legend"]][["width"]], fig.height=main.cll.tpll[["legend"]][["height"]], out.width=300, out.height=150----
#FIG# 2A
grid.draw(main.cll.tpll[["legend"]][["plot"]])
## --------------------------------------------------------------------------
# select the data
mtcll = as.data.frame(t(giveViabMatrix(diag="CLL",
screen="main",
chnnl="viaraw.4_5")))
colnames(mtcll) = drugs[colnames(mtcll),"name"]
# function which plots the scatter plot
scatdr = function(drug1, drug2, coldot, mtNEW, min){
dataNEW = mtNEW[,c(drug1, drug2)]
colnames(dataNEW) = c("A", "B")
p = ggplot(data=dataNEW, aes(A, B)) + geom_point(size=3, col=coldot, alpha=0.8) +
labs(x = drug1, y = drug2) + ylim(c(min, 1.35)) + xlim(c(min, 1.35)) +
theme(panel.background = element_blank(),
axis.text = element_text(size = 15),
axis.title = element_text(size = rel(1.5)),
axis.line.x = element_line(colour = "black", size = 0.5),
axis.line.y = element_line(colour = "black", size = 0.5),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
geom_smooth(method=lm) +
geom_text(x=1, y=min+0.1,
label=paste0("Pearson-R = ",
round(cor(dataNEW$A, dataNEW$B ), 2)),
size = 5)
return(p)
}
## ----cor_scatter, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), , fig.width=6, fig.height=4, warning=FALSE, out.width=420, out.height=280----
#FIG# 2A
scatdr("ibrutinib", "spebrutinib", coldot="deeppink1", mtNEW=mtcll, min=0.4)
scatdr("ibrutinib", "PRT062607 HCl", coldot="deeppink1", mtNEW=mtcll, min=0.4)
scatdr("ibrutinib", "idelalisib", coldot="deeppink1", mtNEW=mtcll, min=0.4)
scatdr("venetoclax", "navitoclax", coldot="goldenrod2", mtNEW=mtcll, min=0.2)
scatdr("SD51", "MIS-43", coldot="dodgerblue3", mtNEW=mtcll, min=0.2)
## --------------------------------------------------------------------------
main.tpll = BloodCancerMultiOmics2017:::makeCorrHeatmap(
mt=giveViabMatrix(diag="T-PLL", screen="main", chnnl="viaraw.4_5"),
colsc=palette.cor1, concNo="one")
## ----main.T-PLL, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=main.tpll[["figure"]][["width"]], fig.height=main.tpll[["figure"]][["height"]]----
#FIG# S6 B
grid.draw(main.tpll[["figure"]][["plot"]])
## ----main.tpll.legend, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=main.tpll[["legend"]][["width"]], fig.height=main.tpll[["legend"]][["height"]], out.width=300, out.height=150----
#FIG# S6 B
grid.draw(main.tpll[["legend"]][["plot"]])
## --------------------------------------------------------------------------
main.mcl = BloodCancerMultiOmics2017:::makeCorrHeatmap(
mt=giveViabMatrix(diag="MCL", screen="main", chnnl="viaraw.4_5"),
colsc=palette.cor1, concNo="one")
## ----main.MCL, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=main.mcl[["figure"]][["width"]], fig.height=main.mcl[["figure"]][["height"]]----
#FIG# S6 A
grid.draw(main.mcl[["figure"]][["plot"]])
## ----main.mcl.legend, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=main.mcl[["legend"]][["width"]], fig.height=main.mcl[["legend"]][["height"]], out.width=300, out.height=150----
#FIG# S6 A
grid.draw(main.mcl[["legend"]][["plot"]])
## --------------------------------------------------------------------------
data(list=c("lpdAll", "conctab", "patmeta"))
## --------------------------------------------------------------------------
#Select rows contain drug response data
lpdSub <- lpdAll[fData(lpdAll)$type == "viab",]
#Only use samples with complete values
lpdSub <- lpdSub[,complete.cases(t(Biobase::exprs(lpdSub)))]
#Transformation of the values
Biobase::exprs(lpdSub) <- log(Biobase::exprs(lpdSub))
Biobase::exprs(lpdSub) <- t(scale(t(Biobase::exprs(lpdSub))))
#annotation for drug ID
anno <- sprintf("%s(%s)",fData(lpdSub)$name,fData(lpdSub)$subtype)
names(anno) <- rownames(lpdSub)
## --------------------------------------------------------------------------
tsneRun <- function(distMat,perplexity=10,theta=0,max_iter=5000, seed = 1000) {
set.seed(seed)
tsneRes <- Rtsne(distMat, perplexity = perplexity, theta = theta,
max_iter = max_iter, is_distance = TRUE, dims =2)
tsneRes <- tsneRes$Y
rownames(tsneRes) <- labels(distMat)
colnames(tsneRes) <- c("x","y")
tsneRes
}
## --------------------------------------------------------------------------
colDiagFill = c(`CLL` = "grey80",
`U-CLL` = "grey80",
`B-PLL`="grey80",
`T-PLL`="#cc5352",
`Sezary`="#cc5352",
`PTCL-NOS`="#cc5352",
`HCL`="#b29441",
`HCL-V`="mediumaquamarine",
`AML`="#addbaf",
`MCL`="#8e65ca",
`MZL`="#c95e9e",
`FL`="darkorchid4",
`LPL`="#6295cd",
`hMNC`="pink")
colDiagBorder <- colDiagFill
colDiagBorder["U-CLL"] <- "black"
## --------------------------------------------------------------------------
annoDiagNew <- function(patList, lpdObj = lpdSub) {
Diagnosis <- pData(lpdObj)[patList,c("Diagnosis","IGHV Uppsala U/M")]
DiagNew <- c()
for (i in seq(1:nrow(Diagnosis))) {
if (Diagnosis[i,1] == "CLL") {
if (is.na(Diagnosis[i,2])) {
DiagNew <- c(DiagNew,"CLL")
} else if (Diagnosis[i,2] == "U") {
DiagNew <- c(DiagNew,sprintf("%s-%s",Diagnosis[i,2],Diagnosis[i,1]))
} else if (Diagnosis[i,2] == "M") {
DiagNew <- c(DiagNew,"CLL")
}
} else DiagNew <- c(DiagNew,Diagnosis[i,1])
}
DiagNew
}
## --------------------------------------------------------------------------
#prepare distance matrix
distLpd <- dist(t(Biobase::exprs(lpdSub)))
#run t-SNE
plotTab <- data.frame(tsneRun(distLpd,perplexity=25, max_iter=5000, seed=338))
#annotated patient sample
plotTab$Diagnosis <- pData(lpdSub[,rownames(plotTab)])$Diagnosis
plotTab$Diagnosis <- annoDiagNew(rownames(plotTab,lpdSub)) #consider IGHV status
plotTab$Diagnosis <- factor(plotTab$Diagnosis,levels = names(colDiagFill))
## ----tSNE, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=10, fig.height=8, warning=FALSE, out.width=600, out.height=480----
#FIG# 2 C
p <- (ggplot(plotTab, aes(x=x,y=y)) +
geom_point(size=3, shape= 21, aes(col = Diagnosis, fill = Diagnosis)) +
theme_classic() +
theme(axis.ticks=element_line(color="black",size=0.5),
text=element_text(size=20),
axis.line.x = element_line(color="black",size=0.5),
axis.line.y = element_line(color="black",size=0.5),
legend.position="right") +
scale_fill_manual(values = colDiagFill) +
scale_color_manual(values = colDiagBorder) +
xlab("Component 1") + ylab("Component 2")) +
coord_cartesian(xlim = c(-20,20),ylim=c(-20,20))
print(p)
## --------------------------------------------------------------------------
lpdPlot <- lpdAll[fData(lpdAll)$type == "viab",]
concList <- c()
for (drugID in rownames(fData(lpdPlot))) {
concIndex <- as.character(fData(lpdPlot)[drugID,"subtype"])
concSplit <- unlist(strsplit(as.character(concIndex),":"))
ID <- substr(drugID,1,5)
if (length(concSplit) == 1) {
realConc <- conctab[ID,as.integer(concSplit)]
concList <- c(concList,realConc)
} else {
realConc <- sprintf("%s:%s",
conctab[ID,as.integer(concSplit[1])],
conctab[ID,as.integer(concSplit[2])])
concList <- c(concList,realConc)
}
}
fData(lpdPlot)$concValue <- concList
lpdPlot <- lpdPlot[,complete.cases(t(Biobase::exprs(lpdPlot)))]
## --------------------------------------------------------------------------
patDiag <- c("CLL","T-PLL","HCL","MCL")
drugID <- c("D_012_5","D_017_4","D_039_3","D_040_5","D_081_4","D_083_5")
lpdBee <- lpdPlot[drugID,pData(lpdPlot)$Diagnosis %in% patDiag]
## --------------------------------------------------------------------------
lpdCurve <-
lpdPlot[fData(lpdPlot)$name %in% fData(lpdBee)$name,
pData(lpdPlot)$Diagnosis %in% patDiag]
lpdCurve <- lpdCurve[fData(lpdCurve)$subtype %in% seq(1,5),]
dataCurve <- data.frame(Biobase::exprs(lpdCurve))
dataCurve <- cbind(dataCurve,fData(lpdCurve)[,c("name","concValue")])
tabCurve <- melt(dataCurve,
id.vars = c("name","concValue"), variable.name = "patID")
tabCurve$Diagnosis <- factor(pData(lpdCurve[,tabCurve$patID])$Diagnosis,
levels = patDiag)
tabCurve$value <- tabCurve$value
tabCurve$concValue <- as.numeric(tabCurve$concValue)
# set order
tabCurve$name <- factor(tabCurve$name, levels = fData(lpdBee)$name)
#calculate the mean and mse for each drug+cencentration in different disease
tabGroup <- group_by(tabCurve,name,concValue,Diagnosis)
tabSum <- summarise(tabGroup,meanViab = mean(value))
## ----viabilityCurve, fig.path=plotDir, dev=c("png", "pdf"), fig.width=4, fig.height=3----
#FIG# 2 C
tconc = expression("Concentration [" * mu * "M]")
fmt_dcimals <- function(decimals=0){
# return a function responpsible for formatting the
# axis labels with a given number of decimals
function(x) as.character(round(x,decimals))
}
for (drugName in unique(tabSum$name)) {
tabDrug <- filter(tabSum, name == drugName)
p <- (ggplot(data=tabDrug, aes(x=concValue,y=meanViab, col=Diagnosis)) +
geom_line() + geom_point(pch=16, size=4) +
scale_color_manual(values = colDiagFill[patDiag])
+ theme_classic() +
theme(panel.border=element_blank(),
axis.line.x=element_line(size=0.5,
linetype="solid", colour="black"),
axis.line.y = element_line(size = 0.5,
linetype="solid", colour="black"),
legend.position="none",
plot.title = element_text(hjust = 0.5, size=20),
axis.text = element_text(size=15),
axis.title = element_text(size=20)) +
ylab("Viability") + xlab(tconc) + ggtitle(drugName) +
scale_x_log10(labels=fmt_dcimals(2)) +
scale_y_continuous(limits = c(0,1.3), breaks = seq(0,1.3,0.2)))
plot(p)
}
## ----viabilityBee, fig.path=plotDir, dev=c("png", "pdf"), fig.width=5, fig.height=10, warning=FALSE----
#FIG# 2 D
lpdDiag <- lpdAll[,pData(lpdAll)$Diagnosis %in% c("CLL", "MCL", "HCL", "T-PLL")]
dr <- c("D_012_5", "D_083_5", "D_081_3", "D_040_4", "D_039_3")
m <- data.frame(t(Biobase::exprs(lpdDiag)[dr, ]), diag=pData(lpdDiag)$Diagnosis)
m <- melt(m)
m$lable <- 1
for (i in 1:nrow(m )) {
m[i, "lable"] <- giveDrugLabel(as.character(m[i, "variable"]), conctab, drugs)
}
ggplot( m, aes(diag, value, color=factor(diag) ) ) +
ylim(0,1.3) + ylab("Viability") +
xlab("") +
geom_boxplot(outlier.shape = NA) +
geom_beeswarm(cex=1.4, size=1.4,alpha=0.5, color="grey80") +
scale_color_manual("diagnosis", values=c(colDiagFill["CLL"], colDiagFill["MCL"],
colDiagFill["HCL"], colDiagFill["T-PLL"])) +
theme_bw() +
theme(legend.position="right") +
theme(
panel.background = element_blank(),
panel.grid.minor.x = element_blank(),
axis.text = element_text(size=15),
axis.title = element_text(size=15),
strip.text = element_text(size=15)
) +
facet_wrap(~ lable, ncol=1)
## ---- include=!exists(".standalone"), eval=!exists(".standalone")----------
# sessionInfo()
## ---- message=FALSE, warning=FALSE, include=FALSE--------------------------
rm(list=ls())
## ---- message=FALSE, include=!exists(".standalone"), eval=!exists(".standalone")----
# library("BloodCancerMultiOmics2017")
# library("readxl")
# library("dplyr")
# library("ggplot2")
# library("reshape2")
# library("xtable")
## ----echo=FALSE------------------------------------------------------------
plotDir = ifelse(exists(".standalone"), "", "part14/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
## ---- echo=FALSE-----------------------------------------------------------
# make trasperent
makeTransparent = function(..., alpha=0.18) {
if(alpha<0 | alpha>1) stop("alpha must be between 0 and 1")
alpha = floor(255*alpha)
newColor = col2rgb(col=unlist(list(...)), alpha=FALSE)
.makeTransparent = function(col, alpha) {
rgb(red=col[1], green=col[2], blue=col[3], alpha=alpha, maxColorValue=255)
}
newColor = apply(newColor, 2, .makeTransparent, alpha=alpha)
return(newColor)
}
## ----load_affinity_data----------------------------------------------------
# AZD7762 binding affinity constants
azd = read_excel(system.file("extdata","TargetProfiling.xlsx",
package="BloodCancerMultiOmics2017"), sheet = 1)
# PF477736 binding affinity constants
pf = read_excel(system.file("extdata","TargetProfiling.xlsx",
package="BloodCancerMultiOmics2017"), sheet = 2)
# BCR tagets Proc Natl Acad Sci U S A. 2016 May 17;113(20):5688-93
pProt = read_excel(system.file("extdata","TargetProfiling.xlsx",
package="BloodCancerMultiOmics2017"),sheet = 3)
## --------------------------------------------------------------------------
p <- full_join(azd, pf )
p <- full_join(p, pProt )
pp <- p[p$BCR_effect=="Yes",]
pp <- data.frame(pp[-which(is.na(pp$BCR_effect)),])
## ----azd-pf, fig.path=plotDir, dev=c('png','pdf'), fig.width=7, fig.height=4----
#FIG# 2B
rownames(pp) <- 1:nrow(pp)
pp <- as.data.frame(pp)
pp <- melt(pp)
colnames(pp)[3] <- "Drugs"
colnames(pp)[4] <- "Score"
ggplot(pp, aes(x= reorder(gene, Score), Score, colour=Drugs ) )+ geom_point(size=3) +
scale_colour_manual(values = c(makeTransparent("royalblue1", alpha = 0.75),
makeTransparent("royalblue4", alpha = 0.75),
makeTransparent("brown1", alpha = 0.55),
makeTransparent("brown3", alpha = 0.35)),
breaks = c("az10", "az2", "pf10", "pf2"),
labels = c("AZD7762 10 µM","AZD7762 2 µM","PF477736 10 µM","PF477736 2 µM") ) +
ylab("Binding affinity") +
theme_bw() + geom_hline(yintercept = 0.5) +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
axis.title.x=element_blank() )
## ----kinobead--------------------------------------------------------------
j <- apply(p[,c("az10", "az2", "pf10", "pf2")], 1, function (x) { min(x, na.rm=FALSE) } )
p <- p[which(j<0.5), ]
p <- unique(p, by = p$gene)
knitr::kable(p)
## ----kinobead_write_file, echo=FALSE, results='hide', eval=FALSE-----------
# write(print(p), file=paste0(plotDir,"kinobead.tex"))
## ---- include=!exists(".standalone"), eval=!exists(".standalone")----------
# sessionInfo()
## ---- message=FALSE, warning=FALSE, include=FALSE--------------------------
rm(list=ls())
## ----setup, message=FALSE, include=!exists(".standalone"), eval=!exists(".standalone")----
# options(error = recover)
# knitr::opts_chunk$set(tidy = FALSE, # cache = TRUE, autodep = TRUE,
# message = FALSE, error = FALSE, warning = TRUE)
#
# library("BloodCancerMultiOmics2017")
# library("Biobase")
# library("dplyr")
# library("RColorBrewer")
# library("dendsort")
# library("nat") # provides nlapply
# library("grid")
# library("magrittr")
## ----echo=FALSE------------------------------------------------------------
plotDir = ifelse(exists(".standalone"), "", "part03/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
## ----load------------------------------------------------------------------
data("conctab", "lpdAll", "drugs", "patmeta")
lpdCLL <- lpdAll[, lpdAll$Diagnosis=="CLL"]
lpdAll = lpdAll[, lpdAll$Diagnosis!="hMNC"]
## --------------------------------------------------------------------------
someMatch <- function(...) {
rv <- match(...)
if (all(is.na(rv)))
stop(sprintf("`match` failed to match any of the following: %s",
paste(x[is.na(rv)], collapse=", ")))
rv
}
## ----geneusage-------------------------------------------------------------
colnames(pData(lpdAll))
gu <- pData(lpdAll)$`IGHV Uppsala gene usage`
tabgu <- sort(table(gu), decreasing = TRUE)
biggestGroups <- names(tabgu)[1:5]
gu[ is.na(match(gu, biggestGroups)) & !is.na(gu) ] <- "other"
pData(lpdAll)$`IGHV gene usage` <- factor(gu, levels = c(biggestGroups, "other"))
## ----targetedDrugNames-----------------------------------------------------
stopifnot(is.null(drugs$id))
drugs$id <- rownames(drugs)
targetedDrugNames <- c("ibrutinib", "idelalisib", "PRT062607 HCl",
"duvelisib", "spebrutinib", "selumetinib", "MK-2206",
"everolimus", "encorafenib")
id1 <- safeMatch(targetedDrugNames, drugs$name)
targetedDrugs <- paste( rep(drugs[id1, "id"], each = 2), 4:5, sep="_" )
chemoDrugNames <- c("fludarabine", "doxorubicine", "nutlin-3")
id2 <- safeMatch(chemoDrugNames, drugs$name)
chemoDrugs <- paste( rep(drugs[id2, "id"], each = 5), 3:5, sep="_" )
tzselDrugNames <- c("ibrutinib", "idelalisib", "duvelisib", "selumetinib",
"AZD7762", "MK-2206", "everolimus", "venetoclax", "thapsigargin",
"AT13387", "YM155", "encorafenib", "tamatinib", "ruxolitinib",
"PF 477736", "fludarabine", "nutlin-3")
id3 <- safeMatch(tzselDrugNames, drugs$name)
tzselDrugs <- unlist(lapply(seq(along = tzselDrugNames), function(i)
paste(drugs[id3[i], "id"],
if (tzselDrugNames[i] %in% c("fludarabine", "nutlin-3")) 2:3 else 4:5,
sep = "_" )))
## ----addBCR----------------------------------------------------------------
bcrDrugs <- c("ibrutinib", "idelalisib", "PRT062607 HCl", "spebrutinib")
everolID <- drugs$id[ safeMatch("everolimus", drugs$name)]
bcrID <- drugs$id[ safeMatch(bcrDrugs, drugs$name)]
is_BCR <-
(fData(lpdAll)$id %in% bcrID) & (fData(lpdAll)$subtype %in% paste(4:5))
is_mTOR <-
(fData(lpdAll)$id %in% everolID) & (fData(lpdAll)$subtype %in% paste(4:5))
myin <- function(x, y) as.numeric( (x %in% y) & !is.na(x) )
weights1 <- data.frame(
hclust = rep(1, nrow(lpdAll)) + 1.75 * is_mTOR,
score = as.numeric( is_BCR ),
row.names = rownames(lpdAll))
weights2 <- data.frame(
row.names = tzselDrugs,
hclust = myin(drugs$target_category[id3], "B-cell receptor") * 0.3 + 0.7,
score = rep(1, length(tzselDrugs)))
## ----badDrugs--------------------------------------------------------------
badDrugs <- c(bortezomib = "D_008", `NSC 74859` = "D_025")
stopifnot(identical(drugs[ badDrugs, "name"], names(badDrugs)))
candDrugs <- rownames(lpdAll)[
fData(lpdAll)$type=="viab" & !(fData(lpdAll)$id %in% badDrugs) &
fData(lpdAll)$subtype %in% paste(2:5)
]
## ----threshs---------------------------------------------------------------
thresh <- list(effectVal = 0.7, effectNum = 4, viab = 0.6, maxval = 1.1)
overallMean <- rowMeans(Biobase::exprs(lpdAll)[candDrugs, ])
nthStrongest <- apply(Biobase::exprs(lpdAll)[candDrugs, ], 1,
function(x) sort(x)[thresh$effectNum])
## ----hists, fig.width = 8, fig.height = 2.6--------------------------------
par(mfrow = c(1, 3))
hist(overallMean, breaks = 30, col = "pink")
abline(v = thresh$viab, col="blue")
hist(nthStrongest, breaks = 30, col = "pink")
abline(v = thresh$effectVal, col="blue")
plot(overallMean, nthStrongest)
abline(h = thresh$effectVal, v = thresh$viab, col = "blue")
## ----seldrugs1-------------------------------------------------------------
seldrugs1 <- candDrugs[ overallMean >= thresh$viab &
nthStrongest <= thresh$effectVal ] %>%
union(targetedDrugs) %>%
union(chemoDrugs)
d1 <- Biobase::exprs(lpdAll[seldrugs1,, drop = FALSE ]) %>%
deckel(lower = 0, upper = thresh$maxval)
d2 <- Biobase::exprs(lpdAll[tzselDrugs,, drop = FALSE ]) %>%
deckel(lower = 0, upper = thresh$maxval)
## ----differentspread, fig.width = 4, fig.height = 4------------------------
spreads <- sapply(list(mad = mad, `Q95-Q05` = function(x)
diff(quantile(x, probs = c(0.05, 0.95)))), function(s) apply(d1, 1, s))
plot( spreads )
jj <- names(which( spreads[, "mad"] < 0.15 & spreads[, "Q95-Q05"] > 0.7))
jj
drugs[ stripConc(jj), "name" ]
## ----transf----------------------------------------------------------------
medianCenter_MadScale <- function(x) {
s <- median(x)
(x - s) / deckel(mad(x, center = s), lower = 0.05, upper = 0.2)
}
scaleDrugResp <- function(x) t(apply(x, 1, medianCenter_MadScale))
scd1 <- scaleDrugResp(d1)
scd2 <- scaleDrugResp(d2)
## ----defdisgrp-------------------------------------------------------------
sort(table(lpdAll$Diagnosis), decreasing = TRUE)
diseaseGroups <- list(
`CLL` = c("CLL"),
`MCL` = c("MCL"),
`HCL` = c("HCL", "HCL-V"),
`other B-cell` = c("B-PLL", "MZL", "LPL", "FL"),
`T-cell` = c("T-PLL", "Sezary", "PTCL-NOS"),
`myeloid` = c("AML"))
stopifnot(setequal(unlist(diseaseGroups), unique(lpdAll$Diagnosis)))
fdg <- factor(rep(NA, ncol(lpdAll)), levels = names(diseaseGroups))
for (i in names(diseaseGroups))
fdg[ lpdAll$Diagnosis %in% diseaseGroups[[i]] ] <- i
lpdAll$`Disease Group` <- fdg
## ----matclust--------------------------------------------------------------
matClust <- function(x,
rowweights,
colgroups = factor(rep("all", ncol(x))),
reorderrows = FALSE) {
stopifnot(is.data.frame(rowweights),
c("hclust", "score") %in% colnames(rowweights),
!is.null(rownames(rowweights)),
!is.null(rownames(x)),
all(rownames(x) %in% rownames(rowweights)),
is.factor(colgroups),
!any(is.na(colgroups)),
length(colgroups) == ncol(x))
wgt <- rowweights[ rownames(x), ]
columnsClust <- function(xk) {
score <- -svd(xk * wgt[, "score"])$v[, 1]
cmns <- colSums(xk * wgt[, "score"])
## make sure that high score = high sensitivity
if (cor(score, cmns) > 0) score <- (-score)
ddraw <- as.dendrogram(hclust(dist(t(xk * wgt[, "hclust"]),
method = "euclidean"),
method = "complete"))
dd <- reorder(ddraw, wts = -score, agglo.FUN = mean)
ord <- order.dendrogram(dd)
list(dd = dd, ord = ord, score = score)
}
sp <- split(seq(along = colgroups), colgroups)
cc <- lapply(sp, function(k) columnsClust(x[, k, drop=FALSE]))
cidx <- unlist(lapply(seq(along = cc), function (i)
sp[[i]][ cc[[i]]$ord ]))
csc <- unlist(lapply(seq(along = cc), function (i)
cc[[i]]$score[ cc[[i]]$ord ]))
rddraw <- as.dendrogram(hclust(dist(x, method = "euclidean"),
method = "complete"))
ridx <- if (reorderrows) {
ww <- (colgroups == "CLL")
stopifnot(!any(is.na(ww)), any(ww))
rowscore <- svd(t(x) * ww)$v[, 1]
dd <- reorder(rddraw, wts = rowscore, agglo.FUN = mean)
order.dendrogram(dd)
} else {
rev(order.dendrogram(dendsort(rddraw)))
}
res <- x[ridx, cidx]
stopifnot(identical(dim(res), dim(x)))
attr(res, "colgap") <- cumsum(sapply(cc, function(x) length(x$score)))
res
}
## --------------------------------------------------------------------------
translation = list(IGHV=c(U=0, M=1),
Methylation_Cluster=c(`LP-CLL`=0, `IP-CLL`=1, `HP-CLL`=2))
## ----annosamp--------------------------------------------------------------
make_pd <- function(cn, ...) {
df <- function(...) data.frame(..., check.names = FALSE)
x <- lpdAll[, cn]
pd <- df(
t(Biobase::exprs(x)[ c("del17p13", "TP53", "trisomy12"), , drop = FALSE]) %>%
`colnames<-`(c("del 17p13", "TP53", "trisomy 12")))
# pd <- df(pd,
# t(Biobase::exprs(x)[ c("SF3B1", "del11q22.3", "del13q14_any"),, drop = FALSE]) %>%
# `colnames<-`(c("SF3B1", "del11q22.3", "del13q14")))
pd <- df(pd,
cbind(as.integer(Biobase::exprs(x)["KRAS",] | Biobase::exprs(x)["NRAS",])) %>%
`colnames<-`("KRAS | NRAS"))
pd <- df(pd,
# IGHV = Biobase::exprs(x)["IGHV Uppsala U/M", ],
`IGHV (%)` = cut(x[["IGHV Uppsala % SHM"]],
breaks = c(0, seq(92, 100, by = 2), Inf), right = FALSE),
`Meth. Cluster` = names(translation$Methylation_Cluster)[
someMatch(paste(Biobase::exprs(x)["Methylation_Cluster", ]),
translation$Methylation_Cluster)],
`Gene usage` = x$`IGHV gene usage`)
if(length(unique(x$Diagnosis)) > 1)
pd <- df(pd, Diagnosis = x$Diagnosis)
pd <- df(pd,
pretreated = ifelse(patmeta[colnames(x),"IC50beforeTreatment"],"no","yes"),
alive = ifelse(patmeta[colnames(x),"died"]>0, "no", "yes"),
sex = factor(x$Gender))
rownames(pd) <- colnames(Biobase::exprs(x))
for (i in setdiff(colnames(pd), "BCR score")) {
if (!is.factor(pd[[i]]))
pd[[i]] <- factor(pd[[i]])
if (any(is.na(pd[[i]]))) {
levels(pd[[i]]) <- c(levels(pd[[i]]), "n.d.")
pd[[i]][ is.na(pd[[i]]) ] <- "n.d."
}
}
pd
}
## ----annokey---------------------------------------------------------------
gucol <- rev(brewer.pal(nlevels(lpdAll$`IGHV gene usage`), "Set3")) %>%
`names<-`(sort(levels(lpdAll$`IGHV gene usage`)))
gucol["IGHV3-21"] <- "#E41A1C"
make_ann_colors <- function(pd) {
bw <- c(`TRUE` = "darkblue", `FALSE` = "#ffffff")
res <- list(
Btk = bw, Syk = bw, PI3K = bw, MEK = bw)
if ("exptbatch" %in% colnames(pd))
res$exptbatch <- brewer.pal(nlevels(pd$exptbatch), "Set2") %>%
`names<-`(levels(pd$exptbatch))
if ("IGHV (%)" %in% colnames(pd))
res$`IGHV (%)` <-
c(rev(colorRampPalette(
brewer.pal(9, "Blues"))(nlevels(pd$`IGHV (%)`)-1)), "white") %>%
`names<-`(levels(pd$`IGHV (%)`))
if ("CD38" %in% colnames(pd))
res$CD38 <- colorRampPalette(
c("blue", "yellow"))(nlevels(pd$CD38)) %>% `names<-`(levels(pd$CD38))
if("Gene usage" %in% colnames(pd))
res$`Gene usage` <- gucol
if("Meth. Cluster" %in% colnames(pd))
res$`Meth. Cluster` <- brewer.pal(9, "Blues")[c(1, 5, 9)] %>%
`names<-`(names(translation$Methylation_Cluster))
res <- c(res, BloodCancerMultiOmics2017:::sampleColors) # from addons.R
if("Diagnosis" %in% colnames(pd))
res$Diagnosis <- BloodCancerMultiOmics2017:::colDiagS[
names(BloodCancerMultiOmics2017:::colDiagS) %in% levels(pd$Diagnosis) ] %>%
(function(x) x[order(names(x))])
for(i in names(res)) {
whnd <- which(names(res[[i]]) == "n.d.")
if(length(whnd)==1)
res[[i]][whnd] <- "#e0e0e0" else
res[[i]] <- c(res[[i]], `n.d.` = "#e0e0e0")
stopifnot(all(pd[[i]] %in% names(res[[i]])))
}
res
}
## ----dm--------------------------------------------------------------------
theatmap <- function(x, cellwidth = 7, cellheight = 11) {
stopifnot(is.matrix(x))
patDat <- make_pd(colnames(x))
bpp <- brewer.pal(9, "Set1")
breaks <- 2.3 * c(seq(-1, 1, length.out = 101)) %>% `names<-`(
colorRampPalette(c(rev(brewer.pal(7, "YlOrRd")),
"white", "white", "white",
brewer.pal(7, "Blues")))(101))
if (!is.null(attr(x, "colgap")))
stopifnot(last(attr(x, "colgap")) == ncol(x))
pheatmapwh(deckel(x, lower = first(breaks), upper = last(breaks)),
cluster_rows = FALSE,
cluster_cols = FALSE,
gaps_col = attr(x, "colgap"),
gaps_row = attr(x, "rowgap"),
scale = "none",
annotation_col = patDat,
annotation_colors = make_ann_colors(patDat),
color = names(breaks),
breaks = breaks,
show_rownames = TRUE,
show_colnames = !TRUE,
cellwidth = cellwidth, cellheight = cellheight,
fontsize = 10, fontsize_row = 11, fontsize_col = 8,
annotation_legend = TRUE, drop_levels = TRUE)
}
## ----doheatmaps------------------------------------------------------------
clscd1 <- matClust(scd1, rowweights = weights1,
colgroups = lpdAll$`Disease Group`)
clscd2 <- matClust(scd2, rowweights = weights2,
colgroups = lpdAll$`Disease Group`, reorderrows = TRUE)
## ----gapsrow---------------------------------------------------------------
setGapPositions <- function(x, gapAt) {
rg <- if (missing(gapAt)) c(0) else {
s <- strsplit(gapAt, split = "--")
stopifnot(all(listLen(s) == 2L))
s <- strsplit(unlist(s), split = ":")
spname <- drugs[safeMatch(sapply(s, `[`, 1), drugs$name), "id"]
spconc <- as.numeric(sapply(s, `[`, 2))
spi <- mapply(function(d, cc) {
i <- which(cc == conctab[d, ])
stopifnot(length(i) == 1)
i
}, spname, spconc)
spdrug <- paste(spname, spi, sep = "_")
mt <- safeMatch(spdrug, rownames(x))
igp <- seq(1, length(mt), by = 2L)
stopifnot(all( mt[igp] - mt[igp + 1] == 1))
#stopifnot(all( mt[igp] - mt[igp + 1] == 1))
c(mt[igp + 1], 0)
}
attr(x, "rowgap") <- rg
x
}
clscd1 %<>% setGapPositions(gapAt = c(
"PF 477736:10--idelalisib:10",
"spebrutinib:2.5--PF 477736:2.5",
"PRT062607 HCl:10--selumetinib:2.5",
"selumetinib:10--MK-2206:2.5",
"MK-2206:0.156--tipifarnib:10",
"AT13387:0.039--encorafenib:10",
"encorafenib:2.5--SD07:1.111",
"doxorubicine:0.016--encorafenib:0.625",
"encorafenib:0.156--rotenone:2.5",
"SCH 900776:0.625--everolimus:0.625",
"everolimus:10--afatinib:1.667",
"arsenic trioxide:1--thapsigargin:5",
"thapsigargin:0.313--fludarabine:0.156"
))
clscd2 %<>% setGapPositions(gapAt = c(
"AT13387:0.039--everolimus:0.156",
"everolimus:0.625--nutlin-3:10",
"fludarabine:10--thapsigargin:0.078",
"thapsigargin:0.313--encorafenib:0.625",
"encorafenib:0.156--ruxolitinib:0.156"
))
## ----Fig3AheatmapV1, dev = c("png", "pdf"), fig.path = plotDir, fig.height = 8.0 + nrow(clscd1) * 0.12, fig.width = 31, fig.wide = TRUE----
#FIG# S8
rownames(clscd1) <- with(fData(lpdAll)[ rownames(clscd1),, drop = FALSE],
paste0(drugs[id, "name"], " ", conctab[cbind(id, paste0("c", subtype))], "uM"))
rownames(clscd1)
theatmap(clscd1)
## ----Fig3AheatmapV2, dev = c("png", "pdf"), fig.path = plotDir, fig.height = 8.0 + nrow(clscd2) * 0.12, fig.width = 31, fig.wide = TRUE----
#FIG# 3A
rownames(clscd2) <- with(fData(lpdAll)[ rownames(clscd2),, drop = FALSE],
paste0(drugs[id, "name"], " ", conctab[cbind(id, paste0("c", subtype))], "uM"))
rownames(clscd2)
theatmap(clscd2)
## ---- include=!exists(".standalone"), eval=!exists(".standalone")----------
# devtools::session_info()
## ---- message=FALSE, warning=FALSE, include=FALSE, eval=exists(".standalone")----
rm(list=ls())
## ---- message=FALSE, include=!exists(".standalone"), eval=!exists(".standalone")----
# library("BloodCancerMultiOmics2017")
# library("Biobase")
# library("RColorBrewer")
# library("colorspace") # hex2RGB
# library("ggtern") # ggtern
# library("reshape2") # melt
# library("limma")
# library("pheatmap")
# library("beeswarm")
# library("dplyr")
# library("ggplot2")
# library("tibble") # as_tibble
# library("survival")
# library("SummarizedExperiment")
# library("DESeq2")
## ----echo=FALSE------------------------------------------------------------
plotDir = ifelse(exists(".standalone"), "", "part10/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
## --------------------------------------------------------------------------
options(stringsAsFactors=TRUE)
## --------------------------------------------------------------------------
data("conctab", "lpdAll", "drugs", "patmeta")
## --------------------------------------------------------------------------
lpdCLL <- lpdAll[, lpdAll$Diagnosis=="CLL"]
## --------------------------------------------------------------------------
makeTransparent = function(..., alpha=0.18) {
if(alpha<0 | alpha>1) stop("alpha must be between 0 and 1")
alpha = floor(255*alpha)
newColor = col2rgb(col=unlist(list(...)), alpha=FALSE)
.makeTransparent = function(col, alpha) {
rgb(red=col[1], green=col[2], blue=col[3], alpha=alpha, maxColorValue=255)
}
newColor = apply(newColor, 2, .makeTransparent, alpha=alpha)
return(newColor)
}
giveColors = function(idx, alpha=1) {
bp = brewer.pal(12, "Paired")
makeTransparent(
sequential_hcl(12, h = coords(as(hex2RGB(bp[idx]), "polarLUV"))[1, "H"])[1],
alpha=alpha)
}
## ----calculate-input-------------------------------------------------------
# calculate (x+c)/(s+3c), (y+c)/(s+3c), (z+c)/(s+3c)
prepareTernaryData = function(lpd, targets, invDrugs) {
# calculate values for ternary
df = sapply(targets, function(tg) {
dr = paste(invDrugs[tg], c(4,5), sep="_")
tmp = 1-Biobase::exprs(lpd)[dr,]
tmp = colMeans(tmp)
pmax(tmp, 0)
})
df = data.frame(df, sum=rowSums(df), max=rowMax(df))
tern = apply(df[,targets], 2, function(x) {
(x+0.005) / (df$sum+3*0.005)
})
colnames(tern) = paste0("tern", 1:3)
# add IGHV status
cbind(df, tern, IGHV=patmeta[rownames(df),"IGHV"],
treatNaive=patmeta[rownames(df),"IC50beforeTreatment"])
}
## ----plottingFunction------------------------------------------------------
makeTernaryPlot = function(td=ternData, targets, invDrugs) {
drn = setNames(drugs[invDrugs[targets],"name"], nm=targets)
plot = ggtern(data=td, aes(x=tern1, y=tern2, z=tern3)) +
#countours
stat_density_tern(geom='polygon', aes(fill=..level..),
position = "identity", contour=TRUE, n=400,
weight = 1, base = 'identity', expand = c(1.5, 1.5)) +
scale_fill_gradient(low='lightblue', high='red', guide = FALSE) +
#points
geom_mask() +
geom_point(size=35*td[,"max"],
fill=ifelse(td[,"treatNaive"],"green","yellow"),
color="black", shape=21) +
#themes
theme_rgbw( ) +
theme_custom(
col.T=giveColors(2),
col.L=giveColors(10),
col.R=giveColors(4),
tern.plot.background="white", base_size = 18 ) +
labs( x = targets[1], xarrow = drn[targets[1]],
y = targets[2], yarrow = drn[targets[2]],
z = targets[3], zarrow = drn[targets[3]] ) +
theme_showarrows() + theme_clockwise() +
# lines
geom_Tline(Tintercept=.5, colour=giveColors(2)) +
geom_Lline(Lintercept=.5, colour=giveColors(10)) +
geom_Rline(Rintercept=.5, colour=giveColors(4))
plot
}
## --------------------------------------------------------------------------
# RUN TERNARY
makeTernary = function(lpd, targets, ighv=NA) {
# list of investigated drugs and their targets
invDrugs = c("PI3K"="D_003", "BTK"="D_002", "SYK"="D_166",
"MTOR"="D_063", "MEK"="D_012")
ternData = prepareTernaryData(lpd, targets, invDrugs)
if(!is.na(ighv)) ternData = ternData[which(ternData$IGHV==ighv),]
print(table(ternData$treatNaive))
ternPlot = makeTernaryPlot(ternData, targets, invDrugs)
ternPlot
}
## ----BCR-tern-CLL, fig.path=plotDir, dev=c('png','pdf'), fig.width = 7, fig.height = 7----
#FIG# 3B
makeTernary(lpdCLL, c("PI3K", "BTK", "SYK"), ighv=NA)
## ----BCR-tern-MCLL, fig.path=plotDir, dev=c('png','pdf'), fig.width = 7, fig.height = 7----
#FIG# 3B
makeTernary(lpdCLL, c("PI3K", "BTK", "SYK"), ighv="M")
## ----BCR-tern-UCLL, fig.path=plotDir, dev=c('png','pdf'), fig.width = 7, fig.height = 7----
#FIG# 3B
makeTernary(lpdCLL, c("PI3K", "BTK", "SYK"), ighv="U")
## ----MEK-mTOR-tern, fig.path=plotDir, dev=c('png','pdf'), fig.width = 7, fig.height = 7----
#FIG# 3BC
makeTernary(lpdCLL, c("MTOR", "BTK", "MEK"), ighv=NA)
## ----MEK-mTOR-tern-MCLL, fig.path=plotDir, dev=c('png','pdf'), fig.width = 7, fig.height = 7----
#FIG# 3BC
makeTernary(lpdCLL, c("MTOR", "BTK", "MEK"), ighv="M")
## ----MEK-mTOR-tern-UCLL, fig.path=plotDir, dev=c('png','pdf'), fig.width = 7, fig.height = 7----
#FIG# 3BC
makeTernary(lpdCLL, c("MTOR", "BTK", "MEK"), ighv="U")
## ----PI3K-MEK-mTOR-tern-CLL, fig.path=plotDir, dev=c('png','pdf'), fig.width = 7, fig.height = 7----
#FIG# S9 left
makeTernary(lpdCLL, c("MTOR", "PI3K", "MEK"), ighv=NA)
## ----SYK-MEK-mTOR-tern-CLL, fig.path=plotDir, dev=c('png','pdf'), fig.width = 7, fig.height = 7----
#FIG# S9 right
makeTernary(lpdCLL, c("MTOR", "SYK", "MEK"), ighv=NA)
## --------------------------------------------------------------------------
data("exprTreat", "drugs")
## --------------------------------------------------------------------------
e <- exprTreat
colnames( pData(e) ) <- sub( "PatientID", "Patient", colnames( pData(e) ) )
colnames( pData(e) ) <- sub( "DrugID", "Drug", colnames( pData(e) ) )
pData(e)$Drug[ is.na(pData(e)$Drug) ] <- "none"
pData(e)$Drug <- relevel( factor( pData(e)$Drug ), "none" )
pData(e)$SampleID <- colnames(e)
colnames(e) <- paste( pData(e)$Patient, pData(e)$Drug, sep=":" )
head( pData(e) )
## --------------------------------------------------------------------------
fData(e) <- fData(e)[ , c( "ProbeID", "Entrez_Gene_ID", "Symbol",
"Cytoband", "Definition" ) ]
## --------------------------------------------------------------------------
pheatmap( cor(Biobase::exprs(e)), symm=TRUE, cluster_rows = FALSE, cluster_cols = FALSE,
color = colorRampPalette(c("gray10","lightpink"))(100) )
## --------------------------------------------------------------------------
mm <- model.matrix( ~ 0 + Patient + Drug, pData(e) )
colnames(mm) <- sub( "Patient", "", colnames(mm) )
colnames(mm) <- sub( "Drug", "", colnames(mm) )
head(mm)
## --------------------------------------------------------------------------
fit <- lmFit( e, mm )
fit <- eBayes( fit )
## --------------------------------------------------------------------------
a <- decideTests( fit, p.value = 0.1 )
t( apply( a[ , grepl( "D_...", colnames(a) ) ], 2,
function(x) table( factor(x,levels=c(-1,0,1)) ) ) )
## ----overlap---------------------------------------------------------------
a <-
sapply( levels(pData(e)$Drug)[-1], function(dr1)
sapply( levels(pData(e)$Drug)[-1], function(dr2)
100*( length( intersect(
unique( topTable( fit, coef=dr1, p.value=0.1,
number=Inf )$`Entrez_Gene_ID` ),
unique( topTable( fit, coef=dr2, p.value=0.1,
number=Inf )$`Entrez_Gene_ID` ) ) ) /
length(unique( topTable( fit, coef=dr1, p.value=0.1,
number=Inf )$`Entrez_Gene_ID`)))
)
)
rownames(a) <-drugs[ rownames(a), "name" ]
colnames(a) <-rownames(a)
a <- a[-4, -4]
a
## ----corrGenes-------------------------------------------------------------
extractGenes = function(fit, coef) {
tmp = topTable(fit, coef=coef, number=Inf )[, c("ProbeID", "logFC")]
rownames(tmp) = tmp$ProbeID
colnames(tmp)[2] = drugs[coef,"name"]
tmp[order(rownames(tmp)),2, drop=FALSE]
}
runExtractGenes = function(fit, drs) {
tmp = do.call(cbind, lapply(drs, function(dr) {
extractGenes(fit, dr)
}))
as.matrix(tmp)
}
mt = runExtractGenes(fit, drs=c("D_002","D_003","D_012","D_063"))
## --------------------------------------------------------------------------
mt <- cbind( mt, median=rowMedians(mt))
mt <- mt[order(mt[,"median"]), ]
mt <- mt[1:2000, ]
mt <- mt[,-5]
## --------------------------------------------------------------------------
(mtcr = cor(mt))
## ----plot-corrGenes, fig.path=plotDir, fig.width = 4, fig.height = 4, dev = c("png", "pdf")----
#FIG# 3D
pheatmap(mtcr, cluster_rows = FALSE, cluster_cols = FALSE,
col=colorRampPalette(c("white", "lightblue","darkblue") )(100))
## --------------------------------------------------------------------------
data("lpdAll", "drugs")
## --------------------------------------------------------------------------
lpdCLL <- lpdAll[ , lpdAll$Diagnosis== "CLL"]
## ----z_factor--------------------------------------------------------------
z_factor <- qnorm(0.05, lower.tail = FALSE)
z_factor
## ----seldrugs--------------------------------------------------------------
drugnames <- c( "ibrutinib", "everolimus", "selumetinib")
ib <- "D_002_4:5"
ev <- "D_063_4:5"
se <- "D_012_4:5"
stopifnot(identical(fData(lpdAll)[c(ib, ev, se), "name"], drugnames))
## --------------------------------------------------------------------------
df <- Biobase::exprs(lpdAll)[c(ib, ev, se), lpdAll$Diagnosis=="CLL"] %>%
t %>% as_tibble %>% `colnames<-`(drugnames)
## --------------------------------------------------------------------------
mdf <- melt(data.frame(df))
## ----scatter2, fig.width = 7, fig.height = 3.5, fig.path=plotDir-----------
grid.arrange(ncol = 2,
ggplot(df, aes(x = 1-ibrutinib, y = 1-everolimus )) + geom_point(),
ggplot(df, aes(x = 1-everolimus, y = 1-selumetinib)) + geom_point()
)
## ----mirror----------------------------------------------------------------
pmdf <- filter(mdf, value >= 1)
ssd <- mean( (pmdf$value - 1) ^ 2 ) ^ 0.5
ssd
## ----dnorm-----------------------------------------------------------------
dn <- tibble(
x = seq(min(mdf$value), max(mdf$value), length.out = 100),
y = dnorm(x, mean = 1, sd = ssd) * 2 * nrow(pmdf) / nrow(mdf)
)
## ----hist1, fig.width = 10, fig.height = 5, dev = c("png", "pdf"), fig.path=plotDir----
#FIG# S10 A
thresh <- 1 - z_factor * ssd
thresh
hh <- ggplot() +
geom_histogram(aes(x = value, y = ..density..),
binwidth = 0.025, data = mdf) +
theme_minimal() +
geom_line(aes(x = x, y = y), col = "darkblue", data = dn) +
geom_vline(col = "red", xintercept = thresh)
hh
## ----hist2, fig.width = 10, fig.height = 5, dev = c("png", "pdf"), fig.path=plotDir----
hh + facet_grid( ~ variable)
## ----dec-------------------------------------------------------------------
thresh
df <- mutate(df,
group = ifelse(ibrutinib < thresh, "BTK",
ifelse(everolimus < thresh, "mTOR",
ifelse(selumetinib < thresh, "MEK", "Non-responder")))
)
## ----scatter3, fig.width = 10, fig.height = 5, dev = c("png", "pdf"), fig.path=plotDir, warning=FALSE----
#FIG# S10 B
mycol <- c(`BTK` = "Royalblue4",
`mTOR` = "chartreuse4",
`MEK` = "mediumorchid4",
`Non-responder` = "grey61")
plots <- list(
ggplot(df, aes(x = 1-ibrutinib, y = 1-everolimus)),
ggplot(filter(df, group != "BTK"), aes(x = 1-selumetinib, y = 1-everolimus))
)
plots <- lapply(plots, function(x)
x + geom_point(aes(col = group), size = 1.5) + theme_minimal() +
coord_fixed() +
scale_color_manual(values = mycol) +
geom_hline(yintercept = 1 - thresh) +
geom_vline(xintercept = 1- thresh) +
ylim(-0.15, 0.32) + xlim(-0.15, 0.32)
)
grid.arrange(ncol = 2, grobs = plots)
## --------------------------------------------------------------------------
sel = defineResponseGroups(lpd=lpdAll)
## ----co-sens-all-----------------------------------------------------------
# colors
c1 = giveColors(2, 0.5)
c2 = giveColors(4, 0.85)
c3 = giveColors(10, 0.75)
# vectors
p <- vector(); d <- vector();
pMTOR <- vector(); pBTK <- vector(); pMEK <- vector(); pNONE <- vector()
dMTOR <- vector(); dBTK <- vector(); dMEK <- vector(); dNONE <- vector()
dMTOR_NONE <- vector(); pMTOR_NONE <- vector()
# groups
sel$grMTOR_NONE <- ifelse(sel$group=="mTOR", "mTOR", NA)
sel$grMTOR_NONE <- ifelse(sel$group=="none", "none", sel$grMTOR_NONE)
sel$grMTOR <- ifelse(sel$group=="mTOR", "mTOR", "rest")
sel$col <- ifelse(sel$group=="mTOR", c3, "grey")
sel$grBTK <- ifelse(sel$group=="BTK", "BTK", "rest")
sel$col <- ifelse(sel$group=="BTK", c1, sel$col)
sel$grMEK <- ifelse(sel$group=="MEK", "MEK", "rest")
sel$col <- ifelse(sel$group=="MEK", c2, sel$col)
sel$grNONE <- ifelse(sel$group=="none", "none", "rest")
for (i in 1: max(which(fData(lpdCLL)$type=="viab"))) {
fit <- aov(Biobase::exprs(lpdCLL)[i, rownames(sel)] ~ sel$group)
p[i] <- summary(fit)[[1]][["Pr(>F)"]][1]
calc_p = function(clmn) {
p.adjust(
t.test(Biobase::exprs(lpdCLL)[i, rownames(sel)] ~ sel[,clmn],
alternative = c("two.sided") )$p.value, "BH" )
}
calc_d = function(clmn) {
diff(
t.test(Biobase::exprs(lpdCLL)[i, rownames(sel)] ~ sel[,clmn])$estimate,
alternative = c("two.sided") )
}
pMTOR_NONE[i] <- calc_p("grMTOR_NONE")
dMTOR_NONE[i] <- calc_d("grMTOR_NONE")
pMTOR[i] <- calc_p("grMTOR")
dMTOR[i] <- calc_d("grMTOR")
pBTK[i] <- calc_p("grBTK")
dBTK[i] <- calc_d("grBTK")
pMEK[i] <- calc_p("grMEK")
dMEK[i] <- calc_d("grMEK")
pNONE[i] <- calc_p("grNONE")
dNONE[i] <- calc_d("grNONE")
# drugnames
d[i] <- rownames(lpdCLL)[i]
}
## ----heatmap_mean_4_5, fig.path=plotDir, fig.width = 6, fig.height = 8, dev = c("png", "pdf")----
#FIG# 3F
#construct data frame
ps <- data.frame(drug=d, pMTOR, pBTK, pMEK, pNONE, p )
ds <- data.frame(dMTOR, dBTK, dMEK, dNONE)
rownames(ps) <- ps[,1]; rownames(ds) <- ps[,1]
# selcet only rows for singel concentrations, set non-sig to zero
ps45 <- ps[rownames(ps)[grep(rownames(ps), pattern="_4:5")],2:5 ]
for (i in 1:4) { ps45[,i] <- ifelse(ps45[,i]<0.05, ps45[,i], 0) }
ds45 <- ds[rownames(ds)[grep(rownames(ds), pattern="_4:5")],1:4 ]
for (i in 1:4) { ds45[,i] <- ifelse(ps45[,i]<0.05, ds45[,i], 0) }
# exclude non-significant rows
selDS <- rownames(ds45)[rowSums(ps45)>0]
selPS <- rownames(ps45)[rowSums(ps45)>0]
ps45 <- ps45[selPS, ]
ds45 <- ds45[selDS, ]
groupMean = function(gr) {
rowMeans(Biobase::exprs(lpdCLL)[rownames(ps45), rownames(sel)[sel$group==gr]])
}
MBTK <- groupMean("BTK")
MMEK <- groupMean("MEK")
MmTOR <- groupMean("mTOR")
MNONE <- groupMean("none")
# create data frame, new colnames
ms <- data.frame(BTK=MBTK, MEK=MMEK, mTOR=MmTOR, NONE=MNONE)
colnames(ms) <- c("BTK", "MEK", "mTOR", "WEAK")
rownames(ms) <- drugs[substr(selPS, 1,5), "name"]
# select rows with effect sizes group vs. rest >0.05
ms <- ms[ which(rowMax(as.matrix(ds45)) > 0.05 ) , ]
# exclude some drugs
ms <- ms[-c(
which(rownames(ms) %in%
c("everolimus", "ibrutinib", "selumetinib", "bortezomib"))),]
pheatmap(ms[, c("MEK", "BTK","mTOR", "WEAK")], cluster_cols = FALSE,
cluster_rows =TRUE, clustering_method = "centroid",
scale = "row",
color=colorRampPalette(
c(rev(brewer.pal(7, "YlOrRd")), "white", "white", "white",
brewer.pal(7, "Blues")))(101)
)
## ----sel-beeswarm, fig.path=plotDir, fig.width = 7, fig.height = 5.5, dev = c("png", "pdf")----
#FIG# 3G
# drug label
giveDrugLabel3 <- function(drid) {
vapply(strsplit(drid, "_"), function(x) {
k <- paste(x[1:2], collapse="_")
paste0(drugs[k, "name"])
}, character(1))
}
groups = sel[,"group", drop=FALSE]
groups[which(groups=="none"), "group"] = "WEAK"
# beeswarm function
beeDrug <- function(xDrug) {
par(bty="l", cex.axis=1.5)
beeswarm(
Biobase::exprs(lpdCLL)[xDrug, rownames(sel)] ~ groups$group,
axes=FALSE, cex.lab=1.5, ylab="Viability", xlab="", pch = 16,
pwcol=sel$col, cex=1,
ylim=c(min( Biobase::exprs(lpdCLL)[xDrug, rownames(sel)] ) - 0.05, 1.2) )
boxplot(Biobase::exprs(lpdCLL)[xDrug, rownames(sel)] ~ groups$group, add = T,
col="#0000ff22", cex.lab=2, outline = FALSE)
mtext(side=3, outer=F, line=0,
paste0(giveDrugLabel3(xDrug) ), cex=2)
}
beeDrug("D_001_4:5")
beeDrug("D_081_4:5")
beeDrug("D_013_4:5")
beeDrug("D_003_4:5")
beeDrug("D_020_4:5")
beeDrug("D_165_3")
## ----surv, fig.path=plotDir, fig.width = 8, fig.height = 8, dev = c("png", "pdf")----
#FIG# S11 A
patmeta[, "group"] <- sel[rownames(patmeta), "group"]
c1n <- giveColors(2)
c2n <- giveColors(4)
c3n <- giveColors(10)
c4n <- "lightgrey"
survplot(Surv(patmeta[ , "T5"],
patmeta[ , "treatedAfter"] == TRUE) ~ patmeta$group ,
lwd=3, cex.axis = 1.2, cex.lab=1.5, col=c(c1n, c2n, c3n, c4n),
data = patmeta,
legend.pos = 'bottomleft', stitle = 'Drug response',
xlab = 'Time (years)', ylab = 'Patients without treatment (%)',
)
## --------------------------------------------------------------------------
data(lpdAll)
## ----selectCLL-------------------------------------------------------------
lpdCLL <- lpdAll[ , lpdAll$Diagnosis== "CLL"]
## --------------------------------------------------------------------------
sel = defineResponseGroups(lpd=lpdAll)
## --------------------------------------------------------------------------
genes <- data.frame(
t(Biobase::exprs(lpdCLL)[fData(lpdCLL)$type %in% c("gen", "IGHV"), rownames(sel)]),
group = factor(sel$group)
)
genes <- genes[!(is.na(rownames(genes))), ]
colnames(genes) %<>%
sub("del13q14_any", "del13q14", .) %>%
sub("IGHV.Uppsala.U.M", "IGHV", .)
Nmut = rowSums(genes[, colnames(genes) != "group"], na.rm = TRUE)
mf <- sapply(c("BTK", "MEK", "mTOR", "none"), function(i)
mean(Nmut[genes$group==i])
)
## ----bp, fig.width = 4, fig.height = 3.5-----------------------------------
barplot(mf, ylab="Total number of mutations/CNVs per patient", col="darkgreen")
## ----mutationTests---------------------------------------------------------
mutsUse <- setdiff( colnames(genes), "group" )
mutsUse <- mutsUse[ colSums(genes[, mutsUse], na.rm = TRUE) >= 4 ]
mutationTests <- lapply(mutsUse, function(m) {
tibble(
mutation = m,
p = fisher.test(genes[, m], genes$group)$p.value)
}) %>% bind_rows %>% mutate(pBH = p.adjust(p, "BH")) %>% arrange(p)
mutationTests <- mutationTests %>% filter(pBH < 0.1)
## ----numMutationTests, fig.width = 3, fig.height = 2-----------------------
nrow(mutationTests)
## ----groupTests------------------------------------------------------------
groupTests <- lapply(unique(genes$group), function(g) {
tibble(
group = g,
p = fisher.test(
colSums(genes[genes$group == g, mutsUse], na.rm=TRUE),
colSums(genes[genes$group != g, mutsUse], na.rm=TRUE),
simulate.p.value = TRUE, B = 10000)$p.value)
}) %>% bind_rows %>% arrange(p)
groupTests
## ----ft--------------------------------------------------------------------
fisher.genes <- function(g, ref) {
stopifnot(length(ref) == 1)
ggg <- ifelse(g$group == ref, ref, "other")
idx <- which(colnames(g) != "group")
lapply(idx, function(i)
if (sum(g[, i], na.rm = TRUE) > 2) {
ft <- fisher.test(ggg, g[, i])
tibble(
p = ft$p.value,
es = ft$estimate,
prop = sum((ggg == ref) & !is.na(g[, i]), na.rm = TRUE),
mut1 = sum((ggg == ref) & (g[, i] != 0), na.rm = TRUE),
gene = colnames(g)[i])
} else {
tibble(p = 1, es = 1, prop = 0, mut1 = 0, gene = colnames(g)[i])
}
) %>% bind_rows
}
## --------------------------------------------------------------------------
pMTOR <- fisher.genes(genes, ref="mTOR")
pBTK <- fisher.genes(genes, ref="BTK")
pMEK <- fisher.genes(genes, ref="MEK")
pNONE <- fisher.genes(genes, ref="none")
p <- cbind(pBTK$p, pMEK$p, pMTOR$p, pNONE$p)
es <- cbind(pBTK$es, pMEK$es, pMTOR$es, pNONE$es)
prop <- cbind(pBTK$prop, pMEK$prop, pMTOR$prop, pNONE$prop)
mut1 <- cbind(pBTK$mut1, pMEK$mut1, pMTOR$mut1, pNONE$mut1)
## ----prepmat---------------------------------------------------------------
p <- ifelse(p < 0.05, 1, 0)
p <- ifelse(es <= 1, p, -p)
rownames(p) <- pMTOR$gene
colnames(p) <- c("BTK", "MEK", "mTOR", "NONE")
pM <- p[rowSums(abs(p))!=0, ]
propM <- prop[rowSums(abs(p))!=0, ]
## --------------------------------------------------------------------------
N <- cbind( paste0(mut1[,1],"/",prop[,1] ),
paste0(mut1[,2],"/",prop[,2] ),
paste0(mut1[,3],"/",prop[,3] ),
paste0(mut1[,4],"/",prop[,4] )
)
rownames(N) <- rownames(p)
## ----heatmap2_1, fig.path=plotDir, fig.width = 7, fig.height = 7, dev = c("png", "pdf")----
#FIG# S11 B
mutationTests <-
mutationTests[which(!(mutationTests$mutation %in%
c("del13q14_bi", "del13q14_mono"))), ]
pMA <- p[mutationTests$mutation, ]
pMA
pheatmap(pMA, cluster_cols = FALSE,
cluster_rows = FALSE, legend=TRUE, annotation_legend = FALSE,
color = c("red", "white", "lightblue"),
display_numbers = N[ rownames(pMA), ]
)
## --------------------------------------------------------------------------
data("dds")
sel = defineResponseGroups(lpd=lpdAll)
sel$group = gsub("none","weak", sel$group)
## --------------------------------------------------------------------------
# select patients with CLL in the main screen data
colnames(dds) <- colData(dds)$PatID
pat <- intersect(colnames(lpdCLL), colnames(dds))
dds_CLL <- dds[,which(colData(dds)$PatID %in% pat)]
# add group label
colData(dds_CLL)$group <- factor(sel[colnames(dds_CLL), "group"])
colData(dds_CLL)$IGHV = factor(patmeta[colnames(dds_CLL),"IGHV"])
## --------------------------------------------------------------------------
vsd <- varianceStabilizingTransformation( assay(dds_CLL) )
colnames(vsd) = colData(dds_CLL)$PatID
rowVariance <- setNames(rowVars(vsd), nm=rownames(vsd))
sortedVar <- sort(rowVariance, decreasing=TRUE)
mostVariedGenes <- sortedVar[1:10000]
dds_CLL <- dds_CLL[names(mostVariedGenes), ]
## --------------------------------------------------------------------------
cb <- combn(unique(colData(dds_CLL)$group), 2)
gg <- list(); ggM <- list(); ggU <- list()
DE <- function(ighv) {
for (i in 1:ncol(cb)) {
dds_sel <- dds_CLL[,which(colData(dds_CLL)$IGHV %in% ighv)]
dds_sel <- dds_sel[,which(colData(dds_sel)$group %in% cb[,i])]
design(dds_sel) = ~group
dds_sel$group <- droplevels(dds_sel$group)
dds_sel$group <- relevel(dds_sel$group, ref=as.character(cb[2,i]) )
dds_sel <- DESeq(dds_sel)
res <- results(dds_sel)
gg[[i]] <- res[order(res$padj, decreasing = FALSE), ]
names(gg)[i] <- paste0(cb[1,i], "_", cb[2,i])
}
return(gg)
}
## ----eval=FALSE------------------------------------------------------------
# ggM <- DE(ighv="M")
# ggU <- DE(ighv="U")
# gg <- DE(ighv=c("M", "U"))
## ----echo=FALSE, eval=FALSE------------------------------------------------
# save(ggM, ggU, gg, file=paste0(plotDir,"gexGroups.RData"))
## ----echo=TRUE, eval=TRUE--------------------------------------------------
load(system.file("extdata","gexGroups.RData", package="BloodCancerMultiOmics2017"))
## ---- eval=FALSE-----------------------------------------------------------
# library("biomaRt")
# # extract all ensembl ids
# allGenes = unique(unlist(lapply(gg, function(x) rownames(x))))
# # get gene ids for ensembl ids
# genSymbols = getBM(filters="ensembl_gene_id",
# attributes=c("ensembl_gene_id", "hgnc_symbol"),
# values=allGenes, mart=mart)
# # select first id if more than one is present
# genSymbols = genSymbols[!duplicated(genSymbols[,"ensembl_gene_id"]),]
# # set rownames to ens id
# rownames(genSymbols) = genSymbols[,"ensembl_gene_id"]
## ----echo=TRUE, eval=TRUE--------------------------------------------------
load(system.file("extdata","genSymbols.RData", package="BloodCancerMultiOmics2017"))
## ----IL10, fig.width=6.5, fig.height=6.5, dev = c("png", "pdf"), fig.path=plotDir----
#FIG# S14
gen="ENSG00000136634" #IL10
drug <- "D_063_4:5"
patsel <- intersect( rownames(sel)[sel$group %in% c("mTOR")], colnames(vsd) )
c <- cor.test( Biobase::exprs(lpdCLL)[drug, patsel], vsd[gen, patsel] )
# get hgnc_symbol for gen
# mart = useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
# genSym = getBM(filters="ensembl_gene_id", attributes="hgnc_symbol",
# values=gen, mart=mart)
genSym = genSymbols[gen, "hgnc_symbol"]
plot(vsd[gen, patsel], Biobase::exprs(lpdCLL)[drug, patsel],
xlab=paste0(genSym, " expression"),
ylab="Viability (everolimus)", pch=19, ylim=c(0.70, 0.92), col="purple",
main = paste0("mTOR-group", "\n cor = ", round(c$estimate, 3),
", p = ", signif(c$p.value,2 )),
cex=1.2)
abline(lm(Biobase::exprs(lpdCLL)[drug, patsel] ~ vsd[gen, patsel]))
## --------------------------------------------------------------------------
c1 = giveColors(2, 0.4)
c2 = giveColors(4, 0.7)
c3 = giveColors(10, 0.6)
## --------------------------------------------------------------------------
sigEx <- function(real) {
ggsig = lapply(real, function(x) {
x = data.frame(x)
x = x[which(!(is.na(x$padj))),]
x = x[x$padj<0.1,]
x = x[order(x$padj, decreasing = TRUE),]
x = data.frame(x[ ,c("padj","log2FoldChange")], ensg=rownames(x) )
x$hgnc1 = genSymbols[rownames(x), "hgnc_symbol"]
x$hgnc2 = ifelse(x$hgnc1=="" | x$hgnc1=="T" | is.na(x$hgnc1),
as.character(x$ensg), x$hgnc1)
x[-(grep(x[,"hgnc2"], pattern="IG")),]
})
return(ggsig)
}
## --------------------------------------------------------------------------
barplot1 <- function(real, tit) {
# process real diff genes
sigExPlus = sigEx(real)
ng <- lapply(sigExPlus, function(x){ cbind(
up=nrow(x[x$log2FoldChange>0, ]),
dn=nrow(x[x$log2FoldChange<0, ]) ) } )
ng = melt(ng)
p <- ggplot(ng, aes(reorder(L1, -value)), ylim(-500:500)) +
geom_bar(data = ng, aes(y = value, fill=Var2), stat="identity",
position=position_dodge() ) +
scale_fill_brewer(palette="Paired", direction = -1,
labels = c("up", "down")) +
ggtitle(label=tit) +
geom_hline(yintercept = 0,colour = "grey90") +
theme(
panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.text.x = element_text(size=14, angle = 60, hjust = 1),
axis.ticks.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_text(size=14, colour="black"),
axis.line.y = element_line(colour = "black",
size = 0.5, linetype = "solid"),
legend.key = element_rect(fill = "white"),
legend.background = element_rect(fill = "white"),
legend.title = element_blank(),
legend.text = element_text(size=14, colour="black"),
panel.background = element_rect(fill = "white", color="white")
)
plot(p)
}
## ----deGroups, fig.width = 10, fig.height = 7, dev = c("png", "pdf"), fig.path=plotDir----
#FIG# S11 C
barplot1(real=gg, tit="")
## --------------------------------------------------------------------------
# beeswarm funtion
beefun <- function(df, sym) {
par(bty="l", cex.axis=1.5)
beeswarm(df$x ~ df$y, axes=FALSE, cex.lab=1.5, col="grey", ylab=sym, xlab="",
pch = 16, pwcol=sel[colnames(vsd),"col"], cex=1.3)
boxplot(df$x ~ df$y, add = T, col="#0000ff22", cex.lab=1.5)
}
## ----groups----------------------------------------------------------------
sel = defineResponseGroups(lpdCLL)
sel[,1:3] = 1-sel[,1:3]
sel$IGHV = pData(lpdCLL)[rownames(sel), "IGHV Uppsala U/M"]
c1 = giveColors(2, 0.5)
c2 = giveColors(4, 0.85)
c3 = giveColors(10, 0.75)
# add colors
sel$col <- ifelse(sel$group=="mTOR", c3, "grey")
sel$col <- ifelse(sel$group=="BTK", c1, sel$col)
sel$col <- ifelse(sel$group=="MEK", c2, sel$col)
## --------------------------------------------------------------------------
cytokines <- c("CXCL2","TGFB1","CCL2","IL2","IL12B","IL4","IL6","IL10","CXCL8",
"TNF")
cyEN = sapply(cytokines, function(i) {
genSymbols[which(genSymbols$hgnc_symbol==i)[1],"ensembl_gene_id"]
})
makeEmpty = function() {
data.frame(matrix(ncol=3, nrow=length(cyEN),
dimnames=list(names(cyEN), c("BTK", "MEK", "mTOR"))) )
}
p = makeEmpty()
ef = makeEmpty()
## ----CYTOKINES, fig.path=plotDir, fig.width=7, fig.height=5.5, dev=c("png", "pdf")----
for (i in 1:length(cyEN) ) {
geneID <- cyEN[i]
df <- data.frame(x=vsd[geneID, ], y=sel[colnames(vsd) ,"group"])
df$y <- as.factor(df$y)
beefun(df, sym=names(geneID))
df <- within(df, y <- relevel(y, ref = "none"))
fit <- lm(x ~y, data=df)
p[i,] <- summary(fit)$coefficients[ 2:4, "Pr(>|t|)"]
abtk = mean(df[df$y=="BTK", "x"], na.rm=TRUE) - mean(df[df$y=="none", "x"],
na.rm=TRUE)
amek = mean(df[df$y=="MEK", "x"], na.rm=TRUE) - mean(df[df$y=="none", "x"],
na.rm=TRUE)
amtor= mean(df[df$y=="mTOR", "x"], na.rm=TRUE) - mean(df[df$y=="none", "x"],
na.rm=TRUE)
ef[i,] <- c(as.numeric(abtk), as.numeric(amek), as.numeric(amtor))
mtext( paste( "pBTK=", summary(fit)$coefficients[ 2, "Pr(>|t|)"],
"\npMEK=", summary(fit)$coefficients[ 3, "Pr(>|t|)"],
"\npMTOR=", summary(fit)$coefficients[ 4, "Pr(>|t|)"],
side=3 ))
}
## ----CYTOKINES-summary, fig.path=plotDir, fig.width=4.5, fig.height=3.5, dev = c("png", "pdf")----
#FIG# S11 D
# log p-values
plog <- apply(p, 2, function(x){-log(x)} )
plog_m <- melt(as.matrix(plog))
ef_m <- melt(as.matrix(ef))
# introduce effect direction
plog_m$value <- ifelse(ef_m$value>0, plog_m$value, -plog_m$value)
rownames(plog_m) <- 1:nrow(plog_m)
# fdr
fdrmin = min( p.adjust(p$mTOR, "fdr") )
### plot ####
colnames(plog_m) <- c("cytokine", "group", "p")
lev = names(sort(tapply(plog_m$p, plog_m$cytokine, function(p) min(p))))
plog_m$cytokine <- factor(plog_m$cytokine, levels=lev)
ggplot(data=plog_m, mapping=aes(x=cytokine, y=p, color=group)) +
scale_colour_manual(values = c(c1, c2, c3)) +
geom_point( size=3 ) +
theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.9),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
axis.line.x=element_line(),
axis.line.y=element_line(),
legend.position="none"
) +
scale_y_continuous(name="-log(p-value)", breaks=seq(-3,7.5,2),
limits=c(-3,7.5)) +
xlab("") +
geom_hline(yintercept = 0) +
geom_hline(yintercept = -log(0.004588897), color="purple", linetype="dashed") +
geom_hline(yintercept = (-log(0.05)), color="grey", linetype="dashed")
## ---- include=!exists(".standalone"), eval=!exists(".standalone")----------
# sessionInfo()
## ---- message=FALSE, warning=FALSE, include=FALSE--------------------------
rm(list=ls())
## ---- message=FALSE, include=!exists(".standalone"), eval=!exists(".standalone")----
# library("BloodCancerMultiOmics2017")
# library("ggplot2")
# library("tidyverse")
# library("gridExtra")
## ----echo=FALSE------------------------------------------------------------
plotDir = ifelse(exists(".standalone"), "", "part16/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
## --------------------------------------------------------------------------
data("cytokineViab")
## ----cytokinesC0, fig.width=8, fig.height=6, warning=FALSE, fig.path=plotDir, dev=c("png", "pdf")----
cond <- c("IL-2", "IL-4", "IL-10", "IL-21", "LPS", "IgM")
for (i in cond){
plot = ggplot(
filter(cytokineViab, Duplicate%in%c("1"), Stimulation==i, Timepoint=="48h"),
aes(x=as.factor(Cytokine_Concentration2), y=Normalized_DMSO, colour=mtor,
group=interaction(Patient))) +
ylab("viability") + xlab("c(stimulation)") + ylim(c(0.8, 1.4)) +
geom_line() + geom_point() + ggtitle(i) + theme_bw() + guides(color="none")
assign(paste0("p",i), plot)
}
grid.arrange(`pIL-2`,`pIL-10`,`pIL-4`,`pIL-21`,pLPS, pIgM, nrow=2)
## ---- include=!exists(".standalone"), eval=!exists(".standalone")----------
# sessionInfo()
## ---- message=FALSE, warning=FALSE, include=FALSE--------------------------
rm(list=ls())
## ---- message=FALSE, include=!exists(".standalone"), eval=!exists(".standalone")----
# library("BloodCancerMultiOmics2017")
# library("DESeq2")
# library("reshape2")
# library("dplyr")
# library("tibble")
# library("Biobase")
# library("SummarizedExperiment")
# library("genefilter")
# library("piano") # loadGSC
# library("ggplot2")
# library("gtable")
# library("grid")
## ----echo=FALSE------------------------------------------------------------
plotDir = ifelse(exists(".standalone"), "", "part07/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
## --------------------------------------------------------------------------
options(stringsAsFactors=FALSE)
## --------------------------------------------------------------------------
data(list=c("dds", "lpdAll"))
gmts = list(H=system.file("extdata","h.all.v5.1.symbols.gmt",
package="BloodCancerMultiOmics2017"),
C6=system.file("extdata","c6.all.v5.1.symbols.gmt",
package="BloodCancerMultiOmics2017"))
## --------------------------------------------------------------------------
patGroup = defineResponseGroups(lpd=lpdAll)
## ----subset----------------------------------------------------------------
lpdCLL <- lpdAll[fData(lpdAll)$type=="viab",
pData(lpdAll)$Diagnosis %in% c("CLL")]
ddsCLL <- dds[,colData(dds)$PatID %in% colnames(lpdCLL)]
## ----pat_group-------------------------------------------------------------
ddsCLL <- ddsCLL[,colData(ddsCLL)$PatID %in% rownames(patGroup)]
#remove genes without gene symbol annotations
ddsCLL <- ddsCLL[!is.na(rowData(ddsCLL)$symbol),]
ddsCLL <- ddsCLL[rowData(ddsCLL)$symbol != "",]
#add drug sensitivity annotations to coldata
colData(ddsCLL)$group <- factor(patGroup[colData(ddsCLL)$PatID, "group"])
## ----remove_low_counts-----------------------------------------------------
#only keep genes that have counts higher than 10 in any sample
keep <- apply(counts(ddsCLL), 1, function(x) any(x >= 10))
ddsCLL <- ddsCLL[keep,]
dim(ddsCLL)
## ----filtering-------------------------------------------------------------
ddsCLL <- estimateSizeFactors(ddsCLL)
sds <- rowSds(counts(ddsCLL, normalized = TRUE))
sh <- shorth(sds)
ddsCLL <- ddsCLL[sds >= sh,]
## ---- cache=TRUE-----------------------------------------------------------
ddsCLL.norm <- varianceStabilizingTransformation(ddsCLL)
## ----DE, cache=TRUE--------------------------------------------------------
DEres <- list()
design(ddsCLL) <- ~ group
rnaRaw <- DESeq(ddsCLL, betaPrior = FALSE)
#extract results for different comparisons
# responders versus weak-responders
DEres[["BTKnone"]] <- results(rnaRaw, contrast = c("group","BTK","none"))
DEres[["MEKnone"]] <- results(rnaRaw, contrast = c("group","MEK","none"))
DEres[["mTORnone"]] <- results(rnaRaw, contrast = c("group","mTOR","none"))
## --------------------------------------------------------------------------
pCut = 0.05
## ----gsea_function---------------------------------------------------------
runGSEA <- function(inputTab, gmtFile, GSAmethod="gsea", nPerm=1000){
inGMT <- loadGSC(gmtFile,type="gmt")
#re-rank by score
rankTab <- inputTab[order(inputTab[,1],decreasing = TRUE),,drop=FALSE]
if (GSAmethod == "gsea"){
#readin geneset database
#GSEA analysis
res <- runGSA(geneLevelStats = rankTab,
geneSetStat = GSAmethod,
adjMethod = "fdr", gsc=inGMT,
signifMethod = 'geneSampling', nPerm = nPerm)
GSAsummaryTable(res)
} else if (GSAmethod == "page"){
res <- runGSA(geneLevelStats = rankTab,
geneSetStat = GSAmethod,
adjMethod = "fdr", gsc=inGMT,
signifMethod = 'nullDist')
GSAsummaryTable(res)
}
}
## --------------------------------------------------------------------------
runGSE = function(gmt) {
Res <- list()
for (i in names(DEres)) {
dataTab <- data.frame(DEres[[i]])
dataTab$ID <- rownames(dataTab)
#filter using pvalues
dataTab <- filter(dataTab, pvalue <= pCut) %>%
arrange(pvalue) %>%
mutate(Symbol = rowData(ddsCLL[ID,])$symbol)
dataTab <- dataTab[!duplicated(dataTab$Symbol),]
statTab <- data.frame(row.names = dataTab$Symbol, stat = dataTab$stat)
resTab <- runGSEA(inputTab=statTab, gmtFile=gmt, GSAmethod="page")
Res[[i]] <- arrange(resTab,desc(`Stat (dist.dir)`))
}
Res
}
## ----get_genes_function----------------------------------------------------
getGenes <- function(inputTab, gmtFile){
geneList <- loadGSC(gmtFile,type="gmt")$gsc
enrichedUp <- lapply(geneList, function(x)
intersect(rownames(inputTab[inputTab[,1] >0,,drop=FALSE]),x))
enrichedDown <- lapply(geneList, function(x)
intersect(rownames(inputTab[inputTab[,1] <0,,drop=FALSE]),x))
return(list(up=enrichedUp, down=enrichedDown))
}
## ----intersection_heatmap--------------------------------------------------
plotSetHeatmap <-
function(geneTab, enrichTab, topN, gmtFile, tittle="",
asterixList = NULL, anno=FALSE) {
if (nrow(enrichTab) < topN) topN <- nrow(enrichTab)
enrichTab <- enrichTab[seq(1,topN),]
geneList <- getGenes(geneTab,gmtFile)
geneList$up <- geneList$up[enrichTab[,1]]
geneList$down <- geneList$down[enrichTab[,1]]
#form a table
allGenes <- unique(c(unlist(geneList$up),unlist(geneList$down)))
allSets <- unique(c(names(geneList$up),names(geneList$down)))
plotTable <- matrix(data=NA,ncol = length(allSets),
nrow = length(allGenes),
dimnames = list(allGenes,allSets))
for (setName in names(geneList$up)) {
plotTable[geneList$up[[setName]],setName] <- 1
}
for (setName in names(geneList$down)) {
plotTable[geneList$down[[setName]],setName] <- -1
}
if(is.null(asterixList)) {
#if no correlation table specified, order by the number of
# significant gene
geneOrder <- rev(
rownames(plotTable[order(rowSums(plotTable, na.rm = TRUE),
decreasing =FALSE),]))
} else {
#otherwise, order by the p value of correlation
asterixList <- arrange(asterixList, p)
geneOrder <- filter(
asterixList, symbol %in% rownames(plotTable))$symbol
geneOrder <- c(
geneOrder, rownames(plotTable)[! rownames(plotTable) %in% geneOrder])
}
plotTable <- melt(plotTable)
colnames(plotTable) <- c("gene","set","value")
plotTable$gene <- as.character(plotTable$gene)
if(!is.null(asterixList)) {
#add + if gene is positivily correlated with sensitivity, else add "-"
plotTable$ifSig <- asterixList[
match(plotTable$gene, asterixList$symbol),]$coef
plotTable <- mutate(plotTable, ifSig =
ifelse(is.na(ifSig) | is.na(value), "",
ifelse(ifSig > 0, "-", "+")))
}
plotTable$value <- replace(plotTable$value,
plotTable$value %in% c(1), "Up")
plotTable$value <- replace(plotTable$value,
plotTable$value %in% c(-1), "Down")
allSymbols <- plotTable$gene
geneSymbol <- geneOrder
if (anno) { #if add functional annotations in addition to gene names
annoTab <- tibble(symbol = rowData(ddsCLL)$symbol,
anno = sapply(rowData(ddsCLL)$description,
function(x) unlist(strsplit(x,"[[]"))[1]))
annoTab <- annoTab[!duplicated(annoTab$symbol),]
annoTab$combine <- sprintf("%s (%s)",annoTab$symbol, annoTab$anno)
plotTable$gene <- annoTab[match(plotTable$gene,annoTab$symbol),]$combine
geneOrder <- annoTab[match(geneOrder,annoTab$symbol),]$combine
geneOrder <- rev(geneOrder)
}
plotTable$gene <- factor(plotTable$gene, levels =geneOrder)
plotTable$set <- factor(plotTable$set, levels = enrichTab[,1])
g <- ggplot(plotTable, aes(x=set, y = gene)) +
geom_tile(aes(fill=value), color = "black") +
scale_fill_manual(values = c("Up"="red","Down"="blue")) +
xlab("") + ylab("") + theme_classic() +
theme(axis.text.x=element_text(size=7, angle = 60, hjust = 0),
axis.text.y=element_text(size=7),
axis.ticks = element_line(color="white"),
axis.line = element_line(color="white"),
legend.position = "none") +
scale_x_discrete(position = "top") +
scale_y_discrete(position = "right")
if(!is.null(asterixList)) {
g <- g + geom_text(aes(label = ifSig), vjust =0.40)
}
# construct the gtable
wdths = c(0.05, 0.25*length(levels(plotTable$set)), 5)
hghts = c(2.8, 0.1*length(levels(plotTable$gene)), 0.05)
gt = gtable(widths=unit(wdths, "in"), heights=unit(hghts, "in"))
## make grobs
ggr = ggplotGrob(g)
## fill in the gtable
gt = gtable_add_grob(gt, gtable_filter(ggr, "panel"), 2, 2)
gt = gtable_add_grob(gt, ggr$grobs[[5]], 1, 2) # top axis
gt = gtable_add_grob(gt, ggr$grobs[[9]], 2, 3) # right axis
return(list(list(plot=gt,
width=sum(wdths),
height=sum(hghts),
genes=geneSymbol)))
}
## --------------------------------------------------------------------------
statTab = setNames(lapply(c("mTORnone","BTKnone","MEKnone"), function(gr) {
dataTab <- data.frame(DEres[[gr]])
dataTab$ID <- rownames(dataTab)
#filter using pvalues
dataTab <- filter(dataTab, pvalue <= pCut) %>%
arrange(pvalue) %>%
mutate(Symbol = rowData(ddsCLL[ID,])$symbol) %>%
filter(log2FoldChange > 0)
dataTab <- dataTab[!duplicated(dataTab$Symbol),]
data.frame(row.names = dataTab$Symbol, stat = dataTab$stat)
}), nm=c("mTORnone","BTKnone","MEKnone"))
## ----run_GSE_hallmark, warning=FALSE, message= FALSE-----------------------
hallmarkRes = runGSE(gmt=gmts[["H"]])
## ----run_GSE_C6, warning=FALSE, message= FALSE-----------------------------
c6Res = runGSE(gmt=gmts[["C6"]])
## --------------------------------------------------------------------------
ddsCLL.mTOR <- ddsCLL.norm[,ddsCLL.norm$group %in% "mTOR"]
viabMTOR <- Biobase::exprs(lpdCLL["D_063_4:5", ddsCLL.mTOR$PatID])[1,]
stopifnot(all(ddsCLL.mTOR$PatID == colnames(viabMTOR)))
## --------------------------------------------------------------------------
#only keep genes that have counts higher than 10 in any sample
keep <- apply(assay(ddsCLL.mTOR), 1, function(x) any(x >= 10))
ddsCLL.mTOR <- ddsCLL.mTOR[keep,]
dim(ddsCLL.mTOR)
## --------------------------------------------------------------------------
corResult <- tibble(ID = rownames(ddsCLL.mTOR),
symbol = rowData(ddsCLL.mTOR)$symbol,
coef = NA,
p = NA)
for (i in seq(nrow(corResult))) {
res <- cor.test(viabMTOR, assay(ddsCLL.mTOR[i,])[1,], method = "pearson")
corResult[i,c(3,4)] <- c(res$estimate[[1]], res$p.value)
}
corResult <- arrange(corResult, p) %>% mutate(p.adj = p.adjust(p, method="BH"))
## --------------------------------------------------------------------------
pCut = 0.05
corResult.sig <- filter(corResult, p <= pCut)
c6Plot <- plotSetHeatmap(geneTab=statTab[["mTORnone"]],
enrichTab=c6Res[["mTORnone"]],
topN=5, gmtFile=gmts[["C6"]],
#add asterix in front of the overlapped genes
asterixList = corResult.sig,
anno=TRUE, i)
## ----fig_mTOR_C6_asterix, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.height=c6Plot[[1]][["height"]], fig.width=c6Plot[[1]][["width"]]----
#FIG# S13 B left
grid.draw(c6Plot[[1]]$plot)
## --------------------------------------------------------------------------
hallmarkPlot <- plotSetHeatmap(geneTab=statTab[["mTORnone"]],
enrichTab=hallmarkRes[["mTORnone"]],
topN=5, gmtFile=gmts[["H"]],
asterixList = corResult.sig,
anno=TRUE, i)
## ----fig_mTOR_H_asterix, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.height=hallmarkPlot[[1]][["height"]], fig.width=hallmarkPlot[[1]][["width"]]----
#FIG# S13 B right
grid.draw(hallmarkPlot[[1]]$plot)
## --------------------------------------------------------------------------
ddsCLL.BTK <- ddsCLL.norm[,ddsCLL.norm$group %in% "BTK"]
viabBTK <- Biobase::exprs(lpdCLL["D_002_4:5", ddsCLL.BTK$PatID])[1,]
stopifnot(all(ddsCLL.BTK$PatID == colnames(viabBTK)))
## --------------------------------------------------------------------------
#only keep genes that have counts higher than 10 in any sample
keep <- apply(assay(ddsCLL.BTK), 1, function(x) any(x >= 10))
ddsCLL.BTK <- ddsCLL.BTK[keep,]
dim(ddsCLL.BTK)
## --------------------------------------------------------------------------
corResult <- tibble(ID = rownames(ddsCLL.BTK),
symbol = rowData(ddsCLL.BTK)$symbol,
coef = NA,
p = NA)
for (i in seq(nrow(corResult))) {
res <- cor.test(viabBTK, assay(ddsCLL.BTK[i,])[1,], method = "pearson")
corResult[i,c(3,4)] <- c(res$estimate[[1]], res$p.value)
}
corResult <- arrange(corResult, p) %>% mutate(p.adj = p.adjust(p, method="BH"))
## --------------------------------------------------------------------------
pCut = 0.05
corResult.sig <- filter(corResult, p <= pCut)
c6Plot <- plotSetHeatmap(geneTab=statTab[["BTKnone"]],
enrichTab=c6Res[["BTKnone"]],
topN=5, gmtFile=gmts[["C6"]],
#add asterix in front of the overlapped genes
asterixList = corResult.sig,
anno=TRUE, i)
## ----fig_BTK_C6_asterix, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.height=c6Plot[[1]][["height"]], fig.width=c6Plot[[1]][["width"]]----
#FIG# S12 left
grid.draw(c6Plot[[1]]$plot)
## --------------------------------------------------------------------------
hallmarkPlot <- plotSetHeatmap(geneTab=statTab[["BTKnone"]],
enrichTab=hallmarkRes[["BTKnone"]],
topN=5, gmtFile=gmts[["H"]],
asterixList = corResult.sig,
anno=TRUE, i)
## ----fig_BTK_H_asterix, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.height=hallmarkPlot[[1]][["height"]], fig.width=hallmarkPlot[[1]][["width"]]----
#FIG# S12 right
grid.draw(hallmarkPlot[[1]]$plot)
## --------------------------------------------------------------------------
ddsCLL.MEK <- ddsCLL.norm[,ddsCLL.norm$group %in% "MEK"]
viabMEK <- Biobase::exprs(lpdCLL["D_012_4:5", ddsCLL.MEK$PatID])[1,]
stopifnot(all(ddsCLL.MEK$PatID == colnames(viabMEK)))
## --------------------------------------------------------------------------
#only keep genes that have counts higher than 10 in any sample
keep <- apply(assay(ddsCLL.MEK), 1, function(x) any(x >= 10))
ddsCLL.MEK <- ddsCLL.MEK[keep,]
dim(ddsCLL.MEK)
## --------------------------------------------------------------------------
corResult <- tibble(ID = rownames(ddsCLL.MEK),
symbol = rowData(ddsCLL.MEK)$symbol,
coef = NA,
p = NA)
for (i in seq(nrow(corResult))) {
res <- cor.test(viabMEK, assay(ddsCLL.MEK[i,])[1,], method = "pearson")
corResult[i,c(3,4)] <- c(res$estimate[[1]], res$p.value)
}
corResult <- arrange(corResult, p) %>% mutate(p.adj = p.adjust(p, method="BH"))
## --------------------------------------------------------------------------
pCut = 0.05
corResult.sig <- filter(corResult, p <= pCut)
c6Plot <- plotSetHeatmap(geneTab=statTab[["MEKnone"]],
enrichTab=c6Res[["MEKnone"]],
topN=5, gmtFile=gmts[["C6"]],
anno=TRUE, i)
## ----fig_MEK_C6_asterix, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.height=c6Plot[[1]][["height"]], fig.width=c6Plot[[1]][["width"]]----
#FIG# S13 A left
grid.draw(c6Plot[[1]]$plot)
## --------------------------------------------------------------------------
hallmarkPlot <- plotSetHeatmap(geneTab=statTab[["MEKnone"]],
enrichTab=hallmarkRes[["MEKnone"]],
topN=5, gmtFile=gmts[["H"]],
asterixList = corResult.sig,
anno=TRUE, i)
## ----fig_MEK_H_asterix, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.height=hallmarkPlot[[1]][["height"]], fig.width=hallmarkPlot[[1]][["width"]]----
#FIG# S13 A right
grid.draw(hallmarkPlot[[1]]$plot)
## ---- include=!exists(".standalone"), eval=!exists(".standalone")----------
# sessionInfo()
## ---- message=FALSE, warning=FALSE, include=FALSE--------------------------
rm(list=ls())
## ---- message=FALSE, include=!exists(".standalone"), eval=!exists(".standalone")----
# library("BloodCancerMultiOmics2017")
# library("Biobase")
# library("ggplot2")
# library("grid")
## ----echo=FALSE------------------------------------------------------------
plotDir = ifelse(exists(".standalone"), "", "part04/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
## --------------------------------------------------------------------------
options(stringsAsFactors=FALSE)
## --------------------------------------------------------------------------
data(list=c("drpar", "patmeta", "drugs", "mutCOM", "conctab"))
## --------------------------------------------------------------------------
testFactors = function(msrmnt, factors, test="student", batch=NA) {
# cut out the data
tmp = colnames(factors)
factors = data.frame(factors[rownames(msrmnt),], check.names=FALSE)
colnames(factors) = tmp
for(cidx in 1:ncol(factors))
factors[,cidx] = factor(factors[,cidx], levels=c(0,1))
# calculate the group size
groupSizes = do.call(rbind, lapply(factors, function(tf) {
tmp = table(tf)
data.frame(n.0=tmp["0"], n.1=tmp["1"])
}))
# remove the factors with less then 2 cases per group
factors = factors[,names(which(apply(groupSizes, 1,
function(i) all(i>2)))), drop=FALSE]
# calculate the effect
effect = do.call(rbind, lapply(colnames(factors), function(tf) {
tmp = aggregate(msrmnt ~ fac, data=data.frame(fac=factors[,tf]), mean)
rownames(tmp) = paste("mean", tmp$fac, sep=".")
tmp = t(tmp[2:ncol(tmp)])
data.frame(TestFac=tf,
DrugID=rownames(tmp),
FacDr=paste(tf, rownames(tmp), sep="."),
n.0=groupSizes[tf,"n.0"], n.1=groupSizes[tf,"n.1"],
tmp, WM=tmp[,"mean.0"]-tmp[,"mean.1"])
}))
# do the test
T = if(test=="student") {
do.call(rbind, lapply(colnames(factors), function(tf) {
tmp = do.call(rbind, lapply(colnames(msrmnt), function(dr) {
res = t.test(msrmnt[,dr] ~ factors[,tf], var.equal=TRUE)
data.frame(DrugID=dr, TestFac=tf,
pval=res$p.value, t=res$statistic,
conf1=res$conf.int[1], conf2=res$conf.int[2])
}))
tmp
}))
} else if(test=="anova") {
do.call(rbind, lapply(colnames(factors), function(tf) {
tmp = do.call(rbind, lapply(colnames(msrmnt), function(dr) {
# make sure that the order in batch is the same as in msrmnt
stopifnot(identical(rownames(msrmnt), names(batch)))
res = anova(lm(msrmnt[,dr] ~ factors[,tf]+batch))
data.frame(DrugID=dr, TestFac=tf, pval=res$`Pr(>F)`[1],
f=res$`F value`[1], meanSq1=res$`Mean Sq`[1],
meanSq2=res$`Mean Sq`[2])
}))
tmp
}))
} else {
NA
}
enhanceObject = function(obj) {
# give nice drug names
obj$Drug = giveDrugLabel(obj$DrugID, conctab, drugs)
# combine the testfac and drug id
obj$FacDr = paste(obj$TestFac, obj$DrugID, sep=".")
# select just the drug name
obj$DrugID2 = substring(obj$DrugID, 1, 5)
obj
}
list(effect=effect, test=enhanceObject(T))
}
## --------------------------------------------------------------------------
## VIABILITIES
## list of matrices; one matrix per screen/day
## each matrix contains all CLL patients
measurements=list()
### Main Screen
patM = colnames(drpar)[which(patmeta[colnames(drpar),"Diagnosis"]=="CLL")]
measurements[["main"]] =
do.call(cbind,
lapply(list("viaraw.1","viaraw.2","viaraw.3","viaraw.4","viaraw.5"),
function(viac) {
tmp = t(assayData(drpar)[[viac]][,patM])
colnames(tmp) = paste(colnames(tmp), conctab[colnames(tmp),
paste0("c",substring(viac,8,8))], sep="-")
tmp
}))
pats = sort(unique(patM))
## TESTING FACTORS
testingFactors = list()
# ighv
ighv = setNames(patmeta[pats, "IGHV"], nm=pats)
# mutations
tmp = cbind(IGHV=ifelse(ighv=="U",1,0), assayData(mutCOM)$binary[pats,])
testingFactors[["mutation"]] = tmp[,-grep("Chromothripsis", colnames(tmp))]
# BATCHES
j = which(pData(drpar)[patM, "ExpDate"] < as.Date("2014-01-01"))
k = which(pData(drpar)[patM, "ExpDate"] < as.Date("2014-08-01") &
pData(drpar)[patM, "ExpDate"] > as.Date("2014-01-01"))
l = which(pData(drpar)[patM, "ExpDate"] > as.Date("2014-08-01"))
measurements[["main"]] = measurements[["main"]][c(patM[j], patM[k], patM[l]),]
batchvec = factor(
setNames(c(rep(1, length(j)), rep(2, length(k)), rep(3, length(l))),
nm=c(patM[j], patM[k], patM[l])))
# LABELS FOR GROUPING
beelabs = t(sapply(colnames(testingFactors[["mutation"]]), function(fac) {
if(fac=="IGHV")
c(`0`="IGHV mut", `1`="IGHV unmut")
else if(grepl("[[:upper:]]",fac)) # if all letters are uppercase
c(`0`=paste(fac, "wt"),`1`=paste(fac, "mt"))
else
c(`0`="wt",`1`=fac)
}))
## --------------------------------------------------------------------------
allresultsT = testFactors(msrmnt=measurements[["main"]],
factors=testingFactors[["mutation"]],
test="student", batch=NA)
resultsT = allresultsT$test
resultsT$adj.pval = p.adjust(resultsT$pval, method="BH")
## --------------------------------------------------------------------------
allresultsA = testFactors(msrmnt=measurements[["main"]],
factors=testingFactors[["mutation"]],
test="anova", batch=batchvec)
resultsA = allresultsA$test
resultsA$adj.pval = p.adjust(resultsA$pval, method="BH")
## ----batchEffect, results='asis', echo=FALSE, fig.path=plotDir, dev=c('png','pdf'), fig.height=20, fig.width=14----
#FIG# S30
xylim = 1e-8
tmp = merge(resultsT[,c("FacDr","Drug","DrugID","DrugID2","FacDr","pval")],
resultsA[,c("FacDr","pval")], by.x="FacDr", by.y="FacDr")
tmp$DrugName = toCaps(drugs[tmp$DrugID2, "name"])
tmp$Shape = ifelse(tmp$pval.x < xylim | tmp$pval.y < xylim, "tri", "dot")
tmp$pval.cens.x = ifelse(tmp$pval.x < xylim, xylim, tmp$pval.x)
tmp$pval.cens.y = ifelse(tmp$pval.y < xylim, xylim, tmp$pval.y)
ggplot(tmp) + geom_abline(intercept=0, slope=1, colour="hotpink") +
geom_point(aes(x=pval.cens.x, y=pval.cens.y, shape=Shape), alpha=0.6) +
facet_wrap(~DrugName, ncol=7) +
scale_x_log10(labels=scientific_10, limits=c(xylim,1)) +
scale_y_log10(labels=scientific_10, limits=c(xylim,1)) +
theme_bw() + coord_fixed() + xlab("P-value from Student t-test") +
ylab("P-value from ANOVA (including batch group factor)") +
theme(axis.text.x=element_text(angle=0, hjust=1, vjust=0.5, size=rel(1.5)),
axis.title=element_text(size=18)) + guides(shape=FALSE)
## --------------------------------------------------------------------------
badrugs = c("D_008", "D_025")
measurements = lapply(measurements, function(drres) {
drres[,-grep(paste(badrugs, collapse="|"), colnames(drres))]
})
## --------------------------------------------------------------------------
allresults1 = lapply(measurements, function(measurement) {
testFactors(msrmnt=measurement, factors=testingFactors[["mutation"]],
test="student", batch=NA)
})
effects1 = lapply(allresults1, function(res) res[["effect"]])
results1 = lapply(allresults1, function(res) res[["test"]])
results1 = lapply(results1, function(res) {
res$adj.pval = p.adjust(res$pval, method="BH")
res
})
## --------------------------------------------------------------------------
measurements1 = measurements
testingFactors1 = testingFactors
beelabs1 = beelabs
## ----echo=FALSE------------------------------------------------------------
## CREATE THE PLOTS
plotmp01 = BloodCancerMultiOmics2017:::run.ggvolcGr2(results=results1, effects=effects1,
screen="main", mts="IGHV", fdr=0.1, maxX=0.75,
maxY=NA, expY=0.05, hghBox=0.15, axisMarkY=4,
breaksX=c(-0.75,-0.5,-0.25,0,0.25,0.5,0.75),
arrowLength=0.5, Xhang=0.3, minConc=3)
plotmp02 = BloodCancerMultiOmics2017:::run.ggvolcGr2(results=results1, effects=effects1,
screen="main", mts="trisomy12", fdr=0.1,
maxX=0.75, maxY=NA, expY=0.05, hghBox=0.15,
axisMarkY=4,
breaksX=c(-0.75,-0.5,-0.25,0,0.25,0.5,0.75),
arrowLength=0.5, Xhang=0.3, minConc=1)
## ----volc_IGHV, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp01$IGHV$figure$width, fig.height=plotmp01$IGHV$figure$height----
#FIG# 4B
grid.draw(plotmp01$IGHV$figure$plot)
## ----volc_IGHV_legend, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp01$IGHV$legend$width, fig.height=plotmp01$IGHV$legend$height, out.height=300, out.width=150----
#FIG# 4B legend
grid.draw(plotmp01$IGHV$legend$plot)
## ----volc_trisomy12, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp02$trisomy12$figure$width, fig.height=plotmp02$trisomy12$figure$height----
#FIG# 4C
grid.draw(plotmp02$trisomy12$figure$plot)
## ----volc_trisomy12_legend, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp02$trisomy12$legend$width, fig.height=plotmp02$trisomy12$legend$height, out.height=300, out.width=150----
#FIG# 4C legend
grid.draw(plotmp02$trisomy12$legend$plot)
## --------------------------------------------------------------------------
fac2test = lapply(measurements, function(mea) {
tf = testingFactors[["mutation"]][rownames(mea),]
names(which(apply(tf,2,function(i) {
if(length(table(i,tf[,1]))!=4)
FALSE
else
all(table(i,tf[,1])>2)
})))
})
## --------------------------------------------------------------------------
measurements2 = setNames(lapply(names(measurements), function(mea) {
ig = testingFactors[["mutation"]][rownames(measurements[[mea]]),"IGHV"]
patU = names(which(ig==1))
patM = names(which(ig==0))
list(U=measurements[[mea]][patU,], M=measurements[[mea]][patM,])
}), nm=names(measurements))
## --------------------------------------------------------------------------
allresults2 = setNames(lapply(names(measurements2), function(mea) {
list(U = testFactors(msrmnt=measurements2[[mea]]$U,
factors=testingFactors[["mutation"]][
rownames(measurements2[[mea]]$U),fac2test[[mea]]]),
M = testFactors(msrmnt=measurements2[[mea]]$M,
factors=testingFactors[["mutation"]][
rownames(measurements2[[mea]]$M),fac2test[[mea]]]))
}), nm=names(measurements2))
## --------------------------------------------------------------------------
results2 = lapply(allresults2, function(allres) {
list(U=allres[["U"]][["test"]], M=allres[["M"]][["test"]])
})
effects2 = lapply(allresults2, function(allres) {
list(U=allres[["U"]][["effect"]], M=allres[["M"]][["effect"]])
})
## --------------------------------------------------------------------------
results2 = lapply(results2, function(res) {
tmp = p.adjust(c(res$U$pval,res$M$pval), method="BH")
l = length(tmp)
res$U$adj.pval = tmp[1:(l/2)]
res$M$adj.pval = tmp[(l/2+1):l]
res
})
## --------------------------------------------------------------------------
testingFactors2 = testingFactors
beelabs2 = beelabs
## ----echo=FALSE------------------------------------------------------------
## CREATE THE PLOTS
plotmp03 = BloodCancerMultiOmics2017:::run.ggvolcGr2(results=results2$main, effects=effects2$main,
screen="U", mts="trisomy12", fdr=0.1,
maxX=0.75, maxY=7, expY=0.05, hghBox=NA,
axisMarkY=4,
breaksX=c(-0.75,-0.5,-0.25,0,0.25,0.5,0.75),
arrowLength=0.5, Xhang=0.3, minConc=1,
fixedHght=6)
plotmp04 = BloodCancerMultiOmics2017:::run.ggvolcGr2(results=results2$main, effects=effects2$main,
screen="M", mts="trisomy12", fdr=0.1,
maxX=0.75, maxY=7, expY=0.05, hghBox=NA,
axisMarkY=4,
breaksX=c(-0.75,-0.5,-0.25,0,0.25,0.5,0.75),
arrowLength=0.5, Xhang=0.3, minConc=1,
fixedHght=6)
## ----volc_trisomy12_U, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp03$trisomy12$figure$width, fig.height=plotmp03$trisomy12$figure$height----
#FIG# S21 right
grid.draw(plotmp03$trisomy12$figure$plot)
## ----volc_trisomy12_U_legend, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp03$trisomy12$legend$width, fig.height=plotmp03$trisomy12$legend$height, out.height=300, out.width=150----
#FIG# S21 right legend
grid.draw(plotmp03$trisomy12$legend$plot)
## ----volc_trisomy12_M, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp04$trisomy12$figure$width, fig.height=plotmp04$trisomy12$figure$height----
#FIG# S21 left
grid.draw(plotmp04$trisomy12$figure$plot)
## ----volc_trisomy12_M_legend, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp04$trisomy12$legend$width, fig.height=plotmp04$trisomy12$legend$height, out.height=300, out.width=150----
#FIG# S21 left legend
grid.draw(plotmp04$trisomy12$legend$plot)
## --------------------------------------------------------------------------
## VIABILITIES
### main
pats = colnames(drpar)
# make the big matrix with viabilities
measureTMP = do.call(cbind,
lapply(list("viaraw.1","viaraw.2","viaraw.3",
"viaraw.4","viaraw.5"), function(viac) {
tmp = t(assayData(drpar)[[viac]][,pats])
colnames(tmp) = paste(colnames(tmp),
conctab[colnames(tmp),
paste0("c",substring(viac,8,8))], sep="-")
tmp
}))
# select diagnosis to work on
pats4diag = tapply(pats, patmeta[pats,"Diagnosis"], function(i) i)
diags = names(which(table(patmeta[pats,"Diagnosis"])>2))
diags = diags[-which(diags=="CLL")]
# there will be two lists: one with CLL and the second with other diagnosis
# (first one is passed as argument to the createObjects function)
pats4diag2 = pats4diag[diags]
# function that creates testingFactors, measurements and beelabs
createObjects = function(pats4diag1, beefix="") {
measurements=list()
testingFactors=list()
# make the list for testing
for(m in names(pats4diag1)) {
for(n in names(pats4diag2)) {
p1 = pats4diag1[[m]]
p2 = pats4diag2[[n]]
measurements[[paste(m,n,sep=".")]] = measureTMP[c(p1, p2),]
testingFactors[[paste(m,n,sep=".")]] = setNames(c(rep(0,length(p1)),
rep(1,length(p2))),
nm=c(p1,p2))
}
}
# reformat testingFactors to the df
pats=sort(unique(c(unlist(pats4diag1),unlist(pats4diag2))))
testingFactors = as.data.frame(
do.call(cbind, lapply(testingFactors, function(tf) {
setNames(tf[pats], nm=pats)
})))
# Labels for beeswarms
beelabs = t(sapply(colnames(testingFactors), function(fac) {
tmp = unlist(strsplit(fac, ".", fixed=TRUE))
c(`0`=paste0(tmp[1], beefix),`1`=tmp[2])
}))
return(list(msrmts=measurements, tf=testingFactors, bl=beelabs))
}
# all CLL together
res = createObjects(pats4diag1=pats4diag["CLL"])
measurements3 = res$msrmts
testingFactors3 = res$tf
beelabs3 = res$bl
## ---- results='hide'-------------------------------------------------------
allresults3 = setNames(lapply(names(measurements3), function(mea) {
tmp = data.frame(testingFactors3[,mea])
colnames(tmp) = mea
rownames(tmp) = rownames(testingFactors3)
testFactors(msrmnt=measurements3[[mea]], factors=tmp)
}), nm=names(measurements3))
## --------------------------------------------------------------------------
results3 = lapply(allresults3, function(res) res[["test"]])
effects3 = lapply(allresults3, function(res) res[["effect"]])
## --------------------------------------------------------------------------
results3 = lapply(results3, function(res) {
res$adj.pval = p.adjust(res$pval, method="BH")
res
})
## ----echo=FALSE------------------------------------------------------------
tmpheat = BloodCancerMultiOmics2017:::ggheat(results3, effects3)
## ----cll.diag, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=tmpheat$figure$width, fig.height=tmpheat$figure$height----
#FIG# S7 plot
grid.draw(tmpheat$figure$plot)
## ----cll.diag.legend, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=tmpheat$legend$width, fig.height=tmpheat$legend$height----
#FIG# S7 legend
grid.draw(tmpheat$legend$plot)
## --------------------------------------------------------------------------
data(drugs, lpdAll, mutCOM, conctab)
## --------------------------------------------------------------------------
lpdCLL = lpdAll[ , lpdAll$Diagnosis %in% "CLL"]
## ----beesMutMain, fig.path=plotDir, dev=c("png", "pdf"), fig.width=8, fig.height=10, out.width=280, out.height=350----
#FIG# 4D
BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_010_2", "TP53", cs=T, y1=0, y2=1.2, custc=T)
BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_006_3", "TP53", cs=T,y1=0, y2=1.2, custc=T)
BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_063_5", "CREBBP", cs=T, y1=0, y2=1.2, custc=T)
BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_056_5", "PRPF8", cs=T, y1=0, y2=1.2, custc=T)
BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_012_5", "trisomy12", cs=F,y1=0.6, y2=1.2, custc=T)
## ----beesMutSupp, fig.path=plotDir, fig.width = 18, fig.height = 20, dev = c("png", "pdf")----
#FIG# S17
par(mfrow = c(3,4), mar=c(5,4.5,5,2))
BloodCancerMultiOmics2017:::beeF(diag="CLL", drug="D_159_3", mut="TP53", cs=T, y1=0, y2=1.2, custc=T)
BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_006_2", "del17p13", cs=T, y1=0, y2=1.2, custc=T)
BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_159_3", "del17p13", cs=T, y1=0, y2=1.2, custc=T)
BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_010_2", "del17p13", cs=T, y1=0, y2=1.2, custc=T)
BloodCancerMultiOmics2017:::beeF(diag="MCL", "D_006_2", "TP53", cs=T, y1=0, y2=1.2, custc=T)
BloodCancerMultiOmics2017:::beeF(diag="MCL", "D_010_2", "TP53", cs=T, y1=0, y2=1.2, custc=T)
BloodCancerMultiOmics2017:::beeF(diag=c("HCL", "HCL-V"), "D_012_3", "BRAF", cs=T, y1=0, y2=1.2,
custc=F)
BloodCancerMultiOmics2017:::beeF(diag=c("HCL", "HCL-V"), "D_083_4", "BRAF", cs=T, y1=0, y2=1.2,
custc=F)
BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_012_5", "KRAS", cs=T, y1=0.6, y2=1.2, custc=T)
BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_083_5", "KRAS", cs=T, y1=0.6, y2=1.45, custc=T)
BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_081_4", "UMODL1", cs=T, y1=0, y2=1.2, custc=T)
BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_001_4", "UMODL1", cs=T, y1=0, y2=1.2, custc=T)
## ----colorbar, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=2, fig.height=4, out.height=200, out.width=100----
# the below function comes from the CRAN package named monogeneaGM
colorBar <-function(colpalette=NULL, min, max=-min, nticks=11, ticks=seq(min, max, len=nticks),
tit="") {
if(is.null(colpalette)) {
scale <- (length(colpalette)-1)/(max-min)
rwb <- c("#99000D","#FB6A4A","white","#6BAED6","#084594")
colpalette<- colorRampPalette(rwb, space="Lab")(101)
}
scale <- (length(colpalette)-1)/(max-min)
plot(c(0,10), c(min,max), type="n", bty="n", xaxt="n", xlab="", yaxt="n", ylab="", main=tit)
axis(2, ticks, las=1)
for (i in 1:(length(colpalette)-1)) {
y = (i-1)/scale + min
rect(0,y,10,y+1/scale, col=colpalette[i], border=NA)
}
}
colorBar(colorRampPalette(c('coral1','blue4'))(100), min=0, max = 1,
ticks=c(0,0.5,1))
## ----bee-pretreatment, fig.path=plotDir, fig.width=15, fig.height=10, dev = c("png", "pdf")----
#FIG# S18
par(mfrow = c(2,3), mar=c(5,4.5,2,2))
BloodCancerMultiOmics2017:::beePretreatment(lpdCLL, "D_006_1:5", y1=0.2, y2=1.3, fac="TP53",
val=c(0,1), name="Fludarabine")
BloodCancerMultiOmics2017:::beePretreatment(lpdCLL, "D_006_1:5", y1=0.2, y2=1.3, fac="TP53",
val=c(0), name="p53 wt: Fludarabine")
BloodCancerMultiOmics2017:::beePretreatment(lpdCLL, "D_006_1:5", y1=0.2, y2=1.3, fac="TP53",
val=c(1), name="p53 mut: Fludarabine")
BloodCancerMultiOmics2017:::beePretreatment(lpdCLL, "D_002_4:5", y1=0.4, y2=1.3,
fac="IGHV Uppsala U/M", val=c(0,1), name="Ibrutinib")
BloodCancerMultiOmics2017:::beePretreatment(lpdCLL, "D_002_4:5", y1=0.4, y2=1.3,
fac="IGHV Uppsala U/M", val=c(0), name="U-CLL: Ibrutinib")
BloodCancerMultiOmics2017:::beePretreatment(lpdCLL, "D_002_4:5", y1=0.4, y2=1.3,
fac="IGHV Uppsala U/M", val=c(1), name="M-CLL: Ibrutinib")
## ---- include=!exists(".standalone"), eval=!exists(".standalone")----------
# sessionInfo()
## ---- message=FALSE, warning=FALSE, include=FALSE--------------------------
rm(list=ls())
## ---- message=FALSE, warning=FALSE, include=!exists(".standalone"), eval=!exists(".standalone")----
# library("BloodCancerMultiOmics2017")
# library("Biobase")
# library("dplyr")
# library("tidyr")
# library("broom")
# library("ggplot2")
# library("grid")
# library("gridExtra")
# library("reshape2")
# library("foreach")
# library("doParallel")
# library("scales")
# library("knitr")
# registerDoParallel()
## ----echo=FALSE------------------------------------------------------------
plotDir = ifelse(exists(".standalone"), "", "part05/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
## --------------------------------------------------------------------------
options(stringsAsFactors=FALSE)
## --------------------------------------------------------------------------
# get drug responsee data
get.drugresp <- function(lpd) {
drugresp = t(Biobase::exprs(lpd[fData(lpd)$type == 'viab'])) %>%
dplyr::tbl_df() %>% dplyr::select(-ends_with(":5")) %>%
dplyr::mutate(ID = colnames(lpd)) %>%
tidyr::gather(drugconc, viab, -ID) %>%
dplyr::mutate(drug = drugs[substring(drugconc, 1, 5), "name"],
conc = sub("^D_([0-9]+_)", "", drugconc)) %>%
dplyr::mutate(conc = as.integer(gsub("D_CHK_", "", conc)))
drugresp
}
# extract mutations and IGHV status
get.somatic <- function(lpd) {
somatic = t(Biobase::exprs(lpd[Biobase::fData(lpd)$type == 'gen' |
Biobase::fData(lpd)$type == 'IGHV']))
## rename IGHV Uppsala to 'IGHV' (simply)
colnames(somatic)[grep("IGHV", colnames(somatic))] = "IGHV"
## at least 3 patients should have this mutation
min.samples = which(Matrix::colSums(somatic, na.rm = T) > 2)
somatic = dplyr::tbl_df(somatic[, min.samples]) %>%
dplyr::select(-one_of("del13q14_bi", "del13q14_mono",
"Chromothripsis", "RP11-766F14.2")) %>%
dplyr::rename(del13q14 = del13q14_any) %>%
dplyr::mutate(ID = colnames(lpd)) %>%
tidyr::gather(mutation, mut.value, -ID)
somatic
}
## ----ggTheme---------------------------------------------------------------
t1<-theme(
plot.background = element_blank(),
panel.grid.major = element_line(),
panel.grid.major.x = element_line(linetype = "dotted", colour = "grey"),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(size=.4),
axis.line.x = element_line(),
axis.line.y = element_line(),
axis.text.x = element_text(angle=90, size=12,
face="bold", hjust = 1, vjust = 0.4),
axis.text.y = element_text(size = 14),
axis.ticks.x = element_line(linetype = "dotted"),
axis.ticks.length = unit(0.3,"cm"),
axis.title.x = element_text(face="bold", size=16),
axis.title.y = element_text(face="bold", size=20),
plot.title = element_text(face="bold", size=16, hjust = 0.5)
)
## theme for the legend
t.leg <- theme(legend.title = element_text(face='bold',
hjust = 1, size=11),
legend.position = c(0, 0.76),
legend.key = element_blank(),
legend.text = element_text(size=12),
legend.background = element_rect(color = "black"))
## ----colorPalette----------------------------------------------------------
colors= c("#015872","#3A9C94","#99977D","#ffbf00","#5991C7","#99cc00",
"#D5A370","#801416","#B2221C","#ff5050","#33bbff","#5c5cd6",
"#E394BB","#0066ff","#C0C0C0")
## --------------------------------------------------------------------------
get.pretreat <- function(patmeta, lpd) {
patmeta = patmeta[rownames(patmeta) %in% colnames(lpd),]
data.frame(ID=rownames(patmeta), pretreat=!patmeta$IC50beforeTreatment) %>%
mutate(pretreat = as.factor(pretreat))
}
## --------------------------------------------------------------------------
make.dr <- function(resp, features, patmeta, lpd) {
treat = get.pretreat(patmeta, lpd)
dr = full_join(resp, features) %>%
inner_join(treat)
}
## --------------------------------------------------------------------------
get.medp <- function(drugresp) {
tab = drugresp %>% group_by(drug, conc) %>%
do(data.frame(v = .$viab, ID = .$ID)) %>% spread(ID, v)
med.p = foreach(n=unique(tab$drug), .combine = cbind) %dopar% {
tb = filter(tab, drug == n) %>% ungroup() %>% dplyr::select(-(drug:conc)) %>%
as.matrix %>% `rownames<-`(1:5)
mdp = stats::medpolish(tb)
df = as.data.frame(mdp$col) + mdp$overall
colnames(df) <- n
df
}
medp.viab = dplyr::tbl_df(med.p) %>% dplyr::mutate(ID = rownames(med.p)) %>%
tidyr::gather(drug, viab, -ID)
medp.viab
}
## --------------------------------------------------------------------------
get.labels <- function(pvals) {
lev = levels(factor(pvals$mutation))
lev = gsub("^(gain)([0-9]+)([a-z][0-9]+)$", "\\1(\\2)(\\3)", lev)
lev = gsub("^(del)([0-9]+)([a-z].+)$", "\\1(\\2)(\\3)", lev)
lev = gsub("trisomy12", "trisomy 12", lev)
lev
}
## --------------------------------------------------------------------------
get.mutation.order <- function(lev) {
ord = c("trisomy 12", "TP53",
"del(11)(q22.3)", "del(13)(q14)",
"del(17)(p13)",
"gain(8)(q24)",
"BRAF", "CREBBP", "PRPF8",
"KLHL6", "NRAS", "ABI3BP", "UMODL1")
mut.order = c(match(ord, lev),
grep("Other", lev), grep("Below", lev))
mut.order
}
## --------------------------------------------------------------------------
get.drug.order <- function(pvals, drugs) {
## determine drug order by column sums of log-p values
dr.order = pvals %>%
mutate(logp = -log10(p.value)) %>%
group_by(drug) %>% summarise(logsum = sum(logp))
dr.order = inner_join(dr.order, pvals %>%
group_by(drug) %>%
summarise(n = length(unique(mutation)))) %>%
arrange(desc(n), desc(logsum))
dr.order = inner_join(dr.order, drugs %>% rename(drug = name))
dr.order = left_join(dr.order, dr.order %>%
group_by(`target_category`) ) %>%
arrange(`target_category`, drug) %>%
filter(! `target_category` %in% c("ALK", "Angiogenesis", "Other")) %>%
filter(!is.na(`target_category`))
dr.order
}
## --------------------------------------------------------------------------
make.annot <- function(g, dr.order) {
# make a color palette for drug pathways
drug.class = c("#273649", "#647184", "#B1B2C8",
"#A7755D", "#5D2E1C", "#38201C")
pathways = c("BH3","B-cell receptor","DNA damage",
"MAPK", "PI3K", "Reactive oxygen species")
names(pathways) = c("BH3", "BCR inhibitors", "DNA damage",
"MAPK", "PI3K", "ROS")
for (i in 1:6) {
prange = grep(pathways[i], dr.order$`target_category`)
path.grob <- grobTree(rectGrob(gp=gpar(fill=drug.class[i])),
textGrob(names(pathways)[i],
gp = gpar(cex =0.8, col = "white")))
g = g +
annotation_custom(path.grob,
xmin = min(prange) -0.25 - 0.1 * ifelse(i == 2, 1, 0),
xmax = max(prange) + 0.25 + 0.1 * ifelse(i == 2, 1, 0),
ymin = -0.52, ymax = -0.2)
}
g
}
## --------------------------------------------------------------------------
g_legend<-function(a.gplot){
tmp <- ggplot_gtable(ggplot_build(a.gplot))
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
legend
} ## end define
## --------------------------------------------------------------------------
data(list=c("conctab", "drugs", "lpdAll", "patmeta"))
## ----preprocesslpd---------------------------------------------------------
lpdCLL <- lpdAll[ , lpdAll$Diagnosis=="CLL"]
## extract viability data for 5 different concentrations
drugresp = get.drugresp(lpdCLL)
## ----preprocessMuts--------------------------------------------------------
## extract somatic variants
somatic = get.somatic(lpdCLL) %>%
mutate(mut.value = as.factor(mut.value))
## ----medPolish, warning=FALSE, results='hide'------------------------------
## compute median polish patient effects and recalculate p-values
medp.viab = get.medp(drugresp)
dr = make.dr(medp.viab, somatic, patmeta, lpdCLL)
## ----pvalsAndFDR-----------------------------------------------------------
pvals = dr %>% group_by(drug, mutation) %>%
do(tidy(t.test(viab ~ mut.value, data = ., var.equal = T))) %>%
dplyr::select(drug, mutation, p.value)
# compute the FDR threshold
fd.thresh = 10
padj = p.adjust(pvals$p.value, method = "BH")
fdr = max(pvals$p.value[which(padj <= fd.thresh/100)])
## ----subsetPvals-----------------------------------------------------------
# selected mutations
select.mutations = c("trisomy12", "TP53",
"del11q22.3", "del13q14",
"del17p13",
"gain8q24",
"BRAF", "CREBBP", "PRPF8",
"KLHL6", "NRAS", "ABI3BP", "UMODL1")
pvals = filter(pvals, mutation != 'IGHV')
pvals = pvals %>% ungroup() %>%
mutate(mutation = ifelse(p.value > fdr,
paste0("Below ", fd.thresh,"% FDR"), mutation)) %>%
mutate(mutation = ifelse(!(mutation %in% select.mutations) &
!(mutation == paste0("Below ", fd.thresh,"% FDR")),
"Other", mutation)) %>%
filter(drug != "bortezomib" & drug != "NSC 74859")
## ----renameMuts------------------------------------------------------------
## order of mutations
lev = get.labels(pvals)
folge = get.mutation.order(lev)
## ----pvalsForGGplot--------------------------------------------------------
drugs = drugs[,c("name", "target_category")]
# get the drug order
dr.order = get.drug.order(pvals, drugs)
## --------------------------------------------------------------------------
plot.pvalues <- function(pvals, dr.order, folge, colors, shapes) {
g = ggplot(data = filter(pvals, drug %in% dr.order$drug)) +
geom_point(aes(x = factor(drug, levels = dr.order$drug), y = -log10(p.value),
colour = factor(mutation, levels(factor(mutation))[folge]),
shape = factor(mutation, levels(factor(mutation))[folge])),
size=5, show.legend = T) +
scale_color_manual(name = "Mutations",
values = colors,
labels = lev[folge]) +
scale_shape_manual(name = "Mutations",
values = shapes,
labels = lev[folge]) + t1 +
labs(x = "", y = expression(paste(-log[10], "p")), title = "") +
scale_y_continuous(expression(italic(p)*"-value"),
breaks=seq(0,10,5),
labels=math_format(expr=10^.x)(-seq(0,10,5)))
g
}
## ----pvalsMain, fig.path=plotDir, dev=c("png", "pdf"), fig.width=14, fig.height=10----
#FIG# 4A
## plot the p-values
g = plot.pvalues(pvals, dr.order, folge,
colors, shapes = c(rep(16,13), c(1,1)))
## add FDR threshold
g = g + geom_hline(yintercept = -log10(fdr),
linetype="dashed", size=0.3)
g = g +
annotation_custom(grob = textGrob(label = paste0("FDR", fd.thresh, "%"),
hjust = 1, vjust = 1,
gp = gpar(cex = 0.5,
fontface = "bold",
fontsize = 25)),
ymin = -log10(fdr) - 0.2,
ymax = -log10(fdr) + 0.5,
xmin = -1.3, xmax = 1.5) +
theme(legend.position = "none")
# generate pathway/target annotations for certain drug classes
#g = make.annot(g, dr.order)
# legend guide
leg.guides <- guides(colour = guide_legend(ncol = 1,
byrow = TRUE,
override.aes = list(size = 3),
title = "Mutations",
label.hjust = 0,
keywidth = 0.4,
keyheight = 0.8),
shape = guide_legend(ncol = 1,
byrow = TRUE,
title = "Mutations",
label.hjust = 0,
keywidth = 0.4,
keyheight = 0.8))
# create a legend grob
legend = g_legend(g + t.leg + leg.guides)
## arranget the main plot and the legend
# using grid graphics
gt <- ggplot_gtable(ggplot_build(g + theme(legend.position = 'none')))
gt$layout$clip[gt$layout$name == "panel"] <- "off"
grid.arrange(gt, legend,
ncol=2, nrow=1, widths=c(0.92,0.08))
## --------------------------------------------------------------------------
## lm(viab ~ mutation + pretreatment.status)
pvals = dr %>% group_by(drug, mutation) %>%
do(tidy(lm(viab ~ mut.value + pretreat, data = .))) %>%
filter(term == 'mut.value1') %>%
dplyr::select(drug, mutation, p.value)
# compute the FDR threshold
fd.thresh = 10
padj = p.adjust(pvals$p.value, method = "BH")
fdr = max(pvals$p.value[which(padj <= fd.thresh/100)])
pvals = filter(pvals, mutation != 'IGHV')
pvals = pvals %>% ungroup() %>%
mutate(mutation = ifelse(p.value > fdr,
paste0("Below ", fd.thresh,"% FDR"),
mutation)) %>%
mutate(mutation = ifelse(!(mutation %in% select.mutations) &
!(mutation == paste0("Below ",
fd.thresh,"% FDR")),
"Other", mutation)) %>%
filter(drug != "bortezomib" & drug != "NSC 74859")
lev = get.labels(pvals)
folge = get.mutation.order(lev)
# get the drug order
dr.order = get.drug.order(pvals, drugs)
mut.order = folge[!is.na(folge)]
## ----pvalsSupp, fig.path=plotDir, dev=c("png", "pdf"), fig.width=14, fig.height=10----
#FIG# S19
## plot the p-values
g = plot.pvalues(pvals, dr.order, mut.order,
colors[which(!is.na(folge))], shapes = c(rep(16,9), c(1,1)))
## add FDR threshold
g = g + geom_hline(yintercept = -log10(fdr),
linetype="dashed", size=0.3)
g = g +
annotation_custom(grob = textGrob(label = paste0("FDR", fd.thresh, "%"),
hjust = 1, vjust = 1,
gp = gpar(cex = 0.5,
fontface = "bold",
fontsize = 25)),
ymin = -log10(fdr) - 0.2,
ymax = -log10(fdr) + 0.5,
xmin = -1.3, xmax = 1.5) +
theme(legend.position = "none")
# generate pathway/target annotations for certain drug classes
#g = make.annot(g, dr.order)
# legend guide
leg.guides <- guides(colour = guide_legend(ncol = 1,
byrow = TRUE,
override.aes = list(size = 3),
title = "Mutations",
label.hjust = 0,
keywidth = 0.4,
keyheight = 0.8),
shape = guide_legend(ncol = 1,
byrow = TRUE,
title = "Mutations",
label.hjust = 0,
keywidth = 0.4,
keyheight = 0.8))
# create a legend grob
legend = g_legend(g + t.leg + leg.guides)
## arranget the main plot and the legend
# using grid graphics
gt <- ggplot_gtable(ggplot_build(g + theme(legend.position = 'none')))
gt$layout$clip[gt$layout$name == "panel"] <- "off"
grid.arrange(gt, legend,
ncol=2, nrow=1, widths=c(0.92,0.08))
## --------------------------------------------------------------------------
pvals.main = dr %>% group_by(drug, mutation) %>%
do(tidy(t.test(viab ~ mut.value, data = ., var.equal = T))) %>%
dplyr::select(drug, mutation, p.value)
p.main.adj = p.adjust(pvals.main$p.value, method = "BH")
fdr.main = max(pvals.main$p.value[which(p.main.adj <= fd.thresh/100)])
pvals.main = filter(pvals.main, mutation != "IGHV") %>%
rename(p.main = p.value)
## lm(viab ~ mutation + pretreatment.status)
pvals.sup = dr %>% group_by(drug, mutation) %>%
do(tidy(lm(viab ~ mut.value + pretreat, data = .))) %>%
filter(term == 'mut.value1') %>%
dplyr::select(drug, mutation, p.value)
p.sup.adj = p.adjust(pvals.sup$p.value, method = "BH")
fdr.sup = max(pvals.sup$p.value[which(p.sup.adj <= fd.thresh/100)])
pvals.sup = filter(pvals.sup, mutation != "IGHV") %>%
rename(p.sup = p.value)
pvals = inner_join(pvals.main, pvals.sup)
pvals = mutate(pvals, signif = ifelse(p.main > fdr.main,
ifelse(p.sup > fdr.sup,
"Below 10% FDR in both models",
"Significant with pretreatment accounted"),
ifelse(p.sup > fdr.sup,
"Significant without pretreatment in the model",
"Significant in both models")))
t2<-theme(
plot.background = element_blank(),
panel.grid.major = element_line(),
panel.grid.major.x = element_line(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(size=.4),
axis.line.x = element_line(),
axis.line.y = element_line(),
axis.text.x = element_text(size=12),
axis.text.y = element_text(size = 12),
axis.title.x = element_text(face="bold", size=12),
axis.title.y = element_text(face="bold", size=12),
legend.title = element_text(face='bold',
hjust = 1, size=10),
legend.position = c(0.78, 0.11),
legend.key = element_blank(),
legend.text = element_text(size=10),
legend.background = element_rect(color = "black")
)
## ----pvalComparisonScatterplot, fig.path=plotDir, dev=c("png", "pdf"), fig.width=10, fig.height=7----
#FIG# S19
ggplot(pvals, aes(-log10(p.main), -log10(p.sup), colour = factor(signif))) +
geom_point() + t2 + labs(x = expression(paste(-log[10], "p, pretreatment not considered", sep = "")),
y = expression(paste(-log[10], "p, accounting for pretreatment", sep = ""))) +
coord_fixed() +
scale_x_continuous(breaks = seq(0,9,by = 3)) +
scale_y_continuous(breaks = seq(0,9,by = 3)) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
scale_color_manual(name = "Statistical Significance",
values = c("#F1BB7B","#669999", "#FD6467", "#5B1A18"))
## --------------------------------------------------------------------------
signif.in.one = filter(pvals,
signif %in% c("Significant with pretreatment accounted",
"Significant without pretreatment in the model")) %>%
arrange(signif)
kable(signif.in.one, digits = 4,
align = c("l", "l", "c", "c", "c"),
col.names = c("Drug", "Mutation", "P-value (Main)",
"P-value (Supplement)", "Statistical significance"),
format.args = list(width = 14))
## ---- comment=NA, eval=FALSE-----------------------------------------------
# print(kable(signif.in.one, format = "latex", digits = 4,
# align = c("l", "l", "c", "c", "c"),
# col.names = c("Drug", "Mutation", "P-value (Main)",
# "P-value (Supplement)", "Statistical significance")))
## ---- include=!exists(".standalone"), eval=!exists(".standalone")----------
# sessionInfo()
## ---- message=FALSE, warning=FALSE, include=FALSE--------------------------
rm(list=ls())
## ---- message=FALSE, include=!exists(".standalone"), eval=!exists(".standalone")----
# library("BloodCancerMultiOmics2017")
# library("Biobase")
# library("ggbeeswarm")
# library("ggplot2")
# library("gridExtra")
# library("dplyr")
## ----echo=FALSE------------------------------------------------------------
plotDir = ifelse(exists(".standalone"), "", "part15/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
## --------------------------------------------------------------------------
data(list= c("validateExp","lpdAll"))
## --------------------------------------------------------------------------
plotTab <- filter(validateExp, Drug %in% c("Ganetespib", "Onalespib")) %>%
mutate(IGHV = Biobase::exprs(lpdAll)["IGHV Uppsala U/M", patientID]) %>%
filter(!is.na(IGHV)) %>%
mutate(IGHV = as.factor(ifelse(IGHV == 1, "M","U")),
Concentration = as.factor(Concentration))
## --------------------------------------------------------------------------
pTab <- group_by(plotTab, Drug, Concentration) %>%
do(data.frame(p = t.test(viab ~ IGHV, .)$p.value)) %>%
mutate(p = format(p, digits =2, scientific = TRUE))
## ----HSP90confirm, fig.width=14, fig.height=5, warning=FALSE, fig.path=plotDir, dev=c("png", "pdf")----
pList <- group_by(plotTab, Drug) %>%
do(plots = ggplot(., aes(x=Concentration, y = viab)) +
stat_boxplot(geom = "errorbar", width = 0.3,
position = position_dodge(width=0.6),
aes(dodge = IGHV)) +
geom_boxplot(outlier.shape = NA, position = position_dodge(width=0.6),
col="black", width=0.5, aes(dodge = IGHV)) +
geom_beeswarm(size=1,dodge.width=0.6, aes(col=IGHV)) +
theme_classic() +
scale_y_continuous(expand = c(0, 0),breaks=seq(0,1.2,0.20)) +
coord_cartesian(ylim = c(0,1.30)) +
xlab("Concentration (µM)") + ylab("Viability") +
ggtitle(unique(.$Drug)) +
geom_text(data=filter(pTab, Drug == unique(.$Drug)), y = 1.25,
aes(x=Concentration, label=sprintf("p=%s",p)),
size = 4.5) +
theme(axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text = element_text(size=15),
axis.title = element_text(size =15),
legend.text = element_text(size=13),
legend.title = element_text(size=15),
plot.title = element_text(face="bold", hjust=0.5, size=17),
plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")))
grid.arrange(grobs = pList$plots, ncol =2)
## --------------------------------------------------------------------------
plotTab <- filter(validateExp, Drug %in%
c("Cobimetinib","SCH772984","Trametinib")) %>%
mutate(Trisomy12 = Biobase::exprs(lpdAll)["trisomy12", patientID]) %>%
filter(!is.na(Trisomy12)) %>%
mutate(Trisomy12 = as.factor(ifelse(Trisomy12 == 1, "present","absent")),
Concentration = as.factor(Concentration))
## --------------------------------------------------------------------------
pTab <- group_by(plotTab, Drug, Concentration) %>%
do(data.frame(p = t.test(viab ~ Trisomy12, .)$p.value)) %>%
mutate(p = format(p, digits =2, scientific = FALSE))
## ----tris12confirm, fig.width=8, fig.height=12, warning=FALSE, fig.path=plotDir, dev=c("png", "pdf")----
pList <- group_by(plotTab, Drug) %>%
do(plots = ggplot(., aes(x=Concentration, y = viab)) +
stat_boxplot(geom = "errorbar", width = 0.3,
position = position_dodge(width=0.6),
aes(dodge = Trisomy12)) +
geom_boxplot(outlier.shape = NA, position = position_dodge(width=0.6),
col="black", width=0.5, aes(dodge = Trisomy12)) +
geom_beeswarm(size=1,dodge.width=0.6, aes(col=Trisomy12)) +
theme_classic() +
scale_y_continuous(expand = c(0, 0),breaks=seq(0,1.2,0.2)) +
coord_cartesian(ylim = c(0,1.3)) +
xlab("Concentration (µM)") + ylab("Viability") +
ggtitle(unique(.$Drug)) +
geom_text(data=filter(pTab, Drug == unique(.$Drug)), y = 1.25,
aes(x=Concentration, label=sprintf("p=%s",p)), size = 5) +
theme(axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text = element_text(size=15),
axis.title = element_text(size =15),
legend.text = element_text(size=13),
legend.title = element_text(size=15),
plot.title = element_text(face="bold", hjust=0.5, size=17),
plot.margin = unit(c(0.5,0,0.5,0), "cm")))
grid.arrange(grobs = pList$plots, ncol =1)
## ---- include=!exists(".standalone"), eval=!exists(".standalone")----------
# sessionInfo()
## ---- message=FALSE, warning=FALSE, include=FALSE--------------------------
rm(list=ls())
## ---- echo=FALSE, include=!exists(".standalone")---------------------------
knitr::opts_chunk$set(cache = TRUE)
## ---- message=FALSE, include=!exists(".standalone"), eval=!exists(".standalone")----
# library("BloodCancerMultiOmics2017")
# library("DESeq2")
# library("piano")
# library("pheatmap")
# library("genefilter")
# library("grid")
# library("gridExtra")
# library("RColorBrewer")
# library("cowplot")
# library("tidyverse")
## ----echo=FALSE------------------------------------------------------------
plotDir = ifelse(exists(".standalone"), "", "part13/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
## --------------------------------------------------------------------------
data(list=c("dds", "patmeta", "mutCOM"))
#load genesets
gmts = list(
H=system.file("extdata","h.all.v5.1.symbols.gmt",
package="BloodCancerMultiOmics2017"),
C6=system.file("extdata","c6.all.v5.1.symbols.gmt",
package="BloodCancerMultiOmics2017"),
KEGG=system.file("extdata","c2.cp.kegg.v5.1.symbols.gmt",
package="BloodCancerMultiOmics2017"))
## --------------------------------------------------------------------------
#only choose CLL samples
colData(dds)$Diagnosis <- patmeta[match(dds$PatID,rownames(patmeta)),]$Diagnosis
ddsCLL <- dds[,dds$Diagnosis %in% "CLL"]
#add trisomy 12 and IGHV information
colData(ddsCLL)$trisomy12 <-
factor(assayData(mutCOM[ddsCLL$PatID,])$binary[,"trisomy12"])
colData(ddsCLL)$IGHV <- factor(patmeta[ddsCLL$PatID,]$IGHV)
#remove samples that do not have trisomy 12 information
ddsCLL <- ddsCLL[,!is.na(ddsCLL$trisomy12)]
#how many genes and samples we have?
dim(ddsCLL)
## ---- cache=TRUE-----------------------------------------------------------
#remove genes without gene symbol annotations
ddsCLL <- ddsCLL[!is.na(rowData(ddsCLL)$symbol),]
ddsCLL <- ddsCLL[rowData(ddsCLL)$symbol != "",]
#only keep genes that have counts higher than 10 in any sample
keep <- apply(counts(ddsCLL), 1, function(x) any(x >= 10))
ddsCLL <- ddsCLL[keep,]
#Remove transcripts do not show variance across samples
ddsCLL <- estimateSizeFactors(ddsCLL)
sds <- rowSds(counts(ddsCLL, normalized = TRUE))
sh <- shorth(sds)
ddsCLL <- ddsCLL[sds >= sh,]
#variance stabilization
ddsCLL.norm <- varianceStabilizingTransformation(ddsCLL, blind=TRUE)
#how many genes left
dim(ddsCLL)
## ---- cache=TRUE-----------------------------------------------------------
design(ddsCLL) <- ~ trisomy12
ddsCLL <- DESeq(ddsCLL, betaPrior = FALSE)
DEres <- results(ddsCLL)
DEres.shr <- lfcShrink(ddsCLL, contrast = c("trisomy12","1","0"), res = DEres)
## ---- warning=FALSE--------------------------------------------------------
#FIG# S23 A
plotTab <- as.data.frame(DEres)
plotTab$onChr12 <- rowData(ddsCLL)$chromosome == 12
dosePlot <- ggplot(plotTab) +
geom_density(aes(x=log2FoldChange, col=onChr12, fill=onChr12), alpha=0.4) +
xlim( -3, 3 )
dosePlot
## --------------------------------------------------------------------------
#filter genes
fdrCut <- 0.1
fcCut <- 1.5
allDE <- data.frame(DEres.shr) %>%
rownames_to_column(var = "ID") %>%
mutate(Symbol = rowData(ddsCLL[ID,])$symbol,
Chr = rowData(ddsCLL[ID,])$chromosome) %>%
filter(padj <= fdrCut & abs(log2FoldChange) > fcCut) %>%
arrange(pvalue) %>% filter(!duplicated(Symbol)) %>%
mutate(Chr12 = ifelse(Chr == 12, "yes", "no"))
#get the expression matrix
plotMat <- assay(ddsCLL.norm[allDE$ID,])
colnames(plotMat) <- ddsCLL.norm$PatID
rownames(plotMat) <- allDE$Symbol
#sort columns of plot matrix based on trisomy 12 status
plotMat <- plotMat[,order(ddsCLL.norm$trisomy12)]
#calculate z-score and scale
plotMat <- t(scale(t(plotMat)))
plotMat[plotMat >= 4] <- 4
plotMat[plotMat <= -4] <- -4
## ---- trisomy12_heatmap, dev = c("png", "pdf"), fig.path=plotDir, fig.width = 8, fig.height = 10----
#FIG# S23 B
#prepare colums and row annotations
annoCol <- data.frame(row.names=ddsCLL.norm$PatID, Tris12=ddsCLL.norm$trisomy12)
levels(annoCol$Tris12) <- list(wt = 0, mut =1)
annoRow <- data.frame(row.names = allDE$Symbol, Chr12 = allDE$Chr12)
annoColor <- list(Tris12 = c(wt = "grey80", mut = "black"),
Chr12 = c(yes="red", no = "grey80"))
pheatmap(plotMat,
color=colorRampPalette(rev(brewer.pal(n=7, name="RdBu")))(100),
cluster_cols = FALSE,
annotation_row = annoRow, annotation_col = annoCol,
show_colnames = FALSE, fontsize_row = 3,
breaks = seq(-4,4, length.out = 101),
annotation_colors = annoColor, border_color = NA)
## --------------------------------------------------------------------------
runGSEA <- function(inputTab, gmtFile, GSAmethod="gsea", nPerm=1000){
inGMT <- loadGSC(gmtFile,type="gmt")
#re-rank by score
rankTab <- inputTab[order(inputTab[,1],decreasing = TRUE),,drop=FALSE]
if (GSAmethod == "gsea"){
#readin geneset database
#GSEA analysis
res <- runGSA(geneLevelStats = rankTab,
geneSetStat = GSAmethod,
adjMethod = "fdr", gsc=inGMT,
signifMethod = 'geneSampling', nPerm = nPerm)
GSAsummaryTable(res)
} else if (GSAmethod == "page"){
res <- runGSA(geneLevelStats = rankTab,
geneSetStat = GSAmethod,
adjMethod = "fdr", gsc=inGMT,
signifMethod = 'nullDist')
GSAsummaryTable(res)
}
}
## --------------------------------------------------------------------------
plotEnrichmentBar <- function(resTab, pCut=0.05, ifFDR=FALSE,
setName="Signatures") {
pList <- list()
rowNum <- c()
for (i in names(resTab)) {
plotTab <- resTab[[i]]
if (ifFDR) {
plotTab <- dplyr::filter(
plotTab, `p adj (dist.dir.up)` <= pCut | `p adj (dist.dir.dn)` <= pCut)
} else {
plotTab <- dplyr::filter(
plotTab, `p (dist.dir.up)` <= pCut | `p (dist.dir.dn)` <= pCut)
}
if (nrow(plotTab) == 0) {
print("No sets passed the criteria")
next
} else {
#firstly, process the result table
plotTab <- apply(plotTab, 1, function(x) {
statSign <- as.numeric(x[3])
data.frame(Name = x[1],
p = as.numeric(ifelse(statSign >= 0, x[4], x[6])),
geneNum = ifelse(statSign >= 0, x[8], x[9]),
Direction = ifelse(statSign > 0, "Up", "Down"),
stringsAsFactors = FALSE)
}) %>% do.call(rbind,.)
plotTab$Name <- sprintf("%s (%s)",plotTab$Name,plotTab$geneNum)
plotTab <- plotTab[with(plotTab,order(Direction, p, decreasing=TRUE)),]
plotTab$Direction <- factor(plotTab$Direction, levels = c("Down","Up"))
plotTab$Name <- factor(plotTab$Name, levels = plotTab$Name)
#plot the barplot
pList[[i]] <- ggplot(data=plotTab, aes(x=Name, y= -log10(p),
fill=Direction)) +
geom_bar(position="dodge",stat="identity", width = 0.5) +
scale_fill_manual(values=c(Up = "blue", Down = "red")) +
coord_fixed(ratio = 0.5) + coord_flip() + xlab(setName) +
ggtitle(i) + theme_bw() + theme(
plot.title = element_text(face = "bold", hjust =0.5),
axis.title = element_text(size=15))
rowNum <-c(rowNum,nrow(plotTab))
}
}
if (length(pList) == 0) {
print("Nothing to plot")
} else {
rowNum <- rowNum
grobList <- lapply(pList, ggplotGrob)
grobList <- do.call(rbind,c(grobList,size="max"))
panels <- grobList$layout$t[grep("panel", grobList$layout$name)]
grobList$heights[panels] <- unit(rowNum, "null")
}
return(grobList)
}
## ----message=FALSE---------------------------------------------------------
pCut <- 0.05
dataTab <- data.frame(DEres)
dataTab$ID <- rownames(dataTab)
#filter using raw pvalues
dataTab <- filter(dataTab, pvalue <= pCut) %>%
arrange(pvalue) %>%
mutate(Symbol = rowData(ddsCLL[ID,])$symbol)
dataTab <- dataTab[!duplicated(dataTab$Symbol),]
statTab <- data.frame(row.names = dataTab$Symbol, stat = dataTab$stat)
## ----fig.width=8, fig.height=6, message=FALSE------------------------------
hallmarkRes <- list()
#run PAGE
resTab <- runGSEA(statTab, gmts$H ,GSAmethod = "page")
#remove the HALLMARK_
resTab$Name <- gsub("HALLMARK_","",resTab$Name)
hallmarkRes[["Gene set enrichment analysis"]] <-
arrange(resTab,desc(`Stat (dist.dir)`))
hallBar <- plotEnrichmentBar(hallmarkRes, pCut = 0.01, ifFDR = TRUE,
setName = "Hallmark gene sets")
## ---- message=FALSE--------------------------------------------------------
keggRes <- list()
resTab <- runGSEA(statTab,gmts$KEGG,GSAmethod = "page")
#remove the KEGG_
resTab$Name <- gsub("KEGG_","",resTab$Name)
keggRes[["Gene set enrichment analysis"]] <- resTab
keggBar <- plotEnrichmentBar(keggRes, pCut = 0.01, ifFDR = TRUE,
setName = "KEGG gene sets")
## --------------------------------------------------------------------------
#select differentially expressed genes
fdrCut <- 0.05
cytoDE <- data.frame(DEres) %>% rownames_to_column(var = "ID") %>%
mutate(Symbol = rowData(ddsCLL[ID,])$symbol,
Chr=rowData(ddsCLL[ID,])$chromosome) %>%
filter(padj <= fdrCut, log2FoldChange > 0) %>%
arrange(pvalue) %>% filter(!duplicated(Symbol)) %>%
mutate(Chr12 = ifelse(Chr == 12, "yes", "no"))
#get the expression matrix
plotMat <- assay(ddsCLL.norm[cytoDE$ID,])
colnames(plotMat) <- ddsCLL.norm$PatID
rownames(plotMat) <- cytoDE$Symbol
#sort columns of plot matrix based on trisomy 12 status
plotMat <- plotMat[,order(ddsCLL.norm$trisomy12)]
#calculate z-score and sensor
plotMat <- t(scale(t(plotMat)))
plotMat[plotMat >= 4] <- 4
plotMat[plotMat <= -4] <- -4
annoCol <- data.frame(row.names = ddsCLL.norm$PatID,
Tris12 = ddsCLL.norm$trisomy12)
levels(annoCol$Tris12) <- list(wt = 0, mut =1)
annoRow <- data.frame(row.names = cytoDE$Symbol, Chr12 = cytoDE$Chr12)
## --------------------------------------------------------------------------
gsc <- loadGSC(gmts$KEGG)
geneList <- gsc$gsc$KEGG_CHEMOKINE_SIGNALING_PATHWAY
plotMat.chemo <- plotMat[rownames(plotMat) %in% geneList,]
keggHeatmap <- pheatmap(plotMat.chemo,
color = colorRampPalette(
rev(brewer.pal(n=7, name="RdBu")))(100),
cluster_cols = FALSE, clustering_method = "ward.D2",
annotation_row = annoRow, annotation_col = annoCol,
show_colnames = FALSE, fontsize_row = 8,
breaks = seq(-4,4, length.out = 101), treeheight_row = 0,
annotation_colors = annoColor, border_color = NA,
main = "CHEMOKINE_SIGNALING_PATHWAY",silent = TRUE)$gtable
## ----geneEnrichment_result, dev = c("png", "pdf"), fig.path=plotDir, fig.width = 14, fig.height = 13----
#FIG# S24 ABC
ggdraw() +
draw_plot(hallBar, 0, 0.7, 0.5, 0.3) +
draw_plot(keggBar, 0.5, 0.7, 0.5, 0.3) +
draw_plot(keggHeatmap, 0.1, 0, 0.8, 0.65) +
draw_plot_label(c("A","B","C"), c(0, 0.5, 0.1), c(1, 1, 0.7),
fontface = "plain", size=20)
## ---- include=!exists(".standalone"), eval=!exists(".standalone")----------
# sessionInfo()
## ---- message=FALSE, warning=FALSE, include=FALSE--------------------------
rm(list=ls())
## ---- message=FALSE, include=!exists(".standalone"), eval=!exists(".standalone")----
# library("BloodCancerMultiOmics2017")
# library("Biobase")
# library("SummarizedExperiment")
# library("AnnotationDbi")
# library("org.Hs.eg.db")
# library("dplyr")
# library("abind")
# library("reshape2")
# library("RColorBrewer")
# library("glmnet")
# library("ipflasso")
# library("ggplot2")
# library("grid")
# library("DESeq2")
## ----echo=FALSE------------------------------------------------------------
plotDir = ifelse(exists(".standalone"), "", "part11/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
## --------------------------------------------------------------------------
options(stringsAsFactors=FALSE)
## --------------------------------------------------------------------------
data(list=c("conctab", "drpar", "drugs", "patmeta", "lpdAll", "dds", "mutCOM",
"methData"))
## --------------------------------------------------------------------------
e<-dds
colnames(e)<-colData(e)$PatID
#only consider CLL patients
CLLPatients<-rownames(patmeta)[which(patmeta$Diagnosis=="CLL")]
#Methylation Data
methData = t(assay(methData))
#RNA Data
eCLL<-e[,colnames(e) %in% CLLPatients]
###
#filter out genes without gene namce
AnnotationDF<-data.frame(EnsembleId=rownames(eCLL),stringsAsFactors=FALSE)
AnnotationDF$symbol <- mapIds(org.Hs.eg.db,
keys=rownames(eCLL),
column="SYMBOL",
keytype="ENSEMBL",
multiVals="first")
eCLL<-eCLL[AnnotationDF$EnsembleId[!is.na(AnnotationDF$symbol)],]
#filter out low count genes
###
minrs <- 100
rs <- rowSums(assay(eCLL))
eCLL<-eCLL[ rs >= minrs, ]
#variance stabilize the data
#(includes normalizing for library size and dispsersion estimation)
vstCounts<-varianceStabilizingTransformation(eCLL)
vstCounts<-assay(vstCounts)
#no NAs in data
any(is.na(vstCounts))
#filter out low variable genes
ntop<-5000
vstCountsFiltered<-vstCounts[order(apply(vstCounts, 1, var, na.rm=T),
decreasing = T)[1:ntop],]
eData<-t(vstCountsFiltered)
#no NAs
any(is.na(eData))
#genetics
#remove features with less than 5 occurences
mutCOMbinary<-channel(mutCOM, "binary")
mutCOMbinary<-mutCOMbinary[featureNames(mutCOMbinary) %in% CLLPatients,]
genData<-Biobase::exprs(mutCOMbinary)
idx <- which(colnames(genData) %in% c("del13q14_bi", "del13q14_mono"))
genData <- genData[,-idx]
colnames(genData)[which(colnames(genData)=="del13q14_any")] = "del13q14"
minObs <- 5
#remove feutes with less than 5 occurecnes
genData<-genData[,colSums(genData, na.rm=T)>=minObs]
#IGHV
translation <- c(`U` = 0, `M` = 1)
stopifnot(all(patmeta$IGHV %in% c("U","M", NA)))
IGHVData <- matrix(translation[patmeta$IGHV],
dimnames = list(rownames(patmeta), "IGHV"), ncol = 1)
IGHVData<-IGHVData[rownames(IGHVData) %in% CLLPatients,,drop=F]
#remove patiente with NA IGHV status
IGHVData<-IGHVData[!is.na(IGHVData), ,drop=F]
any(is.na(IGHVData))
#demographics (age and sex)
patmeta<-subset(patmeta, Diagnosis=="CLL")
gender <- ifelse(patmeta[,"Gender"]=="m",0,1)
# impute missing values in age by mean
ImputeByMean <- function(x) {x[is.na(x)] <- mean(x, na.rm=TRUE); return(x)}
age<-ImputeByMean(patmeta[,"Age4Main"])
demogrData <- cbind(age=age,gender=gender)
rownames(demogrData) <- rownames(patmeta)
#Pretreatment
pretreated<-patmeta[,"IC50beforeTreatment", drop=FALSE]
##### drug viabilites
summaries <- c(paste("viaraw", 1:5, sep=".") %>% `names<-`(paste(1:5)),
`4:5` = "viaraw.4_5", `1:5` = "viaraw.1_5")
a <- do.call( abind, c( along=3, lapply( summaries,
function(x) assayData(drpar)[[x]])))
dimnames(a)[[3]] <- names(summaries)
names(dimnames(a)) <- c( "drug", "patient", "summary" )
viabData <- acast( melt(a), patient ~ drug + summary )
rownames(viabData)<-c(substr(rownames(viabData),1,4)[1:3],
substr(rownames(viabData),1,5)[4:nrow(viabData)])
## --------------------------------------------------------------------------
# common patients
Xlist<-list(RNA=eData, meth=methData, gen=genData, IGHV=IGHVData,
demographics=demogrData, drugs=viabData, pretreated=pretreated)
PatientsPerOmic<-lapply(Xlist, rownames)
sapply(PatientsPerOmic, length)
allPatients<-Reduce(union, PatientsPerOmic)
PatientOverview<-sapply(Xlist, function(M) allPatients %in% rownames(M))
Patients <- (1:nrow(PatientOverview))
Omics <- (1:ncol(PatientOverview))
image(Patients,Omics, PatientOverview*1, axes=F, col=c("white", "black"),
main="Sample overview across omics")
axis(2, at = 1:ncol(PatientOverview), labels=colnames(PatientOverview), tick=F)
commonPatients<-Reduce(intersect, PatientsPerOmic)
length(commonPatients)
XlistCommon<-lapply(Xlist, function(data) data[commonPatients,, drop=F])
#Take care of missing values (present in genetic data)
ImputeByMean <- function(x) {x[is.na(x)] <- mean(x, na.rm=TRUE); return(x)}
#NAs in genetic
#remove feauters with less 90% completeness
RarlyMeasuredFeautres<-
which(colSums(is.na(XlistCommon$gen))>0.1*nrow(XlistCommon$gen))
XlistCommon$gen<-XlistCommon$gen[,-RarlyMeasuredFeautres]
#remove patients with less than 90% of genetic feautres measured
IncompletePatients<-
rownames(XlistCommon$gen)[
(rowSums(is.na(XlistCommon$gen))>0.1*ncol(XlistCommon$gen))]
commonPatients<-commonPatients[!commonPatients %in% IncompletePatients]
XlistCommon<-lapply(XlistCommon, function(data) data[commonPatients,, drop=F])
#replace remaining NA by mean and round to 0 or 1
XlistCommon$gen<-round(apply(XlistCommon$gen, 2, ImputeByMean))
#NAs in methylation
#remove feauters with less 90% completeness
XlistCommon$meth<-
XlistCommon$meth[,colSums(is.na(XlistCommon$meth))<0.1*nrow(methData)]
#impute remainin missing values by mean for each feautre across patients
XlistCommon$meth<-(apply(XlistCommon$meth, 2, ImputeByMean))
#final dimensions of the data
sapply(XlistCommon, dim)
## ---- fig.width=12, fig.height=10------------------------------------------
pcaMeth<-prcomp(XlistCommon$meth, center=T, scale. = F)
XlistCommon$MethPCs<-pcaMeth$x[,1:20]
colnames(XlistCommon$MethPCs)<-
paste("meth",colnames(XlistCommon$MethPCs), sep="")
pcaExpr<-prcomp(XlistCommon$RNA, center=T, scale. = F)
XlistCommon$RNAPCs<-pcaExpr$x[,1:20]
colnames(XlistCommon$RNAPCs)<-paste("RNA",colnames(XlistCommon$RNAPCs), sep="")
## --------------------------------------------------------------------------
DOI <- c("D_006_1:5", "D_010_1:5", "D_159_1:5","D_002_4:5", "D_003_4:5",
"D_012_4:5", "D_063_4:5", "D_166_4:5")
drugviab<-XlistCommon$drugs
drugviab<-drugviab[,DOI, drop=F]
colnames(drugviab) <- drugs[substr(colnames(drugviab),1,5),"name"]
## --------------------------------------------------------------------------
ZPCs<-list(expression=XlistCommon$RNAPCs,
genetic=XlistCommon$gen,
methylation= XlistCommon$MethPCs,
demographics=XlistCommon$demographics,
IGHV=XlistCommon$IGHV,
pretreated = XlistCommon$pretreated)
ZPCs$all<-do.call(cbind, ZPCs)
ZPCsunscaled<-ZPCs
ZPCsscaled<-lapply(ZPCs, scale)
lapply(ZPCsscaled, colnames)
## --------------------------------------------------------------------------
set1 <- brewer.pal(9,"Set1")
colMod<-c(paste(set1[c(4,1,5,3,2,7)],"88",sep=""), "grey")
names(colMod) <-
c("demographics", "genetic", "IGHV","expression", "methylation", "pretreated",
"all")
## ---- echo=F---------------------------------------------------------------
#Function to calculate Var Explained for Penalized Regression
R2ForPenRegPlusViz<-function(Z, drugviab, nfolds=10, alpha=1, nrep=100,
Parmfrow=c(2,4), ylimMax=0.4, standardize=TRUE){
Zlocal<-Z
set.seed(1030)
seeds<-sample(1:10000000, nrep)
RepeatedR2list<-lapply(1:nrep, function(outer){
#Use same folds for all omics to make comparable
set.seed(seeds[outer])
foldsAssignment<-sample(rep(seq(nfolds), length=nrow(drugviab)))
R2echOmicadj<-sapply(colnames(drugviab), function(dname) {
d<-drugviab[,dname]
sapply(names(Zlocal), function(nameZ) {
pred<-Zlocal[[nameZ]]
#fit a lasso model for omics with more than one features
if(ncol(pred)>1){
fitcv<-cv.glmnet(pred,d,family="gaussian", standardize=standardize,
alpha=alpha, foldid=foldsAssignment)
R2 <- 1-min(fitcv$cvm)/fitcv$cvm[1]
}else {#fit a liner model for single feautres (IGHV)
fitlm<-lm(d~., data.frame(pred))
R2<- summary(fitlm)$r.squared
}
R2
})
}
)
})
RepeatedR2<-RepeatedR2list
#calculate mean and sd across repitions wiht different folds
meanR2<-apply(simplify2array(RepeatedR2), 1:2, mean)
sdR2<-apply(simplify2array(RepeatedR2), 1:2, sd)
par(mfrow=Parmfrow, mar=c(7,5,3,3))
for (i in 1: ncol(meanR2)) {
barc<-barplot(meanR2[,i], main= colnames(meanR2)[i], ylim=c(0,ylimMax),
las=2, col=colMod[rownames(meanR2)], ylab="R2")
segments(barc, meanR2[,i]-sdR2[,i],barc, meanR2[,i]+sdR2[,i])
}
RepeatedR2
}
## ----lasso_main, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=9, fig.height=8----
#FIG# 5A
resultLassoPCA<-R2ForPenRegPlusViz(ZPCsscaled, drugviab, nfold=10, alpha=1,
nrep=100, ylimMax=0.4)
df_resultLassoPCA <- melt(resultLassoPCA)
colnames(df_resultLassoPCA) <- c("omic", "drug", "R2", "run")
summaryR2 <- df_resultLassoPCA %>% group_by(omic, drug) %>%
dplyr::summarise(meanR2=mean(R2),sdR2 = sd(R2), nR2 = length(R2)) %>%
mutate(seR2 = sdR2/sqrt(nR2))
ggplot(summaryR2, aes(x=omic, y=meanR2, fill=omic, group=omic))+
geom_bar(stat = "identity") + scale_fill_manual(values=colMod) +
geom_errorbar(aes(ymax = meanR2 + sdR2,ymin = meanR2 - sdR2), position = "dodge", width = 0.25) +facet_wrap(~drug, ncol=4) +theme_bw(base_size = 18) +ylab(bquote(R^2)) + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.text = element_text(colour = "black")) +xlab("")+guides(fill=guide_legend(title="Data type"))
## ---- eval=FALSE, fig.path='supp_'-----------------------------------------
# nfolds<-10
# nrep<-100
#
# DOI <-
# grepl("1:5",colnames(XlistCommon$drugs)) |
# grepl("4:5",colnames(XlistCommon$drugs))
# drugviabAll<-XlistCommon$drugs
# drugviabAll<-drugviabAll[,DOI]
# colnames(drugviabAll) <-
# paste0(drugs[substr(colnames(drugviabAll),1,5),"name"],
# substr(colnames(drugviabAll),6,9))
#
# R2ForPenReg(Zscaled, drugviabAll, nfold=nfolds, alpha=1, nrep=nrep,
# Parmfrow=c(4,4), ylimMax=0.6)
## ----echo=FALSE------------------------------------------------------------
dataType = c(M="Methylation_Cluster", V="viab", G="gen", I="IGHV", P="pretreat")
## ----echo=FALSE------------------------------------------------------------
coldef<-list()
coldef["I"]<-brewer.pal(9, "Blues")[7]
coldef["M"]<-list(brewer.pal(9, "Blues")[c(1, 5, 9)])
coldef["G"]<-brewer.pal(8, "YlOrRd")[8]
coldef["P"]<-"chocolate4"
## --------------------------------------------------------------------------
lpdCLL = lpdAll[ , lpdAll$Diagnosis=="CLL"]
## --------------------------------------------------------------------------
lpdCLL = lpdAll[ , lpdAll$Diagnosis=="CLL"]
lpdCLL = BloodCancerMultiOmics2017:::prepareLPD(lpd=lpdCLL, minNumSamplesPerGroup=5)
(predictorList = BloodCancerMultiOmics2017:::makeListOfPredictors(lpdCLL))
## --------------------------------------------------------------------------
drs = list("1:5"=c("D_006", "D_010", "D_159"),
"4:5"=c("D_002", "D_003", "D_012", "D_063", "D_166"))
## --------------------------------------------------------------------------
predvar = unlist(BloodCancerMultiOmics2017:::makePredictions(drs=drs,
lpd=lpdCLL,
predictorList=predictorList,
lim=0.15, std=FALSE, adaLasso=TRUE,
colors=coldef),
recursive=FALSE)
## ----echo=FALSE------------------------------------------------------------
details = function(dr, what) {
predvar[[dr]][["plot"]][[what]]
}
## ----prediction-D_006, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=details("D_006","width"), fig.height=details("D_006","height"), eval=TRUE----
#FIG# 5B
grid.draw(details("D_006","plot"))
## ----prediction-D_010, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=details("D_010","width"), fig.height=details("D_010","height"), eval=TRUE----
#FIG# 5B
grid.draw(details("D_010","plot"))
## ----prediction-D_159, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=details("D_159","width"), fig.height=details("D_159","height"), eval=TRUE----
#FIG# 5B
grid.draw(details("D_159","plot"))
## ----prediction-D_002, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=details("D_002","width"), fig.height=details("D_002","height"), eval=TRUE----
#FIG# 5B
grid.draw(details("D_002","plot"))
## ----prediction-D_003, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=details("D_003","width"), fig.height=details("D_003","height"), eval=TRUE----
#FIG# 5B
grid.draw(details("D_003","plot"))
## ----prediction-D_012, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=details("D_012","width"), fig.height=details("D_012","height"), eval=TRUE----
#FIG# 5B
grid.draw(details("D_012","plot"))
## ----prediction-D_063, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=details("D_063","width"), fig.height=details("D_063","height"), eval=TRUE----
#FIG# 5B
grid.draw(details("D_063","plot"))
## ----prediction-D_166, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=details("D_166","width"), fig.height=details("D_166","height"), eval=TRUE----
#FIG# 5B
grid.draw(details("D_166","plot"))
## --------------------------------------------------------------------------
legends = BloodCancerMultiOmics2017:::makeLegends(legendFor=c("G","I","M", "P"),
coldef)
## ----legend, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=legends[["width"]], fig.height=legends[["height"]]----
#FIG# 5B legend
grid.draw(legends[["plot"]])
## --------------------------------------------------------------------------
drs_rot = list("4:5"=c("D_067"))
predvar_rot = unlist(BloodCancerMultiOmics2017:::makePredictions(drs=drs_rot,
lpd=lpdCLL,
predictorList=predictorList,
lim=0.23, std=FALSE, adaLasso=TRUE,
colors=coldef),
recursive=FALSE)
## ----prediction-D_067, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=predvar_rot[["D_067"]][["plot"]][["width"]], fig.height=predvar_rot[["D_067"]][["plot"]][["height"]], eval=TRUE----
#FIG# S26
grid.draw(predvar_rot[["D_067"]][["plot"]][["plot"]])
## --------------------------------------------------------------------------
alldrs = unique(fData(lpdCLL)[fData(lpdCLL)$type=="viab","id"])
drs = list("1:5"=alldrs, "4:5"=alldrs)
## ---- eval=TRUE------------------------------------------------------------
predvar2 = BloodCancerMultiOmics2017:::makePredictions(drs=drs,
lpd=lpdCLL,
predictorList=predictorList,
lim=0.23,
colors=coldef)
## --------------------------------------------------------------------------
givePreatreatSum = function(predNum) {
idx = sapply(predvar2[[predNum]], function(x) length(x)==1)
predvar2[[predNum]] = predvar2[[predNum]][!idx]
# get model coefficients and reshape
coeffs <- do.call(cbind,lapply(predvar2[[predNum]], "[[", 'coeffs'))
coeffs <- coeffs[-1,]
coeffs <- as.matrix(coeffs)
# colnames(coeffs) <- unlist(drs["1:5"])
colnames(coeffs) = names(predvar2[[predNum]])
colnames(coeffs) <- drugs[colnames(coeffs),"name"]
coeffDF <- melt(as.matrix(coeffs))
colnames(coeffDF) <- c("predictor", "drug", "beta")
coeffDF$selected <- coeffDF$beta !=0
#sort by times being selected
coeffDF$predictor <- factor(coeffDF$predictor, level=)
# number of drugs a predictor is chosen for
gg1 <- coeffDF %>% group_by(predictor) %>%
dplyr::summarize(selectedSum = sum(selected)) %>%
mutate(predictor = factor(predictor,
levels=predictor[order(selectedSum)])) %>%
ggplot(aes(x=predictor, y=selectedSum)) +
geom_bar(stat="identity")+ylab("# drugs selected for") +
coord_flip()
# boxplots of non-zero coeffients
orderMedian <- filter(coeffDF, selected) %>% group_by(predictor) %>%
dplyr::summarize(medianBeta = median(abs(beta)))
coeffDF$predictor <- factor(
coeffDF$predictor,
levels=orderMedian$predictor[order(orderMedian$medianBeta)] )
gg2 <- ggplot(filter(coeffDF, selected), aes(x=predictor, y=abs(beta))) +
geom_boxplot() +
coord_flip() + ggtitle("Distribution of non-zero coefficients")
gridExtra::grid.arrange(gg1,gg2, ncol=1)
# coefficeints per drug
ggplot(filter(coeffDF, selected),
aes(x= drug, y=abs(beta), col= predictor=="Pretreatment")) +
geom_point() +
coord_flip()
#drugs pretreatment is selected for
as.character(filter(coeffDF, predictor=="Pretreatment" & beta!=0)$drug)
PselDrugs <- as.character(
filter(coeffDF, predictor=="Pretreatment" & beta!=0)$drug)
length(PselDrugs)
# length(drs[[1]])
}
## ----pretreatment_c1-5, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=7, fig.height=10, eval=TRUE----
givePreatreatSum(predNum=1)
## ----pretreatment_c4-5, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=7, fig.height=10, eval=TRUE----
givePreatreatSum(predNum=2)
## ---- include=!exists(".standalone"), eval=!exists(".standalone")----------
# sessionInfo()
## ---- message=FALSE, warning=FALSE, include=FALSE--------------------------
rm(list=ls())
## ---- message=FALSE, include=!exists(".standalone"), eval=!exists(".standalone")----
# library("BloodCancerMultiOmics2017")
# library("Biobase")
# library("RColorBrewer")
# library("grid")
# library("ggplot2")
# library("survival")
# library("gtable")
# library("forestplot")
# library("xtable")
# library("maxstat")
## ----echo=FALSE------------------------------------------------------------
plotDir = ifelse(exists(".standalone"), "", "part08/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
## --------------------------------------------------------------------------
options(stringsAsFactors=FALSE)
## --------------------------------------------------------------------------
data(lpdAll, patmeta, drugs)
## --------------------------------------------------------------------------
lpdCLL <- lpdAll[ , lpdAll$Diagnosis=="CLL" ]
# data rearrangements
survT = patmeta[colnames(lpdCLL),]
survT[which(survT[,"IGHV"]=="U") ,"IGHV"] = 0
survT[which(survT[,"IGHV"]=="M") ,"IGHV"] = 1
survT$IGHV = as.numeric(survT$IGHV)
colnames(survT) = gsub("Age4Main", "age", colnames(survT))
survT$ibr45 <- 1-Biobase::exprs(lpdCLL)[ "D_002_4:5", rownames(survT) ]
survT$ide45 <- 1-Biobase::exprs(lpdCLL)[ "D_003_4:5", rownames(survT) ]
survT$prt45 <- 1-Biobase::exprs(lpdCLL)[ "D_166_4:5", rownames(survT) ]
survT$selu45 <- 1-Biobase::exprs(lpdCLL)[ "D_012_4:5", rownames(survT) ]
survT$ever45 <- 1-Biobase::exprs(lpdCLL)[ "D_063_4:5", rownames(survT) ]
survT$nut15 <- 1-Biobase::exprs(lpdCLL)[ "D_010_1:5", rownames(survT) ]
survT$dox15 <- 1-Biobase::exprs(lpdCLL)[ "D_159_1:5", rownames(survT) ]
survT$flu15 <- 1-Biobase::exprs(lpdCLL)[ "D_006_1:5", rownames(survT) ]
survT$SF3B1 <- Biobase::exprs(lpdCLL)[ "SF3B1", rownames(survT) ]
survT$NOTCH1 <- Biobase::exprs(lpdCLL)[ "NOTCH1", rownames(survT) ]
survT$BRAF <- Biobase::exprs(lpdCLL)[ "BRAF", rownames(survT) ]
survT$TP53 <- Biobase::exprs(lpdCLL)[ "TP53", rownames(survT) ]
survT$del17p13 <- Biobase::exprs(lpdCLL)[ "del17p13", rownames(survT) ]
survT$del11q22.3 <- Biobase::exprs(lpdCLL)[ "del11q22.3", rownames(survT) ]
survT$trisomy12 <- Biobase::exprs(lpdCLL)[ "trisomy12", rownames(survT) ]
survT$IGHV_cont <- patmeta[ rownames(survT) ,"IGHV Uppsala % SHM"]
# competinting risk endpoint fpr
survT$compE <- ifelse(survT$treatedAfter == TRUE, 1, 0)
survT$compE <- ifelse(survT$treatedAfter == FALSE & survT$died==TRUE,
2, survT$compE )
survT$T7 <- ifelse(survT$compE == 1, survT$T5, survT$T6 )
## ----forest----------------------------------------------------------------
forest <- function(Time, endpoint, title, sdrugs, split, sub) {
stopifnot(is.character(Time), is.character(title), is.character(split),
is.character(endpoint),
all(c(Time, split, endpoint) %in% colnames(survT)),
is.logical(survT[[endpoint]]),
is.character(sdrugs), !is.null(names(sdrugs)))
clrs <- fpColors(box="royalblue",line="darkblue", summary="royalblue")
res <- lapply(sdrugs, function(g) {
drug <- survT[, g] * 10
suse <- if (identical(sub, "none"))
rep(TRUE, nrow(survT))
else
(survT[[split]] == sub)
stopifnot(sum(suse, na.rm = TRUE) > 1)
surv <- coxph(Surv(survT[,Time], survT[,endpoint]) ~ drug, subset=suse)
sumsu <- summary(surv)
c(p = sumsu[["coefficients"]][, "Pr(>|z|)"],
coef = sumsu[["coefficients"]][, "exp(coef)"],
lower = sumsu[["conf.int"]][, "lower .95"],
higher = sumsu[["conf.int"]][, "upper .95"])
})
s <- do.call(rbind, res)
rownames(s) <- names(sdrugs)
tabletext <- list(c(NA, rownames(s)), append(list("p-value"),
sprintf("%.4f", s[,"p"])))
forestplot(tabletext,
rbind(
rep(NA, 3),
s[, 2:4]),
page = new,
clip = c(0.8,20),
xlog = TRUE, xticks = c(0.5,1, 1.5), title = title,
col = clrs,
txt_gp = fpTxtGp(ticks = gpar(cex=1) ),
new_page = TRUE)
}
## ----forest-together-------------------------------------------------------
com <- function( Time, endpoint, scaleX, sub, d, split, drug_names) {
res <- lapply(d, function(g) {
drug <- survT[,g] * scaleX
## all=99, M-CLL=1, U-CLL=0
if(sub==99) { surv <- coxph(Surv(survT[,paste0(Time)],
survT[,paste0(endpoint)] == TRUE) ~ drug)}
if(sub<99) { surv <- coxph(Surv(survT[,paste0(Time)],
survT[,paste0(endpoint)] == TRUE) ~ drug,
subset=survT[,paste0(split)]==sub)}
c(summary(surv)[[7]][,5], summary(surv)[[7]][,2],
summary(surv)[[8]][,3],
summary(surv)[[8]][,4])
})
s <- do.call(rbind, res)
colnames(s) <- c("p", "HR", "lower", "higher")
rownames(s) <- drug_names
s
}
fp <- function( sub, title, d, split, drug_names, a, b, scaleX) {
ttt <- com(Time="T5", endpoint="treatedAfter", sub=sub, d=d,
split=split, drug_names=drug_names, scaleX=scaleX)
rownames(ttt) <- paste0(rownames(ttt), "_TTT")
os <- com(Time="T6", endpoint="died", sub=sub, d=d, split=split,
drug_names=drug_names, scaleX=scaleX)
rownames(os) <- paste0(rownames(os), "_OS")
n <- c( p=NA, HR=NA, lower=NA, higher=NA )
nn <- t( data.frame( n ) )
for (i in 1:(nrow(ttt)-1) ) { nn <-rbind(nn, n ) }
rownames(nn) <- drug_names
od <- order( c(seq(nrow(nn)), seq(nrow(ttt)), seq(nrow(os)) ))
s <- data.frame( rbind(nn, ttt, os)[od, ] )
s$Name <- rownames(s)
s$x <- 1:nrow(s)
s$col <- rep(c("white", "black", "darkgreen"), nrow(ttt) )
s$Endpoint <- factor( c(rep("nn", nrow(nn) ), rep("TTT", nrow(ttt) ),
rep("OS", nrow(os) ) )[od] )
s$features <- ""; s[ which(s$Endpoint=="OS"),"features"] <- drug_names
s[which(s$Endpoint=="nn"), "Endpoint"] <- ""
s <- rbind(s, rep(NA, 8))
p <- ggplot(data=s ,aes(x=x, y=HR, ymin=lower, ymax=higher,
colour=Endpoint)) + geom_pointrange() +
theme(legend.position="top", legend.text = element_text(size = 20) ) +
scale_x_discrete(limits=s$x, labels=s$features ) +
expand_limits(y=c(a,b)) +
scale_y_log10(breaks=c(0.01,0.1,0.5,1,2,5,10),
labels=c(0.01,0.1,0.5,1,2,5,10)) +
theme(
panel.grid.minor = element_blank(),
axis.title.x = element_text(size=16),
axis.text.x = element_text(size=16, colour="black"),
axis.title.y = element_blank(),
axis.text.y = element_text(size=12, colour="black"),
legend.key = element_rect(fill = "white"),
legend.background = element_rect(fill = "white"),
legend.title = element_blank(),
panel.background = element_rect(fill = "white", color="black"),
panel.grid.major = element_blank(),
axis.ticks.y = element_blank()
) +
coord_flip() +
scale_color_manual(values=c("OS"="darkgreen", "TTT"="black"),
labels=c("OS", "TTT", "")) +
geom_hline(aes(yintercept=1), colour="black", size=1.5,
linetype="dashed", alpha=0.3) +
annotate("text", x = 1:nrow(s)+0.5, y = s$HR+0.003,
label = ifelse( s$p<0.001, paste0("p<","0.001"),
paste0("p=", round(s$p,3) ) ), colour=s$col)
plot(p)
}
## ----Fig5A, fig.path=plotDir, fig.width=5.55, fig.height=( (1+8*1.2) ), dev = c("png", "pdf"), warning=FALSE----
#FIG# S27
d <- c("SF3B1", "NOTCH1", "BRAF", "TP53", "del17p13", "del11q22.3",
"trisomy12", "IGHV")
drug_names <- c("SF3B1", "NOTCH1", "BRAF", "TP53", "del17p13", "del11q22.3",
"Trisomy12" ,"IGHV")
fp(sub=99, d=d, drug_names=drug_names, split="IGHV", title="", a=0, b=10,
scaleX=1)
## ----Fig5B, fig.path=plotDir, fig.width=6.0, fig.height=( (1+8*1.2) ), dev = c("png", "pdf"), warning=FALSE----
#FIG# 6A
d <- c("flu15", "nut15", "dox15", "ibr45", "ide45", "prt45", "selu45",
"ever45")
drug_names <- c("Fludarabine", "Nutlin-3", "Doxorubicine", "Ibrutinib",
"Idelalisib", "PRT062607 HCl", "Selumetinib" ,"Everolimus")
fp(sub=99, d=d, drug_names=drug_names, split="TP53", title="", a=0, b=5,
scaleX=10)
## ----SFig_genetics1, fig.path=plotDir, fig.width = 8.46, fig.height = 4.5, dev = c("png", "pdf")----
#FIG# S27 left (top+bottom)
par(mfcol=c(1,2))
for (fac in paste(c("IGHV", "TP53"))) {
survplot( Surv(survT$T5, survT$treatedAfter == TRUE) ~ as.factor(survT[,fac]),
snames=c("wt", "mut"),
lwd=1.5, cex.axis = 1, cex.lab=1, col= c("darkmagenta", "dodgerblue4"),
show.nrisk = FALSE,
legend.pos = FALSE, stitle = "", hr.pos= "topright",
main = paste(fac),
xlab = 'Time (Years)', ylab = 'Time to treatment')
}
## ----SFig_genetics2, fig.path=plotDir, fig.width = 8.46, fig.height = 4.5, dev = c("png", "pdf")----
#FIG# 6B left
#FIG# S27 right (top+bottom)
par(mfcol=c(1,2))
for (fac in paste(c("IGHV", "TP53"))) {
survplot( Surv(survT$T6, survT$died == TRUE) ~ as.factor(survT[,fac]),
snames=c("wt", "mut"),
lwd=1.5, cex.axis = 1.0, cex.lab=1.0, col= c("darkmagenta", "dodgerblue4"),
show.nrisk = FALSE,
legend.pos = FALSE, stitle = "", hr.pos= "bottomleft",
main = paste(fac),
xlab = 'Time (Years)', ylab = 'Overall survival')
}
## ----KM, echo=FALSE--------------------------------------------------------
km <- function(drug, split, title, t, hr, c) {
stopifnot(is.character(drug), length(drug)==1, is.character(title),
length(title)==3)
surv <- survT[!(is.na(survT[,split])), ]
k <- Biobase::exprs(lpdCLL)[ drug, rownames(surv) ]
ms5 <- maxstat.test(Surv(T5, treatedAfter) ~ k,
data = surv,
smethod = "LogRank",
minprop = 0.2,
maxprop = 0.8,
alpha = NULL)
ms6 <- maxstat.test(Surv(T6, died) ~ k,
data = surv,
smethod = "LogRank",
minprop = 0.2,
maxprop = 0.8,
alpha = NULL)
# median & TTT
if (c=="med" & t=="TTT") {
surv$cutA <- ifelse(
k >= median( k[which(!(is.na(surv$T5)) ) ] ), "weak", "good")
surv$cutM <- ifelse(
k >= median( k[ which( surv[,paste0(split)]==1 & !(is.na(surv$T5)) ) ],
na.rm=TRUE ), "weak", "good")
surv$cutU <- ifelse(
k >= median( k[ which( surv[,paste0(split)]==0 & !(is.na(surv$T5)) ) ],
na.rm=TRUE ), "weak", "good")
}
# median & OS
if (c=="med" & t=="OS") {
surv$cutA <- ifelse(k >= median(k), "weak", "good")
surv$cutM <- ifelse(k >= median( k[ which( surv[,paste0(split)]==1 ) ] ),
"weak", "good")
surv$cutU <- ifelse(k >= median( k[ which( surv[,paste0(split)]==0 ) ] ),
"weak", "good")
}
#TTT & maxstat
if (c=="maxstat" & t=="TTT") {
surv$cutA <- surv$cut5 <- ifelse(k >= ms5$estimate, "weak", "good")
surv$cutM <- surv$cut5 <- ifelse(k >= ms5$estimate, "weak", "good")
surv$cutU <- surv$cut5 <- ifelse(k >= ms5$estimate, "weak", "good")
}
#OS & maxstat
if (c=="maxstat" & t=="OS") {
surv$cutA <- surv$cut5 <- ifelse(k >= ms6$estimate, "weak", "good")
surv$cutM <- surv$cut5 <- ifelse(k >= ms6$estimate, "weak", "good")
surv$cutU <- surv$cut5 <- ifelse(k >= ms6$estimate, "weak", "good")
}
drName <- toCaps(drugs[stripConc(drug), "name"])
sp <- function(...)
survplot(...,
lwd = 3, cex.axis = 1.2, cex.lab = 1.5,
col= c("royalblue", "darkred"), show.nrisk = FALSE,
legend.pos = FALSE, stitle = "",
hr.pos=ifelse(hr=="bl", "bottomleft", "topright" ),
xlab = 'Time (Years)')
if (t=="TTT") {
yl <- "Fraction w/o treatment"
if (c=="med"){
cat(sprintf("%s median-cutpoint for TTT: %5.2g\n",
drName, median(k) ) ) } else
{ cat(sprintf("%s cutpoint for TTT: %5.2g\n", drName,
ms5$estimate )) }
sp(Surv(surv$T5, surv$treatedAfter) ~ surv$cutA,
subset = rep(TRUE, nrow(surv)), ylab = yl, main = drName)
sp(Surv(surv$T5, surv$treatedAfter) ~ surv$cutM,
subset = surv[, split]==1, ylab = yl,
main = paste(drName, title[1], title[3]))
sp(Surv(surv$T5, surv$treatedAfter) ~ surv$cutU,
subset = surv[ ,split]==0, ylab = yl,
main = paste(drName, title[1], title[2])) }
# OS
else {
yl <- "Fraction overall survival"
if (c=="med"){
cat(sprintf("%s median-cutpoint for OS: %5.2g\n",
drName, median(k) ) ) } else {
cat(sprintf("%s cutpoint for OS: %5.2g\n",
drName, ms6$estimate ))}
sp(Surv(surv$T6, surv$died) ~ surv$cutA,
subset = rep(TRUE, nrow(surv)), ylab = yl, main = drName)
sp(Surv(surv$T6, surv$died) ~ surv$cutM,
subset = surv[, split]==1, ylab = yl,
main = paste(drName, title[1], title[3]))
sp(Surv(surv$T6, surv$died) ~ surv$cutU,
subset = surv[ ,split]==0, ylab = yl,
main = paste(drName, title[1], title[2]))
}
}
## ----KM-TTT-maxstat, fig.path=plotDir, fig.width = 10, fig.height = 3.3, dev = c("png", "pdf")----
par(mfrow=c(1,3), mar=c(5,5,2,0.9))
km(drug = "D_006_1:5", split = "TP53", t="TTT",
title=c("(TP53", "wt)", "mut)"), hr="tr", c="maxstat")
km(drug = "D_159_1:5", split = "TP53", t="TTT",
title=c("(TP53", "wt)", "mut)"), hr="tr", c="maxstat")
km(drug = "D_010_1:5", split = "TP53", t="TTT",
title=c("(TP53", "wt)", "mut)"), hr="tr", c="maxstat")
km(drug = "D_002_4:5", split = "IGHV", t="TTT",
title=c("(IGHV", "wt)" , "mut)"), hr="tr", c="maxstat" )
km(drug = "D_003_4:5", split = "IGHV", t="TTT",
title=c("(IGHV", "wt)" , "mut)"), hr="tr", c="maxstat" )
km(drug = "D_166_4:5", split = "IGHV", t="TTT",
title=c("(IGHV", "wt)" , "mut)"), hr="tr", c="maxstat" )
## ----KM-OS-maxstat, fig.path=plotDir, fig.width = 10, fig.height = 3.3, dev = c("png", "pdf")----
par(mfrow=c(1,3), mar=c(5,5,2,0.9))
km(drug = "D_006_1:5", split = "TP53", t="OS",
title=c("(TP53", "wt)", "mut)"), hr="bl", c="maxstat")
#FIG# 6B right
#FIG# 6C
km(drug = "D_159_1:5", split = "TP53", t="OS", # doxorubicine
title=c("(TP53", "wt)", "mut)"), hr="bl", c="maxstat" )
#FIG# 6B middle
km(drug = "D_010_1:5", split = "TP53", t="OS", # nutlin-3
title=c("(TP53", "wt)", "mut)"), hr="bl", c="maxstat" )
km(drug = "D_002_4:5", split = "IGHV", t="OS",
title=c("(IGHV", "wt)" , "mut)"), hr="bl", c="maxstat" )
km(drug = "D_003_4:5", split = "IGHV", t="OS",
title=c("(IGHV", "wt)" , "mut)"), hr="bl", c="maxstat" )
km(drug = "D_166_4:5", split = "IGHV", t="OS",
title=c("(IGHV", "wt)" , "mut)"), hr="bl", c="maxstat" )
## ----extract---------------------------------------------------------------
extractSome <- function(x) {
sumsu <- summary(x)
data.frame(
`p-value` =
sprintf("%6.3g", sumsu[["coefficients"]][, "Pr(>|z|)"]),
`HR` =
sprintf("%6.3g", signif( sumsu[["coefficients"]][, "exp(coef)"], 2) ),
`lower 95% CI` =
sprintf("%6.3g", signif( sumsu[["conf.int"]][, "lower .95"], 2) ),
`upper 95% CI` =
sprintf("%6.3g", signif( sumsu[["conf.int"]][, "upper .95"], 2),
check.names = FALSE) )
}
## ----covariates, echo=FALSE------------------------------------------------
survT$age <- survT$age/10
survT$IC50beforeTreatment <- ifelse(survT$IC50beforeTreatment==TRUE, 1, 0)
survT$IGHVwt <- ifelse(survT$IGHV==1, 0, 1)
survT$flu15 <- survT$flu15*10
survT$dox15 <- survT$dox15*10
survT$ibr45 <- survT$ibr45*10
survT$ide45 <- survT$ide45*10
survT$prt45 <- survT$prt45*10
## --------------------------------------------------------------------------
surv1 <- coxph(
Surv(T6, died) ~
age +
as.factor(IC50beforeTreatment) +
as.factor(trisomy12) +
as.factor(del11q22.3) +
as.factor(del17p13) +
as.factor(TP53) +
IGHVwt +
flu15, # continuous
#dox15 + # continuous
#flu15:TP53,
#TP53:dox15,
data = survT )
extractSome(surv1)
cat(sprintf("%s patients considerd in the model; number of events %1g\n",
summary(surv1)$n, summary(surv1)[6] ) )
## ----echo=FALSE, results='hide', eval=FALSE--------------------------------
# write(print(xtable(extractSome(surv1))), file=paste0(plotDir,"flu_MD.tex"))
## --------------------------------------------------------------------------
surv2 <- coxph(
Surv(T6, died) ~ #as.factor(survT$TP53) , data=survT )
age +
as.factor(IC50beforeTreatment) +
as.factor(trisomy12) +
as.factor(del11q22.3) +
as.factor(del17p13) +
as.factor(TP53) +
IGHVwt +
#flu15 + # continuous
dox15 , # continuous
#flu15:TP53 ,
#TP53:dox15,
data = survT )
extractSome(surv2)
cat(sprintf("%s patients considerd in the model; number of events %1g\n",
summary(surv2)$n, summary(surv2)[6] ) )
## ----echo=FALSE, results='hide', eval=FALSE--------------------------------
# write(print(xtable(extractSome(surv2))), file=paste0(plotDir,"dox_MD.tex"))
## --------------------------------------------------------------------------
surv4 <- coxph(
Surv(T5, treatedAfter) ~
age +
as.factor(IC50beforeTreatment) +
as.factor(trisomy12) +
as.factor(del11q22.3) +
as.factor(del17p13) +
IGHVwt +
ibr45 +
IGHVwt:ibr45,
data = survT )
extractSome(surv4)
cat(sprintf("%s patients considerd in the model; number of events %1g\n",
summary(surv4)$n, summary(surv4)[6] ) )
## ----echo=FALSE, results='hide', eval=FALSE--------------------------------
# write(print(xtable(extractSome(surv4))), file=paste0(plotDir,"ibr_TTT.tex"))
## --------------------------------------------------------------------------
surv6 <- coxph(
Surv(T5, treatedAfter) ~
age +
as.factor(IC50beforeTreatment) +
as.factor(trisomy12) +
as.factor(del11q22.3) +
as.factor(del17p13) +
IGHVwt +
ide45 +
IGHVwt:ide45,
data = survT )
extractSome(surv6)
cat(sprintf("%s patients considerd in the model; number of events %1g\n",
summary(surv6)$n, summary(surv6)[6] ) )
## ----echo=FALSE, results='hide', eval=FALSE--------------------------------
# write(print(xtable(extractSome(surv6))), file=paste0(plotDir,"ide_TTT.tex"))
## --------------------------------------------------------------------------
surv8 <- coxph(
Surv(T5, treatedAfter) ~
age +
as.factor(IC50beforeTreatment) +
as.factor(trisomy12) +
as.factor(del11q22.3) +
as.factor(del17p13) +
IGHVwt +
prt45 +
IGHVwt:prt45,
data = survT )
extractSome(surv8)
cat(sprintf("%s patients considerd in the model; number of events %1g\n",
summary(surv8)$n, summary(surv8)[6] ) )
## ----echo=FALSE, results='hide', eval=FALSE--------------------------------
# write(print(xtable(extractSome(surv8))), file=paste0(plotDir,"prt_TTT.tex"))
## ---- include=!exists(".standalone"), eval=!exists(".standalone")----------
# sessionInfo()
## ---- message=FALSE, warning=FALSE, include=FALSE--------------------------
rm(list=ls())
## --------------------------------------------------------------------------
sessionInfo()