单细胞测序实践

前言

1. mRNA中的polyA尾巴作用

polyA尾的作用是维持mRNA作为翻译模板的活性并增加其mRNA本身的稳定性(防止核酸外切酶)。
自然界中极大多数真核细胞核内转录的mRNA，会在其3`末端修饰“加尾”成约有180-200个单一的腺苷酸残基构成的多聚腺苷酸链。
多聚腺苷酸的产生不依赖与DNA序列，而是与转录的终止同时进行。

2. 外周血

骨髓是人体的造血组织，外周血是除骨髓之外的血液，临床上常用一些方法将骨髓中的造血干细胞释放到血液中，再从血液中提取分离得到造血干细胞，我们把这样得到的干细胞称为外周血干细胞。
外周血单个核细胞（pbmc）是外周血中具有单个核的细胞，包括淋巴细胞和单核细胞。

3. tsv与csv的区别

tsv数据是以Tab制表符‘\t’作为分隔符，而csv数据以’,’ 作为分隔符

4. R中make.unique 函数

通过对重复项追加序号，使字符向量的元素唯一
参数为sep，即选择分隔符，默认为’.’

5. R中ReadLines函数

整行读取，每行作为一个字符串，返回向量

6. R中 Matrix::readMM函数

基因表达矩阵Matrix.mtx文件使用Matrix Market

6.1 Matrix Market格式

Matrix Market(MM) 提供了一个简单的机制来促进矩阵数据的转换
在其最初的规范中，定义了以下两种矩阵格式：

• 坐标格式(coordinate):
  适用于存储稀疏矩阵，以（行号，列号，非零元素值）的形式来存储矩阵中的非零元素值
• 数组格式
  用于存储稠密矩阵和稠密向量，存储矩阵中的所有元素的值(默认按列存储)

6.2 支持的数据类型和矩阵结构

• 支持数据类型：实数（real）、复数（complex）、整数（integer）、二进制（pattern）
• 支持的矩阵结构：通用general（正常存储所有非零元）、Hermitian（下三角）、对称symmetry（只存储下三角）、斜对称矩阵skew-symmetric（存储下三角但没有对角线，因为斜对称矩阵为矩阵加上转置矩阵为0，所以斜对称矩阵对角线元素为0，不必存储）

6.3 MM格式存储方式

MM格式使用“.mtx”的ASCII文本文件存储矩阵数据，其中第一行对存储的矩阵类型和存储格式进行了说明，例如：
%%MatrixMarket matrix coordinate real general
其中，coordinate表示为坐标格式存储，存储数据为实数，存储所有非零元
“%”号开头为文件头，包含了矩阵的信息，之后第一行会给出矩阵的行数、列数和非零元素个数，之后才是矩阵数据

6.4 R中读取mtx文件

使用Matrix::readMM函数

7. R中as函数

使用Matrix::readMM中读取的矩阵格式为dgTMatrix，即coordinate三元组（行，列，数据值）,使用as函数将dgTMatrix对象转换为dgCMatrix（按列压缩的稀疏矩阵格式）

7.1 dgRMatrix（CSR Matrix）存储格式

CSR即按行压缩稀疏矩阵的存储格式，采用三类数据来表示矩阵：数值(values)，列号(column indices)，行偏移(row offsets)。数值和列号的含义与COOR(坐标三元组表示)一致，表示元素的列号和元素值。行偏移值为重点，如上图中，为0 2 4 7 9 ，其中 0的含义表示第一行的首个元素在数值列表（values）和列号列表（column indices）上的index，2的含义为第二行的首个元素在数值列表和列号列表上的index，依次类推，第三行首个元素的偏移值为4，第四行首个元素的偏移值为7，最后在行偏移列表（row offsets）最后补上矩阵总的元素个数。CSC（dgCMatrix）是CSR相对应的一种方式，即按列压缩稀疏矩阵

一、 原理

1. 10X Genomics单细胞测序原理

上图最左边为带有标记的凝胶微珠，凝胶微珠逐渐向右流去，在第一个十字路口，放进去细胞悬液和酶，然后凝胶微珠和细胞悬液和酶一起向右流动，到第二个十字路口遇到油，变为油包水结构，正常情况下为一个油包水结构里面含有一个凝胶微珠和一个细胞以及酶，然后若干个油包水结构被收集起来。

油包水结构中的细胞在酶的作用下破裂，释放出细胞内的RNA，这些RNA会连接到凝胶微珠上。
之所以RNA会连接到凝胶微珠上，是因为在凝胶微珠表面覆盖了一层短的核酸序列，这些段的核酸序列如上图所示，都是由4部分结构组成的：R1，10X Barcode，UMI，Poly(dT)VN。其中poly(dT)VN部分是富集T(胸腺嘧啶)的序列，我们知道在RNA的3`末端有一个polyA的尾巴，所以当RNA被释放出来的时候，它的polyA的尾巴刚好和凝胶微珠上的poly(dT)VN序列结合，一个RNA连接到凝胶微珠上的一个序列(一个凝胶微珠上有很多短的核酸序列，一个RNA只结合一条短核酸序列)。
当拿到带有RNA的凝胶微珠之后，就可以在反转录酶的作用下生成cDNA，然后就是cDNA的扩增，就和Bulk RNA-seq里的扩增一样了。

