Seurat软件使用操作教程

Posted

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了Seurat软件使用操作教程相关的知识，希望对你有一定的参考价值。

参考技术A Peripheral Blood Mononuclear Cells (PBMC) 是 10X Genomics dataset page 提供的一个数据，包含2700个单细胞，出自Illumina NextSeq 500平台。

PBMCs是来自健康供体具有相对少量RNA（around 1pg RNA/cell）的原代细胞。在Illumina NextSeq 500平台，检测到2700个单细胞，每个细胞获得69000 reads。

matrix.mtx：matrix.mtx 是 MatrixMarket格式文件；更多内容见： http://math.nist.gov/MatrixMarket/formats.html

总共得到2700个细胞和13714个基因。
count matrix长什么样呢？我们可以看

其中的.表示0，即no molecules detected。当然，这个地方还有另外一种含义就是这个基因是真的没有表达。

由于单细胞测序数据中大多数的值都为0，因此，seurat使用一个稀疏矩阵来保存测序得到的count matrix，这样有利于数据存储空间的节省。

我们来看看使用稀疏矩阵和使用0来存储两种方式的大小对比。

dense为转换为0后的matrix存储大小，709591472 bytes，
sparse为.即稀疏矩阵的大小，29905192 bytes。
两者比值为23.7 倍，即使用系数矩阵来存储单细胞水平的基因表达值非常节省空间。
使用pbmc数据初始化Seurat对象

预处理主要包括基于QC指标的细胞和基因过滤，数据标准化和归一化，高变基因选择。

一般筛选条件有三个：

我们将使用以下标准进行基因过滤：
We filter cells that have unique feature counts over 2,500 or less than 200
We filter cells that have >5% mitochondrial counts
过滤前三个指标可视化

我们就出了小提请图。
查看基因数目, 线粒体基因占比与UMI数目的关系

筛选检测到基因数目超过2500或低于200的细胞
单个细胞中线粒体基因数目占比超过>5%

默认使用数据标准化方法是LogNormalize, 每个细胞总的表达量都标准化到10000，然后log取对数；结果存放于pbmc[["RNA"]]@data。
标准化前，每个细胞总的表达量

鉴定在细胞间表达高度变化的基因，后续研究需要集中于这部分基因。Seurat内置的FindVariableFeatures()函数，首先计算每一个基因的均值和方差，并且直接模拟其关系。默认返回2000个基因。

线性转换缩放数据，ScaleData()函数可以实现此功能。

最终每个基因均值为0，方差为1。

结果存放于pbmc[["RNA"]]@scale.data。

设置参数features是因为ScaleData默认处理前面鉴定的差异基因。这一步怎么做都不会影响到后续pca和聚类，但是会影响做热图。

移除影响方差的因素。

对缩放后的数据进行PCA分析，默认使用前面鉴定表达变化大的基因。使用features参数可以重新定义数据集。

VizDimReduction, DimPlot, 和 DimHeatmap可以从基因或细胞角度可视化pca结果
查看对每个主成分影响比较大的基因集

查看前五的可变基因

可视化对每个主成分影响比较大的基因集

两个主成分的展示

DimHeatmap绘制基于单个主成分的热图，细胞和基因的排序都是基于他们的主成分分数。对于数据异质性的探索是很有帮助的，可以帮助用户选择用于下游分析的主成分维度。

为了避免单个基因影响，Seurat聚类细胞时使用pca结果。首先需要确定的是使用多少个主成分用于后续分析。常用有两种方法，一种是基于零分布的统计检验方法，这种方法耗时且可能不会返回明确结果。另一种是主成分分析常用的启发式评估。
JackStraw()
在JackStraw()函数中, 使用基于零分布的置换检验方法。随机抽取一部分基因(默认1%)然后进行pca分析得到pca分数，将这部分基因的pca分数与先前计算的pca分数进行比较得到显著性p-Value，。根据主成分（pc）所包含基因的p-value进行判断选择主成分。最终的结果是每个基因与每个主成分的关联的p-Value。保留下来的主成分是那些富集小的p-Value基因的主成分。
处理大数据时会花费大量时间；ElbowPlot()内置了一些其它的方法可以减少运行时间。

JackStrawPlot()函数提供可视化方法，用于比较每一个主成分的p-value的分布，虚线是均匀分布；显著的主成分富集有小p-Value基因，实线位于虚线左上方。下图表明保留10个pca主成分用于后续分析是比较合理的。

ElbowPlot

启发式评估方法生成一个Elbow plot图。在图中展示了每个主成分对数据方差的解释情况（百分比表示），并进行排序。根据自己需要选择主成分，图中发现第9个主成分是一个拐点，后续的主成分(PC)变化都不大了。

注意：鉴别数据的真实维度不是件容易的事情；除了上面两种方法，Serat官当文档还建议将主成分(数据异质性的相关来源有关)与GSEA分析相结合。Dendritic cell 和 NK aficionados可能识别的基因与主成分 12 和 13相关，定义了罕见的免疫亚群 (i.e. MZB1 is a marker for plasmacytoid DCs)。如果不是事先知道的情况下，很难发现这些问题。

