锐客网

  1. 首页 > it技术 > >

单细胞数据分析神器——Seurat

近年来，单细胞技术日益火热，并且有着愈演愈烈的趋势。在2015年至2017年，甚至对某细胞群体或组织进行单细胞测序，解析其细胞成分就能发一篇CNS级别的文章。近两三年，单细胞技术从最开始的基因组，转录组测序，发展成现在的单细胞DNA甲基化，单细胞ATAC-seq等等。测序手段也从早期的10X Genomics、 Drop-seq等，发展为现在的多种多样个性化的方法。研究内容更不仅仅局限于解析细胞群体的成分，而是向研究细胞功能和生物学特性发展。今天小编向大家简单一个实用并且易上手的单细胞数据分析软件——Seurat,大家躺在床上为国家做贡献的同时也能get新技能。
单细胞数据分析神器——Seurat
文章图片
介绍一下今天的主角，Seurat是由New York Genome Center, Satija Lab开发的单细胞数据分析集成软件包。其功能不仅包含基本的数据分析流程，如质控，细胞筛选，细胞类型鉴定，特征基因选择，差异表达分析，数据可视化等等。同时也包括一些高级功能，如时序单细胞数据分析，不同组学单细胞数据整合分析等。今天，小编以官网中提供的单细胞基因表达数据为例，为大家简单介绍一下Seurat软件包中的基础分析流程，希望能够抛砖引玉，祝大家在科研的道路上越走越远。
第一步，数据集导入
在本教程中，我们将分析从10X基因组学免费获得的外周血单个核细胞（PBMC）数据集，来源于Illumina NextSeq 500测得的2700个单细胞转录组数据。首先，我们需要把数据集存储成Seurat可识别的数据格式，
读入的数据可以是一个矩阵，行代表基因，列代表细胞。
library(dplyr)library(Seurat)# Load the PBMC datasetpbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")# Initialize the Seurat object with the raw (non-normalized data).pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)pbmc## An object of class Seurat ## 13714 features across 2700 samples within 1 assay## Active assay: RNA (13714 features)
数据导入成功以后，我们可以看到pbmc对象中包含了一个13714（基因数）X2700（细胞数）的矩阵，其实在数据导入的时候，数据集中测到的少于200个基因的细胞（min.features = 200），和少于3个细胞覆盖的基因（min.cells = 3），就已经被过滤掉了。


第二步，数据质控
质控的参数主要有两个：1.每个细胞测到的unique feature数目（unique feature代表一个细胞检测到的基因的数目，可以根据数据的质量进行调整）2.每个细胞检测到的线粒体基因的比例，理论上线粒体基因组与核基因组相比，只占很小一部分。所以线粒体基因表达比例过高的细胞会被过滤。
在质控前，我们需要目测一下数据集的基本情况。
# The [[ operator can add columns to object metadata. This is a great place to stash QC statspbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")# Visualize QC metrics as a violin plotVlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
代码运行后会得到下图所示结果，nFeature_RNA代表每个细胞测到的基因数目，nCount代表每个细胞测到所有基因的表达量之和，percent.mt代表测到的线粒体基因的比例。
单细胞数据分析神器——Seurat
文章图片
这里，我们过滤掉上图所示的离群点，并对基因的表达量进行log转换。
第三步，鉴定细胞间表达量高变的基因（feature selection）
这一步的目的是鉴定出细胞与细胞之间表达量相差很大的基因，用于后续鉴定细胞类型，我们使用默认参数，即“vst”方法选取2000个高变基因。
pbmc<- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)# Identify the 10 most highly variable genestop10<- head(VariableFeatures(pbmc), 10)# plot variable features with and without labelsplot1<- VariableFeaturePlot(pbmc)plot2<- LabelPoints(plot = plot1, points = top10, repel = TRUE)CombinePlots(plots=list(plot1, plot2))
代码运行后会得到下图所示的结果。
单细胞数据分析神器——Seurat
文章图片
第四步，细胞分类
1) 分类前首先要对数据集进行降维
#Scaling the dataall.genes <- rownames(pbmc)pbmc<- ScaleData(pbmc, features =all.genes)#Perform linear dimensional reductionpbmc<- RunPCA(pbmc, features = VariableFeatures(object = pbmc))# Examine and visualize PCA results a few different waysprint(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
## PC_ 1 ## Positive:CST3, TYROBP, LST1, AIF1, FTL ## Negative:MALAT1, LTB, IL32, IL7R, CD2 ## PC_ 2 ## Positive:CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1 ## Negative:NKG7, PRF1, CST7, GZMB, GZMA ## PC_ 3 ## Positive:HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1 ## Negative:PPBP, PF4, SDPR, SPARC, GNG11 ## PC_ 4 ## Positive:HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1 ## Negative:VIM, IL7R, S100A6, IL32, S100A8 ## PC_ 5 ## Positive:GZMB, NKG7, S100A8, FGFBP2, GNLY ## Negative:LTB, IL7R, CKB, VIM, MS4A7
官网中提供了不同的方法可视化降维的结果，此处不再赘述。感兴趣的读者可以自行探索降维的结果，代码如下。
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")DimPlot(pbmc, reduction = "pca")DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
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 timepbmc <- JackStraw(pbmc, num.replicate = 100)pbmc <- ScoreJackStraw(pbmc, dims = 1:20)JackStrawPlot(pbmc, dims = 1:15)ElbowPlot(pbmc)
JackStraw（左图）和Elbow（右图）都可以决定数据的“维度”。但是Elbow比较直观，我们选择Elbow结果进行解读。可以看到，主成分（PC）7到10之间，数据的标准差基本不在下降。所以我们需要在7到10之间进行选择，为了尊重官网的建议，我们选取10，即前10个主成分用于细胞的分类。


