-
Notifications
You must be signed in to change notification settings - Fork 0
/
FigS5B.R
121 lines (95 loc) · 5.09 KB
/
FigS5B.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
################################################################################
#################### Setup #####################################################
# loc
setwd("~/git/CAIRO5_Public")
# libraries
library(readxl)
library(data.table)
library(ggplot2)
library(ggallin)
library(scales)
# Data
BloodSampleCharacteristics <- as.data.table(read_xlsx("Data/CAIRO5_Supplemental.xlsx", sheet=2, skip=1))
SummaryCopyNumber <- as.data.table(read_xlsx("Data/CAIRO5_Supplemental.xlsx", sheet=5, skip=1))
################################################################################
################# pair down data to required cols ##############################
# Blood Samples
keep <- c('Subject ID', 'Sample Time Point', 'DELFI-TF', 'MAF', 'ichorCNA')
BloodSampleCharacteristics <- BloodSampleCharacteristics[,..keep]
setkeyv(BloodSampleCharacteristics, c('Subject ID', 'Sample Time Point'))
# CNA data
setkeyv(SummaryCopyNumber, c('Subject ID', 'Sample Time Point'))
# merge data together
finalData <- merge(SummaryCopyNumber, BloodSampleCharacteristics)
################################################################################
################ format data ###################################################
# melt the data first by technology and then by CNA feature
finalData <- melt(finalData, measure.vars=c('DELFI-TF', 'MAF'), variable.name="Technology", value.name="Value")
finalData <- melt(finalData, measure.vars=c('PLCG1 Copy Number', 'MBD1 Copy Number'), variable.name="Gene", value.name="CN")
# rename copy number value names
finalData[Gene == 'PLCG1 Copy Number', Gene := 'PLCG1 (Chr20q)']
finalData[Gene == 'MBD1 Copy Number', Gene := 'MBD1 (Chr18q)']
# rename MAF for plot
finalData[Technology == 'MAF', Technology := 'RAS/BRAF MAF']
################################################################################
########### run correlation ###################################################
# switch statement to re-define p-values
relabel_p <- function(x){
if(x < .0001) return('p < 0.0001')
if(x < .001) return('p < 0.001')
if(x < .01) return('p < 0.01')
if(x < .05) return('p < 0.05')
return(paste0('p = ', signif(x, 3)))
}
# run correlations
MAF_PLCG1 <- cor.test(finalData[Technology == 'RAS/BRAF MAF' & Gene == 'PLCG1 (Chr20q)']$Value,
finalData[Technology == 'RAS/BRAF MAF' & Gene == 'PLCG1 (Chr20q)']$CN)
MAF_PLCG1 <- paste0("r=", signif(MAF_PLCG1$estimate, 2), "\n", relabel_p(MAF_PLCG1$p.value))
delfi_PLCG1 <- cor.test(finalData[Technology == 'DELFI-TF' & Gene == 'PLCG1 (Chr20q)']$Value,
finalData[Technology == 'DELFI-TF' & Gene == 'PLCG1 (Chr20q)']$CN)
delfi_PLCG1 <- paste0("r=", signif(delfi_PLCG1$estimate, 2), "\n", relabel_p(delfi_PLCG1$p.value))
MAF_mbd1 <- cor.test(finalData[Technology == 'RAS/BRAF MAF' & Gene == 'MBD1 (Chr18q)']$Value,
finalData[Technology == 'RAS/BRAF MAF' & Gene == 'MBD1 (Chr18q)']$CN)
MAF_mbd1 <- paste0("r=", signif(MAF_mbd1$estimate, 2), "\n", relabel_p(MAF_mbd1$p.value))
delfi_mbd1 <- cor.test(finalData[Technology == 'DELFI-TF' & Gene == 'MBD1 (Chr18q)']$Value,
finalData[Technology == 'DELFI-TF' & Gene == 'MBD1 (Chr18q)']$CN)
delfi_mbd1 <- paste0("r=", signif(delfi_mbd1$estimate, 2), "\n", relabel_p(delfi_mbd1$p.value))
# make statistics DT
statDT_mbd1 <- data.table(Technology=c('DELFI-TF', 'RAS/BRAF MAF'),
Gene=c('MBD1 (Chr18q)'),
x=c(.12, .15),
y=c(-1),
label=c(delfi_mbd1, MAF_mbd1))
statDT_PLCG1 <- data.table(Technology=c('DELFI-TF', 'RAS/BRAF MAF'),
Gene=c('PLCG1 (Chr20q)'),
x=c(.12, .15),
y=c(.85),
label=c(delfi_PLCG1, MAF_PLCG1))
statDT <- rbindlist(list(statDT_mbd1, statDT_PLCG1))
################################################################################
############## Make final plots ################################################
# construct the plot
FigS5B <- ggplot() +
geom_point(data=finalData, mapping=aes(x=Value, y=CN, color=Technology), alpha=.5, size=3) +
geom_smooth(data=finalData, mapping=aes(x=Value, y=CN, color=Technology), method='glm') +
scale_color_manual("",values=c("DELFI-TF"="#046C9A", "RAS/BRAF MAF"="#D69C4E")) +
scale_x_continuous(labels=scales::percent) +
geom_text(data=statDT, mapping=aes(x=x, y=y, label=label), size=12) +
facet_grid(Gene~Technology, scales="free") +
guides(color = guide_legend(override.aes = list(size = 6))) +
ylab("cfDNA Log2 Copy Number Ratio") +
xlab("Tumor Fraction") +
theme_bw() +
theme(panel.grid=element_blank()) +
theme(strip.text=element_text(color="white", size=42)) +
theme(strip.background = element_rect(fill="#1B465A")) +
theme(legend.position="none") +
theme(axis.text=element_text(size=36), axis.title=element_text(size=42)) +
theme(legend.text=element_text(size=36))
# save the results
png(filename="Figures/FigS5B.png", height=14, width=16, units="in", res=150)
FigS5B
dev.off()
cairo_pdf(file="Figures/FigS5B.pdf", height=14, width=16)
FigS5B
dev.off()