analysis.diff.survival.TCGA <- function(interest.gene, diff.gene.pro, exp.data.process, clin.data, group.by = "median", EnhancedVolcano.plot = T, Box.plot = T, main = NULL, meta.signature = F, single.signature = T){
require(ggpubr)
# require interest.gene to be data.frame
if(EnhancedVolcano.plot){
require(EnhancedVolcano)
gene.overlap <- intersect(rownames(diff.gene.pro), rownames(interest.gene))
diff.gene.pro.sig.gene <- diff.gene.pro[gene.overlap,]
interest.gene.sig <- interest.gene[gene.overlap,]
up.gene <- which(diff.gene.pro.sig.gene$log2FoldChange>0)
if(length(up.gene)>0){
down.gene <- setdiff(1:nrow(diff.gene.pro.sig.gene), up.gene)
interest.gene.sig$label <- rep("down", nrow(diff.gene.pro.sig.gene))
interest.gene.sig$label[up.gene] <- "up"
}else{
down.gene <- 1:nrow(diff.gene.pro.sig.gene)
interest.gene.sig$label <- rep("down", nrow(diff.gene.pro.sig.gene))
}
interest.gene.sig$TCGA_log2FoldChange <- diff.gene.pro.sig.gene$log2FoldChange
interest.gene.sig$TCGA_padj <- diff.gene.pro.sig.gene$padj
p <- EnhancedVolcano(diff.gene.pro.sig.gene,
lab = rownames(diff.gene.pro.sig.gene),
selectLab = rownames(diff.gene.pro.sig.gene),
x = 'log2FoldChange',
y = 'padj', pCutoff = 0.05, drawConnectors = TRUE, FCcutoff = 1,
widthConnectors = 0.2,
title = "Tumor VS Normal", xlab = "", ylab = "-Log10(padj)")
print(p)
if(length(gene.overlap) < nrow(interest.gene)){
a <- setdiff(rownames(interest.gene), gene.overlap)
cat("gene ", paste(a, sep = ";"), " not in TCGA!\n")
}else{
cat("All gene in TCGA data!\n")
}
}else{
#Require interest.gene to be a character vector
gene.overlap <- intersect(interest.gene, rownames(diff.gene.pro))
diff.gene.pro.sig.gene <- diff.gene.pro[gene.overlap,]
interest.gene.sig <- data.frame(Gene = gene.overlap)
if(length(gene.overlap) < length(interest.gene)){
a <- setdiff(interest.gene, gene.overlap)
cat("gene ", paste(a, sep = ";"), " not in TCGA!\n")
}else{
cat("All gene in TCGA data!\n")
}
}
##Draw the expression barplot of each gene
exp.data.interest <- exp.data.process[gene.overlap,]
sample.id <- colnames(exp.data.process)
sample.id <- substr(sample.id, 1, 16)
sample.type <- substr(sample.id, 14, 15)
normal.index <- which(sample.type == "11") #72
cancer.index <- which(sample.type == "01") #533
patient.type <- rep("Tumor", ncol(exp.data.process))
patient.type[normal.index] <- "Normal"
if(Box.plot){
cat("Normal sample: ", length(normal.index), "\n")
cat("Tumor sample: ", length(cancer.index), "\n")
box.res <- sapply(gene.overlap, function(x){
interest.gene.matrix <- exp.data.process[x,]
genes <- data.frame(type = patient.type, expressionValue = interest.gene.matrix)
p <- ggboxplot(genes, x = "type", y = "expressionValue", ylab = "Normalized Count", color = "type", palette = "jco", add = "jitter", title = x) + stat_compare_means()
print(p)
})
}
##对应表达和生存数据
sample.id <- substr(sample.id, 1, 15)
patient.overlap <- intersect(clin.data$Sample, sample.id) # 533 tumor samples
exp.data.process.pro <- exp.data.process[,match(patient.overlap, sample.id)]
clin.info <- clin.data[match(patient.overlap, clin.data$Sample),]
#all gene as the signature
if(meta.signature){
interest.matrix <- exp.data.process.pro[gene.overlap,]
signature.score <- colSums(interest.matrix)/nrow(interest.matrix)
median.score <- median(signature.score)
high.group <- which(signature.score > median.score)
sample.label <- rep("Low group", length(signature.score))
sample.label[high.group] <- "High group"
OS.data <- data.frame(Patient_ID = patient.overlap, event = clin.info$OS, time = clin.info$OS_time, sample.label = sample.label)
DFS.data <- data.frame(Patient_ID = patient.overlap, event = clin.info$DFS, time = clin.info$DFS_time, sample.label = sample.label)
p1 <- plot.surv(OS.data, risk.table = T, HR = T, ylab = "Overall Survival", main = main, surv.median.line = "hv", xlab = "Time (Month)")
print(p1)
p2 <- plot.surv(DFS.data, risk.table = T, HR = T, ylab = "Disease-Free Survival", main = main, surv.median.line = "hv", xlab = "Time (Month)")
print(p2)
}
## single gene model
if(single.signature){
suv.res <- sapply(gene.overlap, function(x){
gene.expression <- as.numeric(exp.data.process.pro[x,])
if(group.by == "top33"){
quantile.value <- quantile(gene.expression, seq(0,1,0.33))
top33 <- as.numeric(quantile.value[3])
bottom33 <- as.numeric(quantile.value[2])
index1 <- which(gene.expression >= top33)
index2 <- which(gene.expression <= bottom33)
sample.label <- c(rep("Low group", length(index1)), rep("High group", length(index2)))
clin.info <- clin.info[c(index2,index1),]
OS.data <- data.frame(Patient_ID = patient.overlap[c(index2,index1)], event = clin.info$OS, time = clin.info$OS_time, sample.label = sample.label)
DFS.data <- data.frame(Patient_ID = patient.overlap[c(index2,index1)], event = clin.info$DFS, time = clin.info$DFS_time, sample.label = sample.label)
}else{
# median as threshold
med.exp <- median(gene.expression)
high.group <- which(gene.expression>med.exp)
sample.label <- rep("Low group", length(gene.expression))
sample.label[high.group] <- "High group"
OS.data <- data.frame(Patient_ID = patient.overlap, event = clin.info$OS, time = clin.info$OS_time, sample.label = sample.label)
DFS.data <- data.frame(Patient_ID = patient.overlap, event = clin.info$DFS, time = clin.info$DFS_time, sample.label = sample.label)
}
p1 <- plot.surv(OS.data, risk.table = T, HR = T, ylab = "Overall Survival", main = x, surv.median.line = "hv", xlab = "Time (Month)")
print(p1)
p2 <- plot.surv(DFS.data, risk.table = T, HR = T, ylab = "Disease-Free Survival", main = x, surv.median.line = "hv", xlab = "Time (Month)")
print(p2)
OS.obj <- coxph(Surv(time, event)~sample.label, data=OS.data)
DFS.obj <- coxph(Surv(time, event)~sample.label, data=DFS.data)
return(c(summary(OS.obj)$logtest[3], summary(DFS.obj)$logtest[3]))
})
interest.gene.sig$OS_logrank_p <- suv.res[1,]
interest.gene.sig$DFS_logrank_p <- suv.res[2,]
}
return(interest.gene.sig)
}
plot.surv <- function(clinical.data, upper.time = NULL, xscale = 1, xlab = "Time", median.time = TRUE,
surv.median.line = "none", HR = FALSE, risk.table = TRUE, pval = TRUE,
conf.int = FALSE, main = NULL, ylab = "Survival probability", colors = c("#D95319", "#3B6793","#EA4335","#4285F4","#34A853","#000000")) {
#load R package
require(survival)
require(survminer)
require(RColorBrewer)
require(gridExtra)
# Determine the type of event and the unit of time
# survival.event <- survival.event[1];
# unit.xlabel <- unit.xlabel[1];
# If upper.time is set, the samples whose survival time exceeds upper.time will be removed
if (!is.null(upper.time)) clinical.data <- clinical.data[clinical.data$time <= upper.time,]
# #date format
# xSL <- data.frame(xScale=c(1,7,30,365.25),xLab=c("Days","Weeks","Months","Years"), stringsAsFactors=FALSE)
# switch(unit.xlabel, year={xScale <- 365.25;}, month={xScale <- 30;}, week={xScale <- 7;}, day={xScale <- 1})
# xLab <- xSL[which(xSL[,1]==xScale),2];
# color
if (!is.factor(clinical.data$sample.label))
clinical.data$sample.label <- as.factor(clinical.data$sample.label)
t.name <- levels(clinical.data$sample.label)
if (length(t.name) > 6) stop("Sample grouping >6, exceeding the function acceptance range!")
t.col <- colors[1:length(t.name)]
# 构造生存对象
km.curves <- survfit(Surv(time, event)~sample.label, data=clinical.data)
# HR and 95%CI
if (length(t.name) == 2) {
if (HR) {
cox.obj <- coxph(Surv(time, event)~sample.label, data=clinical.data)
tmp <- summary(cox.obj)
HRs <- round(tmp$coefficients[ ,2], digits = 2)
HR.confint.lower <- round(tmp$conf.int[,"lower .95"], 2)
HR.confint.upper <- round(tmp$conf.int[,"upper .95"], 2)
HRs <- paste0(HRs, " (", HR.confint.lower, "-", HR.confint.upper, ")")
}
}
# Construct the legend display text in the survival image
legend.content <- substr(names(km.curves$strata), start = 14, stop = 1000)
# x-axis scale unit conversion
if (is.numeric(xscale) | (xscale %in% c("d_m", "d_y", "m_d", "m_y", "y_d", "y_m"))) {
xscale = xscale
} else {
stop('xscale should be numeric or one of c("d_m", "d_y", "m_d", "m_y", "y_d", "y_m").')
}
# conversion of survival time units
.format_xticklabels <- function(labels, xscale){
# 1 year = 365.25 days
# 1 month = 365.25/12 = 30.4375 days
if (is.numeric(xscale)) xtrans <- 1/xscale
else
xtrans <- switch(xscale,
d_m = 12/365.25,
d_y = 1/365.25,
m_d = 365.25/12,
m_y = 1/12,
y_d = 365.25,
y_m = 12,
1
)
round(labels*xtrans,2)
}
# Add the median survival time and its 95% CI in the figure and place it in the subtitle position
subtitle <- NULL
if (median.time) {
if (is.numeric(xscale)) {
median.km.obj = km.curves
} else if (xscale %in% c("d_m", "d_y", "m_d", "m_y", "y_d", "y_m")) {
clinical.data$time <- .format_xticklabels(labels = clinical.data$time, xscale = xscale)
median.km.obj <- survfit(Surv(time, event)~sample.label, data=clinical.data)
}
survival.time.info <- NULL
survival.time.info <- rbind(survival.time.info, summary(median.km.obj)$table)
median.survival <- round(survival.time.info[!duplicated(survival.time.info[,7:9]),7:9], digits = 2) # 注意：这里取得的置信区间上界可能为NA
if (length(levels(clinical.data$sample.label)) == 1) {
tmp1 <- levels(clinical.data$sample.label)
} else {
tmp1 <- do.call(rbind,strsplit(rownames(summary(median.km.obj)$table), split = "="))[,2]
}
tmp2 <- paste(median.survival[,1], "(", median.survival[,2], "-", median.survival[,3], ")")
subtitle <- paste(tmp1, tmp2, sep = ":", collapse = "\n")
}
# ggsurvplot
ggsurv <- ggsurvplot(km.curves, # survfit object with calculated statistics.
data = clinical.data, # data used to fit survival curves.
palette = t.col,
risk.table = risk.table, # show risk table.
pval = pval, # show p-value of log-rank test.
surv.median.line = surv.median.line, # add the median survival pointer.
title = main, #main title
subtitle = subtitle, #sub title
font.main = 15,
xlab = xlab, # customize X axis label.
ylab = ylab, # customize Y axis label
xscale = xscale,
#legend
legend.title = "",
legend.labs = legend.content,
legend = c(0.8,0.9),
font.legend = 9,
#risk table
tables.theme = theme_cleantable(),
risk.table.title = "No. at risk:",
risk.table.y.text.col = T,
risk.table.y.text = FALSE,
tables.height = 0.15,
risk.table.fontsize = 3
)
if (length(t.name) == 2) {
if (HR)
ggsurv$plot <- ggsurv$plot + ggplot2::annotate("text", x = max(km.curves$time)/12,
y = 0.15, size = 5, label = paste("HR=", HRs))
}
ggsurv$plot <- ggsurv$plot + theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(size = 10),
plot.margin = unit(c(5.5, 5.5, 5.5, 50), "points"))
ggsurv$table <- ggsurv$table + theme(plot.title = element_text(hjust = -0.04),
plot.margin = unit(c(5.5, 5.5, 5.5, 50), "points"))
if(length(t.name) > 2) {
# pairwise: log rank P value
res <- pairwise_survdiff(Surv(time, event)~sample.label, data=clinical.data);
pairwise.pvalue <- round(res$p.value, digits = 4);
pairwise.pvalue[which(pairwise.pvalue < 0.0001)] <- "<0.0001";
pairwise.pvalue[is.na(pairwise.pvalue)] <- "-"
# add table
tt <- ttheme_minimal(core = list(fg_params = list(col = "black"),bg_params = list(fill = NA, col = "black")),
colhead = list(fg_params = list(col = NA),bg_params = list(fill = t.col, col = "black")),
rowhead = list(fg_params = list(col = NA, hjust = 1),bg_params = list(fill = c("white",t.col[-1]), col = "black"))
)
pairwise.table <- tableGrob(pairwise.pvalue, theme = tt)
ggsurv <- ggarrange(ggarrange(ggsurv$plot, ggsurv$table, nrow=2, heights=c(2,0.5)),
pairwise.table, nrow=2, heights = c(2,0.5),
labels = c("","p from pairwise comparisons"),
hjust = 0, font.label = list(size = 15, face = "plain"))
}
ggsurv
}
Datasets
Models
