处理单细胞? Bioconductor就够用了

最近，bioconductor团队出版了一本电子书，其整合了其网站上关于单细胞的R包并制定了一套常规的分析流程包括分析，可视化，导入导出。不仅如此，前三章还分别教你如何下载使用R，使用bioconductor网站以及如何设计单细胞实验，对初学者很友好了，哪怕你对R语言一窍不懂，也能跟着走完流程。

《Orchestrating Single-Cell Analysis with Bioconductor》

1.简介

1.1 能学到什么？

能看这本书的都是对单细胞测序有所需求或这有这个意愿去学习相关知识的。这本书主要是整合目前常见的单细胞分析流程并尽可能详细的解释这些流程的每一个步骤，包括原理，所使用的工具以及给出栗子。所以可以根据自己的实际需求取选取合适自己的workflow与其中的步骤。

1.2 不能学到什么？

当然是与单细胞测序分析无关的其他分析，毕竟一口不能吃成大胖子。

2.R与Bioconductor(for beginner)

因为这本书关于单细胞分析的软件都是用R写的。如果对R完全不了解，你也不知道做出来的到底什么鬼，所以还是简单的介绍一下必备的背景吧。（最有效率的办法还是跟着生信技能树的公众号去学习！）

• R 是一门高水平的编程语言，主要用于数据分析
• 网上有很多关于R如何使用的教程例如 Learn R series，R for Data Science以及tidyverse等，但都是全英文的。看不懂的话，就跟着Jimmy老师团队学习吧，保证入股不亏。
• 本地安装 R 就跟常规软件安装是一样的,但是最好还是安装 RStudio ，因为界面更加人性化。想想你的大脑要处理复杂的编程语言，为什么不让它轻松一点呢。
• R语言现在有很多大神编写好的R包，对应各种各样的功能。
• 目前主流的R包三大下载软件：

1. CRAN
2. github(。。。大陆的孩子一言难尽)
3. 以及bioconductor

P.S 最近bioconductor容易出现联网失败的问题，具体的解决方法见Jimmy老师的推文【紧急通知】下载R包却联网失败？初学者的痛

3.R进阶

• 这里推荐几本电子书：What They Forgot to Teach You About R ，R Inferno
• 以及Jimmy老师推荐的实体书《R语言实战》
• 记住一本书看5遍or看5本书

4.SingleCellExperiment

这里有个坑注意一下，这本书的R包要求bioconductor version至少要3.10以上才能用，所以可以一键升级所有R包

#####
rm(list = ls())
options(stringsAsFactors = F)
options(download.file.method = 'libcurl')
options(url.method='libcurl')
chooseBioCmirror()###记得选bioc镜像，清华or中科大
BiocManager::valid() ###测试R包是否过期
BiocManager::install(version = "3.10")###升级bioconductor版本

但是时间比较久，因为要升级全部的R包，建议中午睡个午觉，下午起来就搞定了。

SingleCellExperiment

SingleCellExperiment是一个能够完成常用的70种单细胞测序R包的数据交换的R包，因为它能存储一个sc-seq数据的所有数据在一个class里，可以放单个细胞数据，也可以放不同细胞间基因表达量数据，还可以放不同批次的细胞等等，一个数据就是一个assay。class中分门别类的放好，而且方便调用，一个函数就能去你想要得数据。举个生活中的例子，class就相当于支付宝中的余额宝，数据就相当于你的钱，随时可以从不同银行（不同R包）转入转出，调用的函数就等于你的密码。

本章用的6个功能

1586590300873

一个实战的例子

这里选用书中最简单的列子来说明，注意的是SingleCellExperiment输入必须是矩阵，就像你存在支付宝的钱必须是人民币一样

4.1.存入数据

rm(list = ls())
options(stringsAsFactors = F)
options(download.file.method = 'libcurl')
options(url.method='libcurl')
library(scater)
library(SingleCellExperiment)
####
###建立一个矩阵，行名是基因，列名是细胞，注意的是输入必须是矩阵格式
counts_matrix <- data.frame(cell_1 = rpois(10, 10), 
                            cell_2 = rpois(10, 10), 
                            cell_3 = rpois(10, 30))
