Chapter 3 Target Evaluation
3.1 Load packages
library(openxlsx)
library(dplyr)
library(DT)
library(nlme)
library(DESeq2)
library(SummarizedExperiment)
library(openxlsx)
library(dendextend)
library(matrixStats)
library(data.table)
library(ggpubr)
library(limma)
library(DESeq2)
library(pheatmap)
library(gplots)
library(survminer)
library(survival)
library(enrichR)
library(stringr)
library(ggthemes)
library(fmsb)
# color theme
tp53.mutant <- c("WT","L130V","L130F","L130R","P142L","V143M","A161T","M237V","G245D","F270I","F270V","V272A","V272M","R248Q","R273H","null")
tp53.mutant.color <- c("red",rep("#494949","14"),"#E3E3E3")
names(tp53.mutant.color) <- tp53.mutant
anno.color <- list(
Mutant = tp53.mutant.color)
3.2 Merge profiling
We provided the merged profiling, which could be downloaded via https://github.com/NRCTM-bioinfo/p53_compounds.
tp53.targets <- read.xlsx("01.Profiling/p53_targets.xlsx")
rownames(tp53.targets) <- tp53.targets$Gene
tp53.target.g32 <- tp53.targets$Gene[which(tp53.targets$LiteratureNum >= 10)]
tp53.target.g116 <- tp53.targets$Gene
merge.meta.data <- read.xlsx("01.Profiling/p53_meta.xlsx", sheet = 1)
merge.compare.data <- read.xlsx("01.Profiling/Table1.xlsx", sheet = 1, startRow = 2)
merge.compare.data$Sample <- merge.compare.data$`Sample.name.(after.treatment)`
merge.compare.data$Control <- merge.compare.data$`Sample.name.(untreated.samples)`
merge.compare.data$Group <- str_replace(merge.compare.data$Sample, regex("_rep.$"), "")
rownames(merge.compare.data) <- merge.compare.data$`Sample.name.(after.treatment)`
merge.mat <- read.xlsx("01.Profiling/Table2.xlsx", sheet = 1, startRow = 3, rowNames = T)
3.3 Normalization
# Median normalization
merge.nor.mat <- merge.mat
sub.1 <- aggregate(colMedians(as.matrix(merge.nor.mat)),
list(Cell = merge.meta.data$Dataset), min)
for (i in 1:ncol(merge.nor.mat)) {
merge.nor.mat[, i] <- merge.nor.mat[, i] - sub.1$x[which(sub.1$Cell == merge.meta.data$Dataset[i])]
}
# Generation of fold change (Log2-transformed for FPKM) of cell line samples after treatment compared to before treatment
compare.mat <- NULL
for (i in 1:nrow(merge.compare.data)) {
if (length(grep("; ", merge.compare.data$Control[i])) > 0) {
sub.1 <- unlist(str_split(merge.compare.data$Control[i], "; "))
sub <- merge.nor.mat[, merge.compare.data$Sample[i]] - rowMeans(merge.nor.mat[, sub.1])
} else {
sub <- merge.nor.mat[, merge.compare.data$Sample[i]] - merge.nor.mat[, merge.compare.data$Control[i]]
}
compare.mat <- cbind(compare.mat, sub)
}
colnames(compare.mat) <- merge.compare.data$Sample
head(compare.mat, 3)
## U937_WT U937_L130V_ATO U937_L130F_ATO U937_L130R_ATO U937_P142L_ATO U937_V143M_ATO U937_M237V_ATO
## CDKN1A 3.1916457 3.080478 1.957199 2.0917992 0.7774611 2.1667723 1.9077609
## RRM2B 0.8391627 1.001993 1.187001 0.6459767 0.4785670 0.8683713 0.5825729
## MDM2 2.7245150 2.028156 1.302068 1.1567997 0.8098108 1.5695778 1.3455218
## U937_F270I_ATO U937_F270V_ATO U937_V272A_ATO U937_V272M_ATO U937_V272M_PAT
## CDKN1A 3.947027 2.909291 2.035912 3.675682 3.269377
## RRM2B 1.860397 1.479516 0.708025 1.567485 1.557309
## MDM2 2.642251 2.444877 1.671464 3.256886 2.245890
## Kasumi-1_R248Q_Kevetrin_6h_rep1 Kasumi-1_R248Q_Kevetrin_6h_rep2 Kasumi-1_R248Q_Kevetrin_6h_rep3
## CDKN1A -0.00507247 -0.03563892 0.17165276
## RRM2B 0.66265896 0.57210179 0.51180155
## MDM2 0.04900353 -0.04690794 -0.01286232
## SKM1_R248Q_APR-246_rep1 SKM1_R248Q_APR-246_rep2 SKM1_R248Q_APR-246_rep3 staET7.1_R273H_APR-246_rep1
## CDKN1A -0.07281667 0.02994333 -0.26504667 0.4066695
## RRM2B 0.13078333 0.05577333 -0.27296667 0.6471832
## MDM2 0.18355667 -0.07633333 -0.09720333 0.2011240
## staET7.1_R273H_APR-246_rep2 staET7.2_R273H_APR-246_rep1 staET7.3_R273H_APR-246_rep1
## CDKN1A -0.31610198 -0.1298085 -0.41130580
## RRM2B 0.12121419 0.3475557 0.07482002
## MDM2 -0.04862277 0.1818754 0.07292721
## staET7.3_R273H_APR-246_rep2 staET7.2_R273H_APR-246_rep2 Saos-2_R273H_APR-246_6h_rep1
## CDKN1A -0.091652792 1.1191052 -0.02540054
## RRM2B 0.356190359 0.5005976 1.05526980
## MDM2 0.001335775 0.2663554 -0.08762767
## Saos-2_R273H_APR-246_6h_rep2 Saos-2_R273H_APR-246_6h_rep3 Saos-2_R273H_APR-246_12h_rep1
## CDKN1A -0.1177016 -0.534773959 -0.62120069
## RRM2B 1.1814193 0.431825568 1.03123586
## MDM2 0.3947754 -0.008325736 0.06401986
## Saos-2_R273H_APR-246_12h_rep2 Saos-2_R273H_APR-246_12h_rep3 Kasumi-1_R248Q_Kevetrin_48h_rep1
## CDKN1A -0.10077979 -0.3485730 0.3656561
## RRM2B 0.75826079 0.6167059 1.0807021
## MDM2 0.02695217 0.0326328 -0.4937894
## Kasumi-1_R248Q_Kevetrin_48h_rep2 Kasumi-1_R248Q_Kevetrin_48h_rep3 PCI13_G245D_COTI-2_rep1
## CDKN1A 0.6150592 0.5977127 0.7164330
## RRM2B 1.0029696 1.2170216 0.1130768
## MDM2 -0.5052966 -0.2238308 0.2377859
## PCI13_G245D_COTI-2_rep2 PCI13_G245D_COTI-2_rep3 U266_A161T_PRIMA-1_rep1 U266_A161T_PRIMA-1_rep2
## CDKN1A 0.6274429 0.647518846 0.1974164 0.4389011
## RRM2B -0.0679887 -0.006853216 -0.2733003 -0.3027616
## MDM2 0.2552303 0.341083919 -0.2199022 -0.2003108
## U937_R273H_ATO U937_R273H_PAT PCI13_G245D_CDDP_rep1 PCI13_G245D_CDDP_rep2 PCI13_G245D_CDDP_rep3
## CDKN1A 0.3927590 0.1909721 -0.21575276 -0.29477320 -0.3104395
## RRM2B 0.2643679 -0.1862437 -0.28507761 -0.56260138 -0.4041368
## MDM2 -0.5397466 -0.4929459 0.06734607 0.07106408 0.1046652
## 8266R5_null_PRIMA-1_rep1 8266R5_null_PRIMA-1_rep2 Saos-2_null_APR-246_6h_rep1 Saos-2_null_APR-246_6h_rep2
## CDKN1A -0.1462286 -0.1858493 1.0204064 -0.2922208
## RRM2B -0.2320481 0.2372455 0.2411493 0.4065886
## MDM2 0.2465478 0.3169312 -0.6160103 -0.6468282
## Saos-2_null_APR-246_6h_rep3 Saos-2_null_APR-246_12h_rep1 Saos-2_null_APR-246_12h_rep2
## CDKN1A 0.16556354 -0.252144504 -0.2329088
## RRM2B -0.31879809 0.490890719 -0.4246623
## MDM2 -0.05262877 0.001146008 -0.1160399
3.4 32 confident p53 targets
The heatmap plot was used for visualization of the log2-transformed fold change of 32 confident p53 targets.
my.breaks <- c(seq(-1.5, -0.01, by = 0.001), seq(0.01, 1.5, by = 0.001) )
my.colors <- c(colorRampPalette(colors = c("#00599F","#00599F","#287dba","#62aee5","#a8dcff","white"))(length(my.breaks)/2),
colorRampPalette(colors = c("white","#ffca79","#ef6e51","#db1b18","#b50600","#b50600"))(length(my.breaks)/2))
plot.mat <- compare.mat[tp53.target.g32, ]
anno.data <- data.frame(row.names = merge.compare.data$Sample,
Mutant = merge.compare.data$p53.status)
# Visualization of the 32 confident p53 targets
p <- pheatmap(plot.mat, scale = "none",
color = my.colors, breaks = my.breaks,
cluster_row = F, cluster_col = F, border_color = NA,
annotation_col = anno.data,
annotation_colors = anno.color,
clustering_method = "ward.D2",
clustering_distance_rows = "manhattan",
clustering_distance_cols = "manhattan",
fontsize_col = 6,
fontsize_row = 6)
p
# Visualization of the other 18,264 genes (exclude the 32 confident p53 targets)
plot.mat <- compare.mat[!rownames(compare.mat) %in% tp53.target.g32, ]
p <- pheatmap(plot.mat, scale = "none",
color = my.colors, breaks = my.breaks,
cluster_row = F, cluster_col = F, border_color = NA,
annotation_col = anno.data,
annotation_colors = anno.color,
clustering_method = "ward.D2",
clustering_distance_rows = "manhattan",
clustering_distance_cols = "manhattan",
fontsize_col = 6,
fontsize_row = 0.1)
p
3.5 116 p53 targets
The heatmap plot was used for visualization of the log2-transformed fold change of the 116 p53 targets.
my.breaks <- c(seq(-1.5, -0.01, by = 0.001), seq(0.01, 1.5, by = 0.001) )
my.colors <- c(colorRampPalette(colors = c("#00599F","#00599F","#287dba","#62aee5","#a8dcff","white"))(length(my.breaks)/2),
colorRampPalette(colors = c("white","#ffca79","#ef6e51","#db1b18","#b50600","#b50600"))(length(my.breaks)/2))
plot.mat <- compare.mat[tp53.target.g116, ]
anno.data <- data.frame(row.names = merge.compare.data$Sample, Mutant = merge.compare.data$p53.status)
# Visualization of the 116 p53 targets
p <- pheatmap(plot.mat, scale = "none",
color = my.colors, breaks = my.breaks,
cluster_row = F, cluster_col = F, border_color = NA,
annotation_col = anno.data,
annotation_colors = anno.color,
clustering_method = "ward.D2",
clustering_distance_rows = "manhattan",
clustering_distance_cols = "manhattan",
fontsize_col = 6,
fontsize_row = 2)
p
# Visualization of the other 18,180 genes (exclude the 116 p53 targets)
plot.mat <- compare.mat[!rownames(compare.mat) %in% tp53.target.g116, ]
p <- pheatmap(plot.mat, scale = "none",
color = my.colors, breaks = my.breaks,
cluster_row = F, cluster_col = F, border_color = NA,
annotation_col = anno.data,
annotation_colors = anno.color,
clustering_method = "ward.D2",
clustering_distance_rows = "manhattan",
clustering_distance_cols = "manhattan",
fontsize_col = 8,
fontsize_row = 0.1)
p
3.6 p53 transactivation
plot.order <- c("U937_WT","U937_L130V_ATO","U937_L130F_ATO","U937_L130R_ATO","U937_P142L_ATO","U937_V143M_ATO","U937_M237V_ATO","U937_F270I_ATO","U937_F270V_ATO","U937_V272A_ATO","U937_V272M_ATO","U937_V272M_PAT","Kasumi-1_R248Q_Kevetrin_6h","SKM1_R248Q_APR-246","staET7.1_R273H_APR-246","staET7.2_R273H_APR-246","staET7.3_R273H_APR-246","Saos-2_R273H_APR-246_6h","Saos-2_R273H_APR-246_12h","Kasumi-1_R248Q_Kevetrin_48h","PCI13_G245D_COTI-2","U266_A161T_PRIMA-1","U937_R273H_ATO","U937_R273H_PAT","PCI13_G245D_CDDP","8266R5_null_PRIMA-1","Saos-2_null_APR-246_6h","Saos-2_null_APR-246_12h")
total.act <- NULL
total.test.cutoff <- c(1.2, 1.5, 2.0)
for (i in 1:length(total.test.cutoff)) {
cutoff_tmp <- round(log2(total.test.cutoff[i]), 2)
plot.data <- data.frame(Sample = colnames(compare.mat),
UP_116 = colSums(compare.mat[tp53.target.g116, ] > cutoff_tmp) ,
UP_32 = colSums(compare.mat[tp53.target.g32, ] > cutoff_tmp),
UP_other_116 = colSums(compare.mat[!rownames(compare.mat) %in% tp53.target.g116, ] > cutoff_tmp) ,
UP_other_32 = colSums(compare.mat[!rownames(compare.mat) %in% tp53.target.g32, ] > cutoff_tmp) )
plot.data$Mutant <- merge.compare.data$p53.status[match(plot.data$Sample, merge.compare.data$Sample)]
plot.data$MutantType <- "MUT"
plot.data$MutantType[which(plot.data$Mutant == "WT")] <- "WT"
plot.data$Act_32 <- (plot.data$UP_32) / (plot.data$UP_other_32 + plot.data$UP_32) * 100
plot.data$Act_116 <- (plot.data$UP_116) / (plot.data$UP_other_116 + plot.data$UP_116) * 100
plot.data$Group <- merge.compare.data$Group[match(plot.data$Sample, merge.compare.data$Sample)]
plot.data$Act_32[which(is.na(plot.data$Act_32))] = 1
plot.data$Act_116[which(is.na(plot.data$Act_116))] = 1
plot.data.group <- aggregate(plot.data[, c("Act_32","Act_116")], list(Sample = plot.data$Group), mean)
plot.data.group <- plot.data.group[match(plot.order, plot.data.group$Sample), ]
plot.data.group$Mutant <- merge.compare.data$p53.status[match(plot.data.group$Sample, merge.compare.data$Group)]
plot.data.group$MutantType <- "MUT"
plot.data.group$MutantType[which(plot.data.group$Mutant == "null")] <- "null_other"
plot.data.group$MutantType[which(plot.data.group$Sample %in% c("U937_R273H_ATO","U937_R273H_PAT","PCI13_G245D_CDDP"))] <- "null_other"
plot.data.group$MutantType[which(plot.data.group$Mutant == "WT")] <- "WT"
if (cutoff_tmp != 1) {
group.tag <- paste0("FC_", total.test.cutoff[i] )
} else {
group.tag <- "FC_2.0"
}
plot.data.group$Group <- group.tag
total.act <- rbind(total.act, plot.data.group)
}
total.act$Group <- factor(as.character(total.act$Group), levels = c("FC_1.5","FC_1.2","FC_2.0"))
p <- ggbarplot(total.act, x = "Sample", y = "Act_32",
color = "MutantType", fill = "MutantType",
sorting = "none", size = 0.5,
order = plot.order,
main = paste0("32 confident p53 targets"), xlab = "", ylab = "",
palette = c("black","#A3A3A3","red"), width = 0.6 ) + theme_few()
p <- p + geom_hline(yintercept = 0, linetype = 2, color = "black")
p <- p + theme_base() + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 4))
p <- p + theme(plot.background = element_blank())
p <- p + scale_y_continuous(limits = c(0,11) )
p <- p + facet_wrap( ~ Group, ncol = 3)
p
p <- ggbarplot(total.act, x = "Sample", y = "Act_116",
color = "MutantType", fill = "MutantType",
sorting = "none", size = 0.5,
order = plot.order,
main = paste0("116 p53 targets"), xlab = "", ylab = "",
palette = c("black","#A3A3A3","red"), width = 0.6 ) + theme_few()
p <- p + geom_hline(yintercept = 0, linetype = 2, color = "black")
p <- p + theme_base() + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 4))
p <- p + theme(plot.background = element_blank())
p <- p + scale_y_continuous(limits = c(0,15) )
p <- p + facet_wrap( ~ Group, ncol = 3)
p