---
title: "geomeTriD Vignette: Plot data in 3D"
author: "Jianhong Ou"
date: "`r BiocStyle::doc_date()`"
package: "`r BiocStyle::pkg_ver('geomeTriD')`"
abstract: >
  Visualize epigeomic data in 2D or 3D plots.
vignette: >
  %\VignetteIndexEntry{geomeTriD Vignette: Plot data in 3D}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
output:
  html_document:
    theme: simplex
    toc: true
    toc_float: true
    toc_depth: 4
    fig_caption: true
---

```{r, echo=FALSE, results="hide", warning=FALSE}
suppressPackageStartupMessages({
  options(rgl.useNULL=TRUE)
  library(geomeTriD)
  library(TxDb.Hsapiens.UCSC.hg19.knownGene)
  library(org.Hs.eg.db)
  library(InteractionSet)
  library(rgl)
})
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
```

# Introduction

The chromatin interactions is involved in precise quantitative and spatiotemporal
control of gene expression. The development of high-throughput experimental 
techniques, such as HiC-seq, HiCAR-seq, and InTAC-seq, for analyzing both the 
higher-order structure of chromatin and the interactions between protein and
their nearby and remote regulatory elements has been developed to reveal how
gene expression is controlled in genome-wide. 
The geomeTriD package enables users to visualize epigenomic data in both 2D and
3D.

# Installation

```{r installation, eval=FALSE}
if (!require("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}

BiocManager::install("geomeTriD")
```

# Quick start

There are three major steps to view the epigenomic data by geomeTriD in 3D.

1. Load data in the `GRanges` object with metadata of positions.

2. Parse the data into a list of `threeJsGeometry` objects.

