Figure 4

Figure 4 showcases the geomeTriD package, demonstrating how it presents 3D models generated from Dip-C.

Load Libraries

library(geomeTriD)
library(geomeTriD.documentation)
library(GenomicRanges)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(org.Hs.eg.db)
library(clusterCrit)
library(colorRamps)

present single cell 3D structure for human data

Load and plot data

This data were downloaded from GEO with accession GSE117874.

## set the data folder, all data are available in the extdata folder of this package
extdata <- system.file('extdata', 'GSE117874', package='geomeTriD.documentation')
hickit_3dg <- dir(extdata, '3dg', full.names = TRUE)
hickit <- import3dg(hickit_3dg, parental_postfix=c('a', 'b'))[[1]]
## split it into maternal and paternal 3D structures.
hickit.a <- hickit[hickit$parental=='a'] ## mat
hickit.b <- hickit[hickit$parental=='b'] ## pat
## delete the parental information
hickit.a$parental <- hickit.b$parental <- NULL
## set data range
range <- GRanges('X:1-155270560')
## prepare the coordinates for the four superloops. 
supperloops <- GRanges(c('X:56800000-56850000',
                      'X:75350000-75400000',
                      'X:115000000-115050000',
                      'X:130850000-130900000'))
names(supperloops) <- c('ICCE', 'x75', 'DXZ4', 'FIRRE')
supperloops$label <- names(supperloops)
supperloops$col <- 2:5 ## set colors for each element
supperloops$type <- 'gene' ## set it as gene
## subset the data by given range
hickit.a <- subsetByOverlaps(hickit.a, range)
hickit.b <- subsetByOverlaps(hickit.b, range)
hickit.a <- alignCoor(hickit.a, hickit.b)
## prepare the backbone colors
resolution <- 3
backbone_colors <- matlab.like2(n=resolution*length(hickit.a))
## add the superloops as segments
c1 <- view3dStructure(hickit.a,
                      feature.gr=supperloops,
                      renderer = 'none',
                      region = range,
                      resolution=resolution,
                      show_coor=FALSE,
                      lwd.backbone = 0.25,
                      col.backbone = backbone_colors,
                      lwd.gene=6)
c2 <- view3dStructure(hickit.b,
                      feature.gr=supperloops,
                      renderer = 'none',
                      region = range,
                      resolution=resolution,
                      show_coor=FALSE,
                      lwd.backbone = 0.25,
                      col.backbone = backbone_colors,
                      lwd.gene=6)
c2 <- lapply(c2, function(.ele) { ## put pat to right pannel
  .ele$side = 'right'
  .ele
})
## view the data
threeJsViewer(c1, c2, title = c('GM12878 cell 3 mat', 'GM12878 cell 3 pat'))
#widget <-threeJsViewer(c1, c2, title = c('mat', 'pat'))
#tempfile <- 'Fig4.html'
#htmlwidgets::saveWidget(widget, file=tempfile)
#utils::browseURL(tempfile)
## view the superloops as spheres
addFeaturesAsSphere <- function(obj, features, ...){
  backbone <- extractBackbonePositions(obj)
  spheres <- createTADGeometries(features, backbone, ...)
  c(obj, spheres)
}
mat <- view3dStructure(hickit.a, renderer = 'none', region = range,
                       resolution=resolution, show_coor=FALSE,
                       lwd.backbone = 0.25,
                       col.backbone = backbone_colors)
mat <- addFeaturesAsSphere(mat, supperloops, alpha = 0.5)
pat <- view3dStructure(hickit.b, renderer = 'none', region = range,
                       resolution=resolution, show_coor=FALSE,
                       lwd.backbone = 0.25,
                       col.backbone = backbone_colors)
mat <- addFeaturesAsSphere(mat, supperloops, alpha = 0.5)
pat <- addFeaturesAsSphere(pat, supperloops, alpha = 0.5)
showPairs(mat, pat, title = c('GM12878 cell 3 mat', 'GM12878 cell 3 pat'), height = NULL)

Plot all cells to show the cell heterogeneity

## load the processed data.
hickit.a <- readRDS(file.path(extdata, 'hickit.a.rds'))
hickit.b <- readRDS(file.path(extdata, 'hickit.b.rds'))
widgets <- mapply(function(a, b, i){
  mat <- view3dStructure(a, renderer = 'none', region = range,
                         resolution=resolution, show_coor=FALSE, lwd.backbone = 0.25,
                      col.backbone = backbone_colors)
  mat <- addFeaturesAsSphere(mat, supperloops, alpha = 0.5)
  pat <- view3dStructure(b, renderer = 'none', region = range,
                         resolution=resolution, show_coor=FALSE, lwd.backbone = 0.25,
                      col.backbone = backbone_colors)
  pat <- addFeaturesAsSphere(pat, supperloops, alpha = 0.5)
  showPairs(mat, pat, title = paste('cell', i, c('mat', 'pat')),
            background = c("#11111188",
                           "#222222DD",
                           "#222222DD",
                           "#11111188"),
            height = NULL)
}, hickit.a, hickit.b, names(hickit.a), SIMPLIFY = FALSE)