扩增之后的序列如上图所示，两端的P5和P7是为了测序加上的，它们和测序芯片上的寡聚核苷酸链进行互补结合；Read1和Read2为测序引物；中间灰色的横线为要检测的核酸序列；Sample Index是加进去的样本ID，因为虽然建库的时候是分开的，但是测序的时候是多个样本一起的，所以加进去样本ID以作区分；
剩下的10xBarcode和UMI为关键部分：每个凝胶微珠上有很多短序列，但每条短序列上的10xBarcode是一样的，即每个凝胶微珠上只携带一种10xBarcode序列，用来标记细胞（测序完成之后区分哪些序列是属于同一个细胞，就看它们的Barcode是否相同）；UMI全程是Unique Molecular Index，一个凝胶微珠上的UMI是和不相同的，它的作用是对RNA的绝对表达水平进行统计，细胞破裂之后一个RNA分子与一个凝胶微珠上的短序列结合，当对比确定RNA属于哪个基因后，只需要把来自相同基因的RNA聚到一起，统计一下UMI的种类即可，注意是UMI的种类，并不是UMI的个数，避免扩增的偏好性的问题（某个基因表达水平高，转录的mRNA分子数量多，则细胞破裂时RNA分子结合凝胶微珠上的短序列条数数量多，即UMI种类多，因为凝胶微珠上每条短序列上的UMI各不相同）

二、数据集

数据集选择pbmc
数据集下载地址：戳！
数据集下载后为三个文件：gene.tsv、barcodes.tsv、matrix.mtx

• gene.tsv为基因数据，第二列为基因名
  
• barcodes.tsv为细胞数据，每一行为一个细胞
  
• matrix.mtx为基因表达矩阵，共32738个基因，2700个细胞，表达矩阵中共有2286884个非0值，其余每行为非0值的坐标和表达量
  

三、读取数据集

1. 普通读取

# 设置文件路径
dir.10x <- "D://RData/seurat.analysis/data/filtered_gene_bc_matrices/"
genes <- read.table(paste0(dir.10x,"genes.tsv"),stringsAsFactors = F,header = F)$V2
# 标记重复，如a.1,a.2
genes <- make.unique(genes,sep = '.')
# 整行读取，一行为一个字符串,一行代表一个unique细胞
barcodes <- readLines(paste0(dir.10x,"barcodes.tsv"))
# 基因数量 细胞数量 非0元素
# 32738    2700     2286884
# 行数     列数   基因表达量
# 32709     1         4
# 32707     1         1
# 32706     1        10
mtx <- Matrix::readMM(paste0(dir.10x,"matrix.mtx"))
# 强制将mtx(dgTMatrix)转换为dgCMatrix对象
# dgCMatrix：按列压缩的稀疏矩阵
mtx <- as(mtx,"dgCMatrix")
# 细胞为列
colnames(mtx) <- barcodes
# 基因为行
rownames(mtx) <- genes
# 参数筛选基因和细胞
# min.cells : 对于一个基因来说必须最少在n个基因中表达
# min.features : 对于一个细胞来说必须最少200个基因表达量不为0
pbmc <- CreateSeuratObject(counts = mtx, project = "pbmc3k",
                           min.cells = 3, min.features = 200)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

2. 直接读取

dir.10x <- "D://RData/seurat.analysis/data/filtered_gene_bc_matrices/"
pbmc.data <- Read10X(data.dir = dir.10x)
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k",
                           min.cells = 3, min.features = 200)
1
2
3
4

3. Seurat 对象

3.1 Seurat 整体结构

