Figure 3
Figure 3 showcases the geomeTriD package, demonstrating
how it presents the distance changes between a cis-regulatory element
and its target gene.
Load Libraries
library(trackViewer)
library(geomeTriD)
library(InteractionSet)
library(org.Dr.eg.db)
library(GenomicFeatures)
library(GenomeInfoDb)
library(GenomicRanges)
library(ChIPpeakAnno)Prepare the annotation and genomic data
This data were downloaded from GEO with accession GSE231771.
## prepare the TxDb object,
## or load the refseq gene annotation by library(TxDb.Drerio.UCSC.danRer11.refGene)
txdb <- txdbmaker::makeTxDbFromEnsembl(release=109, organism='Danio rerio')
## set the silencer coordinates
silencer <- GRanges("chr5", IRanges(21795000, 21800000))
## set data range
range <- GRanges("chr5", IRanges(21100000, 22000000))
## set the plot range
fig3PlotRange <- GRanges('chr5', IRanges(21300000, 2.2e+07))
## prepare gene annotation
genes <- GenomicFeatures::genes(txdb)
GenomeInfoDb::seqlevelsStyle(genes) <- "UCSC"
genes <- subsetByOverlaps(genes, range, ignore.strand=TRUE)
id <- genes$gene_id
eid <- ChIPpeakAnno::xget(id, org.Dr.egENSEMBL2EG, output = "first")
symbol <- id
symbol[!is.na(eid)] <- ChIPpeakAnno::xget(eid[!is.na(eid)],
org.Dr.egSYMBOL,
output = "first")
names(symbol) <- id
feature.gr <- subsetByOverlaps(genes, range, ignore.strand=TRUE)
feature.gr$label <- symbol[feature.gr$gene_id]
# remove genes without symbols to clean the plot
feature.gr <- feature.gr[!grepl('ENSDARG', feature.gr$label)]
feature.gr$type <- 'gene'
feature.gr$col <- 1 # set black color to uninterested genes
feature.gr$col[feature.gr$label %in% c('tenm1', 'smarca1', 'sh2d1aa')] <-
c(2, 4, 5) ## set the colors for the given genes
silencer$label <- 'silencer'
silencer$type <- 'cRE' ## cis regulatory element
silencer$col <- 3 ## set the color to green
silencer$pch <- 7 ## set the point as a cross in a square
feature.gr <- c(feature.gr, silencer)Import ATAC-seq signals
## set the data folder, all data are available in the extdata folder of this package
extdata <- system.file('extdata', 'GSE231771', package='geomeTriD.documentation')
## list all ATAC-seq bigwig files
bws <- dir(extdata, "bigWig", full.names = TRUE)
## import the ATAC-seq signales
## the ATAC-seq signals are HiCAR R2 read coverage
bwsigs <- lapply(bws, trackViewer::importScore, format="BigWig", ranges=range)
## create names for the data for each sample
(names(bwsigs) <-sub('.bigWig', ' ATAC-seq', basename(bws)))## [1] "0dpa ATAC-seq" "4dpa ATAC-seq"
## get the hic files. Here we use the prepared GInteraction objects
hics <- dir(extdata, "10000.ginteractions.rds", full.names = TRUE)
## set the sample names
(names(hics) <- sub("MAPS_signals_(.*?).10000.ginteractions.rds", "\\1 heatmap",
basename(hics)))## [1] "0dpa heatmap" "4dpa heatmap"
gis <- lapply(hics, function(.ele){
## The RDS files were generated by trackViewer function `importGInteractions`
## from the HIC file. Because the MAPS signal are values normalized by
## fitting models, we set the matrixType to 'observed' and normalization to 'NONE'
## range.ext are ranges with extension to include all the intrachromosal interactions
# range.ext <- range
# start(range.ext) <- start(range.ext) - 10000000
# end(range.ext) <- end(range.ext) + 10000000
# importGInteractions(.ele, ranges = range.ext,
# format = "hic",
# resolution = 10000,
# out = "GInteractions",
# normalization = "NONE",
# matrixType = "observed")
.ele <- readRDS(.ele)
})Plot ShRec3D results
ShRec3D can handler very sparse matrix, here we use the MAPS normalized matrix as input contact matrix.
## list all ShRec3D output files
superRec <- dir(extdata, pattern="SuperRec.*.chr5.rds", full.names = TRUE)
names(superRec) <-
sub('SuperRec.(.*).chr5.rds', '\\1',basename(superRec))
superRec.gr <- lapply(superRec, readRDS)
superRec.gr <- lapply(superRec.gr,
function(.ele){
subsetByOverlaps(.ele, range, type = 'within')
})
#feature.gr.s <- feature.gr[feature.gr$label %in% c('tenm1', 'smarca1', 'silencer')]
#feature.gr.s$col <- c('cyan', 'blue', 'green')
## assign colors for the interactions
gis.sub <- lapply(gis, subsetByOverlaps, ranges = feature.gr[feature.gr$label=='silencer'], minoverlap=2500)
gis.sub <- lapply(gis.sub, function(.ele){ ## remove the short distance interactions
.ele[distance(first(.ele), second(.ele))>5000]
})
gis.sub.scores <- range(unlist(lapply(gis.sub, function(.ele) .ele$score)))
gis.sub.scores <- seq(gis.sub.scores[1], gis.sub.scores[2],
length.out = 8)
colors <- colorRampPalette(c('green', 'darkred'))(8)
gis.sub <- lapply(gis.sub, function(.ele){
.ele$color <- colors[findInterval(.ele$score, vec = gis.sub.scores,
all.inside = TRUE)]
.ele
})
obj0 <- view3dStructure(superRec.gr[['0dpa']],
k=3, feature.gr=feature.gr,
genomicSigs=gis.sub[['0dpa heatmap']],
renderer = 'none', topN=5)
obj4 <- view3dStructure(superRec.gr[['4dpa']],
k=3, feature.gr=feature.gr,
genomicSigs=gis.sub[['4dpa heatmap']],
renderer = 'none', topN=5)
obj4 <- lapply(obj4, function(.ele) {
.ele$side = 'right'
.ele
})
threeJsViewer(obj0, obj4, title = c('0dpa', '4dpa'))
#widget <- threeJsViewer(obj0, obj4, title = c('0dpa', '4dpa'))
#tempfile <- 'Fig3.html'
#htmlwidgets::saveWidget(widget, file=tempfile)
#utils::browseURL(tempfile)
## list all interaction BEDPE files
(interactions <- dir(extdata,
"10k.2.sig3Dinteractions.bedpe",
full.names = TRUE))## [1] "/tmp/Rtmpkea0Nw/temp_libpath10fa667a28b4/geomeTriD.documentation/extdata/GSE231771/0dpa.10k.2.sig3Dinteractions.bedpe"
## [2] "/tmp/Rtmpkea0Nw/temp_libpath10fa667a28b4/geomeTriD.documentation/extdata/GSE231771/4dpa.10k.2.sig3Dinteractions.bedpe"
## import the loop for `loopBouquetPlot`
loops_loopBouquet <- lapply(interactions, function(.ele){
## read bedpe as text file
.loop <- read.delim(.ele, header = FALSE)
## create GInteractions object
.loop <- with(.loop, InteractionSet::GInteractions(
## narrow the interaction region for better visualization effect
GRanges(V1, IRanges(V2+2501, V3-2500)),
GRanges(V4, IRanges(V5+2501, V6-2500)),
## set the score to log2 value to decrease the difference of signals
score = log2(V8)
))
## subset the data with given range
.loop <- subsetByOverlaps(.loop, range)
## remove the interchromosal interactions
.loop <- .loop[seqnames(first(.loop))==seqnames(range)[1] &
seqnames(second(.loop))==seqnames(range)[1]]
})
## create names for the data for each sample
(names(loops_loopBouquet) <- sub(".10k.*$", " interactions",
basename(interactions)))## [1] "0dpa interactions" "4dpa interactions"Plot the data by loopBouquetPlot for subfigure
c and d
set.seed(1) ## set seed to control the plot layout
for(i in seq_along(loops_loopBouquet)){
## pdf(paste0(names(loops_loopBouquet)[i], '.pdf')) # if using terminal
geomeTriD::loopBouquetPlot(loops_loopBouquet[[i]],
fig3PlotRange,
feature.gr,
genomicSigs=bwsigs[[i]],
lwd.backbone=2,
col.backbone='gray10',
col.backbone_background = 'gray70',
lwd.gene = 2,
show_coor = TRUE,
coor_tick_unit = 1e4,
coor_mark_interval=2e5,
length.arrow=unit(0.1, 'inch'),
show_edges = FALSE)
grid::grid.text(label=names(loops_loopBouquet)[i],
x=.15, y=.9)
## dev.off() # if using terminal
}
Subfigure e
There are several steps to repeat the heatmaps and links plot:
Prepare the annotations
Import the ATAC-seq signals
Import the called interactions for link tracks
Import the interactions for heatmap tracks
Plot the data by
viewTracksfunction.
Import data for subfigure e
## create gene track object
trs <- geneTrack(id, txdb, symbol, asList = FALSE)
seqlevelsStyle(trs$dat) <- "UCSC" ## make sure the seqnames are in same style
## import the interactions as 'track' objects
loops <- lapply(interactions, function(.ele){
.loop <- read.delim(.ele, header = FALSE)
.loop <- with(.loop, GInteractions(
GRanges(V1, IRanges(V2+2501, V3-2500)),
GRanges(V4, IRanges(V5+2501, V6-2500)),
score = log2(V8)
))
.loop <- gi2track(.loop)
## set the tracktype to 'link' to plot as links
setTrackStyleParam(.loop, "tracktype", "link")
.loop
})
names(loops) <- sub(".10k.*$", " interactions", basename(interactions))
## check the order of the samples
stopifnot(identical(sub(" .*$", "", names(gis)),
sub(" .*$", "", names(loops))))
## create heatmap signals
heatmaps <- mapply(gis, interactions, FUN=function(.gi, .ele){
.loop <- read.delim(.ele, header = FALSE)
.loop <- with(.loop, GInteractions(
GRanges(V1, IRanges(V2+2501, V3-2500)),
GRanges(V4, IRanges(V5+2501, V6-2500)),
score = log2(V8)
))
## highlight the called interactions by add color the border
.gi$border_color <- NA
ol <- findOverlaps(.gi, .loop)
.gi$border_color[queryHits(ol)] <- "black"
## create 'track' object from GInteraction object
gi2track(.gi)
})Plot the subfigure e by viewTracks
## create the plot style by default algorithm
op <- optimizeStyle(
trackList = trackList(trs,
rev(bwsigs), ## reverse the signal
rev(loops), ## to keep 0 dpa samples
rev(heatmaps), ## in top
## set the height for each track
heightDist = c(2, length(bwsigs),
length(loops),
length(heatmaps)*8)))
tL <- op$tracks
sty <- op$style
for(i in which(grepl('ATAC-seq', names(tL)))){#R2 ATAC-seq signal
## set the ylim. The raw range is [0,1]
setTrackStyleParam(tL[[i]], "ylim", c(0.15, .8))
## set the y-axis to the right side
setTrackYaxisParam(tL[[i]], "main", FALSE)
## plot the y-axis lables
setTrackYaxisParam(tL[[i]], "label", TRUE)
## draw the y-axis
setTrackYaxisParam(tL[[i]], "draw", TRUE)
}
## set color keys for the link tracks
for(i in which(grepl('interactions', names(tL)))){
setTrackStyleParam(tL[[i]], "breaks",
## set the breaks from min to max
c(seq(from=2, to=6.5, by=.5), 10))
setTrackStyleParam(tL[[i]], "color",
c("lightblue", "yellow", "red"))
setTrackYaxisParam(tL[[i]], "draw", TRUE)
}
## set color keys for the heatmaps
for(i in which(grepl('heatmap', names(tL)))){
setTrackStyleParam(tL[[i]], "breaks",
c(seq(from=0, to=5, by=.5), 8.5))
setTrackStyleParam(tL[[i]], "color",
c("#FDF8F5", "#EA7556", "#4D0E0B"))
}
##
## pdf('Fig1a.pdf') # if using terminal
trackViewer::viewTracks(tL, viewerStyle = sty,
gr=fig3PlotRange, autoOptimizeStyle = FALSE)
## set the TADs to plot with the heatmap
tad <- GRanges("chr5",
IRanges(c(21060000, 21220000, 21326000, 21890000, 22000000),
c(21130000, 21290000, 21860000, 21960000, 22090000)))
## add the annotations to the heatmap
for(i in which(grepl('heatmap', names(tL)))){
## add TADs as dashed line (lty=2)
trackViewer::addInteractionAnnotation(tad, i, grid.lines,
gp=gpar(col="black",
lwd=1, lty=2))
## add the guide line for the silencer with slope -1 with lwd=3
trackViewer::addInteractionAnnotation(-21795000, i,
gp=gpar(col = "#009E73",
lwd=3, lty=3))
}
## dev.off() # if using terminalSessionInfo
## R version 4.5.0 (2025-04-11)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] grid stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ChIPpeakAnno_3.43.0 GenomicFeatures_1.61.3
## [3] org.Dr.eg.db_3.21.0 AnnotationDbi_1.71.0
## [5] InteractionSet_1.37.0 SummarizedExperiment_1.39.0
## [7] Biobase_2.69.0 MatrixGenerics_1.21.0
## [9] matrixStats_1.5.0 geomeTriD_1.3.11
## [11] trackViewer_1.45.0 GenomicRanges_1.61.0
## [13] GenomeInfoDb_1.45.4 IRanges_2.43.0
## [15] S4Vectors_0.47.0 BiocGenerics_0.55.0
## [17] generics_0.1.4
##
## loaded via a namespace (and not attached):
## [1] splines_4.5.0 BiocIO_1.19.0 bitops_1.0-9
## [4] filelock_1.0.3 tibble_3.3.0 graph_1.87.0
## [7] XML_3.99-0.18 rpart_4.1.24 lifecycle_1.0.4
## [10] httr2_1.1.2 pwalign_1.5.0 aricode_1.0.3
## [13] globals_0.18.0 lattice_0.22-7 ensembldb_2.33.0
## [16] MASS_7.3-65 backports_1.5.0 magrittr_2.0.3
## [19] Hmisc_5.2-3 sass_0.4.10 rmarkdown_2.29
## [22] jquerylib_0.1.4 yaml_2.3.10 plotrix_3.8-4
## [25] Gviz_1.53.0 DBI_1.2.3 RColorBrewer_1.1-3
## [28] abind_1.4-8 purrr_1.0.4 AnnotationFilter_1.33.0
## [31] biovizBase_1.57.0 RCurl_1.98-1.17 rgl_1.3.18
## [34] nnet_7.3-20 VariantAnnotation_1.55.0 rappdirs_0.3.3
## [37] grImport_0.9-7 listenv_0.9.1 parallelly_1.45.0
## [40] pkgdown_2.1.3 codetools_0.2-20 DelayedArray_0.35.1
## [43] xml2_1.3.8 tidyselect_1.2.1 futile.logger_1.4.3
## [46] UCSC.utils_1.5.0 farver_2.1.2 universalmotif_1.27.2
## [49] BiocFileCache_2.99.5 base64enc_0.1-3 GenomicAlignments_1.45.0
## [52] jsonlite_2.0.0 multtest_2.65.0 progressr_0.15.1
## [55] Formula_1.2-5 survival_3.8-3 systemfonts_1.2.3
## [58] dbscan_1.2.2 tools_4.5.0 progress_1.2.3
## [61] ragg_1.4.0 strawr_0.0.92 Rcpp_1.0.14
## [64] glue_1.8.0 gridExtra_2.3 SparseArray_1.9.0
## [67] xfun_0.52 dplyr_1.1.4 formatR_1.14
## [70] fastmap_1.2.0 latticeExtra_0.6-30 rhdf5filters_1.21.0
## [73] digest_0.6.37 R6_2.6.1 textshaping_1.0.1
## [76] colorspace_2.1-1 jpeg_0.1-11 dichromat_2.0-0.1
## [79] biomaRt_2.65.0 RSQLite_2.4.1 tidyr_1.3.1
## [82] data.table_1.17.4 rtracklayer_1.69.0 prettyunits_1.2.0
## [85] httr_1.4.7 htmlwidgets_1.6.4 S4Arrays_1.9.1
## [88] regioneR_1.41.2 pkgconfig_2.0.3 gtable_0.3.6
## [91] blob_1.2.4 XVector_0.49.0 htmltools_0.5.8.1
## [94] RBGL_1.85.0 ProtGenerics_1.41.0 clue_0.3-66
## [97] scales_1.4.0 png_0.1-8 knitr_1.50
## [100] lambda.r_1.2.4 rstudioapi_0.17.1 rjson_0.2.23
## [103] checkmate_2.3.2 curl_6.3.0 cachem_1.1.0
## [106] rhdf5_2.53.1 stringr_1.5.1 parallel_4.5.0
## [109] foreign_0.8-90 restfulr_0.0.15 desc_1.4.3
## [112] pillar_1.10.2 vctrs_0.6.5 RANN_2.6.2
## [115] dbplyr_2.5.0 cluster_2.1.8.1 htmlTable_2.4.3
## [118] evaluate_1.0.3 VennDiagram_1.7.3 cli_3.6.5
## [121] compiler_4.5.0 futile.options_1.0.1 Rsamtools_2.25.0
## [124] rlang_1.1.6 crayon_1.5.3 future.apply_1.20.0
## [127] interp_1.1-6 fs_1.6.6 stringi_1.8.7
## [130] deldir_2.0-4 BiocParallel_1.43.3 txdbmaker_1.5.4
## [133] Biostrings_2.77.1 lazyeval_0.2.2 Matrix_1.7-3
## [136] BSgenome_1.77.0 hms_1.1.3 bit64_4.6.0-1
## [139] future_1.58.0 ggplot2_3.5.2 Rhdf5lib_1.31.0
## [142] KEGGREST_1.49.0 igraph_2.1.4 memoise_2.0.1
## [145] bslib_0.9.0 bit_4.6.0