plot44 · 代码 0 · 图 0 · 笔记 5 · 代码库专题
HTAN Lung：TP53 相关肺腺癌细胞与空间图谱
论文：Zhao et al., “A cellular and spatial atlas of TP53-associated tissue remodeling defines a multicellular tumor ecosystem in lung adenocarcinoma”, Nature Cancer 6, 1857–1879 (2025). [论文原文](https://ww
HTAN Lung：TP53 相关肺腺癌细胞与空间图谱（自动生成占位图）
图源：（自动生成占位图）
分析代码
8 个源码文件，按目录顺序；本机绝对路径已替换为 ~。 代码可复制，但需要按自己的数据改路径与参数。
figure_catalog/extended/bulk__3a/plot_code.r
# EXTENDED FIG 3A, D, E
v.s3ade <- c("p53_targets", "AT2", "CC.G2M", "`Glycolysis/Hypox`", "pEMT") # row 1
levels.tmp <- levels(luad.2018@meta.data$comut)
createBoxplots(plotVector = v.s3ade,
xVar = "comut",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x))[-6],
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "3a")
figure_catalog/extended/bulk__3a/source_original.r
######################################################
#################### TSANKOV LAB #####################
######################################################
# ----------------------------------------------------
# load the libraries
# ----------------------------------------------------
library(Seurat)
library(tidyverse)
library(ggpubr)
# ----------------------------------------------------
# set directories
# ----------------------------------------------------
data.path <- "../data/"
figures.dir <- "../../figures/"
system(paste0("mkdir -p ", figures.dir))
# ----------------------------------------------------
# load the objects
# ----------------------------------------------------
luad.2018 <- readRDS(file.path(data.path, "luad.2018.processed.rds"))
cptac.luad.prot <- readRDS(file.path(data.path, "cptac.luad.prot.processed.rds"))
# ----------------------------------------------------
# load the dataframe
# ----------------------------------------------------
df.pvals.tcga <- read_csv(file.path(data.path, "tcga.regression_plots.p53_kras_egfr_stk11_CancerProp.filtered.csv")) |>
select(-"...1")
# ----------------------------------------------------
# custom function
# ----------------------------------------------------
createBoxplots <- function(plotVector,
plotDF,
outfileName,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
df.pvals = NULL,
numCol = 4,
w = 3.5,
h = 3) {
plotlist <- list()
for (plotVar in plotVector) {
plotVar.clean <- str_remove_all(string = plotVar, pattern = "`")
if (!plotVar.clean %in% colnames(plotDF)) {
message(paste0(plotVar.clean, " not found in plotDF!"))
next
}
message(paste0("Adding plot and saving csv file for ", plotVar, "..."))
saveDF <- plotDF |>
drop_na(all_of(xVar)) |>
select(all_of(c(xVar, plotVar.clean)))
# write.csv(saveDF, file.path(paste0(data.path, outfileName, ".", make.names(plotVar.clean), ".csv")))
if (length(levels(plotDF[[xVar]])) == 2) {
if (is.null(df.pvals)) {
# for cptac
plotlist[[plotVar]] <- plotDF |>
drop_na(all_of(xVar)) |>
ggboxplot(x = xVar,
y = plotVar,
color = xVar,
palette = c("#0091CA", "#D8423D"),
title = plotVar) +
stat_compare_means(comparisons = comparisonList) +
scale_x_discrete(labels = function(x) paste0(x, "\n(n = ", table(plotDF[[xVar]])[x], ")")) +
rotate_x_text(0) +
xlab(NULL) +
theme(plot.title = element_text(size = 12, face = "bold"), legend.position = "none")
} else {
# for tcga
pvals_vector <- setNames(df.pvals$pval, df.pvals$program)
pvals_manual <- pvals_vector[[plotVar.clean]]
plotlist[[plotVar]] <- plotDF |>
drop_na(all_of(xVar)) |>
ggboxplot(x = xVar,
y = plotVar,
color = xVar,
palette = c("#0091CA", "#D8423D"),
title = plotVar) +
# stat_compare_means(comparisons = comparisonList) +
scale_x_discrete(labels = function(x) paste0(x, "\n(n = ", table(plotDF[[xVar]])[x], ")")) +
rotate_x_text(0) +
xlab(NULL) +
theme(plot.title = element_text(size = 12, face = "bold"), legend.position = "none") +
ggsignif::geom_signif(comparisons = list(c("WT", "mut")), # Replace "WT" and "mut" with actual levels if different
annotations = paste0("p = ", signif(pvals_manual, digits = 3)),
y_position = max(plotDF[[plotVar.clean]], na.rm = TRUE) * 1.1, # Adjust as needed
tip_length = 0.02, # Adjust the length of the horizontal line tips
textsize = 3)
}
} else {
plotlist[[plotVar]] <- plotDF |>
drop_na(all_of(xVar)) |>
ggboxplot(x = xVar,
y = plotVar,
color = xVar,
title = plotVar) +
stat_compare_means(comparisons = comparisonList) +
scale_x_discrete(labels = function(x) paste0(x, "\n(n = ", table(plotDF[[xVar]])[x], ")")) +
rotate_x_text(90) +
xlab(NULL) +
theme(plot.title = element_text(size = 12, face = "bold"), legend.position = "none")
}
}
ggsave(filename = file.path(paste0(figures.dir, outfileName, ".pdf")),
plot = patchwork::wrap_plots(plotlist, ncol = numCol),
device = "pdf",
width = w * numCol,
height = h * ceiling(length(plotlist) / 4))
}
# ----------------------------------------------------
# generate figures
# ----------------------------------------------------
# FIG 1F
v.1f <- c("Endothelial", "Pericytes")
createBoxplots(plotVector = v.1f, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 2,
plotDF = luad.2018@meta.data,
outfileName = "1f")
# FIG 2E
v.2e <- c("p53_targets", "AT2", "CC.G2M", "`Glycolysis/Hypox`", "pEMT")
createBoxplots(plotVector = v.2e, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "2e")
# EXTENDED FIG 3A, D, E
v.s3ade <- c("p53_targets", "AT2", "CC.G2M", "`Glycolysis/Hypox`", "pEMT") # row 1
levels.tmp <- levels(luad.2018@meta.data$comut)
createBoxplots(plotVector = v.s3ade,
xVar = "comut",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x))[-6],
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "3a")
levels.tmp <- levels(luad.2018@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s3ade,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "3d")
levels.tmp <- levels(luad.2018@meta.data$impact)
createBoxplots(plotVector = v.s3ade,
xVar = "impact",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "3e")
# FIG 5D
v.5d <- c("SPP1")
createBoxplots(plotVector = v.5d, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "5d")
# EXTENDED FIG 1D
v.s1d <- c("B.cells", "Fibroblast", "Cancer", "Mast", "`T.CD4+`", # row 1
"Myeloid", "NK.cells", "Smooth.muscle", "T.cells", "`T.CD8+`") # row 2
createBoxplots(plotVector = v.s1d, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "s1d")
# EXTENDED FIG 3B
v.s3b <- c("p53_targets", "AT2", "CC.G2M", "Glycolysis.Hypox", "pEMT")
createBoxplots(plotVector = v.s3b,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 5,
plotDF = cptac.luad.prot@meta.data,
outfileName = "s3b")
# EXTENDED FIG 3C
v.s3c <- c("p53_targets", "AT2", "CC.G2M")
levels.tmp <- levels(cptac.luad.prot@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s3c, w = 3.5, h = 5,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 3,
plotDF = cptac.luad.prot@meta.data,
outfileName = "s3cA")
levels.tmp <- levels(cptac.luad.prot@meta.data$impact)
createBoxplots(plotVector = v.s3c, w = 3.5, h = 5,
xVar = "impact",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 3,
plotDF = cptac.luad.prot@meta.data,
outfileName = "s3cB")
# EXTENDED FIG 4A
v.s4a <- c("Endothelial")
levels.tmp <- levels(luad.2018@meta.data$comut)
createBoxplots(plotVector = v.s4a, w = 3, h = 5,
xVar = "comut",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x))[-6],
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s4aA")
levels.tmp <- levels(luad.2018@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s4a, w = 3, h = 5,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s4aB")
# EXTENDED FIG 4E
v.s4e <- c("Aerocyte", "Arterial")
createBoxplots(plotVector = v.s4e, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 2,
plotDF = luad.2018@meta.data,
outfileName = "s4e")
# EXTENDED FIG 5G
v.s5g <- c("CAF.ADH1B")
createBoxplots(plotVector = v.s5g, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s5g")
# EXTENDED FIG 5H
v.s5h <- c("Pericytes")
levels.tmp <- levels(luad.2018@meta.data$comut)
createBoxplots(plotVector = v.s5h, w = 3, h = 5,
xVar = "comut",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x))[-6],
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s5hA")
levels.tmp <- levels(luad.2018@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s5h, w = 3, h = 5,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s5hB")
# EXTENDED FIG 6B
v.s6b <- c("TAM.CXCL")
createBoxplots(plotVector = v.s6b, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s6b")
# EXTENDED FIG 6C
v.s6c <- c("CXCL9", "CXCL10", "CXCL11")
createBoxplots(plotVector = v.s6c, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 3,
plotDF = luad.2018@meta.data,
outfileName = "s6c")
# EXTENDED FIG 6D,E
v.s6de <- c("SPP1", "CXCL9", "CXCL10", "CXCL11")
levels.tmp <- levels(luad.2018@meta.data$comut)
createBoxplots(plotVector = v.s6de, w = 3, h = 5,
xVar = "comut",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x))[-6],
numCol = 4,
plotDF = luad.2018@meta.data,
outfileName = "s6d")
levels.tmp <- levels(luad.2018@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s6de, w = 3, h = 5,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 4,
plotDF = luad.2018@meta.data,
outfileName = "s6eA")
levels.tmp <- levels(luad.2018@meta.data$impact)
createBoxplots(plotVector = v.s6de, w = 3, h = 5,
xVar = "impact",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 4,
plotDF = luad.2018@meta.data,
outfileName = "s6eB")
# EXTENDED FIG 9D
v.s9d <- c("T.Exhausted", "CD8.GZMK")
createBoxplots(plotVector = v.s9d, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 2,
plotDF = luad.2018@meta.data,
outfileName = "s9d")
# EXTENDED FIG 7F
v.s7f <- c("CD274", "CD86", "PVR", "PDCD1", "CTLA4", "TIGIT")
createBoxplots(plotVector = v.s7f, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 6,
plotDF = luad.2018@meta.data,
outfileName = "s7f")
# EXTENDED FIG 7G
v.s7g <- c("PVR", "CD274")
createBoxplots(plotVector = v.s7g, w = 3, h = 4,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 2,
plotDF = cptac.luad.prot@meta.data,
outfileName = "s7g")
# ----------------------------------------------------
# smoking pack years
# ----------------------------------------------------
smoking_comparison_p53 <- read_csv(file.path(data.path, "smoking_comparison_p53.csv"))
p53.filt <- c("MGH1172", "MGH1173", "MGH1175", "MGH1182", "BWH16", "BWH11", "BWH14", "BWH06")
WT.filt <- c("P14", "MGH1174", "MGH1176", "MGH1183", "BWH01", "BWH04", "BWH05", "BWH09", "BWH19", "BWH23") # no MGH1170
xVar <- "p53_status"
plotVar <- "pack_yrs"
comparisonList = list(c("WT", "mut"))
plotDF <- smoking_comparison_p53 |>
# mutate(p53_status = case_when(Patient %in% c(p53.filt, "P17", "P19", "MGH1181") ~ "mut", Patient %in% c(WT.filt, "MGH1179") ~ "WT", .default = NA)) |>
mutate(p53_status = case_when(Patient %in% p53.filt ~ "mut", Patient %in% WT.filt ~ "WT", .default = NA)) |>
rename("pack_yrs" = "Pack yrs")
p.smoking_pack_yrs <- plotDF |>
drop_na(all_of(xVar)) |>
ggboxplot(x = xVar,
y = plotVar,
color = xVar,
palette = c("#0091CA", "#D8423D"),
add = "jitter",
title = plotVar) +
stat_compare_means(comparisons = comparisonList) +
scale_x_discrete(labels = function(x) paste0(x, "\n(n = ", table(plotDF[[xVar]])[x], ")")) +
rotate_x_text(0) +
xlab(NULL) +
ylab("Pack years") +
theme(plot.title = element_text(size = 12, face = "bold"), legend.position = "none")
p.smoking_pack_yrs
ggsave(filename = file.path(paste0(figures.dir, "smoking.pack_years.pdf")),
plot = p.smoking_pack_yrs,
device = "pdf",
width = 3,
height = 4)
# plotDF |>
# drop_na(all_of(xVar)) |>
# select(Patient, pack_yrs, p53_status) |>
# write.csv(file.path(data.path, "smoking.pack_years.csv"), row.names = F)
# # ----------------------------------------------------
figure_catalog/extended/bulk__3d/plot_code.r
levels.tmp <- levels(luad.2018@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s3ade,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "3d")
figure_catalog/extended/bulk__3d/source_original.r
######################################################
#################### TSANKOV LAB #####################
######################################################
# ----------------------------------------------------
# load the libraries
# ----------------------------------------------------
library(Seurat)
library(tidyverse)
library(ggpubr)
# ----------------------------------------------------
# set directories
# ----------------------------------------------------
data.path <- "../data/"
figures.dir <- "../../figures/"
system(paste0("mkdir -p ", figures.dir))
# ----------------------------------------------------
# load the objects
# ----------------------------------------------------
luad.2018 <- readRDS(file.path(data.path, "luad.2018.processed.rds"))
cptac.luad.prot <- readRDS(file.path(data.path, "cptac.luad.prot.processed.rds"))
# ----------------------------------------------------
# load the dataframe
# ----------------------------------------------------
df.pvals.tcga <- read_csv(file.path(data.path, "tcga.regression_plots.p53_kras_egfr_stk11_CancerProp.filtered.csv")) |>
select(-"...1")
# ----------------------------------------------------
# custom function
# ----------------------------------------------------
createBoxplots <- function(plotVector,
plotDF,
outfileName,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
df.pvals = NULL,
numCol = 4,
w = 3.5,
h = 3) {
plotlist <- list()
for (plotVar in plotVector) {
plotVar.clean <- str_remove_all(string = plotVar, pattern = "`")
if (!plotVar.clean %in% colnames(plotDF)) {
message(paste0(plotVar.clean, " not found in plotDF!"))
next
}
message(paste0("Adding plot and saving csv file for ", plotVar, "..."))
saveDF <- plotDF |>
drop_na(all_of(xVar)) |>
select(all_of(c(xVar, plotVar.clean)))
# write.csv(saveDF, file.path(paste0(data.path, outfileName, ".", make.names(plotVar.clean), ".csv")))
if (length(levels(plotDF[[xVar]])) == 2) {
if (is.null(df.pvals)) {
# for cptac
plotlist[[plotVar]] <- plotDF |>
drop_na(all_of(xVar)) |>
ggboxplot(x = xVar,
y = plotVar,
color = xVar,
palette = c("#0091CA", "#D8423D"),
title = plotVar) +
stat_compare_means(comparisons = comparisonList) +
scale_x_discrete(labels = function(x) paste0(x, "\n(n = ", table(plotDF[[xVar]])[x], ")")) +
rotate_x_text(0) +
xlab(NULL) +
theme(plot.title = element_text(size = 12, face = "bold"), legend.position = "none")
} else {
# for tcga
pvals_vector <- setNames(df.pvals$pval, df.pvals$program)
pvals_manual <- pvals_vector[[plotVar.clean]]
plotlist[[plotVar]] <- plotDF |>
drop_na(all_of(xVar)) |>
ggboxplot(x = xVar,
y = plotVar,
color = xVar,
palette = c("#0091CA", "#D8423D"),
title = plotVar) +
# stat_compare_means(comparisons = comparisonList) +
scale_x_discrete(labels = function(x) paste0(x, "\n(n = ", table(plotDF[[xVar]])[x], ")")) +
rotate_x_text(0) +
xlab(NULL) +
theme(plot.title = element_text(size = 12, face = "bold"), legend.position = "none") +
ggsignif::geom_signif(comparisons = list(c("WT", "mut")), # Replace "WT" and "mut" with actual levels if different
annotations = paste0("p = ", signif(pvals_manual, digits = 3)),
y_position = max(plotDF[[plotVar.clean]], na.rm = TRUE) * 1.1, # Adjust as needed
tip_length = 0.02, # Adjust the length of the horizontal line tips
textsize = 3)
}
} else {
plotlist[[plotVar]] <- plotDF |>
drop_na(all_of(xVar)) |>
ggboxplot(x = xVar,
y = plotVar,
color = xVar,
title = plotVar) +
stat_compare_means(comparisons = comparisonList) +
scale_x_discrete(labels = function(x) paste0(x, "\n(n = ", table(plotDF[[xVar]])[x], ")")) +
rotate_x_text(90) +
xlab(NULL) +
theme(plot.title = element_text(size = 12, face = "bold"), legend.position = "none")
}
}
ggsave(filename = file.path(paste0(figures.dir, outfileName, ".pdf")),
plot = patchwork::wrap_plots(plotlist, ncol = numCol),
device = "pdf",
width = w * numCol,
height = h * ceiling(length(plotlist) / 4))
}
# ----------------------------------------------------
# generate figures
# ----------------------------------------------------
# FIG 1F
v.1f <- c("Endothelial", "Pericytes")
createBoxplots(plotVector = v.1f, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 2,
plotDF = luad.2018@meta.data,
outfileName = "1f")
# FIG 2E
v.2e <- c("p53_targets", "AT2", "CC.G2M", "`Glycolysis/Hypox`", "pEMT")
createBoxplots(plotVector = v.2e, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "2e")
# EXTENDED FIG 3A, D, E
v.s3ade <- c("p53_targets", "AT2", "CC.G2M", "`Glycolysis/Hypox`", "pEMT") # row 1
levels.tmp <- levels(luad.2018@meta.data$comut)
createBoxplots(plotVector = v.s3ade,
xVar = "comut",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x))[-6],
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "3a")
levels.tmp <- levels(luad.2018@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s3ade,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "3d")
levels.tmp <- levels(luad.2018@meta.data$impact)
createBoxplots(plotVector = v.s3ade,
xVar = "impact",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "3e")
# FIG 5D
v.5d <- c("SPP1")
createBoxplots(plotVector = v.5d, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "5d")
# EXTENDED FIG 1D
v.s1d <- c("B.cells", "Fibroblast", "Cancer", "Mast", "`T.CD4+`", # row 1
"Myeloid", "NK.cells", "Smooth.muscle", "T.cells", "`T.CD8+`") # row 2
createBoxplots(plotVector = v.s1d, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "s1d")
# EXTENDED FIG 3B
v.s3b <- c("p53_targets", "AT2", "CC.G2M", "Glycolysis.Hypox", "pEMT")
createBoxplots(plotVector = v.s3b,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 5,
plotDF = cptac.luad.prot@meta.data,
outfileName = "s3b")
# EXTENDED FIG 3C
v.s3c <- c("p53_targets", "AT2", "CC.G2M")
levels.tmp <- levels(cptac.luad.prot@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s3c, w = 3.5, h = 5,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 3,
plotDF = cptac.luad.prot@meta.data,
outfileName = "s3cA")
levels.tmp <- levels(cptac.luad.prot@meta.data$impact)
createBoxplots(plotVector = v.s3c, w = 3.5, h = 5,
xVar = "impact",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 3,
plotDF = cptac.luad.prot@meta.data,
outfileName = "s3cB")
# EXTENDED FIG 4A
v.s4a <- c("Endothelial")
levels.tmp <- levels(luad.2018@meta.data$comut)
createBoxplots(plotVector = v.s4a, w = 3, h = 5,
xVar = "comut",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x))[-6],
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s4aA")
levels.tmp <- levels(luad.2018@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s4a, w = 3, h = 5,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s4aB")
# EXTENDED FIG 4E
v.s4e <- c("Aerocyte", "Arterial")
createBoxplots(plotVector = v.s4e, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 2,
plotDF = luad.2018@meta.data,
outfileName = "s4e")
# EXTENDED FIG 5G
v.s5g <- c("CAF.ADH1B")
createBoxplots(plotVector = v.s5g, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s5g")
# EXTENDED FIG 5H
v.s5h <- c("Pericytes")
levels.tmp <- levels(luad.2018@meta.data$comut)
createBoxplots(plotVector = v.s5h, w = 3, h = 5,
xVar = "comut",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x))[-6],
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s5hA")
levels.tmp <- levels(luad.2018@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s5h, w = 3, h = 5,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s5hB")
# EXTENDED FIG 6B
v.s6b <- c("TAM.CXCL")
createBoxplots(plotVector = v.s6b, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s6b")
# EXTENDED FIG 6C
v.s6c <- c("CXCL9", "CXCL10", "CXCL11")
createBoxplots(plotVector = v.s6c, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 3,
plotDF = luad.2018@meta.data,
outfileName = "s6c")
# EXTENDED FIG 6D,E
v.s6de <- c("SPP1", "CXCL9", "CXCL10", "CXCL11")
levels.tmp <- levels(luad.2018@meta.data$comut)
createBoxplots(plotVector = v.s6de, w = 3, h = 5,
xVar = "comut",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x))[-6],
numCol = 4,
plotDF = luad.2018@meta.data,
outfileName = "s6d")
levels.tmp <- levels(luad.2018@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s6de, w = 3, h = 5,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 4,
plotDF = luad.2018@meta.data,
outfileName = "s6eA")
levels.tmp <- levels(luad.2018@meta.data$impact)
createBoxplots(plotVector = v.s6de, w = 3, h = 5,
xVar = "impact",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 4,
plotDF = luad.2018@meta.data,
outfileName = "s6eB")
# EXTENDED FIG 9D
v.s9d <- c("T.Exhausted", "CD8.GZMK")
createBoxplots(plotVector = v.s9d, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 2,
plotDF = luad.2018@meta.data,
outfileName = "s9d")
# EXTENDED FIG 7F
v.s7f <- c("CD274", "CD86", "PVR", "PDCD1", "CTLA4", "TIGIT")
createBoxplots(plotVector = v.s7f, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 6,
plotDF = luad.2018@meta.data,
outfileName = "s7f")
# EXTENDED FIG 7G
v.s7g <- c("PVR", "CD274")
createBoxplots(plotVector = v.s7g, w = 3, h = 4,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 2,
plotDF = cptac.luad.prot@meta.data,
outfileName = "s7g")
# ----------------------------------------------------
# smoking pack years
# ----------------------------------------------------
smoking_comparison_p53 <- read_csv(file.path(data.path, "smoking_comparison_p53.csv"))
p53.filt <- c("MGH1172", "MGH1173", "MGH1175", "MGH1182", "BWH16", "BWH11", "BWH14", "BWH06")
WT.filt <- c("P14", "MGH1174", "MGH1176", "MGH1183", "BWH01", "BWH04", "BWH05", "BWH09", "BWH19", "BWH23") # no MGH1170
xVar <- "p53_status"
plotVar <- "pack_yrs"
comparisonList = list(c("WT", "mut"))
plotDF <- smoking_comparison_p53 |>
# mutate(p53_status = case_when(Patient %in% c(p53.filt, "P17", "P19", "MGH1181") ~ "mut", Patient %in% c(WT.filt, "MGH1179") ~ "WT", .default = NA)) |>
mutate(p53_status = case_when(Patient %in% p53.filt ~ "mut", Patient %in% WT.filt ~ "WT", .default = NA)) |>
rename("pack_yrs" = "Pack yrs")
p.smoking_pack_yrs <- plotDF |>
drop_na(all_of(xVar)) |>
ggboxplot(x = xVar,
y = plotVar,
color = xVar,
palette = c("#0091CA", "#D8423D"),
add = "jitter",
title = plotVar) +
stat_compare_means(comparisons = comparisonList) +
scale_x_discrete(labels = function(x) paste0(x, "\n(n = ", table(plotDF[[xVar]])[x], ")")) +
rotate_x_text(0) +
xlab(NULL) +
ylab("Pack years") +
theme(plot.title = element_text(size = 12, face = "bold"), legend.position = "none")
p.smoking_pack_yrs
ggsave(filename = file.path(paste0(figures.dir, "smoking.pack_years.pdf")),
plot = p.smoking_pack_yrs,
device = "pdf",
width = 3,
height = 4)
# plotDF |>
# drop_na(all_of(xVar)) |>
# select(Patient, pack_yrs, p53_status) |>
# write.csv(file.path(data.path, "smoking.pack_years.csv"), row.names = F)
# # ----------------------------------------------------
figure_catalog/extended/bulk__3e/plot_code.r
levels.tmp <- levels(luad.2018@meta.data$impact)
createBoxplots(plotVector = v.s3ade,
xVar = "impact",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "3e")
figure_catalog/extended/bulk__3e/source_original.r
######################################################
#################### TSANKOV LAB #####################
######################################################
# ----------------------------------------------------
# load the libraries
# ----------------------------------------------------
library(Seurat)
library(tidyverse)
library(ggpubr)
# ----------------------------------------------------
# set directories
# ----------------------------------------------------
data.path <- "../data/"
figures.dir <- "../../figures/"
system(paste0("mkdir -p ", figures.dir))
# ----------------------------------------------------
# load the objects
# ----------------------------------------------------
luad.2018 <- readRDS(file.path(data.path, "luad.2018.processed.rds"))
cptac.luad.prot <- readRDS(file.path(data.path, "cptac.luad.prot.processed.rds"))
# ----------------------------------------------------
# load the dataframe
# ----------------------------------------------------
df.pvals.tcga <- read_csv(file.path(data.path, "tcga.regression_plots.p53_kras_egfr_stk11_CancerProp.filtered.csv")) |>
select(-"...1")
# ----------------------------------------------------
# custom function
# ----------------------------------------------------
createBoxplots <- function(plotVector,
plotDF,
outfileName,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
df.pvals = NULL,
numCol = 4,
w = 3.5,
h = 3) {
plotlist <- list()
for (plotVar in plotVector) {
plotVar.clean <- str_remove_all(string = plotVar, pattern = "`")
if (!plotVar.clean %in% colnames(plotDF)) {
message(paste0(plotVar.clean, " not found in plotDF!"))
next
}
message(paste0("Adding plot and saving csv file for ", plotVar, "..."))
saveDF <- plotDF |>
drop_na(all_of(xVar)) |>
select(all_of(c(xVar, plotVar.clean)))
# write.csv(saveDF, file.path(paste0(data.path, outfileName, ".", make.names(plotVar.clean), ".csv")))
if (length(levels(plotDF[[xVar]])) == 2) {
if (is.null(df.pvals)) {
# for cptac
plotlist[[plotVar]] <- plotDF |>
drop_na(all_of(xVar)) |>
ggboxplot(x = xVar,
y = plotVar,
color = xVar,
palette = c("#0091CA", "#D8423D"),
title = plotVar) +
stat_compare_means(comparisons = comparisonList) +
scale_x_discrete(labels = function(x) paste0(x, "\n(n = ", table(plotDF[[xVar]])[x], ")")) +
rotate_x_text(0) +
xlab(NULL) +
theme(plot.title = element_text(size = 12, face = "bold"), legend.position = "none")
} else {
# for tcga
pvals_vector <- setNames(df.pvals$pval, df.pvals$program)
pvals_manual <- pvals_vector[[plotVar.clean]]
plotlist[[plotVar]] <- plotDF |>
drop_na(all_of(xVar)) |>
ggboxplot(x = xVar,
y = plotVar,
color = xVar,
palette = c("#0091CA", "#D8423D"),
title = plotVar) +
# stat_compare_means(comparisons = comparisonList) +
scale_x_discrete(labels = function(x) paste0(x, "\n(n = ", table(plotDF[[xVar]])[x], ")")) +
rotate_x_text(0) +
xlab(NULL) +
theme(plot.title = element_text(size = 12, face = "bold"), legend.position = "none") +
ggsignif::geom_signif(comparisons = list(c("WT", "mut")), # Replace "WT" and "mut" with actual levels if different
annotations = paste0("p = ", signif(pvals_manual, digits = 3)),
y_position = max(plotDF[[plotVar.clean]], na.rm = TRUE) * 1.1, # Adjust as needed
tip_length = 0.02, # Adjust the length of the horizontal line tips
textsize = 3)
}
} else {
plotlist[[plotVar]] <- plotDF |>
drop_na(all_of(xVar)) |>
ggboxplot(x = xVar,
y = plotVar,
color = xVar,
title = plotVar) +
stat_compare_means(comparisons = comparisonList) +
scale_x_discrete(labels = function(x) paste0(x, "\n(n = ", table(plotDF[[xVar]])[x], ")")) +
rotate_x_text(90) +
xlab(NULL) +
theme(plot.title = element_text(size = 12, face = "bold"), legend.position = "none")
}
}
ggsave(filename = file.path(paste0(figures.dir, outfileName, ".pdf")),
plot = patchwork::wrap_plots(plotlist, ncol = numCol),
device = "pdf",
width = w * numCol,
height = h * ceiling(length(plotlist) / 4))
}
# ----------------------------------------------------
# generate figures
# ----------------------------------------------------
# FIG 1F
v.1f <- c("Endothelial", "Pericytes")
createBoxplots(plotVector = v.1f, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 2,
plotDF = luad.2018@meta.data,
outfileName = "1f")
# FIG 2E
v.2e <- c("p53_targets", "AT2", "CC.G2M", "`Glycolysis/Hypox`", "pEMT")
createBoxplots(plotVector = v.2e, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "2e")
# EXTENDED FIG 3A, D, E
v.s3ade <- c("p53_targets", "AT2", "CC.G2M", "`Glycolysis/Hypox`", "pEMT") # row 1
levels.tmp <- levels(luad.2018@meta.data$comut)
createBoxplots(plotVector = v.s3ade,
xVar = "comut",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x))[-6],
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "3a")
levels.tmp <- levels(luad.2018@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s3ade,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "3d")
levels.tmp <- levels(luad.2018@meta.data$impact)
createBoxplots(plotVector = v.s3ade,
xVar = "impact",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "3e")
# FIG 5D
v.5d <- c("SPP1")
createBoxplots(plotVector = v.5d, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "5d")
# EXTENDED FIG 1D
v.s1d <- c("B.cells", "Fibroblast", "Cancer", "Mast", "`T.CD4+`", # row 1
"Myeloid", "NK.cells", "Smooth.muscle", "T.cells", "`T.CD8+`") # row 2
createBoxplots(plotVector = v.s1d, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 5,
plotDF = luad.2018@meta.data,
outfileName = "s1d")
# EXTENDED FIG 3B
v.s3b <- c("p53_targets", "AT2", "CC.G2M", "Glycolysis.Hypox", "pEMT")
createBoxplots(plotVector = v.s3b,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 5,
plotDF = cptac.luad.prot@meta.data,
outfileName = "s3b")
# EXTENDED FIG 3C
v.s3c <- c("p53_targets", "AT2", "CC.G2M")
levels.tmp <- levels(cptac.luad.prot@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s3c, w = 3.5, h = 5,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 3,
plotDF = cptac.luad.prot@meta.data,
outfileName = "s3cA")
levels.tmp <- levels(cptac.luad.prot@meta.data$impact)
createBoxplots(plotVector = v.s3c, w = 3.5, h = 5,
xVar = "impact",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 3,
plotDF = cptac.luad.prot@meta.data,
outfileName = "s3cB")
# EXTENDED FIG 4A
v.s4a <- c("Endothelial")
levels.tmp <- levels(luad.2018@meta.data$comut)
createBoxplots(plotVector = v.s4a, w = 3, h = 5,
xVar = "comut",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x))[-6],
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s4aA")
levels.tmp <- levels(luad.2018@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s4a, w = 3, h = 5,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s4aB")
# EXTENDED FIG 4E
v.s4e <- c("Aerocyte", "Arterial")
createBoxplots(plotVector = v.s4e, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 2,
plotDF = luad.2018@meta.data,
outfileName = "s4e")
# EXTENDED FIG 5G
v.s5g <- c("CAF.ADH1B")
createBoxplots(plotVector = v.s5g, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s5g")
# EXTENDED FIG 5H
v.s5h <- c("Pericytes")
levels.tmp <- levels(luad.2018@meta.data$comut)
createBoxplots(plotVector = v.s5h, w = 3, h = 5,
xVar = "comut",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x))[-6],
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s5hA")
levels.tmp <- levels(luad.2018@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s5h, w = 3, h = 5,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s5hB")
# EXTENDED FIG 6B
v.s6b <- c("TAM.CXCL")
createBoxplots(plotVector = v.s6b, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 1,
plotDF = luad.2018@meta.data,
outfileName = "s6b")
# EXTENDED FIG 6C
v.s6c <- c("CXCL9", "CXCL10", "CXCL11")
createBoxplots(plotVector = v.s6c, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 3,
plotDF = luad.2018@meta.data,
outfileName = "s6c")
# EXTENDED FIG 6D,E
v.s6de <- c("SPP1", "CXCL9", "CXCL10", "CXCL11")
levels.tmp <- levels(luad.2018@meta.data$comut)
createBoxplots(plotVector = v.s6de, w = 3, h = 5,
xVar = "comut",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x))[-6],
numCol = 4,
plotDF = luad.2018@meta.data,
outfileName = "s6d")
levels.tmp <- levels(luad.2018@meta.data$DNE_LOFclass)
createBoxplots(plotVector = v.s6de, w = 3, h = 5,
xVar = "DNE_LOFclass",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 4,
plotDF = luad.2018@meta.data,
outfileName = "s6eA")
levels.tmp <- levels(luad.2018@meta.data$impact)
createBoxplots(plotVector = v.s6de, w = 3, h = 5,
xVar = "impact",
comparisonList = lapply(setdiff(levels.tmp, "WT"), function(x) c("WT", x)),
numCol = 4,
plotDF = luad.2018@meta.data,
outfileName = "s6eB")
# EXTENDED FIG 9D
v.s9d <- c("T.Exhausted", "CD8.GZMK")
createBoxplots(plotVector = v.s9d, w = 3, h = 4, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 2,
plotDF = luad.2018@meta.data,
outfileName = "s9d")
# EXTENDED FIG 7F
v.s7f <- c("CD274", "CD86", "PVR", "PDCD1", "CTLA4", "TIGIT")
createBoxplots(plotVector = v.s7f, df.pvals = df.pvals.tcga,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 6,
plotDF = luad.2018@meta.data,
outfileName = "s7f")
# EXTENDED FIG 7G
v.s7g <- c("PVR", "CD274")
createBoxplots(plotVector = v.s7g, w = 3, h = 4,
xVar = "p53_status",
comparisonList = list(c("WT", "mut")),
numCol = 2,
plotDF = cptac.luad.prot@meta.data,
outfileName = "s7g")
# ----------------------------------------------------
# smoking pack years
# ----------------------------------------------------
smoking_comparison_p53 <- read_csv(file.path(data.path, "smoking_comparison_p53.csv"))
p53.filt <- c("MGH1172", "MGH1173", "MGH1175", "MGH1182", "BWH16", "BWH11", "BWH14", "BWH06")
WT.filt <- c("P14", "MGH1174", "MGH1176", "MGH1183", "BWH01", "BWH04", "BWH05", "BWH09", "BWH19", "BWH23") # no MGH1170
xVar <- "p53_status"
plotVar <- "pack_yrs"
comparisonList = list(c("WT", "mut"))
plotDF <- smoking_comparison_p53 |>
# mutate(p53_status = case_when(Patient %in% c(p53.filt, "P17", "P19", "MGH1181") ~ "mut", Patient %in% c(WT.filt, "MGH1179") ~ "WT", .default = NA)) |>
mutate(p53_status = case_when(Patient %in% p53.filt ~ "mut", Patient %in% WT.filt ~ "WT", .default = NA)) |>
rename("pack_yrs" = "Pack yrs")
p.smoking_pack_yrs <- plotDF |>
drop_na(all_of(xVar)) |>
ggboxplot(x = xVar,
y = plotVar,
color = xVar,
palette = c("#0091CA", "#D8423D"),
add = "jitter",
title = plotVar) +
stat_compare_means(comparisons = comparisonList) +
scale_x_discrete(labels = function(x) paste0(x, "\n(n = ", table(plotDF[[xVar]])[x], ")")) +
rotate_x_text(0) +
xlab(NULL) +
ylab("Pack years") +
theme(plot.title = element_text(size = 12, face = "bold"), legend.position = "none")
p.smoking_pack_yrs
ggsave(filename = file.path(paste0(figures.dir, "smoking.pack_years.pdf")),
plot = p.smoking_pack_yrs,
device = "pdf",
width = 3,
height = 4)
# plotDF |>
# drop_na(all_of(xVar)) |>
# select(Patient, pack_yrs, p53_status) |>
# write.csv(file.path(data.path, "smoking.pack_years.csv"), row.names = F)
# # ----------------------------------------------------
figure_catalog/extended/bulk__ext_fig8b_left/plot_code.r
##################### Exteneded 8B (left)
##################### adeno + all treatments
se = which(clin$Histology_Harmonized == "Adeno") ## adeno + all treatments
clin.subset = clin[se,]
### mut of common genes
clin.subset1 = clin.subset[,c("EGFR_status", "STK11_status", "KEAP1_status", "KRAS_status",
"TP53_status", "OS_Years", "PFS_Years", "OS_Status", "PFS_Status")]
sub = colnames(clin.subset1)[1:5]
survdiff.pvalOS = coxph.pvalOS = survdiff.pvalPFS = coxph.pvalPFS = rep(0, length(sub))
hrOS = hrPFS = lbOS = ubOS = lbPFS = ubPFS = num_high = num_low = rep(0, length(sub))
for (k in 1:length(sub))
{
cat("\r", k)
label = rep(NA, nrow(clin.subset1))
se = which(clin.subset1[,k] == "Wildtype")
label[se] = 0
se = which(clin.subset1[,k] == "Mutants")
label[se] = 1
xx = cbind(clin.subset1, label)
se = !is.na(xx$label)
xx = xx[se,]
a = table(label)
num_high[k] = as.numeric(a[names(a)==1])
num_low[k] = as.numeric(a[names(a)==0])
mycox = coxph(Surv(OS_Years, OS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalOS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrOS[k] = tmp[1]
lbOS[k] = tmp[3]
ubOS[k] = tmp[4]
mycox = coxph(Surv(PFS_Years, PFS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalPFS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrPFS[k] = tmp[1]
lbPFS[k] = tmp[3]
ubPFS[k] = tmp[4]
}
res.surv = data.frame(sub, num_high, num_low, coxph.pvalOS, hrOS, lbOS, ubOS, coxph.pvalPFS, hrPFS, lbPFS, ubPFS)
res.surv1 = res.surv
### mut of p53 specific subsets
se.list = list()
se.list[[1]] = which(clin.subset$Impact1 == 1)
se.list[[2]] = which(clin.subset$Impact2 == 1)
se.list[[3]] = which(clin.subset$DNE_LOFclass == "DNE_LOF")
se.list[[4]] = which(clin.subset$DNE_LOFclass == "notDNE_LOF")
se.list[[5]] = which(clin.subset$DNE_LOFclass == "notDNE_notLOF")
se.list[[6]] = grep("Missense_Mutation", clin.subset$TP53.all)
se.list[[7]] = grep("Nonsense_Mutation", clin.subset$TP53.all)
se.list[[8]] = grep("Splice_Site", clin.subset$TP53.all)
se.list[[9]] = which(clin.subset$TP53.all == "Splice_Site")
se.list[[10]] = grep("Frame_Shift_Ins|Frame_Shift_Del", clin.subset$TP53.all)
sub = c("Impact1", "Impact2", "DNE_LOF", "notDNE_LOF", "notDNE_notLOF",
"Missense", "Nonsense", "Splice_Site_All", "Splice_Site_Unique", "Frameshift")
survdiff.pvalOS = coxph.pvalOS = survdiff.pvalPFS = coxph.pvalPFS = rep(0, length(sub))
hrOS = hrPFS = lbOS = ubOS = lbPFS = ubPFS = num_high = num_low = rep(0, length(sub))
for (k in 1:length(sub))
{
cat("\r", k)
label = rep(NA, nrow(clin.subset))
se = which(clin.subset$TP53.all == "Wildtype")
label[se] = 0
se = se.list[[k]]
label[se] = 1
xx = cbind(clin.subset, label)
se = !is.na(xx$label)
xx = xx[se,]
a = table(label)
num_high[k] = as.numeric(a[names(a)==1])
num_low[k] = as.numeric(a[names(a)==0])
mycox = coxph(Surv(OS_Years, OS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalOS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrOS[k] = tmp[1]
lbOS[k] = tmp[3]
ubOS[k] = tmp[4]
mycox = coxph(Surv(PFS_Years, PFS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalPFS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrPFS[k] = tmp[1]
lbPFS[k] = tmp[3]
ubPFS[k] = tmp[4]
}
res.surv = data.frame(sub, num_high, num_low, coxph.pvalOS, hrOS, lbOS, ubOS, coxph.pvalPFS, hrPFS, lbPFS, ubPFS)
res.surv2 = res.surv
###
res.surv = rbind(res.surv1, res.surv2)
##
coxph.pvalPFS = res.surv$coxph.pvalPFS
lbPFS = res.surv$lbPFS
hrPFS = res.surv$hrPFS
ubPFS = res.surv$ubPFS
sub = res.surv$sub
num_high = res.surv$num_high
num_low = res.surv$num_low
coxph.pvalPFS.ann = NULL
for (i in 1:length(coxph.pvalPFS))
{
tmp = coxph.pvalPFS[i]
coxph.pvalPFS.ann[i] = as.character(cut(tmp,
breaks = c(-Inf, 0.001, 0.1, Inf),
labels = c(paste("p =", formatC(tmp, 1, format = "e")),
paste("p =", formatC(signif(tmp, 2), 2, flag = "#")),
"p > 0.1")))
}
sub2 = c("EGFR", "STK11", "KEAP1", "KRAS", "TP53", " Impactful 1", " Impactful 2", " DNE,LOF", " notDNE,LOF", " notDNE,notLOF",
" Missense", " Nonsense", " Splice Site", " Unique Splice Site", " Frameshift")
sub2 = paste0(sub2, " (", num_high, ")")
res1 = data.frame(lbPFS, hrPFS, ubPFS)
res2 = data.frame(sub2, coxph.pvalPFS.ann)
##############
pdf(paste0(figures.dir, "Ext_Fig8B_left.pdf"), useDingbats = F, 9, 9)
forestplot(res2, as.matrix(res1), graph.pos = 2, new_page = F,
xlab = "Hazard ratio",
xticks = c(0, 1, 2, 3, 4),
graphwidth = unit(90, "mm"),
vertices = TRUE,
col = fpColors(lines = "black"),
boxsize = 0.3,
txt_gp = fpTxtGp(xlab = gpar(cex = 2.25),
ticks = gpar(cex = 1.75),
title = gpar(cex = 2.5, fontface = "bold"),
cex = 1.5),
align = c("l", "l"),
zero = 1)
dev.off()
figure_catalog/extended/bulk__ext_fig8b_left/source_original.r
#####################@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#####################@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
##################### FIGURE 7E AND EXTENDED 8B
##################### SURVIVAL ANALYSES OF SU2C-MARK COHORT
library(readxl)
library(survival)
library(forestplot)
# my outputs
data.path <- "../data/"
figures.dir <- "../../figures/"
#####################
#####################
clin = read_excel(paste0(data.path, "su2c_clinical_info.xlsx"))
clin = as.data.frame(clin)
clin = clin[-1,]
colnames(clin) = clin[1,]
clin = clin[-1,]
clin$OS_Years = as.numeric(clin$Harmonized_OS_Days) / 365
clin$PFS_Years = as.numeric(clin$Harmonized_PFS_Days) / 365
clin$OS_Status = as.numeric(clin$Harmonized_OS_Event)
clin$PFS_Status = as.numeric(clin$Harmonized_PFS_Event)
clin.bkup1 = clin
se = which(clin$WES_All == 1)
clin = clin[se,]
rownames(clin) = clin$Harmonized_SU2C_WES_Tumor_Sample_ID_v2
########
p53 = get(load(paste0(data.path, "su2c.mut.info.p53.Rda")))
mymut = get(load(paste0(data.path, "su2c.mut.info.selectgenes1.Rda")))
old.clin = clin
all(rownames(clin) == rownames(p53))
all(rownames(clin) == rownames(mymut))
clin = cbind(clin, mymut[,-c(3,4)], p53)
######## convert silent mutations to WT
se = which(clin$EGFR.all == "Silent")
clin$EGFR_status[se] = "Wildtype"
se = which(clin$KEAP1.all == "Silent")
clin$KEAP1_status[se] = "Wildtype"
se = which(clin$KRAS.all == "Silent")
clin$KRAS_status[se] = "Wildtype"
se = which(clin$STK11.all == "Silent") ## 0
se = which(clin$TP53.all == "Silent") ## 0
##################### 7E (left)
##################### adeno + anti-pd1
se = which(clin$Agent_PD1 %in% c("Nivolumab", "Pembrolizumab") & clin$Histology_Harmonized == "Adeno") ## adeno+ anti-pd-1
clin.subset = clin[se,]
### mut of common genes
clin.subset1 = clin.subset[,c("EGFR_status", "STK11_status", "KEAP1_status", "KRAS_status",
"TP53_status", "OS_Years", "PFS_Years", "OS_Status", "PFS_Status")]
sub = colnames(clin.subset1)[1:5]
survdiff.pvalOS = coxph.pvalOS = survdiff.pvalPFS = coxph.pvalPFS = rep(0, length(sub))
hrOS = hrPFS = lbOS = ubOS = lbPFS = ubPFS = num_high = num_low = rep(0, length(sub))
for (k in 1:length(sub))
{
cat("\r", k)
label = rep(NA, nrow(clin.subset1))
se = which(clin.subset1[,k] == "Wildtype")
label[se] = 0
se = which(clin.subset1[,k] == "Mutants")
label[se] = 1
xx = cbind(clin.subset1, label)
se = !is.na(xx$label)
xx = xx[se,]
a = table(label)
num_high[k] = as.numeric(a[names(a)==1])
num_low[k] = as.numeric(a[names(a)==0])
mycox = coxph(Surv(OS_Years, OS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalOS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrOS[k] = tmp[1]
lbOS[k] = tmp[3]
ubOS[k] = tmp[4]
mycox = coxph(Surv(PFS_Years, PFS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalPFS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrPFS[k] = tmp[1]
lbPFS[k] = tmp[3]
ubPFS[k] = tmp[4]
}
res.surv = data.frame(sub, num_high, num_low, coxph.pvalOS, hrOS, lbOS, ubOS, coxph.pvalPFS, hrPFS, lbPFS, ubPFS)
res.surv1 = res.surv
### mut of p53 specific subsets
se.list = list()
se.list[[1]] = which(clin.subset$Impact1 == 1)
se.list[[2]] = which(clin.subset$Impact2 == 1)
se.list[[3]] = which(clin.subset$DNE_LOFclass == "DNE_LOF")
se.list[[4]] = which(clin.subset$DNE_LOFclass == "notDNE_LOF")
se.list[[5]] = which(clin.subset$DNE_LOFclass == "notDNE_notLOF")
se.list[[6]] = grep("Missense_Mutation", clin.subset$TP53.all)
se.list[[7]] = grep("Nonsense_Mutation", clin.subset$TP53.all)
se.list[[8]] = grep("Splice_Site", clin.subset$TP53.all)
se.list[[9]] = which(clin.subset$TP53.all == "Splice_Site")
se.list[[10]] = grep("Frame_Shift_Ins|Frame_Shift_Del", clin.subset$TP53.all)
sub = c("Impact1", "Impact2", "DNE_LOF", "notDNE_LOF", "notDNE_notLOF",
"Missense", "Nonsense", "Splice_Site_All", "Splice_Site_Unique", "Frameshift")
survdiff.pvalOS = coxph.pvalOS = survdiff.pvalPFS = coxph.pvalPFS = rep(0, length(sub))
hrOS = hrPFS = lbOS = ubOS = lbPFS = ubPFS = num_high = num_low = rep(0, length(sub))
for (k in 1:length(sub))
{
cat("\r", k)
label = rep(NA, nrow(clin.subset))
se = which(clin.subset$TP53.all == "Wildtype")
label[se] = 0
se = se.list[[k]]
label[se] = 1
xx = cbind(clin.subset, label)
se = !is.na(xx$label)
xx = xx[se,]
a = table(label)
num_high[k] = as.numeric(a[names(a)==1])
num_low[k] = as.numeric(a[names(a)==0])
mycox = coxph(Surv(OS_Years, OS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalOS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrOS[k] = tmp[1]
lbOS[k] = tmp[3]
ubOS[k] = tmp[4]
mycox = coxph(Surv(PFS_Years, PFS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalPFS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrPFS[k] = tmp[1]
lbPFS[k] = tmp[3]
ubPFS[k] = tmp[4]
}
res.surv = data.frame(sub, num_high, num_low, coxph.pvalOS, hrOS, lbOS, ubOS, coxph.pvalPFS, hrPFS, lbPFS, ubPFS)
res.surv2 = res.surv
###
res.surv = rbind(res.surv1, res.surv2)
##
coxph.pvalPFS = res.surv$coxph.pvalPFS
lbPFS = res.surv$lbPFS
hrPFS = res.surv$hrPFS
ubPFS = res.surv$ubPFS
sub = res.surv$sub
num_high = res.surv$num_high
num_low = res.surv$num_low
coxph.pvalPFS.ann = NULL
for (i in 1:length(coxph.pvalPFS))
{
tmp = coxph.pvalPFS[i]
coxph.pvalPFS.ann[i] = as.character(cut(tmp,
breaks = c(-Inf, 0.001, 0.1, Inf),
labels = c(paste("p =", formatC(tmp, 1, format = "e")),
paste("p =", formatC(signif(tmp, 2), 2, flag = "#")),
"p > 0.1")))
}
sub2 = c("EGFR", "STK11", "KEAP1", "KRAS", "TP53", " Impactful 1", " Impactful 2", " DNE,LOF", " notDNE,LOF", " notDNE,notLOF",
" Missense", " Nonsense", " Splice Site", " Unique Splice Site", " Frameshift")
sub2 = paste0(sub2, " (", num_high, ")")
res1 = data.frame(lbPFS, hrPFS, ubPFS)
res2 = data.frame(sub2, coxph.pvalPFS.ann)
##############
pdf(paste0(figures.dir, "Fig7E_left.pdf"), useDingbats = F, 9, 9)
forestplot(res2, as.matrix(res1), graph.pos = 2, new_page = F,
xlab = "Hazard ratio",
xticks = c(0, 1, 2, 3, 4),
graphwidth = unit(90, "mm"),
vertices = TRUE,
col = fpColors(lines = "black"),
boxsize = 0.3,
txt_gp = fpTxtGp(xlab = gpar(cex = 2.25),
ticks = gpar(cex = 1.75),
title = gpar(cex = 2.5, fontface = "bold"),
cex = 1.5),
align = c("l", "l"),
zero = 1)
dev.off()
##################### Exteneded 8B (left)
##################### adeno + all treatments
se = which(clin$Histology_Harmonized == "Adeno") ## adeno + all treatments
clin.subset = clin[se,]
### mut of common genes
clin.subset1 = clin.subset[,c("EGFR_status", "STK11_status", "KEAP1_status", "KRAS_status",
"TP53_status", "OS_Years", "PFS_Years", "OS_Status", "PFS_Status")]
sub = colnames(clin.subset1)[1:5]
survdiff.pvalOS = coxph.pvalOS = survdiff.pvalPFS = coxph.pvalPFS = rep(0, length(sub))
hrOS = hrPFS = lbOS = ubOS = lbPFS = ubPFS = num_high = num_low = rep(0, length(sub))
for (k in 1:length(sub))
{
cat("\r", k)
label = rep(NA, nrow(clin.subset1))
se = which(clin.subset1[,k] == "Wildtype")
label[se] = 0
se = which(clin.subset1[,k] == "Mutants")
label[se] = 1
xx = cbind(clin.subset1, label)
se = !is.na(xx$label)
xx = xx[se,]
a = table(label)
num_high[k] = as.numeric(a[names(a)==1])
num_low[k] = as.numeric(a[names(a)==0])
mycox = coxph(Surv(OS_Years, OS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalOS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrOS[k] = tmp[1]
lbOS[k] = tmp[3]
ubOS[k] = tmp[4]
mycox = coxph(Surv(PFS_Years, PFS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalPFS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrPFS[k] = tmp[1]
lbPFS[k] = tmp[3]
ubPFS[k] = tmp[4]
}
res.surv = data.frame(sub, num_high, num_low, coxph.pvalOS, hrOS, lbOS, ubOS, coxph.pvalPFS, hrPFS, lbPFS, ubPFS)
res.surv1 = res.surv
### mut of p53 specific subsets
se.list = list()
se.list[[1]] = which(clin.subset$Impact1 == 1)
se.list[[2]] = which(clin.subset$Impact2 == 1)
se.list[[3]] = which(clin.subset$DNE_LOFclass == "DNE_LOF")
se.list[[4]] = which(clin.subset$DNE_LOFclass == "notDNE_LOF")
se.list[[5]] = which(clin.subset$DNE_LOFclass == "notDNE_notLOF")
se.list[[6]] = grep("Missense_Mutation", clin.subset$TP53.all)
se.list[[7]] = grep("Nonsense_Mutation", clin.subset$TP53.all)
se.list[[8]] = grep("Splice_Site", clin.subset$TP53.all)
se.list[[9]] = which(clin.subset$TP53.all == "Splice_Site")
se.list[[10]] = grep("Frame_Shift_Ins|Frame_Shift_Del", clin.subset$TP53.all)
sub = c("Impact1", "Impact2", "DNE_LOF", "notDNE_LOF", "notDNE_notLOF",
"Missense", "Nonsense", "Splice_Site_All", "Splice_Site_Unique", "Frameshift")
survdiff.pvalOS = coxph.pvalOS = survdiff.pvalPFS = coxph.pvalPFS = rep(0, length(sub))
hrOS = hrPFS = lbOS = ubOS = lbPFS = ubPFS = num_high = num_low = rep(0, length(sub))
for (k in 1:length(sub))
{
cat("\r", k)
label = rep(NA, nrow(clin.subset))
se = which(clin.subset$TP53.all == "Wildtype")
label[se] = 0
se = se.list[[k]]
label[se] = 1
xx = cbind(clin.subset, label)
se = !is.na(xx$label)
xx = xx[se,]
a = table(label)
num_high[k] = as.numeric(a[names(a)==1])
num_low[k] = as.numeric(a[names(a)==0])
mycox = coxph(Surv(OS_Years, OS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalOS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrOS[k] = tmp[1]
lbOS[k] = tmp[3]
ubOS[k] = tmp[4]
mycox = coxph(Surv(PFS_Years, PFS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalPFS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrPFS[k] = tmp[1]
lbPFS[k] = tmp[3]
ubPFS[k] = tmp[4]
}
res.surv = data.frame(sub, num_high, num_low, coxph.pvalOS, hrOS, lbOS, ubOS, coxph.pvalPFS, hrPFS, lbPFS, ubPFS)
res.surv2 = res.surv
###
res.surv = rbind(res.surv1, res.surv2)
##
coxph.pvalPFS = res.surv$coxph.pvalPFS
lbPFS = res.surv$lbPFS
hrPFS = res.surv$hrPFS
ubPFS = res.surv$ubPFS
sub = res.surv$sub
num_high = res.surv$num_high
num_low = res.surv$num_low
coxph.pvalPFS.ann = NULL
for (i in 1:length(coxph.pvalPFS))
{
tmp = coxph.pvalPFS[i]
coxph.pvalPFS.ann[i] = as.character(cut(tmp,
breaks = c(-Inf, 0.001, 0.1, Inf),
labels = c(paste("p =", formatC(tmp, 1, format = "e")),
paste("p =", formatC(signif(tmp, 2), 2, flag = "#")),
"p > 0.1")))
}
sub2 = c("EGFR", "STK11", "KEAP1", "KRAS", "TP53", " Impactful 1", " Impactful 2", " DNE,LOF", " notDNE,LOF", " notDNE,notLOF",
" Missense", " Nonsense", " Splice Site", " Unique Splice Site", " Frameshift")
sub2 = paste0(sub2, " (", num_high, ")")
res1 = data.frame(lbPFS, hrPFS, ubPFS)
res2 = data.frame(sub2, coxph.pvalPFS.ann)
##############
pdf(paste0(figures.dir, "Ext_Fig8B_left.pdf"), useDingbats = F, 9, 9)
forestplot(res2, as.matrix(res1), graph.pos = 2, new_page = F,
xlab = "Hazard ratio",
xticks = c(0, 1, 2, 3, 4),
graphwidth = unit(90, "mm"),
vertices = TRUE,
col = fpColors(lines = "black"),
boxsize = 0.3,
txt_gp = fpTxtGp(xlab = gpar(cex = 2.25),
ticks = gpar(cex = 1.75),
title = gpar(cex = 2.5, fontface = "bold"),
cex = 1.5),
align = c("l", "l"),
zero = 1)
dev.off()
##################### 7E (right)
##################### squamous + anti-pd1
se = which(clin$Agent_PD1 %in% c("Nivolumab", "Pembrolizumab") & clin$Histology_Harmonized == "Squamous") ## sq+ anti-pd-1
clin.subset = clin[se,]
### mut of common genes
clin.subset1 = clin.subset[,c("EGFR_status", "STK11_status", "KEAP1_status", "KRAS_status",
"TP53_status", "OS_Years", "PFS_Years", "OS_Status", "PFS_Status")]
sub = colnames(clin.subset1)[1:5]
survdiff.pvalOS = coxph.pvalOS = survdiff.pvalPFS = coxph.pvalPFS = rep(0, length(sub))
hrOS = hrPFS = lbOS = ubOS = lbPFS = ubPFS = num_high = num_low = rep(0, length(sub))
for (k in 2:length(sub)) ## egfr has no muts (pdl1), only 1 mut (all treatments)
{
cat("\r", k)
label = rep(NA, nrow(clin.subset1))
se = which(clin.subset1[,k] == "Wildtype")
label[se] = 0
se = which(clin.subset1[,k] == "Mutants")
label[se] = 1
xx = cbind(clin.subset1, label)
se = !is.na(xx$label)
xx = xx[se,]
a = table(label)
num_high[k] = as.numeric(a[names(a)==1])
num_low[k] = as.numeric(a[names(a)==0])
mycox = coxph(Surv(OS_Years, OS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalOS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrOS[k] = tmp[1]
lbOS[k] = tmp[3]
ubOS[k] = tmp[4]
mycox = coxph(Surv(PFS_Years, PFS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalPFS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrPFS[k] = tmp[1]
lbPFS[k] = tmp[3]
ubPFS[k] = tmp[4]
}
res.surv = data.frame(sub, num_high, num_low, coxph.pvalOS, hrOS, lbOS, ubOS, coxph.pvalPFS, hrPFS, lbPFS, ubPFS)
res.surv1 = res.surv
### mut of p53 specific subsets
se.list = list()
se.list[[1]] = which(clin.subset$Impact1 == 1)
se.list[[2]] = which(clin.subset$Impact2 == 1)
se.list[[3]] = which(clin.subset$DNE_LOFclass == "DNE_LOF")
se.list[[4]] = which(clin.subset$DNE_LOFclass == "notDNE_LOF")
se.list[[5]] = which(clin.subset$DNE_LOFclass == "notDNE_notLOF")
se.list[[6]] = grep("Missense_Mutation", clin.subset$TP53.all)
se.list[[7]] = grep("Nonsense_Mutation", clin.subset$TP53.all)
se.list[[8]] = grep("Splice_Site", clin.subset$TP53.all)
se.list[[9]] = which(clin.subset$TP53.all == "Splice_Site")
se.list[[10]] = grep("Frame_Shift_Ins|Frame_Shift_Del", clin.subset$TP53.all)
sub = c("Impact1", "Impact2", "DNE_LOF", "notDNE_LOF", "notDNE_notLOF",
"Missense", "Nonsense", "Splice_Site_All", "Splice_Site_Unique", "Frameshift")
survdiff.pvalOS = coxph.pvalOS = survdiff.pvalPFS = coxph.pvalPFS = rep(0, length(sub))
hrOS = hrPFS = lbOS = ubOS = lbPFS = ubPFS = num_high = num_low = rep(0, length(sub))
for (k in 1:length(sub))
{
cat("\r", k)
label = rep(NA, nrow(clin.subset))
se = which(clin.subset$TP53.all == "Wildtype")
label[se] = 0
se = se.list[[k]]
label[se] = 1
xx = cbind(clin.subset, label)
se = !is.na(xx$label)
xx = xx[se,]
a = table(label)
num_high[k] = as.numeric(a[names(a)==1])
num_low[k] = as.numeric(a[names(a)==0])
mycox = coxph(Surv(OS_Years, OS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalOS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrOS[k] = tmp[1]
lbOS[k] = tmp[3]
ubOS[k] = tmp[4]
mycox = coxph(Surv(PFS_Years, PFS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalPFS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrPFS[k] = tmp[1]
lbPFS[k] = tmp[3]
ubPFS[k] = tmp[4]
}
res.surv = data.frame(sub, num_high, num_low, coxph.pvalOS, hrOS, lbOS, ubOS, coxph.pvalPFS, hrPFS, lbPFS, ubPFS)
res.surv2 = res.surv
###
res.surv = rbind(res.surv1, res.surv2)
##
coxph.pvalPFS = res.surv$coxph.pvalPFS
lbPFS = res.surv$lbPFS
hrPFS = res.surv$hrPFS
ubPFS = res.surv$ubPFS
sub = res.surv$sub
num_high = res.surv$num_high
num_low = res.surv$num_low
coxph.pvalPFS.ann = NULL
for (i in 1:length(coxph.pvalPFS))
{
tmp = coxph.pvalPFS[i]
coxph.pvalPFS.ann[i] = as.character(cut(tmp,
breaks = c(-Inf, 0.001, 0.1, Inf),
labels = c(paste("p =", formatC(tmp, 1, format = "e")),
paste("p =", formatC(signif(tmp, 2), 2, flag = "#")),
"p > 0.1")))
}
sub2 = c("EGFR", "STK11", "KEAP1", "KRAS", "TP53", " Impactful 1", " Impactful 2", " DNE,LOF", " notDNE,LOF", " notDNE,notLOF",
" Missense", " Nonsense", " Splice Site", " Unique Splice Site", " Frameshift")
sub2 = paste0(sub2, " (", num_high, ")")
res1 = data.frame(lbPFS, hrPFS, ubPFS)
res2 = data.frame(sub2, coxph.pvalPFS.ann)
##
res1 = res1[-1,]
res2 = res2[-1,]
## set range hr >17 -> 17
res1.mod = res1
se = which(res1.mod[,"ubPFS"] > 17)
res1.mod[se,"ubPFS"] = 17
##############
pdf(paste0(figures.dir, "Fig7E_right.pdf"), useDingbats = F, 9, 9)
forestplot(res2, as.matrix(res1.mod), graph.pos = 2, new_page = F,
xlab = "Hazard ratio",
xticks = c(0, 1, 5, 10, 15, 17),
graphwidth = unit(90, "mm"),
vertices = TRUE,
col = fpColors(lines = "black"),
boxsize = 0.3,
txt_gp = fpTxtGp(xlab = gpar(cex = 2.25),
ticks = gpar(cex = 1.75),
title = gpar(cex = 2.5, fontface = "bold"),
cex = 1.5),
align = c("l", "l"),
zero = 1)
dev.off()
##################### Exteneded 8B (right))
#####################
se = which(clin$Histology_Harmonized == "Squamous") ## sq + all treatments
clin.subset = clin[se,]
### mut of common genes
clin.subset1 = clin.subset[,c("EGFR_status", "STK11_status", "KEAP1_status", "KRAS_status",
"TP53_status", "OS_Years", "PFS_Years", "OS_Status", "PFS_Status")]
sub = colnames(clin.subset1)[1:5]
survdiff.pvalOS = coxph.pvalOS = survdiff.pvalPFS = coxph.pvalPFS = rep(0, length(sub))
hrOS = hrPFS = lbOS = ubOS = lbPFS = ubPFS = num_high = num_low = rep(0, length(sub))
for (k in 2:length(sub)) ## egfr has no muts (pdl1), only 1 mut (all treatments)
{
cat("\r", k)
label = rep(NA, nrow(clin.subset1))
se = which(clin.subset1[,k] == "Wildtype")
label[se] = 0
se = which(clin.subset1[,k] == "Mutants")
label[se] = 1
xx = cbind(clin.subset1, label)
se = !is.na(xx$label)
xx = xx[se,]
a = table(label)
num_high[k] = as.numeric(a[names(a)==1])
num_low[k] = as.numeric(a[names(a)==0])
mycox = coxph(Surv(OS_Years, OS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalOS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrOS[k] = tmp[1]
lbOS[k] = tmp[3]
ubOS[k] = tmp[4]
mycox = coxph(Surv(PFS_Years, PFS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalPFS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrPFS[k] = tmp[1]
lbPFS[k] = tmp[3]
ubPFS[k] = tmp[4]
}
res.surv = data.frame(sub, num_high, num_low, coxph.pvalOS, hrOS, lbOS, ubOS, coxph.pvalPFS, hrPFS, lbPFS, ubPFS)
res.surv1 = res.surv
### mut of p53 specific subsets
se.list = list()
se.list[[1]] = which(clin.subset$Impact1 == 1)
se.list[[2]] = which(clin.subset$Impact2 == 1)
se.list[[3]] = which(clin.subset$DNE_LOFclass == "DNE_LOF")
se.list[[4]] = which(clin.subset$DNE_LOFclass == "notDNE_LOF")
se.list[[5]] = which(clin.subset$DNE_LOFclass == "notDNE_notLOF")
se.list[[6]] = grep("Missense_Mutation", clin.subset$TP53.all)
se.list[[7]] = grep("Nonsense_Mutation", clin.subset$TP53.all)
se.list[[8]] = grep("Splice_Site", clin.subset$TP53.all)
se.list[[9]] = which(clin.subset$TP53.all == "Splice_Site")
se.list[[10]] = grep("Frame_Shift_Ins|Frame_Shift_Del", clin.subset$TP53.all)
sub = c("Impact1", "Impact2", "DNE_LOF", "notDNE_LOF", "notDNE_notLOF",
"Missense", "Nonsense", "Splice_Site_All", "Splice_Site_Unique", "Frameshift")
survdiff.pvalOS = coxph.pvalOS = survdiff.pvalPFS = coxph.pvalPFS = rep(0, length(sub))
hrOS = hrPFS = lbOS = ubOS = lbPFS = ubPFS = num_high = num_low = rep(0, length(sub))
for (k in 1:length(sub))
{
cat("\r", k)
label = rep(NA, nrow(clin.subset))
se = which(clin.subset$TP53.all == "Wildtype")
label[se] = 0
se = se.list[[k]]
label[se] = 1
xx = cbind(clin.subset, label)
se = !is.na(xx$label)
xx = xx[se,]
a = table(label)
num_high[k] = as.numeric(a[names(a)==1])
num_low[k] = as.numeric(a[names(a)==0])
mycox = coxph(Surv(OS_Years, OS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalOS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrOS[k] = tmp[1]
lbOS[k] = tmp[3]
ubOS[k] = tmp[4]
mycox = coxph(Surv(PFS_Years, PFS_Status)~label, xx)
mycox = summary(mycox)
coxph.pvalPFS[k] = mycox$coefficients[5]
tmp = mycox$conf.int
hrPFS[k] = tmp[1]
lbPFS[k] = tmp[3]
ubPFS[k] = tmp[4]
}
res.surv = data.frame(sub, num_high, num_low, coxph.pvalOS, hrOS, lbOS, ubOS, coxph.pvalPFS, hrPFS, lbPFS, ubPFS)
res.surv2 = res.surv
###
res.surv = rbind(res.surv1, res.surv2)
##
coxph.pvalPFS = res.surv$coxph.pvalPFS
lbPFS = res.surv$lbPFS
hrPFS = res.surv$hrPFS
ubPFS = res.surv$ubPFS
sub = res.surv$sub
num_high = res.surv$num_high
num_low = res.surv$num_low
coxph.pvalPFS.ann = NULL
for (i in 1:length(coxph.pvalPFS))
{
tmp = coxph.pvalPFS[i]
coxph.pvalPFS.ann[i] = as.character(cut(tmp,
breaks = c(-Inf, 0.001, 0.1, Inf),
labels = c(paste("p =", formatC(tmp, 1, format = "e")),
paste("p =", formatC(signif(tmp, 2), 2, flag = "#")),
"p > 0.1")))
}
sub2 = c("EGFR", "STK11", "KEAP1", "KRAS", "TP53", " Impactful 1", " Impactful 2", " DNE,LOF", " notDNE,LOF", " notDNE,notLOF",
" Missense", " Nonsense", " Splice Site", " Unique Splice Site", " Frameshift")
sub2 = paste0(sub2, " (", num_high, ")")
res1 = data.frame(lbPFS, hrPFS, ubPFS)
res2 = data.frame(sub2, coxph.pvalPFS.ann)
##
res1 = res1[-1,]
res2 = res2[-1,]
## set range hr >18 -> 18
res1.mod = res1
se = which(res1.mod[,"ubPFS"] > 18)
res1.mod[se,"ubPFS"] = 18
##############
pdf(paste0(figures.dir, "Ext_Fig8B_right.pdf"), useDingbats = F, 9, 9)
forestplot(res2, as.matrix(res1.mod), graph.pos = 2, new_page = F,
xlab = "Hazard ratio",
xticks = c(0, 1, 5, 10, 15, 18),
graphwidth = unit(90, "mm"),
vertices = TRUE,
col = fpColors(lines = "black"),
boxsize = 0.3,
txt_gp = fpTxtGp(xlab = gpar(cex = 2.25),
ticks = gpar(cex = 1.75),
title = gpar(cex = 2.5, fontface = "bold"),
cex = 1.5),
align = c("l", "l"),
zero = 1)
dev.off()