widgets[[1]]
widgets[[2]]

present single cell 3D structure for mouse data

Figure 4 showcases the geomeTriD package, demonstrating how it presents 3D models generated from HiRES.

Load Libraries

library(geomeTriD)
library(geomeTriD.documentation)
library(GenomicRanges)
library(colorRamps)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(org.Mm.eg.db)
library(geometry)
## HiRES
extdata <- system.file('extdata', 'GSE223917', package = 'geomeTriD.documentation')
HiRES <- readRDS(file.path(extdata, 'HiRES.radial_glias.G1.chrX.rds'))
exprs <- readRDS(file.path(extdata, 'expr.radial_glias.G1.chrX.rds'))
pairs <- readRDS(file.path(extdata, 'sel.imput.pairs.chrX.rds'))

### supperloop
supperloops <- GRanges(c('chrX:50555744-50635321',
                      'chrX:75725458-75764699', # 4933407K13RiK, NR_029443
                      'chrX:103422010-103484957',
                      'chrX:105040854-105117090')) # 5530601H04RiK, NR_015467 and Pbdc1
names(supperloops) <- c('Firre', 'Dxz4', 'Xist/Tsix', 'x75')
supperloops$label <- names(supperloops)
supperloops$col <- 2:5 ## set colors for each element
supperloops$type <- 'gene' ## set it as gene

## plot region
range <- as(seqinfo(TxDb.Mmusculus.UCSC.mm10.knownGene)['chrX'], 'GRanges')

## annotations
genes <- genes(TxDb.Mmusculus.UCSC.mm10.knownGene)
##   66 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.
genes_symbols <- mget(genes$gene_id, org.Mm.egSYMBOL, ifnotfound = NA)
genes$symbols <- sapply(genes_symbols, `[`, i=1)
geneX <- genes[seqnames(genes)=='chrX']
geneX$label <- geneX$symbols
geneX <- geneX[geneX$symbols %in% rownames(exprs)]

## get data for female, must have mat and pat info for chromosome X.
HiRES.mat <- lapply(HiRES, function(.ele) {
  .ele <- .ele[.ele$parental=='(mat)']
  .ele$parental <- NULL
  .ele
})
l.mat <- lengths(HiRES.mat)
HiRES.pat <- lapply(HiRES, function(.ele) {
  .ele <- .ele[.ele$parental=='(pat)']
  .ele$parental <- NULL
  .ele
})
l.pat <- lengths(HiRES.pat)
k <- l.mat>0 & l.pat>0
HiRES.mat <- HiRES.mat[k]
HiRES.pat <- HiRES.pat[k]
# Make all GRanges object in a list same length by filling with NA
HiRES <- paddingGRangesList(c(HiRES.mat, HiRES.pat))
HiRES.mat.xyzs <- HiRES$xyzs[seq_along(HiRES.mat)]
HiRES.pat.xyzs <- HiRES$xyzs[-seq_along(HiRES.mat)]
## 
HiRES.mat <- lapply(HiRES.mat.xyzs, function(.ele){
  gr <- HiRES$gr
  mcols(gr) <- .ele[, c('x', 'y', 'z')]
  gr
})
HiRES.pat <- lapply(HiRES.pat.xyzs, function(.ele){
  gr <- HiRES$gr
  mcols(gr) <- .ele[, c('x', 'y', 'z')]
  gr
})

## backbone color
resolution <- 3
backbone_colors <- matlab.like2(n=resolution*length(HiRES.mat[[1]]))
backbone_bws <- backbone_colors
backbone_bws[1001:(length(backbone_colors)-1000)] <- 'gray'

getV <- function(points){
  vol <- convhulln(points, options='Fa')$vol
}

## calculate the Root Mean Square Deviation (RMSD)
RMSD_mat_pat <- mapply(function(mat, pat){
  mat <- as.data.frame(mcols(mat))
  pat <- as.data.frame(mcols(pat))
  mat <- fill_NA(mat)
  pat <- fill_NA(pat)
  pat <- alignCoor(pat, mat)
  mat.center <- colMeans(mat, na.rm = TRUE)
  pat.center <- colMeans(pat, na.rm = TRUE)
  mat <- t(t(mat) - mat.center)
  pat <- t(t(pat) - pat.center)
  mean(sqrt(rowSums((mat - pat)^2)), na.rm = TRUE)
}, HiRES.mat, HiRES.pat)
XlinkExpr <- data.frame(Xist=exprs['Xist', ], total=colSums(exprs), RMSD=RMSD_mat_pat)
XlinkExpr <- XlinkExpr[order(XlinkExpr$total), ]

