背景

关注生信方向的朋友都知道，目前单细胞各个细分领域比较热门。相较于常规的Bulk测序，单细胞能发现细胞层面的异质性。觉得很有前景，故打算整理半年前参加线上培训的内容，跟大家分享单细胞的常规分析思路，出一个小系列。文章最后会贴出线上培训的视频链接。

单细胞扫盲

数据获取

单细胞测序目前价格较贵，测一个样品需要好几万，入门可以从公共数据库下载，以GSE159677为例说明常规数据获取渠道。本教程实际数据见文末“原始数据及Code”节的百度网盘。

• 进入NCBI官网
• 选择GEO DataSets
• 输入从文献中得到的数据序号：GSE159677
• 选择包含该数据的文献
• 找到Supplementary file中的链接进行下载

• 文件下载成功后依次解压
• 发现该数据有6个样品，选择其中一个解压之后可得三个文件
• 该文件就是处理的起始文件

文件解读

上述解压得到的三个文件分别为barcodes(区分各个细胞)、features（或genes，包含基因信息）、matrix（包含表达量信息，通过坐标定位具体表达量归属）。根据三个文件整合可以得到基因的表达矩阵。可能有的朋友会疑惑，直接用一个文件直接保存表达矩阵不就可以了，为什么需要三个文件？这是因为单细胞的表达矩阵是大部分为0的稀疏矩阵，分为三个文件可以节省内存。

genes.tsv文件（有时也叫features.tsv文件）

文件内容：有两列，第一列为基因ID，第二列为基因Symbol ID。

barcodes.tsv文件

文件内容：有一列，内容为测序时为了区分各个细胞的标记信息，称为Barcodes

matrix.mtx

文件内容：有三列，数字的第一行是测序的汇总信息。第一行的第一个为测序的总基因数（说明genes.tsv文件中有32738行），第二个为测序的总细胞数（说明barcodes.tsv文件中有2700行），第三个为测序的总reads数。 除第一行，其余数字记录矩阵信息，前两列代表坐标，第三列表示具体某个细胞某个基因的表达量。

表达量矩阵

根据上述三个文件可以得到表达量的稀疏矩阵。如右下角所示。

过程详解

读取表达矩阵

library(Seurat)                                #加载Seurat
library(dplyr)                                 #加载dplyr包
library(SingleR)                               #加载SingleR包，细胞注释包
library(celldex)                               #加载celldex包

#读取矩阵文件
pbmc.counts <- Read10X(data.dir = "hg19/")    #读取预处理的表达文件,需要根据自己情况设置文件位置
pbmc <- CreateSeuratObject(counts = pbmc.counts,min.cells = 3,min.features = 200)                    #创建Seurat对象

注：

• 脚本与“hg19”处于同级目录下，上述Read10X（）函数中参数为相对路径，也可添加绝对路径（即完整路径）。“hg19”中为genes/barcodes/matrix文件。
• CreateSeuratObject（）函数参数意义：

  • min.cells每个基因至少在3个细胞中表达（即，表达量不为0）。
  • min.features每个细胞中至少包含200个基因才会保留。

细胞过滤

#计算线粒体基因表达比例
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")                #计算线粒体基因表达比例
# 注意MT大小写，小写表示另一个物种

VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)   #绘制Feature,Count,mt比例图

#细胞过滤
pbmc <- subset(pbmc,subset = nFeature_RNA <=2000 & nCount_RNA <= 5000 & percent.mt <= 5) #过滤Feature,Count,mt
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),group.by = 'orig.ident', ncol = 3)        #绘制过滤后的图片

注：

• PercentageFeatureSet（）参数意义:

  • pattern = "^MT-"为正则表达式写法，意为匹配"MT-"开头的基因，即线粒体基因。区分MT大小写，此处如果小写则代表另一个物种。
  • 为什么将线粒体基因作为一个过滤条件？“因为线粒体是参与细胞凋亡启动和执行的主要细胞器之一。线粒体基因高表达意为着样品质量差，导致大量细胞凋亡或裂解。以及特定样本的生物学，例如肿瘤活检，可能由于代谢活动和/或坏死而增加线粒体基因表达。”
