In lujunyan1118/seahorseCLL: Characterising CLL bioenergetics through Seahorse extracellular flux assay
Association between the clinical phenotype and energy metabolic features
library(seahorseCLL)
library(BloodCancerMultiOmics2017)
library(SummarizedExperiment)
library(grid)
library(survival)
library(gridExtra)
library(maxstat)
library(xtable)
library(DESeq2)
library(cowplot)
library(ggbeeswarm)
library(survminer)
library(tidyverse)
plotDir = ifelse(exists(".standalone"), "", "section05/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)
options(stringsAsFactors=FALSE)
Load datasets
data("patmeta", "seaCombat","lpdAll", "pretreat","doublingTime")
Correlation between seahorse measurements and pretreatment status
Prepare data table for t-test
testTab <- assays(seaCombat)$seaMedian %>% data.frame() %>% rownames_to_column("Measurement") %>%
gather(key = "patID", value = "value", -Measurement) %>%
mutate(pretreated = pretreat[patID,]) %>%
mutate(IGHV = factor(Biobase::exprs(lpdAll)["IGHV Uppsala U/M",patID])) %>%
as.tibble()
Performing t-test for each measurement
tTest <- function(x, y) {
noNA <- !is.na(x) & !is.na(y)
res <- t.test(x[noNA] ~ y[noNA], equal.var = TRUE)
tibble(p = res$p.value,
dm = res$estimate[[2]] - res$estimate[[1]])
}
pTab <- group_by(testTab, Measurement) %>% do(tTest(.$value, .$pretreated)) %>% ungroup() %>%
mutate(p.adj = p.adjust(p, method = "BH"))
Perform ANOVA test, including IGHV as a blocking factor
aovTest <- function(x,y,block) {
noNA <- !is.na(x) & !is.na(y) & !is.na(block)
dataTab <- data.frame(x = x[noNA], y = y[noNA], block = block[noNA])
res <- anova(lm(x ~ block + y))
tibble(p = res$`Pr(>F)`[2])
}
pTab.IGHV <- group_by(testTab, Measurement) %>% do(aovTest(.$value, .$pretreated, .$IGHV)) %>%
ungroup() %>% mutate(p.adj = p.adjust(p, method = "BH"))
A table as output
outTable <- tibble(`Seahorse mearuement` = formatSea(pTab$Measurement),
`p value` = format(pTab$p, digits = 2),
`adjusted p` = format(pTab$p.adj, digits = 3),
`p value (IGHV blocked)` = format(pTab.IGHV$p, digits = 2),
`adjusted p (IGHV blocked)` = format(pTab.IGHV$p.adj, digits = 3))
write(
print(xtable(outTable, digits = 3,
caption = "Student's t-test and ANOVA test (IGHV blocked) results of energy metabolic measurements related to pretreatment status"),
include.rownames=FALSE,
caption.placement = "top")
,file = paste0(plotDir,"tTest_SeahorseVSpretreat.tex"))
Association between IGHV status and pretreatment
IGHVtab <- filter(testTab, !duplicated(patID))
chiRes <- chisq.test(IGHVtab$pretreated, IGHVtab$IGHV)
chiRes
Plot significant associations as beeswarms plot (for supplementary figures)
Plot
seaList <- c("glycolytic.capacity", "glycolytic.reserve")
plotList <- lapply(seaList, function(seaName) {
#prepare table for plotting
plotTab <- filter(testTab, Measurement == seaName) %>%
mutate(Measurement = formatSea(Measurement),
pretreated = ifelse(pretreated,
sprintf("pretreated (n=%s)", sum(pretreated == 1)),
sprintf("untreated (n=%s)", sum(pretreated == 0))),
IGHV = ifelse(is.na(IGHV), "unknown",
ifelse(IGHV == 1, "mutated", "unmutated"))) %>%
mutate(IGHV = factor(IGHV, levels = c("mutated","unmutated","unknown"))) %>%
filter(!is.na(value))
#p value for annotation
pval <- filter(pTab, Measurement == seaName)$p
#color for IGHV status
colorList <- c(mutated = "red", unmutated = "black", unknown = "grey80")
#formatted seaname
seaTitle <- unique(plotTab$Measurement)
#title
plotTitle <- paste(sprintf("'%s (p = '~",seaTitle),
sciPretty(pval, digits = 2),"*')'")
ggplot(plotTab, aes(x=pretreated, y = value)) +
stat_boxplot(geom = "errorbar", width = 0.3) +
geom_boxplot(outlier.shape = NA, col="black", width=0.4) +
geom_beeswarm(cex=2, size =1, aes(col = IGHV)) + theme_classic() +
xlab("") + ylab("ECAR (pMol/min)") + ggtitle(parse(text=plotTitle)) +
scale_color_manual(values = colorList) +
theme(axis.line.x = element_blank(), axis.ticks.x = element_blank(),
axis.title = element_text(size=12, face="bold"),
axis.text = element_text(size=12),
plot.title = element_text(hjust=0.5,size=15),
axis.title.x = element_text(face="bold"))
})
Save the plot
grid.arrange(grobs = plotList, ncol =2)
Correlation between seahorse measurement and lymphocyte doubling time
Correlation test
Pre-processing data
testTab <- assays(seaCombat)$seaMedian %>% data.frame() %>% rownames_to_column("Measurement") %>%
gather(key = "patID", value = "value", -Measurement) %>%
mutate(pretreated = pretreat[patID,],
doubling.time = LDT[patID,]) %>%
mutate(pretreated = factor(pretreated),
IGHV = factor(Biobase::exprs(lpdAll)["IGHV Uppsala U/M",patID])) %>%
filter(!is.na(doubling.time), !is.na(value)) %>%
as.tibble()
Correlation test between seahorse measurement and doubling time
corTest <- function(x, y, block = NULL, method = "pearson") {
if (is.null(block)) {
res <- cor.test(x,y, method = method)
tibble(p = res$p.value, coef = res$estimate[[1]])
} else {
tab <- data.frame(x = x, y = y, block = block)
res <- summary(lm( y ~ x + block, tab)) #how much y can be explained by x and block
tibble(p = res$coefficients[2,4],
coef = cor(x,y))
}
}
corRes <- group_by(testTab, Measurement) %>% do(corTest(.$value, .$doubling.time)) %>%
ungroup() %>% mutate(p.adj = p.adjust(p, method = "BH"))
Correlations between lymphocyte doubling time and IGHV stauts
pRes <- filter(testTab, !duplicated(patID)) %>%
do(data.frame(p = t.test(doubling.time ~ IGHV,.)$p.value))
pRes
Correlation test between seahorse measurement and doubling time
(Considering IGHV as cofactor)
corRes.aov <- group_by(testTab, Measurement) %>% filter(!is.na(IGHV)) %>%
do(corTest(x = .$value, y = .$doubling.time, block = .$IGHV)) %>%
ungroup() %>% mutate(p.adj = p.adjust(p, method = "BH"))
Correlation test between seahorse measurement and doubling time
(within M-CLL)
corRes.M <- group_by(testTab, Measurement) %>% filter(IGHV %in% 1) %>%
do(corTest(x = .$value, y = .$doubling.time)) %>%
ungroup() %>% mutate(p.adj = p.adjust(p, method = "BH"))
Correlation test between seahorse measurement and doubling time
(within U-CLL)
corRes.U <- group_by(testTab, Measurement) %>% filter(IGHV %in% 0) %>%
do(corTest(x = .$value, y = .$doubling.time)) %>%
ungroup() %>% mutate(p.adj = p.adjust(p, method = "BH"))
A table for output
outTable <- tibble(`Seahorse mearuement` = formatSea(corRes$Measurement),
`p value` = format(corRes$p, digits = 2),
`adjusted p` = format(corRes$p.adj, digits = 3),
`p value (IGHV blocked)` = format(corRes.aov$p, digits = 3),
`adjusted p (IGHV blocked)` = format(corRes.aov$p.adj, digits = 3))
write(
print(xtable(outTable, digits = 3,
caption = "Correlation tests between each Seahorse measurements and lymphocyte doubling time"),
include.rownames=FALSE,
caption.placement = "top")
,paste0(plotDir,"tTest_SeahorseVSldt.tex"))
Plot correlations (for supplementary figure)
seaList <- c("glycolysis", "glycolytic.capacity")
plotList.cor <- lapply(seaList, function(seaName) {
#prepare table for plotting
plotTab <- filter(testTab, Measurement == seaName) %>%
mutate(Measurement = formatSea(Measurement),
IGHV = ifelse(is.na(IGHV), "unknown",
ifelse(IGHV == 1, "mutated", "unmutated"))) %>%
mutate(IGHV = factor(IGHV, levels = c("mutated","unmutated","unknown"))) %>%
filter(!is.na(value), !is.na(doubling.time))
#p value for annotation
pval <- filter(corRes, Measurement == seaName)$p
coef <- filter(corRes, Measurement == seaName)$coef
#color for IGHV status
colorList <- c(mutated = "red", unmutated = "black", unknown = "grey80")
#formatted seaname
seaTitle <- unique(plotTab$Measurement)
#prepare correlation test annotations
annoText <- paste("'coefficient ='~",format(coef,digits = 2),"*","', p ='~",sciPretty(pval,digits=2))
ggplot(plotTab, aes(x=value,y=doubling.time)) +
geom_point(size=1, aes(col = IGHV)) +
geom_smooth(method = "lm", se= FALSE) +
xlab("ECAR (pMol/min)") + ylab("Lymphocyte doubling time (days)") +
theme_bw() + ggtitle(seaTitle) +
scale_color_manual(values = colorList) +
annotate("text", x = -Inf, y = Inf, label = annoText,
vjust=1, hjust=0, size = 5, parse = TRUE, col= "darkred") +
theme(panel.grid = element_blank(),
axis.title = element_text(size=13, face="bold"),
axis.text = element_text(size=12),
plot.title = element_text(face="bold", hjust=0.5, size = 15 ),
legend.position = c(0.9,0.1),
axis.title.x = element_text(face="bold"),
plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm"))
})
Show the plot
grid.arrange(grobs = plotList.cor, ncol =2)
Correlation between seahorse measurement and survival
Prepare dataset
seaTable <- assays(seaCombat)$seaMedian
survT = patmeta[colnames(seaTable),]
survT[which(survT[,"IGHV"]=="U") ,"IGHV"] = 0
survT[which(survT[,"IGHV"]=="M") ,"IGHV"] = 1
survT$IGHV = as.numeric(survT$IGHV)
survT$pretreatment <- pretreat[rownames(survT),]
colnames(survT) = gsub("Age4Main", "age", colnames(survT))
survT$treatment <- survT$treatedAfter | survT$pretreatment
#add seahorse measurement information
survT <- cbind(survT, t(seaTable))
Univariate survival analysis
Caclulation correlations
Function to calculate correlations
com <- function( Time, endpoint, scaleX, sub, d, split, drug_names) {
res <- lapply(d, function(g) {
#drug <- survT[,g] * scaleX
drug <- survT[,g]
#drug <- (drug - mean(drug, na.rm = TRUE))/sd(drug, na.rm=TRUE)
## 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
}
All samples
d <- rownames(seaTable)
sea_names <- formatSea(d)
ttt <- com(Time="T5", endpoint="treatedAfter", sub=99, d=d,
split="", drug_names=sea_names, scaleX=1)
os <- com(Time="T6", endpoint="died", sub=99, d=d, split="",
drug_names=sea_names, scaleX=1)
#TTT result
ttt
#OS result
os
M-CLL samples
d <- rownames(seaTable)
sea_names <- formatSea(d)
ttt.M <- com(Time="T5", endpoint="treatedAfter", sub=1, d=d,
split="IGHV", drug_names=sea_names, scaleX=1)
os.M <- com(Time="T6", endpoint="died", sub=1, d=d,
split="IGHV",drug_names=sea_names, scaleX=1)
#TTT result
ttt.M
#OS result
os.M
U-CLL samples
d <- rownames(seaTable)
sea_names <- formatSea(d)
ttt.U <- com(Time="T5", endpoint="treatedAfter", sub=0, d=d,
split="IGHV", drug_names=sea_names, scaleX=1)
os.U <- com(Time="T6", endpoint="died", sub=0, d=d,
split="IGHV",drug_names=sea_names, scaleX=1)
#TTT result
ttt.U
#OS result
os.U
KM plot
Function for KM plot
ggkm <- function(response, time, endpoint, titlePlot = "KM plot", pval = NULL, stat = "median") {
#function for km plot
survS <- data.frame(time = time,
endpoint = endpoint)
if (stat == "maxstat") {
ms <- maxstat.test(Surv(time, endpoint) ~ response,
data = survS,
smethod = "LogRank",
minprop = 0.2,
maxprop = 0.8,
alpha = NULL)
survS$group <- factor(ifelse(response >= ms$estimate, "high", "low"))
} else if (stat == "median") {
med <- median(response, na.rm = TRUE)
survS$group <- factor(ifelse(response >= med, "high", "low"))
} else if (stat == "binary") {
survS$group <- factor(response)
}
if (is.null(pval)) pval = TRUE
p <- ggsurvplot(survfit(Surv(time, endpoint) ~ group, data = survS),
data = survS, pval = TRUE, conf.int = TRUE,
ylab = "Fraction", xlab = "Time (years)", title = titlePlot,
ggtheme = theme_bw() + theme(plot.title = element_text(hjust =0.5)))$plot
p
}
#a fucntion to assign unit, either ECAR or OCR or none
giveUnit <- function(x) {
ocrList <- c("basal.respiration","ATP.production","proton.leak","maximal.respiration",
"spare.respiratory.capacity","OCR")
ecarList <- c("glycolysis","glycolytic.capacity","glycolytic.reserve","ECAR")
if (x %in% ocrList)
return("OCR")
else if (x %in% ecarList)
return("ECAR")
else return("")
}
Stratrified by metabolism
TTT
tttList <- lapply(rownames(seaCombat), function(n) {
survS <- data.frame(time = survT$T5,
endpoint = survT$treatedAfter)
val <- survT[[n]]
ms <- maxstat.test(Surv(time, endpoint) ~ val,
data = survS,
smethod = "LogRank",
minprop = 0.2,
maxprop = 0.8,
alpha = NULL)
plotTab <- survS %>% mutate(response = val) %>%
filter(!is.na(response), !is.na(endpoint),!is.na(time)) %>%
mutate(group = ifelse(response >= ms$estimate,
sprintf("high %s (%s >= %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate),
sprintf("low %s (%s < %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate)))
nTab <- table(plotTab$group)
plotTab <- mutate(plotTab, group = paste0(group,nTab[group],")"))
pval <- sprintf("p = %1.3f",ttt[formatSea(n),"p"])
p <- ggsurvplot(survfit(Surv(time, endpoint) ~ group, data = plotTab),
data = plotTab, pval = pval, conf.int = FALSE,
legend = c(0.6, 0.12),pval.coord = c(0.1,0.25),
legend.labs = sort(unique(plotTab$group)),
legend.title = "group", palette = "Dark2",
ylab = "Fraction treatment free", xlab = "Time (years)", title = formatSea(n),
ggtheme = theme_classic() + theme(axis.title = element_text(size=13, face="bold"),
axis.text = element_text(size=12),
plot.title = element_text(face="bold", hjust=0.5, size = 15 ),
axis.title.x = element_text(face="bold"),
plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")))$plot
})
names(tttList) <- names(seaCombat)
OS
osList <- lapply(rownames(seaCombat), function(n) {
survS <- data.frame(time = survT$T6,
endpoint = survT$died)
val <- survT[[n]]
ms <- maxstat.test(Surv(time, endpoint) ~ val,
data = survS,
smethod = "LogRank",
minprop = 0.2,
maxprop = 0.8,
alpha = NULL)
plotTab <- survS %>% mutate(response = val) %>%
filter(!is.na(response), !is.na(endpoint),!is.na(time)) %>%
mutate(group = ifelse(response >= ms$estimate,
sprintf("high %s (%s >= %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate),
sprintf("low %s (%s < %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate)))
nTab <- table(plotTab$group)
plotTab <- mutate(plotTab, group = paste0(group,nTab[group],")"))
pval <- sprintf("p = %1.3f",os[formatSea(n),"p"])
p <- ggsurvplot(survfit(Surv(time, endpoint) ~ group, data = plotTab),
data = plotTab, pval = pval, conf.int = FALSE,
legend = c(0.6,0.12),pval.coord = c(0.1,0.25),
legend.labs = sort(unique(plotTab$group)),
legend.title = "group", palette = "Dark2",
ylab = "Fraction overall survival", xlab = "Time (years)", title = formatSea(n),
ggtheme = theme_classic() + theme(axis.title = element_text(size=13, face="bold"),
axis.text = element_text(size=12),
plot.title = element_text(face="bold", hjust=0.5, size = 15 ),
axis.title.x = element_text(face="bold"),
plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")))$plot
})
names(osList) <- names(seaCombat)
KM plot for significant assocations in single variante cox model
p<-plot_grid(tttList$glycolytic.reserve, tttList$maximal.respiration,
tttList$spare.respiratory.capacity,
osList$glycolytic.capacity, osList$glycolytic.reserve,
labels = NULL)
p
Straitified by both metabolism and IGHV status
TTT
tttKm <- lapply(rownames(seaCombat), function(n) {
survS <- data.frame(time = survT$T5,
endpoint = survT$treatedAfter)
val <- survT[[n]]
ms <- maxstat.test(Surv(time, endpoint) ~ val,
data = survS,
smethod = "LogRank",
minprop = 0.2,
maxprop = 0.8,
alpha = NULL)
plotTab <- survS %>% mutate(response = val, IGHV = survT$IGHV) %>%
filter(!is.na(response),!is.na(IGHV), !is.na(endpoint),!is.na(time)) %>%
mutate(group = ifelse(response >= ms$estimate,
sprintf("high %s (%s >= %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate),
sprintf("low %s (%s < %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate))) %>%
mutate(group = ifelse(IGHV==1, paste0("M-CLL with ",group), paste0("U-CLL with ",group)))
nTab <- table(plotTab$group)
plotTab <- mutate(plotTab, group = paste0(group,nTab[group],")"))
p <- ggsurvplot(survfit(Surv(time, endpoint) ~ group, data = plotTab),
data = plotTab, pval = TRUE, conf.int = FALSE,
legend = c(0.6,0.2),legend.labs = sort(unique(plotTab$group)),
pval.coord = c(0.1,0.4),
legend.title = "group", palette = "Dark2",
ylab = "Fraction treatment free", xlab = "Time (years)", title = formatSea(n),
ggtheme = theme_classic() + theme(axis.title = element_text(size=13, face="bold"),
axis.text = element_text(size=12),
plot.title = element_text(face="bold", hjust=0.5, size = 15 ),
axis.title.x = element_text(face="bold"),
plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")))$plot
})
names(tttKm) <- names(seaCombat)
OS
osKm <- lapply(rownames(seaCombat), function(n) {
survS <- data.frame(time = survT$T6,
endpoint = survT$died)
val <- survT[[n]]
ms <- maxstat.test(Surv(time, endpoint) ~ val,
data = survS,
smethod = "LogRank",
minprop = 0.2,
maxprop = 0.8,
alpha = NULL)
plotTab <- survS %>% mutate(response = val, IGHV = survT$IGHV) %>%
filter(!is.na(response),!is.na(IGHV), !is.na(endpoint),!is.na(time)) %>%
mutate(group = ifelse(response >= ms$estimate,
sprintf("high %s (%s >= %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate),
sprintf("low %s (%s < %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate))) %>%
mutate(group = ifelse(IGHV==1, paste0("M-CLL with ",group), paste0("U-CLL with ",group)))
nTab <- table(plotTab$group)
plotTab <- mutate(plotTab, group = paste0(group,nTab[group],")"))
p <- ggsurvplot(survfit(Surv(time, endpoint) ~ group, data = plotTab),
data = plotTab, pval = TRUE, conf.int = FALSE,
legend = c(0.6,0.2),legend.labs = sort(unique(plotTab$group)),
pval.coord = c(0.1,0.4),
legend.title = "group", palette = "Dark2",
ylab = "Fraction overall survival", xlab = "Time (years)", title = formatSea(n),
ggtheme = theme_classic() + theme(axis.title = element_text(size=13, face="bold"),
axis.text = element_text(size=12),
plot.title = element_text(face="bold", hjust=0.5, size = 15 ),
axis.title.x = element_text(face="bold"),
plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")))$plot
})
names(osKm) <- names(seaCombat)
Multi-variate cox model
Prepare dataset
#add genetic variants to survival table
survT$SF3B1 <- Biobase::exprs(lpdAll)[ "SF3B1", rownames(survT) ]
survT$NOTCH1 <- Biobase::exprs(lpdAll)[ "NOTCH1", rownames(survT) ]
survT$BRAF <- Biobase::exprs(lpdAll)[ "BRAF", rownames(survT) ]
survT$TP53 <- Biobase::exprs(lpdAll)[ "TP53", rownames(survT) ]
survT$del17p13 <- Biobase::exprs(lpdAll)[ "del17p13", rownames(survT) ]
survT$del11q22.3 <- Biobase::exprs(lpdAll)[ "del11q22.3", rownames(survT) ]
survT$trisomy12 <- Biobase::exprs(lpdAll)[ "trisomy12", rownames(survT) ]
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$IGHVwt <- ifelse(survT$IGHV==1, 0, 1)
survT$treatment <- ifelse(survT$treatment, 1, 0)
TTT
glycolytic reserve
surv1 <- coxph(
Surv(T5, treatedAfter) ~
age +
as.factor(trisomy12) +
as.factor(del11q22.3) +
as.factor(del17p13) +
as.factor(TP53) +
as.factor(IGHVwt) +
glycolytic.reserve, # continuous
data = survT )
colFactor <- data.frame(factor = c("age",
"trisomy12", "del11q22.3",
"del17p13","TP53","U-CLL",
"glycolytic reserve"))
outTab <- cbind(colFactor,extractSome(surv1))
write(print(xtable(outTab,
caption = "Multivariate Cox regression model for time to treatment with glycolytic reserve as a covariate"),
include.rownames=FALSE,
caption.placement = "top"), file = paste0(plotDir,"glyRes_TTT.tex"))
cat(sprintf("%s patients considerd in the model; number of events %1g\n",
summary(surv1)$n, summary(surv1)[6] ) )
maximal respiration
surv1 <- coxph(
Surv(T5, treatedAfter) ~
age +
as.factor(trisomy12) +
as.factor(del11q22.3) +
as.factor(del17p13) +
as.factor(TP53) +
IGHVwt +
maximal.respiration, # continuous
data = survT )
colFactor <- data.frame(factor = c("age",
"trisomy12", "del11q22.3",
"del17p13","TP53","U-CLL",
"maximal respiration"))
outTab <- cbind(colFactor,extractSome(surv1))
write(print(xtable(outTab,
caption = "Multivariate Cox regression model for time to treatment with maximal respiration as a covariate"),
include.rownames=FALSE,
caption.placement = "top"), file = paste0(plotDir,"maxRes_TTT.tex"))
cat(sprintf("%s patients considerd in the model; number of events %1g\n",
summary(surv1)$n, summary(surv1)[6] ) )
spare respiratory capacity
surv1 <- coxph(
Surv(T5, treatedAfter) ~
age +
as.factor(trisomy12) +
as.factor(del11q22.3) +
as.factor(del17p13) +
as.factor(TP53) +
IGHVwt +
spare.respiratory.capacity, # continuous
data = survT )
colFactor <- data.frame(factor = c("age",
"trisomy12", "del11q22.3",
"del17p13","TP53","U-CLL",
"spare.respiratory.capacity"))
outTab <- cbind(colFactor,extractSome(surv1))
write(print(xtable(outTab,
caption = "Multivariate Cox regression model for time to treatment with glycolytic reserve as a covariate"),
include.rownames=FALSE,
caption.placement = "top"), file = paste0(plotDir,"spResCap_TTT.tex"))
cat(sprintf("%s patients considerd in the model; number of events %1g\n",
summary(surv1)$n, summary(surv1)[6] ) )
OS
glycolytic reserve
surv1 <- coxph(
Surv(T6, died) ~
age + as.factor(treatment) +
as.factor(trisomy12) +
as.factor(del11q22.3) +
as.factor(del17p13) +
as.factor(TP53) +
IGHVwt +
glycolytic.reserve, # continuous
data = survT )
colFactor <- data.frame(factor = c("age", "treatment",
"trisomy12", "del11q22.3",
"del17p13","TP53","U-CLL",
"glycolytic reserve"))
outTab <- cbind(colFactor,extractSome(surv1))
write(print(xtable(outTab,
caption = "Multivariate Cox regression model for overall survival with glycolytic reserve as a covariate"),
include.rownames=FALSE,
caption.placement = "top"), file = paste0(plotDir,"glyRes_OS.tex"))
cat(sprintf("%s patients considerd in the model; number of events %1g\n",
summary(surv1)$n, summary(surv1)[6] ) )
write.csv2(outTab, "glycolytic.reserve_VS_OS.csv")
glycolytic capacity
surv1 <- coxph(
Surv(T6, died) ~
age + as.factor(treatment) +
as.factor(trisomy12) +
as.factor(del11q22.3) +
as.factor(del17p13) +
as.factor(TP53) +
IGHVwt +
glycolytic.capacity, # continuous
data = survT )
colFactor <- data.frame(factor = c("age", "treatment",
"trisomy12", "del11q22.3",
"del17p13","TP53","U-CLL",
"glycolytic capacity"))
outTab <- cbind(colFactor,extractSome(surv1))
write(print(xtable(outTab,
caption = "Multivariate Cox regression model for overall survival with glycolytic capacity as a covariate"),
include.rownames=FALSE,
caption.placement = "top"), file = paste0(plotDir,"glyCap_OS.tex"))
cat(sprintf("%s patients considerd in the model; number of events %1g\n",
summary(surv1)$n, summary(surv1)[6] ))
write.csv2(outTab, "glycolytic.capacity_VS_OS.csv")
Association between bioenergetic features and CD38/CD49d expression
Prepare data
Query CD38 and CD49d (ITGA4) expressions from RNAseq data
data("dds")
geneList <- c("CD38","ITGA4")
ddsSea <- dds[,dds$PatID %in% colnames(seaOri)]
ddsSea <- estimateSizeFactors(ddsSea)
ddsSea <- varianceStabilizingTransformation(ddsSea)
ddsSea <- ddsSea[rowData(ddsSea)$symbol %in% geneList,]
countSea <- assay(ddsSea)
rownames(countSea) <- rowData(ddsSea)$symbol
colnames(countSea) <- ddsSea$PatID
countSea <- data.frame(t(countSea)) %>% rownames_to_column("patID")
Combine the phenotype table and subset for patients with Seahorse measurement
phenoTab <- countSea[match(colnames(seaCombat), countSea$patID),] %>%
filter(!is.na(patID)) %>%
mutate(IGHV = Biobase::exprs(lpdAll)["IGHV Uppsala U/M",patID]) %>%
mutate(IGHV = ifelse(is.na(IGHV),NA, ifelse(IGHV ==1, "M","U")))
Correlation test
seaTab <- assay(seaCombat)[,phenoTab$patID]
stopifnot(all(colnames(seaTab) == phenoTab$patID))
corRes <- lapply(rownames(seaTab), function(seaName) {
lapply(colnames(phenoTab)[2:3], function(phenoName){
testTab <- data.frame(sea = seaTab[seaName,], pheno = phenoTab[[phenoName]], IGHV = as.factor(phenoTab$IGHV))
lmRes <- summary(lm(pheno ~ sea, data = testTab, na.action = na.omit))
lmRes.IGHV <- summary(lm(pheno ~ sea + IGHV, data = testTab, na.action = na.omit))
data.frame(measurement = seaName, phenotype = phenoName,
p = lmRes$coefficients[2,4], p.IGHV = lmRes.IGHV$coefficients[2,4],
stringsAsFactors = FALSE)
}) %>% bind_rows()
}) %>% bind_rows() %>% arrange(p) %>% mutate(p.adj = p.adjust(p, method = "BH"),
p.IGHV.adj = p.adjust(p.IGHV, method = "BH"))
Save a table for supplementary table
tabOut <- filter(corRes, p.adj <= 0.05) %>% mutate(measurement = formatSea(measurement)) %>%
dplyr::rename(Measurement = measurement, Gene = phenotype,
`p value` = p, `p value (IGHV blocked)` = p.IGHV,
`adjusted p value` = p.adj, `adjusted p value (IGHV blocked)`= p.IGHV.adj)
write(print(xtable(tabOut, digits = 3,
caption = "Correlation test results between energy metabolic measurements and CD38/IGTA4(CD49d) expression"),
include.rownames=FALSE,
caption.placement = "top"), file = paste0(plotDir,"corTest_seaVSCD38.tex"))
Plot significant correlations
corRes.sig <- filter(corRes, p.IGHV.adj < 0.05)
pList.all <- lapply(seq(nrow(corRes.sig)), function(i) {
seaName <- corRes.sig$measurement[i]
phenoName <- corRes.sig$phenotype[i]
plotTab <- data.frame(sea = seaTab[seaName,],
pheno = phenoTab[[phenoName]],
IGHV = phenoTab$IGHV) %>%
mutate(IGHV = ifelse(is.na(IGHV),"unknown",ifelse(
IGHV == "M", "mutated","unmutated")))
corRes <- cor.test(plotTab$sea, plotTab$pheno)
annoText <- paste("'coefficient ='~",format(corRes$estimate,digits = 2),"*","', p ='~",sciPretty(corRes$p.value,digits=2))
ggplot(plotTab, aes(x=sea, y = pheno)) +
geom_point(aes(color = IGHV)) + geom_smooth(method = "lm", se=FALSE ) +
theme_bw() + ylab(paste0(phenoName," (normalized expression)")) + xlab("ECAR (pMol/min)") +
scale_color_manual(values = c(mutated = "red", unmutated = "black", unknown = "grey80")) + ggtitle(formatSea(seaName)) +
annotate("text", x = -Inf, y = Inf, label = annoText,
vjust=1, hjust=0, size = 5, parse = TRUE, col= "darkred") +
theme(panel.grid = element_blank(),
axis.title = element_text(size=13, face="bold"),
axis.text = element_text(size=12),
plot.title = element_text(face="bold", hjust=0.5, size = 15 ),
legend.position = c(0.85,0.15),
legend.background = element_rect(fill = NA),
axis.title.x = element_text(face="bold"),
plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm"))
})
grid.arrange(grobs = pList.all, ncol=2)
pList.U <- lapply(seq(nrow(corRes.sig)), function(i) {
seaName <- corRes.sig$measurement[i]
phenoName <- corRes.sig$phenotype[i]
plotTab <- data.frame(sea = seaTab[seaName,],
pheno = phenoTab[[phenoName]],
IGHV = phenoTab$IGHV) %>%
filter(IGHV == "U")
corRes <- cor.test(plotTab$sea, plotTab$pheno)
annoText <- paste("'coefficient ='~",format(corRes$estimate,digits = 2),"*","', p ='~",sciPretty(corRes$p.value,digits=2))
ggplot(plotTab, aes(x=sea, y = pheno)) +
geom_point(color = "black") + geom_smooth(method = "lm", se=FALSE ) +
theme_bw() + ylab(paste0(phenoName," (normalized expression)")) + xlab("ECAR (pMol/min)") +
scale_color_manual(values = c(mutated = "red", unmutated = "black", unknown = "grey80")) + ggtitle(paste0(formatSea(seaName)," (U-CLL samples only)")) +
annotate("text", x = -Inf, y = Inf, label = annoText,
vjust=1, hjust=0, size = 5, parse = TRUE, col= "darkred") +
theme(panel.grid = element_blank(),
axis.title = element_text(size=13, face="bold"),
axis.text = element_text(size=12),
plot.title = element_text(face="bold", hjust=0.5, size = 15 ),
legend.position = c(0.9,0.1),
axis.title.x = element_text(face="bold"),
plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm"))
})
pList.M <- lapply(seq(nrow(corRes.sig)), function(i) {
seaName <- corRes.sig$measurement[i]
phenoName <- corRes.sig$phenotype[i]
plotTab <- data.frame(sea = seaTab[seaName,],
pheno = phenoTab[[phenoName]],
IGHV = phenoTab$IGHV) %>%
filter(IGHV == "M")
corRes <- cor.test(plotTab$sea, plotTab$pheno)
annoText <- paste("'coefficient ='~",format(corRes$estimate,digits = 2),"*","', p ='~",sciPretty(corRes$p.value,digits=2))
ggplot(plotTab, aes(x=sea, y = pheno)) +
geom_point(color = "red") + geom_smooth(method = "lm", se=FALSE ) +
theme_bw() + ylab(paste0(phenoName," (normalized expression)")) + xlab("ECAR (pMol/min)") +
scale_color_manual(values = c(mutated = "red", unmutated = "black", unknown = "grey80")) + ggtitle(paste0(formatSea(seaName)," (M-CLL samples only)")) +
annotate("text", x = -Inf, y = Inf, label = annoText,
vjust=1, hjust=0, size = 5, parse = TRUE, col= "darkred") +
theme(panel.grid = element_blank(),
axis.title = element_text(size=13, face="bold"),
axis.text = element_text(size=12),
plot.title = element_text(face="bold", hjust=0.5, size = 15 ),
legend.position = c(0.9,0.1),
axis.title.x = element_text(face="bold"),
plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm"))
})
grid.arrange(grobs = c(pList.U,pList.M), ncol=2)
Organize a figure for the clinical part
figOut <- c(osKm$glycolytic.capacity, osKm$glycolytic.reserve, pList.all[[1]],pList.all[[2]])
title = ggdraw() + draw_figure_label("Figure 5", fontface = "bold", position = "top.left",size=20)
p<-plot_grid(osKm$glycolytic.capacity, osKm$glycolytic.reserve,
pList.all[[1]],
pList.all[[2]],
labels = c("A","B","C","D"), label_size = 22)
plot_grid(title, p, rel_heights = c(0.05,0.95), ncol = 1)
lujunyan1118/seahorseCLL documentation built on May 12, 2019, 6:16 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.
Association between the clinical phenotype and energy metabolic features
library(seahorseCLL) library(BloodCancerMultiOmics2017) library(SummarizedExperiment) library(grid) library(survival) library(gridExtra) library(maxstat) library(xtable) library(DESeq2) library(cowplot) library(ggbeeswarm) library(survminer) library(tidyverse)
plotDir = ifelse(exists(".standalone"), "", "section05/") if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir) options(stringsAsFactors=FALSE)
Load datasets
data("patmeta", "seaCombat","lpdAll", "pretreat","doublingTime")
Correlation between seahorse measurements and pretreatment status
Prepare data table for t-test
testTab <- assays(seaCombat)$seaMedian %>% data.frame() %>% rownames_to_column("Measurement") %>% gather(key = "patID", value = "value", -Measurement) %>% mutate(pretreated = pretreat[patID,]) %>% mutate(IGHV = factor(Biobase::exprs(lpdAll)["IGHV Uppsala U/M",patID])) %>% as.tibble()
Performing t-test for each measurement
tTest <- function(x, y) { noNA <- !is.na(x) & !is.na(y) res <- t.test(x[noNA] ~ y[noNA], equal.var = TRUE) tibble(p = res$p.value, dm = res$estimate[[2]] - res$estimate[[1]]) } pTab <- group_by(testTab, Measurement) %>% do(tTest(.$value, .$pretreated)) %>% ungroup() %>% mutate(p.adj = p.adjust(p, method = "BH"))
Perform ANOVA test, including IGHV as a blocking factor
aovTest <- function(x,y,block) { noNA <- !is.na(x) & !is.na(y) & !is.na(block) dataTab <- data.frame(x = x[noNA], y = y[noNA], block = block[noNA]) res <- anova(lm(x ~ block + y)) tibble(p = res$`Pr(>F)`[2]) } pTab.IGHV <- group_by(testTab, Measurement) %>% do(aovTest(.$value, .$pretreated, .$IGHV)) %>% ungroup() %>% mutate(p.adj = p.adjust(p, method = "BH"))
A table as output
outTable <- tibble(`Seahorse mearuement` = formatSea(pTab$Measurement), `p value` = format(pTab$p, digits = 2), `adjusted p` = format(pTab$p.adj, digits = 3), `p value (IGHV blocked)` = format(pTab.IGHV$p, digits = 2), `adjusted p (IGHV blocked)` = format(pTab.IGHV$p.adj, digits = 3)) write( print(xtable(outTable, digits = 3, caption = "Student's t-test and ANOVA test (IGHV blocked) results of energy metabolic measurements related to pretreatment status"), include.rownames=FALSE, caption.placement = "top") ,file = paste0(plotDir,"tTest_SeahorseVSpretreat.tex"))
Association between IGHV status and pretreatment
IGHVtab <- filter(testTab, !duplicated(patID)) chiRes <- chisq.test(IGHVtab$pretreated, IGHVtab$IGHV) chiRes
Plot significant associations as beeswarms plot (for supplementary figures)
Plot
seaList <- c("glycolytic.capacity", "glycolytic.reserve") plotList <- lapply(seaList, function(seaName) { #prepare table for plotting plotTab <- filter(testTab, Measurement == seaName) %>% mutate(Measurement = formatSea(Measurement), pretreated = ifelse(pretreated, sprintf("pretreated (n=%s)", sum(pretreated == 1)), sprintf("untreated (n=%s)", sum(pretreated == 0))), IGHV = ifelse(is.na(IGHV), "unknown", ifelse(IGHV == 1, "mutated", "unmutated"))) %>% mutate(IGHV = factor(IGHV, levels = c("mutated","unmutated","unknown"))) %>% filter(!is.na(value)) #p value for annotation pval <- filter(pTab, Measurement == seaName)$p #color for IGHV status colorList <- c(mutated = "red", unmutated = "black", unknown = "grey80") #formatted seaname seaTitle <- unique(plotTab$Measurement) #title plotTitle <- paste(sprintf("'%s (p = '~",seaTitle), sciPretty(pval, digits = 2),"*')'") ggplot(plotTab, aes(x=pretreated, y = value)) + stat_boxplot(geom = "errorbar", width = 0.3) + geom_boxplot(outlier.shape = NA, col="black", width=0.4) + geom_beeswarm(cex=2, size =1, aes(col = IGHV)) + theme_classic() + xlab("") + ylab("ECAR (pMol/min)") + ggtitle(parse(text=plotTitle)) + scale_color_manual(values = colorList) + theme(axis.line.x = element_blank(), axis.ticks.x = element_blank(), axis.title = element_text(size=12, face="bold"), axis.text = element_text(size=12), plot.title = element_text(hjust=0.5,size=15), axis.title.x = element_text(face="bold")) })
Save the plot
grid.arrange(grobs = plotList, ncol =2)
Correlation between seahorse measurement and lymphocyte doubling time
Correlation test
Pre-processing data
testTab <- assays(seaCombat)$seaMedian %>% data.frame() %>% rownames_to_column("Measurement") %>% gather(key = "patID", value = "value", -Measurement) %>% mutate(pretreated = pretreat[patID,], doubling.time = LDT[patID,]) %>% mutate(pretreated = factor(pretreated), IGHV = factor(Biobase::exprs(lpdAll)["IGHV Uppsala U/M",patID])) %>% filter(!is.na(doubling.time), !is.na(value)) %>% as.tibble()
Correlation test between seahorse measurement and doubling time
corTest <- function(x, y, block = NULL, method = "pearson") { if (is.null(block)) { res <- cor.test(x,y, method = method) tibble(p = res$p.value, coef = res$estimate[[1]]) } else { tab <- data.frame(x = x, y = y, block = block) res <- summary(lm( y ~ x + block, tab)) #how much y can be explained by x and block tibble(p = res$coefficients[2,4], coef = cor(x,y)) } } corRes <- group_by(testTab, Measurement) %>% do(corTest(.$value, .$doubling.time)) %>% ungroup() %>% mutate(p.adj = p.adjust(p, method = "BH"))
Correlations between lymphocyte doubling time and IGHV stauts
pRes <- filter(testTab, !duplicated(patID)) %>% do(data.frame(p = t.test(doubling.time ~ IGHV,.)$p.value)) pRes
Correlation test between seahorse measurement and doubling time (Considering IGHV as cofactor)
corRes.aov <- group_by(testTab, Measurement) %>% filter(!is.na(IGHV)) %>% do(corTest(x = .$value, y = .$doubling.time, block = .$IGHV)) %>% ungroup() %>% mutate(p.adj = p.adjust(p, method = "BH"))
Correlation test between seahorse measurement and doubling time (within M-CLL)
corRes.M <- group_by(testTab, Measurement) %>% filter(IGHV %in% 1) %>% do(corTest(x = .$value, y = .$doubling.time)) %>% ungroup() %>% mutate(p.adj = p.adjust(p, method = "BH"))
Correlation test between seahorse measurement and doubling time (within U-CLL)
corRes.U <- group_by(testTab, Measurement) %>% filter(IGHV %in% 0) %>% do(corTest(x = .$value, y = .$doubling.time)) %>% ungroup() %>% mutate(p.adj = p.adjust(p, method = "BH"))
A table for output
outTable <- tibble(`Seahorse mearuement` = formatSea(corRes$Measurement), `p value` = format(corRes$p, digits = 2), `adjusted p` = format(corRes$p.adj, digits = 3), `p value (IGHV blocked)` = format(corRes.aov$p, digits = 3), `adjusted p (IGHV blocked)` = format(corRes.aov$p.adj, digits = 3)) write( print(xtable(outTable, digits = 3, caption = "Correlation tests between each Seahorse measurements and lymphocyte doubling time"), include.rownames=FALSE, caption.placement = "top") ,paste0(plotDir,"tTest_SeahorseVSldt.tex"))
Plot correlations (for supplementary figure)
seaList <- c("glycolysis", "glycolytic.capacity") plotList.cor <- lapply(seaList, function(seaName) { #prepare table for plotting plotTab <- filter(testTab, Measurement == seaName) %>% mutate(Measurement = formatSea(Measurement), IGHV = ifelse(is.na(IGHV), "unknown", ifelse(IGHV == 1, "mutated", "unmutated"))) %>% mutate(IGHV = factor(IGHV, levels = c("mutated","unmutated","unknown"))) %>% filter(!is.na(value), !is.na(doubling.time)) #p value for annotation pval <- filter(corRes, Measurement == seaName)$p coef <- filter(corRes, Measurement == seaName)$coef #color for IGHV status colorList <- c(mutated = "red", unmutated = "black", unknown = "grey80") #formatted seaname seaTitle <- unique(plotTab$Measurement) #prepare correlation test annotations annoText <- paste("'coefficient ='~",format(coef,digits = 2),"*","', p ='~",sciPretty(pval,digits=2)) ggplot(plotTab, aes(x=value,y=doubling.time)) + geom_point(size=1, aes(col = IGHV)) + geom_smooth(method = "lm", se= FALSE) + xlab("ECAR (pMol/min)") + ylab("Lymphocyte doubling time (days)") + theme_bw() + ggtitle(seaTitle) + scale_color_manual(values = colorList) + annotate("text", x = -Inf, y = Inf, label = annoText, vjust=1, hjust=0, size = 5, parse = TRUE, col= "darkred") + theme(panel.grid = element_blank(), axis.title = element_text(size=13, face="bold"), axis.text = element_text(size=12), plot.title = element_text(face="bold", hjust=0.5, size = 15 ), legend.position = c(0.9,0.1), axis.title.x = element_text(face="bold"), plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) })
Show the plot
grid.arrange(grobs = plotList.cor, ncol =2)
Correlation between seahorse measurement and survival
Prepare dataset
seaTable <- assays(seaCombat)$seaMedian survT = patmeta[colnames(seaTable),] survT[which(survT[,"IGHV"]=="U") ,"IGHV"] = 0 survT[which(survT[,"IGHV"]=="M") ,"IGHV"] = 1 survT$IGHV = as.numeric(survT$IGHV) survT$pretreatment <- pretreat[rownames(survT),] colnames(survT) = gsub("Age4Main", "age", colnames(survT)) survT$treatment <- survT$treatedAfter | survT$pretreatment #add seahorse measurement information survT <- cbind(survT, t(seaTable))
Univariate survival analysis
Caclulation correlations
Function to calculate correlations
com <- function( Time, endpoint, scaleX, sub, d, split, drug_names) { res <- lapply(d, function(g) { #drug <- survT[,g] * scaleX drug <- survT[,g] #drug <- (drug - mean(drug, na.rm = TRUE))/sd(drug, na.rm=TRUE) ## 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 }
All samples
d <- rownames(seaTable) sea_names <- formatSea(d) ttt <- com(Time="T5", endpoint="treatedAfter", sub=99, d=d, split="", drug_names=sea_names, scaleX=1) os <- com(Time="T6", endpoint="died", sub=99, d=d, split="", drug_names=sea_names, scaleX=1) #TTT result ttt #OS result os
M-CLL samples
d <- rownames(seaTable) sea_names <- formatSea(d) ttt.M <- com(Time="T5", endpoint="treatedAfter", sub=1, d=d, split="IGHV", drug_names=sea_names, scaleX=1) os.M <- com(Time="T6", endpoint="died", sub=1, d=d, split="IGHV",drug_names=sea_names, scaleX=1) #TTT result ttt.M #OS result os.M
U-CLL samples
d <- rownames(seaTable) sea_names <- formatSea(d) ttt.U <- com(Time="T5", endpoint="treatedAfter", sub=0, d=d, split="IGHV", drug_names=sea_names, scaleX=1) os.U <- com(Time="T6", endpoint="died", sub=0, d=d, split="IGHV",drug_names=sea_names, scaleX=1) #TTT result ttt.U #OS result os.U
KM plot
Function for KM plot
ggkm <- function(response, time, endpoint, titlePlot = "KM plot", pval = NULL, stat = "median") { #function for km plot survS <- data.frame(time = time, endpoint = endpoint) if (stat == "maxstat") { ms <- maxstat.test(Surv(time, endpoint) ~ response, data = survS, smethod = "LogRank", minprop = 0.2, maxprop = 0.8, alpha = NULL) survS$group <- factor(ifelse(response >= ms$estimate, "high", "low")) } else if (stat == "median") { med <- median(response, na.rm = TRUE) survS$group <- factor(ifelse(response >= med, "high", "low")) } else if (stat == "binary") { survS$group <- factor(response) } if (is.null(pval)) pval = TRUE p <- ggsurvplot(survfit(Surv(time, endpoint) ~ group, data = survS), data = survS, pval = TRUE, conf.int = TRUE, ylab = "Fraction", xlab = "Time (years)", title = titlePlot, ggtheme = theme_bw() + theme(plot.title = element_text(hjust =0.5)))$plot p } #a fucntion to assign unit, either ECAR or OCR or none giveUnit <- function(x) { ocrList <- c("basal.respiration","ATP.production","proton.leak","maximal.respiration", "spare.respiratory.capacity","OCR") ecarList <- c("glycolysis","glycolytic.capacity","glycolytic.reserve","ECAR") if (x %in% ocrList) return("OCR") else if (x %in% ecarList) return("ECAR") else return("") }
Stratrified by metabolism
TTT
tttList <- lapply(rownames(seaCombat), function(n) { survS <- data.frame(time = survT$T5, endpoint = survT$treatedAfter) val <- survT[[n]] ms <- maxstat.test(Surv(time, endpoint) ~ val, data = survS, smethod = "LogRank", minprop = 0.2, maxprop = 0.8, alpha = NULL) plotTab <- survS %>% mutate(response = val) %>% filter(!is.na(response), !is.na(endpoint),!is.na(time)) %>% mutate(group = ifelse(response >= ms$estimate, sprintf("high %s (%s >= %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate), sprintf("low %s (%s < %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate))) nTab <- table(plotTab$group) plotTab <- mutate(plotTab, group = paste0(group,nTab[group],")")) pval <- sprintf("p = %1.3f",ttt[formatSea(n),"p"]) p <- ggsurvplot(survfit(Surv(time, endpoint) ~ group, data = plotTab), data = plotTab, pval = pval, conf.int = FALSE, legend = c(0.6, 0.12),pval.coord = c(0.1,0.25), legend.labs = sort(unique(plotTab$group)), legend.title = "group", palette = "Dark2", ylab = "Fraction treatment free", xlab = "Time (years)", title = formatSea(n), ggtheme = theme_classic() + theme(axis.title = element_text(size=13, face="bold"), axis.text = element_text(size=12), plot.title = element_text(face="bold", hjust=0.5, size = 15 ), axis.title.x = element_text(face="bold"), plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")))$plot }) names(tttList) <- names(seaCombat)
OS
osList <- lapply(rownames(seaCombat), function(n) { survS <- data.frame(time = survT$T6, endpoint = survT$died) val <- survT[[n]] ms <- maxstat.test(Surv(time, endpoint) ~ val, data = survS, smethod = "LogRank", minprop = 0.2, maxprop = 0.8, alpha = NULL) plotTab <- survS %>% mutate(response = val) %>% filter(!is.na(response), !is.na(endpoint),!is.na(time)) %>% mutate(group = ifelse(response >= ms$estimate, sprintf("high %s (%s >= %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate), sprintf("low %s (%s < %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate))) nTab <- table(plotTab$group) plotTab <- mutate(plotTab, group = paste0(group,nTab[group],")")) pval <- sprintf("p = %1.3f",os[formatSea(n),"p"]) p <- ggsurvplot(survfit(Surv(time, endpoint) ~ group, data = plotTab), data = plotTab, pval = pval, conf.int = FALSE, legend = c(0.6,0.12),pval.coord = c(0.1,0.25), legend.labs = sort(unique(plotTab$group)), legend.title = "group", palette = "Dark2", ylab = "Fraction overall survival", xlab = "Time (years)", title = formatSea(n), ggtheme = theme_classic() + theme(axis.title = element_text(size=13, face="bold"), axis.text = element_text(size=12), plot.title = element_text(face="bold", hjust=0.5, size = 15 ), axis.title.x = element_text(face="bold"), plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")))$plot }) names(osList) <- names(seaCombat)
KM plot for significant assocations in single variante cox model
p<-plot_grid(tttList$glycolytic.reserve, tttList$maximal.respiration, tttList$spare.respiratory.capacity, osList$glycolytic.capacity, osList$glycolytic.reserve, labels = NULL) p
Straitified by both metabolism and IGHV status
TTT
tttKm <- lapply(rownames(seaCombat), function(n) { survS <- data.frame(time = survT$T5, endpoint = survT$treatedAfter) val <- survT[[n]] ms <- maxstat.test(Surv(time, endpoint) ~ val, data = survS, smethod = "LogRank", minprop = 0.2, maxprop = 0.8, alpha = NULL) plotTab <- survS %>% mutate(response = val, IGHV = survT$IGHV) %>% filter(!is.na(response),!is.na(IGHV), !is.na(endpoint),!is.na(time)) %>% mutate(group = ifelse(response >= ms$estimate, sprintf("high %s (%s >= %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate), sprintf("low %s (%s < %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate))) %>% mutate(group = ifelse(IGHV==1, paste0("M-CLL with ",group), paste0("U-CLL with ",group))) nTab <- table(plotTab$group) plotTab <- mutate(plotTab, group = paste0(group,nTab[group],")")) p <- ggsurvplot(survfit(Surv(time, endpoint) ~ group, data = plotTab), data = plotTab, pval = TRUE, conf.int = FALSE, legend = c(0.6,0.2),legend.labs = sort(unique(plotTab$group)), pval.coord = c(0.1,0.4), legend.title = "group", palette = "Dark2", ylab = "Fraction treatment free", xlab = "Time (years)", title = formatSea(n), ggtheme = theme_classic() + theme(axis.title = element_text(size=13, face="bold"), axis.text = element_text(size=12), plot.title = element_text(face="bold", hjust=0.5, size = 15 ), axis.title.x = element_text(face="bold"), plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")))$plot }) names(tttKm) <- names(seaCombat)
OS
osKm <- lapply(rownames(seaCombat), function(n) { survS <- data.frame(time = survT$T6, endpoint = survT$died) val <- survT[[n]] ms <- maxstat.test(Surv(time, endpoint) ~ val, data = survS, smethod = "LogRank", minprop = 0.2, maxprop = 0.8, alpha = NULL) plotTab <- survS %>% mutate(response = val, IGHV = survT$IGHV) %>% filter(!is.na(response),!is.na(IGHV), !is.na(endpoint),!is.na(time)) %>% mutate(group = ifelse(response >= ms$estimate, sprintf("high %s (%s >= %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate), sprintf("low %s (%s < %1.1f, n=",formatSea(n),giveUnit(n),ms$estimate))) %>% mutate(group = ifelse(IGHV==1, paste0("M-CLL with ",group), paste0("U-CLL with ",group))) nTab <- table(plotTab$group) plotTab <- mutate(plotTab, group = paste0(group,nTab[group],")")) p <- ggsurvplot(survfit(Surv(time, endpoint) ~ group, data = plotTab), data = plotTab, pval = TRUE, conf.int = FALSE, legend = c(0.6,0.2),legend.labs = sort(unique(plotTab$group)), pval.coord = c(0.1,0.4), legend.title = "group", palette = "Dark2", ylab = "Fraction overall survival", xlab = "Time (years)", title = formatSea(n), ggtheme = theme_classic() + theme(axis.title = element_text(size=13, face="bold"), axis.text = element_text(size=12), plot.title = element_text(face="bold", hjust=0.5, size = 15 ), axis.title.x = element_text(face="bold"), plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")))$plot }) names(osKm) <- names(seaCombat)
Multi-variate cox model
Prepare dataset
#add genetic variants to survival table survT$SF3B1 <- Biobase::exprs(lpdAll)[ "SF3B1", rownames(survT) ] survT$NOTCH1 <- Biobase::exprs(lpdAll)[ "NOTCH1", rownames(survT) ] survT$BRAF <- Biobase::exprs(lpdAll)[ "BRAF", rownames(survT) ] survT$TP53 <- Biobase::exprs(lpdAll)[ "TP53", rownames(survT) ] survT$del17p13 <- Biobase::exprs(lpdAll)[ "del17p13", rownames(survT) ] survT$del11q22.3 <- Biobase::exprs(lpdAll)[ "del11q22.3", rownames(survT) ] survT$trisomy12 <- Biobase::exprs(lpdAll)[ "trisomy12", rownames(survT) ]
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$IGHVwt <- ifelse(survT$IGHV==1, 0, 1) survT$treatment <- ifelse(survT$treatment, 1, 0)
TTT
glycolytic reserve
surv1 <- coxph( Surv(T5, treatedAfter) ~ age + as.factor(trisomy12) + as.factor(del11q22.3) + as.factor(del17p13) + as.factor(TP53) + as.factor(IGHVwt) + glycolytic.reserve, # continuous data = survT ) colFactor <- data.frame(factor = c("age", "trisomy12", "del11q22.3", "del17p13","TP53","U-CLL", "glycolytic reserve")) outTab <- cbind(colFactor,extractSome(surv1)) write(print(xtable(outTab, caption = "Multivariate Cox regression model for time to treatment with glycolytic reserve as a covariate"), include.rownames=FALSE, caption.placement = "top"), file = paste0(plotDir,"glyRes_TTT.tex")) cat(sprintf("%s patients considerd in the model; number of events %1g\n", summary(surv1)$n, summary(surv1)[6] ) )
maximal respiration
surv1 <- coxph( Surv(T5, treatedAfter) ~ age + as.factor(trisomy12) + as.factor(del11q22.3) + as.factor(del17p13) + as.factor(TP53) + IGHVwt + maximal.respiration, # continuous data = survT ) colFactor <- data.frame(factor = c("age", "trisomy12", "del11q22.3", "del17p13","TP53","U-CLL", "maximal respiration")) outTab <- cbind(colFactor,extractSome(surv1)) write(print(xtable(outTab, caption = "Multivariate Cox regression model for time to treatment with maximal respiration as a covariate"), include.rownames=FALSE, caption.placement = "top"), file = paste0(plotDir,"maxRes_TTT.tex")) cat(sprintf("%s patients considerd in the model; number of events %1g\n", summary(surv1)$n, summary(surv1)[6] ) )
spare respiratory capacity
surv1 <- coxph( Surv(T5, treatedAfter) ~ age + as.factor(trisomy12) + as.factor(del11q22.3) + as.factor(del17p13) + as.factor(TP53) + IGHVwt + spare.respiratory.capacity, # continuous data = survT ) colFactor <- data.frame(factor = c("age", "trisomy12", "del11q22.3", "del17p13","TP53","U-CLL", "spare.respiratory.capacity")) outTab <- cbind(colFactor,extractSome(surv1)) write(print(xtable(outTab, caption = "Multivariate Cox regression model for time to treatment with glycolytic reserve as a covariate"), include.rownames=FALSE, caption.placement = "top"), file = paste0(plotDir,"spResCap_TTT.tex")) cat(sprintf("%s patients considerd in the model; number of events %1g\n", summary(surv1)$n, summary(surv1)[6] ) )
OS
glycolytic reserve
surv1 <- coxph( Surv(T6, died) ~ age + as.factor(treatment) + as.factor(trisomy12) + as.factor(del11q22.3) + as.factor(del17p13) + as.factor(TP53) + IGHVwt + glycolytic.reserve, # continuous data = survT ) colFactor <- data.frame(factor = c("age", "treatment", "trisomy12", "del11q22.3", "del17p13","TP53","U-CLL", "glycolytic reserve")) outTab <- cbind(colFactor,extractSome(surv1)) write(print(xtable(outTab, caption = "Multivariate Cox regression model for overall survival with glycolytic reserve as a covariate"), include.rownames=FALSE, caption.placement = "top"), file = paste0(plotDir,"glyRes_OS.tex")) cat(sprintf("%s patients considerd in the model; number of events %1g\n", summary(surv1)$n, summary(surv1)[6] ) ) write.csv2(outTab, "glycolytic.reserve_VS_OS.csv")
glycolytic capacity
surv1 <- coxph( Surv(T6, died) ~ age + as.factor(treatment) + as.factor(trisomy12) + as.factor(del11q22.3) + as.factor(del17p13) + as.factor(TP53) + IGHVwt + glycolytic.capacity, # continuous data = survT ) colFactor <- data.frame(factor = c("age", "treatment", "trisomy12", "del11q22.3", "del17p13","TP53","U-CLL", "glycolytic capacity")) outTab <- cbind(colFactor,extractSome(surv1)) write(print(xtable(outTab, caption = "Multivariate Cox regression model for overall survival with glycolytic capacity as a covariate"), include.rownames=FALSE, caption.placement = "top"), file = paste0(plotDir,"glyCap_OS.tex")) cat(sprintf("%s patients considerd in the model; number of events %1g\n", summary(surv1)$n, summary(surv1)[6] )) write.csv2(outTab, "glycolytic.capacity_VS_OS.csv")
Association between bioenergetic features and CD38/CD49d expression
Prepare data
Query CD38 and CD49d (ITGA4) expressions from RNAseq data
data("dds") geneList <- c("CD38","ITGA4") ddsSea <- dds[,dds$PatID %in% colnames(seaOri)] ddsSea <- estimateSizeFactors(ddsSea) ddsSea <- varianceStabilizingTransformation(ddsSea) ddsSea <- ddsSea[rowData(ddsSea)$symbol %in% geneList,] countSea <- assay(ddsSea) rownames(countSea) <- rowData(ddsSea)$symbol colnames(countSea) <- ddsSea$PatID countSea <- data.frame(t(countSea)) %>% rownames_to_column("patID")
Combine the phenotype table and subset for patients with Seahorse measurement
phenoTab <- countSea[match(colnames(seaCombat), countSea$patID),] %>% filter(!is.na(patID)) %>% mutate(IGHV = Biobase::exprs(lpdAll)["IGHV Uppsala U/M",patID]) %>% mutate(IGHV = ifelse(is.na(IGHV),NA, ifelse(IGHV ==1, "M","U")))
Correlation test
seaTab <- assay(seaCombat)[,phenoTab$patID] stopifnot(all(colnames(seaTab) == phenoTab$patID)) corRes <- lapply(rownames(seaTab), function(seaName) { lapply(colnames(phenoTab)[2:3], function(phenoName){ testTab <- data.frame(sea = seaTab[seaName,], pheno = phenoTab[[phenoName]], IGHV = as.factor(phenoTab$IGHV)) lmRes <- summary(lm(pheno ~ sea, data = testTab, na.action = na.omit)) lmRes.IGHV <- summary(lm(pheno ~ sea + IGHV, data = testTab, na.action = na.omit)) data.frame(measurement = seaName, phenotype = phenoName, p = lmRes$coefficients[2,4], p.IGHV = lmRes.IGHV$coefficients[2,4], stringsAsFactors = FALSE) }) %>% bind_rows() }) %>% bind_rows() %>% arrange(p) %>% mutate(p.adj = p.adjust(p, method = "BH"), p.IGHV.adj = p.adjust(p.IGHV, method = "BH"))
Save a table for supplementary table
tabOut <- filter(corRes, p.adj <= 0.05) %>% mutate(measurement = formatSea(measurement)) %>% dplyr::rename(Measurement = measurement, Gene = phenotype, `p value` = p, `p value (IGHV blocked)` = p.IGHV, `adjusted p value` = p.adj, `adjusted p value (IGHV blocked)`= p.IGHV.adj) write(print(xtable(tabOut, digits = 3, caption = "Correlation test results between energy metabolic measurements and CD38/IGTA4(CD49d) expression"), include.rownames=FALSE, caption.placement = "top"), file = paste0(plotDir,"corTest_seaVSCD38.tex"))
Plot significant correlations
corRes.sig <- filter(corRes, p.IGHV.adj < 0.05) pList.all <- lapply(seq(nrow(corRes.sig)), function(i) { seaName <- corRes.sig$measurement[i] phenoName <- corRes.sig$phenotype[i] plotTab <- data.frame(sea = seaTab[seaName,], pheno = phenoTab[[phenoName]], IGHV = phenoTab$IGHV) %>% mutate(IGHV = ifelse(is.na(IGHV),"unknown",ifelse( IGHV == "M", "mutated","unmutated"))) corRes <- cor.test(plotTab$sea, plotTab$pheno) annoText <- paste("'coefficient ='~",format(corRes$estimate,digits = 2),"*","', p ='~",sciPretty(corRes$p.value,digits=2)) ggplot(plotTab, aes(x=sea, y = pheno)) + geom_point(aes(color = IGHV)) + geom_smooth(method = "lm", se=FALSE ) + theme_bw() + ylab(paste0(phenoName," (normalized expression)")) + xlab("ECAR (pMol/min)") + scale_color_manual(values = c(mutated = "red", unmutated = "black", unknown = "grey80")) + ggtitle(formatSea(seaName)) + annotate("text", x = -Inf, y = Inf, label = annoText, vjust=1, hjust=0, size = 5, parse = TRUE, col= "darkred") + theme(panel.grid = element_blank(), axis.title = element_text(size=13, face="bold"), axis.text = element_text(size=12), plot.title = element_text(face="bold", hjust=0.5, size = 15 ), legend.position = c(0.85,0.15), legend.background = element_rect(fill = NA), axis.title.x = element_text(face="bold"), plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) }) grid.arrange(grobs = pList.all, ncol=2)
pList.U <- lapply(seq(nrow(corRes.sig)), function(i) { seaName <- corRes.sig$measurement[i] phenoName <- corRes.sig$phenotype[i] plotTab <- data.frame(sea = seaTab[seaName,], pheno = phenoTab[[phenoName]], IGHV = phenoTab$IGHV) %>% filter(IGHV == "U") corRes <- cor.test(plotTab$sea, plotTab$pheno) annoText <- paste("'coefficient ='~",format(corRes$estimate,digits = 2),"*","', p ='~",sciPretty(corRes$p.value,digits=2)) ggplot(plotTab, aes(x=sea, y = pheno)) + geom_point(color = "black") + geom_smooth(method = "lm", se=FALSE ) + theme_bw() + ylab(paste0(phenoName," (normalized expression)")) + xlab("ECAR (pMol/min)") + scale_color_manual(values = c(mutated = "red", unmutated = "black", unknown = "grey80")) + ggtitle(paste0(formatSea(seaName)," (U-CLL samples only)")) + annotate("text", x = -Inf, y = Inf, label = annoText, vjust=1, hjust=0, size = 5, parse = TRUE, col= "darkred") + theme(panel.grid = element_blank(), axis.title = element_text(size=13, face="bold"), axis.text = element_text(size=12), plot.title = element_text(face="bold", hjust=0.5, size = 15 ), legend.position = c(0.9,0.1), axis.title.x = element_text(face="bold"), plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) }) pList.M <- lapply(seq(nrow(corRes.sig)), function(i) { seaName <- corRes.sig$measurement[i] phenoName <- corRes.sig$phenotype[i] plotTab <- data.frame(sea = seaTab[seaName,], pheno = phenoTab[[phenoName]], IGHV = phenoTab$IGHV) %>% filter(IGHV == "M") corRes <- cor.test(plotTab$sea, plotTab$pheno) annoText <- paste("'coefficient ='~",format(corRes$estimate,digits = 2),"*","', p ='~",sciPretty(corRes$p.value,digits=2)) ggplot(plotTab, aes(x=sea, y = pheno)) + geom_point(color = "red") + geom_smooth(method = "lm", se=FALSE ) + theme_bw() + ylab(paste0(phenoName," (normalized expression)")) + xlab("ECAR (pMol/min)") + scale_color_manual(values = c(mutated = "red", unmutated = "black", unknown = "grey80")) + ggtitle(paste0(formatSea(seaName)," (M-CLL samples only)")) + annotate("text", x = -Inf, y = Inf, label = annoText, vjust=1, hjust=0, size = 5, parse = TRUE, col= "darkred") + theme(panel.grid = element_blank(), axis.title = element_text(size=13, face="bold"), axis.text = element_text(size=12), plot.title = element_text(face="bold", hjust=0.5, size = 15 ), legend.position = c(0.9,0.1), axis.title.x = element_text(face="bold"), plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")) }) grid.arrange(grobs = c(pList.U,pList.M), ncol=2)
Organize a figure for the clinical part
figOut <- c(osKm$glycolytic.capacity, osKm$glycolytic.reserve, pList.all[[1]],pList.all[[2]]) title = ggdraw() + draw_figure_label("Figure 5", fontface = "bold", position = "top.left",size=20) p<-plot_grid(osKm$glycolytic.capacity, osKm$glycolytic.reserve, pList.all[[1]], pList.all[[2]], labels = c("A","B","C","D"), label_size = 22) plot_grid(title, p, rel_heights = c(0.05,0.95), ncol = 1)
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.