Homec4science
Diffusion PhD_thesis (master) Browse

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

Last Change
Raw File

figure_ctcf_sp1_myc_ebf1_footprint.R
View Options

setwd(file.path("/", "local", "groux", "scATAC-seq"))
# libraries
library(RColorBrewer)
# functions
source(file.path("scripts", "functions.R"))
data.dir = file.path("/", "local", "groux", "scATAC-seq", "data", "10xgenomics_PBMC_5k_motifs")
dest.dir = file.path("/", "local", "groux", "phd_thesis", "scATAC-seq", "figures")
# CTCF data
## open chromatin
ctcf.open.1.atac = as.matrix(read.table(file.path(data.dir, "ctcf_motifs_10e-6_open_bin1bp_read_atac.mat")))
## nucleosomes
ctcf.nucl.1.cent = as.matrix(read.table(file.path(data.dir, "ctcf_motifs_10e-6_nucleosomes_bin1bp_fragment_center.mat")))
ctcf.open.1.atac = colMeans(ctcf.open.1.atac) / max(colMeans(ctcf.open.1.atac))
ctcf.nucl.1.cent = colMeans(ctcf.nucl.1.cent) / max(colMeans(ctcf.nucl.1.cent))
# SP1 data
## open chromatin
sp1.open.1.atac = as.matrix(read.table(file.path(data.dir, "sp1_motifs_10e-7_open_bin1bp_read_atac.mat")))
## nucleosomes
sp1.nucl.1.cent = as.matrix(read.table(file.path(data.dir, "sp1_motifs_10e-7_nucleosomes_bin1bp_fragment_center.mat")))
sp1.open.1.atac = colMeans(sp1.open.1.atac) / max(colMeans(sp1.open.1.atac))
sp1.nucl.1.cent = colMeans(sp1.nucl.1.cent) / max(colMeans(sp1.nucl.1.cent))
# myc data
## open chromatin
myc.open.1.atac = as.matrix(read.table(file.path(data.dir, "myc_motifs_10e-6_open_bin1bp_read_atac.mat")))
## nucleosomes
myc.nucl.1.cent = as.matrix(read.table(file.path(data.dir, "myc_motifs_10e-6_nucleosomes_bin1bp_fragment_center.mat")))
myc.open.1.atac = colMeans(myc.open.1.atac) / max(colMeans(myc.open.1.atac))
myc.nucl.1.cent = colMeans(myc.nucl.1.cent) / max(colMeans(myc.nucl.1.cent))
# EBF1 data
## open chromatin
ebf1.open.1.atac = as.matrix(read.table(file.path(data.dir, "ebf1_motifs_10e-6_open_bin1bp_read_atac.mat")))
## nucleosomes
ebf1.nucl.1.cent = as.matrix(read.table(file.path(data.dir, "ebf1_motifs_10e-6_nucleosomes_bin1bp_fragment_center.mat")))
ebf1.open.1.atac = colMeans(ebf1.open.1.atac) / max(colMeans(ebf1.open.1.atac))
ebf1.nucl.1.cent = colMeans(ebf1.nucl.1.cent) / max(colMeans(ebf1.nucl.1.cent))
# colors
col = brewer.pal(4, "Set1")
# display center only
idx = 200:600
# x-axis
axis.at = seq(-200, 200, length.out=3) + 200
axis.lab = seq(-200, 200, by=200)
y.lim = c(0,1)
# X11(width=18, height=9)
png(filename=file.path(dest.dir, "ctcf_sp1_myc_ebf1_footprint.png"),
units="in", res=480, width=18, height=9)
m = matrix(nrow=2, ncol=2,
data=c(1,3,
2,4),
byrow=T)
l = layout(mat=m, widths=c(1,1), heights=c(1,1))
p = par(mar=c(5.1, 5.1, 4.1, 2.1))
# CTCF
plot(ctcf.open.1.atac[idx], col=col[1], lwd=3, type='l',
main="CTCF motif", xlab="pos[bp]", ylab="Prop max. signal", xaxt='n', ylim=y.lim,
cex.axis=2, cex.lab=2, cex.main=1.8)
lines(ctcf.nucl.1.cent[idx], col=col[2], lwd=3)
axis(side=1, at=axis.at, labels=axis.lab, cex.axis=1.8)
abline(v=191, lty=2, lwd=3)
abline(v=210, lty=2, lwd=3)
# SP1
plot(sp1.open.1.atac[idx], col=col[1], lwd=3, type='l',
main="SP1 motif", xlab="pos[bp]", ylab="Prop max. signal", xaxt='n', ylim=y.lim,
cex.axis=2, cex.lab=2, cex.main=1.8)
lines(sp1.nucl.1.cent[idx], col=col[2], lwd=3)
axis(side=1, at=axis.at, labels=axis.lab, cex.axis=1.8)
abline(v=194, lty=2, lwd=3)
abline(v=210, lty=2, lwd=3)
# myc
plot(myc.open.1.atac[idx], col=col[1], lwd=3, type='l',
main="myc motif", xlab="pos[bp]", ylab="Prop max. signal", xaxt='n', ylim=y.lim,
cex.axis=2, cex.lab=2, cex.main=1.8)
lines(myc.nucl.1.cent[idx], col=col[2], lwd=3)
axis(side=1, at=axis.at, labels=axis.lab, cex.axis=1.8)
abline(v=195, lty=2, lwd=3)
abline(v=209, lty=2, lwd=3)
# EBF1
plot(ebf1.open.1.atac[idx], col=col[1], lwd=3, type='l',
main="EBF1 motif", xlab="pos[bp]", ylab="Prop max. signal", xaxt='n', ylim=y.lim,
cex.axis=2, cex.lab=2, cex.main=1.8)
lines(ebf1.nucl.1.cent[idx], col=col[2], lwd=3)
axis(side=1, at=axis.at, labels=axis.lab, cex.axis=1.8)
abline(v=197, lty=2, lwd=3)
abline(v=206, lty=2, lwd=3)
dev.off()
c4science · Help