---
title: "R Package for Analyzing Tomo-seq Data"
author: "Ryosuke Matsuzawa"
output: BiocStyle::html_document
vignette: >
    %\VignetteIndexEntry{tomoseqr}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8}
---

# Introduction

Tomo-seq is a genome-wide RNA tomography method that combines combining
high-throughput RNA sequencing with cryosectioning for spatially resolved
transcriptomics.
In 2014, [Junker et al.](https://doi.org/10.1016/j.cell.2014.09.038) published a method for reconstructing 3D gene
expression patterns from tomoseq data (iterative proportional fitting ([Fienberg, 1970](https://doi.org/10.1214/aoms/1177696968))). tomoseqr is a package to run
this method in R. tomoseqr also includes a 3D silhouette maker
for samples and tools to visualize the reconstruction results.

# Installation and loading

Install the `tomoseqr` package with the following command.
```{r eval = FALSE}
# for bioconductor
if (!require("BiocManager"))
    install.packages("BiocManager")

BiocManager::install("tomoseqr")
```

Load the `tomoseqr` package with the following command.

```{r setup, message = FALSE}
library(tomoseqr)
```

# Data preparation

`tomoseqr` requires Tomo-seq data in three directions (x, y, and z-axes) and a mask.

## Tomo-seq data preparation

Tomo-seq data must meet the following conditions.


1. It is a data frame with a header.
1. The first column is the gene ID (string), and the second and subsequent
columns are the expression levels (numeric) at each sections.
1. The order of the second and subsequent rows is the same as the order of
the sections.
1. The second and subsequent rows are ordered from left to right
(not right to left) of the intercept for each axis.
1. One row corresponds to one gene.

The `tomoseqr` package provides artificial tomoseq data as an example.

```{r}
data(testx, testy, testz)
head(testx)
```

## Tomo-seq data from Junker 2014

You can download part of the Tomo-seq data published by [Junker et al.](https://www.sciencedirect.com/science/article/pii/S0092867414012264#mmc3) with following commands. 

```{r, eval = FALSE}
tomoCache <- downloadJunker2014()
junker2014 <- doadJunker2014(tomoCache)

sheldAV <- junker2014[["sheldAV"]]
sheldVD <- junker2014[["sheldVD"]]
sheldLR <- junker2014[["sheldLR"]]
mask <- junker2014[["mask"]]
```

## Mask data preparation

Mask is a 3-dimensional matrix that represents the shape of the sample.
The value of each element indicates whether a voxel is included in the
sample or not. For example, if the voxel at coordinates (2, 3, 4) is
contained in the material, then The value of the element (2, 3, 4) in
the mask is 1. Conversely, the voxel at coordinates (10, 11, 12) is not
included in the sample, then the value of the mask element (10, 11, 12)
is 0.

The `tomoseqr` package provides artificial mask data as an example.

```{r}
data(mask)
```

## How to use Masker (tomoseqr mask maker)

The `tomoseqr` package provides a tool to create mask data,
which can be run with the following command.

```{r eval = FALSE}
masker()
```

![Interface of masker](masker_image.png)

You can create mask data using masker by create cross sections along
certain axes.

1. To create a new mask, select `new mask` under `How input`.
1. Set the number of sections for each axis in the
`Number of intercept` field.
1. In `width of cells`, set the thickness of the intercept for each axis
as a relative value.
(This is only related to the appearance of the mask creation screen, so
it does not have to be strict.)
1. Set along which axis you create the mask in the `Along which axis?`
field.
1. Click `make table`, and a mask creation screen appears.
1. On the mask creation screen, create a cross-sectional view
as if you were creating pixel art.
1. Click `Download as .rda` to output mask data as rda file.

# Reconstruction of 3D expression patterns

Use the `estimate3dExpressions` function to reconstruct the 3D expression
pattern of the gene specified by `query` from Tomo-seq data.

```{r}
tomoObj <- estimate3dExpressions(
    testx,
    testy,
    testz,
    mask = mask,
    query = c("gene1", "gene2", "gene3")
)
```

# Visualize result of reconstruction

Visualize the reconstruction results using `imageViewer`.

```{r eval = FALSE}
imageViewer(tomoObj)
```

This command starts a shiny app `imageViewer`.

![Interface of ImageViewer (2D view)](2Dimage.png)

![Interface of ImageViewer (3D view)](3Dimage.png)

1. Select gene ID.
1. Edit labels of each axes.
1. Set the width of sections for each axis. it is used to calculate the
aspect ratio of figures.
1. Use this tab to switch between 2D and 3D views.
1. Select which axis to look at the cross section perpendicular to.
1. Move the slider to select the location you want to see. Press the Auto Play button to view each section of the image
frame by frame.
1. Edit figure titles and labels.
1. Specify the range of expression levels to be displayed (default is
the full range).
1. Rendering area.
1. Obtain a cross-sectional view (png) or a frame-by-frame image (gif)
of each fragment being displayed.
1. Threshold of expression level. Expression points with levels lower than
this value are not plotted.
1. Toggle whether or not to plot points representing the mask.

# Find axial peak gene

Detect genes with significant differences in expression levels between
sections from Tomo-seq data.

```{r}
axialGeneTable <- findAxialGenes(testx)
print(axialGeneTable)
```
The results include these columns:

- GeneID: Gene ID.  
- max: The maximum expression level of the gene.  
- meanExeptMax: Average of expression levels excluding the maximum.  
- index: Number of the section containing the maximum expression level.

You can also specify a vector of candidate gene IDs.

```{r}
axialGeneTable <- findAxialGenes(testx, genes = c("gene1", "gene3"))
print(axialGeneTable)
```

# SessionInfo

```{r}
sessionInfo()
```