单细胞数据分析神器——Seurat
文章图片


3) 细胞分类
pbmc <- FindNeighbors(pbmc, dims = 1:10)pbmc <- FindClusters(pbmc, resolution = 0.5)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck## ## Number of nodes: 2638## Number of edges: 96033## ## Running Louvain algorithm...## Maximum modularity in 10 random starts: 0.8720## Number of communities: 9## Elapsed time: 0 seconds
这里我们设置了dims = 1:10 即选取前10个主成分来分类细胞。分类的结果如下,可以看到，细胞被分为9个类别。
# Look at cluster IDs of the first 5 cellshead(Idents(pbmc), 5)## AAACATACAACCAC AAACATTGAGCTAC AAACATTGATCAGC AAACCGTGCTTCCG AAACCGTGTATGCG## 1 3 1 2 6## Levels: 0 1 2 3 4 5 6 7 8
4) 可视化分类结果
TSNE和UMAP两种方法经常被用于可视化细胞类别。
# If you haven't installed UMAP, you can do so via reticulate::py_install(packages = 'umap-learn')pbmc <- RunUMAP(pbmc, dims = 1:10)# note that you can set `label = TRUE` or use the LabelClusters function to help label individual clustersDimPlot(pbmc, reduction = "umap")saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")
单细胞数据分析神器——Seurat
文章图片
第五步，提取各个细胞类型的marker gene
利用 FindMarkers 命令，可以找到找到各个细胞类型中于其他类别的差异表达基因，作为该细胞类型的生物学标记基因。其中ident.1参数设置待分析的细胞类别，min.pct表示该基因表达数目占该类细胞总数的比例
# find all markers of cluster 1cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)head(cluster1.markers, n = 5)
单细胞数据分析神器——Seurat
文章图片
利用 DoHeatmap 命令可以可视化marker基因的表达
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)DoHeatmap(pbmc, features = top10$gene) + NoLegend()
单细胞数据分析神器——Seurat
文章图片
第六步，探索感兴趣的基因
Seurat提供了小提琴图和散点图两种方法，使我们能够方便的探索感兴趣的基因在各个细胞类型中的表达情况
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
单细胞数据分析神器——Seurat
文章图片
我们能够看到，MS4A1和CD79A两个基因在细胞群体3中特异性表达。
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
单细胞数据分析神器——Seurat
文章图片
这种展示方法把基因表达量映射到UMAP结果中，同样可以直观的看到基因表达的特异性。
第七步，利用先验知识定义细胞类型
通过对比我们鉴定的marker gene与已发表的细胞类型特意的基因表达marker，可以定义我们划分出来的细胞类群。
单细胞数据分析神器——Seurat
文章图片
最后，给我们定义好的细胞类群加上名称
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()
单细胞数据分析神器——Seurat
文章图片
以上就是小编学习到的基本的单细胞转录组数据分析流程，整理出来供大家参考，本文主要参考了Seurat官网给出的单细胞转录组数据分析教程，若文中有小编解读错误之处，还请广大读者批评指正。
Seurat官网地址
【单细胞数据分析神器——Seurat】https://satijalab.org/seurat/v3.1/pbmc3k_tutorial.html

    推荐阅读


    上一篇：青少年宫六激教学体系学习1

    下一篇：中医原来是这样啊！！！