Homec4science
Diffusion scATAC-seq (master) Browse

scATAC-seq/scripts/bulk_sequencing/analysis_cluster_ctcf_mnase_k562.R
1bfc0a93ba65master

Last Change
Raw File

analysis_cluster_ctcf_mnase_k562.R
View Options

setwd(file.path("/", "local", "groux", "scATAC-seq"))
# libraries
library(RColorBrewer)
# functions
source(file.path("scripts", "functions.R"))
# data
data.1 = read.read.models(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_1class_ref.mat"))
ref.1 = data.1$models
prob.1 = data.1$prob
# aic.1 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_1class_aic.txt")))
data.1 = NULL
data.2 = read.read.models(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_2class_ref.mat"))
ref.2 = data.2$models
prob.2 = data.2$prob
# aic.2 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_2class_aic.txt")))
data.2 = NULL
data.3 = read.read.models(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_3class_ref.mat"))
ref.3 = data.3$models
prob.3 = data.3$prob
# aic.3 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_3class_aic.txt")))
data.3 = NULL
data.4 = read.read.models(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_4class_ref.mat"))
ref.4 = data.4$models
prob.4 = data.4$prob
# aic.4 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_4class_aic.txt")))
data.4 = NULL
data.5 = read.read.models(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_5class_ref.mat"))
ref.5 = data.5$models
prob.5 = data.5$prob
# aic.5 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_5class_aic.txt")))
data.5 = NULL
data.6 = read.read.models(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_6class_ref.mat"))
ref.6 = data.6$models
prob.6 = data.6$prob
# aic.6 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_6class_aic.txt")))
data.6 = NULL
data.7 = read.read.models(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_7class_ref.mat"))
ref.7 = data.7$models
prob.7 = data.7$prob
# aic.7 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_7class_aic.txt")))
data.7 = NULL
data.8 = read.read.models(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_8class_ref.mat"))
ref.8 = data.8$models
prob.8 = data.8$prob
# aic.8 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_8class_aic.txt")))
data.8 = NULL
data.9 = read.read.models(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_9class_ref.mat"))
ref.9 = data.9$models
prob.9 = data.9$prob
# aic.9 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_9class_aic.txt")))
data.9 = NULL
data.10 = read.read.models(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_10class_ref.mat"))
ref.10 = data.10$models
prob.10 = data.10$prob
# aic.10 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_mnase_k562_10class_aic.txt")))
data.10 = NULL
ref = list(ref.10, ref.9, ref.8, ref.7, ref.6, ref.5, ref.4, ref.3, ref.2, ref.1)
prob = list(prob.10, prob.9, prob.8, prob.7, prob.6, prob.5, prob.4, prob.3, prob.2, prob.1)
# aic = c(aic.15, aic.14, aic.13, aic.12, aic.11, aic.10, aic.9, aic.8, aic.7, aic.6, aic.5, aic.4, aic.3, aic.2, aic.1)
aic = rep(0, length(ref))
# number of runs
n_run = length(ref)
# number of different classes overall
n_class_tot = sum(unlist(lapply(ref, nrow)))
# max value of K
n_class_max = max(unlist(lapply(ref, nrow)))
# some colors
colors = rep(brewer.pal(9, "Set1")[2], n_class_max)
# construct a matrix with all discovered references on the rows
references = matrix(nrow=n_class_tot, ncol=ncol(ref[[1]]))
run_value = vector(length=n_class_tot)
k_value = vector(length=n_class_tot)
probabilities = vector(length=n_class_tot)
k = 1
for(i in 1:n_run)
{
for(j in 1:nrow(ref[[i]]))
{ references[k,] = ref[[i]][j,]
probabilities[k] = prob[[i]][j]
run_value[k] = i
k_value[k] = j
k = k + 1
}
}
# distance matrix between all references
distances = distance.model(references)
rownames(distances) = 1:nrow(distances)
colnames(distances) = 1:ncol(distances)
plot.references(file.path("results","bulk_sequencing", "ctcf_mnase.png"),
references, probabilities, colors, aic, distances, n_run, run_value, n_class_max)
c4science · Help