scMaSigPro: Bridging maSigPro Analysis
Guide to applying maSigPro workflow with scMaSigPro
22 febrero 2024
Abstract
scMaSigPro is an R package designed for the serial analysis of single-cell
RNA-seq (scRNA-Seq) data along inferred pseudotime. It builds upon the maSigPro
Bioconductor package to identify genes with significant expression changes across different
branching paths in a pseudotime-ordered scRNA-Seq dataset. This vignette illustrates
the workflow of the original maSigPro Bioconductor package of using scMaSigPro.
Package
scMaSigPro 0.0.1
Contents
- 0.1 Introduction
- 0.2 Overview of
maSigPro - 0.3 Adaptation for Bulk-RNA Seq Data
- 0.4 Setup and Installation
- 0.5 Creating objects for
scMaSigPro - 0.6 Creatin of Regression Matrix
- 0.7 Detecting Non-Flat Profiles with
maSigPro::p.vector() - 0.8 Detecting Significant Polynomial Terms with
maSigPro::T.fit() - 0.9 R-Square Filter with
maSigPro::get.siggenes() - 0.10 Visualizing the trends and Intersection
- 0.11 Session Info
0.1 Introduction
scMaSigPro is a polynomial regression-based approach inspired by the maSigPro
Bioconductor package tailored for scRNA-Seq data. It first discretizes single cell
expression along inferred pseudotime while preserving order. Afterwards it applies
the maSigPro model
to pinpoint genes exhibiting significant expression profile differences among
branching paths and pseudotime.
0.2 Overview of maSigPro
maSigPro (microarray significant profiles) is a bioconductor package. Originally designed for microarray data analysis, maSigPro focuses on identifying genes with significant expression profiles that vary among different experimental groups, particularly in time-course experiments.
0.3 Adaptation for Bulk-RNA Seq Data
In 2016, the maSigPro package was adapted for bulk-RNA seq data, a crucial development given the distinct differences between microarray and RNA-Seq technologies. Unlike the distributions typically seen in microarray data, bulk-RNA Seq data generally follows a Negative Binomial Distribution, which is more apt for count data like RNA-Seq, characterized by variance often exceeding the mean.
This adaptation brought in explicit theta parameterization, critical
for RNA-Seq data analysis as the theta parameter represents the dispersion in
the Negative Binomial Distribution. Such parameterization enhances the accuracy
in modeling gene expression variability, a key factor in identifying significant
changes in gene expression, thereby marking a significant step forward in the
application of maSigPro to the more complex RNA-Seq data.
0.4 Setup and Installation
0.4.1 Installing scMaSigPro
scMaSigPro is currently available on GitHub. Follow these steps to install it:
# Install devtools if not already installed
if (!requireNamespace("devtools", quietly = TRUE)) {
install.packages("devtools")
}
# Install scMaSigPro
devtools::install_github("BioBam/scMaSigPro",
ref = "main",
build_vignettes = FALSE,
build_manual = TRUE,
upgrade = "never",
force = TRUE,
quiet = TRUE,
)0.4.2 Installing maSigPro
You can easily install maSigPro from Bioconductor using the following steps:
# Install Bioconductor if not already installed
if (!require("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
# install maSigPro
BiocManager::install("maSigPro")0.4.3 Setup
First, we need to load both scMaSigPro and maSigPro. Additionally, we will load
testthat, a user-friendly testing framework for R that integrates well with
existing workflows. We will use testthat to verify if scMaSigPro produces
the same results as maSigPro when used in bulk mode.
# load scMaSigPro
library(scMaSigPro)
# load maSigPro
library(maSigPro)
# We will also keep checking the output
library(testthat)0.4.4 Setting up the data
In this tutorial, we will use datasets from maSigPro and aim to replicate the
same analysis pipeline using scMaSigPro. For more detailed information about
the data analysis, refer to the original maSigPro vignette-userguide
# load counts
data(data.abiotic)
# load Design file
data(edesign.abiotic)
# save with different names for scMaSigPro
count <- data.abiotic
metadata <- as.data.frame(edesign.abiotic)Let’s briefly examine each dataset:
# View a subset of count data
print(data.abiotic[c(1:5), c(1:5)])## Control_3H_1 Control_3H_2 Control_3H_3 Control_9H_1 Control_9H_2
## STMDF90 0.13735714 -0.3653065 -0.15329448 0.44754535 0.287476796
## STMCJ24 NA NA NA NA NA
## STMJH42 0.07864449 0.1002328 -0.17365488 -0.25279484 0.184855409
## STMDE66 0.22976991 0.4740975 0.46930716 0.37101059 -0.004992029
## STMIX74 0.14407618 -0.4801864 -0.07847999 0.05692331 0.013045420# View a subset of the design file
print(edesign.abiotic[c(1:5), c(1:5)])## Time Replicate Control Cold Heat
## Control_3H_1 3 1 1 0 0
## Control_3H_2 3 1 1 0 0
## Control_3H_3 3 1 1 0 0
## Control_9H_1 9 2 1 0 0
## Control_9H_2 9 2 1 0 0In edesign.abiotic, there are four groups, including one control group. The experiment features experimental capture times, which are distinct from Pseudotime, and includes temporal replicates.
To simplify the analysis, we’ll convert the binarized groups into a single column.
scMaSigPro automatically binarizes groups internally, reducing the likelihood of errors.
# Add group column
metadata$Group <- apply(metadata[, c(3:6)], 1, FUN = function(x) {
return(names(x[x == 1]))
})
# Remove binary columns
metadata <- metadata[, !(colnames(metadata) %in%
c("Control", "Cold", "Heat", "Salt")),
drop = FALSE
]
# View
print(metadata[c(1:5), ])## Time Replicate Group
## Control_3H_1 3 1 Control
## Control_3H_2 3 1 Control
## Control_3H_3 3 1 Control
## Control_9H_1 9 2 Control
## Control_9H_2 9 2 Control0.5 Creating objects for scMaSigPro
scMaSigPro utilizes S4 classes, and therefore, the initial step involves creating
an object of the scMaSigPro class, also known as scmpObject. To do this, we use
the create_scmp() function. We will use the same object throughout the tutorial
to show the concept of S4 classes. On the other hand for maSigPro we will be making
different S3 objects as shown in the `original tutorial.
As we are using scMaSigPro in bulk mode, we will set use_as_bin = TRUE. This
directs the data to the dense slot of the scmpObject, allowing us to proceed with
the analysis without using sc.squeeze(), which is specifically designed for scRNA-Seq
pseudo-bulking.
scmp_ob <- create_scmp(
counts = count,
cell_data = metadata,
ptime_col = "Time",
path_col = "Group",
use_as_bin = TRUE
)
# Print Console output
scmp_ob## Class: ScMaSigProClass
## nCells: 36
## nFeatures: 1000
## Pseudotime Range: 3 27
## Branching Paths: Control, Cold, Heat, Salt
## Binned Pseudotime: 3-27(Range), 13(Mean),
## Number of bins-> Cold: 9 Control: 9 Heat: 9 Salt: 9
## Average bin Size-> Cold: 1 Control: 1 Heat: 1 Salt: 1In the console output, we can observe that the data has been transferred to the
dense slot. We can visualize this using plotBinTile().
plotBinTile(scmp_ob)
Next, we need to verify if scMaSigPro has created the appropriate count matrix.
For this, we will use testthat::expect_identical(). This function highlights any
mismatches; if everything is identical, it will not throw an error.
The eDense() function retrieves the dense expression file from the scMaSigPro object.
expect_identical(eDense(scmp_ob), expected = data.abiotic)No errors were thrown, indicating that our setup is correct and we can proceed with the analysis.
0.6 Creatin of Regression Matrix
The function maSigPro::make.design.matrix() is used to create a regression matrix
for fitting the full regression model. In scMaSigPro, this task is accomplished
using sc.set.poly(), where the regression matrix is referred to as “predictor_matrix”.
We will also verify that scMaSigPro has correctly performed binarization, as
specified in edesign.abiotic.
0.6.1 maSigPro::make.design.matrix() and scMaSigPro::sc.set.poly()
# Using maSigPro
design <- make.design.matrix(edesign.abiotic, degree = 2)
# Using scMaSigPro
scmp_ob <- sc.set.poly(scmp_ob, poly_degree = 2)
# Comparing binarization results
expect_equal(pathAssign(scmp_ob), expected = edesign.abiotic)
# Comparing regression matrices
expect_identical(scmp_ob@Design@predictor_matrix, expected = design$dis)No errors were thrown, so we can proceed with the analysis.
0.7 Detecting Non-Flat Profiles with maSigPro::p.vector()
To compute a regression fit for each gene, maSigPro uses p.vector(), while scMaSigPro
uses sc.p.vector(). The main difference is the offset parameter in scMaSigPro,
which accounts for variations in bin sizes. Since we are not dealing with
scRNA-Seq data, we will disable this parameter.
# Using maSigPro
gc <- capture_output(fit <- p.vector(data.abiotic, design,
Q = 0.05,
MT.adjust = "BH", min.obs = 20
))
# Using ScMaSigPro
scmp_ob <- sc.p.vector(scmp_ob,
min_na = 20,
verbose = FALSE,
link = "identity",
offset = FALSE,
max_it = 25,
epsilon = 0.00001,
family = gaussian()
)Now, we’ll compare the S3 object fit frommaSigPro::p.vector()with the results fromsc.p.vector()`:
# Compare p-values
expect_identical(
matrix(scmp_ob@Profile@p_values,
dimnames = list(names(scmp_ob@Profile@p_values), "p.value")
),
expected = as.matrix(fit$p.vector[, 1, drop = FALSE])
)
# Compare adjusted p-values
pad <- scmp_ob@Profile@adj_p_values
names(pad) <- NULL
expect_identical(pad, expected = fit$p.adjusted)0.8 Detecting Significant Polynomial Terms with maSigPro::T.fit()
Next, we apply a variable selection procedure to identify significant polynomial
terms for each gene, using T.fit() and sc.t.fit().
# Using maSigPro
gc <- capture_output(tstep <- T.fit(fit, step.method = "backward", alfa = 0.05))
# Using scMaSigPro
scmp_ob <- sc.t.fit(scmp_ob,
offset = FALSE,
verbose = FALSE,
epsilon = 0.00001,
family = gaussian()
)We will compare the s3 object tstep from maSigPro::T.fit() with the results from sc.t.fit():
# Solutions
expect_identical(showSol(scmp_ob), expected = tstep$sol)
# Coefficients
expect_identical(showCoeff(scmp_ob), expected = tstep$coefficients)
# Group Coefficients
expect_identical(showGroupCoeff(scmp_ob), expected = tstep$group.coeffs)
# tscore
expect_identical(as.data.frame(showTS(scmp_ob)), expected = tstep$t.score)0.9 R-Square Filter with maSigPro::get.siggenes()
The next step is to generate lists of significant genes based on desired biological question.
This is achieved using the get.siggenes() function in maSigPro and sc.filter() in scMaSigPro.
# maSigPro
sigs <- get.siggenes(tstep, rsq = 0.6, vars = "groups")
# scMaSigPro
scmp_ob <- sc.filter(scmp_ob, rsq = 0.6, vars = "groups")We’ll compare the s3 object (sigs) from `maSigPro::get.siggenes()’ with the results from sc.filter()’:
# Compare Cold vs Control
expect_identical(scmp_ob@Significant@genes$ColdvsControl,
expected = sigs$summary$ColdvsControl
)
# Compare Heat vs Control
expect_identical(scmp_ob@Significant@genes$HeatvsControl,
expected = sigs$summary$HeatvsControl[sigs$summary$HeatvsControl != " "]
)
# Compare Salt vs Control
expect_identical(scmp_ob@Significant@genes$SaltvsControl,
expected = sigs$summary$SaltvsControl[sigs$summary$SaltvsControl != " "]
)
# Compare Control
expect_identical(scmp_ob@Significant@genes$Control,
expected = sigs$summary$Control[sigs$summary$Control != " "]
)All comparisons indicate equality, confirming that the results are reproducible up to this point.
0.10 Visualizing the trends and Intersection
Having confirmed that the numerical values are identical for both methods, our next step is to assess whether the actual trends are also consistent.
0.10.1 Visualizing Intersection maSigPro::suma2Venn()
# Venn Diagram of maSigPro
suma2Venn(sigs$summary[, c(1:4)])
# Upset plot of scMaSigPro
plotIntersect(
scmpObj = scmp_ob,
min_intersection_size = 0
)
0.10.2 Visualizing Single trend with maSigPro::PlotGroups()
The PlotGroups() function in maSigPro generates a plot of gene expression profiles
with time as the x-axis and gene expression on the y-axis. In these plots, gene
expression from the same experimental group is represented in the same color, and
lines are drawn to join the averages of each time-group to visualize the trend of
each group over time. For scMaSigPro, this functionality is achieved using
plotTrend().
# Extracting gene "STMDE66" from data
STMDE66 <- data.abiotic[rownames(data.abiotic) == "STMDE66", ]
# Plotting with maSigPro
PlotGroups(STMDE66,
edesign = edesign.abiotic, show.fit = TRUE,
dis = design$dis, groups.vector = design$groups.vector
)
# Plotting the same gene with scMaSigPro
plotTrend(scmp_ob, "STMDE66",
logs = FALSE, pseudoCount = 0,
smoothness = 0.01, significant = FALSE,
summary_mode = "mean"
)
0.10.3 Visualizing Cluster trend with maSigPro::PlotProfiles()
# Plot clustered Trend
gc <- capture_output(
res <- see.genes(sigs$sig.genes$ColdvsControl,
show.fit = TRUE, dis = design$dis,
cluster.method = "hclust", cluster.data = 1, k = 9
)
)
# Compute Clusters
scmp_ob <- sc.cluster.trend(scmp_ob,
geneSet = "ColdvsControl",
cluster_by = "counts"
)0.10.4 Plot cluster Trends
plotTrendCluster(
scmpObj = scmp_ob,
plot = "coeff",
verbose = FALSE
)
0.11 Session Info
sessionInfo(package = "scMaSigPro")## R version 4.3.0 (2023-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=es_ES.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=es_ES.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=es_ES.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Madrid
## tzcode source: system (glibc)
##
## attached base packages:
## character(0)
##
## other attached packages:
## [1] scMaSigPro_0.0.1
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-7 remotes_2.4.2.1
## [3] testthat_3.2.1 rlang_1.1.3
## [5] magrittr_2.0.3 matrixStats_1.2.0
## [7] e1071_1.7-14 compiler_4.3.0
## [9] mgcv_1.9-1 roxygen2_7.2.3
## [11] venn_1.11 vctrs_0.6.5
## [13] stringr_1.5.1 profvis_0.3.8
## [15] pkgconfig_2.0.3 crayon_1.5.2
## [17] fastmap_1.1.1 magick_2.8.2
## [19] XVector_0.42.0 ellipsis_0.3.2
## [21] labeling_0.4.3 utf8_1.2.4
## [23] promises_1.2.1 rmarkdown_2.25
## [25] grDevices_4.3.0 sessioninfo_1.2.2
## [27] waldo_0.5.2 purrr_1.0.2
## [29] xfun_0.42 zlibbioc_1.48.0
## [31] cachem_1.0.8 graphics_4.3.0
## [33] GenomeInfoDb_1.38.6 jsonlite_1.8.8
## [35] highr_0.10 later_1.3.2
## [37] DelayedArray_0.28.0 parallel_4.3.0
## [39] R6_2.5.1 bslib_0.6.1
## [41] stringi_1.8.3 ComplexUpset_1.3.3
## [43] parallelly_1.37.0 pkgload_1.3.4
## [45] brio_1.1.4 jquerylib_0.1.4
## [47] GenomicRanges_1.54.1 bookdown_0.37
## [49] Rcpp_1.0.12 assertthat_0.2.1
## [51] SummarizedExperiment_1.32.0 knitr_1.45
## [53] usethis_2.2.2 IRanges_2.36.0
## [55] splines_4.3.0 httpuv_1.6.14
## [57] Matrix_1.6-5 igraph_1.6.0
## [59] tidyselect_1.2.0 yaml_2.3.8
## [61] rstudioapi_0.15.0 abind_1.4-5
## [63] miniUI_0.1.1.1 admisc_0.34
## [65] pkgbuild_1.4.3 lattice_0.22-5
## [67] tibble_3.2.1 Biobase_2.62.0
## [69] shiny_1.8.0 withr_3.0.0
## [71] evaluate_0.23 base_4.3.0
## [73] desc_1.4.3 proxy_0.4-27
## [75] urlchecker_1.0.1 xml2_1.3.6
## [77] mclust_6.0.1 BiocManager_1.30.22
## [79] pillar_1.9.0 MatrixGenerics_1.14.0
## [81] stats4_4.3.0 RColorConesa_1.0.0
## [83] plotly_4.10.4 generics_0.1.3
## [85] rprojroot_2.0.4 RCurl_1.98-1.14
## [87] S4Vectors_0.40.2 ggplot2_3.4.4
## [89] munsell_0.5.0 scales_1.3.0
## [91] BiocStyle_2.30.0 stats_4.3.0
## [93] xtable_1.8-4 class_7.3-22
## [95] glue_1.7.0 lazyeval_0.2.2
## [97] tools_4.3.0 datasets_4.3.0
## [99] data.table_1.15.0 fs_1.6.3
## [101] grid_4.3.0 utils_4.3.0
## [103] tidyr_1.3.1 methods_4.3.0
## [105] devtools_2.4.5 colorspace_2.1-0
## [107] SingleCellExperiment_1.24.0 nlme_3.1-164
## [109] patchwork_1.1.3 GenomeInfoDbData_1.2.11
## [111] cli_3.6.2 fansi_1.0.6
## [113] S4Arrays_1.2.0 viridisLite_0.4.2
## [115] dplyr_1.1.4 gtable_0.3.4
## [117] sass_0.4.8 digest_0.6.34
## [119] BiocGenerics_0.48.1 SparseArray_1.2.4
## [121] farver_2.1.1 htmlwidgets_1.6.4
## [123] memoise_2.0.1 maSigPro_1.74.0
## [125] entropy_1.3.1 htmltools_0.5.7
## [127] lifecycle_1.0.4 httr_1.4.7
## [129] mime_0.12 MASS_7.3-60sessionInfo(package = "maSigPro")## R version 4.3.0 (2023-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=es_ES.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=es_ES.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=es_ES.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Madrid
## tzcode source: system (glibc)
##
## attached base packages:
## character(0)
##
## other attached packages:
## [1] maSigPro_1.74.0
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-7 remotes_2.4.2.1
## [3] testthat_3.2.1 rlang_1.1.3
## [5] magrittr_2.0.3 matrixStats_1.2.0
## [7] e1071_1.7-14 compiler_4.3.0
## [9] mgcv_1.9-1 roxygen2_7.2.3
## [11] venn_1.11 vctrs_0.6.5
## [13] stringr_1.5.1 profvis_0.3.8
## [15] pkgconfig_2.0.3 crayon_1.5.2
## [17] fastmap_1.1.1 magick_2.8.2
## [19] XVector_0.42.0 ellipsis_0.3.2
## [21] labeling_0.4.3 utf8_1.2.4
## [23] promises_1.2.1 rmarkdown_2.25
## [25] grDevices_4.3.0 sessioninfo_1.2.2
## [27] waldo_0.5.2 purrr_1.0.2
## [29] xfun_0.42 zlibbioc_1.48.0
## [31] cachem_1.0.8 graphics_4.3.0
## [33] GenomeInfoDb_1.38.6 jsonlite_1.8.8
## [35] highr_0.10 later_1.3.2
## [37] DelayedArray_0.28.0 parallel_4.3.0
## [39] R6_2.5.1 bslib_0.6.1
## [41] stringi_1.8.3 ComplexUpset_1.3.3
## [43] parallelly_1.37.0 pkgload_1.3.4
## [45] brio_1.1.4 jquerylib_0.1.4
## [47] GenomicRanges_1.54.1 bookdown_0.37
## [49] Rcpp_1.0.12 scMaSigPro_0.0.1
## [51] assertthat_0.2.1 SummarizedExperiment_1.32.0
## [53] knitr_1.45 usethis_2.2.2
## [55] IRanges_2.36.0 splines_4.3.0
## [57] httpuv_1.6.14 Matrix_1.6-5
## [59] igraph_1.6.0 tidyselect_1.2.0
## [61] yaml_2.3.8 rstudioapi_0.15.0
## [63] abind_1.4-5 miniUI_0.1.1.1
## [65] admisc_0.34 pkgbuild_1.4.3
## [67] lattice_0.22-5 tibble_3.2.1
## [69] Biobase_2.62.0 shiny_1.8.0
## [71] withr_3.0.0 evaluate_0.23
## [73] base_4.3.0 desc_1.4.3
## [75] proxy_0.4-27 urlchecker_1.0.1
## [77] xml2_1.3.6 mclust_6.0.1
## [79] BiocManager_1.30.22 pillar_1.9.0
## [81] MatrixGenerics_1.14.0 stats4_4.3.0
## [83] RColorConesa_1.0.0 plotly_4.10.4
## [85] generics_0.1.3 rprojroot_2.0.4
## [87] RCurl_1.98-1.14 S4Vectors_0.40.2
## [89] ggplot2_3.4.4 munsell_0.5.0
## [91] scales_1.3.0 BiocStyle_2.30.0
## [93] stats_4.3.0 xtable_1.8-4
## [95] class_7.3-22 glue_1.7.0
## [97] lazyeval_0.2.2 tools_4.3.0
## [99] datasets_4.3.0 data.table_1.15.0
## [101] fs_1.6.3 grid_4.3.0
## [103] utils_4.3.0 tidyr_1.3.1
## [105] methods_4.3.0 devtools_2.4.5
## [107] colorspace_2.1-0 SingleCellExperiment_1.24.0
## [109] nlme_3.1-164 patchwork_1.1.3
## [111] GenomeInfoDbData_1.2.11 cli_3.6.2
## [113] fansi_1.0.6 S4Arrays_1.2.0
## [115] viridisLite_0.4.2 dplyr_1.1.4
## [117] gtable_0.3.4 sass_0.4.8
## [119] digest_0.6.34 BiocGenerics_0.48.1
## [121] SparseArray_1.2.4 farver_2.1.1
## [123] htmlwidgets_1.6.4 memoise_2.0.1
## [125] entropy_1.3.1 htmltools_0.5.7
## [127] lifecycle_1.0.4 httr_1.4.7
## [129] mime_0.12 MASS_7.3-60