详细视频教程：Seurat对象结构

Seurat对象大概包含assays、meta.data、active.assay、active.ident、reductions、version和commands
其中，

• assays表示分析对象，返回包含一个或多个的Assay对象列表，例如多组学分析时会有多个Assay。Seurat对象调用Assay时使用"[[]]" ，例如pbmc[['RNA']]
• meta.data为dataframe格式，行表示细胞，列是细胞的某些属性，下游分析的结果也通常会加在meta.data上
• active.assay表示当前激活的分析对象。一个Seurat对象中assays中可以有多个assay，但某一时刻只能有一个激活的assay
• active.ident表示细胞的类型注释，同active.assay。一个Seurat对象中可能有多种不同方式得到的细胞注释，但某一时刻只能有一种
• reduction返回一个或多个DimReduct对象的列表,每个DimReduct对象表示一种降维分析后得到的结果，常见的有pca、umap、T-sne
• version表示创建对象时所用的Seurat版本
• commands返回一个列表，表示workflow中的每个步骤所使用的命令和参数以及时间等信息
  

Seurat对象还涉及到两个子对象：Assays和DimReduct

3.2 Assays对象

以Seurat对象实例pbmc为例，可以使用 pbmc[['RNA']]来访问Assay对象
Assay有以下slot：

• counts ： 表示未经处理的原始数据
• data : 表示标准化后的数据
• scale.data : 表示经过scale（z-score）后的数据
• key：每个Assay对象都有一个key，例如“_rna”
• var.features: 表示高表达变异的特征名，可以使用VariableFeatures函数来获取
• meta.features : 类似meta.data。也是dataframe类型，行表示features（对于不同的assay有不同的含义，例如对于单细胞测序rna_seq则是表示基因），列表示features的某些属性

3.3 DimReduct对象

以Seurat对象实例pbmc为例，可以使用 pbmc[['pca']]来访问对象
DimReduct有以下slot：

• cell.embeddings : 表示每个细胞在PC轴上的坐标
• feature.loadings : 每个基因低对每个PC轴上的贡献度
• feature.loadings.projected : 投影特征加载矩阵
• assay.used : 用于降维的assay对象名称
• stdev ： 每个维度的标准差
• key : 设置pc轴的名称

四、QC质量控制

####################################################
# QC
# Seurat对象：
# 1. meta.data：orig.ident存储细胞的样本来源
#               nCount_RNA为每个细胞的UMI数量
#               nFeature_RNA为基因数目
# 实际上下游分析的处理结果也会存储在meta.data中
####################################################
# 计算线粒体比例，过滤掉线粒体含量过高的细胞
# 细胞线粒体含量比较高，则可能为死细胞或损伤细胞
pbmc$percent.mt <- PercentageFeatureSet(pbmc,pattern = "^MT-")
pbmc$log10GenePerUMI <- log10(pbmc$nFeature_RNA) / log10(pbmc$nCount_RNA)

VlnPlot(pbmc,features = c("nFeature_RNA","nCount_RNA","log10GenePerUMI","percent.mt"),ncol = 4)
plot1 <- FeatureScatter(pbmc,feature1 = "nCount_RNA",feature2 = "percent.mt" )
plot2 <- FeatureScatter(pbmc,feature1 = "nCount_RNA",feature2 = "nFeature_RNA" )
plot1 + plot2

# 过滤基因数量小于200和大于2500和线粒体比例大于5%的细胞
pbcm <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

4.1 meta.data

Seurat会自动为每个细胞创建一些元数据，meta.data 格式如下，

其中，每行表示一个细胞，列包括：

• orig.ident : 样品名，默认为创建Seurat对象时赋予的project参数的值
• nCount_RNA : 每个细胞的UMI数目
• nFeature_RNA : 每个细胞所检测到的基因数目

此外，需要根据以上信息来计算另外一些指标来判断数据集，例如

• 每个UMI所检测到的基因数目： 用来了解数据集的复杂性，每个UMI所检测到的基因数目越多，则数据集越复杂。
  pbmc$log10GenePerUMI <- log10(pbmc$nFeature_RNA) / log10(pbmc$nCount_RNA),我们取lg转换结果，以便更好的进行样本之间的对比