fit <- lm(RMSD ~ total, data = XlinkExpr)
plot(XlinkExpr$total, XlinkExpr$RMSD,
     xlab='Total expression level of X-lined gene', 
     ylab='RMSD between maternal and paternal')
abline(fit, col = "blue", lwd=2)

widgets <- lapply(c(head(rownames(XlinkExpr), n=2),
                    tail(rownames(XlinkExpr), n=2)), function(cell_id){
    exprSig <- geneX
    exprSig$score <- exprs[geneX$symbols, cell_id]
    mat_cell <- HiRES.mat[[cell_id]]
    mcols(mat_cell) <- fill_NA(as.data.frame(mcols(mat_cell)))
    pat_cell <- HiRES.pat[[cell_id]]
    mcols(pat_cell) <- fill_NA(as.data.frame(mcols(pat_cell)))
    ## check the volumn of the chrX, 
    ## bigger one is Xa
    ## condensed one is Xi
    v_mat <- getV(as.matrix(mcols(mat_cell)))
    v_pat <- getV(as.matrix(mcols(pat_cell)))
    mat_only_pairs <- pairs$mat[[cell_id]]
    mat_only_pairs$color <- 'black'
    mat_only_pairs$lwd <- 4
    mat_cell <- view3dStructure(mat_cell,
                                feature.gr=supperloops,
                                lwd.gene = 4,
                                renderer = 'none',
                                region = range,
                                resolution=resolution,
                                genomicSigs = if(v_mat>v_pat) {
                                  list(mat_rna_reads=exprSig, mat_pairs=mat_only_pairs)
                                } else {list(mat_pairs=mat_only_pairs)},
                                signalTransformFun = c,
                                reverseGenomicSigs = FALSE,
                                show_coor=FALSE,
                                lwd.backbone = 0.25,
                                col.backbone = if(v_mat>v_pat) backbone_bws else backbone_colors)
    pat_only_pairs <- pairs$pat[[cell_id]]
    pat_only_pairs$color <- 'black'
    pat_only_pairs$lwd <- 4
    pat_cell <- view3dStructure(pat_cell,
                                feature.gr=supperloops,
                                lwd.gene = 4,
                                renderer = 'none',
                                region = range,
                                resolution=resolution,
                                genomicSigs = if(v_mat<=v_pat) {
                                  list(pat_rna_reads=exprSig, pat_pairs=pat_only_pairs)
                                } else {list(pat_pairs=pat_only_pairs)},
                                signalTransformFun = c,
                                reverseGenomicSigs = FALSE,
                                show_coor=FALSE,
                                lwd.backbone = 0.25,
                                col.backbone = if(v_mat>v_pat) backbone_colors else backbone_bws)
    # widget <-showPairs(mat_cell, pat_cell, title = paste(c('mat', 'pat'), cell_id))
    # tempfile <- paste0('cell.', cell_id, '.html')
    # htmlwidgets::saveWidget(widget, file=tempfile, selfcontained = FALSE, libdir = 'js')
    showPairs(mat_cell, pat_cell,
              title = paste(c('mat', 'pat'), cell_id),
              height=NULL)
})
## low expression of X-linked genes
widgets[[1]]
widgets[[2]]
## high expression of X-linked genes
widgets[[3]]
widgets[[4]]

SessionInfo

