scRNA-seq 教程

本教程，我们将为大家介绍如何分析单细胞测序数据。高通量单细胞RNA-Seq产品依托于Droplet技术（液滴法），可实现上万个单细胞的同时进行分离和建库，大幅的降低建库成本的同时实现全局性的组织细胞转录图谱。

单细胞目前面临的挑战主要是：扩增量目前最多是一百万个细胞；有可能一个基因在一个细胞中能检测到中等表达量，但是在另一个细胞中却检测不到，这种现象叫做"gene dropouts"。

方法主要分为两大部分：定量与捕获。

定量包括两种类型：全长以及基于标签(tag)。前者对每个转录本都试图获得一致的read覆盖度，后者只捕获5‘或者3’端的RNA。Smart-Seq方法以其可实现完整mRNA扩增的优势，被广泛的应用。

捕获的技术决定了细胞如何被筛选、获取怎样的测序外的补充信息、数据产量，三种最常用的方法是基于微孔microwell-(流式细胞术FACS分选到孔板)、微液流microfluidic-(Fluidigm's C1)、微滴droplet-(Drop seq, 10x genomics)

TED视频: Drop seq to study brains, which is a data driven discovery.

Human Cell Atlas一项旨在为人体中每种细胞类型创建详细分类的全球性计划。

[推荐] 单细胞技术的历史、分类、优劣扫盲。

现在开始学习单细胞测序分析常用的R包 Seurat: QC以及后续分析探索数据

现在开始学习单细胞测序分析常用的R包 monocle: 伪时间分析 Monocle performs differential expression and time-series analysis for single-cell expression experiments.

现在开始学习单细胞测序分析常用的R包 SingleR: 全自动细胞类型注释R包

最佳实践

paper: 2019 MSB paper

基础扫盲

【实用技能】单细胞测序平台有哪些？

单细胞RNA-seq介绍

10x scRNA-seq

10X genomics资料: 10x 官方数据集 | 10x 文章 | 单细胞技术扫盲 |

教程 Analysis of single cell RNA-seq data: 一步一步的操作步骤。

公众号教程: 基于Seurat的基础分析、以及基于clusterProfiler、Monocle、SingleR等R包的延伸分析。

单细胞课程： 初级8讲和高级分析8讲, 单细胞分析十八般武艺, 单细胞分析环境,
4细胞类型鉴定: CellMarker库 /PanglaoDB库， SingleR包;
5细胞再聚类: 细胞子集的提取使用subset函数，主要依据meta.data的分类项目选择，也可以自定义细胞barcode提取。
6伪时间分析: Monocle假设轨迹是树状结构，一端是“根”，另一端是“叶”。Monocle2使用DDTree降维图，Monocle3使用UMAP降维图。
7差异基因富集分析: seurat聚类分析得到的按cluster分类，singleR分析得到的按细胞类型分类，monocle分析得到的按拟时状态（state）分类。不同的细胞类型之间，有哪些表达差异基因呢，这些差异基因有特别的意义吗？
8可视化工具总结: RidgePlot山脊图 /VInPlot小提琴图 /FeaturePlot特征图 /DotPlot点图 /DoHeatmap热图 /FeatureScatter散点图 /DimPlot降维图
公众号专题：单细胞进阶

谁适合阅读？

本教程介绍的是高通量测序分析知识，适合从事高通量测序数据分析和数据可视化的人员阅读。

需要具备的知识？

如果你熟悉主流计算机编程语言，有一定的Linux、R和python基础。

生物信息与R语言QQ群: 187923577

表达矩阵

scRNA-seq数据的分析以其表达矩阵为起点，表达矩阵中的每一行代表一个基因，每列代表一个细胞。每个条目代表给定细胞中特定基因的表达水平。

1、表达矩阵中表达值的数据类型

counts: Raw count data, e.g., number of reads or transcripts for a particular gene.
normcounts: Normalized values on the same scale as the original counts. For example, counts divided by cell-specific size factors that are centred at unity.
logcounts: Log-transformed counts or count-like values. In most cases, this will be defined as log-transformed normcounts, e.g., using log base 2 and a pseudo-count of 1.
cpm: Counts-per-million. This is the read count for each gene in each cell, divided by the library size of each cell in millions.
tpm: Transcripts-per-million. This is the number of transcripts for each gene in each cell, divided by the total number of transcripts in that cell (in millions).

一般取 apply(df, 2, function(x){ log2( x/sum(x) +1 } )

2. 独特的分子标识符（UMI）

unique Molecular identifiers 是在逆转录过程中添加到转录本中的短随机条码（4-10bp）。它们可以将测序读序列分配给单个转录分子，从而消除scRNASeq数据中的放大噪声和偏差。当对包含UMI的数据进行测序时，技术通常只对包含UMI的转录本的末端（通常是3’端）进行特定的测序。

3. 过滤矩阵

去除在任何细胞中都不表达的基因；
具有少量reads/molecules的细胞，很可能是已经被破坏或捕获细胞失败，应该移除这类细胞；
手动过滤细胞：大多数细胞检测到的基因在7000-10000之间，这对于高深度scRNAseq来说是正常的。但是这也取决于实验protocol和测序深度。比如说基于droplet或其他测序深度较低的方法，通常其每个细胞检测到的基因较少。如果细胞检出率相等，则分布应近似正常，因此我们移除分布尾部的那些细胞（检测到的基因少于7000的细胞）；
自动过滤细胞：自动异常检测来识别可能存在问题的细胞：scater包创建一个矩阵（行代表cell，列代表不同的QC度量），然后PAC提供了按照QC度量排序的cells的2D表示，然后使用来自mvoutlier包的方法检测异常值。
过滤基因：移除表达水平被认为是undetectable的基因。We define a gene as detectable if at least two cells contain more than 1 transcript from the gene. If we were considering read counts rather than UMI counts a reasonable threshold is to require at least five reads in at least two cells. 然而很多时候阈值取决于测序深度，而且很重要的一点是基因过滤一定要在细胞过滤之后，因为一些基因可能仅仅被检测到只存在在低质量的细胞里。

处理混杂因素

技术上的混杂因素包括识别和移除表达数据中与感兴趣的生物学信号无关（即混杂因素）的变异源。目前存在各种各样的方法，其中一些使用spike-in or housekeeping genes,一些使用endogenous genes。

数据可视化

1. PCA

PC对应于协方差矩阵的特征向量，特征向量按照特征值排序，因此第一主成分尽可能的考虑数据的可变性，并且每个后续成分在与前面的成分正交的约束下具有最高的方差。

注意：对数转换对数据是有益的：减少了第一主成分的方差，并且已经分离了一些生物效应。而且使表达值的分布更正常。

2. tSNE（tSNE已死，请用 UMAP）

tSNE（t-Distributed Stochastic Neighbor Embedding）将维数降低与最近邻网络上的随机游走相结合，将高维数据（例如14214维表达矩阵）映射到二维空间，同时保持细胞间的局部距离。与PCA相比，tSNE是一种随机算法，在统一数据集上运行多次该方法会得到不同的图。由于tSNE的非线性和随机性，导致很难解释。为了确保再现性（reproducibility）,我们可以固定随机数产生器的种子，以保证总是得到相同的图。

此外，tSNE需要提供一个perplexity值，该值反应了用于构建最近邻网络的邻居数量。该值越高表明产生了一个密集网络，将细胞聚集在一起；该值越低表明网路更加稀疏，允许细胞群彼此分离。在scater包中，perplexity的默认值是：cells总数除以5。

3. UMAP

2018年McInnes提出了算法，UMAP（Uniform Manifold Approximation and Projection for Dimension Reduction，一致的流形逼近和投影以进行降维）。一致的流形近似和投影（UMAP）是一种降维技术，类似于t-SNE，可用于可视化，但也可用于一般的非线性降维。

UMAP应该说是目前最好的降维算法了，能最大程度的保留原始数据的特征同时大幅度的降低特征维数。处理单细胞数据时，UMAP是显著优于tSNE的。样本量大时 tSNE/FltSNE 很慢且太消耗内存。tSNE没有保留全局结构，也就是组内的相似，但是组间的距离远并不反映相似程度的小。另外，tSNE只能嵌入到2维或3维中。

论文：McInnes, L, Healy, J, UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction, ArXiv e-prints 1802.03426, 2018
https://pypi.org/project/umap-learn/

 #手写数字聚类，每个数字 64 维。
from sklearn.datasets import load_digits
import matplotlib.pyplot as plt

digits = load_digits()
fig, ax_array = plt.subplots(20, 20)
axes = ax_array.flatten()
for i, ax in enumerate(axes):
    ax.imshow(digits.images[i], cmap='gray_r')
plt.setp(axes, xticks=[], yticks=[], frame_on=False)
plt.tight_layout(h_pad=0.5, w_pad=0.01)
plt.show() # 显示手写数字图像


# 使用umap降至2维并绘制散点图
import umap
import numpy as np

reducer = umap.UMAP(random_state=42)
embedding = reducer.fit_transform(digits.data)
print(embedding.shape)
# 可视化
plt.scatter(embedding[:, 0], embedding[:, 1], c=digits.target, cmap='Spectral', s=5)
plt.gca().set_aspect('equal', 'datalim')
plt.colorbar(boundaries=np.arange(11)-0.5).set_ticks(np.arange(10))
plt.title('UMAP projection of the Digits dataset')
plt.show() 


# 更多参数
embedding = umap.UMAP(n_neighbors=5,
                      min_dist=0.3,
                      metric='correlation', random_state=8).fit_transform(digits.data)
# 可视化
plt.scatter(embedding[:, 0], embedding[:, 1], c=digits.target, cmap='Spectral', s=5)
plt.gca().set_aspect('equal', 'datalim')
plt.colorbar(boundaries=np.arange(11)-0.5).set_ticks(np.arange(10))
plt.title('UMAP projection of the Digits dataset')
plt.show() 


附录: install UMAP on Ubuntu 20.04
$ sudo pip3 install -U llvmlite -i https://pypi.douban.com/simple/
$ sudo pip3 install -U numba -i https://pypi.douban.com/simple/

$ pip3 uninstall llvmlite
$ pip3 install -U --ignore-installed numba

$ pip install umap-learn
ref: https://zhuanlan.zhihu.com/p/109584077