In MalgorzataOles/BloodCancerMultiOmics2017: "Drug-perturbation-based stratification of blood cancer" by Dietrich S, Oles M, Lu J et al. - experimental data and complete analysis
library("BloodCancerMultiOmics2017")
library("Biobase")
library("RColorBrewer")
library("grid")
library("ggplot2")
library("survival")
library("gtable")
library("forestplot")
library("xtable")
library("maxstat")
plotDir = ifelse(exists(".standalone"), "", "part08/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
options(stringsAsFactors=FALSE)
Survival analysis
Load the data.
data(lpdAll, patmeta, drugs)
Prepare survival data.
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 )
Univariate survival analysis
Define forest functions
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)
}
Combine OS and TTT in one plot.
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)
}
Forest plot for genetic factors
#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)
Forest plot for drug responses
#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)
Kaplan-Meier curves
Genetics factors
#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')
}
#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')
}
Drug responses
Drug responses were dichotomized using maximally selected rank statistics. The analysis is also perforemd within subgroups: TP53 wt/ mut. and IGHV wt/ mut.
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]))
}
}
Time to next treatment (maxstat).
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" )
Overall survival (maxstat).
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" )
Multivariate Cox-model
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) )
}
Define covariates and effects.
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
Chemotherapies
Fludarabine
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] ) )
write(print(xtable(extractSome(surv1))), file=paste0(plotDir,"flu_MD.tex"))
Doxorubicine
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] ) )
write(print(xtable(extractSome(surv2))), file=paste0(plotDir,"dox_MD.tex"))
Targeted therapies
Ibrutinib TTT
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] ) )
write(print(xtable(extractSome(surv4))), file=paste0(plotDir,"ibr_TTT.tex"))
Idelalisib TTT
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] ) )
write(print(xtable(extractSome(surv6))), file=paste0(plotDir,"ide_TTT.tex"))
PRT062607 HCl TTT
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] ) )
write(print(xtable(extractSome(surv8))), file=paste0(plotDir,"prt_TTT.tex"))
sessionInfo()
rm(list=ls())
MalgorzataOles/BloodCancerMultiOmics2017 documentation built on March 29, 2024, 2:29 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library("BloodCancerMultiOmics2017") library("Biobase") library("RColorBrewer") library("grid") library("ggplot2") library("survival") library("gtable") library("forestplot") library("xtable") library("maxstat")
plotDir = ifelse(exists(".standalone"), "", "part08/") if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
options(stringsAsFactors=FALSE)
Survival analysis
Load the data.
data(lpdAll, patmeta, drugs)
Prepare survival data.
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 )
Univariate survival analysis
Define forest functions
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) }
Combine OS and TTT in one plot.
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) }
Forest plot for genetic factors
#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)
Forest plot for drug responses
#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)
Kaplan-Meier curves
Genetics factors
#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') }
#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') }
Drug responses
Drug responses were dichotomized using maximally selected rank statistics. The analysis is also perforemd within subgroups: TP53 wt/ mut. and IGHV wt/ mut.
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])) } }
Time to next treatment (maxstat).
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" )
Overall survival (maxstat).
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" )
Multivariate Cox-model
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) ) }
Define covariates and effects.
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
Chemotherapies
Fludarabine
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] ) )
write(print(xtable(extractSome(surv1))), file=paste0(plotDir,"flu_MD.tex"))
Doxorubicine
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] ) )
write(print(xtable(extractSome(surv2))), file=paste0(plotDir,"dox_MD.tex"))
Targeted therapies
Ibrutinib TTT
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] ) )
write(print(xtable(extractSome(surv4))), file=paste0(plotDir,"ibr_TTT.tex"))
Idelalisib TTT
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] ) )
write(print(xtable(extractSome(surv6))), file=paste0(plotDir,"ide_TTT.tex"))
PRT062607 HCl TTT
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] ) )
write(print(xtable(extractSome(surv8))), file=paste0(plotDir,"prt_TTT.tex"))
sessionInfo()
rm(list=ls())
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.