# Roadmap to prepare the input matrix for `scTensor`

Koki Tsuyuzaki1, Manabu Ishii1 and Itoshi Nikaido1*

1Laboratory for Bioinformatics Research, RIKEN Center for Biosystems Dynamics Reseach

*k.t.the-answer@hotmail.co.jp

#### 6 July 2020

#### Package

scTensor 1.4.3

# Contents

# 1 Introduction

We explain the way to create a matrix, in which the row names are NCBI Gene ID (ENTREZID), for specifying an input of scTensor. Typical scTensor workflow can be described like below.

    library("scTensor")
    library("LRBase.XXX.eg.db") # e.g. LRBase.Hsa.eg.db
    # Input matrix
    input <- ...
    sce <- SingleCellExperiment(assays=list(counts = input))
    # Color vector
    color <- ...
    # Celltype vector
    label <- ...
    cellCellSetting(sce, LRBase.XXX.eg.db, color, label)

In scTensor, the row names of the input matrix is limited only NCBI Gene ID to cooperate with the other R packages (cf. data(“Germline”)). Since the user has many different types of the data matrix, here we introduce some situations and the way to convert the row names of the user’s matrix as NCBI Gene ID.

# 2 Step.1 : Create gene-level expression matrix

First of all, we have to prepare the expression matrix (gene \(\times\) cell). There are many types of single-cell RNA-Seq (scRNA-Seq) technologies and the situation will be changed by the used experimental methods and quantification tools described below.

## 2.1 Case I: Gene-level quantification

In Plate-based scRNA-Seq experiment (i.e. Smart-Seq2, Quart-Seq2, CEL-Seq2, MARS-Seq,…etc), FASTQ file is generated in each cell. After the mapping of reads in each FASTQ file to reference genome, the same number of BAM files will be generated. By using some quantification tools such as featureCounts, or HTSeq-count, user can get the expression matrix and used as the input of scTensor. These tools simply count the number of reads in union-exon in each gene. However, these tools do not take “multimap” of not unique reads into consideration and the quantification is not accurate. Therefore, we recommend the transcript-level quantification and gene-level summarization explained below.

## 2.2 Case II : Transcript-level quantification

Some quantification tools such as RSEM, Sailfish, Salmon, Kallisto, and StringTie use the reference transcriptome instead of genome, and quantify the expression level in each transcript. After the quantification, the transcript-level expression can be summarized to gene-level by using `summarizeToGene` function of tximport. The paper of tximport showed that the gene-level summalization is more accurate than featureCounts.

Note that if you use the reference transcriptome of GENCODE, this step becomes slightly complicated. Although the number of transcripts of GENCODE and that of Ensembl is almost the same, and actually, most of the transcript is duplicated in these two databases, the gene identifier used in GENCODE looks complicated like “ENST00000456328.2|ENSG00000223972.5|OTTHUMG00000000961.2|OTTHUMT00000362751.1|DDX11L1-202|DDX11L1|1657|processed_transcript|”. In such case, firstly only Ensembl Transcript ID should be extracted (e.g. ENST00000456328.2), removed the version (e.g. ENST00000456328), summarized to Ensembl Gene ID by tximport (e.g. ENSG00000223972), and then converted to NCBI Gene ID (e.g. 100287102).

## 2.3 Case III: UMI-count

In the droplet-based scRNA-Seq experiment (i.e. Drop-Seq, inDrop RNA-Seq, 10X Chromium), unique molecular identifier (UMI) is introduced for avoiding the bias of PCR amplification, and after multiplexing by cellular barcode, digital gene expression (DGE) matrix is generated by counting the number of types of UMI mapped in each gene.

When user perform Drop-seq, Drop-Seq tool can generate the DGE matrix.

Another tool Alevin, which is a subcommand of Salmon is also applicable to Drop-seq data. In such case [tximport] with option “type = ‘alevin’” can import the result of Alevin into R and summarize the DGE matrix.

When the user performs 10X Chromium, using Cell Ranger developed by 10X Genomics is straightforward.

Although Cell Ranger is implemented by Python, starting from the files generated by Cell Ranger (e.g. filtered_gene_bc_matrices/{hg19,mm10}/{barcodes.tsv,genes.tsv,matrix.mtx}), Seurat can import these files to an R object.

