Using ctgGEM
2021-09-05
Contents
1 Introduction
Single-cell gene expression profiling methods are well suited for determining transcriptional heterogeneity across a cell population. In many cases, this heterogeneity may reflect dynamic changes in the transcriptional landscape due to biological processes such as differentiation or tumorigenesis. The overall trajectory of this change may be unified or branched and is partially obscured by the noise of the transcriptional programs within each cell. To address this computational challenge, a number of packages have been created to calculate pseudotemporal ordering of cells to reflect cell-state hierarchies from gene expression data. Each of these packages have unique inputs, outputs, and visualization schemes, which makes comparisons across these methods time-consuming. Here we present Cell Tree Generator for Gene Expression Matrices (ctgGEM), a package to streamline the building of cell-state hierarchies from single-cell gene expression data across multiple existing tools for improved comparability and reproducibility. ctgGEM provides the user with simplified way to build trees with these packages using one function call with a single dataset and the desired visualization name as a parameter. Results are also stored in the SIF file format for use in downstream analysis workflows or input into Cytoscape. The packages currently supported by ctgGEM are:
2 Installing ctgGEM
ctgGEM requires R version 4.0 or higher (available here) and the most recent version of Bioconductor. For more information on using Bioconductor, please see their website at https:// bioconductor.org. The following code will install Bioconductor, ctgGEM and the following CRAN packages and their dependencies: ggm, ggplot2, igraph, irlba, maptpx, VGAM
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ctgGEM")After installing, attach the ctgGEM package with the following command:
library(ctgGEM)3 Preparing a ctgGEMset object
The ctgGEM workflow is based around a single data class,
ctgGEMset, that extends the Biobase SummarizedExperiment
class, which provides a common interface familiar to those who have analyzed
microarray experiments with Bioconductor. The ctgGEMset class requires the
following three inputs:
exprsData, a numeric matrix of expression values, where rows are genes, and columns are cells/samplesphenoData, a data frame, where rows are cells/samples, and columns are cell/sample attributes (cell type, culture condition, day captured, etc.)featureData, a data frame, where rows are features (e.g. genes), and columns are gene attributes (gene identifiers, gc content, etc.)
The expression value matrix must have the same number of columns as the
phenoData has rows, and it must have the same number of rows as the featureData
data frame has rows. Row names of the phenoData object must match the column
names of the expression matrix. Row names of the featureData object must
matchrow names of the expression matrix. Details for loading the data for and
constructing an example ctgGEMset object suitable for this vignette can be
found in the following section.
3.1 Loading the Data
A ctgGEMset object that will support all of the tree types in
ctgGEM requires a gene expression matrix, which for the monocle
and TSCAN tree types must contain strictly non-negative values.
For this vignette, we will construct a toy object, using
the data provided in the HSMMSingleCell package. For more
information on this data set, please see the documentation available at
https://doi.org/doi:10.18129/B9.bioc.HSMMSingleCell
# load HSMMSingleCell package
library(HSMMSingleCell)
# load the data
data(HSMM_expr_matrix)
data(HSMM_sample_sheet)
data(HSMM_gene_annotation)3.2 Constructing a ctgGEMset object
Using the data loaded in the previous step, we can construct the ctgGEMset
object for this vignette as follows:
toyGEMset <- ctgGEMset(exprsData = HSMM_expr_matrix,
phenoData = HSMM_sample_sheet,
featureData = HSMM_gene_annotation)3.2.1 Additional Info for the “monocle” method
Using the “monocle” method requires a pair of parameters prepared with a column
name of gene identifiers in featureData() that corresponds to the gene short
names as the first item to be set, followed by the data type. The data type
must be one of “UMI”, “TC”, “FPKM”, “TPM”, “LTFPKM”, or “LTTPM”, where “UMI” is
UMI counts, “TC” is transcript counts, “FPKM” is FPKM, and “TPM” is TPM.
Monocle works best with untransformed data, but if you want to use
log-transformed FPKM or TPM, use data type “LTFPKM” or “LTTPM”, respectively.
Here we will use the gene_short_name column and the data type “FPKM”.
monocleInfo(toyGEMset, "gene_id") <- "gene_short_name"
monocleInfo(toyGEMset, "ex_type") <- "FPKM"The “monocle” method can be used in either semi-supervised or unsupervised mode. To run in semi-supervised mode, the names of two known classifying marker genes from the column specified as “gene_id” need to be set. If these two marker genes are not set, unsupervised mode will run. The two genes we will use for this demonstration are MYF5, a known marker for myoblasts, and ANPEP, a known marker for fibroblasts.
# Set two marker genes to use semi-supervised mode
# Alternatively, omit the next two lines to run in unsupervised mode
monocleInfo(toyGEMset, "cell_id_1") <- "MYF5" # marks myoblasts
monocleInfo(toyGEMset, "cell_id_2") <- "ANPEP" # marks fibroblasts3.2.2 Additional Info for the “TSCAN” method
In addition to its primary plot, the “TSCAN” method can also generate a single gene vs. pseudotime plot. To generate this plot, we need to supply the rowname of a single gene row in exprs(), and store it in TSCANinfo. Here we will use the “ENSG00000000003.10” gene.
TSCANinfo(toyGEMset) <- "ENSG00000000003.10"3.2.3 Additional Info for the “sincell” method
The “sincell” method can be used with a variety of parameters to control the
distance method used, which type (if any) dimensionality reduction to be used,
and which clustering method to use. If no options are specified, PCA will be
applied, with KNN clustering. Additional details concerning these parameters
can be found in the sincell package documentation. To use a distance method
with no dimensionality reduction, set the “method” parameter using the
sincellInfo() function to one of the following: “euclidean”, “L1” (Manhattan
distance), “cosine”, “pearson”, “spearman”, or “MI” (Mutual Information).
sincellInfo(toyGEMset, "method") <- "pearson" To use dimensionality reduction with the “sincell” method, set the “method” parameter to “PCA” for Principal Component Analysis, “ICA” for Independent Component Analysis, “tSNE” for t-Distributed Stochastic Neighbor Embedding, “classical-MDS” for classical Multidimensional Scaling, or “nonmetric-MDS” for non-metric Multidimensional Scaling. If using “classical-MDS” or “nonmetric- MDS”, we can also select the distance method to use with the “MDS.distance” parameter. The options for “MDS.distance” are the same as those listed above for no dimensionality reduction, with the exception the that the “cosine” method cannot be selected for “MDS.distance”.
sincellInfo(toyGEMset, "method") <- "classical-MDS"
sincellInfo(toyGEMset, "MDS.distance") <- "spearman"The final optional parameter for the “sincell” method is “clust.method”.
Acceptable values are “max-distance”, “percent”, “knn”, “k-medoids”, “ward.D”,
“ward.D2”, “single”, “complete”, “average”, “mcquitty”, “median”, or
“centroid”. In this example, we will set it to “k-medoids”.
sincellInfo(toyGEMset, "clust.method") <- "k-medoids"4 Using ctgGEM
To use ctgGEM, call the generate_tree function with the
desired tree method and ctgGEMset object. The user can specify the desired output
directory using the ‘outputDir’ parameter, or default to the temporary directory
returned by tempdir.
4.1 Using the “monocle” method
To use our example ctgGEMset and the “monocle” method, we would type the
following:
# Normally monocle will generate an internal warning for a condition with
# length > 1. If your R environment has the _R_CHECK_LENGTH_1_CONDITION_ set
# to verbose and/or abort, this will raise an error
if(!any(c("abort","verbose") %in% Sys.getenv("_R_CHECK_LENGTH_1_CONDITION_"))){
toyGEMset <- generate_tree(dataSet = toyGEMset, treeType = "monocle")
}This stores the final trees in the originalTrees list within the ctgGEMset
object, a simplified igraph version of the tree in the
treeList within the ctgGEMset object, and if necessary creates a new
folder, called “CTG-Output”, that contains a folder called “SIFs” containing
the .SIF text file for the final tree, and a folder called “Plots” containing a
.png image of the following plot:
4.2 Using the “sincell” method
To use our example ctgGEMset and the “sincell” method, we would type the
following:
toyGEMset <- generate_tree(dataSet = toyGEMset, treeType = "sincell")This stores the final trees in the originalTrees list within the ctgGEMset
object, a simplified igraph version of the tree in the
treeList within the ctgGEMset object, and if necessary creates a new
folder, called “CTG-Output”, that contains a folder called “Plots” containing
.png images of the following plots:
4.3 Using the “TSCAN” method
To use our example ctgGEMset and the “TSCAN”" method, we would type the
following:
toyGEMset <- generate_tree(dataSet = toyGEMset, treeType = "TSCAN")This stores the final trees in the originalTrees list within the ctgGEMset
object, a simplified igraph version of the tree in the
treeList within the ctgGEMset object, and if necessary creates a new
folder, called “CTG-Output”, that contains a folder called “SIFs” containing
the .SIF text file for the final tree, and a folder called “Plots” containing
.png images of the following plots:
5 Replotting Trees
If at some point we wish to view the plot of a tree generated after it’s been
created, but don’t want to have to regenerate it and all its files, ctgGEM has a function named plotOriginalTree(), that will
reproduce a plot stored in a ctgGEMset object. To use this function, we must
know the name of the tree we wish to plot. We can view the names of the trees
in our toyGEMset object using the names() function.
names(originalTrees(toyGEMset))## [1] "sincellMST" "sincellSST" "sincellIMC" "TSCANclustering"
## [5] "TSCANsingleGene"Once we have the names, we can choose a tree to plot. Let’s plot the “sincellIMC” tree again.
plotOriginalTree(toyGEMset, "sincellIMC")
Using this function eliminates the need to regenerate the tree to view a plot that was already created, thereby saving time for trees that require extensive computations to generate.
To store your analysis session result for later use, you can use the .Rda format.
save(toyGEMset, file = "toyGEMset.Rda")6 Session Information
sessionInfo()## R version 4.1.1 (2021-08-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows Server x64 (build 17763)
##
## 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] splines stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] HSMMSingleCell_1.13.0 ctgGEM_1.5.2
## [3] SummarizedExperiment_1.23.4 GenomicRanges_1.45.0
## [5] GenomeInfoDb_1.29.8 IRanges_2.27.2
## [7] S4Vectors_0.31.3 MatrixGenerics_1.5.4
## [9] matrixStats_0.60.1 monocle_2.21.1
## [11] DDRTree_0.1.5 irlba_2.3.3
## [13] VGAM_1.1-5 ggplot2_3.3.5
## [15] Biobase_2.53.0 BiocGenerics_0.39.2
## [17] Matrix_1.3-4 BiocStyle_2.21.3
##
## loaded via a namespace (and not attached):
## [1] Rtsne_0.15 colorspace_2.0-2
## [3] ellipsis_0.3.2 mclust_5.4.7
## [5] XVector_0.33.0 proxy_0.4-26
## [7] farver_2.1.0 ggrepel_0.9.1
## [9] fansi_0.5.0 codetools_0.2-18
## [11] docopt_0.7.1 knitr_1.33.5
## [13] spam_2.7-0 jsonlite_1.7.2
## [15] entropy_1.3.0 cluster_2.1.2
## [17] pheatmap_1.0.12 shiny_1.6.0
## [19] BiocManager_1.30.16 compiler_4.1.1
## [21] assertthat_0.2.1 fastmap_1.1.0
## [23] limma_3.49.4 later_1.3.0
## [25] htmltools_0.5.2 tools_4.1.1
## [27] igraph_1.2.6 dotCall64_1.0-1
## [29] gtable_0.3.0 glue_1.4.2
## [31] GenomeInfoDbData_1.2.6 RANN_2.6.1
## [33] reshape2_1.4.4 dplyr_1.0.7
## [35] maps_3.3.0 Rcpp_1.0.7
## [37] slam_0.1-48 jquerylib_0.1.4
## [39] TSCAN_1.31.0 vctrs_0.3.8
## [41] nlme_3.1-152 iterators_1.0.13
## [43] xfun_0.25 stringr_1.4.0
## [45] mime_0.11 lifecycle_1.0.0
## [47] TrajectoryUtils_1.1.2 gtools_3.9.2
## [49] statmod_1.4.36 zlibbioc_1.39.0
## [51] MASS_7.3-54 scales_1.1.1
## [53] TSP_1.1-10 promises_1.2.0.1
## [55] parallel_4.1.1 RColorBrewer_1.1-2
## [57] fields_12.5 SingleCellExperiment_1.15.2
## [59] yaml_2.2.1 gridExtra_2.3
## [61] sass_0.4.0 fastICA_1.2-2
## [63] stringi_1.7.4 highr_0.9
## [65] foreach_1.5.1 sincell_1.25.0
## [67] caTools_1.18.2 densityClust_0.3
## [69] rlang_0.4.11 pkgconfig_2.0.3
## [71] bitops_1.0-7 qlcMatrix_0.9.7
## [73] evaluate_0.14 lattice_0.20-44
## [75] purrr_0.3.4 labeling_0.4.2
## [77] tidyselect_1.1.1 plyr_1.8.6
## [79] magrittr_2.0.1 bookdown_0.24
## [81] R6_2.5.1 magick_2.7.3
## [83] gplots_3.1.1 generics_0.1.0
## [85] combinat_0.0-8 DelayedArray_0.19.2
## [87] DBI_1.1.1 pillar_1.6.2
## [89] withr_2.4.2 mgcv_1.8-36
## [91] scatterplot3d_0.3-41 RCurl_1.98-1.4
## [93] tibble_3.1.4 crayon_1.4.1
## [95] KernSmooth_2.23-20 utf8_1.2.2
## [97] rmarkdown_2.10 viridis_0.6.1
## [99] grid_4.1.1 FNN_1.1.3
## [101] sparsesvd_0.2 digest_0.6.27
## [103] xtable_1.8-4 httpuv_1.6.2
## [105] munsell_0.5.0 viridisLite_0.4.0
## [107] bslib_0.3.07 References
7.1 monocle
Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, Lennon NJ, Livak KJ, Mikkelsen TS and Rinn JL (2014). “The dynamics and regulators of cell fate decisions are revealed by pseudo-temporal ordering of single cells.” Nature Biotechnology.
monocle package URL: https://bioconductor.org/packages/release/bioc/html/monocle.html
7.2 sincell
Juliá M, Telenti A, Rausell A (2014): Sincell: R package for the statistical assessment of cell state hierarchies from single-cell RNA-seq data. bioRxiv preprint
sincell package URL: https://bioconductor.org/packages/release/bioc/html/sincell.html
7.3 TSCAN
Zhicheng Ji and Hongkai Ji (2015). TSCAN: TSCAN: Tools for Single-Cell ANalysis. R package version 1.12.0.
TSCAN package URL: https://bioconductor.org/packages/release/bioc/html/TSCAN.html