## ----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)

## ----installation, eval=FALSE-------------------------------------------------
# if (!require("BiocManager", quietly = TRUE)) {
#   install.packages("BiocManager")
# }
# 
# BiocManager::install("geomeTriD")

## ----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)

## ----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)

## ----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")
}

## ----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"
  )
}

## ----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"
  )
}

## ----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)

## ----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)

## ----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))
}

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

## ----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)

## ----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
# )

## -----------------------------------------------------------------------------
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"
)

## ----eval=FALSE---------------------------------------------------------------
# ## not run to reduce the package size.
# library(manipulateWidget)
# clear3d() # Remove the earlier display
# rglViewer(cellobjs, background = "white")
# rglwidget()

## ----threejs------------------------------------------------------------------
threeJsViewer(cellobjs)

## -----------------------------------------------------------------------------
# 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"
  )
}

## -----------------------------------------------------------------------------
for(i in seq_along(objs)){
  if(grepl('tick|arrow|cRE|tss', names(objs)[i])){
    objs[[i]]$layer <- 'bottom'
  }
}
if(interactive()){
  threeJsViewer(objs)
}

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

## -----------------------------------------------------------------------------
## 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))

## ----sessionInfo, results='asis'----------------------------------------------
sessionInfo()