---
title: "simPIC: simulating single-cell ATAC-seq data"
author: "Sagrika Chugh"
package: simPIC
date: "`r Sys.Date()`"
output:
BiocStyle::html_document:
toc: true
toc_float: true
vignette: >
\usepackage[utf8]{inputenc}
%\VignetteIndexEntry{simPIC: simulating single-cell ATAC-seq data}
%\VignetteEngine{knitr::rmarkdown}
editor_options:
markdown:
wrap: 80
---
```{r knitr-options, echo = FALSE, message = FALSE, warning = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
options(max.print = 30)
```
# 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.
# Installation
simPIC can be installed from Bioconductor:
```{r install-bioc, eval = FALSE, include=TRUE}
BiocManager::install("simPIC")
```
# 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.
```{r quickstart}
# 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)
# Simulate data using estimated parameters
sim <- simPICsimulate(est)
```
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.
# 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 `r BiocStyle::Biocpkg("PICsnATAC")`
[vignette][PIC-vignette]. Briefly, you need three input files for PIC -
1. cell barcodes with metadata (`singlecell.csv`)
2. list of peak regions (`peaks.bed`)
3. fragment files (`fragment.tsv.gz`)
```{r pic, eval=FALSE, include=TRUE}
pic_mat <- PIC_counting(
cells = cells,
fragment_tsv_gz_file_location = fragment_tsv_gz_file_location,
peak_sets = peak_sets
)
```
# 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](https://en.wikipedia.org/wiki/Weibull_distribution) with default
settings. Users have the flexibility to choose from `gamma`, `lognormal-gamma` 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](https://en.wikipedia.org/wiki/Poisson_distribution).
# 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.
```{r simPICparams}
sim.params <- newsimPICcount()
```
we can see the default values for the `simPICcount` object parameters. These
values are based on provided test data.
```{r params}
sim.params
```
## Getting and setting parameters
To get a particular parameter, for e.g., number of peaks, we can use `simPICget`
function:
```{r getParam}
simPICget(sim.params, "nPeaks")
```
Alternatively, to give a parameter a new value we can use the
`setsimPICparameters` function:
```{r setParam}
sim.params <- setsimPICparameters(sim.params, nPeaks = 2000)
simPICget(sim.params, "nPeaks")
```
To get or set multiple parameters use `simPICget` or `setsimPICparameters`
functions:
```{r getParams-setParams}
# 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
```
# 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.
```{r simPICestimate}
# Get the counts from test data
count <- readRDS(system.file("extdata", "test.rds", package = "simPIC"))
# Check that counts is a dgCMatrix
class(count)
typeof(count)
# Check the dimensions, each row is a peak, each column is a cell
dim(count)
# Show the first few entries
count[1:5, 1:5]
new <- newsimPICcount()
new <- simPICestimate(count)
## 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:
1. Mean parameters are estimated by fitting a weibull distribution (default) to
the peak means.
2. Library size parameters are estimated by fitting a log-normal distribution
to the library sizes.
3. Sparsity parameter is estimated by fitting a Bernoulli distribution.
For more details of the estimation procedures see `?simPICestimate`.
# 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:
```{r simPICsimulate}
sim <- simPICsimulate(new, nCells = 1000)
sim
## 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 `r nrow(sim)` features (peaks) and
`r ncol(sim)` 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.
```{r SCE}
# Access the counts
counts(sim)[1:5, 1:5]
# Information about peaks
head(rowData(sim))
# Information about cells
head(colData(sim))
# Peak by cell matrices
names(assays(sim))
```
For more details about the `r BiocStyle::Biocpkg("SingleCellExperiment")`
object refer to the[vignette][SCE-vignette].
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`.
# 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.
```{r comparison}
sim1 <- simPICsimulate(nPeaks = 20000, nCells = 1000)
sim2 <- simPICsimulate(nPeaks = 20000, nCells = 1000)
comparison <- simPICcompare(list(real = sim1, simPIC = sim2))
names(comparison)
names(comparison$Plots)
```
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:
```{r comparison-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.
# Citing simPIC
If you use simPIC in your work please cite our paper:
```{r citation}
citation("simPIC")
```
# Session information {.unnumbered}
```{r sessionInfo}
sessionInfo()
```
[SCE-vignette]:https://bioconductor.org/packages/devel/bioc/vignettes/SingleCellExperiment/inst/doc/intro.html
[PIC-vignette]:https://github.com/Zhen-Miao/PICsnATAC/blob/main/vignettes/Run_PIC_counting_on_pbmc_3k_data.ipynb).