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')
## 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 Under development (unstable) (2025-10-13 r88918)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sequoia 15.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.6-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## 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.2
## [3] org.Mm.eg.db_3.22.0
## [4] TxDb.Mmusculus.UCSC.mm39.knownGene_3.21.0
## [5] GenomicFeatures_1.61.6
## [6] AnnotationDbi_1.71.2
## [7] Biobase_2.69.1
## [8] GenomicRanges_1.61.5
## [9] Seqinfo_0.99.2
## [10] IRanges_2.43.5
## [11] S4Vectors_0.47.4
## [12] BiocGenerics_0.55.4
## [13] generics_0.1.4
## [14] geomeTriD.documentation_0.0.6
## [15] geomeTriD_1.3.19
##
## loaded via a namespace (and not attached):
## [1] strawr_0.0.92 rstudioapi_0.17.1
## [3] jsonlite_2.0.0 magrittr_2.0.4
## [5] farver_2.1.2 rmarkdown_2.30
## [7] fs_1.6.6 BiocIO_1.19.0
## [9] ragg_1.5.0 vctrs_0.6.5
## [11] memoise_2.0.1 Rsamtools_2.25.3
## [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_7.0.0 Rhdf5lib_1.31.1
## [21] rhdf5_2.53.6 SparseArray_1.9.1
## [23] Formula_1.2-5 sass_0.4.10
## [25] parallelly_1.45.1 bslib_0.9.0
## [27] htmlwidgets_1.6.4 desc_1.4.3
## [29] Gviz_1.53.1 httr2_1.2.1
## [31] cachem_1.1.0 GenomicAlignments_1.45.4
## [33] igraph_2.2.0 lifecycle_1.0.4
## [35] pkgconfig_2.0.3 Matrix_1.7-4
## [37] R6_2.6.1 fastmap_1.2.0
## [39] MatrixGenerics_1.21.0 future_1.67.0
## [41] aricode_1.0.3 clue_0.3-66
## [43] digest_0.6.37 colorspace_2.1-2
## [45] textshaping_1.0.4 Hmisc_5.2-4
## [47] RSQLite_2.4.3 filelock_1.0.3
## [49] progressr_0.17.0 httr_1.4.7
## [51] abind_1.4-8 compiler_4.6.0
## [53] withr_3.0.2 bit64_4.6.0-1
## [55] backports_1.5.0 htmlTable_2.4.3
## [57] S7_0.2.0 BiocParallel_1.43.4
## [59] DBI_1.2.3 R.utils_2.13.0
## [61] biomaRt_2.65.16 MASS_7.3-65
## [63] rappdirs_0.3.3 DelayedArray_0.35.3
## [65] rjson_0.2.23 tools_4.6.0
## [67] foreign_0.8-90 future.apply_1.20.0
## [69] nnet_7.3-20 R.oo_1.27.1
## [71] glue_1.8.0 restfulr_0.0.16
## [73] dbscan_1.2.3 InteractionSet_1.37.1
## [75] rhdf5filters_1.21.4 checkmate_2.3.3
## [77] cluster_2.1.8.1 gtable_0.3.6
## [79] BSgenome_1.77.2 R.methodsS3_1.8.2
## [81] ensembldb_2.33.2 data.table_1.17.8
## [83] hms_1.1.4 XVector_0.49.1
## [85] RANN_2.6.2 pillar_1.11.1
## [87] stringr_1.5.2 dplyr_1.1.4
## [89] BiocFileCache_2.99.6 lattice_0.22-7
## [91] deldir_2.0-4 rtracklayer_1.69.1
## [93] bit_4.6.0 biovizBase_1.57.1
## [95] tidyselect_1.2.1 Biostrings_2.77.2
## [97] knitr_1.50 gridExtra_2.3
## [99] ProtGenerics_1.41.0 SummarizedExperiment_1.39.2
## [101] xfun_0.53 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.5 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.3.1 jquerylib_0.1.4
## [115] dichromat_2.0-0.1 Rcpp_1.1.0
## [117] GenomeInfoDb_1.45.12 globals_0.18.0
## [119] grImport_0.9-7 dbplyr_2.5.1
## [121] png_0.1-8 XML_3.99-0.19
## [123] parallel_4.6.0 pkgdown_2.1.3
## [125] rgl_1.3.24 ggplot2_4.0.0
## [127] blob_1.2.4 prettyunits_1.2.0
## [129] jpeg_0.1-11 latticeExtra_0.6-31
## [131] AnnotationFilter_1.33.0 bitops_1.0-9
## [133] txdbmaker_1.5.6 listenv_0.9.1
## [135] VariantAnnotation_1.55.2 scales_1.4.0
## [137] purrr_1.1.0 crayon_1.5.3
## [139] rlang_1.1.6 KEGGREST_1.49.2