rownames(counts_matrix) <- paste0("gene_", 1:10)
counts_matrix <- as.matrix(counts_matrix) # must be a matrix object!
head(counts_matrix)
#        cell_1 cell_2 cell_3
# gene_1     10      9     29
# gene_2      8     13     36
# gene_3      4      9     36
# gene_4      6     11     42
# gene_5     10      8     26
# gene_6      4      9     35
###之后把矩阵放入class中，记得给list命名，这样方便提取
sce <- SingleCellExperiment(assays = list(counts = counts_matrix))
sce
# class: SingleCellExperiment  是什么class
# dim: 10 3 维度
# metadata(0): 
# assays(1): counts 
# rownames(10): gene_1 gene_2 ... gene_9 gene_10 行名
# rowData names(0): 基因的元数据的行名
# colnames(3): cell_1 cell_2 cell_3 列名
# colData names(0): 细胞的元数据的列名
# reducedDimNames(0): 
# spikeNames(0):
# altExpNames(0): 

上面就是存入，对应的就是取出

4.2.取出数据

有2种方法，一种是完整函数，一种是更便捷的办法（类似于密码与指纹支付）

assay(sce, "counts")##取出assay中名为counts的矩阵
counts(sce)##取出assay中名为counts的矩阵

第二种方法名字一定要一致，例如”counts“只能取”counts“，就像别人的指纹不可能转走你的钱，对吧。

4.3.编辑数据

SingleCellExperiment支持直接编辑数据，就像你不需要把余额宝的钱转支付宝，就可以在淘宝购物一样。但有一点不一样的是，不会覆盖你的原始数据，会作为一个新的assay存放在class中

sce <- scater::logNormCounts(sce)##对sce的数据进行log化和正态化
# class: SingleCellExperiment 
# dim: 10 3 
# metadata(0): 
# assays(2): counts logcounts 可以看到这里多了个名字为logcounts的assay
# rownames(10): gene_1 gene_2 ... gene_9 gene_10
# rowData names(0):
# colnames(3): cell_1 cell_2 cell_3
# colData names(0):
# reducedDimNames(0):
# spikeNames(0):
# altExpNames(0):
###提取log数据
logcounts(sce)
#           cell_1   cell_2   cell_3
# gene_1  4.473275 4.146048 3.961101
# gene_2  4.167494 4.651269 4.254920
# gene_3  3.245625 4.146048 4.254920
# gene_4  3.778973 4.420662 4.466477
# gene_5  4.473275 3.986273 3.814209
# gene_6  3.245625 4.146048 4.216435
# gene_7  4.328471 3.986273 4.006919
# gene_8  3.778973 4.146048 3.591729
# gene_9  3.986273 3.806574 4.254920
# gene_10 4.836720 3.361811 4.136242
###可以看到数据已经log+正态化

4.4.元数据

元数据，作为一种注释的数据，注释试验批次，实验条件等等，类似于你的工资条，告诉你你的工资构成。注意元数据是以数据框的形式放入class里

添加细胞的元数据

cell_metadata <- data.frame(batch = c(1, 1, 2))
rownames(cell_metadata) <- paste0("cell_", 1:3)
cell_metadata
#        batch
# cell_1     1
# cell_2     1
# cell_3     2

sce <- SingleCellExperiment(assays = list(counts = counts_matrix),
    colData = cell_metadata) ##为counts数据附上注释信息
sce
# class: SingleCellExperiment 
# dim: 10 3 
# metadata(0):
# assays(1): counts
# rownames(10): gene_1 gene_2 ... gene_9 gene_10
# rowData names(0):
# colnames(3): cell_1 cell_2 cell_3
# colData names(1): batch  ###细胞的元数据的列名为batch
# reducedDimNames(0):
# spikeNames(0):
# altExpNames(0):
##元数据也可以用2种方法取出，一种取出来的是数据框，另一种则是向量
colData(sce)
# DataFrame with 3 rows and 1 column
#            batch
#        <numeric>
# cell_1         1
# cell_2         1
# cell_3         2
sce$batch
# [1] 1 1 2

当然也可以用函数添加更多的元数据例如质控数据，总数等等。

