.Rhistory

cbttc.clin[pts %in% cbttc.clin$Patient.ID,c("Patient.ID", "OS_MONTHS", "OS_STATUS")]
cbttc.clin[which(pts %in% cbttc.clin$Patient.ID),c("Patient.ID", "OS_MONTHS", "OS_STATUS")]
pts
table(droplevels(na.surv)$Kids_First_Participant_ID)
pts <- na.surv$Kids_First_Participant_ID
pts
na.surv$Kids_First_Participant_ID<- droplevels(na.surv$Kids_First_Participant_ID)
pts <- na.surv$Kids_First_Participant_ID
pts
cbttc.clin[which(pts %in% cbttc.clin$Patient.ID),c("Patient.ID", "OS_MONTHS", "OS_STATUS")]
pts
cbttc.clin[na.surv$Kids_First_Participant_ID %in% cbttc.clin$Patient.ID,c("Patient.ID", "OS_MONTHS", "OS_STATUS")]
na.surv$Kids_First_Participant_ID
cbttc.clin[which(na.surv$Kids_First_Participant_ID %in% cbttc.clin$Patient.ID),c("Patient.ID", "OS_MONTHS", "OS_STATUS")]
na.surv$Kids_First_Participant_ID
na.surv
View(cbttc.clin)
load("~/Documents/GitHub/pptc-pdx-RNA-Seq-clustering/PPTC_GSEA_Input.rda", verbose = T)
View(pptc_gsea)
View(clinical_data)
filt.ann <- read.delim("~/Box Sync/D3B-share/projects/pbta-fusions/fusion_filt_pipeline_v16_all_PBTA/results/Filtered_Annotated_Fusions.txt",
sep  = "\t", header = T, check.names = F)
View(filt.ann)
driver <- read.delim("~/Box Sync/D3B-share/projects/pbta-fusions/fusion_filt_pipeline_v16_all_PBTA/results/PutativeDriver_Annotated_Fusions.txt",
sep  = "\t", header = T, check.names = F)
View(driver)
melted_cormat <- melt(cormat)
melted_cormat
# working directory (created with git clone)
mainDir <- "~/Documents/GitHub/pptc-pdx-classifier-analysis/"
dataDir <- paste0(mainDir,"data/")
# set path to your git cloned repo
script.folder <- paste0(mainDir, "R/")
##create directory for output files
ifelse(!dir.exists(file.path(paste0(mainDir, "results/"))), dir.create(file.path(paste0(mainDir, "results/"))),
"Directory exists!")
results.folder <- paste0(mainDir, "results/")
##set wd
setwd(mainDir)
####Install package Dependencies if needed
source(paste0(script.folder, "install.packages.R"))
####Source plot theme script
source(paste0(script.folder, "theme.R"))
####Load histology color function
source(paste0(script.folder, "mutation-color-function.R"))
###ROC CURVE
roc_file <- file.path(dataDir, "full_roc_threshold_results.tsv")
roc_df <- readr::read_tsv(roc_file,
col_types = readr::cols(
.default = readr::col_double(),
gene = readr::col_character(),
shuffled = readr::col_character()
)
)
curve_colors <- c(
"#313695",
"#35978f",
"#c51b7d",
"#313695",
"#35978f",
"#c51b7d"
)
curve_labels <- c(
"TP53 False" = "TP53 (AUROC = 0.89)",
"Ras False" = "Ras (AUROC = 0.55)",
"NF1 False" = "NF1 (AUROC = 0.77)",
"TP53 True" = "TP53 Shuffled (AUROC = 0.46)",
"Ras True" = "Ras Shuffled (AUROC = 0.55)",
"NF1 True" = "NF1 Shuffled (AUROC = 0.43)"
)
roc_df$model_groups <- paste(roc_df$gene, roc_df$shuffled)
roc_df$model_groups <- factor(roc_df$model_groups, levels = names(curve_labels))
ggplot(roc_df,
aes(x = fpr,
y = tpr,
color = model_groups)) +
coord_fixed() +
geom_step(aes(linetype = shuffled), size = 0.7) +
geom_segment(aes(x = 0,
y = 0,
xend = 1,
yend = 1),
linetype = "dashed",
color = "black") +
scale_y_continuous(labels = scales::percent) +
scale_x_continuous(labels = scales::percent) +
scale_color_manual(name = "Classifier",
values = curve_colors,
labels = curve_labels) +
scale_linetype_manual(name = "Data",
values = c("solid",
"dashed"),
labels = c("True" = "Shuffled",
"False" = "Real")) +
xlab("False Positive Rate") +
ylab("True Positive Rate") +
theme_Publication()
##read in variants
variants <- read.delim(paste0(dataDir, "ras-tp53-nf1-alterations.txt"), sep = "\t", header = T, as.is = T)
##subset tp53 variants
p53.var <- subset(variants, classifier == "TP53")
##load and merge with clinical file
clin = read.delim(paste0(dataDir, "pptc-pdx-clinical-web.txt"), sep = "\t", header = T, as.is = T)
##select only those models with RNA-Seq
rna <- subset(clin, RNA.Part.of.PPTC == "yes")
##read in scores
scores <- read.delim(paste0(dataDir, "classifier_scores.tsv"),
sep = "\t", header = T, as.is = T)
colnames(scores)[colnames(scores) == "sample_id"] <- "Model"
scores.clin <- merge(scores, clin[,c("Model", "Histology.Detailed")], all.x = T)
##select only tp53 scores to work with
tp53scores <- scores.clin[,c("Model", "Histology.Detailed", "tp53_score")]
###add tp53 variants to tp53 scores df
var.score <- merge(tp53scores, p53.var[,c("Model", "Hugo_Symbol", "classifier", "Variant_Classification", "cDNA_Change", "Protein_Change", "Fusion")], all.x = T)
##add classification to WT samples
var.score$Variant_Classification <- ifelse(is.na(var.score$Hugo_Symbol), "wild-type", var.score$Variant_Classification)
var.score$Hugo_Symbol <- ifelse(is.na(var.score$Hugo_Symbol), "wild-type", var.score$Hugo_Symbol)
var.score$classifier <- ifelse(is.na(var.score$classifier), "unclassified", var.score$classifier)
###reorder
var.score$Hugo_Symbol <- factor(var.score$Hugo_Symbol, levels = c("wild-type", "TP53", "CDKN2A", "MDM2", "MDM4", "GORAB", "ATM", "ATR", "RB1", "CHEK1", "CHEK2"))
###plot for figure 4D
pdf(paste(results.folder, Sys.Date(), "-tp53-scoresbyvarianttype.pdf", sep = ""), width = 20, height = 3)
ggplot(var.score, aes(x = Variant_Classification, y = tp53_score, color = Variant_Classification, alpha = 0.5)) +
geom_boxplot(outlier.shape = 21, fill = 'white') +
geom_point(pch = 21, position = position_jitterdodge())+
geom_hline(aes(yintercept = 0.5), colour = 'black', linetype = "longdash")+
stat_boxplot(geom ='errorbar', width = 0.5) +
facet_grid(~Hugo_Symbol, scales = "free_x")+
guides(alpha = FALSE, fill = FALSE) +
theme_Publication() + xlab("") + ylab('TP53 Classifier Score') +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank()) +
ylim(c(0,1))+
scale_color_manual(values = colores)
dev.off()
###subset only osteosarcoma scores for supplemental figures
ost <- subset(var.score, Histology.Detailed == "Osteosarcoma")
ost$Protein_Change <- ifelse(ost$Variant_Classification == "Splice_Site", "splice", ost$Protein_Change)
ost$Alteration <- ifelse(is.na(ost$Protein_Change) & is.na(ost$Fusion),
paste(ost$Hugo_Symbol, ost$Variant_Classification, sep = " "),
ifelse(is.na(ost$Protein_Change), paste0(ost$Fusion, " Fusion"),
paste(ost$Hugo_Symbol, ost$Protein_Change, sep = " ")))
table(ost$Alteration)
ost.gather <- ost[,c("Model", "tp53_score", "Alteration")]
united <- ost.gather %>%
group_by(Model, tp53_score) %>%
summarise_all(funs(trimws(paste(., collapse = ', '))))
print(ggplot(ost, aes(x = Variant_Classification, y = tp53_score, color = Variant_Classification, alpha = 0.5)) +
geom_boxplot(outlier.shape = 21, fill = 'white') +
geom_point(pch = 21, position = position_jitterdodge())+
geom_hline(aes(yintercept = 0.5), colour = 'black', linetype = "longdash")+
stat_boxplot(geom ='errorbar', width = 0.5) +
facet_grid(~Hugo_Symbol, scales = "free_x")+
guides(alpha = FALSE, fill = FALSE) +
theme_Publication() + xlab("") + ylab('TP53 Classifier Score') +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank()) +
ylim(c(0,1))+
scale_color_manual(values = colores))
####remove osteosarcoma to replot
no.ost <- subset(var.score, Histology.Detailed != "Osteosarcoma")
print(ggplot(no.ost, aes(x = Variant_Classification, y = tp53_score, color = Variant_Classification, alpha = 0.5)) +
geom_boxplot(outlier.shape = 21, fill = 'white') +
geom_point(pch = 21, position = position_jitterdodge())+
geom_hline(aes(yintercept = 0.5), colour = 'black', linetype = "longdash")+
facet_grid(~Hugo_Symbol, scales = "free_x")+
stat_boxplot(geom ='errorbar', width = 0.5) +
guides(alpha = FALSE, fill = FALSE) +
theme_Publication() + xlab("") + ylab('TP53 Classifier Score') +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank()) +
ylim(c(0,1))+
scale_color_manual(values = colores))
###define CNVs
cnvs <- c("Hemizygous_Deletion", "Amplification", "Homozygous_Deletion")
####add CNV as a variant type for comparison by type
var.score$type <- ifelse(var.score$Variant_Classification %in% cnvs, "CNV",
ifelse (var.score$Variant_Classification == "Fusion", "Fusion",
ifelse (var.score$Variant_Classification == "wild-type", "wild-type", "SNV")))
##relevel
var.score$type <- factor(var.score$type, levels = c("wild-type", "SNV", "CNV", "Fusion"))
print(ggplot(var.score, aes(x = Hugo_Symbol, y = tp53_score, alpha = 0.5)) +
geom_boxplot(outlier.shape = 21, fill = 'white') +
geom_jitter(position=position_jitter(width=.1, height=0), shape = 21) +
stat_boxplot(geom ='errorbar', width = 0.5) +
geom_hline(aes(yintercept = 0.5), colour = 'black', linetype = "longdash")+
guides(alpha = FALSE, fill = FALSE) +
theme_Publication() + xlab("Gene") + ylab('TP53 Classifier Score') +
ylim(c(0,1))+
scale_color_manual(values = colores)+
stat_n_text() +
scale_fill_manual(values = colores))
mycomps <- list(c("wild-type", "SNV"), c("wild-type", "CNV"), c("wild-type", "Fusion"))
ggviolin(var.score, x = "type", y = "tp53_score", fill = "type",alpha = 0.5,
palette =colores, add.params = list(fill = "white"),
shape = 1,
add = "boxplot",
bxp.errorbar = T)+
stat_compare_means(comparisons = mycomps, label = "p.format", method = "wilcox.test", ref.group = "wild-type", step.increase = 0.2) + ##posthoc p-values
stat_compare_means(label.y = 1.7) + #global p value
geom_hline(aes(yintercept = 0.5), colour = 'black', linetype = "longdash")+
theme_Publication() +
xlab("Type of Alteration") +
ylab('TP53 Classifier Score') +
ylim(c(0,1.7))
print(ggplot(var.score, aes(x = type, y = tp53_score, color = type, alpha = 0.5)) +
geom_boxplot(outlier.shape = 21, fill = 'white') +
geom_jitter(position=position_jitter(width=.1, height=0), shape = 21) +
geom_hline(aes(yintercept = 0.5), colour = 'black', linetype = "longdash")+
facet_wrap(~Hugo_Symbol, scales = "free_x", nrow = 2)+
stat_boxplot(geom ='errorbar', width = 0.5) +
guides(alpha = FALSE) +
theme_Publication() + xlab("Type of Alteration") + ylab('TP53 Classifier Score') +
ylim(c(0,1))+
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())+ #+
scale_color_manual(values = colores)+
scale_fill_manual(values = colores))
###overall classifier status
var.score$status <- ifelse(var.score$Variant_Classification == "wild-type", "wild-type", "altered")
var.score$status <- factor(var.score$status, levels = c("wild-type", "altered"))
####add only tp53 alterations
var.score$p53status <- ifelse(var.score$Hugo_Symbol == "TP53" & var.score$status != "wild-type", "altered", "wild-type")
var.score$p53status <- factor(var.score$p53status, levels = c("wild-type", "altered"))
var.score.noost <- subset(var.score, Histology.Detailed != "Osteosarcoma")
####figure 4A
table(var.score$status) ###duplicate samples included due to multiple TP53 alterations, so remove those
tp53 <- subset(var.score, Hugo_Symbol == "TP53")
#unique altered models
un.tp53 <- unique(tp53[,c("Model", "status")])
###merge with all scores
un.tp53<- merge(un.tp53, rna[,c("Model", "Histology.Detailed")], all.y = T)
un.tp53$status <-ifelse(is.na(un.tp53$status), paste0("wild-type"), paste0(un.tp53$status))
un.tp53.scores <- merge(un.tp53, tp53scores[,c("Model", "tp53_score")])
table(un.tp53.scores$status) ###new Ns for figure 4A
print(ggviolin(un.tp53.scores, x = "status", y = "tp53_score", fill = "status",
palette = colores,alpha = 0.5,
add = "boxplot",
add.params = list(fill = "white"))+
stat_compare_means(comparisons = list(c("wild-type", "altered")), label = "p.signif")+ # Add significance levels
stat_compare_means(label.y = 1.5) + ##global p
theme_Publication() +
xlab("Mutation Status") + ylab('TP53 Classifier Score'))
wilcox.test(tp53_score ~ status, un.tp53.scores)
median(as.numeric(unlist(subset(un.tp53.scores, status=="wild-type", select=tp53_score))))
median(as.numeric(unlist(subset(un.tp53.scores, status=="altered", select=tp53_score))))
###subset for models not osteosarcoma
new.noost <- subset(un.tp53.scores, Histology.Detailed != "Osteosarcoma")
table(new.noost$status) ##Ns
length(subset(clin, Histology.Detailed == "Osteosarcoma" & RNA.Part.of.PPTC == "yes"))
print(ggviolin(new.noost, x = "status", y = "tp53_score", fill = "status",
palette = colores,alpha = 0.5,
add = "boxplot",
add.params = list(fill = "white"))+
stat_compare_means(comparisons = list(c("wild-type", "altered")), label = "p.signif")+ # Add significance levels
stat_compare_means(label.y = 1.5) + ##global p
theme_Publication() +
xlab("Mutation Status") + ylab('TP53 Classifier Score'))
wilcox.test(tp53_score ~ status, un.tp53.scores)
scores
# working directory (created with git clone)
mainDir <- "~/Documents/GitHub/pptc-pdx-classifier-analysis/"
#read in tmb data
tmb <- read.delim(paste0(dataDir, "mutations-per-model.txt"), sep = "\t",
header = T, as.is = T)
#read in signatures
sigs <- read.delim(paste0(dataDir, "allpdx-signature-weights.txt"),
sep = "\t", header = T, as.is = T)
unique(sigs$Model)
sigs$Model <- NULL
##remove sigs with no values and unknown sig
Summat <-  colSums(sigs,na.rm = TRUE)
nonzero.sigs <- sigs[, Summat!=0]
nonzero.sigs$Unknown <- NULL
nonzero.sigs$Model <- rownames(nonzero.sigs)
##merge breakpoints and scores
bp.scores <- merge(bkpt[,c("Model","n.breakpoints")], scores, all = T)
#read in tmb data
bp <- read.delim(paste0(dataDir, "Breakpoints_rawdata.txt"), sep = "\t",
header = T, as.is = T)
tmb <- read.delim(paste0(dataDir, "mutations-per-model.txt"), sep = "\t",
header = T, as.is = T)
#read in signatures
sigs <- read.delim(paste0(dataDir, "allpdx-signature-weights.txt"),
sep = "\t", header = T, as.is = T)
sigs$Model <- NULL
##remove sigs with no values and unknown sig
Summat <-  colSums(sigs,na.rm = TRUE)
nonzero.sigs <- sigs[, Summat!=0]
nonzero.sigs$Unknown <- NULL
nonzero.sigs$Model <- rownames(nonzero.sigs)
##merge breakpoints and scores
bp.scores <- merge(bp[,c("Model","n.breakpoints")], scores, all = T)
summary(lm(bp.scores$n.breakpoints ~ bp.scores$tp53_score))
View(bp.scores)
##add tmb to the merge
bp.scor.tmb <- merge(bp.scores, tmb[,c("Model", "MutperMB")], all = T)
View(bp.scor.tmb)
bp.scor.tmb$tmb <- log(bp.scor.tmb$MutperMB, 2)
View(bp.scor.tmb)
###add signatures to merge
all.mer <- merge(bp.scor.tmb, nonzero.sigs, all = T)
View(all.mer)
nonzero.sigs
#read in signatures
sigs <- read.delim(paste0(dataDir, "allpdx-signature-weights.txt"),
sep = "\t", header = T, as.is = T)
rownames(sigs) <- sigs$Model
sigs$Model <- NULL
##remove sigs with no values and unknown sig
Summat <-  colSums(sigs,na.rm = TRUE)
nonzero.sigs <- sigs[, Summat!=0]
nonzero.sigs$Unknown <- NULL
nonzero.sigs$Model <- rownames(nonzero.sigs)
##merge breakpoints and scores
bp.scores <- merge(bp[,c("Model","n.breakpoints")], scores, all = T)
summary(lm(bp.scores$n.breakpoints ~ bp.scores$tp53_score))
##add tmb to the merge
bp.scor.tmb <- merge(bp.scores, tmb[,c("Model", "MutperMB")], all = T)
##log tmb
bp.scor.tmb$tmb <- log(bp.scor.tmb$MutperMB, 2)
###add signatures to merge
all.mer <- merge(bp.scor.tmb, nonzero.sigs, all = T)
View(all.mer)
View(all.mer)
##move model to rownames for correlation plots
rownames(bp.scor.tmb) <- bp.scor.tmb$Model
###create first cor matrix
first <- bp.scor.tmb[,c("tp53_score", "n.breakpoints", "tp53_shuffle", "tmb")]
###remove ost models to see if correlation still holds
no.ost.models <- var.score.noost$Model
first.no.ost <- first[no.ost.models,]
cor.test(first$tp53_score, first$n.breakpoints, method = "pearson")
cor.test(first.no.ost$tp53_score, first.no.ost$n.breakpoints, method = "pearson") ##yes
cor.test(first$tmb, first$n.breakpoints, method = "pearson")
#remove model column
cor.df <- all.mer[,2:ncol(all.mer)]
####create correlation matrix
cormat1 <- round(cor(first, use = "pairwise", method = "pearson"),2)
# Get lower triangle of the correlation matrix
get_lower_tri<-function(cormat1){
cormat[upper.tri(cormat1)] <- NA
return(cormat1)
}
# Get upper triangle of the correlation matrix
get_upper_tri <- function(cormat1){
cormat[lower.tri(cormat1)]<- NA
return(cormat1)
}
source(paste0(script.folder, "cormat.R"))
####create correlation matrix
cormat <- round(cor(first, use = "pairwise", method = "pearson"),2)
source(paste0(script.folder, "cormat.R"))
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "#005789", high = "#D55E00", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_Publication()+ # minimal theme
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
coord_fixed()
# Print the heatmap with labels
print(ggheatmap +
geom_text(aes(Var2, Var1, label = value), color = "black", size = 6) +
theme_Publication()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
labs(y ="", x = ""))
#######second cor matrix:
second <- cor.df[,c("tp53_score", "n.breakpoints", "tp53_shuffle", "nf1_score", "nf1_shuffle",
"Signature.1", "Signature.13", "Signature.2", "Signature.15", "tmb",
"Signature.3", "Signature.6", "Signature.20", "Signature.26", "Signature.9",
"Signature.10", "Signature.21")]
View(cor.df)
names(cor.df)
#remove model column
cor.df <- all.mer[,2:ncol(all.mer)]
names(cor.df)
#read in signatures
sigs <- read.delim(paste0(dataDir, "allpdx-signature-weights.txt"),
sep = "\t", header = T, as.is = T)
View(sigs)
#######second cor matrix: - sig 9 not in cor.df anymore
second <- cor.df[,c("tp53_score", "n.breakpoints", "tp53_shuffle", "nf1_score", "nf1_shuffle",
"Signature.1", "Signature.13", "Signature.2", "Signature.15", "tmb",
"Signature.3", "Signature.6", "Signature.20", "Signature.26",
"Signature.10", "Signature.21")]
names(cor.df)
cor.test(second$Signature.1, second$Signature.6, method = "pearson")
cor.test(second$Signature.1, second$Signature.3, method = "pearson")
cor.test(second$Signature.2, second$Signature.13, method = "pearson")
cormat <- round(cor(second, use = "pairwise", method = "pearson"),2)
source(paste0(script.folder, "cormat.R"))
# Create a ggheatmap
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "#005789", high = "#D55E00", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_Publication()+ # minimal theme
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
coord_fixed()
# Print the heatmap
print(ggheatmap +
geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
theme_Publication()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
labs(y ="", x = ""))
cor.test(second$Signature.1, second$Signature.6, method = "pearson")
#read in signatures
sigs <- read.delim(paste0(dataDir, "allpdx-signature-weights.txt"),
sep = "\t", header = T, as.is = T)
rownames(sigs) <- sigs$Model
sigs$Model <- NULL
##remove sigs with no values and unknown sig
Summat <-  colSums(sigs,na.rm = TRUE)
nonzero.sigs <- sigs[, Summat!=0]
nonzero.sigs$Unknown <- NULL
nonzero.sigs$Model <- rownames(nonzero.sigs)
##merge breakpoints and scores
bp.scores <- merge(bp[,c("Model","n.breakpoints")], scores, all = T)
summary(lm(bp.scores$n.breakpoints ~ bp.scores$tp53_score))
##add tmb to the merge
bp.scor.tmb <- merge(bp.scores, tmb[,c("Model", "MutperMB")], all = T)
##log tmb
bp.scor.tmb$tmb <- log(bp.scor.tmb$MutperMB, 2)
###add signatures to merge
all.mer <- merge(bp.scor.tmb, nonzero.sigs, all = T)
##move model to rownames for correlation plots
rownames(bp.scor.tmb) <- bp.scor.tmb$Model
###create first cor matrix
first <- bp.scor.tmb[,c("tp53_score", "n.breakpoints", "tp53_shuffle", "tmb")]
###remove ost models to see if correlation still holds
no.ost.models <- var.score.noost$Model
first.no.ost <- first[no.ost.models,]
cor.test(first$tp53_score, first$n.breakpoints, method = "pearson")
cor.test(first.no.ost$tp53_score, first.no.ost$n.breakpoints, method = "pearson") ##yes
cor.test(first$tmb, first$n.breakpoints, method = "pearson")
#remove model column
cor.df <- all.mer[,2:ncol(all.mer)]
####create correlation matrix
cormat <- round(cor(first, use = "pairwise", method = "pearson"),2)
source(paste0(script.folder, "cormat.R"))
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "#005789", high = "#D55E00", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_Publication()+ # minimal theme
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
coord_fixed()
# Print the heatmap with labels
print(ggheatmap +
geom_text(aes(Var2, Var1, label = value), color = "black", size = 6) +
theme_Publication()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
labs(y ="", x = ""))
#######second cor matrix: - sig 9 not in cor.df anymore
second <- cor.df[,c("tp53_score", "n.breakpoints", "tp53_shuffle", "nf1_score", "nf1_shuffle",
"Signature.1", "Signature.13", "Signature.2", "Signature.15", "tmb",
"Signature.3", "Signature.6", "Signature.20", "Signature.26", "Signature.9",
"Signature.10", "Signature.21")]
names(cor.df)
cor.test(second$Signature.1, second$Signature.6, method = "pearson")
cor.test(second$Signature.1, second$Signature.3, method = "pearson")
cor.test(second$Signature.2, second$Signature.13, method = "pearson")
cormat <- round(cor(second, use = "pairwise", method = "pearson"),2)
source(paste0(script.folder, "cormat.R"))
# Create a ggheatmap
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "#005789", high = "#D55E00", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_Publication()+ # minimal theme
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
coord_fixed()
# Print the heatmap
print(ggheatmap +
geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
theme_Publication()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
labs(y ="", x = ""))
###third corplot
third <- cor.df[,c("tp53_score", "n.breakpoints", "tp53_shuffle", "nf1_score", "nf1_shuffle",
"Signature.1", "Signature.13", "Signature.2", "Signature.15", "tmb",
"Signature.3", "Signature.6")]
cormat <- round(cor(third, use = "pairwise", method = "pearson"),2)
source(paste0(script.folder, "cormat.R"))
# Create a ggheatmap
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "#005789", high = "#D55E00", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_Publication()+ # minimal theme
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
coord_fixed()
# Print the heatmap
print(ggheatmap +
geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
theme_Publication()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
labs(y ="", x = ""))