Homec4science
Diffusion PhD_thesis (master) Browse

R9524/scripts/ch_atac-seq/figure_sp1_6classes.R
9c4e9be36fc2master

Last Change
Raw File

figure_sp1_6classes.R
View Options

setwd(file.path("/", "local", "groux", "scATAC-seq"))
# libraries
library(RColorBrewer)
# functions
source(file.path("scripts", "functions.R"))
# path to the images for the logo
path.a = file.path("res/A.png")
path.c = file.path("res/C.png")
path.g = file.path("res/G.png")
path.t = file.path("res/T.png")
# paths
data.dir = file.path("/", "local", "groux", "scATAC-seq", "results")
dest.dir = file.path("/", "local", "groux", "phd_thesis", "images", "ch_atac-seq")
# colors
col = brewer.pal(3, "Set1")
##################################### sp1 with flip only #####################################
# open chromatin
data = read.read.models(file.path(data.dir,
"10xgenomics_PBMC_5k_motifs_classification_0",
"sp1_motifs_10e-7_open_bin1bp_read_atac_6class_model.mat"))
model.open = data$models
model.prob = data$prob
data = NULL
# nucleosomes
model.nucl = read.read.models(file.path(data.dir,
"10xgenomics_PBMC_5k_motifs_classification_0",
"sp1_motifs_10e-7_open_bin1bp_read_atac_6class_nucleosomes_fragment_center_model.mat"))$models
# sequence
model.seq = read.sequence.models(file.path(data.dir,
"10xgenomics_PBMC_5k_motifs_classification_0",
"sp1_motifs_10e-7_open_bin1bp_read_atac_6class_sequences_model.mat"))$models
png(filename=file.path(dest.dir,
"sp1_motifs_6class_noshift_flip.png"),
units="in", res=480, width=18, height=8)
m = matrix(1:6, nrow=3, ncol=2, byrow=F)
layout(m)
# order from most to least probable class
ord = order(model.prob, decreasing=T)
ref.open = model.open[ord,, drop=F]
ref.nucl = model.nucl[ord,, drop=F]
ref.seq = model.seq[,,ord, drop=F]
prob = model.prob[ord]
class = c(1:nrow(ref.open))[ord]
for(i in 1:nrow(ref.open))
{ # plot logo
plot.logo(ref.seq[,,i], path.a, path.c, path.g, path.t,
main=sprintf("class %d (p=%.2f)", class[i], prob[i]))
# x-axis
x.lab = seq(-ncol(ref.open), ncol(ref.open), length.out=3)
x.at = (x.lab + ncol(ref.open)) / 2
axis(1, at=x.at, labels=x.lab)
# y-axis is [0,2] because these are bits but
# label it [0,1] for min/max signal
x.at = seq(0, 2, 1)
axis(2, at=x.at, labels=0.5*x.at)
# plot signal (multiplies by 2 because the y-axis goes to 2 bits)
lines(2*(ref.open[i,] / max(ref.open[i,])), lwd=2, col=col[1])
lines(2*(ref.nucl[i,] / max(ref.nucl[i,])), lwd=2, col=col[2])
}
row_n = 1 # row counter
col_n = 1 # column counter
for(i in 1:nrow(ref.open))
{ # plot logo center
right = 0.5*col_n - 0.01
left = right - 0.2
bottom = 1-(row_n*(1/3))+0.1
top = bottom + 0.2
par(fig=c(left, right, bottom, top), new=T)
idx = 380:420
plot.logo(ref.seq[,idx,i], path.a, path.c, path.g, path.t)
# plot signal (multiplies by 2 because the y-axis goes to 2 bits)
lines(2*(ref.open[i,idx] / max(ref.open[i,])), lwd=2, col=col[1])
lines(2*(ref.nucl[i,idx] / max(ref.nucl[i,])), lwd=2, col=col[2])
# xaxis
x.at = 1:length(idx)
axis(1, at=x.at, labels=x.at)
# yaxis
x.at = seq(0, 2, by=1)
axis(2, at=x.at, labels=x.at)
row_n = row_n + 1
if(i %% 3 == 0)
{ col_n = col_n + 1
row_n = 1
}
}
dev.off()
##################################### sp1 with flip and shift #####################################
# open chromatin
data = read.read.models(file.path(data.dir,
"10xgenomics_PBMC_5k_motifs_classification_1",
"sp1_motifs_10e-7_open_bin1bp_read_atac_6class_model.mat"))
model.open = data$models
model.prob = data$prob
data = NULL
# nucleosomes
model.nucl = read.read.models(file.path(data.dir,
"10xgenomics_PBMC_5k_motifs_classification_1",
"sp1_motifs_10e-7_open_bin1bp_read_atac_6class_1nucl_fragment_center_model.mat"))$models
# sequence
model.seq = read.sequence.models(file.path(data.dir,
"10xgenomics_PBMC_5k_motifs_classification_1",
"sp1_motifs_10e-7_open_bin1bp_read_atac_6class_sequences_model.mat"))$models
png(filename=file.path(dest.dir,
"sp1_motifs_6class_shift_flip.png"),
units="in", res=480, width=18, height=8)
m = matrix(1:6, nrow=3, ncol=2, byrow=F)
layout(m)
# order from most to least probable class
ord = order(model.prob, decreasing=T)
ref.open = model.open[ord,, drop=F]
ref.nucl = model.nucl[ord,, drop=F]
ref.seq = model.seq[,,ord, drop=F]
prob = model.prob[ord]
class = c(1:nrow(ref.open))[ord]
for(i in 1:nrow(ref.open))
{ # plot logo
plot.logo(ref.seq[,,i], path.a, path.c, path.g, path.t,
main=sprintf("class %d (p=%.2f)", class[i], prob[i]))
# x-axis
x.lab = seq(-ncol(ref.open), ncol(ref.open), length.out=3)
x.at = (x.lab + ncol(ref.open)) / 2
axis(1, at=x.at, labels=x.lab)
# y-axis is [0,2] because these are bits but
# label it [0,1] for min/max signal
x.at = seq(0, 2, 1)
axis(2, at=x.at, labels=0.5*x.at)
# plot signal (multiplies by 2 because the y-axis goes to 2 bits)
lines(2*(ref.open[i,] / max(ref.open[i,])), lwd=2, col=col[1])
lines(2*(ref.nucl[i,] / max(ref.nucl[i,])), lwd=2, col=col[2])
}
row_n = 1 # row counter
col_n = 1 # column counter
for(i in 1:nrow(ref.open))
{ # plot logo center
right = 0.5*col_n - 0.01
left = right - 0.2
bottom = 1-(row_n*(1/3))+0.1
top = bottom + 0.2
par(fig=c(left, right, bottom, top), new=T)
idx = 380:420
plot.logo(ref.seq[,idx,i], path.a, path.c, path.g, path.t)
# plot signal (multiplies by 2 because the y-axis goes to 2 bits)
lines(2*(ref.open[i,idx] / max(ref.open[i,])), lwd=2, col=col[1])
lines(2*(ref.nucl[i,idx] / max(ref.nucl[i,])), lwd=2, col=col[2])
# xaxis
x.at = 1:length(idx)
axis(1, at=x.at, labels=x.at)
# yaxis
x.at = seq(0, 2, by=1)
axis(2, at=x.at, labels=x.at)
row_n = row_n + 1
if(i %% 3 == 0)
{ col_n = col_n + 1
row_n = 1
}
}
dev.off()
c4science · Help