Serat官当文档因此鼓励用户使用不同数量的PC（10、15，甚至50）重复下游分析。其实也将观察到的，结果通常没有显著差异。因此，在选择此参数时，可以尽量选大一点的维度，维度太小的话对结果会产生不好的影响。

Seurat v3应用基于图形的聚类方法，例如KNN方法。具有相似基因表达模式的细胞之间绘制边缘，然后将他们划分为一个内联群体。

在PhenoGraph中，首先基于pca维度中（先前计算的pca数据）计算欧式距离（the euclidean distance），然后根据两个细胞在局部的重合情况（Jaccard 相似系数）优化两个细胞之间的边缘权值。此步骤内置于FindNeighbors函数，输入时先前确定的pc数据。

为了聚类细胞，接下来应用模块化优化技术迭代将细胞聚集在一起。（the Louvain algorithm (default) or SLM [SLM, Blondel et al., Journal of Statistical Mechanics]），FindClusters函数实现这一功能，其中需要注意resolution参数，该参数设置下游聚类分析的“granularity”，更大的resolution会导致跟多的细胞类群。3000左右的细胞量，设置resolution为0.4-1.2是比较合适的。细胞数据集越大，需要更大的resolution参数, 会获得更多的细胞聚类。
查看细胞属于那个类群可以使用函数Idents。

添加细胞标签

此时可以保存数据，方便下次直接导入数据修改图形。

Seurat可以通过差异表达分析寻找不同细胞类群的标记基因。FindMarkers函数可以进行此操作，但是默认寻找单个类群(参数ident.1)与其他所有类群阳性和阴性标记基因。FindAllMarkers函数会自动寻找每个类群和其他每个类群之间的标记基因。

min.pct参数:设定在两个细胞群中任何一个被检测到的百分比，通过此设定不检测很少表达基因来缩短程序运行时间。默认0.1

thresh.test参数：设定在两个细胞群中基因差异表达量。可以设置为0 ，程序运行时间会更长。

max.cells.per.ident参数：每个类群细胞抽样设置；也可以缩短程序运行时间。

Seurat有多种方法可视化标记基因的方法
VlnPlot: 基于细胞类群的基因表达概率分布
FeaturePlot：在tSNE 或 PCA图中画出基因表达情况
RidgePlot，CellScatter，DotPlot

Seurat有多种方法可视化标记基因的方法
VlnPlot: 基于细胞类群的基因表达概率分布
FeaturePlot：在tSNE 或 PCA图中画出基因表达情况
RidgePlot，CellScatter，DotPlot

DoHeatmap为指定的细胞和基因花表达热图。每个类群默认展示top 20标记基因。

Cluster ID Markers Cell Type
0 IL7R, CCR7 Naive CD4+ T
1 IL7R, S100A4 Memory CD4+
2 CD14, LYZ CD14+ Mono
3 MS4A1 B
4 CD8A CD8+ T
5 FCGR3A, MS4A7 FCGR3A+ Mono
6 GNLY, NKG7 NK
7 FCER1A, CST3 DC
8 PPBP Platelet

12 最后保存这个运行过得文件

自学Seurat的一点记录

教程来自ANALYSIS OF SINGLE CELL RNA-SEQ DATAANALYSIS OF SINGLE CELL RNA-SEQ DATAhttps://broadinstitute.github.io/KrumlovSingleCellWorkshop2020/index.html

本文主要针对Seurat官网教程上一些难懂的部分进行分享。

 scRNA-Seq和bulk RNA-Seq两个关键的不同是drop-out和potential for quality control (QC) metrics，前者是指单细胞中有限的RNA探测导致很多的0值，也正因如此，使用sparse matrix而不是dense matrix进行存储，这样节省存储空间，因为sparse matrix默认大多数为0（用 . 表示），只存储非零值。

识别分析细胞亚群之前，质量控制非常重要

counts <- Read10X(data.dir = counts_matrix_filename)  # Seurat function to read in 10x count data

 将文件读入后得到sparse Matrix of class "dgCMatrix"

rownames是不同的gene名（feature），colnames是不同的16 nt 长barcode序列（细胞），可以通过dim(counts)等函数查看matrix的行列数，即总细胞数和基因数。

过滤低质量的细胞

通过 Matrix::colSums(counts) 计算每一列的count和，即counts per cell

通过 Matrix::rowSums(counts) 计算每一行的count和，即counts per gene

而 Matrix::colSums(counts > 0)则是分别计算每一列中counts大于0的有多少行，即细胞中有多少个基因被表达了，colSums 和 rowSums函数返回的是向量vector形式。

通过探测的基因数量（文库复杂性）对细胞排序

通过plot图可以观察到failed libraries (lower end outliers) 或cell doublets (higher end outliers).