• 线粒体比例 ： 过滤掉线粒体基因含量较高的细胞（可能为死细胞）
  Seurat中PercentageFeatureSet() 函数将使用一个模式(pattern) 来搜索基因标识符。对于每一个细胞，将选取特征基因的计数之和，除以所有基因的计数之和，再乘以100。
  pbmc$percent.mt <- PercentageFeatureSet(pbmc,pattern = "^MT-") : 模式 ^MT- 适用于人类基因名，要根据物种进行调整，如果不使用基因名作为基因ID，则不起作用

所有指标如下，

接下来将可视化指标，删除质量低的细胞。

4.2 质量管理

细胞QC通常基于 三个QC协变量进行，UMI数目、基因数目、线粒体基因比例。检查QC协变量的分布，通过阈值过滤掉异常值。这些异常的细胞可能对应垂死的细胞、膜破裂的细胞或双胞(doublet)。

4.2.1 可视化

# 对meta.data 每列画小提琴图
VlnPlot(pbmc,features = c("nFeature_RNA","nCount_RNA","log10GenePerUMI","percent.mt"),ncol = 4)
plot1 <- FeatureScatter(pbmc,feature1 = "nCount_RNA",feature2 = "percent.mt" )
plot2 <- FeatureScatter(pbmc,feature1 = "nCount_RNA",feature2 = "nFeature_RNA" )
plot1 + plot2
1
2
3
4
5

四个变量的小提琴图观察数据分布

• nCount_RNA和percent.mt 的散点图可以观察细胞的线粒体基因的占比水平。
• nCount_RNA和nFeature_RNA的散点图中，每个点表示一个细胞，counts和features一般呈线性关系，斜率越大，每个UMI能检测出更多的gene，说明这批细胞的 基因丰度（基因组中该基因的拷贝数量） 较低，反之斜率较小则说明这批细胞的细胞丰度较高。有时候也不是一条可以拟合的曲线，或者是两条可以拟合的曲线，这也反映出细胞的异质性。

4.2.2 数据筛选

根据项目选择合适的阈值去过滤数据
pbcm <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

五、数据归一化

在获得高质量的单细胞之后，下一步是进行聚类。不过为了进行聚类，需要确定细胞之间表达最不同的基因。然后，使用这些基因来确定哪些相关基因造成细胞之间表达的最大差异。所以先进行数据归一化。
数据归一化对于比较细胞之间的基因表达至关重要，归一化是对原始计数值进行缩放的过程，通过归一化，细胞之间的表达水平更具可比性。
总的来说，归一化目的是为了处理每个细胞的总Count不同的问题，消除文库大小的影响。

####################################################
# Normalization
####################################################
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
1
2
3
4

5.1 log1p

log1p() 是 R 中的一个内置函数，用于计算 1 + x 的精确自然对数，其中 x 是指定值，对于 0 抛出无穷大，对于负值抛出 NaN。即：
$$\log 1p(x)=ln(1+x)$$

5.2 归一化原理

默认方法方式是LogNormalize : 对每个细胞来说，每个基因的特征计数除以该细胞的特征总计数，再乘以scale.factor(默认为10000)，，然后对结果进行log1p对数转换。归一化的数值存储在Assay@data中。

• scale.factor = 10000 的原因是：
  
• 进行log转换的原因
  数值变换范围过大，需要进行log转换，将数据压缩在一个区间内。其次也可以改变数据分布，测序数据本身不符合正态分布，log转换能够让数据趋于正态分布，以便于进行下一步分析。

六、基因筛选

筛选出细胞间高表达变异基因，这些基因在在某些细胞中高表达，而在其他细胞中低表达。这些基因有助于突出单细胞数据集中的生物信号。
筛选基因的两个标准：

• 基因在所有细胞中的平均表达量
• 基因的表达量在不同细胞之间的方差大小

####################################################
# Feature selection
# 依据：
# 1. 基因在所有细胞中的平均表达量
# 2. 基因的表达量在不同细胞之间方差大小
####################################################
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# 查看前十个高表达变异基因
top10 <- head(VariableFeatures(pbmc), 10)
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot2

all.gene <- rownames(pbmc)
# 本质为z-score，将数据变为均值为0，方差为1
pbmc <- ScaleData(pbmc, features = all.gene)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

6.1 筛选高表达变异基因

• 根据Seurat包中的 FindVariableFeatures 函数，默认使用的方法为vst

• 使用 VariableFeature函数获取高表达变异基因，为向量形式
  

• 使用VariableFeaturePlot 函数可视化并使用 LabelPoints 来标记前十个高表达变异基因，横坐标为平均表达水平，纵坐标为标准差
  