For example, according to the tutorial of Seurat (Seurat - Guided Clustering Tutorial), PBMC data of Cell Ranger can be imported by `Read10X` function and DGE matrix of UMI-count is available by the output of `CreateSeuratObject` function.

    library("Seurat")

    # Load the PBMC dataset
    pbmc.data <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19/")

    # Initialize the Seurat object with the raw (non-normalized data).
    pbmc <- CreateSeuratObject(counts = pbmc.data,
      project = "pbmc3k", min.cells = 3, min.features = 200)

Note that the matrix is formatted as a sparse matrix of Matrix package (MM: Matrix market), but the scTensor assumes dense matrix. By using `as.matrix` function, the sparse matrix is easily converted to dense matrix as follows.

    # Sparse matrix to dense matrix
    for_sc <- as.matrix(pbmc.data)

# 3 Step.2: Convert the row names of a matrix as NCBI Gene ID (ENTREZID)

Even after creating the gene-level expression matrix in Step.1, many kinds of gene-level gene identifiers can be assigned as row names of the matrix such as Ensembl Gene ID, RefSeq, or Gene Symbol. Again, only NCBI Gene ID can be used as row names of input matrix of scTensor. To do such a task, we originally implemented a function `convertToNCBIGeneID`. The only user has to prepare for using this function is the 1. input matrix (or data.frame) filled with only numbers, 2. current gene-level gene identifier in each row of the input matrix, and 3. corresponding table containing current gene-level gene identifier (left) and corresponding NCBI Gene ID (right). The usage of this function is explained below.

## 3.1 Case I: Ensembl Gene ID to NCBI Gene ID

In addition to 1. and 2., the user has to prepare the 3. corresponding table. Here we introduce two approaches to assign user’s Ensembl Gene ID to NCBI Gene ID. First approarch is using Organism DB packages such as Homo.sapiens, Mus.musculus, and Rattus.norvegicus.

Using the `select` function of Organism DB, the corresponding table can be retrieved like below.

    suppressPackageStartupMessages(library("scTensor"))
    suppressPackageStartupMessages(library("Homo.sapiens"))

    # 1. Input matrix
    input <- matrix(1:20, nrow=4, ncol=5)
    # 2. Gene identifier in each row
    rowID <- c("ENSG00000204531", "ENSG00000181449",
      "ENSG00000136997", "ENSG00000136826")
    # 3. Corresponding table
    LefttoRight <- select(Homo.sapiens,
      column=c("ENSEMBL", "ENTREZID"),
      keytype="ENSEMBL", keys=rowID)

    ## 'select()' returned 1:1 mapping between keys and columns

    # ID conversion
    (input <- convertToNCBIGeneID(input, rowID, LefttoRight))

    ##      [,1] [,2] [,3] [,4] [,5]
    ## 5460    1    5    9   13   17
    ## 6657    2    6   10   14   18
    ## 4609    3    7   11   15   19
    ## 9314    4    8   12   16   20

Second approarch is using AnnotationHub package.

Although only three Organism DB packages are explicitly developed, even if the data is generated from other species (e.g. Zebrafish, Arabidopsis thaliana), similar database is also available from AnnotationHub, and `select` function can be performed like below.

    suppressPackageStartupMessages(library("AnnotationHub"))

    # 1. Input matrix
    input <- matrix(1:20, nrow=4, ncol=5)
    # 3. Corresponding table
    ah <- AnnotationHub()

    ## snapshotDate(): 2020-04-27

    # Database of Human
    hs <- query(ah, c("OrgDb", "Homo sapiens"))[[1]]

    ## loading from cache

    LefttoRight <- select(hs,
      column=c("ENSEMBL", "ENTREZID"),
      keytype="ENSEMBL", keys=rowID)

    ## 'select()' returned 1:1 mapping between keys and columns

    (input <- convertToNCBIGeneID(input, rowID, LefttoRight))

    ##      [,1] [,2] [,3] [,4] [,5]
    ## 5460    1    5    9   13   17
    ## 6657    2    6   10   14   18
    ## 4609    3    7   11   15   19
    ## 9314    4    8   12   16   20

## 3.2 Case II: Gene Symbol to NCBI Gene ID

