# 1.2 总览 | 从实验到分析

## 一张非常重要的图

下面的这两张可以说是贯穿始终

![scRNA-Workflow](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-10-30-080845.png)

![贯穿始终的重点](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-10-26-094429.png)

## 实验设计环节

在正式分析之前，关于实验问题的探讨是很有必要的，**最重要的一个就是技术的选择：**

* **Droplet-based**: 10X Genomics, inDrop, Drop-seq
* **Plate-based** with unique molecular identifiers (**UMIs**): CEL-seq, MARS-seq
* Plate-based with **reads**: Smart-seq2
* Other: sci-RNA-seq, Seq-Well

每种方法都有优劣（[Mereu et al. 2019](https://pubmed.ncbi.nlm.nih.gov/32518403/); [Ziegenhain et al. 2017](https://pubmed.ncbi.nlm.nih.gov/28212749/))，目前以10X为代表的droplet-based方法由于高通量和低细胞成本成为了约定俗成的技术；Plate-based方法可以捕获其他的一些表型信息（如细胞形态），另外可以根据实验目的进行调整，灵活性比较好；Read-based方法可以覆盖全转录本，在分析可变剪切、外显子突变等方面很有用；UMI-based方法可以减轻PCR扩增偏差。

下图来自文章的评测：[Benchmarking single-cell RNA-sequencing protocols for cell atlas projects](https://pubmed.ncbi.nlm.nih.gov/32518403/)

![文章的评测](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-25-044413.png)

### **下一个问题就是：到底要捕获多少细胞？测序要测多深？**

答案言简意赅：**As much as you can afford to spend.**

如果再补充一下这个答案就是：想要发现罕见细胞类群，就要多获得细胞；想要探索潜在的微小差异，就要加大测序深度。目前常用的droplet-based仪器可以捕获1万到10万细胞，测序深度是每个细胞1000到10000 UMIs，在经济条件一定的前提下，它们之间一般是成反比。另外它还要权衡高细胞捕获通量和影响捕获效率的“双细胞（doublets）比例”。

### **实验设计和常规转录组类似**

也是要考虑一个实验条件下多个生物重复，而且实验条件最好不要混杂批次。需要注意的是：生物重复不是指的单个细胞，而是指的提供细胞的供体（donors）或者细胞培养体系（cultures）

## 获得表达矩阵（count matrix）

和常规转录组一样，单细胞转录组也是需要得到表达矩阵，才能进行下游分析。表达矩阵包含的信息就是：每个细胞中比对到每个基因的UMIs或者reads数。有一点需要注意：**它的定量方法和具体的实验技术相关**

* 10X的数据：使用`CellRanger` 软件，基于STAR比对到参考基因组，然后统计每个基因的UMIs数量
* Pseudo-alignment方法（如`alevin`）：就像之前用的salmon、kallisto意思一样，不需要比对参考基因组，节省时间、内存
* 对于一些高度multiplexed的方法：可以使用[scPipe](http://www.bioconductor.org/packages/release/bioc/vignettes/scPipe/inst/doc/scPipe_tutorial.html) 包：提供了一套综合的分析流程，利用`Rsubread`比对，然后统计每个基因的UMIs数量

  > multiplexed：翻译叫做”多路复用“，即：large numbers of libraries to be pooled and sequenced simultaneously during a single run，可以节省成本和时间
* CEL-seq、CEL-seq2数据：[scruff ](https://bioconductor.org/packages/release/bioc/vignettes/scruff/inst/doc/scruff.html)包可以专门分析
* read-based方法：可以使用常规bulk 转录组定量的流程（比如smartseq2就可以用hisat2+featureCounts）
* 任何包含spike-in转录本的数据：spike-in序列都要在比对、定量之前加到参考基因组中

定量结束后，一般是先导入表达矩阵然后创建一个`SingleCellExperiment`对象（例如：`read.table() + SingleCellExperiment()`）。除此以外，还有一些特定的文件格式需要用特定的包，比如[DropletUtils](https://bioconductor.org/packages/release/bioc/vignettes/DropletUtils/inst/doc/DropletUtils.html)可以分析10X数据，[tximport](https://bioconductor.org/packages/release/bioc/vignettes/tximport/inst/doc/tximport.html)/[tximeta ](https://bioconductor.org/packages/release/bioc/vignettes/tximeta/inst/doc/tximeta.html)可以分析pseudo-alignment数据

### **需要注意**

* 如果分析的是人类数据并且加入了ERCC，**我们很多时候直接用`^ERCC`在行名中进行正则匹配，但是这时要小心**，因为**ERCC基因家族在人类基因组注释中确实存在**，很有可能将真的基因作为外源转录本进行分析。这个问题可以通过将表达矩阵的行名设置为Ensembl,或Entrez来解决
* 一些定量工具会统计表达矩阵中的reads比对率，会存在一些未必对的情况。尽管这些信息可以用作质控，但这些数值如果被误认为是表达量信息，那么就会干扰下游分析。因此在进行下游分析之前，这部分信息可以去掉或者保存在`colData`中

## 数据处理与下游分析

1. **首先进行质控**：去掉低质量细胞。这些细胞可能在建库环节被破坏，可能没有被有效捕获（这就是所谓的“dropout”）。一般会统计：每个细胞的全部count数、spike-in或线粒体reads比例、检测到基因的数量
2. **表达矩阵归一化**：为了减小细胞文库的偏差（可能由于细胞捕获效率不同、测序深度的差异而造成文库大小差异），把细胞们放在同一起跑线上，才能进行下面的细胞相似性比较，后面再根据相似性进行细胞分群。一般是基于log转换（当然有的函数也涉及了一些size factor的计算），从而对均值-方差进行校正
3. **挑选一些特征基因**（一般是高变化基因HVGs，Highly Variable Genes）进行下游分析。原理是根据每个基因在细胞之间的差异构建变化模型，然后找那些变化差异大的基因。使用HVGs不用全部基因的原因一是为了减少计算量，二是减少不感兴趣基因（比如在细胞之间没什么差异）对分析产生的噪音
4. **降维处理**：让数据更“紧凑”，一般是线性降维PCA+非线性降维tSNE/umap。PCA一般是先获得初步的低维数据（可能会挑出几十个主成分），然后传给t-SNE进一步压缩，进行可视化
5. **细胞聚类：**根据细胞归一化后的表达量相似性分成组，然后根据每个组marker基因（可理解为这一群细胞的标志性基因）的差异表达对分群进行生物学定义

比如用来自[scRNAseq](https://bioconductor.org/packages/release/data/experiment/vignettes/scRNAseq/inst/doc/scRNAseq.html)的一个droplet-based的视网膜数据【[Macosko et al. (2015)](https://pubmed.ncbi.nlm.nih.gov/26000488/)】，就从原始矩阵得到了分群结果，可以看到这里不使用Seurat也能做质控、挑高变化基因等等

```r
library(scRNAseq)
sce <- MacoskoRetinaData()

# 质控
library(scater)
is.mito <- grepl("^MT-", rownames(sce))
qcstats <- perCellQCMetrics(sce, subsets=list(Mito=is.mito))
filtered <- quickPerCellQC(qcstats, percent_subsets="subsets_Mito_percent")
sce <- sce[, !filtered$discard]

# 归一化
sce <- logNormCounts(sce)

# 挑高变化基因
library(scran)
dec <- modelGeneVar(sce)
hvg <- getTopHVGs(dec, prop=0.1)

# 降维
set.seed(1234)
sce <- runPCA(sce, ncomponents=25, subset_row=hvg)
sce <- runUMAP(sce, dimred = 'PCA', external_neighbors=TRUE)

# 聚类
g <- buildSNNGraph(sce, use.dimred = 'PCA')
sce$clusters <- factor(igraph::cluster_louvain(g)$membership)

# 可视化
plotUMAP(sce, colour_by="clusters")
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-10-30-090848.png)

最后注意这里的分群并不一定是真正有生物学意义的，根据不同的参数可以得到不同的分群结果，而且这里看到的多个小群也有可能是同属一个大群。最后的分群需要计算+生物知识共同实现。


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://jieandze1314.osca.top/01/01-2.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.