• subset（）参数意义：

  • nFeature_RNA <=2000 & nCount_RNA <= 5000 & percent.mt <= 5：根据过滤前小提琴图所示，nFeature_RNA（基因个数），nCount_RNA（count个数），percent.mt（染色体基因比例）主要分别在2000,5000,5以下，故作为相应阈值。
• VlnPlot（）参数意义：

  • ncol=3 代表小提琴图每行3个图

数据标准化

#数据标准化（Normalizing the data）

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

• 数据标准化的意义：

  • 去除测序深度带来的影响：
  • 例：
  • 细胞A中总reads数10，基因a的reads数为2；细胞B中总reads数1000，基因a的reads数为20。并不能说基因a在细胞B中表达量大于A，故需要进行标准化。
  • 标准化原则：
  • 每个细胞的每个基因的count数除以该细胞总count数，然后乘以因子（10000），再进行log(n+1)转换
  • 标准化后的数据保存在pbmc@assays RNA@counts中

识别高变基因

#识别高变基因
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)


# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(pbmc), 10)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)

plot1 + plot2

• FindVariableFeatures（）参数意义：

  • FindVariableFeatures 函数有 3 种选择高表达变异基因的方法，可以通过 selection.method参数来选择，它们分别是： vst（默认值）， mean.var.plot 和 dispersion。 nfeatures 参数的默认值是 2000，可以改变。如果 selection.method 参数选择的是 mean.var.plot，就不需要人为规定高表达变异基因的数目，算法会自动选择合适的数目。 建议使用完 FindVariableFeatures 函数后，用 VariableFeaturePlot 对这些高表达变异基因再做个可视化，看看使用默认值 2000 会不会有问题。

降维聚类

#数据标准化 Scale

#对所有基因进行标准化
# all.genes <- rownames(pbmc)
# pbmc <- ScaleData(pbmc, features = all.genes)    

#只对上述2000个高变基因进行标准化，通常仅对高变基因进行标准化和降维
pbmc <- ScaleData(pbmc)                          

#PCA降维
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc) )

DimPlot(pbmc, reduction = "pca")
ElbowPlot(pbmc)

str(pbmc)

# 细胞聚类
# 有UMAP和tSNE聚类两种
# 都是基于PCA降维的结果进行聚类
# tSNE聚类与UMAP聚类选择一个执行即可
# ！！！此处维度选择10是根据上述“ElbowPlot(pbmc)”生成的主成分图确定的，图中拐点为10，表示10个维度足以反映高维信息
pbmc <- FindNeighbors(pbmc, dims = 1:10)        
pbmc <- FindClusters(pbmc, resolution = 0.5) # 分辨率与最终得到的cluster数量成正比

#UMAP降维
pbmc <- RunUMAP(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap")

#tSNE降维
#pbmc <- RunTSNE(pbmc, dims = 1:10)
#DimPlot(pbmc, reduction = "tsne")

• ScaleData()标准化函数作用：

  • 为后续PCA降维做准备。
  • PCA降维要求数据为正态分布，即平均值为0，方差为1。

• DimPlot（）函数生成主成分分析结果图1

• ElbowPlot（）函数生成主成分分析结果图2

  • 发现该图中拐点在10附近，故确定10作为最终主成分。

• FindNeighbors（）参数意义：

  • dims = 1:10，此处的维度由上述主成分分析2图得到。
• FindClusters() 参数意义：

  • resolution = 0.5，此参数决定了后续所得Clusters的数目，分辨率与最终得到的cluster数量成正比。
• 细胞聚类选择UMAP与tSNE其中一个执行即可，后续以UMAP举例。

  • 当前只能知道所有细胞被分为0~7,8种类型细胞，但是不知道每个Cluster具体属于哪个生物类别，需要经过后续注释才可获悉。

寻找差异表达基因

pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.1, logfc.threshold = 0.25)

