The Figure 2 is the showcase for geomeTriD package to
present multiple genomic signals along with 3D models for a
microcompartment/loop.
library(geomeTriD)
library(geomeTriD.documentation)
library(GenomicRanges)
library(TxDb.Mmusculus.UCSC.mm39.knownGene)
library(org.Mm.eg.db)
library(trackViewer)
library(RColorBrewer)The Region Capture Micro-C data were downloaded from GSE207225. The genomic signals of ChIP-seq were downloaded from GSE178982 and remapped to mm39 genome.
## import Genomic interaction data
geo_acc <- c("DMSO"='GSM6281851',
'IAA'='GSM6281852')
url <- 'https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM6281nnn';
urls <- mapply(geo_acc, names(geo_acc), FUN=function(.f, .cond){
file.path(url, .f, 'suppl',
paste0(.f, '_RCMC_BR1_merged_allCap_', .cond,
'_mm39.merged.50.mcool'),
fsep='/')
})
## import the genomic interactions
range_chr8_wid <- GRanges('chr8:84000000-92000000')
gi <- importGInteractionsFromUrl(urls=urls, resolution=500, range=range_chr8_wid,
format='cool', normalization='balanced')## adding rname 'https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM6281nnn/GSM6281851/suppl/GSM6281851_RCMC_BR1_merged_allCap_DMSO_mm39.merged.50.mcool'## adding rname 'https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM6281nnn/GSM6281852/suppl/GSM6281852_RCMC_BR1_merged_allCap_IAA_mm39.merged.50.mcool'
## create 3d model by mdsPlot function
range_chr8 <- GRanges('chr8:85550000-85800000')
model3d <- lapply(gi, mdsPlot, range = range_chr8, k=3, render = 'granges')## Warning in checkGI(gi, fixedBin = TRUE): There are NA values in the gi score.
## It will be removed.## Some region are missing from the input gi.## initial value 17.000699
## iter 5 value 14.159797
## iter 10 value 13.013113
## iter 15 value 12.575326
## final value 12.445963
## converged## Warning in checkGI(gi, fixedBin = TRUE): There are NA values in the gi score.
## It will be removed.## Some region are missing from the input gi.## initial value 17.114933
## iter 5 value 13.719667
## iter 10 value 12.444197
## iter 15 value 12.102132
## final value 11.986303
## converged
### create gene annotations
#### get all genes
feature.gr <- genes(TxDb.Mmusculus.UCSC.mm39.knownGene)## 3 genes were dropped because they have exons located on both strands of the
## same reference sequence or on more than one reference sequence, so cannot be
## represented by a single genomic range.
## Use 'single.strand.genes.only=FALSE' to get all the genes in a GRangesList
## object, or use suppressMessages() to suppress this message.
#### subset the data by target viewer region
feature.gr <- subsetByOverlaps(feature.gr, range_chr8)
#### assign symbols for each gene
symbols <- mget(feature.gr$gene_id, org.Mm.egSYMBOL, ifnotfound = NA)
feature.gr$label[lengths(symbols) == 1] <-
unlist(symbols[lengths(symbols) == 1])
#### assign colors for each gene
feature.gr$col <- sample(1:7, length(feature.gr), replace = TRUE)
### plot for the detailed region ii to show the open or close loop.
subregion <- GRanges('chr8:85690000-85730000')
model3d.sub <- lapply(model3d, subsetByOverlaps, ranges=subregion)
model3d.sub[['IAA']] <- alignCoor(model3d.sub[['IAA']], model3d.sub[['DMSO']])
### import cohesion ctcf, and yy1 signals, realigned to mm39 for GSE178982
pf <- system.file('extdata', package = 'geomeTriD.documentation')
bws_files <- dir(file.path(pf, 'GSE178982'), '.bw', full.names = TRUE)
(names(bws_files) <- sub('^.*?_(IAA|UT)_(.*?).CPM.*$', '\\2', bws_files))## [1] "CTCF" "RAD21" "SMC1A" "SMC3" "YY1" "CTCF" "RAD21" "SMC1A" "SMC3"
## [10] "YY1"
bw_UT <- bws_files[grepl('_UT_', bws_files)]
bw_IAA <- bws_files[grepl('_IAA_', bws_files)]
signals_UT <- lapply(bw_UT, importScore, format='BigWig', ranges=subregion)
signals_IAA <- lapply(bw_IAA[names(bw_UT)], importScore, format='BigWig', ranges=subregion)
colorSets <- c(CTCF="cyan",YY1="yellow",
RAD21="red", SMC1A="green", SMC3="blue")
for(i in seq_along(signals_UT)){
setTrackStyleParam(signals_UT[[i]], "color", c("gray30", colorSets[names(signals_UT)[i]]))
setTrackStyleParam(signals_IAA[[i]], "color", c("gray30", colorSets[names(signals_UT)[i]]))
}
## set maximal lwd for UT and IAA samples according their max
genomicScoreRanges <- lapply(signals_UT, function(.ele) range(.ele$dat$score))geomeTriD
### create the structure
dmso <- view3dStructure(model3d.sub[['DMSO']], feature.gr = feature.gr,
genomicSigs=signals_UT,
reverseGenomicSigs = FALSE,
genomicScoreRange = genomicScoreRanges,
lwd.maxGenomicSigs = 20,
k = 3, renderer = 'none')## Feature type is missing. Set as gene.
iaa <- view3dStructure(model3d.sub[['IAA']], feature.gr = feature.gr,
genomicSigs=signals_IAA,
reverseGenomicSigs = FALSE,
genomicScoreRange = genomicScoreRanges,
lwd.maxGenomicSigs = 20,
k=3, renderer = 'none')## Feature type is missing. Set as gene.
iaa <- lapply(iaa, function(.ele){
.ele$side <- 'right'
.ele
})
threeJsViewer(dmso, iaa, title = c('DMSO control', 'IAA 3h'))
#widget <- threeJsViewer(dmso, iaa, title = c('DMSO control', 'IAA 3h'), background = 'white')
#tempfile <- 'Fig2.part1.html'#tempfile(fileext = '.html', pattern = 'RCMC_BR1_IAA_vs_DMSO_Klf1_II.3jsViewer.')
#htmlwidgets::saveWidget(widget, file=tempfile)
#utils::browseURL(tempfile)The following plot clearly highlights the top 10 interaction events for DMSO and IAA samples. These interactions shift from the long-distance regions of Hook2 to more localized regions.
## extract backbone coordinates, which will be used as the bone for RCMC data
xyz_dmso <- extractBackbonePositions(dmso)
xyz_iaa <- extractBackbonePositions(iaa)
DMSO <- gi[['DMSO']]
IAA <- gi[['IAA']]
DMSO.subset_II <- subsetByOverlaps(DMSO, subregion)
IAA.subset_II <- subsetByOverlaps(IAA, subregion)
DMSO.subset_II <- DMSO.subset_II[distance(first(DMSO.subset_II), second(DMSO.subset_II))>5000]
IAA.subset_II <- IAA.subset_II[distance(first(IAA.subset_II), second(IAA.subset_II))>5000]
hic_dmso_II <- create3dGenomicSignals(
DMSO.subset_II,
xyz_dmso,
name='dmsoII', # name prefix for the geometry
tag='dmsoII', # name for the layer in the scene
color = c('white', brewer.pal(9, 'YlOrRd')),
topN=10, # only plot the top 10 events ordered by the scores.
lwd.maxGenomicSigs = 3,
alpha=0.5
)
hic_iaa_II <- create3dGenomicSignals(
IAA.subset_II,
xyz_iaa,
name='iaaII', # name prefix for the geometry
tag='iaaII', # name for the layer in the scene
color = c('white', brewer.pal(9, 'YlOrRd')),
topN=10, # only plot the top 10 events ordered by the scores.
lwd.maxGenomicSigs = 3,
alpha=0.5
)
hic_iaa_II <- lapply(hic_iaa_II, function(.ele){
.ele$side <- 'right'
.ele
})
threeJsViewer(dmso, iaa, hic_dmso_II, hic_iaa_II, title = c('DMSO', 'IAA 3h'), background = c('gray10', 'gray20', 'gray20', 'gray10'))
#widget <- threeJsViewer(dmso, iaa, hic_dmso_II, hic_iaa_II, title = c('DMSO', 'IAA 3h'), background = c('gray10', 'gray20', 'gray20', 'gray10'))
#tempfile <- 'Fig2.part2.html'#tempfile(fileext = '.html', pattern = 'RCMC_BR1_IAA_vs_DMSO_Klf1_II.3jsViewer.')
#htmlwidgets::saveWidget(widget, file=tempfile)
#utils::browseURL(tempfile)## 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] RColorBrewer_1.1-3
## [2] trackViewer_1.45.0
## [3] org.Mm.eg.db_3.21.0
## [4] TxDb.Mmusculus.UCSC.mm39.knownGene_3.21.0
## [5] GenomicFeatures_1.61.3
## [6] AnnotationDbi_1.71.0
## [7] Biobase_2.69.0
## [8] GenomicRanges_1.61.0
## [9] GenomeInfoDb_1.45.4
## [10] IRanges_2.43.0
## [11] S4Vectors_0.47.0
## [12] BiocGenerics_0.55.0
## [13] generics_0.1.4
## [14] geomeTriD.documentation_0.0.3
## [15] geomeTriD_1.3.11
##
## loaded via a namespace (and not attached):
## [1] strawr_0.0.92 rstudioapi_0.17.1
## [3] jsonlite_2.0.0 magrittr_2.0.3
## [5] farver_2.1.2 rmarkdown_2.29
## [7] fs_1.6.6 BiocIO_1.19.0
## [9] ragg_1.4.0 vctrs_0.6.5
## [11] memoise_2.0.1 Rsamtools_2.25.0
## [13] RCurl_1.98-1.17 base64enc_0.1-3
## [15] htmltools_0.5.8.1 S4Arrays_1.9.1
## [17] progress_1.2.3 plotrix_3.8-4
## [19] curl_6.3.0 Rhdf5lib_1.31.0
## [21] rhdf5_2.53.1 SparseArray_1.9.0
## [23] Formula_1.2-5 sass_0.4.10
## [25] parallelly_1.45.0 bslib_0.9.0
## [27] htmlwidgets_1.6.4 desc_1.4.3
## [29] Gviz_1.53.0 httr2_1.1.2
## [31] cachem_1.1.0 GenomicAlignments_1.45.0
## [33] igraph_2.1.4 lifecycle_1.0.4
## [35] pkgconfig_2.0.3 Matrix_1.7-3
## [37] R6_2.6.1 fastmap_1.2.0
## [39] MatrixGenerics_1.21.0 future_1.58.0
## [41] aricode_1.0.3 clue_0.3-66
## [43] digest_0.6.37 colorspace_2.1-1
## [45] textshaping_1.0.1 Hmisc_5.2-3
## [47] RSQLite_2.4.1 filelock_1.0.3
## [49] progressr_0.15.1 httr_1.4.7
## [51] abind_1.4-8 compiler_4.5.0
## [53] withr_3.0.2 bit64_4.6.0-1
## [55] backports_1.5.0 htmlTable_2.4.3
## [57] BiocParallel_1.43.3 DBI_1.2.3
## [59] R.utils_2.13.0 biomaRt_2.65.0
## [61] MASS_7.3-65 rappdirs_0.3.3
## [63] DelayedArray_0.35.1 rjson_0.2.23
## [65] tools_4.5.0 foreign_0.8-90
## [67] future.apply_1.20.0 nnet_7.3-20
## [69] R.oo_1.27.1 glue_1.8.0
## [71] restfulr_0.0.15 dbscan_1.2.2
## [73] InteractionSet_1.37.0 rhdf5filters_1.21.0
## [75] checkmate_2.3.2 cluster_2.1.8.1
## [77] gtable_0.3.6 BSgenome_1.77.0
## [79] R.methodsS3_1.8.2 ensembldb_2.33.0
## [81] data.table_1.17.4 hms_1.1.3
## [83] xml2_1.3.8 XVector_0.49.0
## [85] RANN_2.6.2 pillar_1.10.2
## [87] stringr_1.5.1 dplyr_1.1.4
## [89] BiocFileCache_2.99.5 lattice_0.22-7
## [91] deldir_2.0-4 rtracklayer_1.69.0
## [93] bit_4.6.0 biovizBase_1.57.0
## [95] tidyselect_1.2.1 Biostrings_2.77.1
## [97] knitr_1.50 gridExtra_2.3
## [99] ProtGenerics_1.41.0 SummarizedExperiment_1.39.0
## [101] xfun_0.52 matrixStats_1.5.0
## [103] stringi_1.8.7 UCSC.utils_1.5.0
## [105] lazyeval_0.2.2 yaml_2.3.10
## [107] evaluate_1.0.3 codetools_0.2-20
## [109] interp_1.1-6 tibble_3.3.0
## [111] cli_3.6.5 rpart_4.1.24
## [113] systemfonts_1.2.3 jquerylib_0.1.4
## [115] dichromat_2.0-0.1 Rcpp_1.0.14
## [117] globals_0.18.0 grImport_0.9-7
## [119] dbplyr_2.5.0 png_0.1-8
## [121] XML_3.99-0.18 parallel_4.5.0
## [123] pkgdown_2.1.3 rgl_1.3.18
## [125] ggplot2_3.5.2 blob_1.2.4
## [127] prettyunits_1.2.0 jpeg_0.1-11
## [129] latticeExtra_0.6-30 AnnotationFilter_1.33.0
## [131] bitops_1.0-9 txdbmaker_1.5.4
## [133] listenv_0.9.1 VariantAnnotation_1.55.0
## [135] scales_1.4.0 purrr_1.0.4
## [137] crayon_1.5.3 rlang_1.1.6
## [139] KEGGREST_1.49.0