sce <- scater::addPerCellQC(sce)
colData(sce)[, 1:5]
## DataFrame with 3 rows and 5 columns
##            batch       sum  detected percent_top_50 percent_top_100
##        <numeric> <integer> <integer>      <numeric>       <numeric>
## cell_1         1        86        10            100             100
## cell_2         1       104        10            100             100
## cell_3         2       280        10            100             100

添加基因的元数据

跟上面大同小异，调用的是不同的函数而已

sce <- scater::addPerFeatureQC(sce)##添加基因质控元数据
rowData(sce)
# DataFrame with 3 rows and 9 columns
#            batch       sum  detected percent_top_50 percent_top_100 percent_top_200
#        <numeric> <integer> <integer>      <numeric>       <numeric>       <numeric>
# cell_1         1        77        10            100             100             100
# cell_2         1        88        10            100             100             100
# cell_3         2       325        10            100             100             100

4.5.降维数据

通过PCA t-sne等方法计算counts数据的降维后的数据，一般是二维或者三维。

###PCA
sce <- scater::logNormCounts(sce)
sce <- scater::runPCA(sce) ###通过PCA算法计算降维数据
sce
# class: SingleCellExperiment 
# dim: 10 3 
# metadata(0):
# assays(2): counts logcounts
# rownames(10): gene_1 gene_2 ... gene_9 gene_10
# rowData names(2): mean detected
# colnames(3): cell_1 cell_2 cell_3
# colData names(16): batch sum ... percent_top_500 total
# reducedDimNames(1): PCA 可以看到这里多了通过PCA算法计算的降维数据
# spikeNames(0):
# altExpNames(0):
reducedDim(sce, "PCA")
#              PC1        PC2
# cell_1 -1.3487780 -0.1053741
# cell_2  0.8693045 -0.4941895
# cell_3  0.4794735  0.5995637
# attr(,"percentVar")
# [1] 82.02114 17.97886
###tsne
sce <- scater::runTSNE(sce, perplexity = 0.1)
# Perplexity should be lower than K!
sce
# class: SingleCellExperiment 
# dim: 10 3 
# metadata(0):
# assays(2): counts logcounts
# rownames(10): gene_1 gene_2 ... gene_9 gene_10
# rowData names(2): mean detected
# colnames(3): cell_1 cell_2 cell_3
# colData names(16): batch sum ... percent_top_500 total
# reducedDimNames(2): PCA TSNE 多了通过TSNE算法计算的降维数据
# spikeNames(0):
# altExpNames(0):
reducedDim(sce, "TSNE")
#             [,1]      [,2]
# cell_1 -5395.251 -1806.471
# cell_2  4278.781 -3769.183
# cell_3  1116.470  5575.653

除了自动添加，也可以手动添加降维数据

u <- uwot::umap(t(logcounts(sce)), n_neighbors = 2)  
###添加通过Uniform Manifold Approximation and Projection (UMAP) method 计算的降维数据
reducedDim(sce, "UMAP_uwot") <- u
reducedDims(sce) # Now stored in the object.

4.6.可选实验

什么叫可选实验？书中给的定义是 data for a distinct set of features but the same set of samples/cells，就是相同细胞或者样本的不同的features，常见的有spike-in RNA

spike-in RNA

spike-in RNA：人工添加的定量的外源RNA（spike-in RNA），它是不存在第一种因素（真实生物差异）的，所以它们的变化可以直接反映技术噪音，然后将整体的内源基因平均表达量变化与spike-in进行拟合。大致了解一下，后面的步骤会仔细讲解。

###建立一个spike_se的class，只含有feature矩阵
spike_counts <- cbind(cell_1 = rpois(5, 10), 
    cell_2 = rpois(5, 10), 
    cell_3 = rpois(5, 30))
rownames(spike_counts) <- paste0("spike_", 1:5)
spike_se <- SummarizedExperiment(list(counts=spike_counts))
spike_se
# class: SummarizedExperiment 
# dim: 5 3 
# metadata(0):
# assays(1): counts
# rownames(5): spike_1 spike_2 spike_3 spike_4 spike_5
# rowData names(0):
# colnames(3): cell_1 cell_2 cell_3
# colData names(0):

