scTensor
scTensor 2.17.0
This vignette has been changed in BioC 3.14, when each data package (LRBase.XXX.eg.db) is deprecated and the way to provide LRBase data has changed to AnnotationHub-style.
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 as below.
library("scTensor")
library("AnnotationHub")
library("LRBaseDbi")
# Input matrix
input <- ...
sce <- SingleCellExperiment(assays=list(counts = input))
# Celltype vector
label <- ...
# LRBase.XXX.eg.db
ah <- AnnotationHub()
dbfile <- query(ah, c("LRBaseDb", "Homo sapiens", "v002"))[[1]]
LRBase.Hsa.eg.db <- LRBaseDbi::LRBaseDb(dbfile)
# Setting
cellCellSetting(sce, LRBase.Hsa.eg.db, 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.
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.
In Plate-based scRNA-Seq experiment (i.e. Smart-Seq2, Quart-Seq2, CEL-Seq2, MARS-Seq,…etc), the FASTQ file is generated in each cell. After the mapping of reads in each FASTQ file to the 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. One downside of these tools is that such 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.
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 transcript-level expression summalized to gene-level is more accurate than the gene-level expression calculated by 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 a 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) by each organism package such as Homo.sapiens.
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 the Read10X
function and DGE matrix of UMI-count is available by the output of CreateSeuratObject
function.
if(!require(Seurat)){
BiocManager::install("Seurat")
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 for now.
By using as.matrix
function,
the sparse matrix is easily converted to a dense matrix as follows.
# Sparse matrix to dense matrix
for_sc <- as.matrix(pbmc.data)
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 the input matrix of scTensor.
To do such a task, we originally implemented a function convertRowID
function of scTGIF.
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.
In addition to 1. and 2., the user has to prepare the 3. corresponding table. Here we introduce two approaches to assign the 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"))
if(!require(Homo.sapiens)){
BiocManager::install("Homo.sapiens")
suppressPackageStartupMessages(library(Homo.sapiens))
}
## Loading required package: Homo.sapiens
## Loading required package: AnnotationDbi
## Loading required package: OrganismDbi
## Loading required package: GenomicFeatures
## Loading required package: GO.db
## Loading required package: org.Hs.eg.db
##
## Loading required package: TxDb.Hsapiens.UCSC.hg19.knownGene
if(!require(scTGIF)){
BiocManager::install("scTGIF")
suppressPackageStartupMessages(library(scTGIF))
}
## Loading required package: scTGIF
# 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 <- convertRowID(input, rowID, LefttoRight))
##
|
| | 0%
|
|================== | 25%
|
|=================================== | 50%
|
|==================================================== | 75%
|
|======================================================================| 100%
## Input matrix: 4x5
## Output matrix: 4x5
## Some gene expression vectors were collapsed into single vector
## by sum rule
## $output
## [,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
##
## $ctable
## Left Right
## [1,] "ENSG00000204531" "5460"
## [2,] "ENSG00000181449" "6657"
## [3,] "ENSG00000136997" "4609"
## [4,] "ENSG00000136826" "9314"
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()
# 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 <- convertRowID(input, rowID, LefttoRight))
##
|
| | 0%
|
|================== | 25%
|
|=================================== | 50%
|
|==================================================== | 75%
|
|======================================================================| 100%
## Input matrix: 4x5
## Output matrix: 4x5
## Some gene expression vectors were collapsed into single vector
## by sum rule
## $output
## [,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
##
## $ctable
## Left Right
## [1,] "ENSG00000204531" "5460"
## [2,] "ENSG00000181449" "6657"
## [3,] "ENSG00000136997" "4609"
## [4,] "ENSG00000136826" "9314"
When using cellranger or Seurat to quantify UMI-count (cf. Step1, Case III), the row names of the 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 <- convertRowID(input, rowID, LefttoRight))
##
|
| | 0%
|
|================== | 25%
|
|=================================== | 50%
|
|==================================================== | 75%
|
|======================================================================| 100%
## Input matrix: 4x5
## Output matrix: 4x5
## Some gene expression vectors were collapsed into single vector
## by sum rule
## $output
## [,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
##
## $ctable
## Left Right
## [1,] "POU5F1" "5460"
## [2,] "SOX2" "6657"
## [3,] "MYC" "4609"
## [4,] "KLF4" "9314"
# 1. Input matrix
input <- matrix(1:20, nrow=4, ncol=5)
# 3. Corresponding table
ah <- AnnotationHub()
# 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 <- convertRowID(input, rowID, LefttoRight))
##
|
| | 0%
|
|================== | 25%
|
|=================================== | 50%
|
|==================================================== | 75%
|
|======================================================================| 100%
## Input matrix: 4x5
## Output matrix: 4x5
## Some gene expression vectors were collapsed into single vector
## by sum rule
## $output
## [,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
##
## $ctable
## Left Right
## [1,] "POU5F1" "5460"
## [2,] "SOX2" "6657"
## [3,] "MYC" "4609"
## [4,] "KLF4" "9314"
Finally, we introduce some situations to perform some normalization methods of gene expression matrix.
If a user converts a Seurat object to a SingleCellExperient object by using as.SingleCellExperiment
,
the result of the NormalizeData
function (log counts) 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).
if(!require(scater)){
BiocManager::install("scater")
library(scater)
}
sce <- SingleCellExperiment(assays=list(counts = input))
cpm(sce) <- calculateCPM(sce)
sce <- normalize(sce)
assayNames(sce) # counts, normcounts, logcounts, cpm
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.
## R Under development (unstable) (2024-10-21 r87258)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.21-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.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] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] scTGIF_1.21.0
## [2] Homo.sapiens_1.3.1
## [3] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [4] org.Hs.eg.db_3.20.0
## [5] GO.db_3.20.0
## [6] OrganismDbi_1.49.0
## [7] GenomicFeatures_1.59.0
## [8] AnnotationDbi_1.69.0
## [9] SingleCellExperiment_1.29.0
## [10] SummarizedExperiment_1.37.0
## [11] Biobase_2.67.0
## [12] GenomicRanges_1.59.0
## [13] GenomeInfoDb_1.43.0
## [14] IRanges_2.41.0
## [15] S4Vectors_0.45.0
## [16] MatrixGenerics_1.19.0
## [17] matrixStats_1.4.1
## [18] scTensor_2.17.0
## [19] RSQLite_2.3.7
## [20] LRBaseDbi_2.17.0
## [21] AnnotationHub_3.15.0
## [22] BiocFileCache_2.15.0
## [23] dbplyr_2.5.0
## [24] BiocGenerics_0.53.0
## [25] BiocStyle_2.35.0
##
## loaded via a namespace (and not attached):
## [1] fs_1.6.4 bitops_1.0-9 enrichplot_1.27.0
## [4] httr_1.4.7 webshot_0.5.5 RColorBrewer_1.1-3
## [7] Rgraphviz_2.51.0 tools_4.5.0 backports_1.5.0
## [10] utf8_1.2.4 R6_2.5.1 lazyeval_0.2.2
## [13] withr_3.0.2 prettyunits_1.2.0 graphite_1.53.0
## [16] gridExtra_2.3 schex_1.21.0 fdrtool_1.2.18
## [19] cli_3.6.3 TSP_1.2-4 entropy_1.3.1
## [22] sass_0.4.9 genefilter_1.89.0 meshr_2.13.0
## [25] Rsamtools_2.23.0 yulab.utils_0.1.7 txdbmaker_1.3.0
## [28] gson_0.1.0 DOSE_4.1.0 R.utils_2.12.3
## [31] MeSHDbi_1.43.0 AnnotationForge_1.49.0 nnTensor_1.3.0
## [34] plotrix_3.8-4 maps_3.4.2 visNetwork_2.1.2
## [37] generics_0.1.3 gridGraphics_0.5-1 GOstats_2.73.0
## [40] BiocIO_1.17.0 dplyr_1.1.4 dendextend_1.18.1
## [43] Matrix_1.7-1 fansi_1.0.6 abind_1.4-8
## [46] R.methodsS3_1.8.2 lifecycle_1.0.4 yaml_2.3.10
## [49] qvalue_2.39.0 SparseArray_1.7.0 grid_4.5.0
## [52] blob_1.2.4 misc3d_0.9-1 crayon_1.5.3
## [55] ggtangle_0.0.3 lattice_0.22-6 msigdbr_7.5.1
## [58] cowplot_1.1.3 annotate_1.85.0 KEGGREST_1.47.0
## [61] magick_2.8.5 pillar_1.9.0 knitr_1.48
## [64] fgsea_1.33.0 tcltk_4.5.0 rjson_0.2.23
## [67] codetools_0.2-20 fastmatch_1.1-4 glue_1.8.0
## [70] outliers_0.15 ggfun_0.1.7 data.table_1.16.2
## [73] vctrs_0.6.5 png_0.1-8 treeio_1.31.0
## [76] spam_2.11-0 rTensor_1.4.8 gtable_0.3.6
## [79] assertthat_0.2.1 cachem_1.1.0 xfun_0.48
## [82] S4Arrays_1.7.0 mime_0.12 tidygraph_1.3.1
## [85] survival_3.7-0 seriation_1.5.6 iterators_1.0.14
## [88] tinytex_0.53 fields_16.3 nlme_3.1-166
## [91] Category_2.73.0 ggtree_3.15.0 bit64_4.5.2
## [94] progress_1.2.3 filelock_1.0.3 bslib_0.8.0
## [97] colorspace_2.1-1 DBI_1.2.3 tidyselect_1.2.1
## [100] bit_4.5.0 compiler_4.5.0 curl_5.2.3
## [103] httr2_1.0.5 graph_1.85.0 xml2_1.3.6
## [106] DelayedArray_0.33.0 plotly_4.10.4 bookdown_0.41
## [109] rtracklayer_1.67.0 checkmate_2.3.2 scales_1.3.0
## [112] hexbin_1.28.4 RBGL_1.83.0 plot3D_1.4.1
## [115] rappdirs_0.3.3 stringr_1.5.1 digest_0.6.37
## [118] rmarkdown_2.28 ca_0.71.1 XVector_0.47.0
## [121] htmltools_0.5.8.1 pkgconfig_2.0.3 highr_0.11
## [124] fastmap_1.2.0 rlang_1.1.4 htmlwidgets_1.6.4
## [127] UCSC.utils_1.3.0 farver_2.1.2 jquerylib_0.1.4
## [130] jsonlite_1.8.9 BiocParallel_1.41.0 GOSemSim_2.33.0
## [133] R.oo_1.26.0 RCurl_1.98-1.16 magrittr_2.0.3
## [136] GenomeInfoDbData_1.2.13 ggplotify_0.1.2 dotCall64_1.2
## [139] patchwork_1.3.0 munsell_0.5.1 Rcpp_1.0.13
## [142] babelgene_22.9 ape_5.8 viridis_0.6.5
## [145] stringi_1.8.4 tagcloud_0.6 ggraph_2.2.1
## [148] zlibbioc_1.53.0 MASS_7.3-61 plyr_1.8.9
## [151] parallel_4.5.0 ggrepel_0.9.6 Biostrings_2.75.0
## [154] graphlayouts_1.2.0 splines_4.5.0 hms_1.1.3
## [157] igraph_2.1.1 biomaRt_2.63.0 reshape2_1.4.4
## [160] BiocVersion_3.21.1 XML_3.99-0.17 evaluate_1.0.1
## [163] BiocManager_1.30.25 foreach_1.5.2 tweenr_2.0.3
## [166] tidyr_1.3.1 purrr_1.0.2 polyclip_1.10-7
## [169] heatmaply_1.5.0 ggplot2_3.5.1 ReactomePA_1.51.0
## [172] ggforce_0.4.2 xtable_1.8-4 restfulr_0.0.15
## [175] reactome.db_1.89.0 tidytree_0.4.6 viridisLite_0.4.2
## [178] tibble_3.2.1 aplot_0.2.3 ccTensor_1.0.2
## [181] GenomicAlignments_1.43.0 memoise_2.0.1 registry_0.5-1
## [184] cluster_2.1.6 concaveman_1.1.0 GSEABase_1.69.0