Homec4science
Diffusion scATAC-seq (master) Browse

scATAC-seq/scripts/10xgenomics_PBMC_5k_peaks_classification_2/analysis_test_sampled.R
1bfc0a93ba65master

Last Change
Raw File

analysis_test_sampled.R
View Options

setwd(file.path("/", "local", "groux", "scATAC-seq"))
# libraries
library(RColorBrewer)
# functions
source(file.path("scripts", "functions.R"))
# the number of classes searched
n.classes = c(10, 20, 30)
# 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 ctcf motifs ##################
for(k in n.classes)
{
# sequence
data = read.sequence.models(file.path("results", "10xgenomics_PBMC_5k_peaks_classification_2",
sprintf("peaks_rmsk_sampled_sequences_1kb_%dclass_model.mat", k)))
model.seq = data$models
model.prob = data$prob
data = NULL
# open chromatin
model.open = read.read.models(file.path("results", "10xgenomics_PBMC_5k_peaks_classification_2",
sprintf("peaks_rmsk_sampled_openchromatin_1kb_read_atac_%dclass_model.mat", k)))$models
# nucleosomes
model.nucl = read.read.models(file.path("results", "10xgenomics_PBMC_5k_peaks_classification_2",
sprintf("peaks_rmsk_sampled_nucleosomes_1kb_fragment_center_%dclass_model.mat", k)))$models
# plot classes
col = brewer.pal(3, "Set1")
# X11(width=26, height=12)
png(filename=file.path("results", "10xgenomics_PBMC_5k_peaks_classification_2",
sprintf("peaks_rmsk_sampled_sequences_%dclass.png", k)),
units="in", res=720, width=18, height=12)
m = matrix(1:30, nrow=6, ncol=5, 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)-1)/2, (ncol(ref.open)-1)/2, length.out=3)
x.at = seq(1, ncol(ref.open), length.out=length(x.lab))
axis(1, at=x.at, labels=x.lab)
# y-axis is [0,1] for min/max signal
y.at = seq(0, 2, length.out=2)
y.lab = c("min", "max")
axis(2, at=y.at, labels=y.lab)
# 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])
}
# inlets with center
# 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 = (391-1-20):(391+1+20)
# 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 = seq(1, length(idx), length.out = 3)
# x.lab = seq(-(ncol(ref.open)-1)/2, (ncol(ref.open)-1)/2)[idx][x.at]
# axis(1, at=x.at, labels=x.lab)
# # yaxis
# axis(2, at=y.at, labels=y.lab)
# row_n = row_n + 1
# if(i %% 5 == 0)
# { col_n = col_n + 1
# row_n = 1
# }
# }
dev.off()
}
c4science · Help