单细胞测序基础知识

最新推荐文章于 2024-08-15 02:16:41 发布

zhengxj_

最新推荐文章于 2024-08-15 02:16:41 发布

阅读量4.9k

点赞数 6

分类专栏：转载文章标签：人工智能

原文链接：https://mp.weixin.qq.com/s/_xs4NagBiJ44S7p_-TVaKw

版权

转载专栏收录该内容

2 篇文章 0 订阅

订阅专栏

单细胞测序基础知识

一、单细胞测序

学习scRNA-seq先了解一下Bulk RNA-seq

Bulk RNA-seq：测量一个大的细胞群体中每一个基因的平均表达水平，对比较转录组学（例如比较不同物种的相同组织）是有帮助的，但对于研究异质性系统（例如，复杂的组织如脑）还是不够的，对于基因表达的本质研究不够深入。
scRNA-seq：
- 2009年由Tang首创，随着测序技术的发展及测序成本的的下降，2014年才逐渐流行开来。它测定的是细胞种群中每个基因的表达量分布，对于研究特定细胞转录组的变化是重要的。
- 通量：测定的细胞量也由最初的10个左右，达到了现在的百万并且还在不断递增，向着高通量的方向发展。
- 测序流程：现在主流的主要10X Genomics Chromium（较多细胞），SAMRT-seq2（较多基因）和Fluidigm C1等。当然还有其他的如：CELL-seq、Drop-seq、mas-seq和Wafergen ICELL8等。
  图片

测序方法

二、单细胞测序流程

测序流程

三、方法主要分为两大部分：定量与分离。

单细胞定量

包括两种类型：全长以及基于标签(tag)。前者对每个转录本都试图获得一致的read覆盖度，后者只捕获5‘或者3’端的RNA。定量方法的选择也影响了后续分析的方法选择。理论上，全长的方法应该得到转录本的平均覆盖度，但是实际上，覆盖度经常是有偏差的。基于标签的方法能够利用特异性分子标记(Unique Molecular Identifiers, UMIs)提高定量准确度，但是呢，这种限制了转录组一端的方法有降低了转录本的可拼接性，让以后的isoform识别变得困难。
单细胞分离（Isolation of Single Cells）
单细胞分离方法：包括有限稀释法（Limiting dilution technique）显微操作法（micromanipulation）、流式细胞分选(FACS)、激光捕获显微分离(LCM)、微孔（microwell）、微流控（microfluidics）和微滴（droplet）。

有限稀释技术（Limiting dilution technique）是利用移液管稀释分离细胞，这种方法的主要缺点是效率低下，成功率20%左右。

显微操作法（micromanipulation）是一种经典的方法，用于从少量细胞样本中提取细胞，如早期胚胎或未培养的微生物，缺点是耗时和低通量的。

流式细胞分选(FACS)广泛用于分离单个细胞，在悬浮状态下需要较大的起始体积(>10,000个细胞)。

激光捕获显微分离(LCM)是利用计算机辅助的激光系统将单个细胞从固体组织中分离出来。

微孔microwell 优势是可以与荧光触发细胞分离技术(FACS)结合，实现基于表面标记的细胞选择。对于想要分离特定细胞群的研究很有效；另外一个优势就是，可以对细胞进行拍照，帮助确定孔中是否存在损伤的或者重复的细胞。缺点是通量低，对每个细胞进行操作需要很大的工作量。

微流控（microfluidics）相对微孔，它的通量更高。缺点是只有10%的细胞能被捕获到，加入处理的细胞类型比较稀少，那么能被芯片捕获的就更少了，通量不高，成本高。

微液控制技术 例如Fluidigm's C1

微滴技术 是将单个细胞包裹在µl级别的液滴中，液滴被搭载到建库所用的酶上，每个微滴包含一个独特的条码(barcode)，由那个被包装好的细胞产生的所有reads都被贴上了该条码，也是为了之后对于不同细胞reads的分辨。一般微滴技术的通量最大，查资料得知一秒能包装7万个以上的液滴，并且建库费用合适--0.05美元/细胞。但是高额测序费用又成了他的短板，鉴定出的一般只有一千多个差异转录本，覆盖度还是比较低的

这些单细胞分离方案具有各自的优点，可以根据研究目的来选择合适的分离方法。

四、分析流程

黄色部分对于高通量数据的处理都是差不多的流程；
橙色的部分需要整合多个转录组分析流程以及显著性分析，来解决单细胞测序的技术误差；
蓝色是下游表达量、通路、互作网络等分析，需要使用针对单细胞研发的方法。

五、应用

单细胞技术揭示细胞的异质性、细胞分化、发育过程中不同的生理状态下的变化、疾病发生病变等生物学过程。

单细胞转录组测序主要应用方向

1.大规模细胞图谱构建

特定组织裂解后通过单细胞测序获得单细胞转录组图谱，并基于每个细胞基因表达谱数据进行细胞类型聚类，分析研究复杂器官中不同细胞亚型的功能，了解细胞间的差异以及各种细胞群体间的协作关系。

2.细胞亚群细化&稀有细胞类型鉴定

在单细胞类型聚类基础上，依照已知细胞类型标志基因表达情况和新基因表达情况，进一步细分细胞类型并发现新的细胞亚群，分析各个小亚群细胞差异，小亚群细胞和稀有细胞类型在生物学过程参与的功能。

3.肿瘤方向研究

肿瘤细胞异质性
肿瘤微环境
肿瘤耐药性
肿瘤干细胞分化

4.干细胞发育及分化研究

通过单细胞测序全面解析干细胞异质性，鉴定干细胞分化不同阶段的不同类型分化过程细胞，鉴定不同亚型干细胞Marker基因，绘制干细胞发育和分化轨迹，解析干细胞发育和分化的分子机理。

5.免疫方向研究

