一、
背景
单细胞转录组测序（Single cell RNA sequencing，scRNA-seq）是在单细胞水平对转录组进行测序的新技术。它能研究单个细胞内的基因表达情况，解决用组织样本测序无法解决的细胞异质性难题，并实现以下目的：
-
-
-
-
识别在特定条件下（例如治疗或疾病）特定细胞类型中差异表达的基因
-
在整合空间、调控和/或蛋白质信息的同时，探索细胞类型中表达的变化
（Mol Syst Biol. 2019. Current best practices in single‐cell RNA‐seq analysis: a tutorial）
虽然完整的单细胞数据分析流程应该是CellRanger+Seurat，其中CellRanger是由10x genomic公司提供的专门用于其单细胞转录组数据分析的软件包，能够将fastq测序数据比对到参考基因组上，然后进行基因表达定量，生成细胞-基因表达矩阵。
由于GEO数据库上最便捷的数据已经是CellRanger处理后的数据（filtered_barcodes_bc_matrices），所以本文先不介绍CellRanger的下载和使用方法（我们以后会涉及到），而是介绍从下载完GEO数据之后进行分析的基础流程（下图中QC及之后的内容）。
(https://github.com/hbctraining/scRNA-seq/blob/master/lessons/03_SC_quality_control-setup.md）
1.生成计数矩阵：格式化读取数据，映射到参考基因组上并定量。
3.对过滤后的计数进行聚类：基于转录组的相似性对细胞进行聚类（细胞类型 = 不同的聚类）。
（1.属于CellRanger部分的内容，本文主要介绍的是2.3.4.部分）
本文使用的是10X Genomics公司提供的外周血单个核细胞(PBMC)数据集。通过使用 Illumina NextSeq 500测序得到了2700个细胞的单细胞转录组数据。
测试数据可以从Seurat官网
https://satijalab.org/seurat/
获得。
（Seurat官网：https://satijalab.org/seurat/）
前面说到，Seurat的分析是从
单细胞计数矩阵
开始的。单细胞计数矩阵是由
cellranger
流程对下机数据进行处理而产生的，里面存储的是UMI计数矩阵，矩阵中的值表示的是
在每个细胞（列）中检测到的每个特征（即基因，行）的分子数
。
00 加载分析所需要的包
01 构建Seurat对象
（
1
）
加载数据集
测试数据集中包含barcodes.tsv、genes.tsv、matrix.mtx三个文件。
在读取10×这个平台文件的时候，不管如何命名，都需要把三个文件整理成如下的形式。
1.barcodes.tsv存储的是用于
标记不同细胞
的条形码序列，即barcode。该文件的行数对应的就是
细胞数量
，由行数统计结果可知细胞数量为2700个；
2.genes.tsv存储的是
基因
的信息，分别为gene id和gene symbol。与barcodes.tsv文件同理，该文件的行数对应的是
基因数量
，基因数量为32738个；
3.matrix.mtx存储的主要是
有表达量
的基因的counts数量，前两行是注释信息不用关注。第三行跟我们上面统计的每个文件的行数对应，这三个数字对应的便分别是
基因数量、细胞数量、有表达量的基因的数量
（任一基因在任一细胞中出现都算一次计数），去掉前三行，刚好是2286884行。
这种排列方式可以高效地存储数据，例如“32704 1 1”这一行，代表genes.tsv文件中的第
32706
个基因在barcodes.tsv文件中的第
1
个细胞中的counts数为
1
。
（
2
）
创建Seurat对象
CreateSeuratObject的几个主要参数：
-
-
-
min.cells：规定基因的表达范围，即每个基因至少在多少个细胞中表达方可保留
-
min.features：规定细胞表达的基因数范围，即每个细胞至少表达多少个基因方可保留
可以看到在经过过滤以后，数据集变成了一个2700个细胞，13714个基因的矩阵。
然后我们来看一下我们构建的Seurat对象是什么样子的（其中assays存储的为表达矩阵；meta.data存储的为细胞信息；active.ident存储的为表达矩阵细胞名；active.assay存储的为表达矩阵的名）
02 标准预处理
-
nFeature_RNA：
每个细胞中检测到的基因的数量
。
低质量的细胞或空液滴通常只有很低的计数；
细胞多胞体可能表现出异常高的计数。
-
nCount_RNA：
细胞内检测到的分子总数
。即检测到的基因的总数量（一个基因可能会因为有多个转录本而被检测到多次），因此分子总数与基因数量密切相关。
（上面这两个指标由于在CreateSeuratObject()过程中会自动计算基因的数量nFeature_RNA和分子总数nCount_RNA，因此在创建对象的时候，他们就已经被存储在对象的meta.data中了。）
-
percent.mt：
对应线粒体基因组的序列的百分比
。低质量/濒死的细胞通常表现出广泛的线粒体污染，当然也有可能是特殊细胞（心肌细胞或者神经细胞）。我们通常使用PercentageFeatureSet()函数来计算线粒体基因相关的reads的百分比。代码如下：
"[[]]"运算符可以向对象pbmc的元数据meta.data中添加列。meta.data是一个存储质控数据的好地方。
可以看到多出了一列“percent.mt”数据，表示每个细胞当中线粒体基因的百分比。
（
1
）
设置质控阈值
质控的目的是挑选质量更好的细胞进行进一步的分析。质控标准的方法很简单，小提琴图VlnPlot和散点图FeatureScatter ，我们就按照官网上的来：
-
第一张图为所有细胞检测到的基因的表达数量，横坐标为样品名，纵坐标为每个细胞中包含的基因数量，其中一个点代表一个细胞；
-
第二张图为所有细胞检测到的RNA表达的数量，横坐标为样品名，纵坐标为每个细胞中的分子总数；
-
-
第一张图，是样本中RNA的表达数量（分子总数）和线粒体基因比例的关系，通常两者之间没什么关系，这个才是正常的样本。
-
第二章图，基因的数量和RNA的数量（分子总数）呈正相关，即测序深度和基因数量之间的关系。只要排除右上方的离群点，也就是由于双胞体或多胞体造成的分子数量过大的那部分数据。
抛开教程说一句，虽然上面两张图可以给我们提供制订质控标准的依据，但完全凭借个人经验去设置标准风险还是太大了，质控不严容易被诟病，质控太严格又容易丢失一些关键信息。最稳妥的办法还是参考研究领域内相似实验设计的大文章，其Methods部分往往会有单细胞数据分析的细节内容，在此基础上制订质控标准效率会更高。如下：
（Cancer Discovery. 2023. Stromal reprogramming by FAK inhibition overcomes radiation resistance to allow for immune priming and response to checkpoint blockade)
（Cancer Cell. 2023. Fungal mycobiome drives IL-33 secretion and type 2 immunity in pancreatic cancer）
（
2
）
过滤数据
为了和官方教程相统一，这一次我们还是采用官方的过滤阈值
过滤之后，细胞数量由2700减少到了2638，基因数量不变。
03 数据标准化
在这一步默认使用的是LogNormalize，是一种全局缩放的标准化方法（当然也有其他的标准化方法，
尤其是SCTtransform
，我们以后会慢慢涉及）。数据标准化可以排除由于技术误差带来的异质性，更好地展现细胞自身基因的表达差异。
标准化之后的结果会被储存在pbmc[["RNA"]]$data中
可以看到，原始数据中表达量为1的数据在标准化之后发生了改变（“.”表示表达量为0，用这样的表示方法可以节省空间）。
04 鉴定高可变基因
接下来要寻找在不同细胞间有显著差异的基因子集（即它们在一些细胞中高表达，在另一些细胞中低表达）。这样可以节省计算资源，降低计算成本。默认使用前2000个基因进行下游的分析。
05 数据缩放
数据缩放是进行降维处理之前的一个标准步骤，作用如下：
-
调整每个基因的表达量，使得细胞间的平均表达量为0；
-
-
目的：这一步骤赋予每个基因在下游分析中同样的权重，不至于使高表达基因占据主导地位；
-
默认情况下，只对高可变基因进行缩放，可以通过指定features来对其他基因进行缩放。
06 PCA降维
接下来我们对缩放数据执行 PCA线性降维。默认情况下，只对先前确定的高可变基因进行处理。
关于PCA结果的展示，官网其实给了好几种方法，包括VizDimReduce()、 DimPlot()和DimHeatmap()。效果如下：
虽然感觉文章一下子可以多好几张图，但实际上这些图对于我们后续的分析帮助并不大，很多高分的文献也不会把这些结果放到文章里，哪怕是副图。真正会用到的是下面的步骤。
07 确定PCA维度
-
在scRNA-seq数据中，每个单独的特征（基因表达）可能包含大量的技术噪声。为了克服这个问题，Seurat（一种流行的单细胞数据分析工具）使用主成分分析（PCA）来聚类细胞。PCA是一种降维技术，它将原始数据转换到一个新的坐标系统，这个坐标系统由主成分（PCs）组成，每个主成分都是原始特征的线性组合。
-
在这个过程中，每个主成分可以被视为一个“元特征”（metafeature），它结合了一组相关基因。这样，PCA分数（即每个细胞在各个主成分上的坐标）可以帮助我们理解细胞之间的差异，而不仅仅是单个基因的表达。
这边难以理解没有关系，我们只需要知道，确定主成分数量的意义就是选择合适的主成分数量来捕捉数据中的关键信息，同时避免过多的技术噪声。PCA太少，会丢失数据中的一些信息；PCA太多，又会引入很多琐碎的技术干扰。
此处官网提供了一种名为肘形图（拐点图）ElbowPlot()的方法：
这个图我们主要看的是拐点，个人理解就是在PC为10的时候图形趋于平稳，所以这里我们选择PC为10作为我们后续分析的一个参数。
当然官网鼓励我们使用不同的PCs进行多次尝试（向更多PCs去尝试），其实差别不会特别大。
08 细胞聚类及其可视化
终于到了能出图的环节了，很多文章单细胞分析的第一张图甚至很多文章整篇的第一个图就是细胞聚类及可视化结果。
细胞聚类就是把一个整体当中的细胞分为多个亚群，在具有相似表达模式的细胞之间绘制边界，然后划分为高度相关的不同簇（cluster）。
官网推荐分辨率从0.4-1.2都可以尝试，我们可以多次调整分辨率，（在你已经对组织的细胞数量及比例有一个心理预期以后）保证每一个细胞亚群不会太大或者太小。
聚类的可视化方法有两种，都属于
非线性降维
，一种是UMAP，一种是tSNE，两者没有孰好孰坏的区别，在文章当中甚至可以全都放上去。
09 寻找差异表达基因
每一个cluster都会有属于自己的高表达基因和低表达基因，Seurat可以通过差异表达(DE)的基因来定义这些cluster。
差异表达基因可视化的方法也有很多，也许以后我们会专门讲一下这部分（待定），在这里先简单介绍几种方法。
这部分的结果可以用来在一些在线网站或者软件当中辅助细胞注释，所以要保存下来。
同时，由于细胞注释的结果往往不会尽如人意，需要多次的修改，且跟上面的这些流程都相对独立，所以也停在这里保存一下RDS文件。
10 细胞注释
