【原】单细胞初探（seurat基础流程）（2021公开课配套笔记）

健明 2021-07-14

展开全文

新课发布在B站了，马上有热心的粉丝看完后写了配套笔记：

下面是粉丝linbo的笔记投稿

前言

自学生信半载有余，跌跌撞撞，不敢和大佬同称萌新，勉强算得上菜鸡。根据课题组进展，马上要接手一个单细胞课题，适逢jimmy老师在B站发布了《单细胞公开课2021》，结合课程以及公众号推文，上手了一下单细胞基础分析流程。
先说一下感悟吧：安装R包最为痛苦😂😂😂，相同代码在不同版本的R包上运行结果可能不同。

以Seruat官方教程和数据为练习，地址：https:///seurat/articles/pbmc3k_tutorial.html，那下面我们就开始吧。

首先在官网下载原始数据（方法如下）

下载解压后得到如下3个文件，这3个文件就是分析要用的原始数据。

基础分析流程

参考：可视化单细胞亚群的标记基因的5个方法

读取数据

这里非常感谢健明老师让我学会了使用project管理R代码，非常方便。将解压的文件拷贝到project根目录的data目录下。

## 魔幻清空
rm(list=ls()) 
## 加载R包
library(dplyr)
library(Seurat)
library(patchwork)
## 读取数据
pbmc.data <- Read10X(data.dir = "data/pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19")
## 创建Seruat对象
pbmc <- CreateSeuratObject(counts = pbmc.data, 
                           project = "pbmc3k", 
                           min.cells = 3, # min.cell：每个feature至少在多少个细胞中表达(feature=gene)
                           min.features = 200) # min.features：每个细胞中至少有多少个feature被检测到

这个 pbmc 就是全部单细胞数据分析所需要的，查看变量pbmc 如下：

pbmc
# An object of class Seurat 
# 13714 features across 2700 samples within 1 assay 
# Active assay: RNA (13714 features, 0 variable features)

可见一共有13714个基因，2700个细胞。有了这个构建好的Seurat 对象，也可参考官网代码对数据进行探索（此处不再赘述），后续分析代码相对来说非常标准，如下:

数据质控

通常使用的QC指标：

每个细胞被检测到的unique genes数量（低质量的细胞或空的液滴含有unique genes较少；细胞双重态或多重态检测到异常高unique genes）；
细胞内被检测到的features的总数目（与unique genes高度相关）；
线粒体基因占比（低质量/将要死去的细胞经常出现线粒体基因占比过高）。

接下来，根据每个细胞线粒体基因占比、检测到的基因数进行简单的过滤。先计算细胞内线粒体基因占比，类似的核糖体基因(大多为RP开头)也可以这样计算，但是注意线粒体基因的MT-，其中-不能少，不然容易把别的基因也计算在内。

## 计算细胞中线粒体基因比例
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
## 使用小提琴图可视化QC指标
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

## FeatureScatter于可视化特征-特征关系
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2

## 过滤
## 官方推荐过滤掉独特特征计数超过2500或小于200的细胞，或者过滤掉超过5%线粒体基因比例的细胞
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc
# An object of class Seurat 
# 13714 features across 2638 samples within 1 assay 
# Active assay: RNA (13714 features, 0 variable features)

可见过滤掉了62个细胞，继续后续分析

数据标准化

各种组学的原始数据普遍存在数据量不统一，数据变化范围过大，数据变化幅度不统一等问题。各种测序数据的分析流程都要对原始数据进行“标准化”，以符合下游分析的需求，单细胞数据也不例外。

## 标准化
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 1e4)

识别高变基因

高变基因：在一些细胞中表达高，另一些细胞中表达低的基因。单细胞表达矩阵为稀疏矩阵（很多0，且为了压缩文件大小，0用.表示），选高变基因可以找到包含信息最多的基因

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000) #返回两千个高变基因