###通过altExp函数将spike_se这个class放入sce这个class中去并命名为spike
altExp(sce, "spike") <- spike_se
altExps(sce)
# List of length 1
# names(1): spike 已经成功添加了名为spike的altExps

基本的class部分就含有以上6个部分，还有很多可选项可以根据自己的需求选择（https://osca.bioconductor.org/data-infrastructure.html）

5.scRNA分析流程

单细胞测序

单细胞测序（英语：Single cell sequencing）采取优化的下一代DNA测序技术（NGS）检测单细胞的序列，可以获得特定微环境下的细胞序列差异以方便研究其功能差异等[1]。对个体细胞的DNA测序可以帮助我们了解例如在癌症中的小范围细胞的变异；对其进行RNA测序可以帮助我们了解和鉴别不同的细胞类型与其表现的基因，对研究发育生物学等有较大裨益。[2] 来源自wiki百科

个人理解：就是把原来的组织更近一步拆分成以细胞为单位，更加精准区别遗传信息的异质性，更有助于人们了解不同细胞间的差异。

更详细介绍可以见单细胞测序扫盲：是什么？为什么？怎么做？ 以及公众号单细胞天地

常规的scRNA workflow可以分为一下6个步骤，分别为：实验设计，预处理，导入R，数据处理，统计分析与可视化

5.1 实验设计

开始做单细胞测序之前，需要想好用什么用什么技术去测序，因为现有的测序技术实在是太多了，wiki介绍已经有100多种技术以及被发表（https://en.wikipedia.org/wiki/List_of_single_cell_omics_methods）

1586655684033

所有的scRNA-seq技术可以根据Transcript coverage大致分为以下2种：

1. 获取全长 ，代表技术scSmart-seq2
2. 3端或5端，代表公司10xgenome，旗下技术有Drop-seq，Indrop 以及Chromium

表格来源自文章Single-Cell RNA-Seq Technologies and Related Computational Data Analysis

下面介绍一下这2种代表性技术

代表性seq

代表性seq分为scsmart-seq与10xgenome2种，表格来源于[文章](doi: 10.1093/bfgp/elx035 )

1586530477265

两者之间主要的一点不同mRNA逆转录成cDNA，smart用的是全长，而10x只通过3-tag的方法。而选取哪一种主要看你的目的：主要是做isoform or gene fusion，选全长；如果更关注转录组异质性，选10x（因为会有更多细胞）图片来源自文章Single-Cell RNA-Seq Technologies and Related Computational Data Analysis

不好理解的话可以看哔哩哔哩视频【陈巍学基因】视频56：10X分析单细胞表达;以及【陈巍学基因】视频11：单细胞mRNA测序

An external file that holds a picture, illustration, etc. Object name is fgene-10-00317-g001.jpg

5.2 获取表达矩阵+载入R

表达矩阵

根据书的推荐

• 10X Genomics 数据, 用 CellRanger 软件的pipeline（10x公司的配套工具）去获取表达矩阵+ STAR 比对，具体实战例子可以见单细胞天地系列推文[https://mp.weixin.qq.com/s/fP8f4HboMM7m2Nd7AIljlg]
• 想要更快更有效率，可以用alevin软件的pipeline
• 也可以用scPipe 软件的pipeline+Rsubread比对
• 对于CEL-seq or CEL-seq2 data, 用scruff的pipeline
• 对于全长数据，用传统的bulk流程就可以了

载入R

• 可以用常规的read.table
• 也可以用DropletUtils(for 10X data)

后面几部是sc-seq的重头戏

5.3 质控

为什么要做质控？

因为不同原因（例如细胞损伤，PCR扩增效率等）造成低质量的库（含有较低的counts，较少的表达基因以及高比例的线粒体或者核RNA）会导致下游分析结果不准确：

1. 最常见的是会形成特有的簇（通常是线粒体或者核RNA高比例造成）
2. 破坏种群间的异质性，导致分组不准确（通常是low counts）
3. 正态化中，异常高表达（含有污染的转录本）

怎么做质控？

1. 根据counts的总和（total sum）筛选低表达文库
2. 所有细胞内的基因表达量都不为0
3. 筛除高比例spike-in和线粒体的数据

所有筛选都是通过阈值来进行筛选

大致步骤：获取QC矩阵-根据自己的需求确定阈值-移除低质量数据

5.4 正态化

为什么要做正态化数据？

因为技术原因（例如PCR扩增效率等）造成的数据偏差会导致下游分析结果不准确，所以需要进行正态化校正

正态化等于移除批次效应么？

答案是不，正态化只能考虑技术偏好性，而批次效应还需要考虑不同批次间生物样本的差异

怎么正态化数据？

通过每个细胞的size factor来处理数据

size factor:每个细胞都有自己的偏好性（例如PCR扩增效率不同），通过该细胞的均值比而等比例缩放所有基因。

选取什么软件？

不同的工具有不同的计算size factor的方法：

1. 最简单的就是将所有细胞的基因表达量作为库的total counts，可以用scater包的librarySizeFactors
2. 通过反卷积处理数据，可以用scran包的calculateSumFactors
3. 通过spike-in处理数据，可以用scRNAseq包的computeSpikeFactors

最后通过scater包的logNormCounts来log+正态化数据

5.5 挑选高度差异的基因

highly variable genes (HVGs)

如何挑选高度差异基因？

1. 最简单的方法是选top基因（根据log值与spike-in），采取scran包的modelGeneVar+getTopHVGs
2. 也可以根据变异系数（CV2=ratio of S.D to mean）来选取top基因，采取scran包的modelGeneCV2+getTopHVGs
3. 根据显著性FDR来筛选
4. 还可以结合FDR+top基因的方式来筛选

还有一种是基于先验信息的方法（如已知细胞的亚型）。比如通过SCDE软件鉴定已知不同细胞亚型间的差异表达基因，然后再基于差异表达基因来聚类分析等。

5.6 降维

为什么要降维？

许多scRNA-seq分析都通过细胞在多个基因中的表达量来比较细胞。例如，聚类通过计算跨基因的欧几里德距离来鉴定具有相似转录组表达丰度的细胞。在这些应用中，每个单独的基因代表数据的一个维度。也就是说，如果一个细胞只有2个基因，我们可以通过二维的图来区分，但是一般细胞中都含有多个基因，就会有多个维度，而人是无法看见3维以上的数据。

所以我们需要将多个特征压缩到一个维中，这样既可以节约工作量，又可以进行可视化。

怎么去降维？

虽然对于算法来说，推荐10-50个维度，但是不利于人类理解，所以一般都会降到2维或者3维进行可视化。对于降维，通常选用以下3种算法：

PCA 主成分分析

• 在生物信息学的实际应用情况中，通常是得到了成百上千个基因的信息，这些基因相互之间会有影响，通过主成分分析后，得到有限的几个主成分就可以代表它们的基因了。也就是所谓的降维。更具体的理解见Jimmy老师的推文一文看懂主成分分析。可以用scater的runPCA函数来实现。
• 但是PCA是线性，也就是说每个PC仅捕获高维空间中沿线的变化，所以top2的主成分无法完全展示差异

1586696993027

t-stochastic neighbor embedding t分布随机邻近嵌入

• 也就是常说的tsne,是一种常用的非线性降维方法，非常适用于高维数据降维到2维或者3维，从而进行可视化.可以用scater的runTSNE函数来降维。
• 但是一方面tsne的计算量比较大，另一个是需要调试参数，例如困惑度等参数。还有就是t-SNE不会保留非相邻群集的相对位置，因此我们无法使用它们的位置来确定远程群集之间的关系。

1586697021310

uniform manifold approximation and projection，统一流形逼近与投影

• UMAP，与t-sne相似但是其还尝试找到一种低维表示形式，该维数表示形式保留了高维空间中邻居之间的关系，从而可以反映细胞群体之间分化的连续性和组织性，而且运算速度更快，书中预测是未来做降维的主流。可以用scater的runUMAP函数来进行降维。
• 但是UMAP旨在保留更多的全局结构，但这必然会降低每个视觉集群中的分辨率。就是图不如t-sne清晰

1586697058676

小结：t-sne与UMAP是2种主流降维方法，各有利弊。根据自己的需求做选择。

5.7 聚类

为什么要做聚类？

• 聚类是数据挖掘中的概念，就是按照某个特定标准(如距离)把一个数据集分割成不同的类或簇，使得同一个簇内的数据对象的相似性尽可能大，同时不在同一个簇中的数据对象的差异性也尽可能地大。也即聚类后同一类的数据尽可能聚集到一起，不同类数据尽量分离。具体可以见推文超级干货 ：一文读懂聚类算法
• 其是一种无监督的学习程序，可用于scRNA-seq数据分析以凭经验定义具有相似表达谱的细胞群。会通过人类理解的方式来演示数据，使我们能够通过易于理解的离散标记来描述种群异质性

如何聚类？

聚类不等同于细胞类型，是人类通过自身的经验来聚集的概念，所以实际上不存在“真正的聚类”这个说法。所以需要通过不同的聚类算法去探索最合适的分类。

基于图像的聚类（需要再去理解）

在这里插入图片描述

• 先构建一个图，其中每个节点表示一个cell，并为边赋权重。
• 例如用R包Seurat，根据PCA空间中的欧式距离构造 KNN 图，然后基于其局部邻域中的共享重叠来细化任意两个细胞之间的距离权重。最后根据参数优化模型。
• 书里推荐 scran的pipeline+rank-based weight+igraph包的walktrap算法或者 Seurat 的pipeline +基于Jaccard-based weight+Louvain 进行聚类
• 最后通过clusterModularity() 函数，设定 as.ratio=TRUE，会得到配对集群之间观测到的权重与预期权重之比。最后可视化结果得到log后不同cluster之间的观测到的权重与预期权重之比，颜色越深说明聚类效果越好

1586697205540

K-means聚类

在这里插入图片描述

• k-means算法是非监督聚类最常用的一种方法，因其算法简单和很好的适用于大样本数据

1. 随机地选择k个对象，每个对象初始地代表了一个簇的中心；
2. 将细胞，根据其与各簇中心的距离，将它赋给最近的簇；
3. 重新计算每个簇的平均值，更新为新的簇中心；
4. 不断重复2、3，直到准则函数收敛。具体可以见推文超级干货 ：一文读懂聚类算法

• 缺点：1. k_means算法假设数据是各向同性的，即不同簇类的协方差是相等的，通俗讲就是样本数据落在各个方向的概率是相等的。也就是说假设每个细胞到簇的距离是一致的，是一个等半径的球星簇。但是实际上分区是长短不一，形状不一致的2.它取决于随机选择的初始坐标。这需要多次运行以验证群集是否稳定3. k值必须提前设定，但是设定k值会导致只能聚这么多簇，最终可能会导致一些簇合并在一起。
• k-means通常可结合PCA和t-SNE等来使用
• 书中推荐用cluster包的clusGap与maxSE函数计算最佳的K值，再结合PCA降维与t-sne降维可视化。

1586698621562

• 最后通过群集的群集内平方和within-cluster sum of squares (WCSS) 计算均方根偏差the root-mean-squared deviation (RMSD) ,如果RMSD越低，说明分簇效果越好

1586699915066

层次化聚类

• 又称树聚类算法，是一种古老的技术，旨在生成包含样本层次的树状图透过一种层次架构方式，反复将数据进行分裂或聚合将样本聚集成簇，然后将这些聚簇聚集成更大的簇，依此类推，直到所有样本都属于一个聚簇为止

1. 将每个对象看作一类，计算两两之间的最小距离；
2. 将距离最小的两个类合并成一个新类；
3. 重新计算新类与所有类之间的距离；
4. 重复1、2，直到所有类最后合并成一类。算法详情来自于推文超级干货 ：一文读懂聚类算法

• 缺点很明显，就是运算起来太慢，但如果数据集比较小，适用于这种方法

在这里插入图片描述

• 采用R包的dynamicTreeCut的cutreeDynamic功能进行层次聚类，可以看到可以合适的区分不同的簇

1586700066037

• 通过轮廓宽度silhouette width检查聚类的分离，如果为正数且越大，说明分簇效果越好

1586700537027

subclustering

• 监督的方法。比如基于特定细胞亚型的已知marker基因来聚类分析
• 例如通过记忆T细胞簇证明该簇是根据几种标记物鉴定的。

1586701030900

5.8 寻找marker 基因

why？

• 为了解释前面的聚类结果，我们确定每个簇的marker基因。这些基因使我们能够根据其功能注释为每个簇赋予生物学意义。这些marker大多数情况是我们通过先验实验得到的这些与细胞类型相关的基因。也可选用不同簇差异最为明显的基因作为潜在的marker。

鉴定marker

• 可以通过不同测验方法进行DE检测，然后整合这些差异结果将基因进行排名

配对t测验

实际上我们只需要单边t-test，因为只需要考虑上调的基因做为潜在的marker，通过scran包的findMarkers函数去寻找potential marker。另外，可以调整pval.type参数调整筛选的严格程度。

1586702333552

Wilcoxon秩和检验

• Wilcoxon秩和检验（也称为Wilcoxon-Mann-Whitney检验，或WMW检验）是另一种广泛用于观察组之间成对比较的方法。
• 还是通过scran包的findMarkers函数去寻找potential marker，test参数需要选用Wilcoxon。通过 area-under-the-curve(AUC)判断，<0.5为下调，>0.5为上调，通常选择0.7-0.8之间的基因

1586702730628

二项式检验

• 更为严格的检测办法，检测不同簇细胞中的基因表达比例，主要适用于，基因在一个簇中活跃而在另一个簇中沉默这种情况
• 还是通过scran包的findMarkers函数去寻找potential marker，test参数需要选用binom

1586703276441

• 最后使用multiMarkerStats整合多种检验结果

处理已知且不感兴趣的变异因素

• 例如不同批次以及性别等并不感兴趣的factors，可以通过block以及design去除或者作为唯一因素的design矩阵（类似于DE分析的design矩阵）

5.9 细胞注释（最难的部分它来啦！）

• 获得细胞簇非常简单，但是要确定每个簇代表的生物学状态biological state 很困难，现在大部分是通过肉眼辨别后，根据先验经验去判别，但是这种方法并不适合电脑计算。所以现有以下几种方法去计算细胞标签

通过参考数据库

• 通过机器学习，从注释好的数据中，寻找最相似的参考样本，然后根据相似性，给现有样本分配标签。但是很大的问题就是其正确性取决于之前样本的作者的专业程度。
• 本书推荐用singleR包的SingleR函数去分配标签，通过列子可以看到对单核细胞和B细胞分组是木有异议的，但是对部分NK细胞以及cd4+ CD8+ T细胞不能有效地区分

1586744581864

通过基因集

• 通过识别在每个单个细胞中高度表达的marker所在的基因集来分配细胞标签
• 书中推荐AUCell 包去做找topAUC值或者δAUC，作为细胞标签

通过富集分析

• 将marker基因进行功能注释（limma包的goana功能），根据GO分析富集到的GOterms，去确认是否富集到的terms正确。如果正确，再根据terms分配细胞标签。

5.10 批次效应

移除批次效应

什么是批次效应？

• 大型的单细胞测序项目一般都会产生许多细胞，这些样本制备过程很难保持时间一致、试剂一致，另外上机测序的时候也不一定在同一个测序仪上
• 线性可以用 limma包的removeBatchEffect() 以及 sva包的 comBat() ,前提是假定细胞群体组成在批次中是已知或相同的。
• 书里推荐batchelor包中的rescaleBatches()来移除线性的批次效果，因为相当于对每个基因的对数表达值进行线性回归，并进行一些调整以提高性能和效率。对于每个基因，按比例缩小每个批次中的平均表达，直到等于所有批次中的最低均值。我们选择按比例缩小所有表达式的值，因为当批次位于均值方差趋势的不同位置时，这可以缓解方差差异。
• 另外还有一种MNN校正，它基于高维表达空间中最近邻居（MNN）检测的批量修正策略，不需要预先定义或者已知全部细胞群体组成，它只需要在批次之间有关联的一小部分群体。书里推荐用batchelo包的 fastMNN()来进行MNN校正