• 通过 rownames(pbmc) 来获取所有基因
  通过dim、ncol、nrow、colnames、rownames、dimnames 等函数来获取Seurat对象中子集Assay中数据的信息

• 其中dimnames函数用于给矩阵赋行名或者列明，输入为一个列表

  dimnames(mymatrix) <- list(c("row1","row2"),c("col1","col2","col3")) #为矩阵赋予行和列的名称
  #dimnames必须为list
  1
  2

  dimnames和其他函数的区别如下，
  

6.2 数据缩放

数据缩放是在PCA降维之前的必须步骤，因为PCA默认数据是服从正态分布的，使用ScaleData函数进行数据缩放,其本质是对数据进行z-score变换 : pbmc <- ScaleData(pbmc, features = all.gene)，features参数选择全部基因
z-score 标准化:

• 总体数据均值 μ
• 总体数据标准差 $\sigma$
• 代入公式： $z=\frac{x-\mu }{\sigma }$

七、降维

7.1 线性降维PCA

对筛选出来的高表达变异基因进行PCA线性降维

####################################################
# Reduction
####################################################
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
DimPlot(pbmc, reduction = "pca")
# 看拐点
ElbowPlot(pbmc, ndims = 50)
1
2
3
4
5
6
7

• DimPlot 函数图像（未聚类划分细胞类型）
  
  • ElbowPlot 函数图像：看拐点确定数据维度（下图大概为10）
    

八、细胞聚类

Seurat 使用KNN进行聚类

####################################################
# Cluster
# louvain cluster , graph based
####################################################
# 构建graph，选取前10个主成分进行聚类
pbmc <- FindNeighbors(pbmc, dims = 1:10)
# Cluster
pbmc <- FindClusters(pbmc, resolution = 0.5)
# 查看前五个细胞的聚类ID 
head(Idents(pbmc))

# 使用umap降维
pbmc <- RunUMAP(pbmc, dims = 1:10)
# 查看"FCGR3A","CD14"两个基因的降维图像
# FeaturePlot(pbmc, features = c("FCGR3A","CD14"), reduction = "umap")
DimPlot(pbmc, reduction = "umap", group.by = "seurat_clusters", label = T)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

• 聚类结果
  

• 使用Umap降维且以聚类类别分组

• 

九、寻找差异表达基因

利用 FindMarkers 和FindAllMarkers命令，可以找到各个细胞类型中与其他类别的差异表达基因，作为该细胞类型的生物学标记基因

• ident.1 : 待分析的细胞类别
• min.pct : 在两组细胞中的任何一组中检测到的最小百分比

pbmc.markers <- FindAllMarkers(pbmc, only.pos = T, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>%
  group_by(cluster) %>%
  slice_max(n = 2, order_by = avg_log2FC)
# 单独选择实验组和对照组进行差异分析
CD4.mem.DEGs <- FindMarkers(pbmc, ident.1 = "Memory CD4 T", 
                            ident.2 = "Naive CD4 T", min.pct = 0.25)
1
2
3
4
5
6
7

十、注释细胞类型

人工注释

# 根据Mark对细胞类型进行注释
FeaturePlot(pbmc, features = c("MS4A1","TYROBP","CD14","FCGR3A","FCER1A",
                               "CCR7","IL7R","PPBP","CD8A"))
# 人工注释
new.cluster.ids <- c("Naive CD4 T","CD14+ Mono","Memory CD4 T","B","CD8 T","FCGR3A+ Mono",
                     "NK","DC","Platelet")
names(new.cluster.ids) <- levels(pbmc)
# 更改细胞类型名称
pbmc <- RenameIdents(pbmc, new.cluster.ids)
# 修改Idents之后
DimPlot(pbmc, reduction = "umap", label = T)
1
2
3
4
5
6
7
8
9
10
11

• 更改Ident
  

十一、尾

11.1 流程

11.2 Seurat绘图函数

Seurat 绘图函数总结： 绘图函数总结

十二、REF

Seurat 结构 https://www.bilibili.com/video/BV1Qf4y1s7h1p=1&vd_source=4dd8c7cbe5dff7c9c8e4bdc377718621
单细胞测序原理：https://www.cnblogs.com/ZhengAbel/p/13894907.html
绘图总结：https://www.jianshu.com/p/95e61f7e834d
分析流程： https://www.jianshu.com/p/104f9290af8b
矩阵市场 ： https://www.jianshu.com/p/4d2a95ce3055