head(pbmc.markers)

#每个cluster的top10基因
top10 <- pbmc.markers %>%
  group_by(cluster) %>%
  top_n(n = 10, wt = avg_log2FC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()

• FindAllMarkers（）参数意义：

  • only.pos = TRUE：只寻找上调的基因
  • min.pct = 0.1：某基因在细胞中表达的细胞数占相应cluster细胞数最低10%
  • logfc.threshold = 0.25 ：fold change倍数为0.25
  • 该函数是决速步，执行比较耗时。
• 将每个Cluster中Top10的差异基因以热图的形式展现出来

  • 横坐标表示基因，纵坐标代表cluster, cluster由细胞组成，热图颜色表示基因表达量

原始数据及code

原始数据

注：

• 本教程使用的数据不是开头举例的GSE159677，而是下述网盘中的测试数据，更小，利于学习。
• 链接：https://pan.baidu.com/s/1Yi8M3rnKpZzJiiEcUwTqMw
• 提取码：h5td

All_Code

library(Seurat)                                #加载Seurat
library(dplyr)                                 #加载dplyr包
library(SingleR)                               #加载SingleR包，细胞注释包
library(celldex)                               #加载celldex包

#读取矩阵文件
#读取预处理的表达文件,需要根据自己情况设置文件位置
pbmc.counts <- Read10X(data.dir = "hg19/")

# 创建Seurat对象
pbmc <- CreateSeuratObject(counts = pbmc.counts,min.cells = 3,min.features = 200)           

#获取Metadata数据
head(pbmc@meta.data)   #得到metadata信息

#计算线粒体基因表达比例
# 注意MT大小写，小写表示另一个物种
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

#绘制Feature,Count,mt比例图
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)   

#细胞过滤
pbmc <- subset(pbmc,subset = nFeature_RNA <=2000 & nCount_RNA <= 5000 & percent.mt <= 5) 

#绘制过滤后的图片
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),group.by = 'orig.ident', ncol = 3)

#数据标准化（Normalizing the data）
# 去除测序深度带来的影响
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

#识别高变基因
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(pbmc), 10)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)

plot1 + plot2


#数据标准化 Scale
# PCA降维要求数据为正态分布，即平均值为0，方差为1。

#对所有基因进行标准化
# all.genes <- rownames(pbmc)
# pbmc <- ScaleData(pbmc, features = all.genes)

#只对上述2000个高变基因进行标准化，通常仅对高变基因进行标准化和降维
pbmc <- ScaleData(pbmc)                       

#PCA降维
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc) )

DimPlot(pbmc, reduction = "pca")
ElbowPlot(pbmc)

str(pbmc)

# 细胞聚类
pbmc <- FindNeighbors(pbmc, dims = 1:10)        
pbmc <- FindClusters(pbmc, resolution = 0.5) # 分辨率与最终得到的cluster数量成正比

#UMAP降维
pbmc <- RunUMAP(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap")

#tSNE降维
pbmc <- RunTSNE(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "tsne")

#寻找差异表达基因
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.1, logfc.threshold = 0.25)

#每个cluster的top10基因
top10 <- pbmc.markers %>%
  group_by(cluster) %>%
  top_n(n = 10, wt = avg_log2FC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()

其他

• 后续将更新细胞注释及通路富集分析相关内容，敬请期待。

• 该单细胞系列教程参考视频，详见链接 单细胞教程。

• 本人也是在逐步学习单细胞相关内容，上述为培训后根据个人理解整理所得，如有问题欢迎评论区讨论。

• 推文多平台同步发布，公众号内容食用更佳

• 更多内容，请关注微信公众号“生信矿工”

• 如有意见或建议可以在评论区讨论