单细胞初探(seurat基础流程)(2021公开课配套笔记)
新课发布在B站了,马上有热心的粉丝看完后写了配套笔记:
前言
自学生信半载有余,跌跌撞撞,不敢和大佬同称萌新,勉强算得上菜鸡。根据课题组进展,马上要接手一个单细胞课题,适逢jimmy老师在B站发布了《单细胞公开课2021》,结合课程以及公众号推文,上手了一下单细胞基础分析流程。
先说一下感悟吧:安装R包最为痛苦😂😂😂,相同代码在不同版本的R包上运行结果可能不同。
以Seruat官方教程和数据为练习,地址:https://satijalab.org/seurat/articles/pbmc3k_tutorial.html, 那下面我们就开始吧。
首先在官网下载原始数据(方法如下)
下载解压后得到如下3个文件,这3个文件就是分析要用的原始数据。
基础分析流程
读取数据
这里非常感谢健明老师让我学会了使用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对象结构图。
文末友情推荐
与十万人一起学生信,你值得拥有下面的学习班: