The Figure S2 is the showcase for geomeTriD package to
present multiple genomic signals along with 3D models for Sox2
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_chr3_wid <- GRanges('chr3:30000000-40000000')
gi <- importGInteractionsFromUrl(urls=urls, resolution=10000, range=range_chr3_wid,
format='cool', normalization='balanced')
## create 3d model by mdsPlot function
range_chr3 <- GRanges('chr3:34500000-35500000')
model3d <- lapply(gi, mdsPlot, range = range_chr3, k=3, render = 'granges')## Warning in checkGI(gi, fixedBin = TRUE): There are NA values in the gi score.
## It will be removed.## initial value 12.913756
## iter 5 value 9.532312
## iter 10 value 6.952131
## iter 15 value 4.341098
## iter 20 value 3.464997
## iter 25 value 3.053412
## iter 30 value 2.788117
## iter 35 value 2.574587
## iter 40 value 2.359898
## iter 45 value 2.188246
## iter 50 value 2.113491
## final value 2.113491
## stopped after 50 iterations## Warning in checkGI(gi, fixedBin = TRUE): There are NA values in the gi score.
## It will be removed.## initial value 12.778648
## iter 5 value 9.009644
## iter 10 value 6.409643
## iter 15 value 4.584916
## iter 20 value 3.789602
## iter 25 value 3.296491
## iter 30 value 2.924350
## iter 35 value 2.660350
## iter 40 value 2.460614
## iter 45 value 2.306374
## iter 50 value 2.216844
## final value 2.216844
## stopped after 50 iterations
### 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_chr3)
#### 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])
scr <- GRanges("chr3:34810610-34817610:+")
scr$gene_id <- NA
scr$label <- 'SCR'
feature.gr <- c(feature.gr, scr)
#### 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('chr3:34695000-34840000')
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', 'chr3'), '.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.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.3 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.47.0
## [3] org.Mm.eg.db_3.22.0
## [4] TxDb.Mmusculus.UCSC.mm39.knownGene_3.22.0
## [5] GenomicFeatures_1.62.0
## [6] AnnotationDbi_1.72.0
## [7] Biobase_2.70.0
## [8] GenomicRanges_1.62.0
## [9] Seqinfo_1.0.0
## [10] IRanges_2.44.0
## [11] S4Vectors_0.48.0
## [12] BiocGenerics_0.56.0
## [13] generics_0.1.4
## [14] geomeTriD.documentation_0.0.6
## [15] geomeTriD_1.5.0
##
## loaded via a namespace (and not attached):
## [1] BiocIO_1.20.0 bitops_1.0-9
## [3] filelock_1.0.3 tibble_3.3.0
## [5] R.oo_1.27.1 XML_3.99-0.20
## [7] rpart_4.1.24 lifecycle_1.0.4
## [9] httr2_1.2.1 aricode_1.0.3
## [11] globals_0.18.0 lattice_0.22-7
## [13] ensembldb_2.34.0 MASS_7.3-65
## [15] backports_1.5.0 magrittr_2.0.4
## [17] Hmisc_5.2-4 sass_0.4.10
## [19] rmarkdown_2.30 jquerylib_0.1.4
## [21] yaml_2.3.10 plotrix_3.8-4
## [23] Gviz_1.54.0 DBI_1.2.3
## [25] abind_1.4-8 purrr_1.2.0
## [27] R.utils_2.13.0 AnnotationFilter_1.34.0
## [29] biovizBase_1.58.0 RCurl_1.98-1.17
## [31] rgl_1.3.24 nnet_7.3-20
## [33] VariantAnnotation_1.56.0 rappdirs_0.3.3
## [35] grImport_0.9-7 listenv_0.10.0
## [37] parallelly_1.45.1 pkgdown_2.2.0
## [39] codetools_0.2-20 DelayedArray_0.36.0
## [41] tidyselect_1.2.1 UCSC.utils_1.6.0
## [43] farver_2.1.2 matrixStats_1.5.0
## [45] BiocFileCache_3.0.0 base64enc_0.1-3
## [47] GenomicAlignments_1.46.0 jsonlite_2.0.0
## [49] progressr_0.18.0 Formula_1.2-5
## [51] systemfonts_1.3.1 dbscan_1.2.3
## [53] tools_4.5.2 progress_1.2.3
## [55] ragg_1.5.0 strawr_0.0.92
## [57] Rcpp_1.1.0 glue_1.8.0
## [59] gridExtra_2.3 SparseArray_1.10.1
## [61] xfun_0.54 MatrixGenerics_1.22.0
## [63] GenomeInfoDb_1.46.0 dplyr_1.1.4
## [65] withr_3.0.2 fastmap_1.2.0
## [67] latticeExtra_0.6-31 rhdf5filters_1.22.0
## [69] digest_0.6.38 R6_2.6.1
## [71] textshaping_1.0.4 colorspace_2.1-2
## [73] jpeg_0.1-11 dichromat_2.0-0.1
## [75] biomaRt_2.66.0 RSQLite_2.4.4
## [77] cigarillo_1.0.0 R.methodsS3_1.8.2
## [79] data.table_1.17.8 rtracklayer_1.70.0
## [81] prettyunits_1.2.0 InteractionSet_1.38.0
## [83] httr_1.4.7 htmlwidgets_1.6.4
## [85] S4Arrays_1.10.0 pkgconfig_2.0.3
## [87] gtable_0.3.6 blob_1.2.4
## [89] S7_0.2.0 XVector_0.50.0
## [91] htmltools_0.5.8.1 ProtGenerics_1.42.0
## [93] clue_0.3-66 scales_1.4.0
## [95] png_0.1-8 knitr_1.50
## [97] rstudioapi_0.17.1 rjson_0.2.23
## [99] checkmate_2.3.3 curl_7.0.0
## [101] cachem_1.1.0 rhdf5_2.54.0
## [103] stringr_1.6.0 parallel_4.5.2
## [105] foreign_0.8-90 restfulr_0.0.16
## [107] desc_1.4.3 pillar_1.11.1
## [109] vctrs_0.6.5 RANN_2.6.2
## [111] dbplyr_2.5.1 cluster_2.1.8.1
## [113] htmlTable_2.4.3 evaluate_1.0.5
## [115] cli_3.6.5 compiler_4.5.2
## [117] Rsamtools_2.26.0 rlang_1.1.6
## [119] crayon_1.5.3 future.apply_1.20.0
## [121] interp_1.1-6 fs_1.6.6
## [123] stringi_1.8.7 deldir_2.0-4
## [125] BiocParallel_1.44.0 txdbmaker_1.6.0
## [127] Biostrings_2.78.0 lazyeval_0.2.2
## [129] Matrix_1.7-4 BSgenome_1.78.0
## [131] hms_1.1.4 bit64_4.6.0-1
## [133] future_1.67.0 ggplot2_4.0.0
## [135] Rhdf5lib_1.32.0 KEGGREST_1.50.0
## [137] SummarizedExperiment_1.40.0 igraph_2.2.1
## [139] memoise_2.0.1 bslib_0.9.0
## [141] bit_4.6.0