Homec4science
Diffusion PhD_thesis (master) Browse

R9524/scripts/ch_smile-seq/recreate_figure_s4b.R
0170f77735f4master

Last Change
Raw File

recreate_figure_s4b.R
View Options

# libraries
library(RColorBrewer)
library(gplots)
library(ggplot2)
library(gridExtra)
# paths
data.dir = "/local/groux/SMILe_seq/SMILe_seq_results/PWM_evaluation_25000_3"
res.dir = "/local/groux/phd_thesis/images/ch_smile-seq/"
source(file.path(data.dir, "scripts", "functions.R"))
results.25000.1.500 = file.path(data.dir, "PWM_evaluation_25000_3_1-500")
results.25000.501.1000 = file.path(data.dir, "PWM_evaluation_25000_3_501-1000")
results.25000.1001.1500 = file.path(data.dir, "PWM_evaluation_25000_3_1001-1500")
results.25000.1501.2000 = file.path(data.dir, "PWM_evaluation_25000_3_1501-2000")
results.25000.2001.2500 = file.path(data.dir, "PWM_evaluation_25000_3_2001-2500")
results.25000.2501.3000 = file.path(data.dir, "PWM_evaluation_25000_3_2501-3000")
results.25000.3001.3500 = file.path(data.dir, "PWM_evaluation_25000_3_3001-3500")
results.25000.3501.4000 = file.path(data.dir, "PWM_evaluation_25000_3_3501-4000")
results.25000.4001.4500 = file.path(data.dir, "PWM_evaluation_25000_3_4001-4500")
results.25000.4501.5000 = file.path(data.dir, "PWM_evaluation_25000_3_4501-5000")
results.25000.5001.5500 = file.path(data.dir, "PWM_evaluation_25000_3_5001-5500")
results.25000.5501.6000 = file.path(data.dir, "PWM_evaluation_25000_3_5501-6000")
# reading data
DATA.25000.1.500 = read.table(file.path(results.25000.1.500, "AUC.tsv"), header=T, stringsAsFactors=F)
DATA.25000.501.1000 = read.table(file.path(results.25000.501.1000, "AUC.tsv"), header=T, stringsAsFactors=F)
DATA.25000.1001.1500 = read.table(file.path(results.25000.1001.1500, "AUC.tsv"), header=T, stringsAsFactors=F)
DATA.25000.1501.2000 = read.table(file.path(results.25000.1501.2000, "AUC.tsv"), header=T, stringsAsFactors=F)
DATA.25000.2001.2500 = read.table(file.path(results.25000.2001.2500, "AUC.tsv"), header=T, stringsAsFactors=F)
DATA.25000.2501.3000 = read.table(file.path(results.25000.2501.3000, "AUC.tsv"), header=T, stringsAsFactors=F)
DATA.25000.3001.3500 = read.table(file.path(results.25000.3001.3500, "AUC.tsv"), header=T, stringsAsFactors=F)
DATA.25000.3501.4000 = read.table(file.path(results.25000.3501.4000, "AUC.tsv"), header=T, stringsAsFactors=F)
DATA.25000.4001.4500 = read.table(file.path(results.25000.4001.4500, "AUC.tsv"), header=T, stringsAsFactors=F)
DATA.25000.4501.5000 = read.table(file.path(results.25000.4501.5000, "AUC.tsv"), header=T, stringsAsFactors=F)
DATA.25000.5001.5500 = read.table(file.path(results.25000.5001.5500, "AUC.tsv"), header=T, stringsAsFactors=F)
DATA.25000.5501.6000 = read.table(file.path(results.25000.5501.6000, "AUC.tsv"), header=T, stringsAsFactors=F)
# reformat data into matrices, add a mean column
info_col = c(1:3)
data_col = c(4:8)
INFO = df2matrix(DATA.25000.1.500, info_col)
DATA.25000.1.500 = df2matrix(DATA.25000.1.500, data_col)
DATA.25000.1.500 = matrix.add.mean(DATA.25000.1.500)
DATA.25000.501.1000 = df2matrix(DATA.25000.501.1000, data_col)
DATA.25000.501.1000 = matrix.add.mean(DATA.25000.501.1000)
DATA.25000.1001.1500 = df2matrix(DATA.25000.1001.1500, data_col)
DATA.25000.1001.1500 = matrix.add.mean(DATA.25000.1001.1500)
DATA.25000.1501.2000 = df2matrix(DATA.25000.1501.2000, data_col)
DATA.25000.1501.2000 = matrix.add.mean(DATA.25000.1501.2000)
DATA.25000.2001.2500 = df2matrix(DATA.25000.2001.2500, data_col)
DATA.25000.2001.2500 = matrix.add.mean(DATA.25000.2001.2500)
DATA.25000.2501.3000 = df2matrix(DATA.25000.2501.3000, data_col)
DATA.25000.2501.3000 = matrix.add.mean(DATA.25000.2501.3000)
DATA.25000.3001.3500 = df2matrix(DATA.25000.3001.3500, data_col)
DATA.25000.3001.3500 = matrix.add.mean(DATA.25000.3001.3500)
DATA.25000.3501.4000 = df2matrix(DATA.25000.3501.4000, data_col)
DATA.25000.3501.4000 = matrix.add.mean(DATA.25000.3501.4000)
DATA.25000.4001.4500 = df2matrix(DATA.25000.4001.4500, data_col)
DATA.25000.4001.4500 = matrix.add.mean(DATA.25000.4001.4500)
DATA.25000.4501.5000 = df2matrix(DATA.25000.4501.5000, data_col)
DATA.25000.4501.5000 = matrix.add.mean(DATA.25000.4501.5000)
DATA.25000.5001.5500 = df2matrix(DATA.25000.5001.5500, data_col)
DATA.25000.5001.5500 = matrix.add.mean(DATA.25000.5001.5500)
DATA.25000.5501.6000 = df2matrix(DATA.25000.5501.6000, data_col)
DATA.25000.5501.6000 = matrix.add.mean(DATA.25000.5501.6000)
# some data rows are only NA thus the mean is NaN, change these NaN
# to 0 to avoid later computation issues.
DATA.25000.1.500 = matrix.NaN2n(DATA.25000.1.500, col=6, 0)
DATA.25000.501.1000 = matrix.NaN2n(DATA.25000.501.1000, col=6, 0)
DATA.25000.1001.1500 = matrix.NaN2n(DATA.25000.1001.1500, col=6, 0)
DATA.25000.1501.2000 = matrix.NaN2n(DATA.25000.1501.2000, col=6, 0)
DATA.25000.2001.2500 = matrix.NaN2n(DATA.25000.2001.2500, col=6, 0)
DATA.25000.2501.3000 = matrix.NaN2n(DATA.25000.2501.3000, col=6, 0)
DATA.25000.3001.3500 = matrix.NaN2n(DATA.25000.3001.3500, col=6, 0)
DATA.25000.3501.4000 = matrix.NaN2n(DATA.25000.3501.4000, col=6, 0)
DATA.25000.4001.4500 = matrix.NaN2n(DATA.25000.4001.4500, col=6, 0)
DATA.25000.4501.5000 = matrix.NaN2n(DATA.25000.4501.5000, col=6, 0)
DATA.25000.5001.5500 = matrix.NaN2n(DATA.25000.5001.5500, col=6, 0)
DATA.25000.5501.6000 = matrix.NaN2n(DATA.25000.5501.6000, col=6, 0)
# colors
colors = brewer.pal(n=5, name="Spectral")
colors6 = brewer.pal(n=6, name="Spectral")
# data structures
PEAK.RANGES = c("1-500", "501-1000", "1001-1500", "1501-2000", "2001-2500", "2501-3000", "3001-3500", "3501-4000", "4001-4500", "4501-5000",
"5001-5500", "5501-6000")
PWM.TYPES = c("SMILe-seq", "HT-Selex", "Uniprobe", "Jaspar", "HOCOMOCO", "mean")
NAMES = c("1", "501", "1001", "1501", "2001", "2501", "3001", "3501", "4001", "4501", "5001", "5501")
DATA.LIST = list(DATA.25000.1.500, DATA.25000.501.1000, DATA.25000.1001.1500, DATA.25000.1501.2000, DATA.25000.2001.2500,
DATA.25000.2501.3000, DATA.25000.3001.3500, DATA.25000.3501.4000, DATA.25000.4001.4500, DATA.25000.4501.5000,
DATA.25000.5001.5500, DATA.25000.5501.6000)
DATA = list2array(DATA.LIST, INFO[,3], PEAK.RANGES, PWM.TYPES)
DATA.LIST = INFO.LIST = NULL
# boxplots of AUCs ordered by peak interval (all 1-500, all 501-1000, etc)
#########
# CEBPb #
#########
f = "CEBPb"
index = which(dimnames(DATA)[[1]] == f)
n = length(index)
# create dataframe
SUBDF = matrix2ggplot2.df(DATA[index[1],,])
for(i in 2:length(index))
{ SUBDF = rbind(SUBDF, matrix2ggplot2.df(DATA[index[i],,])) }
# remove lines labelled 'mean' for method
SUBDF = SUBDF[-which(SUBDF$methods == "mean"),]
SUBDF$inter = interaction(SUBDF$intervals, SUBDF$methods, lex.order=T)
SUBDF$xlab = as.factor(unlist(strsplit(as.character(SUBDF$inter), split='\\.'))[c(T,F)])
SUBDF$xlab = factor(SUBDF$xlab, levels(SUBDF$xlab)[c(1,11,2,3,4,5,6,7,8,9,10,12)])
SUBDF = SUBDF[-c(which(SUBDF$methods == levels(SUBDF$methods)[3])),] # remove UniProbe
# color vector
color = rep(colors[-3], 12)
# boxplots
p.cebpb = ggplot(SUBDF, aes(x=xlab, y=AUC, fill=inter)) + geom_boxplot()
p.cebpb = p.cebpb + scale_fill_manual(values=color)
p.cebpb = p.cebpb + theme(plot.title=element_text(hjust=0.5),
axis.text.x=element_text(angle=0, hjust=0.5),
axis.title.y=element_text(face="bold")) + ggtitle(sprintf("%s (%d ENCODE datasets)", f, n)) + labs(x="",y="AUC") + guides(fill=FALSE)
p.cebpb = p.cebpb + geom_vline(xintercept=seq(from=1.5, to=11.5, by=1), linetype="dashed")
SUBDF = color_tmp = n = NULL
########
# CTCF #
########
f = "CTCF"
index = which(dimnames(DATA)[[1]] == f)
n = length(index)
# create dataframe
SUBDF = matrix2ggplot2.df(DATA[index[1],,])
for(i in 2:length(index))
{ SUBDF = rbind(SUBDF, matrix2ggplot2.df(DATA[index[i],,])) }
# remove lines labelled 'mean' for method
SUBDF = SUBDF[-which(SUBDF$methods == "mean"),]
SUBDF$inter = interaction(SUBDF$intervals, SUBDF$methods, lex.order=T)
SUBDF$xlab = as.factor(unlist(strsplit(as.character(SUBDF$inter), split='\\.'))[c(T,F)])
SUBDF$xlab = factor(SUBDF$xlab, levels(SUBDF$xlab)[c(1,11,2,3,4,5,6,7,8,9,10,12)])
SUBDF = SUBDF[-c(which(SUBDF$methods == levels(SUBDF$methods)[3])),] # remove UniProbe
# color vector
color = rep(colors[-3], 12)
# boxplots
p.ctcf = ggplot(SUBDF, aes(x=xlab, y=AUC, fill=inter)) + geom_boxplot()
p.ctcf = p.ctcf + scale_fill_manual(values=color)
p.ctcf = p.ctcf + theme(plot.title=element_text(hjust=0.5),
axis.text.x=element_text(angle=0, hjust=0.5),
axis.title.y=element_text(face="bold")) + ggtitle(sprintf("%s (%d ENCODE datasets)", f, n)) + labs(x="",y="AUC") + guides(fill=FALSE)
p.ctcf = p.ctcf + geom_vline(xintercept=seq(from=1.5, to=11.5, by=1), linetype="dashed")
SUBDF = color_tmp = n = NULL
#########
# TCF7 #
########
f = "TCF7"
index = which(dimnames(DATA)[[1]] == f)
n = length(index)
# create dataframe
SUBDF = matrix2ggplot2.df(DATA[index[1],,])
for(i in 2:length(index))
{ SUBDF = rbind(SUBDF, matrix2ggplot2.df(DATA[index[i],,])) }
# remove lines labelled 'mean' for method
SUBDF = SUBDF[-which(SUBDF$methods == "mean"),]
SUBDF$inter = interaction(SUBDF$intervals, SUBDF$methods, lex.order=T)
SUBDF$xlab = as.factor(unlist(strsplit(as.character(SUBDF$inter), split='\\.'))[c(T,F)])
SUBDF$xlab = factor(SUBDF$xlab, levels(SUBDF$xlab)[c(1,11,2,3,4,5,6,7,8,9,10,12)])
SUBDF = SUBDF[-c(which(SUBDF$methods == levels(SUBDF$methods)[3])),] # remove UniProbe
# color vector
color = rep(colors[-3], 12)
# boxplots
p.tcf7 = ggplot(SUBDF, aes(x=xlab, y=AUC, fill=inter)) + geom_boxplot()
p.tcf7 = p.tcf7 + scale_fill_manual(values=color)
p.tcf7 = p.tcf7 + theme(plot.title=element_text(hjust=0.5),
axis.text.x=element_text(angle=0, hjust=0.5),
axis.title.y=element_text(face="bold")) + ggtitle(sprintf("%s (%d ENCODE datasets)", f, n)) + labs(x="",y="AUC") + guides(fill=FALSE)
p.tcf7 = p.tcf7 + geom_vline(xintercept=seq(from=1.5, to=11.5, by=1), linetype="dashed")
png(filename=file.path(res.dir, "figure_s4b_recreate.png"),
units="in", res=720, width=12, height=10)
grid.arrange(p.cebpb, p.ctcf, p.tcf7, nrow=3)
dev.off()
c4science · Help