## 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] geometry_0.5.2                           
##  [2] org.Mm.eg.db_3.21.0                      
##  [3] TxDb.Mmusculus.UCSC.mm10.knownGene_3.10.0
##  [4] colorRamps_2.3.4                         
##  [5] clusterCrit_1.3.0                        
##  [6] org.Hs.eg.db_3.21.0                      
##  [7] TxDb.Hsapiens.UCSC.hg38.knownGene_3.21.0 
##  [8] GenomicFeatures_1.61.3                   
##  [9] AnnotationDbi_1.71.0                     
## [10] Biobase_2.69.0                           
## [11] GenomicRanges_1.61.0                     
## [12] GenomeInfoDb_1.45.4                      
## [13] IRanges_2.43.0                           
## [14] S4Vectors_0.47.0                         
## [15] BiocGenerics_0.55.0                      
## [16] generics_0.1.4                           
## [17] geomeTriD.documentation_0.0.3            
## [18] geomeTriD_1.3.11                         
## 
## loaded via a namespace (and not attached):
##   [1] BiocIO_1.19.0               bitops_1.0-9               
##   [3] filelock_1.0.3              tibble_3.3.0               
##   [5] R.oo_1.27.1                 XML_3.99-0.18              
##   [7] rpart_4.1.24                lifecycle_1.0.4            
##   [9] httr2_1.1.2                 aricode_1.0.3              
##  [11] globals_0.18.0              lattice_0.22-7             
##  [13] ensembldb_2.33.0            MASS_7.3-65                
##  [15] backports_1.5.0             magrittr_2.0.3             
##  [17] Hmisc_5.2-3                 sass_0.4.10                
##  [19] rmarkdown_2.29              jquerylib_0.1.4            
##  [21] yaml_2.3.10                 plotrix_3.8-4              
##  [23] Gviz_1.53.0                 DBI_1.2.3                  
##  [25] RColorBrewer_1.1-3          abind_1.4-8                
##  [27] R.utils_2.13.0              AnnotationFilter_1.33.0    
##  [29] biovizBase_1.57.0           RCurl_1.98-1.17            
##  [31] rgl_1.3.18                  nnet_7.3-20                
##  [33] VariantAnnotation_1.55.0    rappdirs_0.3.3             
##  [35] grImport_0.9-7              listenv_0.9.1              
##  [37] parallelly_1.45.0           pkgdown_2.1.3              
##  [39] codetools_0.2-20            DelayedArray_0.35.1        
##  [41] xml2_1.3.8                  tidyselect_1.2.1           
##  [43] UCSC.utils_1.5.0            farver_2.1.2               
##  [45] matrixStats_1.5.0           BiocFileCache_2.99.5       
##  [47] base64enc_0.1-3             GenomicAlignments_1.45.0   
##  [49] jsonlite_2.0.0              trackViewer_1.45.0         
##  [51] progressr_0.15.1            Formula_1.2-5              
##  [53] systemfonts_1.2.3           dbscan_1.2.2               
##  [55] tools_4.5.0                 progress_1.2.3             
##  [57] ragg_1.4.0                  strawr_0.0.92              
##  [59] Rcpp_1.0.14                 glue_1.8.0                 
##  [61] gridExtra_2.3               SparseArray_1.9.0          
##  [63] xfun_0.52                   MatrixGenerics_1.21.0      
##  [65] dplyr_1.1.4                 fastmap_1.2.0              
##  [67] latticeExtra_0.6-30         rhdf5filters_1.21.0        
##  [69] digest_0.6.37               R6_2.6.1                   
##  [71] textshaping_1.0.1           colorspace_2.1-1           
##  [73] jpeg_0.1-11                 dichromat_2.0-0.1          
##  [75] biomaRt_2.65.0              RSQLite_2.4.1              
##  [77] R.methodsS3_1.8.2           data.table_1.17.4          
##  [79] rtracklayer_1.69.0          prettyunits_1.2.0          
##  [81] InteractionSet_1.37.0       httr_1.4.7                 
##  [83] htmlwidgets_1.6.4           S4Arrays_1.9.1             
##  [85] pkgconfig_2.0.3             gtable_0.3.6               
##  [87] blob_1.2.4                  XVector_0.49.0             
##  [89] htmltools_0.5.8.1           ProtGenerics_1.41.0        
##  [91] clue_0.3-66                 scales_1.4.0               
##  [93] png_0.1-8                   knitr_1.50                 
##  [95] rstudioapi_0.17.1           rjson_0.2.23               
##  [97] magic_1.6-1                 checkmate_2.3.2            
##  [99] curl_6.3.0                  cachem_1.1.0               
## [101] rhdf5_2.53.1                stringr_1.5.1              
## [103] parallel_4.5.0              foreign_0.8-90             
## [105] restfulr_0.0.15             desc_1.4.3                 
## [107] pillar_1.10.2               vctrs_0.6.5                
## [109] RANN_2.6.2                  dbplyr_2.5.0               
## [111] cluster_2.1.8.1             htmlTable_2.4.3            
## [113] evaluate_1.0.3              cli_3.6.5                  
## [115] compiler_4.5.0              Rsamtools_2.25.0           
## [117] rlang_1.1.6                 crayon_1.5.3               
## [119] future.apply_1.20.0         interp_1.1-6               
## [121] fs_1.6.6                    stringi_1.8.7              
## [123] deldir_2.0-4                BiocParallel_1.43.3        
## [125] txdbmaker_1.5.4             Biostrings_2.77.1          
## [127] lazyeval_0.2.2              Matrix_1.7-3               
## [129] BSgenome_1.77.0             hms_1.1.3                  
## [131] bit64_4.6.0-1               future_1.58.0              
## [133] ggplot2_3.5.2               Rhdf5lib_1.31.0            
## [135] KEGGREST_1.49.0             SummarizedExperiment_1.39.0
## [137] igraph_2.1.4                memoise_2.0.1              
## [139] bslib_0.9.0                 bit_4.6.0