1 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. published a method for reconstructing 3D gene expression patterns from tomoseq data (iterative proportional fitting (Fienberg, 1970)). 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.

2 Installation and loading

Install the tomoseqr package with the following command.

# for bioconductor
if (!require("BiocManager"))
    install.packages("BiocManager")

BiocManager::install("tomoseqr")

Load the tomoseqr package with the following command.

library(tomoseqr)

3 Data preparation

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

3.1 Tomo-seq data preparation

Tomo-seq data must meet the following conditions.

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

The tomoseqr package provides artificial tomoseq data as an example.

data(testx, testy, testz)
head(testx)

## # A tibble: 4 × 51
##   geneID    V1    V2    V3    V4    V5    V6    V7    V8    V9   V10   V11   V12
##   <chr>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 gene1      0     0     0     0     0     0     0     0     0     0     0     0
## 2 gene2    230   470   650   850  1010  1170  1370  1510  1650  1730  1910  2030
## 3 gene3      0     0     0     0     0     0     0     0     0     0   100   900
## 4 gene4  29770 29530 29350 29150 28990 28830 28630 28490 28350 28270 27990 27070
## # ℹ 38 more variables: V13 <dbl>, V14 <dbl>, V15 <dbl>, V16 <dbl>, V17 <dbl>,
## #   V18 <dbl>, V19 <dbl>, V20 <dbl>, V21 <dbl>, V22 <dbl>, V23 <dbl>,
## #   V24 <dbl>, V25 <dbl>, V26 <dbl>, V27 <dbl>, V28 <dbl>, V29 <dbl>,
## #   V30 <dbl>, V31 <dbl>, V32 <dbl>, V33 <dbl>, V34 <dbl>, V35 <dbl>,
## #   V36 <dbl>, V37 <dbl>, V38 <dbl>, V39 <dbl>, V40 <dbl>, V41 <dbl>,
## #   V42 <dbl>, V43 <dbl>, V44 <dbl>, V45 <dbl>, V46 <dbl>, V47 <dbl>,
## #   V48 <dbl>, V49 <dbl>, V50 <dbl>

3.2 Tomo-seq data from Junker 2014

You can download part of the Tomo-seq data published by Junker et al. with following commands.

tomoCache <- downloadJunker2014()
junker2014 <- doadJunker2014(tomoCache)

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

3.3 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.

data(mask)

3.4 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.

masker()

Interface of masker

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

To create a new mask, select new mask under How input.
Set the number of sections for each axis in the Number of intercept field.
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.)
Set along which axis you create the mask in the Along which axis? field.
Click make table, and a mask creation screen appears.
On the mask creation screen, create a cross-sectional view as if you were creating pixel art.
Click Download as .rda to output mask data as rda file.

4 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.

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

5 Visualize result of reconstruction

Visualize the reconstruction results using imageViewer.

imageViewer(tomoObj)

This command starts a shiny app imageViewer.

Interface of ImageViewer (2D view)

Interface of ImageViewer (3D view)

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

6 Find axial peak gene

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

axialGeneTable <- findAxialGenes(testx)
print(axialGeneTable)

## # A tibble: 4 × 4
##   geneID   max meanExeptMax isPeakGene
##   <chr>  <dbl>        <dbl>      <dbl>
## 1 gene1   3100        551.           0
## 2 gene2   5560       2334.           0
## 3 gene3    900         22.4         12
## 4 gene4  29990      26898.           0

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.

axialGeneTable <- findAxialGenes(testx, genes = c("gene1", "gene3"))
print(axialGeneTable)

## # A tibble: 2 × 4
##   geneID   max meanExeptMax isPeakGene
##   <chr>  <dbl>        <dbl>      <dbl>
## 1 gene1   3100        551.           0
## 2 gene3    900         22.4         12

7 SessionInfo

sessionInfo()

## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [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: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] tomoseqr_1.6.0   BiocStyle_2.30.0
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.4         xfun_0.40            bslib_0.5.1         
##  [4] ggplot2_3.4.4        htmlwidgets_1.6.2    tzdb_0.4.0          
##  [7] vctrs_0.6.4          tools_4.3.1          generics_0.1.3      
## [10] curl_5.1.0           tibble_3.2.1         fansi_1.0.5         
## [13] RSQLite_2.3.1        blob_1.2.4           pkgconfig_2.0.3     
## [16] data.table_1.14.8    dbplyr_2.3.4         lifecycle_1.0.3     
## [19] stringr_1.5.0        compiler_4.3.1       munsell_0.5.0       
## [22] httpuv_1.6.12        htmltools_0.5.6.1    sass_0.4.7          
## [25] yaml_2.3.7           lazyeval_0.2.2       plotly_4.10.3       
## [28] pillar_1.9.0         later_1.3.1          jquerylib_0.1.4     
## [31] tidyr_1.3.0          ellipsis_0.3.2       cachem_1.0.8        
## [34] mime_0.12            tidyselect_1.2.0     digest_0.6.33       
## [37] stringi_1.7.12       dplyr_1.1.3          purrr_1.0.2         
## [40] animation_2.7        bookdown_0.36        fastmap_1.1.1       
## [43] grid_4.3.1           colorspace_2.1-0     cli_3.6.1           
## [46] magrittr_2.0.3       utf8_1.2.4           withr_2.5.1         
## [49] readr_2.1.4          filelock_1.0.2       scales_1.2.1        
## [52] promises_1.2.1       bit64_4.0.5          rmarkdown_2.25      
## [55] httr_1.4.7           bit_4.0.5            hms_1.1.3           
## [58] memoise_2.0.1        shiny_1.7.5.1        evaluate_0.22       
## [61] knitr_1.44           BiocFileCache_2.10.0 viridisLite_0.4.2   
## [64] rlang_1.1.1          Rcpp_1.0.11          xtable_1.8-4        
## [67] glue_1.6.2           DBI_1.1.3            BiocManager_1.30.22 
## [70] jsonlite_1.8.7       R6_2.5.1

R Package for Analyzing Tomo-seq Data

24 October 2023

Contents