setwd(file.path("/", "local", "groux", "scATAC-seq")) # libraries library(RColorBrewer) # functions source(file.path("scripts", "functions.R")) ################## aggregations around ebf1 motifs ################## # data # open chromatin data.open.1.frag = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_open_bin1bp_fragment.mat"))) data.open.2.frag = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_open_bin2bp_fragment.mat"))) data.open.10.frag = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_open_bin10bp_fragment.mat"))) data.open.1.read = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_open_bin1bp_read.mat"))) data.open.2.read = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_open_bin2bp_read.mat"))) data.open.10.read = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_open_bin10bp_read.mat"))) data.open.1.atac = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_open_bin1bp_read_atac.mat"))) data.open.2.atac = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_open_bin2bp_read_atac.mat"))) data.open.10.atac = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_open_bin10bp_read_atac.mat"))) # mono-nucleosomes data.1nucl.1.frag = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_1nucl_bin1bp_fragment.mat"))) data.1nucl.2.frag = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_1nucl_bin2bp_fragment.mat"))) data.1nucl.10.frag = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_1nucl_bin10bp_fragment.mat"))) data.1nucl.1.read = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_1nucl_bin1bp_read.mat"))) data.1nucl.2.read = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_1nucl_bin2bp_read.mat"))) data.1nucl.10.read = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_1nucl_bin10bp_read.mat"))) data.1nucl.1.cent = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_1nucl_bin1bp_fragment_center.mat"))) data.1nucl.2.cent = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_1nucl_bin2bp_fragment_center.mat"))) data.1nucl.10.cent = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_1nucl_bin10bp_fragment_center.mat"))) # di-nucleosomes data.2nucl.1.frag = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nucl_bin1bp_fragment.mat"))) data.2nucl.2.frag = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nucl_bin2bp_fragment.mat"))) data.2nucl.10.frag = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nucl_bin10bp_fragment.mat"))) data.2nucl.1.read = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nucl_bin1bp_read.mat"))) data.2nucl.2.read = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nucl_bin2bp_read.mat"))) data.2nucl.10.read = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nucl_bin10bp_read.mat"))) data.2nucl.1.cent = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nucl_bin1bp_fragment_center.mat"))) data.2nucl.2.cent = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nucl_bin2bp_fragment_center.mat"))) data.2nucl.10.cent = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nucl_bin10bp_fragment_center.mat"))) # mono-nucleosomes from di-nucleosome data data.nucls.1.frag = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nuclsplitintwo_bin1bp_fragment.mat"))) data.nucls.2.frag = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nuclsplitintwo_bin2bp_fragment.mat"))) data.nucls.10.frag = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nuclsplitintwo_bin10bp_fragment.mat"))) data.nucls.1.read = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nuclsplitintwo_bin1bp_read.mat"))) data.nucls.2.read = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nuclsplitintwo_bin2bp_read.mat"))) data.nucls.10.read = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nuclsplitintwo_bin10bp_read.mat"))) data.nucls.1.cent = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nuclsplitintwo_bin1bp_fragment_center.mat"))) data.nucls.2.cent = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nuclsplitintwo_bin2bp_fragment_center.mat"))) data.nucls.10.cent = as.matrix(read.table(file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_2nuclsplitintwo_bin10bp_fragment_center.mat"))) # colors col = brewer.pal(4, "Set1") # x-axis axis.at.1 = seq(0, ncol(data.open.1.frag), length.out =5) axis.lab.1 = seq(-400, 400, by=200) axis.at.2 = seq(0, ncol(data.open.2.frag), length.out =5) axis.lab.2 = seq(-400, 400, by=200) axis.at.10 = seq(0, ncol(data.open.10.frag), length.out=5) axis.lab.10 = seq(-1000, 1000, by=500) # X11(width=12, height=12) png(filename=file.path("results/10xgenomics_PBMC_5k/ebf1_motifs_10e-6_aggregations.png"), units="in", res=720, width=12, height=9) m = matrix(nrow=4, ncol=4, data=c(16,13,14,15, 10, 1, 4, 7, 11, 2, 5, 8, 12, 3, 6, 9), byrow=T) l = layout(mat=m, widths=c(0.2, 1, 1, 1), heights=c(0.2, 1, 1, 1)) layout.show(l) p = par(mar=c(5.1, 5.1, 4.1, 2.1)) # 1bp resolution ## entire fragments ylim = c(0,max(max(colMeans(data.open.1.frag)), max(colMeans(data.open.1.frag)), max(colMeans(data.1nucl.1.frag)), max(colMeans(data.2nucl.1.frag)), max(colMeans(data.nucls.1.frag)))) plot(colMeans(data.open.1.frag), col=col[1], lwd=3, type='l', main="", xlab="pos[bp]", ylab="Nb of reads", xaxt='n', ylim=ylim, cex.axis=2, cex.lab=2) lines(colMeans(data.open.1.frag), col=col[1], lwd=3) lines(colMeans(data.1nucl.1.frag), col=col[2], lwd=3) lines(colMeans(data.2nucl.1.frag), col=col[3], lwd=3) lines(colMeans(data.nucls.1.frag), col=col[4], lwd=3) axis(side=1, at=axis.at.1, labels=axis.lab.1, cex.axis=1.8) ## entire reads ylim = c(0,max(max(colMeans(data.open.1.read)), max(colMeans(data.open.1.read)), max(colMeans(data.1nucl.1.read)), max(colMeans(data.2nucl.1.read)), max(colMeans(data.nucls.1.read)))) plot(colMeans(data.open.1.read), col=col[1], lwd=3, type='l', main="", xlab="pos[bp]", ylab="Nb of reads", xaxt='n', ylim=ylim, cex.axis=2, cex.lab=2) lines(colMeans(data.1nucl.1.read), col=col[2], lwd=3) lines(colMeans(data.2nucl.1.read), col=col[3], lwd=3) lines(colMeans(data.nucls.1.read), col=col[4], lwd=3) axis(side=1, at=axis.at.1, labels=axis.lab.1, cex.axis=1.8) ## atac reads and centers plot(colMeans(data.open.1.atac)/max(colMeans(data.open.1.atac)), col=col[1], lwd=3, type='l', xaxt='n', main="", xlab="pos[bp]", ylab="Prop max signal", cex.axis=2, cex.lab=2) lines(colMeans(data.1nucl.1.cent)/max(colMeans(data.1nucl.1.cent)), col=col[2], lwd=3) lines(colMeans(data.2nucl.1.cent)/max(colMeans(data.2nucl.1.cent)), col=col[3], lwd=3) lines(colMeans(data.nucls.1.cent)/max(colMeans(data.nucls.1.cent)), col=col[4], lwd=3) axis(side=1, at=axis.at.1, labels=axis.lab.1, cex.axis=1.8) # 2bp resolution ## entire fragments ylim = c(0,max(max(colMeans(data.open.2.frag)), max(colMeans(data.open.2.frag)), max(colMeans(data.1nucl.2.frag)), max(colMeans(data.2nucl.2.frag)), max(colMeans(data.nucls.2.frag)))) plot(colMeans(data.open.2.frag), col=col[1], lwd=3, type='l', main="", xlab="pos[bp]", ylab="Nb of reads", xaxt='n', ylim=ylim, cex.axis=2, cex.lab=2) lines(colMeans(data.1nucl.2.frag), col=col[2], lwd=3) lines(colMeans(data.2nucl.2.frag), col=col[3], lwd=3) lines(colMeans(data.nucls.2.frag), col=col[4], lwd=3) axis(side=1, at=axis.at.2, labels=axis.lab.2, cex.axis=1.8) ## entire reads ylim = c(0,max(max(colMeans(data.open.2.read)), max(colMeans(data.open.2.read)), max(colMeans(data.1nucl.2.read)), max(colMeans(data.2nucl.2.read)), max(colMeans(data.nucls.2.read)))) plot(colMeans(data.open.2.read), col=col[1], lwd=3, type='l', main="", xlab="pos[bp]", ylab="Nb of reads", xaxt='n', ylim=ylim, cex.axis=2, cex.lab=2) lines(colMeans(data.1nucl.2.read), col=col[2], lwd=3) lines(colMeans(data.2nucl.2.read), col=col[3], lwd=3) lines(colMeans(data.nucls.2.read), col=col[4], lwd=3) axis(side=1, at=axis.at.2, labels=axis.lab.2, cex.axis=1.8) ## atac reads and centers plot(colMeans(data.open.2.atac)/max(colMeans(data.open.2.atac)), col=col[1], lwd=3, type='l', xaxt='n', main="", xlab="pos[bp]", ylab="Prop max signal", cex.axis=2, cex.lab=2) lines(colMeans(data.1nucl.2.cent)/max(colMeans(data.1nucl.2.cent)), col=col[2], lwd=3) lines(colMeans(data.2nucl.2.cent)/max(colMeans(data.2nucl.2.cent)), col=col[3], lwd=3) lines(colMeans(data.nucls.2.cent)/max(colMeans(data.nucls.2.cent)), col=col[4], lwd=3) axis(side=1, at=axis.at.2, labels=axis.lab.2, cex.axis=1.8) # 10bp resolution ## entire fragments ylim = c(0,max(max(colMeans(data.open.10.frag)), max(colMeans(data.open.10.frag)), max(colMeans(data.1nucl.10.frag)), max(colMeans(data.2nucl.10.frag)), max(colMeans(data.nucls.10.frag)))) plot(colMeans(data.open.10.frag), col=col[1], lwd=3, type='l', main="", xlab="pos[bp]", ylab="Nb of reads", xaxt='n', ylim=ylim, cex.axis=2, cex.lab=2) lines(colMeans(data.1nucl.10.frag), col=col[2], lwd=3) lines(colMeans(data.2nucl.10.frag), col=col[3], lwd=3) lines(colMeans(data.nucls.10.frag), col=col[4], lwd=3) axis(side=1, at=axis.at.10, labels=axis.lab.10, cex.axis=1.8) ## entire reads ylim = c(0,max(max(colMeans(data.open.10.read)), max(colMeans(data.open.10.read)), max(colMeans(data.1nucl.10.read)), max(colMeans(data.2nucl.10.read)), max(colMeans(data.nucls.10.read)))) plot(colMeans(data.open.10.read), col=col[1], lwd=3, type='l', main="", xlab="pos[bp]", ylab="Nb of reads", xaxt='n', ylim=ylim, cex.axis=2, cex.lab=2) lines(colMeans(data.1nucl.10.read), col=col[2], lwd=3) lines(colMeans(data.2nucl.10.read), col=col[3], lwd=3) lines(colMeans(data.nucls.10.read), col=col[4], lwd=3) axis(side=1, at=axis.at.10, labels=axis.lab.10, cex.axis=1.8) ## atac reads and centers plot(colMeans(data.open.10.atac)/max(colMeans(data.open.10.atac)), col=col[1], lwd=3, type='l', xaxt='n', main="", xlab="pos[bp]", ylab="Prop max signal", cex.axis=2, cex.lab=2) lines(colMeans(data.1nucl.10.cent)/max(colMeans(data.1nucl.10.cent)), col=col[2], lwd=3) lines(colMeans(data.2nucl.10.cent)/max(colMeans(data.2nucl.10.cent)), col=col[3], lwd=3) lines(colMeans(data.nucls.10.cent)/max(colMeans(data.nucls.10.cent)), col=col[4], lwd=3) axis(side=1, at=axis.at.10, labels=axis.lab.10, cex.axis=1.8) # some legends over the rows and columns p = par(mar=c(0,0,0,0)) plot(0, 0, col=0, main="", xlab="", ylab="", xaxt='n', yaxt='n') text(0, 0, labels="FRAGMENTS", cex=2, srt=90) plot(0, 0, col=0, main="", xlab="", ylab="", xaxt='n', yaxt='n') text(0, 0, labels="READS", cex=2, srt=90) plot(0, 0, col=0, main="", xlab="", ylab="", xaxt='n', yaxt='n') text(0, 0, labels="EDGES/CENTERS", cex=2, srt=90) plot(0, 0, col=0, main="", xlab="", ylab="", xaxt='n', yaxt='n') text(0, 0, labels="+/-400bp by 1bp", cex=2) plot(0, 0, col=0, main="", xlab="", ylab="", xaxt='n', yaxt='n') text(0, 0, labels="+/-400bp by 2bp", cex=2) plot(0, 0, col=0, main="", xlab="", ylab="", xaxt='n', yaxt='n') text(0, 0, labels="+/-1kp by 10bp", cex=2) par(p) dev.off() # footprint # x-axis axis.lab.1 = seq(-200, 200, by=100) axis.at.1 = seq(0, 400, length.out=length(axis.lab.1)) axis.lab.2 = seq(-200, 200, by=100) axis.at.2 = seq(0, 200, length.out=length(axis.lab.2)) axis.lab.10 = seq(-200, 200, by=100) axis.at.10 = seq(0, 41, length.out=length(axis.lab.10)) # X11(width=10, height=12) png(filename=file.path("results", "10xgenomics_PBMC_5k", "ebf1_motifs_10e-6_footprint.png"), units="in", res=720, width=10, height=12) p = par(mfrow=c(3,1), mar=c(5.1, 5.1, 4.1, 2.1)) # 1bp resolution index = 200:600 x = 1:length(index) plot(x, colMeans(data.open.1.atac[,index])/max(colMeans(data.open.1.atac[,index])), type='l', lwd=3, col=col[1], main="EBF1 motif 1bp", xlab="pos[bp]", ylab="Prop max signal", xaxt='n', cex.axis=2, cex.lab=2, cex.main=2) lines(x, colMeans(data.1nucl.1.cent[,index])/max(colMeans(data.1nucl.1.cent[,index])), lwd=3, col=col[2]) lines(x, colMeans(data.nucls.1.cent[,index])/max(colMeans(data.nucls.1.cent[,index])), lwd=3, col=col[4]) abline(v=191, lwd=3, lty=2) abline(v=211, lwd=3, lty=2) axis(side=1, at=axis.at.1, labels=axis.lab.1, cex.axis=1.8) # 2bp resolution index = 100:300 x = 1:length(index) plot(x, colMeans(data.open.2.atac[,index])/max(colMeans(data.open.2.atac[,index])), type='l', lwd=3, col=col[1], main="EBF1 motif 2bp", xlab="pos[bp]", ylab="Prop max signal", xaxt='n', cex.axis=2, cex.lab=2, cex.main=2) lines(x, colMeans(data.1nucl.2.cent[,index])/max(colMeans(data.1nucl.2.cent[,index])), lwd=3, col=col[2]) lines(x, colMeans(data.nucls.2.cent[,index])/max(colMeans(data.nucls.2.cent[,index])), lwd=3, col=col[4]) abline(v=96, lwd=3, lty=2) abline(v=106, lwd=3, lty=2) axis(side=1, at=axis.at.1, labels=axis.lab.1, cex.axis=1.8) # 10bp resolution index = 80:120 x = 1:length(index) plot(x, colMeans(data.open.10.atac[,index])/max(colMeans(data.open.10.atac[,index])), type='l', lwd=3, col=col[1], main="EBF1 motif 10bp", xlab="pos[bp]", ylab="Prop max signal", xaxt='n', cex.axis=2, cex.lab=2, cex.main=2) lines(x, colMeans(data.1nucl.10.cent[,index])/max(colMeans(data.1nucl.10.cent[,index])), lwd=3, col=col[2]) lines(x, colMeans(data.nucls.10.cent[,index])/max(colMeans(data.nucls.10.cent[,index])), lwd=3, col=col[4]) abline(v=20, lwd=3, lty=2) abline(v=22, lwd=3, lty=2) axis(side=1, at=axis.at.10, labels=axis.lab.10, cex.axis=1.8) par(p) dev.off()