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$ColData
names(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