Homec4science
Diffusion scATAC-seq (master) Browse

scATAC-seq/scripts/bulk_sequencing/analysis_cluster_ctcf_dnase_k562.R
faf6f1d50a68master

Last Change
Raw File

analysis_cluster_ctcf_dnase_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.references(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_1class_ref.mat"))
ref.1 = data.1$references
prob.1 = data.1$prob
aic.1 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_1class_aic.txt")))
data.1 = NULL
data.2 = read.references(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_2class_ref.mat"))
ref.2 = data.2$references
prob.2 = data.2$prob
aic.2 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_2class_aic.txt")))
data.2 = NULL
data.3 = read.references(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_3class_ref.mat"))
ref.3 = data.3$references
prob.3 = data.3$prob
aic.3 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_3class_aic.txt")))
data.3 = NULL
data.4 = read.references(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_4class_ref.mat"))
ref.4 = data.4$references
prob.4 = data.4$prob
aic.4 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_4class_aic.txt")))
data.4 = NULL
data.5 = read.references(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_5class_ref.mat"))
ref.5 = data.5$references
prob.5 = data.5$prob
aic.5 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_5class_aic.txt")))
data.5 = NULL
data.6 = read.references(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_6class_ref.mat"))
ref.6 = data.6$references
prob.6 = data.6$prob
aic.6 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_6class_aic.txt")))
data.6 = NULL
data.7 = read.references(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_7class_ref.mat"))
ref.7 = data.7$references
prob.7 = data.7$prob
aic.7 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_7class_aic.txt")))
data.7 = NULL
data.8 = read.references(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_8class_ref.mat"))
ref.8 = data.8$references
prob.8 = data.8$prob
aic.8 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_8class_aic.txt")))
data.8 = NULL
data.9 = read.references(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_9class_ref.mat"))
ref.9 = data.9$references
prob.9 = data.9$prob
aic.9 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_9class_aic.txt")))
data.9 = NULL
data.10 = read.references(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_10class_ref.mat"))
ref.10 = data.10$references
prob.10 = data.10$prob
aic.10 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_10class_aic.txt")))
data.10 = NULL
data.11 = read.references(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_11class_ref.mat"))
ref.11 = data.11$references
prob.11 = data.11$prob
aic.11 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_11class_aic.txt")))
data.11 = NULL
data.12 = read.references(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_12class_ref.mat"))
ref.12 = data.12$references
prob.12 = data.12$prob
aic.12 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_12class_aic.txt")))
data.12 = NULL
data.13 = read.references(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_13class_ref.mat"))
ref.13 = data.13$references
prob.13 = data.13$prob
aic.13 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_13class_aic.txt")))
data.13 = NULL
data.14 = read.references(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_14class_ref.mat"))
ref.14 = data.14$references
prob.14 = data.14$prob
aic.14 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_14class_aic.txt")))
data.14 = NULL
data.15 = read.references(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_15class_ref.mat"))
ref.15 = data.15$references
prob.15 = data.15$prob
aic.15 = as.matrix(read.table(file.path("results", "bulk_sequencing", "ctcf_dnase_k562_15class_aic.txt")))
data.15 = NULL
ref = list(ref.15, ref.14, ref.13, ref.12, ref.11, ref.10, ref.9, ref.8, ref.7, ref.6, ref.5, ref.4, ref.3, ref.2, ref.1)
prob = list(prob.15, prob.14, prob.13, prob.12, prob.11, 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)
# 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")[1], 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.ref(references)
rownames(distances) = 1:nrow(distances)
colnames(distances) = 1:ncol(distances)
plot.references(file.path("results","bulk_sequencing", "ctcf_dnase.png"),
references, probabilities, colors, aic, distances, n_run, run_value, n_class_max)
c4science · Help