When using cellranger or Seurat to quantify UMI-count (cf. Step1, Case III), the row names of input matrix might be Gene Symbol, and have to be converted to NCBI Gene ID. As well as the Case I described above, Organism DB and AnnotationHub will support such a task like below.

    # 1. Input matrix
    input <- matrix(1:20, nrow=4, ncol=5)
    # 2. Gene identifier in each row
    rowID <- c("POU5F1", "SOX2", "MYC", "KLF4")
    # 3. Corresponding table
    LefttoRight <- select(Homo.sapiens,
      column=c("SYMBOL", "ENTREZID"),
      keytype="SYMBOL", keys=rowID)

    ## 'select()' returned 1:1 mapping between keys and columns

    # ID conversion
    (input <- convertToNCBIGeneID(input, rowID, LefttoRight))

    ##      [,1] [,2] [,3] [,4] [,5]
    ## 5460    1    5    9   13   17
    ## 6657    2    6   10   14   18
    ## 4609    3    7   11   15   19
    ## 9314    4    8   12   16   20

    # 1. Input matrix
    input <- matrix(1:20, nrow=4, ncol=5)
    # 3. Corresponding table
    ah <- AnnotationHub()

    ## snapshotDate(): 2020-04-27

    # Database of Human
    hs <- query(ah, c("OrgDb", "Homo sapiens"))[[1]]

    ## loading from cache

    LefttoRight <- select(hs,
      column=c("SYMBOL", "ENTREZID"),
      keytype="SYMBOL", keys=rowID)

    ## 'select()' returned 1:1 mapping between keys and columns

    (input <- convertToNCBIGeneID(input, rowID, LefttoRight))

    ##      [,1] [,2] [,3] [,4] [,5]
    ## 5460    1    5    9   13   17
    ## 6657    2    6   10   14   18
    ## 4609    3    7   11   15   19
    ## 9314    4    8   12   16   20

# 4 Step.3: Normalize the count matrix

Finally, we introduce some situations to perform some normalization methods of gene expression matrix.

If a user convert a Seurat object to a SingleCellExperient object by using `as.SingleCellExperiment`, the result of `NormalizeData` function (logcounts) is inherited to the SingleCellExperient object as follows;

    pbmc2 <- NormalizeData(pbmc, normalization.method = "LogNormalize",
        scale.factor = 10000)
    sce <- as.SingleCellExperiment(pbmc2)
    assayNames(sce) # counts, logcounts

If the user want to use scater package, `calculateCPM` or `normalize` function can calculate the normalized expression values as follows; (see also the vignette of scater).

    library("scater")
    sce <- SingleCellExperiment(assays=list(counts = input))
    cpm(sce) <- calculateCPM(sce)
    sce <- normalize(sce)
    assayNames(sce) # counts, normcounts, logcounts, cpm

Actually, any original normalization can be stored in the sce. For example, we can calculate the value of count per median (CPMED) as follows;

    # User's Original Normalization Function
    CPMED <- function(input){
        libsize <- colSums(input)
        median(libsize) * t(t(input) / libsize)
    }
    # Normalization
    normcounts(sce) <- log10(CPMED(counts(sce)) + 1)

We recommend using the normcounts slot to save such original normalization values. After the normalization, such values can be specified by `assayNames` option in `cellCellRanks` `cellCellDecomp` and `cellCellReport` functions.

# Session information

    ## R version 4.0.0 (2020-04-24)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server 2012 R2 x64 (build 9600)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.1252
    ## [3] LC_MONETARY=English_United States.1252
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.1252
    ##
    ## attached base packages:
    ## [1] parallel  stats4    stats     graphics  grDevices utils     datasets
    ## [8] methods   base
    ##
    ## other attached packages:
    ##  [1] AnnotationHub_2.20.0
    ##  [2] BiocFileCache_1.12.0
    ##  [3] dbplyr_1.4.4
    ##  [4] Homo.sapiens_1.3.1
    ##  [5] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
    ##  [6] org.Hs.eg.db_3.11.4
    ##  [7] GO.db_3.11.4
    ##  [8] OrganismDbi_1.30.0
    ##  [9] GenomicFeatures_1.40.0
    ## [10] AnnotationDbi_1.50.1
    ## [11] MeSH.Mmu.eg.db_1.13.0
    ## [12] LRBase.Mmu.eg.db_1.2.0
    ## [13] MeSH.Hsa.eg.db_1.13.0
    ## [14] MeSHDbi_1.24.0
    ## [15] SingleCellExperiment_1.10.1
    ## [16] SummarizedExperiment_1.18.1
    ## [17] DelayedArray_0.14.0
    ## [18] matrixStats_0.56.0
    ## [19] Biobase_2.48.0
    ## [20] GenomicRanges_1.40.0
    ## [21] GenomeInfoDb_1.24.2
    ## [22] IRanges_2.22.2
    ## [23] S4Vectors_0.26.1
    ## [24] BiocGenerics_0.34.0
    ## [25] scTensor_1.4.3
    ## [26] RSQLite_2.2.0
    ## [27] LRBase.Hsa.eg.db_1.2.0
    ## [28] LRBaseDbi_1.6.0
    ## [29] BiocStyle_2.16.0
    ##
    ## loaded via a namespace (and not attached):
    ##   [1] rsvd_1.0.3                    Hmisc_4.4-0
    ##   [3] ica_1.0-2                     Rsamtools_2.4.0
    ##   [5] foreach_1.5.0                 lmtest_0.9-37
    ##   [7] crayon_1.3.4                  MASS_7.3-51.6
    ##   [9] nlme_3.1-148                  backports_1.1.8
    ##  [11] GOSemSim_2.14.0               rlang_0.4.6
    ##  [13] XVector_0.28.0                ROCR_1.0-11
    ##  [15] irlba_2.3.3                   nnTensor_1.0.5
    ##  [17] GOstats_2.54.0                BiocParallel_1.22.0
    ##  [19] tagcloud_0.6                  bit64_0.9-7
    ##  [21] glue_1.4.1                    sctransform_0.2.1
    ##  [23] dotCall64_1.0-0               DOSE_3.14.0
    ##  [25] tidyselect_1.1.0              fitdistrplus_1.1-1
    ##  [27] XML_3.99-0.3                  tidyr_1.1.0
    ##  [29] zoo_1.8-8                     GenomicAlignments_1.24.0
    ##  [31] xtable_1.8-4                  magrittr_1.5
    ##  [33] evaluate_0.14                 ggplot2_3.3.2
    ##  [35] zlibbioc_1.34.0               rstudioapi_0.11
    ##  [37] rpart_4.1-15                  fastmatch_1.1-0
    ##  [39] ensembldb_2.12.1              maps_3.3.0
    ##  [41] fields_10.3                   shiny_1.5.0
    ##  [43] xfun_0.15                     askpass_1.1
    ##  [45] cluster_2.1.0                 caTools_1.18.0
    ##  [47] tidygraph_1.2.0               TSP_1.1-10
    ##  [49] tibble_3.0.1                  interactiveDisplayBase_1.26.3
    ##  [51] ggrepel_0.8.2                 biovizBase_1.36.0
    ##  [53] ape_5.4                       listenv_0.8.0
    ##  [55] dendextend_1.13.4             Biostrings_2.56.0
    ##  [57] png_0.1-7                     future_1.17.0
    ##  [59] bitops_1.0-6                  ggforce_0.3.2
    ##  [61] RBGL_1.64.0                   plyr_1.8.6
    ##  [63] GSEABase_1.50.1               AnnotationFilter_1.12.0
    ##  [65] pillar_1.4.4                  gplots_3.0.4
    ##  [67] graphite_1.34.0               europepmc_0.4
    ##  [69] vctrs_0.3.1                   ellipsis_0.3.1
    ##  [71] generics_0.0.2                plot3D_1.3
    ##  [73] urltools_1.7.3                MeSH.Aca.eg.db_1.13.0
    ##  [75] outliers_0.14                 tools_4.0.0
    ##  [77] foreign_0.8-80                entropy_1.2.1
    ##  [79] munsell_0.5.0                 tweenr_1.0.1
    ##  [81] fgsea_1.14.0                  fastmap_1.0.1
    ##  [83] compiler_4.0.0                abind_1.4-5
    ##  [85] httpuv_1.5.4                  rtracklayer_1.48.0
    ##  [87] Gviz_1.32.0                   plotly_4.9.2.1
    ##  [89] GenomeInfoDbData_1.2.3        gridExtra_2.3
    ##  [91] lattice_0.20-41               visNetwork_2.0.9
    ##  [93] AnnotationForge_1.30.1        later_1.1.0.1
    ##  [95] dplyr_1.0.0                   jsonlite_1.7.0
    ##  [97] concaveman_1.1.0              scales_1.1.1
    ##  [99] graph_1.66.0                  pbapply_1.4-2
    ## [101] genefilter_1.70.0             lazyeval_0.2.2
    ## [103] promises_1.1.1                MeSH.db_1.13.0
    ## [105] latticeExtra_0.6-29           reticulate_1.16
    ## [107] checkmate_2.0.0               rmarkdown_2.3
    ## [109] cowplot_1.0.0                 schex_1.2.0
    ## [111] MeSH.Syn.eg.db_1.13.0         webshot_0.5.2
    ## [113] Rtsne_0.15                    dichromat_2.0-0
    ## [115] BSgenome_1.56.0               uwot_0.1.8
    ## [117] igraph_1.2.5                  gclus_1.3.2
    ## [119] survival_3.2-3                yaml_2.2.1
    ## [121] plotrix_3.7-8                 htmltools_0.5.0
    ## [123] memoise_1.1.0                 VariantAnnotation_1.34.0
    ## [125] rTensor_1.4.1                 Seurat_3.1.5
    ## [127] seriation_1.2-8               graphlayouts_0.7.0
    ## [129] viridisLite_0.3.0             digest_0.6.25
    ## [131] assertthat_0.2.1              ReactomePA_1.32.0
    ## [133] mime_0.9                      rappdirs_0.3.1
    ## [135] registry_0.5-1                spam_2.5-1
    ## [137] future.apply_1.6.0            misc3d_0.8-4
    ## [139] data.table_1.12.8             blob_1.2.1
    ## [141] cummeRbund_2.30.0             splines_4.0.0
    ## [143] Formula_1.2-3                 ProtGenerics_1.20.0
    ## [145] RCurl_1.98-1.2                hms_0.5.3
    ## [147] colorspace_1.4-1              base64enc_0.1-3
    ## [149] BiocManager_1.30.10           nnet_7.3-14
    ## [151] Rcpp_1.0.4.6                  bookdown_0.20
    ## [153] RANN_2.6.1                    MeSH.PCR.db_1.13.0
    ## [155] enrichplot_1.8.1              R6_2.4.1
    ## [157] grid_4.0.0                    ggridges_0.5.2
    ## [159] lifecycle_0.2.0               acepack_1.4.1
    ## [161] curl_4.3                      MeSH.Bsu.168.eg.db_1.13.0
    ## [163] gdata_2.18.0                  leiden_0.3.3
    ## [165] MeSH.AOR.db_1.13.0            meshr_1.24.1
    ## [167] DO.db_2.9                     Matrix_1.2-18
    ## [169] qvalue_2.20.0                 RcppAnnoy_0.0.16
    ## [171] RColorBrewer_1.1-2            iterators_1.0.12
    ## [173] stringr_1.4.0                 htmlwidgets_1.5.1
    ## [175] polyclip_1.10-0               triebeard_0.3.0
    ## [177] biomaRt_2.44.1                purrr_0.3.4
    ## [179] gridGraphics_0.5-0            reactome.db_1.70.0
    ## [181] globals_0.12.5                openssl_1.4.2
    ## [183] htmlTable_2.0.1               patchwork_1.0.1
    ## [185] codetools_0.2-16              gtools_3.8.2
    ## [187] prettyunits_1.1.1             gtable_0.3.0
    ## [189] tsne_0.1-3                    DBI_1.1.0
    ## [191] highr_0.8                     httr_1.4.1
    ## [193] KernSmooth_2.23-17            stringi_1.4.6
    ## [195] progress_1.2.2                reshape2_1.4.4
    ## [197] farver_2.0.3                  heatmaply_1.1.0
    ## [199] annotate_1.66.0               viridis_0.5.1
    ## [201] hexbin_1.28.1                 fdrtool_1.2.15
    ## [203] Rgraphviz_2.32.0              magick_2.4.0
    ## [205] xml2_1.3.2                    rvcheck_0.1.8
    ## [207] ggplotify_0.0.5               Category_2.54.0
    ## [209] BiocVersion_3.11.1            bit_1.1-15.2
    ## [211] scatterpie_0.1.4              jpeg_0.1-8.1
    ## [213] ggraph_2.0.3                  pkgconfig_2.0.3
    ## [215] knitr_1.29