Splatter是一个模拟单细胞RNA测序计数数据的软件包。它提供了一个简单的界面，用于创建可复制且文档充分的复杂模拟。可以从真实数据估计参数，并提供用于比较真实数据集和模拟数据集的函数。

# if (!require("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# BiocManager::install("splatter")library(splatter)
library(scater)### 1. 快速创建模拟数据
# Usage
# mockSCE(ncells = 200, ngenes = 2000, nspikes = 100)
set.seed(1)
sce <- mockSCE() # class: SingleCellExperiment
counts(sce)### 2. 参数估计
# 从数据（模拟数据或真实的单细胞测序数据）中估计参数
params <- splatEstimate(sce)  # Params object### 3. 自定义参数，创建模拟数据
# 创建参数对象
params <- newSplatParams() #默认10000 Genes，100 Cells
# 查看参数对象中的参数设定
getParam(params, "nGenes")
getParam(params, "nCells")
getParams(params, c("nGenes", "mean.rate", "mean.shape"))
# 设定参数对象中的参数
params <- setParam(params, "nGenes", 5000)
params <- setParam(params, "nCells", 50)
params <- setParams(params, mean.shape = 0.5, de.prob = 0.2)
params <- setParams(params, update = list(nGenes = 8000, mean.rate = 0.5))# 根据参数，创建模拟数据
sim <- splatSimulate(params)### 4. 提取SingleCellExperiment类信息
class(sim)
# Information about genes
head(rowData(sim))
# Information about cells
head(colData(sim))
# Gene by cell matrices
names(assays(sim))
# Example of cell means matrix
assays(sim)$CellMeans[1:5, 1:5]
# 表达矩阵
counts <- counts(sim)
counts[1:3,1:5]class(counts) # "matrix" "array"
typeof(counts)  # [1] "integer"
dim(counts) #[1] 8000  100### 5. 数据标准化，对数化和降维
## {scuttle}包的函数:logNormCounts
## {scater}包的函数:runPCA,plotPCA
#  SingleCellExperiment 对象的counts值计算对数转换的归一化表达式值
sim <- logNormCounts(sim)
counts(sim)[1:3,1:5]
logcounts(sim)[1:3,1:5]
# PCA降维
sim <- runPCA(sim)
plotPCA(sim)### 6. 创建复杂的模拟数据
## groups
sim.groups <- splatSimulate(group.prob = c(0.3, 0.7), method = "groups",verbose = FALSE)
sim.groups <- logNormCounts(sim.groups)
sim.groups <- runPCA(sim.groups)
plotPCA(sim.groups, colour_by = "Group")dim(counts(sim.groups))
rowData(sim.groups)$DEFacGroup1
rowData(sim.groups)$DEFacGroup2
metadata(sim.groups)## paths
sim.paths <- splatSimulate(de.prob = 0.2, nGenes = 1000, method = "paths",verbose = FALSE)
sim.paths <- logNormCounts(sim.paths)
sim.paths <- runPCA(sim.paths)
plotPCA(sim.paths, colour_by = "Step")
colData(sim.paths)$Step## batches
sim.batches <- splatSimulate(batchCells = c(50, 50), verbose = FALSE)
sim.batches <- logNormCounts(sim.batches)
sim.batches <- runPCA(sim.batches)
plotPCA(sim.batches, colour_by = "Batch")
rowData(sim.batches)$BatchFacBatch1
rowData(sim.batches)$BatchFacBatch2### 7.比较SingleCellExperiment对象
## 相互比较
sim1 <- splatSimulate(nGenes = 1000, batchCells = 20, verbose = FALSE)
sim2 <- simpleSimulate(nGenes = 1000, nCells = 20, verbose = FALSE)
comparison <- compareSCEs(list(Splat = sim1, Simple = sim2))names(comparison)
comparison$RowData
comparison$ColDatanames(comparison$Plots)
comparison$Plots$Means
comparison$Plots$LibrarySizes
comparison$Plots$ZerosCell## 与参考进行比较
sim1 <- splatSimulate(nGenes = 1000, batchCells = 100, verbose = FALSE)
sim2 <- splatSimulate(nGenes = 1000, batchCells = c(40, 60), verbose = FALSE)
sim3 <- simpleSimulate(nGenes = 1000, nCells = 100, verbose = FALSE)
difference <- diffSCEs(list(Splat1 = sim1, Splat2 = sim2, Simple = sim3), ref = "Simple")
difference$Plots$Means
difference$QQPlots$Means### 8.批次效应参数设置
library("splatter")
library("scater")
library("ggplot2")# Simulation with small batch effects
sim1 <- splatSimulate(params, batchCells = c(100, 100),batch.facLoc = 0.001, batch.facScale = 0.001,verbose = FALSE)
sim1 <- logNormCounts(sim1)
sim1 <- runPCA(sim1)
plotPCA(sim1, colour_by = "Batch") + ggtitle("Small batch effects")# Simulation with big batch effects
sim2 <- splatSimulate(params, batchCells = c(100, 100),batch.facLoc = 0.5, batch.facScale = 0.5,verbose = FALSE)
sim2 <- logNormCounts(sim2)
sim2 <- runPCA(sim2)
plotPCA(sim2, colour_by = "Batch") + ggtitle("Big batch effects")# 是否消除批次效应设置
sim1 <- splatSimulate(params, batchCells = c(100, 100), batch.rmEffect = FALSE,
verbose = FALSE)
sim1 <- logNormCounts(sim1)
sim1 <- runPCA(sim1)
plotPCA(sim1, colour_by = "Batch") + ggtitle("With batch effects")sim2 <- splatSimulate(params, batchCells = c(100, 100), batch.rmEffect = TRUE,verbose = FALSE)
sim2 <- logNormCounts(sim2)
sim2 <- runPCA(sim2)
plotPCA(sim2, colour_by = "Batch") + ggtitle("Batch effects removed")### 9. 离群值参数设置
## ----outlier-prob-------------------------------------------------------------
# Few outliers
sim1 <- splatSimulate(out.prob = 0.001, verbose = FALSE)
ggplot(as.data.frame(rowData(sim1)),aes(x = log10(GeneMean), fill = OutlierFactor != 1)) +geom_histogram(bins = 100) +ggtitle("Few outliers")# Lots of outliers
sim2 <- splatSimulate(out.prob = 0.2, verbose = FALSE)
ggplot(as.data.frame(rowData(sim2)),aes(x = log10(GeneMean), fill = OutlierFactor != 1)) +geom_histogram(bins = 100) +ggtitle("Lots of outliers")### 10.设置组参数
# One small group, one big group
params.groups <- newSplatParams(batchCells = 500, nGenes = 1000)
sim1 <- splatSimulateGroups(params.groups, group.prob = c(0.9, 0.1),verbose = FALSE)
sim1 <- logNormCounts(sim1)
sim1 <- runPCA(sim1)
plotPCA(sim1, colour_by = "Group") + ggtitle("One small group, one big group")# Five groups
sim2 <- splatSimulateGroups(params.groups,group.prob = c(0.2, 0.2, 0.2, 0.2, 0.2),verbose = FALSE)
sim2 <- logNormCounts(sim2)
sim2 <- runPCA(sim2)
plotPCA(sim2, colour_by = "Group") + ggtitle("Five groups")### 11. 设置差异表达基因参数
## de.prob 参数
# Few DE genes
params.groups <- newSplatParams(batchCells = 500, nGenes = 1000)
sim1 <- splatSimulateGroups(params.groups, group.prob = c(0.5, 0.5),de.prob = 0.01, verbose = FALSE)
sim1 <- logNormCounts(sim1)
sim1 <- runPCA(sim1)
plotPCA(sim1, colour_by = "Group") + ggtitle("Few DE genes")# Lots of DE genes
sim2 <- splatSimulateGroups(params.groups, group.prob = c(0.5, 0.5),de.prob = 0.3, verbose = FALSE)
sim2 <- logNormCounts(sim2)
sim2 <- runPCA(sim2)
plotPCA(sim2, colour_by = "Group") + ggtitle("Lots of DE genes")## de.facLoc参数
# Small DE factors
sim1 <- splatSimulateGroups(params.groups, group.prob = c(0.5, 0.5),de.facLoc = 0.01, verbose = FALSE)
sim1 <- logNormCounts(sim1)
sim1 <- runPCA(sim1)
plotPCA(sim1, colour_by = "Group") + ggtitle("Small DE factors")# Big DE factors
sim2 <- splatSimulateGroups(params.groups, group.prob = c(0.5, 0.5),de.facLoc = 0.3, verbose = FALSE)
sim2 <- logNormCounts(sim2)
sim2 <- runPCA(sim2)
plotPCA(sim2, colour_by = "Group") + ggtitle("Big DE factors")## complex-de---------------------------------------------------------------
# Different DE probs
sim1 <- splatSimulateGroups(params.groups,group.prob = c(0.2, 0.2, 0.2, 0.2, 0.2),de.prob = c(0.01, 0.01, 0.1, 0.1, 0.3),verbose = FALSE)
sim1 <- logNormCounts(sim1)
sim1 <- runPCA(sim1)
plotPCA(sim1, colour_by = "Group") +labs(title = "Different DE probabilities",caption = paste("Groups 1 and 2 have very few DE genes,","Groups 3 and 4 have the default number,","Group 5 has many DE genes"))# Different DE factors
sim2 <- splatSimulateGroups(params.groups,group.prob = c(0.2, 0.2, 0.2, 0.2, 0.2),de.facLoc = c(0.01, 0.01, 0.1, 0.1, 0.2),de.facScale = c(0.2, 0.5, 0.2, 0.5, 0.4),verbose = FALSE)
sim2 <- logNormCounts(sim2)
sim2 <- runPCA(sim2)
plotPCA(sim2, colour_by = "Group") +labs(title = "Different DE factors",caption = paste("Group 1 has factors with small location (value),","and scale (variability),","Group 2 has small location and greater scale.\n","Groups 3 and 4 have greater location with small,","and large scales","Group 5 has bigger factors with moderate","variability"))# Combination of everything
sim3 <- splatSimulateGroups(params.groups,group.prob = c(0.05, 0.2, 0.2, 0.2, 0.35),de.prob = c(0.3, 0.1, 0.2, 0.01, 0.1),de.downProb = c(0.1, 0.4, 0.9, 0.6, 0.5),de.facLoc = c(0.6, 0.1, 0.1, 0.01, 0.2),de.facScale = c(0.1, 0.4, 0.2, 0.5, 0.4),verbose = FALSE)
sim3 <- logNormCounts(sim3)
sim3 <- runPCA(sim3)
plotPCA(sim3, colour_by = "Group") +labs(title = "Different DE factors",caption = paste("Group 1 is small with many very up-regulated DE genes,","Group 2 has the default DE parameters,\n","Group 3 has many down-regulated DE genes,","Group 4 has very few DE genes,","Group 5 is large with moderate DE factors"))### 12. path参数设置
# Linear paths
params.groups <- newSplatParams(batchCells = 500, nGenes = 1000)
sim1 <- splatSimulatePaths(params.groups,group.prob = c(0.25, 0.25, 0.25, 0.25),de.prob = 0.5, de.facLoc = 0.2,path.from = c(0, 1, 2, 3),verbose = FALSE)
sim1 <- logNormCounts(sim1)
sim1 <- runPCA(sim1)
plotPCA(sim1, colour_by = "Group") + ggtitle("Linear paths")# Branching path
sim2 <- splatSimulatePaths(params.groups,group.prob = c(0.25, 0.25, 0.25, 0.25),de.prob = 0.5, de.facLoc = 0.2,path.from = c(0, 1, 1, 3),verbose = FALSE)
sim2 <- logNormCounts(sim2)
sim2 <- runPCA(sim2)
plotPCA(sim2, colour_by = "Group") + ggtitle("Branching path")## ----paths-steps------------------------
# Few steps
sim1 <- splatSimulatePaths(params.groups, path.nSteps = 3,de.prob = 0.5, de.facLoc = 0.2, verbose = FALSE)
sim1 <- logNormCounts(sim1)
sim1 <- runPCA(sim1)
plotPCA(sim1, colour_by = "Step") + ggtitle("Few steps")# Lots of steps
sim2 <- splatSimulatePaths(params.groups, path.nSteps = 1000,de.prob = 0.5, de.facLoc = 0.2, verbose = FALSE)
sim2 <- logNormCounts(sim2)
sim2 <- runPCA(sim2)
plotPCA(sim2, colour_by = "Step") + ggtitle("Lots of steps")## ----paths-skew-------------------
# Skew towards the end
sim1 <- splatSimulatePaths(params.groups, path.skew = 0.1,de.prob = 0.5, de.facLoc = 0.2, verbose = FALSE)
sim1 <- logNormCounts(sim1)
sim1 <- runPCA(sim1)
plotPCA(sim1, colour_by = "Step") + ggtitle("Skewed towards the end")

参考

https://www.bioconductor.org/packages/release/bioc/html/splatter.html

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5596896/