Seurat24节气之5清明---标准流程Guided Clustering Tutorial

清明:天气晴朗,草木繁茂。

Seurat官网上详细的指导完全可以满足Seurat包初级使用。不过该网站是英文的,为了方便大家迅速上手,我来走一遍标准流程。我用的是Windows 10, R4.0。 我走的流程原网站地址

首先我们需要在自己的RStudio中安装Seurat包

install.packages('Seurat')
library('Seurat')
packageVersion("Seurat")
?Seurat

原参考页面中还使用了一些相关的R包,所以我们也需要一并安装上,如果你已经安装了这些包就跳过这一步

install.packages(c('dplyr','patchwork'))

安装好R包之后,我们要Load进来现在的工作环境

library(dplyr)
library(Seurat)
library(patchwork)

示例数据可以在官网下载:https://support.10xgenomics.com/ 读入的数据可以是一个矩阵,行代表基因,列代表细胞。

1.数据导入

list.files('pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19')
pbmc.counts <- Read10X(data.dir = "pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19")

创建Seurat对象

pbmc <- CreateSeuratObject(counts = pbmc.counts)
pbmc
str(pbmc)

数据集中测到的少于200个基因的细胞(min.features = 200)和少于3个细胞覆盖的基因(min.cells = 3)被过滤掉

pbmc <- CreateSeuratObject(counts = pbmc.counts, project = "pbmc3k", min.cells = 3, min.features = 200)

2.数据质控

质控的参数主要有两个: 1.每个细胞测到的unique feature数目(unique feature代表一个细胞检测到的基因的数目,可以根据数据的质量进行调整) 2.每个细胞检测到的线粒体基因的比例,理论上线粒体基因组与核基因组相比,只占很小一部分。所以线粒体基因表达比例过高的细胞会被过滤。

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

nFeature_RNA代表每个细胞测到的基因数目,nCount代表每个细胞测到所有基因的表达量之和,percent.mt代表测到的线粒体基因的比例。

image
# FeatureScatter is typically used to visualize feature-feature relationships, but can be used
# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2

[图片上传失败...(image-5e2fdf-1615627815934)]

去除线粒体基因表达比例过高的细胞,和一些极值细胞。

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

3.标准化

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
#鉴定细胞间表达量高变的基因(feature selection)
#这一步的目的是鉴定出细胞与细胞之间表达量相差很大的基因,用于后续鉴定细胞类型,
#我们使用默认参数,即“vst”方法选取2000个高变基因。
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
image

4.细胞分类

1)分类前首先要对数据集进行降维

#Scaling the data
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
#Perform linear dimensional reduction
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Examine and visualize PCA results a few different ways
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
DimPlot(pbmc, reduction = "pca")
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
image
image
image

2)定义数据集的“维度”

这里我们需要选择出主成分的数目,用于后续细胞分类。这里定义的“维度”并不代表细胞类型的数目,而是对细胞分类时需要用到的一个参数。

#NOTE: This process can take a long time for big datasets, comment out for expediency. 
#More approximate techniques such as those implemented in ElbowPlot() can be used to reduce computation time
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:15)
ElbowPlot(pbmc)
image
image

JackStraw和Elbow都可以决定数据的“维度”。但是Elbow比较直观,我们选择Elbow结果进行解读。可以看到,主成分(PC)7到10之间,数据的标准差基本不在下降。所以我们需要在7到10之间进行选择,为了尊重官网的建议,我们选取10,即前10个主成分用于细胞的分类。

3) 细胞分类

选择不同的resolution值可以获得不同的cluster数目,值越大cluster数目越多,默认值是0.5.

pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5) 
#这里我们设置了dims = 1:10 即选取前10个主成分来分类细胞。分类的结果如下,可以看到,细胞被分为9个类别。
#Look at cluster IDs of the first 5 cells
head(Idents(pbmc), 5)

4) 可视化分类结果

TSNE和UMAP两种方法经常被用于可视化细胞类别。

#UMAP
pbmc <- RunUMAP(pbmc, dims = 1:10, label = T)
head(pbmc@reductions$umap@cell.embeddings) # 提取UMAP坐标值。
#note that you can set `label = TRUE` or use the LabelClusters function to help label individual clusters
p1 <- DimPlot(pbmc, reduction = "umap")
#T-SNE
pbmc <- RunTSNE(pbmc, dims = 1:10)
head(pbmc@reductions$tsne@cell.embeddings)
p2 <- DimPlot(pbmc, reduction = "tsne")
p1 + p2
saveRDS(pbmc, file = "pbmc_tutorial.rds")  #保存数据,用于后续个性化分析
image

5.提取各个细胞类型的marker gene

利用 FindMarkers 命令,可以找到找到各个细胞类型中与其他类别的差异表达基因,作为该细胞类型的生物学标记基因。其中ident.1参数设置待分析的细胞类别,min.pct表示该基因表达数目占该类细胞总数的比例

#find all markers of cluster 1
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
#利用 DoHeatmap 命令可以可视化marker基因的表达
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25)
?FindMarkers
top3 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 3, wt = avg_logFC)
DoHeatmap(pbmc, features = top3$gene) + NoLegend()
image

6.探索感兴趣的基因

Seurat提供了许多方法使我们能够方便的探索感兴趣的基因在各个细胞类型中的表达情况

VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
#我们能够看到,MS4A1和CD79A两个基因在细胞群体3中特异性表达。
#you can plot raw counts as well
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)

FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
#这种展示方法把基因表达量映射到UMAP结果中,同样可以直观的看到基因表达的特异性。
image
image
image

7.利用先验知识定义细胞类型

通过对比我们鉴定的marker gene与已发表的细胞类型特意的基因表达marker,可以定义我们划分出来的细胞类群。最后,给我们定义好的细胞类群加上名称

new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "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()
image

参考原文
Seurat包------标准流程

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 205,132评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,802评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,566评论 0 338
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,858评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,867评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,695评论 1 282
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,064评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,705评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 42,915评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,677评论 2 323
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,796评论 1 333
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,432评论 4 322
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,041评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,992评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,223评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,185评论 2 352
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,535评论 2 343

推荐阅读更多精彩内容