plot(sort(genes_per_cell), xlab='cell', log='y', main='genes per cell (ordered)')

使用Seurat beginning

  • 创建Seurat对象 
    seurat <- CreateSeuratObject(counts = counts, min.cells = 3, min.features = 350, project = "10X_NSCLC")

    留下在3个或更多细胞中表达的基因，表达350个或更多基因的细胞。

seurat对象包括多个slots ，不仅存储原始的count数据，也存储计算后的结果，会随着分析不断更新。如meta.data，是一个数据框，原始变量包括orig.ident， nCount_RNA， nFeature_RNA，可以用$或[引用。

ps：orig.ident指cell identity，细胞身份，此时它的identity是指10X_NSCLC，但在细胞聚类后，则指细胞所属的cluster

  • 获取对象信息
GetAssayData(object = seurat, slot = 'data')[1:10, 1:10]
  • 获取seurat包含的method  列出seurat的所有函数 
    utils::methods(class = 'Seurat')
ls("package:Seurat")

    预处理1：过滤掉低质量细胞

  • 除去被破坏的受伤细胞，当细胞应急凋亡，会导致线粒体泄露和RNA降解，因此会有线粒体来源基因的富集，计算每个细胞中线粒体转录本的比例percent.mt，作小提琴图。
# The [[ operator can add columns to object metadata. This is a great place to stash QC

seurat[["percent.mt"]] <- PercentageFeatureSet(object = seurat, pattern = "^MT-")
VlnPlot(object = seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
  • 管家基因：细胞中管家基因的表达数量可以反映细胞常规过程的普遍活跃水平，风度高，稳定表达，对dropout不敏感。计算每个细胞中有多少管家基因表达。
# Load the the list of house keeping genes
hkgenes <- read.table(paste0(dirname, "housekeepers.txt"), skip = 2)
hkgenes <- as.vector(hkgenes$V1)

# remove hkgenes that were not found
hkgenes.found <- which(toupper(rownames(seurat@assays$RNA@data)) %in% hkgenes)

n.expressed.hkgenes <- Matrix::colSums(seurat@assays$RNA@data[hkgenes.found, ] > 0)
seurat <- AddMetaData(object = seurat, metadata = n.expressed.hkgenes, col.name = "n.exp.hkgenes")
VlnPlot(object = seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "n.exp.hkgenes"), ncol = 2)
  • 过滤 
    seurat <- SubsetData(object = seurat, subset.names = c("nFeature_RNA", "percent.mito","n.exp.hkgenes"), low.thresholds = c(350, -Inf,55), high.thresholds = c(5000, 0.1, Inf))

seurat <- subset(seurat, nFeature_RNA < 5000)
seurat <- subset(seurat, nFeature_RNA > 350)
seuart <- subset(seurat, n.exp.hkgenes > 55)
seuart <- subset(seurat, percent.mt < 10)

    预处理2：表达标准化

  • LogNormalize：用细胞中该基因的counts数除以细胞中的counts总数，再乘scale因子，默认10的4次方，再进行对数转换。可以稳定方差 
    seurat <- NormalizeData(object = seurat, normalization.method = "LogNormalize", scale.factor = 1e4)

    差异基因的检测

  • 计算每个基因在细胞中的平均表达量和细胞间离散程度z-score 
    seurat <- FindVariableFeatures(object = seurat, selection.method = "vst", nfeatures = 2000)

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

seurat <- FindVariableFeatures(object = seurat, selection.method = "vst", nfeatures = 2000, 
                               mean.function = ExpMean, dispersion.function = LogVMR, 
                               num.bin = 40, 
                               mean.cutoff = c(0.0125, 1),
                               dispersion.cutoff = c(0, 0.5))

## Check number of variable genes
length(seurat@assays$RNA@var.features)

    selection.method    
    How to choose top variable features. Choose one of :

    vst: First, fits a line to the relationship of log(variance) and log(mean) using local polynomial regression (loess). Then standardizes the feature values using the observed mean and expected variance (given by the fitted line). Feature variance is then calculated on the standardized values after clipping to a maximum (see clip.max parameter).

    mean.var.plot (mvp): First, uses a function to calculate average expression (mean.function) and dispersion (dispersion.function) for each feature. Next, divides features into num.bin (deafult 20) bins based on their average expression, and calculates z-scores for dispersion within each bin. The purpose of this is to identify variable features while controlling for the strong relationship between variability and average expression.

    dispersion (disp): selects the genes with the highest dispersion values

关于selection的3种方法，下期再叙！ 

以上是关于Seurat软件使用操作教程的主要内容，如果未能解决你的问题，请参考以下文章

Seurat 3.0 实例教程

Seurat4.0系列教程21： 结合Cell Hashing分析双细胞

自学Seurat的一点记录

2019-10-22 R语言Seurat包下游分析-1

【单细胞测序数据分析-1】认识Seurat对象数据结构/数据格式及操作

Seurat中对细胞分群（Cluster）的操作