## 提取前10的高变基因
top10 <- head(VariableFeatures(pbmc), 10)
top10
# [1] "PPBP"   "LYZ"    "S100A9" "IGLL5"  "GNLY"   "FTL"    "PF4"    "FTH1"   "GNG11"  "S100A8"

## 展示高变基因
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

数据归一化

数据归一化是将每个基因在所有细胞中的均值变为0，方差标为1，即零-均值归一化（z-score归一化）。是组学数据常用的一种降维前处理方法。

pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")

降维

PCA降维，默认使用前面2000个高变基因的scale矩阵用于降维。Run开头的函数降维，Find开头的函数聚类。resolution参数表达聚类的分辨率，值越大得到的cluster越多，对于3K细胞的单细胞数据0.4-1.2 通常会得到较好的结果。

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc)) ##默认会输出5个主成分
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)

## 查看前5分析细胞聚类数ID
head(Idents(pbmc), 5)
# AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1 
#                3                2                8                4                3 
# Levels: 0 1 2 3 4 5 6 7 8 9

## 查看每个类别多少个细胞
table(pbmc@meta.data$seurat_clusters)
#   0   1   2   3   4   5   6   7   8   9 
# 598 483 345 324 300 214 159 106  96  13

上述结果可见，分为10类，由0-9，细胞数依次递减。

UMAP/tSNE可视化

降维可以将相似的细胞放置在低维度的空间

pbmc <- RunUMAP(pbmc, dims = 1:10)
p1 <- DimPlot(pbmc, reduction = "umap", label = T, label.size = 5)
pbmc <- RunTSNE(pbmc, dims = 1:10)
p2 <- DimPlot(pbmc, reduction = "tsne", label = T, label.size = 5)
p1+p2

在umap图中，cluster之间的距离更明显，图形也更美观。

Marker鉴定细胞

根据生物学背景知识，将表达相应Marker基因的各个单细胞亚群的重命名，官网给出的marker基因和细胞对应关系如下所示，本人该部分生物学知识较浅，暂不深究，待后续深入了解。

VlnPlot(pbmc, features = c("MS4A1", "CD79A")) #查看B细胞marker基因在不同聚类中表达情况

可视化展示marker基因的坐标映射分布

FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", 
                               "CD14", "FCER1A", "FCGR3A", 
                               "LYZ", "PPBP", "CD8A"))

上图结合UMAP可视化中的聚类对聚类重命名，会得到如下所示的分群情况。

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)
DimPlot(pbmc, reduction = 'umap', 
        label = TRUE, pt.size = 0.5) + NoLegend()

这个时候有5个可视化方法选择，分别是：小提琴图，坐标映射图，峰峦图，气泡图，热图。代码参见可视化单细胞亚群的标记基因的5个方法

## 小提琴图
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))

## 坐标映射图
FeaturePlot(pbmc, features = c("MS4A1", "CD79A"))

## 峰峦图
RidgePlot(pbmc, features = c("MS4A1", "CD79A"), ncol = 1)

features= c('IL7R', 'CCR7','CD14', 'LYZ',  'IL7R', 'S100A4',"MS4A1", "CD8A",'FOXP3',
            'FCGR3A', 'MS4A7', 'GNLY', 'NKG7',
            'FCER1A', 'CST3','PPBP')
## 气泡图
DotPlot(pbmc, features = unique(features)) + RotatedAxis()

我个人比较喜欢这个气泡图，看的很清楚，也可以根据气泡图来重命名聚类名。

## 热图
DoHeatmap(subset(pbmc ), features = features, size = 3)

save(pbmc,file = 'basic.sce.pbmc.Rdata') #保存用于后续分析

至此，单细胞基础流程分析就结束了。接下来根据健明老师解答的问题，来加深一下理解。特别是seurat对象的结构，在R语言中同一变量重复赋值，则会覆盖，而创建的Seruat对象pbmc则不会。查阅资料才知道Seruat对象是S4结构，会记录所执行的计算及其信息。在此献上周运来老师总结的一幅Seruat对象结构图。