3. View the data by `threeJsViewer` function by [`threeJs`](https://threejs.org/). .

```{r quickStart}
library(geomeTriD)
## load data
obj <- readRDS(system.file("extdata", "4DNFI1UEG1HD.chr21.FLAMINGO.res.rds",
  package = "geomeTriD"
))
head(obj, n = 2)
feature.gr <- readRDS(system.file("extdata", "4DNFI1UEG1HD.feature.gr.rds",
  package = "geomeTriD"
))
head(feature.gr, n = 2)
## create threeJsGeometry objects
list3d <- view3dStructure(obj,
  feature.gr = feature.gr,
  renderer = "none"
)
## plot the data
threeJsViewer(list3d)
```

if you have trouble in seeing the results, such as [Firefox snap installation issue](https://github.com/karma-runner/karma-firefox-launcher/issues/183), please try to save it as a file and
open it in an explorer.
```{r eval=FALSE}
widget <- threeJsViewer(list3d)
tmpdir <- 'path/to/a/folder'
dir.create(tmpdir, recursive = TRUE)
tmp <- tempfile(tmpdir = tmpdir, fileext = '.html')
htmlwidgets::saveWidget(widget, file=tmp)
utils::browseURL(tmp)
```

# Prepare the inputs

The genomic contact frequency, eg output from HiC,
can be converted into spatial distances and then 
visualized using optimization-based (such as manifold learning techniques)
or probabilistic approaches (such as Markov Chain Monte Carlo).
Available tools include, but not limited to, ShRec3D, GEM-FISH, 
Hierarchical3DGenome, RPR, SuperRec, ShNeigh, PASTIS, and FLAMINGO.

Above example are data generated by [FLAMINGOrLite](https://github.com/JiaxinYangJX/FLAMINGOrLite) from 
the chromosome 21 of the HiC
contact matrix of human GM12878 cells downloaded from [4DN](https://data.4dnucleome.org/files-processed/4DNFI1UEG1HD/#file-overview).

In this package, a simple method to calculate the spatial positions from bin-based contact matrix by multi-dimensional scaling (MDS) algorithm is provided.

## Generate the 3D coordinates from contact matrix

The data in `.hic` or `.cool` can be imported into `GInteractions` objects by
`trackViewer::importGInteractions` function. 
Here we show how to generate the 3D coordinates using the `mdsPlot` function 
from the bin-based contact matrix by Kruskal's
Non-metric Multidimensional Scaling.

```{r propareData}
library(trackViewer)
library(InteractionSet)
# load the interaction data.
# to import your own data, please refer trackViewer help documents
gi_nij <- readRDS(system.file("extdata", "nij.chr6.51120000.53200000.gi.rds",
  package = "geomeTriD"
))
# the data is a GInteractions object with metadata score
head(gi_nij, n = 2)
# define a range to plot
range_chr6 <- GRanges("chr6", IRanges(51120000, 53200000))
# plot it if interactive
if (interactive()) { ## not run to reduce the package size.
  mdsPlot(gi_nij, range = range_chr6, k = 3, render = "threejs")
}
```

## Add gene informations

Here we show how to add gene information along the strand. 
```{r addGene}
# create gene annotations
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(org.Hs.eg.db)
## get all genes
feature.gr <- genes(TxDb.Hsapiens.UCSC.hg19.knownGene)
## subset the data by target viewer region
feature.gr <- subsetByOverlaps(feature.gr, range(regions(gi_nij)))
## assign symbols for each gene
symbols <- mget(feature.gr$gene_id, org.Hs.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)
## Here we replaced the some gene feature as cis-regulatory elements
## to validate how it will be treated by the plot
feature.gr$type <- sample(c("cRE", "gene"),
  length(feature.gr),
  replace = TRUE,
  prob = c(0.1, 0.9)
)
## set the size the cRE in the plot
feature.gr$pch <- rep(NA, length(feature.gr))
feature.gr$pch[feature.gr$type == "cRE"] <- 0.5
## features header
head(feature.gr, n = 2)
# plot it if interactive
if (interactive()) { ## not run to reduce the package size.
  mdsPlot(gi_nij,
    range = range_chr6, feature.gr = feature.gr,
    k = 3, render = "threejs"
  )
}
```

## Add additional signals along the strand

Here we show how to add CTCF ATAC-seq signals along the strand. 

```{r addCTCF}
# one layer of signals. It is a `track` object
ctcf <- readRDS(system.file("extdata", "ctcf.sample.rds",
  package = "trackViewer"
))
# plot it if interactive
if (interactive()) { ## not run to reduce the package size.
  mdsPlot(gi_nij,
    range = range_chr6, feature.gr = feature.gr,
    genomicSigs = ctcf,
    ## we reverse the ATAC-seq signale to show where is open
    reverseGenomicSigs = TRUE,
    k = 3, render = "threejs"
  )
}
```

Add more signals such as ChIP-seq.

```{r addmore}
# create a random signal for demonstration
set.seed(1)
randomSig <- ctcf$dat
randomSig <- randomSig[sort(sample(seq_along(randomSig), 50))]
randomSig$score <- randomSig$score * 2 * runif(50)
head(randomSig, n = 2)
objs <- mdsPlot(gi_nij,
  range = range_chr6, feature.gr = feature.gr,
  genomicSigs = list(ctcf = ctcf, test = randomSig),
  reverseGenomicSigs = c(TRUE, FALSE),
  k = 3, render = "none"
)
threeJsViewer(objs)
```

## view the data by `rgl`

```{r mdsPlot3d, message=FALSE, eval=FALSE}
## not run to reduce the package size.
library(rgl)
library(manipulateWidget)
rgl::clear3d() # Remove the earlier display
rglViewer(objs, background = "white")
rgl::rglwidget() %>%
  rgl::toggleWidget(tags = "tick_minor") %>%
  toggleWidget(tags = "tick_major") %>%
  toggleWidget(tags = "tick_labels") %>%
  toggleWidget(tags = "ctcf") %>%
  toggleWidget(tags = "test") %>%
  toggleWidget(tags = "backbone") %>%
  toggleWidget(tags = "gene_body") %>%
  toggleWidget(tags = "tss_labels") %>%
  toggleWidget(tags = "gene_labels") %>%
  toggleWidget(tags = "cRE") %>%
  rgl::asRow(last = 10)
```

## view the genomic interaction data

Here we will show how to add the interaction pairs into the 3d plot.
```{r jsviewerHic}
## get the backbone 3d positions,
## it will be used as a position reference for the interaction data.
xyz <- extractBackbonePositions(objs)
# make the example easy to see
gi.subset <- gi_nij[distance(first(gi_nij), second(gi_nij))>50000]
hic <- create3dGenomicSignals(
  gi.subset, 
  xyz,
  name='segment', # name prefix for the geometry
  tag='hic', # name for the layer in the scene
  color = c('green', 'darkred'),
  topN=3, # only plot the top 3 events ordered by the scores.
  lwd.maxGenomicSigs = 3,
  alpha=0.5
  )
hic2 <- create3dGenomicSignals(
  gi.subset,
  xyz,
  name='polygon', # name prefix for the geometry
  tag='hic', # name for the layer in the scene
  color = c('blue', 'darkred'),
  topN=seq(4, 10), # only plot the top 4-10 events ordered by the scores.
  type = 'polygon'
  )
width(regions(gi.subset)) <- 101#for high resolution input (such as reads pairs)
hic3 <- create3dGenomicSignals(
  gi.subset, 
  xyz,
  name='segment', # name prefix for the geometry
  tag='hic', # name for the layer in the scene
  color = rainbow(40),
  topN=length(gi.subset),
  lwd.maxGenomicSigs = 5,
  alpha = 0.1
)
# plot it if interactive
if (interactive()) { ## not run to reduce the package size.
  threeJsViewer(c(objs, hic, hic2, hic3))
}
```


# Plot chromatin interactions data as `loopBouquet` or `MDS` plot


Plot the data in 3D using the `mdsPlot` function to generate the 3D coordinates.
Alternatively, users can plot it in 2d. 
Please note that current 2d plot can only handle one genomic signals.

```{r eval=FALSE}
## not run to reduce the package size.
geomeTriD::mdsPlot(gi_nij, range = range_chr6, feature.gr = feature.gr, genomicSigs = ctcf)
```


Different from most of the available tools, `loopBouquetPlot` try to plot the loops with the 2D structure. The nodes indicate the region with interactions and
the edges indicates the interactions. The size of the nodes are relative to the width of the region.
The features could be the cRE or gene. The cRE are shown as
points with symbol 11.

```{r plotGInteractions}
gi_chr2 <- readRDS(system.file("extdata", "gi.rds", package = "trackViewer"))
range_chr2 <- GRanges("chr2", IRanges(234300000, 235000000))
gene_hg19 <- suppressMessages(genes(TxDb.Hsapiens.UCSC.hg19.knownGene))
feature.gr_chr2 <- subsetByOverlaps(gene_hg19, range(regions(gi_chr2)))
feature.gr_chr2$col <- sample(2:7, length(feature.gr_chr2), replace = TRUE)
feature.gr_chr2$type <- sample(c("cRE", "gene"),
  length(feature.gr_chr2),
  replace = TRUE,
  prob = c(0.1, 0.9)
)
feature.gr_chr2$pch <- rep(NA, length(feature.gr_chr2))
feature.gr_chr2$pch[feature.gr_chr2$type == "cRE"] <- 11
symbol <- mget(feature.gr_chr2$gene_id, org.Hs.egSYMBOL, ifnotfound = NA)
symbol <- unlist(lapply(symbol, function(.ele) .ele[1]))
feature.gr_chr2$label <- symbol
geomeTriD::loopBouquetPlot(gi_chr2, range_chr2, feature.gr_chr2)
```

```{r plotRealData, eval=FALSE}
## filter the links to simulate the data with less interactions
keep <- distance(first(gi_nij), second(gi_nij)) > 5e5 & gi_nij$score > 35
gi_nij <- gi_nij[keep]
# narrow the width of anchors to enhance the plots
reg <- regions(gi_nij)
wr <- floor(width(reg) / 4)
start(reg) <- start(reg) + wr
end(reg) <- end(reg) - wr
regions(gi_nij) <- reg
feature.gr <- subsetByOverlaps(gene_hg19, range(regions(gi_nij)))
feature.gr$col <- sample(2:7, length(feature.gr), replace = TRUE)
feature.gr$type <- sample(c("cRE", "gene"),
  length(feature.gr),
  replace = TRUE,
  prob = c(0.1, 0.9)
)
symbol <- mget(feature.gr$gene_id, org.Hs.egSYMBOL, ifnotfound = NA)
symbol <- unlist(lapply(symbol, function(.ele) .ele[1]))
feature.gr$label <- symbol
feature.gr <- c(
  feature.gr[sample(seq_along(feature.gr), 5)],
  feature.gr[feature.gr$type == "cRE"][1]
)
feature.gr <- unique(sort(feature.gr))
geomeTriD::loopBouquetPlot(gi_nij, range_chr6, feature.gr,
  coor_tick_unit = 5e4,
  coor_mark_interval = 5e5,
  genomicSigs = ctcf
)
```

# Plot single cells in 3d

`geomeTriD` can be used to view data in 3D. Here we will show how to use it 
to plot the single cell data in 3D.
The input data should be a data.frame with x,y,z coordinates.

```{r}
library(geomeTriD)
cells <- readRDS(system.file("extdata", "pbmc_small.3d.rds",
  package = "geomeTriD"
))
head(cells)
## color the cells by 'groups' info
cellobjs <- view3dCells(cells,
  x = "umap_1", y = "umap_2", z = "umap_3",
  color = "groups",
  colorFun = function(x) {
    ifelse(x == "g1", "red", "blue")
  },
  renderer = "none"
)
```

Plot the cells by `rgl`.
```{r eval=FALSE}
## not run to reduce the package size.
library(manipulateWidget)
clear3d() # Remove the earlier display
rglViewer(cellobjs, background = "white")
rglwidget()
```

Plot it by `threeJs`.
```{r threejs}
threeJsViewer(cellobjs)
```

We can also try to color the cells by a numeric column.
```{r}
# plot it if interactive
if (interactive()) {
  view3dCells(cells,
    x = "umap_1", y = "umap_2", z = "umap_3",
    color = "nCount_RNA",
    colorFun = function(x, pal = colorRampPalette(c("black", "red"))(5)) {
      limits <- range(x)
      pal[findInterval(x, seq(limits[1], limits[2],
        length.out = length(pal) + 1
      ),
      all.inside = TRUE
      )]
    },
    renderer = "threejs"
  )
}
```

# Advanced usage of `threeJsViewer`

## layout

By default, the viewer built on three.js allows objects in a scene to be separated into layers based on tags. Users can easily toggle the visibility of these layers using the 'toggle' button in the controller. Additionally, layers can be grouped into two main categories (top and bottom), enabling the option to hide internal or auxiliary objects as needed.

In this example, we will move all the ticks and tss markers into bottom group.
Please note that labels will not be affected by this methods.
```{r}
for(i in seq_along(objs)){
  if(grepl('tick|arrow|cRE|tss', names(objs)[i])){
    objs[[i]]$layer <- 'bottom'
  }
}
if(interactive()){
  threeJsViewer(objs)
}
```


The `threeJsViewer` allows users to view two 3D structures side by side. It is recommended to align the 3D coordinates by `alignCoor` function before generating the `threeJsGeometry` objects for optimal comparison.

```{r}
objs2 <- lapply(objs, function(.ele){
  .ele$side <- 'right'
  .ele
})
if(interactive()){
  threeJsViewer(objs, objs2)
}
```

## Create objects from raw

Starting from the raw data, you will be prompted to provide the backbone to which materials will be attached. The backbone is a `GRanges` object containing the x, y, and z coordinates. Since the current resolution of Hi-C data is relatively low, we will first smooth the coordinates. The materials, also input as a `GRanges` object, will have their 3D positions derived from the backbone.

```{r}
## simulate a backbone
library(GenomicRanges)
library(geomeTriD)
backbone <- GRanges(1, IRanges(seq(1, 10000, by=1000), width = 1000))
## set x, y, z position, we will use the torus knots equation
t <- seq(0, 2*pi, length.out=length(backbone))
p <- 3
q <- 2
r <- cos(p*t) + 2
backbone$x <- 2*r*cos(q*t)
backbone$y <- 2*r*sin(q*t)
backbone$z <- -2*sin(p*t)

## simulate the signals to be add
N1 <- 30
sig1 <- GRanges(1, IRanges(sample(seq.int(10000), N1), width = 10),
                score=runif(N1, min=1, max=10))
N2 <- 8
sig2 <- GRanges(1, IRanges(sample(seq.int(10000), N2), width = 1),
                score=runif(N2, min=0.05, max=1),
                shape=sample(c('dodecahedron', 'icosahedron',
                               'octahedron', 'tetrahedron'),
                             N2, replace = TRUE),
                color=sample(seq.int(7)[-1], N2, replace = TRUE))

## do smooth for the coordinates
backbone <- smooth3dPoints(backbone, resolution=10)

## create the line geometry
geo_backbone <- threeJsGeometry(
  x=c(backbone$x0, backbone$x1[length(backbone)]),
  y=c(backbone$y0, backbone$y1[length(backbone)]),
  z=c(backbone$z0, backbone$y1[length(backbone)]),
  colors = "black",
  type = "line",
  tag = "tagNameHere",
  properties= list(size = 2)
)

## create the sphere geometry for sig1
geo_sig1_sphere <- create3dGenomicSignals(
  GenoSig = sig1,
  targetObj = backbone,
  positionTransformFun = function(x) {
      x$y0 <- x$y0 + 0.5 # offset from y
      x$y1 <- x$y1 + 0.5
      x
    }, reverseGenomicSigs = FALSE,
    type = 'sphere',
    tag = 'head', name = 'head',
  color = c('blue', 'red'),
  radius = .25)
## create the cylinder geometry for sig1
geo_sig1_body <- create3dGenomicSignals(
  GenoSig = sig1,
  targetObj = backbone,
  positionTransformFun = function(x) {
      x$y0 <- x$y0 + 0.25 # offset from y
      x$y1 <- x$y1 + 0.25
      x
    }, reverseGenomicSigs = FALSE,
    type = 'cylinder',
    tag = 'body', name = 'body',
  height = .5, radiusTop = .05, radiusBottom = .15)
## create two small spheres
geo_sig1_ear1 <- create3dGenomicSignals(
  GenoSig = sig1,
  targetObj = backbone,
  positionTransformFun = function(x) {
      x$x0 <- x$x0 - 0.25 # offset from x
      x$x1 <- x$x1 - 0.25
      x$y0 <- x$y0 + 0.75 # offset from y
      x$y1 <- x$y1 + 0.75
      x
    }, reverseGenomicSigs = FALSE,
    type = 'sphere',
    tag = 'ear', name = 'ear1',
  color = c('blue', 'red'),
  radius = .1)
geo_sig1_ear2 <- create3dGenomicSignals(
  GenoSig = sig1,
  targetObj = backbone,
  positionTransformFun = function(x) {
      x$x0 <- x$x0 + 0.25 # offset from x
      x$x1 <- x$x1 + 0.25
      x$y0 <- x$y0 + 0.75 # offset from y
      x$y1 <- x$y1 + 0.75
      x
    }, reverseGenomicSigs = FALSE,
    type = 'sphere',
    tag = 'ear', name = 'ear2',
  color = c('blue', 'red'),
  radius = .1)
geo_sig1_nose <- create3dGenomicSignals(
  GenoSig = sig1,
  targetObj = backbone,
  positionTransformFun = function(x) {
      x$y0 <- x$y0 + 0.5 # offset from y
      x$y1 <- x$y1 + 0.5
      x$z0 <- x$z0 - 0.25 # offset from z
      x$z1 <- x$z1 - 0.25
      x
    }, reverseGenomicSigs = FALSE,
    type = 'capsule',
    tag = 'nose', name = 'nose',
  color = c('red', 'blue'),
  radius = .05, height=.05)

## create the icosahedron geometry for sig2
geo_sig2 <- lapply(seq_along(sig2), function(id){
  create3dGenomicSignals(
    GenoSig = sig2[id],
    targetObj = backbone,
    reverseGenomicSigs = FALSE,
    type=sig2[id]$shape,
    tag = 'sig2', name = as.character(id),
    color = sig2[id]$color,
    radius = sig2[id]$score
  )})

## create the a customized geometry by a json file
geo_json <- create3dGenomicSignals(
  GenoSig = GRanges(1, IRanges(5000, width=1), score=1),
  targetObj = backbone,
  type = 'json',
  tag = 'json', name='monkey',
  color = 'green',
  rotation = c(pi, 0, pi), ## we need to rotate the model to change the facing side
  path = system.file('extdata', 'suzanne_buffergeometry.json',
                     package = 'geomeTriD')
)
threeJsViewer(c(backbone=geo_backbone,
                geo_sig1_sphere, geo_sig1_body,
                geo_sig1_ear1, geo_sig1_ear2, geo_sig1_nose,
                geo_sig2, geo_json))
```


# Session Info
```{r sessionInfo, results='asis'}
sessionInfo()
```