构建免疫图谱
构建血液或肿瘤微环境、特定组织器官或疾病状态下免疫细胞构成，分析免疫细胞异质性，解析机体复杂的免疫机制；发现新的免疫细胞Marker基因，鉴定稀有免疫细胞类型。
免疫细胞分化
鉴定不同疾病状态、诱导条件下免疫细胞基因表达变化，鉴定不同分化方向免疫细胞Marker基因，绘制免疫细胞细胞发育和分化轨迹

6.神经科学研究

研究神经元细胞分群情况以及不同神经元细胞类型的异质性，了解神经干细胞-神经元细胞分化过程的分子调控机制。

7.发育生物学

脑发育研究

可以获得复杂的脑细胞基因表达图谱，构建脑细胞群体，解析脑细胞的异质性以及不同亚群细胞的特定功能。
胚胎细胞发育研究

8.生物标志物/疾病分型

发现并鉴定异常增殖的细胞类型，结合传统病理学特征，辅助疾病分型

单细胞转录组数据处理之上游分析流程

scRNA-seq技术到目前为止也有一百多个了，但主流的可以大致分为以下几种：

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

实际上你需要理解的就是10x数据和Smart-seq2技术啦，最常用而且最常见！上游分析流程我们分开讲解，在群主的7个小时的单细胞转录组视频课程（限时免费）视频里面演示的其实是Smart-seq2技术的单细胞转录组数据处理，而且仅仅是半个小时的教学，其实是需要你有非常多的背景知识才可能看得懂。

首先针对10x数据

目前单细胞是10x的天下，而10x的测序数据，御用软件cellranger其实就是star的包装，关于10X仪器的单细胞转录组数据走cellranger流程，我们在单细胞天地多次分享过流程笔记，大家可以自行前往学习，如下：

而且10X的单细胞转录组数据的文章都已经快1000篇了，但凡你认真一点看看文献也基本上可以看到上游数据分析描述（文献描述肯定是非常简洁的啦）

10X的单细胞转录组数据处理文章描述

关键是要搞清楚你的输出和输入，输入数据当然是测序序列的fastq文件，输出的表达矩阵。值得注意的是，每个10x样本都有3个fastq文件作为输入，然后输出的表达矩阵，也是3个文件哦！！！

大家在下面的文章里面可以搜索到10x单细胞转录组数据的文章公布在geo数据库的链接：

在教程：使用seurat3的merge功能整合8个10X单细胞转录组样本和 seurat3的merge功能和cellranger的aggr整合多个10X单细胞转录组对比我也给出了后续R代码读取10x单细胞转录组数据的3个文件的表达矩阵。

如果是10x的单细胞公共数据

比如 GSE128033 和 GSE135893，就是10x数据集，随便下载其中一个，就能看到每个样本都是走流程拿到10x单细胞转录组数据的3个文件的表达矩阵。

2.2M Mar  8  2019 GSM3660655_SC94IPFUP_barcodes.tsv.gz
259K Mar  8  2019 GSM3660655_SC94IPFUP_genes.tsv.gz
 26M Mar  8  2019 GSM3660655_SC94IPFUP_matrix.mtx.gz
2.2M Mar  8  2019 GSM3660656_SC95IPFLOW_barcodes.tsv.gz
259K Mar  8  2019 GSM3660656_SC95IPFLOW_genes.tsv.gz
 31M Mar  8  2019 GSM3660656_SC95IPFLOW_matrix.mtx.gz
2.2M Mar  8  2019 GSM3660657_SC153IPFLOW_barcodes.tsv.gz
259K Mar  8  2019 GSM3660657_SC153IPFLOW_genes.tsv.gz
 33M Mar  8  2019 GSM3660657_SC153IPFLOW_matrix.mtx.gz
2.2M Mar  8  2019 GSM3660658_SC154IPFUP_barcodes.tsv.gz
259K Mar  8  2019 GSM3660658_SC154IPFUP_genes.tsv.gz
 31M Mar  8  2019 GSM3660658_SC154IPFUP_matrix.mtx.gz

下游处理的时候，一定要保证这3个文件同时存在，而且在同一个文件夹下面。

示例代码是：

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
sce1 <- CreateSeuratObject(Read10X('../10x-results/WT/'),
                          "wt")

重点就是 Read10X 函数读取文件夹路径，比如：../10x-results/WT/ ，保证文件夹下面有3个文件。

然后针对Smart-seq2数据

这个其实就是普通的转录组数据处理流程哦，比如我们看2017-scRNA-seq-primary breast cancer，韩国研究团队是这样描述的：

image-20200206154712495

其实Smart-seq2单细胞最出名的应该是北京大学张泽民教授团队的多种癌症免疫浸润T细胞测序，包括肝癌，结直肠癌，肺癌的，比如 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE98638

这样的表达矩阵，每个样本只有一个文件即可，需要注意的是，提供多种多样数据格式，需要自行阅读量过两万的综述翻译及其细节知识点的补充：

上游分析对大家来说比较困难，除非你有Linux基础，代码如下：

hisat2=/home/jianmingzeng/biosoft/HISAT/hisat2-2.0.4/hisat2
# # 如果使用conda安装的 hisat2，那么 hisat2 命令应该是在环境变量的。
## 索引文件需要自己下载
# https://ccb.jhu.edu/software/hisat2/manual.shtml
# wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz

index=/home/jianmingzeng/reference/index/hisat/mm10/genome
ls raw_fq/*gz |  while read id; do
$hisat2 -p 10 -x $index -U $id  -S ${id%%.*}.hisat.sam
done

ls *.sam|while read id ;do (samtools sort -O bam -@ 5  -o $(basename ${id} ".sam").bam   ${id});done
rm *.sam
ls *.bam |xargs -i samtools index {}

## gtf文件推荐去gencode数据库下载
gtf=/home/jianmingzeng/reference/gtf/gencode/gencode.vM12.annotation.gtf
featureCounts=/home/jianmingzeng/biosoft/featureCounts/subread-1.5.3-Linux-x86_64/bin/featureCounts   
# # # 如果使用conda安装的 subread，那么featureCounts  命令应该是在环境变量的。
$featureCounts -T 5 -p -t exon -g gene_id  -a $gtf -o  all.id.txt  *.bam  1>counts.id.log 2>&1 &
 # 得到的就是 count矩阵

大家可能会觉得奇怪，为什么我给到的代码里面的软件，都不是截图文献里面使用的呢？其实转录组数据处理流派太多了，并没有绝对的权威，反正我们生信技能树的粉丝流程都是从我这里教出去的，走hisat2和featureCounts流程来定量拿到表达矩阵，也有文献这样写，如下：

走hisat2和featureCounts流程

同样的，很大概率你不需要自己处理上游数据，而是公共数据库，比如我通常是下载上面的count矩阵，然后走代码：

rm(list=ls())
options(stringsAsFactors = F)
# install.packages('R.utils')
# install.packages('data.table')
library(data.table)
# 这个表达矩阵其实是10x的，不过可以演示
a=fread('GSE117988_raw.expMatrix_PBMC.csv.gz',header = TRUE)
length(a$V1)
length(unique(a$V1))
hg=a$V1
dat=a[,2:ncol(a)]
rownames(dat)=hg
hg[grepl('^MT-',hg)]
colnames(dat)
rownames(dat)
meta=as.data.frame(colnames(dat))
colnames(meta)=c('cell name')
rownames(meta)=colnames(dat)
head(meta)
## 前面大量的代码，都是数据预处理
library(Seurat)
dat[1:4,1:4]
class(dat)
# 重点是构建 Seurat对象
pbmc  <- CreateSeuratObject(counts = dat,
                            meta.data = meta,
                            min.cells = 3, min.features = 200,project = '10x_PBMC')

pbmc
head(colnames(dat))
head(rownames(dat))  
rownames(GetAssayData(pbmc,'counts'))
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

两种单细胞转录组数据，都是走Seurat流程，只不过是构建 Seurat对象的代码不一样，注意仔细分辨！

很大概率上你并不会需要自己走上游流程

主要是因为对计算资源的消耗，实验室搭建上游流程成本太高，还不如一次性付费让公司做出来表达矩阵给到你后下游慢慢探索。

但是懂上游分析流程有助于你更好的认识你的单细胞数据！

为什么10x单细胞转录组表达矩阵有3个文件

因为10x单细胞转录组表达矩阵里面的0值非常多，所以换成3个文件存储更节省空间。

本质上仍然是一个表达矩阵而已，如果你都有了表达矩阵，就没必要去想那3个文件了。

样本量多少都可以发

但凡看过单细胞转录组相关文献的都知道，样本量跟文章所发的杂志影响因子是正相关的。样本数量多少都可以发的，不同数量的样本，能讲好的生物学故事不一样，自然发表的水平不一样咯。

我在单细胞天地分享过，不同数量样本的单细胞项目文章：

在教程：使用seurat3的merge功能整合8个10X单细胞转录组样本和 seurat3的merge功能和cellranger的aggr整合多个10X单细胞转录组对比其实就指出来了，样本数量再多也可以分析。

补充材料里面对样本量和数据量描述得很清楚

单个样本的单细胞转录组很少见了，现在以2个样本项目居多，一个对照一个处理，如果是常规转录组，两个分组的话每个组通常是3个样本，但是我们说了嘛，单细胞还是很贵，单个10x样本还是在3万左右的费用，大多数课题组就是想尝个鲜，而且学术界也承认单细胞的确很贵大家hold不住，所以2个样本项目也可以发出去的啦。

下面就是一个典型的2个样本项目样本量描述：

如果样本量比较多，复杂的实验设计，不同生物学假设的分组，可以用示意图发方式，如下；

单个病人取样示意图

关于测序数据量

测序数据量，其实就是文库大小，每个细胞的reads的总数。

上面截图的10x的单细胞转录组的两个样本的课题，就写清楚了，每个样本是3500左右细胞，每个细胞是8万的reads数量。

实际上10x单细胞转录组每个样本可以是3000到10000的细胞数量都可以，取决于实验设计。每个细胞平均可以是1到10万的reads，都没有问题。

如果是Smart-seq2技术呢，每个样本细胞数量通常是96的倍数，500个左右就很厉害了，然后每个细胞的reads就可以很多，百万级别都没有问题。

记住这些关键数量即可，如果想知道具体原理，多看几个综述即可！

测序数据量和捕获细胞数量对结果的影响

10X官方有PBMC单细胞测试数据，4000K细胞，每个细胞平均是50K的reads。然后取，500, 1k, 2,5k, 5k, 7.5k, 10k, 15k, 25k的随机抽样子集，同样的，取100，200，400，600，500，1，2，3，4K的随机抽样子集。

发现平均每个细胞的0.5K和86K的reads测序量，检测到细胞数量都有4000，而且极低深度测序，仍然是可以比较清晰可见的区分细胞亚型，哪怕在每个细胞的0.5K这样的reads数量情况下每个细胞仅仅是能检测到160个左右的基因。

当然了，这样的探索，仅仅是计算层面的指导，其实你肯定是在公司测序，那么公司的人一定会推荐你每个样本是3~8K细胞，平均每个细胞15-50K的reads，这样的测序策略。

现在一起看看拿到表达矩阵后有一个很重要的质控过程，就是根据一些阈值来过滤不合格细胞和基因。、细胞的不合格很容易理解，可能是那个细胞检测到的基因数量太少，或者太多，也许是那个细胞的文库大小异常，都是需要谨慎考虑是否管理它。

基因的不合格，就需要大家对人类参考基因组的注释信息有所背景了。走RNA-seq定量流程，拿到的表达矩阵通常是取决于gtf文件的注释程度，人类的gtf里面五万多个基因，不可能都在你的单细胞转录组项目数据里面出现。

首先细胞的取舍主要是看文库大小和检测到的基因数量

单细胞转录组通常使用10x数据，所以细胞数量惊人，每个样本可以是3000到10000的细胞数量都可以，取决于实验设计。每个细胞平均可以是1到10万的reads文库，都没有问题。

那么，每个细胞可以检测到多少基因数量呢，就取决于每个细胞分到的reads总数，1到10万的reads文库对应着200到1000的基因数量。

但是Smart-seq2技术的单细胞数据就不一样了，每个样本细胞数量通常是96的倍数，500个左右就很厉害了，然后每个细胞的reads就可以很多，百万级别都没有问题。所以检测到基因数量也很多，如下：

所以在你设置阈值来过滤不合格细胞的时候，一定要先理解好你的单细胞转录组数据，给全部的细胞检测到的基因数量绘制一个boxplot，看看哪些细胞所检测到的基因数量偏多或者偏少，一般来说，这样的的离群细胞就是你需要去除的！比如我在单细胞转录组教学就介绍过一个代码：

box <- lapply(colnames(sample_ann[,1:19]),function(i) {
    dat <-  sample_ann[,i,drop=F] 
    dat$sample=rownames(dat)
    ## 画boxplot 
   ggplot(dat, aes('all cells', get(i))) +
          geom_boxplot() +
          xlab(NULL)+ylab(i)
})
plot_grid(plotlist=box, ncol=5 )

这样的话细胞的上游指标可以批量检验，如下：

可以看到，每个细胞的测序文库大小和检测到基因数量虽然很重要，上游分析的其它质量指标，比对率，质量分数，duplicate情况，插入片段，内含子误差其实也会影响这个细胞的取舍。

然后是删除一些基因

前面我们提到过，gtf文件的注释程度不一样，拿到的表达矩阵通常是全部的基因，比如人类的gtf里面五万多个基因，你的表达矩阵就是5万多行。

实际上，大量的基因在你所有的细胞里面都是不表达的，这就是表达量全部为0的，肯定是需要去除啦。

还有一些基因仅仅是在5%左右的细胞表达，这个时候就很主观了，可能确实这个基因是那5%的稀有细胞的marker基因，也有可能就是单细胞转录组技术导致这个基因大量的drop-out了。

这个时候，你就需要卡一个阈值，到底一个基因表达多少算是表达呢，到底一个基因在多少个细胞里面有表达，你才保留下去，做下游分析呢？

后记

当然了，过滤细胞或者基因仅仅是单细胞转录组项目质量控制的开始，后面还有线粒体基因含量，核糖体基因含量，高表达量基因比例，细胞周期，批次效应，我们慢慢聊哈。

我们埋下了一个伏笔，就是拿到了差不多干净的表达矩阵后还有一个预处理步骤，但是这个步骤是选修，就是说你不会过滤线粒体核糖体基因也可以不做，没有人会责怪你，因为单细胞转录组数据分析本来就有难度。分析要点无穷无尽，即使你真的准备过滤线粒体核糖体基因你也会发现标准不好把握，需要看很多文献，还得结合你自己项目数据的实际情况啦。而且过滤线粒体核糖体基因并不是质控的终点，你还有细胞周期检查，单细胞活性检查，是否有两个细胞混在一起也是需要检查，非常的累。不仅仅是标准不好把握，而且因项目而异，具体情况具体判断。

处理线粒体核糖体基因的几个策略

粗暴的去除全部线粒体核糖体基因

直接在表达矩阵里面，去除掉属于的那一行表达量即可，有点类似于甲基化芯片数据分析，直接去除性染色体上面的全部探针，如下所示：

粗暴的去除线粒体核糖体基因

选择一个阈值来去除高表达量细胞

如下所示：

image-20200220113653069

这个阈值很大程度上取决于你对自己项目的了解程度，不同器官组织提取的单细胞，本来就线粒体基因平均水平不一样，不能一刀切！

这个只能是靠你看文献来获取认知了，多看多学习！

我的经验是，多次反复查看线粒体核糖体基因的影响，分群前后看，不同batch看，多次质控图表里面显示它，判断它是否是一个主要因素。

人类和小鼠的线粒体核糖体基因有哪些呢

我这里以gencode数据库的gtf文件为标准，在人类的gencode.v32.annotation.gtf 文件里面，可以查找到37个:

gene_id "ENSG00000210049.1"; gene_name "MT-TF"; hgnc_id "HGNC:7481"
gene_id "ENSG00000211459.2"; gene_name "MT-RNR1"; hgnc_id "HGNC:7470"
gene_id "ENSG00000210077.1"; gene_name "MT-TV"; hgnc_id "HGNC:7500"
gene_id "ENSG00000210082.2"; gene_name "MT-RNR2"; hgnc_id "HGNC:7471"
gene_id "ENSG00000209082.1"; gene_name "MT-TL1"; hgnc_id "HGNC:7490"
gene_id "ENSG00000198888.2"; gene_name "MT-ND1"; hgnc_id "HGNC:7455"
gene_id "ENSG00000210100.1"; gene_name "MT-TI"; hgnc_id "HGNC:7488"
gene_id "ENSG00000210107.1"; gene_name "MT-TQ"; hgnc_id "HGNC:7495"
gene_id "ENSG00000210112.1"; gene_name "MT-TM"; hgnc_id "HGNC:7492"
gene_id "ENSG00000198763.3"; gene_name "MT-ND2"; hgnc_id "HGNC:7456"
gene_id "ENSG00000210117.1"; gene_name "MT-TW"; hgnc_id "HGNC:7501"
gene_id "ENSG00000210127.1"; gene_name "MT-TA"; hgnc_id "HGNC:7475"
gene_id "ENSG00000210135.1"; gene_name "MT-TN"; hgnc_id "HGNC:7493"
gene_id "ENSG00000210140.1"; gene_name "MT-TC"; hgnc_id "HGNC:7477"
gene_id "ENSG00000210144.1"; gene_name "MT-TY"; hgnc_id "HGNC:7502"
gene_id "ENSG00000198804.2"; gene_name "MT-CO1"; hgnc_id "HGNC:7419"
gene_id "ENSG00000210151.2"; gene_name "MT-TS1"; hgnc_id "HGNC:7497"
gene_id "ENSG00000210154.1"; gene_name "MT-TD"; hgnc_id "HGNC:7478"
gene_id "ENSG00000198712.1"; gene_name "MT-CO2"; hgnc_id "HGNC:7421"
gene_id "ENSG00000210156.1"; gene_name "MT-TK"; hgnc_id "HGNC:7489"
gene_id "ENSG00000228253.1"; gene_name "MT-ATP8"; hgnc_id "HGNC:7415"
gene_id "ENSG00000198899.2"; gene_name "MT-ATP6"; hgnc_id "HGNC:7414"
gene_id "ENSG00000198938.2"; gene_name "MT-CO3"; hgnc_id "HGNC:7422"
gene_id "ENSG00000210164.1"; gene_name "MT-TG"; hgnc_id "HGNC:7486"
gene_id "ENSG00000198840.2"; gene_name "MT-ND3"; hgnc_id "HGNC:7458"
gene_id "ENSG00000210174.1"; gene_name "MT-TR"; hgnc_id "HGNC:7496"
gene_id "ENSG00000212907.2"; gene_name "MT-ND4L"; hgnc_id "HGNC:7460"
gene_id "ENSG00000198886.2"; gene_name "MT-ND4"; hgnc_id "HGNC:7459"
gene_id "ENSG00000210176.1"; gene_name "MT-TH"; hgnc_id "HGNC:7487"
gene_id "ENSG00000210184.1"; gene_name "MT-TS2"; hgnc_id "HGNC:7498"
gene_id "ENSG00000210191.1"; gene_name "MT-TL2"; hgnc_id "HGNC:7491"
gene_id "ENSG00000198786.2"; gene_name "MT-ND5"; hgnc_id "HGNC:7461"
gene_id "ENSG00000198695.2"; gene_name "MT-ND6"; hgnc_id "HGNC:7462"
gene_id "ENSG00000210194.1"; gene_name "MT-TE"; hgnc_id "HGNC:7479"
gene_id "ENSG00000198727.2"; gene_name "MT-CYB"; hgnc_id "HGNC:7427"
gene_id "ENSG00000210195.2"; gene_name "MT-TT"; hgnc_id "HGNC:7499"
gene_id "ENSG00000210196.2"; gene_name "MT-TP"; hgnc_id "HGNC:7494"

所以你拿到自己的表达矩阵后，其实简单的看看基因名字是否以 MT- 开头即可哦。

所以你会看到seurat包的函数：PercentageFeatureSet 可以用来计算线粒体基因含量。

  sce <- CreateSeuratObject(Read10X('../scRNA/filtered_feature_bc_matrix/'), "sce")
  head(sce@meta.data)  
  GetAssayData(sce,'counts')
  sce[["percent.mt"]] <- PercentageFeatureSet(sce, pattern = "^MT-")
  VlnPlot(sce, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

上面的方法是修改 sce[["percent.mt"]] ，下面我们演示 AddMetaData 函数，同样是可以增加线粒体基因含量信息到我们的seurat对象。

mt.genes <- rownames(sce)[grep("^MT-",rownames(sce))]
C<-GetAssayData(object = sce, slot = "counts")
percent.mito <- Matrix::colSums(C[mt.genes,])/Matrix::colSums(C)*100
sce <- AddMetaData(sce, percent.mito, col.name = "percent.mito")
sce[["percent.mito"]]

也可以是添加核糖体基因含量，同样的你需要知道核糖体基因的名字规则：

rb.genes <- rownames(sce)[grep("^RP[SL]",rownames(sce))]
percent.ribo <- Matrix::colSums(C[rb.genes,])/Matrix::colSums(C)*100
sce <- AddMetaData(sce, percent.ribo, col.name = "percent.ribo")

在人类的gencode.v32.annotation.gtf 文件里面，可以看到核糖体基因数量也不少哦。

核糖体基因在人类列表

如果是小鼠，通常是基因名字大小写替换一下：

小鼠的核糖体基因列表

不仅仅是线粒体核糖体基因

还可以去免疫相关基因，缺氧相关基因，就更加的需要深入到你自己的课题，其实细节是无穷无尽的，但是我们的教学没办法做到如此的个性化，只能是精炼了常规单细胞转录组数据分析主线，就是5大R包， scater,monocle,Seurat,scran,M3Drop，然后10个步骤：

step1: 创建对象
step2: 质量控制
step3: 表达量的标准化和归一化
step4: 去除干扰因素(多个样本整合)
step5: 判断重要的基因
step6: 多种降维算法
step7: 可视化降维结果
step8: 多种聚类算法
step9: 聚类后找每个细胞亚群的标志基因进行亚群命名
step10: 继续分类

其实质量控制三部曲，还有一个很关键的点没有讲解，就是多个样本整合，并且区分批次效应和生物学差异。但是这个点很大程度是依赖于经验，就是说，要想搞清楚，需要写很多自定义的代码去一点一滴的探索，而不仅仅是流程。所以我们就只能介绍到这里，假设大家都拿到了干净的表达矩阵，而且可以很肯定的说这个表达矩阵做下游分析是ok的。

那么我们就来看看，有了干净的表达矩阵后下游分析的第一个分析要点就是：归一化和标准化。

归一化和标准化的区别

实际上口语里面通常是没办法很便捷的区分这两个概念，我查阅了大家的资料，发现基本上都是混淆在z-score转换和0-1转换这两个上面。所以我这里把归一化和标准化替换成为去除样本/细胞效应或者去除基因效应：

首先去除样本/细胞效应：因为不同样本或者细胞的测序数据量不一样，所以同样的一个基因在不同细胞，哪怕你看到的表达量是一样的，但是背后细胞整体测序数据量的差异其实反而说明了这个基因在不同细胞表达量其实是有差异的。在seurat3里面的代码是：
```
sce <- NormalizeData(sce, normalization.method =  "LogNormalize", scale.factor = 10000)
```
默认文库大小是1万，默认会取log值。
然后去除基因效应，这个主要是在绘制热图的时候会需要使用，因为个别基因表达量超级高，在热图里面一枝独秀，实际上我们并不会关心不同基因的表达量高低，我们仅仅是想看指定基因在不同细胞的高低而已，这样的话，就把该基因的表达量在不同细胞的数值，进行z-score这样的转换。在seurat3里面的代码是：
```
sce <- ScaleData(sce) 
```

这样处理后的表达矩阵，就可以进行后续的降维聚类分群啦，我们下期再讲。

取log让表达量离散程度降低

这个时候眼尖的朋友其实看到了，在使用NormalizeData函数的时候，有一个 normalization.method = "LogNormalize" 参数被设置了，这个是为什么呢？

其实很简单，原始的raw counts矩阵是一个离散型的变量，离散程度很高。有的基因表达丰度比较高，counts数为10000，有些低表达的基因counts可能10，甚至在有些样本中为0。因为是单细胞转录组，drop-out现象很严重，其实大量的基因在很多细胞的表达量都是0，如果是10x数据，甚至会出现一个表达矩阵里面97%的数值都是0 的现象。

如果对表达量取一下log10，发现10000变成了4，10变成了1，这样之前离散程度很大的数据就被集中了。

有时当表达量为0时，取log会出现错误，可以log(counts+1)来取log值。当x=1时，所有的log系列函数值都为0。这样原本表达量为0的值，取log后仍为0。

这也就是UCSC的XENA下载到的表达矩阵的形式。

使用GetAssayData函数查看不同形式的表达矩阵

在seurat3里面，很方便进行各种形式的归一化或者标准化，同样的，也很容易查看处理前后的表达矩阵，使用GetAssayData函数即可。

This function can be used to pull information from any of the slots in the Assay class. For example, pull one of the data matrices("counts", "data", or "scale.data").

如下代码：

# 最原始数据
GetAssayData(sce,'counts')[1:6,1:6]
# 去除了细胞测序数据量后
1/(colSums(as.data.frame(GetAssayData(sce,'counts')))[2]/10000)
GetAssayData(sce,'data')[1:6,1:6]
# z-score后
GetAssayData(sce,'scale.data')[1:6,1:6]

如下结果：

> # 最原始数据
> GetAssayData(sce,'counts')[1:6,1:6]
6 x 6 sparse Matrix of class "dgCMatrix"
             6A-11 6A-13 6A-14 6A-15 6A-16 6A-17
DDX11L1          .     .     .     .     .     .
WASH7P           .     1     .     .     .     .
MIR6859-1        .     .     .     .     .     .
RP11-34P13.3     .     .     .     .     .     .
OR4G11P          .     .     .     .     .     .
OR4F5            .     .     .     .     .     .
> # 去除了细胞测序数据量后
> 1/(colSums(as.data.frame(GetAssayData(sce,'counts')))[2]/10000)
      6A-13 
0.009725971 
> GetAssayData(sce,'data')[1:6,1:6]
6 x 6 sparse Matrix of class "dgCMatrix"
             6A-11       6A-13 6A-14 6A-15 6A-16 6A-17
DDX11L1          . .               .     .     .     .
WASH7P           . 0.009678978     .     .     .     .
MIR6859-1        . .               .     .     .     .
RP11-34P13.3     . .               .     .     .     .
OR4G11P          . .               .     .     .     .
OR4F5            . .               .     .     .     .
> # z-score后
> GetAssayData(sce,'scale.data')[1:6,1:6]
                  6A-11      6A-13      6A-14      6A-15      6A-16      6A-17
RP5-857K21.9 -0.8302953 -0.8302953 -0.7327066 -0.7224892 -0.7683106 -0.7692732
RP5-857K21.8 -0.5066862 -0.5688591 -0.5349974 -0.4881690 -0.3650311 -0.4524682
PRKCZ        -0.8350477 -0.5996852 -0.8350477 -0.8350477 -0.5988844 -0.8350477
TNFRSF14     -0.2550752 -0.2550752 -0.2550752 -0.2550752 -0.2550752 -0.2550752
TP73         -0.2761133 -0.2761133 -0.2761133 -0.2761133 -0.2761133 -0.2337521
HES2         -0.4569104 -0.4569104 -0.4569104 -0.4569104 -0.4569104

其实样本/细胞效应不仅仅是文库大小

每个细胞测序数据量的不一致是很容易理解的，但其实细胞之间还有很多其它效应，比如线粒体基因含量，ERCC含量等等，那些处理起来，其实就是深入了解我们讲解seurat里面的NormalizeData和ScaleData函数。

降维聚类分群是一条龙分析

我们并不是开发单细胞数据处理算法，所以大概率上，大家其实会把降维聚类分群一起做了，在seurat3里面的代码是：

sce <- NormalizeData(sce, normalization.method = "LogNormalize", scale.factor = 10000)
sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 2000)
# 步骤 ScaleData 的耗时取决于电脑系统配置（保守估计大于一分钟）
sce <- ScaleData(sce)
sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce))
sce <- FindNeighbors(sce, dims = 1:15)
sce <- FindClusters(sce, resolution = 0.2)
# 这个 resolution 可以调整，值越大，分出来的细胞亚群越多，默认是 0.8
table(sce@meta.data$RNA_snn_res.0.2)

首先看两种降维

简单解释一下，这代码里面的FindVariableFeatures和RunPCA函数，是两种不同策略的降维。

首先FindVariableFeatures是硬过滤，根据一些统计指标，比如sd,mad,vst等等来判断你输入的单细胞表达矩阵里面的2万多个基因里面，最重要的2000个基因，其余的1.8万个基因下游分析就不考虑了。
然后RunPCA函数其实跑完之后2000个基因会转变为2000个维度，但是我们通常看前面的十几个维度就ok了，所以也是一个效率非常高的降维方式。

这个时候，你一定有疑问，为什么FindVariableFeatures是挑选2000个基因而不是其它数量呢？RunPCA函数跑完后我们应该是挑选前多少个维度呢？10个还是15个，还是20个还是50呢。

然后看聚类分群

聚类分群是紧密连接的，细胞可以看做是空间的不同点，如果是二维平面空间，点与点之间的距离很方便计算，距离的远近就决定着细胞是否属于一个类群。

计算距离的公式很多，我们就不一一展开，但是需要注意的是二维平面空间，三维球体空间的细胞距离很方便计算，但是如果是50个维度的空间，计算几万个细胞之间的距离就很可怕了，如果是2000个维度，甚至是2万个维度，基本上个人计算机就可以放弃了。这就是为什么我们前面通常是需要降维的。

下面是一个典型的单细胞转录组项目数据处理的描述：

降维聚类分群

可以看到他们的第一步降维是，选取top 5000的表达量离散度大的基因，第二步降维是选取top20的主成分。使用KNN-graph的聚类，最终定下来了10个细胞亚群。

一般来说，如果单细胞转录组数据仅仅是文章生物学故事的一个环节，就会采取标准的seurat流程，如下所示：

标准的seurat流程

如果你看的文献足够多，还会发现，在降维聚类分群之后，通常是有一个细胞在二维平面的散点图展示，如下所示：

细胞在二维平面的散点图

如果你足够心细，也会发现其实细胞的空间距离排布坐标通常是tSNE和umap来展现。

那我们说的tSNE和umap是怎么回事呢

tSNE(t-distributed stochastic neighbor embedding) 是一种非线性降维的算法，是流形学习的一种，当然，你大概率上是无需理解流形学习的，认识一下其它流形学习的经典方法即可，包括：Isomap，LLE，LE和diffusion maps等。我发现这篇推文介绍的非常好：单细胞中的流形（一）：理解 tSNE中的perplexity，看完你需要记住的是；

困惑度（perplexity）可以表示细胞的邻近个数，在tSNE图上的直观反映是细胞点的分布是否紧凑。perplexity设置越大，细胞分布越紧凑。
tSNE的参数设置：perplexity < (细胞数-1)/3，建议perplexity = 细胞数 / 50；
tSNE倾向于保留数据的局部结构。

本来我还期待着他们继续讲解umap，但是看起来似乎公众号已经停更了，唉，现在搞分享，写原创教程也挺不容易的，很容易就quit了。

我给大家的策略是，反正你得多尝试，umap啥都尝试一遍，最后选择效果好的图表去展现即可。

可以修改的参数

第一次降维，表达矩阵里面通常是2~3万的基因数量，需要筛选到千这个数量级，可以是1000-5000，都是可以

因为参数需要自己摸索和调整，所以其实拿到细胞亚群数量是因而而异的，取决于你前面降维的程度，分群的算法和参数。不过最重要的是拿到了不同细胞亚群后需要对它进行命名，给出生物学的解释。不同的人分析同一个数据集，有略微不同的结果是可以接受的，保证自己的生物学故事圆满即可。

细胞亚群注释依赖标记基因

我给大家的单细胞进阶课程里面，示例文章《Acquired cancer resistance to combination immunotherapy from transcriptional loss of class I HLA》，就是一个使用seurat标准流程对PBMC分群后如下:

seurat标准流程对PBMC分群

得到了这13个细胞类群，而且也细致的决定它们的亚群名字。实际上PBMC的不同细胞亚群的标记基因是比较明确，的比如：在文章：Nucleic Acids Res. 2018 Apr ，dropClust: efficient clustering of ultra-large scRNA-seq data.就写的很清楚：

PBMC细胞亚群的标记基因列表

这些标记基因在不同亚群细胞的表达量热图或者小提琴图展示一下，就明白了为什么它们可以作为标记基因，来对细胞亚群进行命名啦。该系列教程如下：

如果你看完这些教程，就应该知道，并不是所有你拿到细胞亚群，都是有明确生物学功能的，你的文章所要讲述的生物学故事也并不需要把全部的细胞亚群一一记流水账。

标记基因的不同来源

这样的标记基因列表，有一些网页工具会收集，比如cellmarker CellMarker: a manually curated resource of cell markers ，作者：X Zhang - ‎2019 。就需要自行学习了，也有自己查询自己领域内的全部文献，然后整理出来标记基因列表。这个步骤至少耗时2个月，比如2018年发表于science杂志的文章就是自己根据文献进行整理的：Single-cell transcriptomes from human kidneys reveal the cellular identity of renal tumors.

多种策略来进行亚群注释

其实有多种策略来进行亚群注释，但是都很耗费精力。可以选择每个细胞亚群的差异表达基因进行简单的富集分析，也可以查询文献，选取文献里面报道过的细胞亚群特异性高表达基因作为标记基因。

查询文献这个工作量是蛮大的，所以一般来说，作者也会把他们最后总结好的细胞亚群注释使用的标记基因整理成为一个表格，如下：

查询文献后整理的标记基因列表

也可以通过差异基因来进行注释

比如下面的描述，取top100的差异基因，即便是如此，研究者也仍然是找到了不少已知的标记基因来可视化验证。

取top100的差异基因

单个标记基因可视化以小提琴图和热图展现

两个可视化方法被集成在seurat包里面，如下所示

小提琴图或者热图

代码通常是：

markers_df <- FindMarkers(object = sce, ident.1 = 8, min.pct = 0.25)
print(x = head(markers_df))
markers_genes =  rownames(head(x = markers_df, n = 5))
VlnPlot(object = sce, features =markers_genes,log =T )
FeaturePlot(object = sce, features=markers_genes )

就是 VlnPlot 和 FeaturePlot 两个函数，当然，也是可以自己获取细胞的二维平面坐标，以及指定基因的表达值，自定义绘图函数。

多个标记基因的可视化也是热图为主

下面的例子是，把细胞亚群的基因表达值取一个统计量（平均值或者中位值）来作为该基因在该细胞亚群的表达量。所以热图展现多个标记基因在全部细胞亚群的表达量热图，也不会拥挤：（后面我们会讲解这个图如何绘制）

多个标记基因的可视化热图

但是，如果要展现多个标记基因在全部细胞的表达量，而不仅仅是细胞亚群，热图就会很拥挤，如下所示：

多个标记基因在全部细胞的表达量热图

可以看到，每个细胞亚群虽然从注释的角度来说，指定的几个高表达量基因即可定义它，但是它本身特异性高表达的其它基因也会很多。所以作者单独把值得一提的标记基因，在热图坐标高亮出来。

细胞亚群可以分成主要亚群和次要亚群

Impact of the Microenvironment on Breast Immune Cells
(A) Breast immune cell atlas constructed from combining all patient samples (BC1–8) and tissues using Biscuit, projected with t-SNE. Each dot represents a cell, colored by cluster; major cell types are marked according to Figure 2F, H and Table S2, 3.

其实理论上细胞亚群是可以无限划分的，因为世界上没有两个一模一样的细胞，关键是要把握一个度，什么样的差异可以判定为不同细胞亚群，什么样的差异是可以容忍的细胞类群内部异质性。

有一个策略就是找出主要因素和次要因素。主要因素划分为主要亚群，比如外周血里面的T,B细胞当然是不同亚群，但是T细胞里面还可以继续划分：CD4或者CD8的T细胞，甚至继续划分，如下图所示：

T细胞的多层次划分亚群

主要和次要细胞亚群同一个tSNE图展现

现在10x单细胞转录组技术一个样本出来3~8K细胞的数据，很多课题都是十几个以上样本，所以妥妥的十万多个单细胞数量级，它们走聚类分群对技术资源的消耗当然是很可观的，而且，拿到的细胞亚群数量也是惊人。比如下面这个，就分成了近100个细胞亚群。

实际上，细胞二维散点图，是没办法写全部细胞亚群的生物学功能定义的，我们通常也是把主要细胞亚群标记上去。然后每个细胞亚群再进行继续分群。只不过是在同一个散点图上面展示。

主要和次要细胞亚群以标记基因来展现

上面的散点图毕竟展示的细胞数量太多，大多数情况下以炫酷为主，一般人很难看出来不同细胞亚群内部具体如何划分子亚群，以及不同亚群到底以什么标记基因来进行区分。

有一个策略是，把标记基因的表达量在所有亚群及子亚群里面展现，以热图形式：

e, Clustergrammer heatmap showing protein marker expression (top) in each MC (left) and the canonical annotation of these communities (right).

The dendrogram bars (light gray) indicate the clustering of MCs based on the cosine distance method in Clustergrammer.

标记基因在所有亚群及子亚群里面展现

这个其实有点像我们前些天在生信技能树提到的学徒练习题：

后来我升级为了bodymap和Gtex数据库的，指定基因在指定组织里面的表达量热图，就需要使用代码，把同一个亚群的全部细胞表达量综合一下：

mat=do.call(rbind,
        lapply(unique(id$SMTS), function(t){
          rowMeans(dat[,id$SMTS==t])
        })
        )

其中 id$SMTS 里面是每个细胞的亚群属性，而dat是我们的表达矩阵，所以对表达矩阵来说，依次取出每个亚群的表达矩阵子集，然后取 rowMeans，就拿到了每行的基因在该亚群的表达量平均值。

这个综合后的表达量值就可以去绘制上面的标记基因的表达量在所有亚群及子亚群里面热图。

多个样本可以分开走聚类分群流程

比如中山大学的最新研究《一个人的15个器官单细胞测序数据》，链接是：https://www.biorxiv.org/content/10.1101/2020.03.18.996975v1.full.pdf 。这里面首先有不同器官，可以分开独立走单细胞流程。然后全部的八万多细胞，分成了主要细胞群之后，仍然是可以进行每个亚群细致研究：

we performed single-cell transcriptomes of 88,622 cells derived from 15 tissue organs of one adult donor and generated an adult human cell atlas (AHCA).
首先是 a total of 20,494 T cells，可以区分成为 CD4+ (7,122) and CD8+ (13,372) T cells 两个大类，然后继续每个大类细分为11 and 22细胞亚群。
然后是 10,655 B and plasma cells
还有 5,605 myeloid cells ，可以细分为7个monocyte亚群，9个macrophage亚群，一个small dendritic cell (DC) 亚群。
继而是 18,090 epithelial cells ，进而细分为33个亚群，值得注意的是这些上皮细胞来源于9种器官，189 genes with tissue-specific expression (FC 5, pct.1 0.2) 。
然后是7,137个 Endothelial cells (ECs) 细胞, including 6,863 blood endothelial cells (BECs, marked with VWF) and 274 lymphatic endothelial cells (LECs, marked with LYVE1).
最后是 a total of 17,835 fibroblasts and smooth muscle cells ，分成
14 fibroblast clusters (11,767 cells, MMP2),
four smooth muscle cell clusters (3,201 cells, ACTA2),
another five clusters assigned as FibSmo (2,867 cells; marked with MMP2 and ACTA2; )

单细胞测序的肺癌肿瘤免疫微环境中基质细胞

文章是：Phenotype molding of stromal cells in the lung tumor microenvironment。

基本上也是细胞分群后继续分群，文献里面公布的主要亚群标记基因是：

endothelial cells：CLDN18, FOLR1, AQP4 and PEBP4
endothilial：CLDN5, FLT1, CDH5 ，RAMP2
epithelial cells: CAPS, TMEM190, PIFO，SNTN
fibroblasts,COL1A1, Decorin (DCN) Collagen type I alpha 2 (COL1A2) and C1R 9 ; B-cells,：genes encoding the B-cell antigen receptor complex-associated protein alpha (CD79A),Immunoglobulin Kappa Constant (IGKC), Immunoglobulin Lambda Constant 3 (IGLC3) and Immunoglobulin Heavy Constant Gamma 3 (IGHG3)
myeloid cells, genes encoding Lysozyme (LYZ), the Macrophage Receptor With Collagenous Structure (MARCO), CD68 and CD16a, the Fc Fragment Of IgG Receptor IIIa (FCGR3A) 10
T-cells：(TRAC, TRBC1 and TRBC2).