simPIC: simulating single-cell ATAC-seq data
2024-10-29
Package
simPIC 1.3.0
1 Introduction to simPIC
simPIC is an R package for the simple simulation of single-cell RNA sequencing data. This vignette gives an overview and introduction to simPIC’s functionality.
2 Installation
simPIC can be installed from Bioconductor:
BiocManager::install("simPIC")3 Quick start
Assuming you already have a matrix of count data similar to that you wish to simulate there are two simple steps to creating a simulated data set with simPIC.
# Load package
suppressPackageStartupMessages({
library(simPIC)
})
# Load test data
set.seed(12)
counts <- readRDS(system.file("extdata", "test.rds", package = "simPIC"))
# Estimate parameters
est <- simPICestimate(counts)
#> simPIC is:
#> estimating library size parameters...
#> estimating sparsity...
#> estimating peak mean parameters...
#> using weibull distribution for estimating peak mean
# Simulate data using estimated parameters
sim <- simPICsimulate(est)
#> simPIC is:
#> updating parameters...
#> setting up SingleCellExperiment object...
#> Simulating library size...
#> Simulating peak mean...
#> using weibull distribution for simulating peak mean
#> Simulating true counts...
#> Done!!These steps will be explained in detail in the following sections but briefly the first step takes a dataset and estimates simulation parameters from it and the second step takes those parameters and simulates a new dataset.
4 Input data
simPIC recommends to use a Paired-Insertion Count (PIC) matrix for optimal use
of the quantitative information present in scATAC-seq data. You can convert your
own matrix to PIC by following the PICsnATAC
vignette. Briefly, you need three input files for PIC -
- cell barcodes with metadata (
singlecell.csv) - list of peak regions (
peaks.bed) - fragment files (
fragment.tsv.gz)
pic_mat <- PIC_counting(
cells = cells,
fragment_tsv_gz_file_location = fragment_tsv_gz_file_location,
peak_sets = peak_sets
)5 The simPIC simulation
The core of the simPIC model is a gamma-Poisson distribution. This model is used
to generate a uniform quantification based peak-by-cell count matrix. Mean
chromatin accessibility signals for each peak are simulated from a weibull
distribution with default
settings. Users have the flexibility to choose from gamma, lognormal or
pareto distribution as well. Each cell is given an expected library size,
simulated from a log-normal distribution to match to a given dataset. Sparsity
is imposed on counts simulated from a Poisson
distribution.
6 The simPICcount class
All the parameters for the simPIC simulation are stored in a simPICcount
object. A class specifically desgined for storing simPIC scATAC-seq simulation
parameters. Let’s create a new one and see what it looks like.
sim.params <- newsimPICcount()we can see the default values for the simPICcount object parameters. These
values are based on provided test data.
sim.params
#> An object of class "simPICcount"
#> Slot "nPeaks":
#> [1] 5000
#>
#> Slot "nCells":
#> [1] 700
#>
#> Slot "seed":
#> [1] 96612
#>
#> Slot "default":
#> [1] TRUE
#>
#> Slot "pm.distr":
#> [1] "weibull"
#>
#> Slot "lib.size.meanlog":
#> [1] 6.687082
#>
#> Slot "lib.size.sdlog":
#> [1] 0.344361
#>
#> Slot "peak.mean.shape":
#> [1] 0.7909301
#>
#> Slot "peak.mean.rate":
#> [1] 7.100648
#>
#> Slot "peak.mean.scale":
#> [1] 0.09522228
#>
#> Slot "peak.mean.pi":
#> [1] -17.17441
#>
#> Slot "peak.mean.meanlog":
#> [1] -2.825233
#>
#> Slot "peak.mean.sdlog":
#> [1] -1.366378
#>
#> Slot "sparsity":
#> [1] 0.31647915 0.41055728 0.77536917 0.92773272 0.64524861 0.18688732
#> [7] 0.29059723 0.38707305 0.56802433 0.04876861 0.00083383 0.09877044
#> [13] 0.39601450 0.65099842 0.50916994 0.92962847 0.85701174 0.75142452
#> [19] 0.53441732 0.33270524 0.26845395 0.27826479 0.20090212 0.41940649
#> [25] 0.07230980 0.38349619 0.28234327 0.42958552 0.91643195 0.75386006
#> [ reached 'max' / getOption("max.print") -- omitted 199970 entries ]6.1 Getting and setting parameters
To get a particular parameter, for e.g., number of peaks, we can use simPICget
function:
simPICget(sim.params, "nPeaks")
#> [1] 5000Alternatively, to give a parameter a new value we can use the
setsimPICparameters function:
sim.params <- setsimPICparameters(sim.params, nPeaks = 2000)
simPICget(sim.params, "nPeaks")
#> [1] 2000To get or set multiple parameters use simPICget or setsimPICparameters
functions:
# Set multiple parameters at once (using a list)
sim.params <- setsimPICparameters(sim.params,
update = list(nPeaks = 8000, nCells = 500)
)
# Extract multiple parameters as a list
params <- simPICgetparameters(
sim.params,
c("nPeaks", "nCells", "peak.mean.shape")
)
# Set multiple parameters at once (using additional arguments)
params <- setsimPICparameters(sim.params,
lib.size.sdlog = 3.5, lib.size.meanlog = 9.07
)
params
#> An object of class "simPICcount"
#> Slot "nPeaks":
#> [1] 8000
#>
#> Slot "nCells":
#> [1] 500
#>
#> Slot "seed":
#> [1] 96612
#>
#> Slot "default":
#> [1] TRUE
#>
#> Slot "pm.distr":
#> [1] "weibull"
#>
#> Slot "lib.size.meanlog":
#> [1] 9.07
#>
#> Slot "lib.size.sdlog":
#> [1] 3.5
#>
#> Slot "peak.mean.shape":
#> [1] 0.7909301
#>
#> Slot "peak.mean.rate":
#> [1] 7.100648
#>
#> Slot "peak.mean.scale":
#> [1] 0.09522228
#>
#> Slot "peak.mean.pi":
#> [1] -17.17441
#>
#> Slot "peak.mean.meanlog":
#> [1] -2.825233
#>
#> Slot "peak.mean.sdlog":
#> [1] -1.366378
#>
#> Slot "sparsity":
#> [1] 0.31647915 0.41055728 0.77536917 0.92773272 0.64524861 0.18688732
#> [7] 0.29059723 0.38707305 0.56802433 0.04876861 0.00083383 0.09877044
#> [13] 0.39601450 0.65099842 0.50916994 0.92962847 0.85701174 0.75142452
#> [19] 0.53441732 0.33270524 0.26845395 0.27826479 0.20090212 0.41940649
#> [25] 0.07230980 0.38349619 0.28234327 0.42958552 0.91643195 0.75386006
#> [ reached 'max' / getOption("max.print") -- omitted 199970 entries ]7 Estimating parameters
simPIC allows you to estimate many of it’s parameters from a
SingleCellExperiment object containing counts or a counts matrix using the
simPICestimate function.
# Get the counts from test data
count <- readRDS(system.file("extdata", "test.rds", package = "simPIC"))
# Check that counts is a dgCMatrix
class(count)
#> [1] "dgCMatrix"
#> attr(,"package")
#> [1] "Matrix"
typeof(count)
#> [1] "S4"
# Check the dimensions, each row is a peak, each column is a cell
dim(count)
#> [1] 5000 700
# Show the first few entries
count[1:5, 1:5]
#> 5 x 5 sparse Matrix of class "dgCMatrix"
#> TTGCCCACACTCAAGT-1 TATGTGGGTCCAAGAG-1
#> chr22-29837913-29838812 1 2
#> chr12-90076794-90077665 . .
#> chr11-63574191-63575102 . .
#> chr12-54977628-54978479 . .
#> chr1-221715272-221716182 . .
#> CGATGATCAAGTAACA-1 GTCGTAATCAGGAAGC-1
#> chr22-29837913-29838812 . .
#> chr12-90076794-90077665 . .
#> chr11-63574191-63575102 . .
#> chr12-54977628-54978479 . .
#> chr1-221715272-221716182 . .
#> TACTAGGCACAAACAA-1
#> chr22-29837913-29838812 .
#> chr12-90076794-90077665 .
#> chr11-63574191-63575102 .
#> chr12-54977628-54978479 .
#> chr1-221715272-221716182 1
new <- newsimPICcount()
new <- simPICestimate(count)
#> simPIC is:
#> estimating library size parameters...
#> estimating sparsity...
#> estimating peak mean parameters...
#> using weibull distribution for estimating peak mean
## estimating using gamma distribution
## new <- simPICestimate(count, pm.distr = "gamma")Here we estimated parameters from a counts matrix using default parameters. The estimation process has the following steps:
- Mean parameters are estimated by fitting a weibull distribution (default) to the peak means.
- Library size parameters are estimated by fitting a log-normal distribution to the library sizes.
- Sparsity parameter is estimated by fitting a Bernoulli distribution.
For more details of the estimation procedures see ?simPICestimate.
8 Simulating counts
Once we have a set of parameters we are happy with we can use simPICsimulate
to simulate counts. To make adjustments to the parameters provide them as
additional arguments. Alternatively if we don’t supply any parameters defaults
will be used:
sim <- simPICsimulate(new, nCells = 1000)
#> simPIC is:
#> updating parameters...
#> setting up SingleCellExperiment object...
#> Simulating library size...
#> Simulating peak mean...
#> using weibull distribution for simulating peak mean
#> Simulating true counts...
#> Done!!
sim
#> class: SingleCellExperiment
#> dim: 5000 1000
#> metadata(1): Params
#> assays(1): counts
#> rownames(5000): Peak1 Peak2 ... Peak4999 Peak5000
#> rowData names(2): Peak exp.peakmean
#> colnames(1000): Cell1 Cell2 ... Cell999 Cell1000
#> colData names(2): Cell exp.libsize
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
## simulating using gamma distribution
## sim <- simPICsimulate(new, nCells =1000, pm.distr = "gamma")Looking at the output of simPICsimulate we can see that sim is
SingleCellExperiment object with 5000 features (peaks) and
1000 samples (cells). The main part of this object is a features by
samples matrix containing the simulated counts (accessed using counts).
Additionally a SingleCellExperiment contains information about each cell
(accessed using colData) and each peak (accessed using rowData). simPIC uses
these slots, as well as assays, to store information about the intermediate
values of the simulation.
# Access the counts
counts(sim)[1:5, 1:5]
#> Cell1 Cell2 Cell3 Cell4 Cell5
#> Peak1 0 0 0 0 0
#> Peak2 0 0 1 1 0
#> Peak3 1 0 1 1 0
#> Peak4 0 0 0 0 0
#> Peak5 0 0 0 0 0
# Information about peaks
head(rowData(sim))
#> DataFrame with 6 rows and 2 columns
#> Peak exp.peakmean
#> <character> <numeric>
#> Peak1 Peak1 6.75102e-06
#> Peak2 Peak2 2.68961e-04
#> Peak3 Peak3 1.02505e-04
#> Peak4 Peak4 3.80896e-06
#> Peak5 Peak5 4.83026e-05
#> Peak6 Peak6 4.76012e-05
# Information about cells
head(colData(sim))
#> DataFrame with 6 rows and 2 columns
#> Cell exp.libsize
#> <character> <numeric>
#> Cell1 Cell1 1467
#> Cell2 Cell2 575
#> Cell3 Cell3 1122
#> Cell4 Cell4 1105
#> Cell5 Cell5 925
#> Cell6 Cell6 887
# Peak by cell matrices
names(assays(sim))
#> [1] "counts"For more details about the SingleCellExperiment object refer to thevignette.
The simPICsimulate function provides additional simulation details:
* Cell information (colData)
* Cell - Unique cell identifier.
* exp.libsize - The expected library size for that cell. (not obtained from
the final simulated counts)
* Peak information (rowData)
* Peak - Unique peak identifier.
* exp.peakmean - The expected peak means for that peak. (not obtained from
the final simulated counts)
* Peak by cell information (assays)
* counts - The final simulated counts.
For more information on the simulation see ?simPICsimulate.
9 Comparing simulations and real data
simPIC provides a function simPICcompare that aims to make these comparisons
easier. This function takes a list of SingleCellExperiment objects, combines
the datasets and produces comparison plots. Let’s make two small simulations and
see how they compare.
sim1 <- simPICsimulate(nPeaks = 20000, nCells = 1000)
#> simPIC is:
#> updating parameters...
#> setting up SingleCellExperiment object...
#> Simulating library size...
#> Simulating peak mean...
#> using weibull distribution for simulating peak mean
#> Simulating true counts...
#> Done!!
sim2 <- simPICsimulate(nPeaks = 20000, nCells = 1000)
#> simPIC is:
#> updating parameters...
#> setting up SingleCellExperiment object...
#> Simulating library size...
#> Simulating peak mean...
#> using weibull distribution for simulating peak mean
#> Simulating true counts...
#> Done!!
comparison <- simPICcompare(list(real = sim1, simPIC = sim2))
names(comparison)
#> [1] "RowData" "ColData" "Plots"
names(comparison$Plots)
#> [1] "Means" "Variances" "MeanVar" "LibrarySizes" "ZerosPeak"
#> [6] "ZerosCell" "MeanZeros" "PeakMeanNzp"The returned list has three items. The first two are the combined datasets by
peak (RowData) and by cell (ColData) and the third contains some comparison
plots (produced using ggplot2), for example a plot of the distribution of
means:
comparison$Plots$Means
These are only a few of the plots you might want to consider but it should be easy to make more using the returned data.
10 Citing simPIC
If you use simPIC in your work please cite our paper:
citation("simPIC")
#> To cite package 'simPIC' in publications use:
#>
#> Chugh S, McCarthy D, Shim H (2024). _simPIC: simPIC: flexible
#> simulation of paired-insertion counts for single-cell ATAC-sequencing
#> data_. doi:10.18129/B9.bioc.simPIC
#> <https://doi.org/10.18129/B9.bioc.simPIC>, R package version 1.3.0,
#> <https://bioconductor.org/packages/simPIC>.
#>
#> A BibTeX entry for LaTeX users is
#>
#> @Manual{,
#> title = {simPIC: simPIC: flexible simulation of paired-insertion counts for single-cell
#> ATAC-sequencing data},
#> author = {Sagrika Chugh and Davis McCarthy and Heejung Shim},
#> year = {2024},
#> note = {R package version 1.3.0},
#> url = {https://bioconductor.org/packages/simPIC},
#> doi = {10.18129/B9.bioc.simPIC},
#> }Session information
sessionInfo()
#> 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] simPIC_1.3.0 SingleCellExperiment_1.29.0
#> [3] SummarizedExperiment_1.37.0 Biobase_2.67.0
#> [5] GenomicRanges_1.59.0 GenomeInfoDb_1.43.0
#> [7] IRanges_2.41.0 S4Vectors_0.45.0
#> [9] BiocGenerics_0.53.0 MatrixGenerics_1.19.0
#> [11] matrixStats_1.4.1 BiocStyle_2.35.0
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 xfun_0.48 bslib_0.8.0
#> [4] ggplot2_3.5.1 lattice_0.22-6 generics_0.1.3
#> [7] vctrs_0.6.5 tools_4.5.0 parallel_4.5.0
#> [10] fansi_1.0.6 tibble_3.2.1 highr_0.11
#> [13] pkgconfig_2.0.3 Matrix_1.7-1 checkmate_2.3.2
#> [16] actuar_3.3-4 lifecycle_1.0.4 GenomeInfoDbData_1.2.13
#> [19] farver_2.1.2 compiler_4.5.0 tinytex_0.53
#> [22] munsell_0.5.1 codetools_0.2-20 htmltools_0.5.8.1
#> [25] sass_0.4.9 yaml_2.3.10 pillar_1.9.0
#> [28] crayon_1.5.3 jquerylib_0.1.4 MASS_7.3-61
#> [ reached 'max' / getOption("max.print") -- omitted 40 entries ]