setwd(file.path("/", "local", "groux", "scATAC-seq")) # libraries library(RColorBrewer) library(seqLogo) # functions source(file.path("scripts", "functions.R")) # the minimum number of classes searched k.min = 1 # the maximum number of classes searched k.max = 10 # 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") ################## sequence patterns around sp1 motifs ################## for(k in k.min:k.max) { # open chromatin data = read.read.models(file.path("results", "10xgenomics_PBMC_5k_classification_3", sprintf("sp1_motifs_10e-7_open_bin1bp_read_atac_%dclass_model.mat", k))) model.open = data$models model.prob = data$prob data = NULL # nucleosomes model.nucl = read.read.models(file.path("results", "10xgenomics_PBMC_5k_classification_3", sprintf("sp1_motifs_10e-7_nucleosomes_bin1bp_fragment_center_%dclass_model.mat", k)))$models # sequence model.seq = read.sequence.models(file.path("results", "10xgenomics_PBMC_5k_classification_3", sprintf("sp1_motifs_10e-7_sequences_%dclass_model.mat", k)))$models # plot classes col = brewer.pal(3, "Set1") # X11(width=17, height=10) png(filename=file.path("results", "10xgenomics_PBMC_5k_classification_3", sprintf("sp1_motifs_10e-7_classification_sequences_%dclass.png", k)), units="in", res=720, width=18, height=12) m = matrix(1:10, nrow=5, 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,1] for min/max signal x.at = seq(0, 1, 0.5) axis(2, at=x.at, labels=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=1, col=col[1]) lines(2*(ref.nucl[i,] / max(ref.nucl[i,])), lwd=1, 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*(0.2))+0.05 top = bottom + 0.15 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=1, col=col[1]) lines(2*(ref.nucl[i,idx] / max(ref.nucl[i,])), lwd=1, 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 %% 5 == 0) { col_n = col_n + 1 row_n = 1 } } dev.off() }