单细胞测序分析之scater包学习笔记
这次学习的是单细胞测序的第三个比较常用的包：Scater。本文主要参考教程：单细胞转录组学习笔记-15-利用scRNAseq包学习scater来走分析流程。
scater包的官方网站：here
这篇文章写的些许混乱，因为官网里的很多函数与教程里的并不一样，是因为有些函数是最新版scater新更新的，所以我在学习的时候有的部分按照官网的来，有的部分则按照教程里的来单细胞转录组学习笔记-15-利用scRNAseq包学习scater。但是实际上大同小异，主要的流程都是一样的。
这个包像之前介绍的两个包一样，是为了分析单细胞测序分析而开发的，也可以进行质控、可视化、下游分析前的预处理。scater包的流程是这样的：
简单的来说，scater包可以做下面几件事：
（1）QC
（2）从read数据里进行转录本定量（根据伪时间排序）
（3）数据格式标准化
（4）可视化
（5）Seamless integration into the Bioconductor universe
（6）标准化数据
安装scater包
>BiocManager::install("scater")
创建SingleCellExperiment对象
和之前介绍的Monocle2和Seurat包一样，scater包也需要一个对象。
SingleCellExperiment这个对象，画出来大概是这么个意思（图片来自：单细胞转录组学习笔记-15-利用scRNAseq包学习scater)：
在导入数据之前，请将表达数据存为矩阵。如果你的数据是10×Genomics的测序结果，可以用[DropletUtils]包的read10xCounts函数，它可以自动的生成SingleCellExperiment对象。这里我们用scater包自带的示例数据来走一遍流程。
首先，你要把数据导入进来：
第一个要导入的是基因的表达量矩阵，行是基因，列是细胞：
> rm(list = ls())
> options(stringsAsFactors = F)
> data("sc_example_counts") #导入表达矩阵
> dim(sc_example_counts) #看一下表达矩阵里多少行多少列
[1] 2000 40
> sc_example_counts[1:3,1:3] # 简单的看一下表达矩阵长什么样
Cell_001 Cell_002 Cell_003
Gene_0001 0 123 2
Gene_0002 575 65 3
Gene_0003 0 0 0
第二，导入细胞信息，注释信息可以是细胞名称，组织来源，收集样品的时间点，处理条件等等, 行是细胞，列是属性：
> data("sc_example_cell_info")
> sc_example_cell_info[1:3,1:3] # 样本注释信息
Cell Mutation_Status Cell_Cycle
Cell_001 Cell_001 positive S
Cell_002 Cell_002 positive G0
Cell_003 Cell_003 negative G1
有了上述的文件，就可以构建对象了：
> example_sce <- SingleCellExperiment(
assays = list(counts = sc_example_counts),
colData = sc_example_cell_info)
#瞅一眼对象里有啥
> example_sce
class: SingleCellExperiment
dim: 2000 40
metadata(0):
assays(1): counts
rownames(2000): Gene_0001 Gene_0002 ... Gene_1999 Gene_2000
rowData names(0):
colnames(40): Cell_001 Cell_002 ... Cell_039 Cell_040
colData names(4): Cell Mutation_Status Cell_Cycle Treatment
reducedDimNames(0):
spikeNames(0):
上面在构建对象的时候，我们将表达矩阵起了一个名字是“counts”，是为了方便有的时候你如果需要把矩阵里的counts提取出来，用counts函数比较方便：
> str(counts(example_sce))
int [1:2000, 1:40] 0 575 0 0 0 0 0 0 416 12 ...
- attr(*, "dimnames")=List of 2
..$ : chr [1:2000] "Gene_0001" "Gene_0002" "Gene_0003" "Gene_0004" ...
..$ : chr [1:40] "Cell_001" "Cell_002" "Cell_003" "Cell_004" ...
你可以通过下面的方法来新增一列meatdata：
> example_sce$whee <- sample(LETTERS, ncol(example_sce), replace=TRUE)
> colData(example_sce)
DataFrame with 40 rows and 5 columns
Cell Mutation_Status Cell_Cycle Treatment whee
<character> <character> <character> <character> <character>
Cell_001 Cell_001 positive S treat1 Z
Cell_002 Cell_002 positive G0 treat1 I
Cell_003 Cell_003 negative G1 treat1 W
Cell_004 Cell_004 negative S treat1 O
Cell_005 Cell_005 negative G1 treat2 I
... ... ... ... ... ...
Cell_036 Cell_036 negative G0 treat1 C
Cell_037 Cell_037 negative G0 treat1 E
Cell_038 Cell_038 negative G0 treat2 D
Cell_039 Cell_039 negative G1 treat1 F
Cell_040 Cell_040 negative G0 treat2 F
新增一行的方法：
#比如我想给每一行加上一个随机数，随机数生成的函数是runif()
> rowData(example_sce)$stuff <- runif(nrow(example_sce))
> rowData(example_sce)
DataFrame with 2000 rows and 1 column
stuff
<numeric>
Gene_0001 0.0491068968549371
Gene_0002 0.611523448256776
Gene_0003 0.282845608890057
Gene_0004 0.0250186347402632
Gene_0005 0.182038959115744
... ...
Gene_1996 0.633494263514876
Gene_1997 0.813403354492038
Gene_1998 0.987525492208079
Gene_1999 0.204807848902419
Gene_2000 0.369947757339105
还有一些其他的函数，以后如果有需要可以用的上，比如（摘自：单细胞转录组学习笔记-15-利用scRNAseq包学习scater)：
（1）如果你的数据里有spike-in control，可以用isSpike 对spike-in操作
（2）sizeFactors 是进行标准化时对细胞文库大小计算的结果
（3）reducedDim 对降维结果（reduced dimensionality results）操作
其他导入数据的方法
如果你的count矩阵是存在一个CSV文件里的，可以用read.table() 来读取。
如果是大的数据集，可以用"Matrix"包里的readSparseCounts()函数读取。
或者，如果你的数据来自10X Genomics测序，你可以用 read10xCounts函数读取，它会自动产生一个sparse矩阵。
另外，kallisto和Salmon的转录本丰度数据可以使用tximeta包来进行读取。
质量控制
scater可以从三个水平上进行质量控制（QC）：
（1）过滤细胞
（2）过滤基因
（3）实验变量的质控
细胞水平QC
细胞水平指标的计算是利用perCellQCMetrics()函数进行的，它包括：
（1）sum: 细胞里的总count数（文库大小）
（2）detected: 细胞里含有count的基因数
（3）subsets_X_percent: count的百分比（某个基因集）
NOTE这里需要安装最新版的scater(1.14.4)，如果是老版的scater，这一步是不能完成的
由于我的破电脑用了N种方法也没安装上最新版的scater，所以质控这些步骤我只把代码放了上来，自己并没有亲自运行一遍，如果有同学可以安装，最好要试一遍质控的过程，因为像下一步perCellQCMetrics函数是最新版里新增的功能:
library(scater)
per.cell <- perCellQCMetrics(example_sce,
subsets=list(Mito=grep("mt-", rownames(example_sce))))
summary(per.cell$sum)
summary(per.cell$detected)
summary(per.cell$subsets_Mito_percent)
colData(example_sce) <- cbind(colData(example_sce), per.cell)
教程里并没有用这个函数，而是用了另外一种：calculateQCMetrics
这个函数对细胞和基因都计算了很多种的指标，分别存在colData和rowData里，默认情况下对count值进行计算，你也可以通过参数exprs_values进行修改：
> example_sce <- calculateQCMetrics(example_sce)
> colnames(colData(example_sce)) #列是样品
[1] "Cell" "Mutation_Status" "Cell_Cycle"
[4] "Treatment" "whee" "is_cell_control"
[7] "total_features_by_counts" "log10_total_features_by_counts" "total_counts"
[10] "log10_total_counts" "pct_counts_in_top_50_features" "pct_counts_in_top_100_features"
[13] "pct_counts_in_top_200_features" "pct_counts_in_top_500_features"
> colnames(rowData(example_sce)) #行是基因
[1] "stuff" "is_feature_control" "mean_counts" #平均表达量 "log10_mean_counts"
[5] "n_cells_by_counts"#多少个细胞表达了该基因 "pct_dropout_by_counts"#基因在细胞中表达量为0的细胞数占比 "total_counts" "log10_total_counts"
如果你的测序里添加了ERCC，那么这里还可以设置对照：
> example_sce <- calculateQCMetrics(example_sce,
feature_controls = list(ERCC = 1:20, mito = 500:1000),
cell_controls = list(empty = 1:5, damaged = 31:40))
>
> all_col_qc <- colnames(colData(example_sce))
> all_col_qc <- all_col_qc[grep("ERCC", all_col_qc)]
> all_col_qc #这里就是ERCC相关的信息了
[1] "total_features_by_counts_ERCC" "log10_total_features_by_counts_ERCC"
[3] "total_counts_ERCC" "log10_total_counts_ERCC"
[5] "pct_counts_ERCC" "pct_counts_in_top_50_features_ERCC"
[7] "pct_counts_in_top_100_features_ERCC" "pct_counts_in_top_200_features_ERCC"
[9] "pct_counts_in_top_500_features_ERCC"
可视化
（一）细胞水平上的质控可以画一个总基因表达量 比上feature control的百分比，如果一个基因在总基因表达量上的比例多，在ERCC里少，就是正常的细胞，反之则不正常。
> plotColData(example_sce, x = "total_features_by_counts",
y = "pct_counts_feature_control", colour = "Mutation_Status") +
theme(legend.position = "top") +
stat_smooth(method = "lm", se = FALSE, size = 2, fullrange = TRUE)
（二）plotScater（单细胞转录组学习笔记-15-利用scRNAseq包学习scater）
先在表达量最高的基因中选一部分（默认是500个），然后从高到低累加，看看它们对每个细胞文库的贡献值如何。它将不同细胞的不同基因表达分布绘制出来，就像芯片数据或bulk转录组数据每个样本做一个箱线图，来看基因表达量分布。但是由于单细胞数据样本太多，没办法全画出箱线图。这个图在处理不同批次细胞时就可以很清楚地看到它们之间的差异。
为了看不同细胞的表达量差异，可以利用colData中的数据将细胞分类；默认使用counts值，但是只要在assays存在的表达矩阵，就可以在exprs_values参数中自定义:
> plotScater(example_sce, block1 = "Mutation_Status", block2 = "Treatment",
colour_by = "Cell_Cycle", nfeatures = 300, exprs_values = "counts")
image.png
（三）QC结果两两比较（单细胞转录组学习笔记-15-利用scRNAseq包学习scater）
# 比较feature的属性
> plotRowData(example_sce, x = "n_cells_by_counts", y = "mean_counts")
# 比较样本的属性
> p1 <- plotColData(example_sce, x = "total_counts",
y = "total_features_by_counts")
> p2 <- plotColData(example_sce, x = "pct_counts_feature_control",
y = "total_features_by_counts")
> p3 <- plotColData(example_sce, x = "pct_counts_feature_control",
y = "pct_counts_in_top_50_features")
# 最后组合
> multiplot(p1, p2, p3, cols = 3)
过滤低质量细胞
举个栗子，我们只保留有10万个counts和500个基因表达的细胞：
> keep.total <- example_sce$total_counts > 1e5
> keep.n <- example_sce$total_features_by_counts > 500
> filtered <- example_sce[,keep.total & keep.n]
> dim(filtered)
[1] 1973 37
isOutlier 函数提供了一个更为合适的方法来选择阈值。它定义距离中位数一定数量的中位数绝对偏差(MADs)的阈值。超过这个值的细胞认为是离群的，从而被过滤掉：
> keep.total <- isOutlier(example_sce$total_counts, nmads=3,
type="lower", log=TRUE)
> filtered <- example_sce[,keep.total]
使用quickPerCellQC()函数还可以更方便地检测几个常见指标的异常值。该方法使用总count数、检测到的基因数和某个基因集的count百分比(如线粒体基因、spik -in转录本)来确定要舍弃哪些细胞：
qc.stats <- quickPerCellQC(per.cell, percent_subsets="subsets_Mito_percent")
colSums(as.matrix(qc.stats))
filtered <- example_sce[,!qc.stats$discard]
isOutlier的方法需要根据你的实验需求进行调整，如测序深度、spike-in的加入、细胞类型。
你也可以对基因进行过滤：
> keep_feature <- nexprs(example_sce, byrow=TRUE) >= 4
> example_sce <- example_sce[keep_feature,]
> dim(example_sce)
[1] 1753 40
基因水平的质控
指标的定义
使用的perFeatureQCMetrics()函数包括下面几个指标：
（1） mean: 所有的细胞之间基因的平均count
（2）detected: 对于每一个基因都是非零的细胞百分比
（3）subsets_Y_ratio: 某一个子集的细胞和总细胞的平均count的比值
# 这里举例：前10个孔是空孔
per.feat <- perFeatureQCMetrics(example_sce, subsets=list(Empty=1:10))
summary(per.feat$mean)
summary(per.feat$detected)
summary(per.feat$subsets_Empty_ratio)
calculateAverage()函数提供了更精确的平均值计算，它根据取平均值之前的相对文库大小来调整计数：
#这里我用的是过滤前的对象来进行演示的
> ave <- calculateAverage(example_sce)
> summary(ave)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.009 5.529 67.730 197.135 190.286 12428.798
我们还可以直接计算出表达基因的细胞数量：
> summary(nexprs(example_sce, byrow=TRUE))
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 6.00 12.00 14.66 21.00 40.00
可视化：
在默认值下，只显示前50个基因。行是基因，橙色的线(bar)代表该基因在每一个细胞里的表达量。圆圈代表这个基因表达量的中位数。
> plotHighestExprs(example_sce, exprs_values = "counts")
下面这个图是对于每一个基因都是非零的细胞百分比：
> plotExprsFreqVsMean(example_sce)
根据行来取子集
从SingleCellExperiment对象里移除一些基因，比如我们把在所有细胞中都不表达的基因删除：
> keep_feature <- rowSums(counts(example_sce) > 0) > 0
> example_sce <- example_sce[keep_feature,]
> dim(example_sce)
[1] 1973 40
也可以用现有的基因注释进行过滤，比如线粒体基因，通常我们对这些基因不感兴趣，对我们了解细胞异质性也没有多大帮助。
变量水平的质控
用getVarianceExplained()函数(经过标准化之后，请参阅下面的内容)。
example_sce <- logNormCounts(example_sce) # see below.
vars <- getVarianceExplained(example_sce,
variables=c("tissue", "total mRNA mol", "sex", "age"))
head(vars)
plotExplanatoryVariables(vars)
计算表达值
标准化文库大小差异
第一种方法：
最常用的函数是logNormCounts()，结果储存在"logcounts"里。它是用每一个count除以自身的size factor，加一个 pseudo-count 和log-transforming。
example_sce <- logNormCounts(example_sce)
assayNames(example_sce)
summary(librarySizeFactors(example_sce))
第二种方法：
用calculateCPM()函数，结果储存在"cpm"中。相关函数还包括：calculateTPM()和calculateFPKM()。
cpm(example_sce) <- calculateCPM(example_sce)
assay(example_sce, "normed") <- normalizeCounts(example_sce,
size_factors=runif(ncol(example_sce)), pseudo_count=1.5)
上面两种方法是官网的代码，我用的是下面这个，是教程里用的代码：
> cpm(example_sce) <- calculateCPM(example_sce)
> example_sce <- normalize(example_sce)
> assayNames(example_sce)
[1] "counts" "cpm" "logcounts"
> plotExplanatoryVariables(example_sce)
也可以选择一部分进行画图：
> plotExplanatoryVariables(example_sce,
variables = c("total_features_by_counts", "total_counts",
"Mutation_Status", "Treatment", "Cell_Cycle"))
结果看到：total_counts 和 total_features_by_counts 对基因表达变化差异贡献很高（接近10%），真实数据中这两个占比应该在1-5%之间
可视化表达值
plotExpression()函数可以画某一个基因子集的表达值，默认使用"logcounts"里的数据，也可以通过添加exprs_values参数来修改。
> plotExpression(example_sce, rownames(example_sce)[1:6])
这是前6个基因的logcount值的图
你也可以按照分组的情况画基因表达情况，比如在这个包自带的数据info里，有Mutation_Status，Cell_cycle的分组情况，你也可以根据这些来画图：
> plotExpression(example_sce, rownames(example_sce)[1:6],
x = "Mutation_Status", exprs_values = "logcounts")
> plotExpression(example_sce, rownames(example_sce)[1:6],
x = "Cell_Cycle", exprs_values = "logcounts")
你还可以同时将两个变量一起画出来：
> plotExpression(example_sce, rownames(example_sce)[1:6],
colour_by = "Cell_Cycle", shape_by = "Mutation_Status",
size_by = "Gene_0002")
降维
（一）PCA
最常用的是PCA降维的方法，PCA结果将存放在SingleCellExperiment对象的reducedDims里：
> example_sce <- runPCA(example_sce)
> str(reducedDim(example_sce, "PCA"))
num [1:40, 1:2] -11.35 -1.41 2.01 -9.8 6.35 ...
- attr(*, "dimnames")=List of 2
..$ : chr [1:40] "Cell_001" "Cell_002" "Cell_003" "Cell_004" ...
..$ : chr [1:2] "PC1" "PC2"
- attr(*, "percentVar")= num [1:2] 0.0679 0.0422
默认情况下，runPCA()函数用第一个PCs（聚类）里变化量最大的500个基因来进行计算。你也可以用subset_row参数来调整你想要进行计算的基因数，或者基因子集。用name参数可以修改在reducedDims slot里的名称：
> example_sce <- runPCA(example_sce, name="PCA2",
subset_row=rownames(example_sce)[1:1000],
ncomponents=25)
（二）其他降维方法(tSNE, UMAP)
> set.seed(1000)
> example_sce <- runTSNE(example_sce, perplexity=10)
> head(reducedDim(example_sce, "TSNE"))
[,1] [,2]
Cell_001 3.7408876 -0.8718345
Cell_002 0.9961570 -1.0279169
Cell_003 0.3328361 1.1045595
Cell_004 4.9745275 0.2320463
Cell_005 -4.1067009 0.1909082
Cell_006 0.5936069 -2.1781319
这里官方网站推荐多用几个不同的随机种子和perplexity值（表示近邻细胞的数量）来跑tSNE，以确保你的结果是比较稳定的。
一个更常见的模式是用你之前跑的PCA结果，作为tSNE的输入，进行分析。这有助于提高运行速度，并且降低噪音，只集中在几个主要的变量上。下面这个代码是用的PCA前10个dimension的结果来进行tSNE计算：
> set.seed(1000)
> example_sce <- runTSNE(example_sce, perplexity=50,
dimred="PCA", n_dimred=10)
> head(reducedDim(example_sce, "TSNE"))
你也可以用UMAP来计算：
> example_sce <- runUMAP(example_sce)
> head(reducedDim(example_sce, "UMAP"))
[,1] [,2]
Cell_001 2.1841452 -0.4125833
Cell_002 0.2940562 1.0461869
Cell_003 -0.5211948 1.2969395
Cell_004 2.1369329 0.1383032
Cell_005 -2.0110685 -0.4748272
Cell_006 1.0558627 0.7567186
你还可以用diffusion maps来降维：
> example_sce <- runDiffusionMap(example_sce)
Registered S3 method overwritten by 'xts':
method from
as.zoo.xts zoo
> plotDiffusionMap(example_sce, colour_by = "Gene_0001", size_by = "Gene_1000")
降维的可视化
#根据Mutation_Status来画
> plotReducedDim(example_sce, use_dimred = "PCA", colour_by = "Mutation_Status")
还可以利用表达量添加颜色、大小分组（参考单细胞转录组学习笔记-15-利用scRNAseq包学习scater）：
> # 看特定基因在PCA过程中起到的作用
> plotReducedDim(example_sce, use_dimred = "PCA",
colour_by = "Gene_1000", size_by = "Gene_0500")
如果你想同时画2个以上成分的可视化图的时候，可以用下面的代码：
> example_sce <- runPCA(example_sce, ncomponents=20)
> plotPCA(example_sce, ncomponents = 4, colour_by = "Mutation_Status")
图里的对角线上的方块图，显示的是每一个component的密度
