MAGICAL visualzation
MAGICAL offers several ways of visualizing circuits, all of which
make use of the package karyoploteR.
Read the results
Let’s start by reading the results from the MAGICAL analysis, “MAGICAL_selected_regulatory_circuits.txt”. We suggest running this demo with our provided file, because as mentioned before, running MAGICAL can return very similar but slightly different results each time.
Basic visualization
MAGICAL offers 3 ways of visualizing the circuits: by circuit index, gene, or peak. For either of these ways, MAGICAL will plot a certain region of one chromosome which includes the circuit(s). Basically, the following information will be manifested on the plot:
- The chromosome id and coverage of the region
- The starting site of the gene
- The starting site and ending site of the peak
- The top 3 (or less if there are less than 3) TFs that are most likely to bind to the peak.
Also, see the advanced options section for the advanced visualization option of gene and peak track. For this section, we are not showing them.
visualize by circuit index
With function plot_circuits_with_idx, you can visualize
one circuit by its index.
visualize by gene
With function plot_circuits_with_gene, you can visualize
all circuits that contain a specific gene.
# visualize circuits with gene ACOT11
plot_circuits_with_gene(data, gene = "ACOT11", gene_track = F, peak_track = F)
# visualize circuits with gene ASAP3
plot_circuits_with_gene(data, gene = "ASAP3", gene_track = F, peak_track = F)
visualize by peak
With function plot_circuits_with_peak, you can visualize
all circuits that contain a specific peak.
# visualize circuits with peak chr4:184817550-184818412
plot_circuits_with_peak(data, peak_chr = "chr4", peak_start = 184817550, peak_end = 184818412, gene_track = F, peak_track = F)
Advanced options
All the 3 functions offer an advanced options, which allows you to see all the genes and peaks around the circuit region. This is useful if you want to see the specificity of the circuit(s) in the context of the whole region, which is, the other genes/peaks are not involved in the circuit(s).
Gene track
Gene track can be visualized by setting the parameter
gene_track to TRUE (default). After that, you
need another parameter, TxDb, to specify the reference
genome from which the function will query genes. Before that, you must
download the corresponding TXDb package from Bioconductor. For example,
if your reference genome is hg38, you need to download the package
TxDb.Hsapiens.UCSC.hg38.knownGene.
## Here our reference genome is hg38
# if (!require("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#
# BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene")
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
# visualize the first circuit with gene track
plot_circuits_with_idx(data, 1, gene_track = T, TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene, peak_track = F)
Peak track
By setting the parameter gene_track to TRUE
(default), you can see the density of peaks around the region. You also
need another input, peaks, which contains the chromosome,
start, and end of the peaks. File “scATAC peaks.txt” in our demo
data provides an example.
# read the peaks
scATAC_peak_file_path <- "Demo input files/scATAC peaks.txt"
peaks <- read.table(scATAC_peak_file_path, header = F, sep = "\t")
head(peaks)
#> V1 V2 V3 V4 V5 V6
#> 1 1 chr1 9752 10693 942 *
#> 2 2 chr1 15729 16627 899 *
#> 3 3 chr1 17059 17955 897 *
#> 4 4 chr1 180705 181746 1042 *
#> 5 5 chr1 182981 183730 750 *
#> 6 6 chr1 183736 184785 1050 *# visualize the first circuit with peak track
plot_circuits_with_idx(data, 1, gene_track = F, peak_track = T, peaks = peaks[, 2:4])
Now let’s put both gene and peak tracks together.
