Homec4science
Diffusion scATAC-seq (master) Browse

scATAC-seq/scripts/test_dendrogram.R
1bfc0a93ba65master

Last Change
Raw File

test_dendrogram.R
View Options

setwd(file.path("/", "local", "groux", "scATAC-seq"))
# libraries
library(RColorBrewer)
library(motifStack)
library(TFBSTools)
library(MotifDb)
# functions
source(file.path("scripts", "functions.R"))
get_pfm_list = function(motifs, prefix_name)
{ pfm_list = list()
for(i in 1:dim(motifs)[3])
{ pfm_list[[i]] = new("pfm",
mat=motifs[,,i],
name=sprintf("%s class %d", prefix_name, i))
}
return(pfm_list)
}
# number of classes searched in the data
n_classes = c(17, 20, 30)
# load motifs from JASPAR clustering used to initialise the classes
motifs_jaspar_paths = c("data/pwm/jaspar_2018_clustering/cluster_3_node_23_20_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_4_node_31_3_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_4_node_22_4_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_7_node_13_2_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_5_node_20_5_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_31_node_4_5_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_31_node_5_2_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_1_node_74_2_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_66_1_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_33_node_1_2_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_20_4_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_7_node_17_5_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_32_node_STAT6_1_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_3_node_24_4_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_4_node_30_2_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_39_node_1_2_motifs_prob.mat",
"data/pwm/jaspar_2018_clustering/cluster_48_node_ctcf_1_motifs_prob.mat")
motifs_jaspar = lapply(motifs_jaspar_paths, read.table)
motifs_jaspar = lapply(motifs_jaspar, as.matrix)
motifs_jaspar_names = c("jun",
"HIF1A",
"myc",
"PU.1",
"CEBPb",
"Irf4",
"Irf2",
"LHX3",
"Fox1H",
"Sox3",
"Mef2c",
"Elf5",
"STAT6",
"NFE2",
"AHR",
"E2F2",
"CTCF")
tmp = list()
for(i in 1:length(motifs_jaspar))
{ rownames(motifs_jaspar[[i]]) = c("A", "C", "G", "T")
tmp[[i]] = new("pfm",
mat=motifs_jaspar[[i]],
name=motifs_jaspar_names[i])
}
motifs_jaspar = tmp
rm(tmp)
for(n_class in n_classes)
{ # load classes found
motifs_found = get_pfm_list(read.sequence.models(file.path("results",
"10xgenomics_PBMC_5k_peaks_classification_4",
sprintf("peaks_rmsk_sampled_sequences_1kb_%dclass_model.mat", n_class)))$models,
"")
# colors
red = brewer.pal(3, "Set1")[1]
blue = brewer.pal(3, "Set1")[2]
color = c(rep(blue, length(motifs_jaspar)),
rep(red, length(motifs_found)))
# plot logo stack with radial style
# X11(height=12, width=12)
png(filename=file.path(sprintf("test_%dclass.png", n_class)),
units="in", res=720, width=14, height=14)
motifStack(c(motifs_jaspar,
motifs_found),
layout="radialPhylog",
circle=0.3,
cleaves = 0.2,
clabel.leaves = 0.5,
col.bg=color,
col.bg.alpha=0.3,
col.leaves=color,
col.inner.label.circle=color,
inner.label.circle.width=0.05,
col.outer.label.circle=color,
outer.label.circle.width=0.02,
circle.motif=.5,
angle=350)
dev.off()
}
c4science · Help