2019-2-1 来自https://satijalab.org/seurat/
1.启动Seurat
library(Seurat)
library(dplyr)
##上传数据集
pbmc.data <- Read10X(data.dir = "~/Downloads/filtered_gene_bc_matrices/hg19/")
##检测常规和稀疏矩阵的内存。
dense.size <- object.size(x = as.matrix(x = pbmc.data))
dense.size
sparse.size <- object.size(x = pbmc.data)
sparse.size
dense.size/sparse.size
##用rawdata创建object;去除少于3个细胞表达的基因(~0.1% of the data),去除少于200基因的细胞
pbmc <- CreateSeuratObject(raw.data = pbmc.data, min.cells = 3, min.genes = 200,
project = "10X_PBMC")
2.标准预处理流程
包括新建object,基于质量控制进行细胞的筛选和滤过,数据标准化和缩放,以及高度可变基因的检测。Seurat v2之前的版本将这一过程都囊括到了Set up这一步。v2版本则将这一步分开。
3.质控和筛选后续分析的细胞
CreatSeuratObject设置了最小基因数量截点的同时,可以根据技术或者生物学参数滤过细胞。Seurat包可以进行QC标准的简单探索和基于使用者设定标准过滤细胞。
Seurat主要是通过作图剔除outlier:
可视化gene coutns和molecule counts,绘制散点图,除去基因数异常点(潜在multiplets);还可以基于线粒体基因比例除去细胞。
##自动计算genes和UMIs的数量(nGene /nUMI表示) ;无UMI的数据,nUMI代表细胞内非标准化数值的总数
##计算线粒体基因的比例,用AddMetaData保存为percent.mito
##用object@raw.data代表未经转化和log标准化的counts
##定位为线粒体基因的UMI%为常见的QC标准。
具体command:
mito.genes <- grep(pattern = "^MT-", x = rownames(x = pbmc@data), value = TRUE)
percent.mito <- Matrix::colSums(pbmc@raw.data[mito.genes, ])/Matrix::colSums(pbmc@raw.data)
# AddMetaData 增加栏目到object@meta.data,存储QC stat值
pbmc <- AddMetaData(object = pbmc, metadata = percent.mito, col.name = "percent.mito")
#绘制小提琴图
VlnPlot(object = pbmc, features.plot = c("nGene", "nUMI", "percent.mito"), nCol = 3)
#绘制GenePlot
par(mfrow = c(1, 2))
GenePlot(object = pbmc, gene1 = "nUMI", gene2 = "percent.mito")
GenePlot(object = pbmc, gene1 = "nUMI", gene2 = "nGene")
#过滤细胞:除去>2500和少于200基因数的细胞
#用low.thresholds和high.thresholds设立阈值,如果不想要更低或者更高的阈值,可以用-Inf和Inf
pbmc <- FilterCells(object = pbmc, subset.names = c("nGene", "percent.mito"), low.thresholds = c(200, -Inf), high.thresholds = c(2500, 0.05))
4.标准化数据
筛选数据之后下一步是标准化数据。默认使用“LogNormalize”,乘以scale factor(默认10000),在进行log转换。
pbmc <- NormalizeData(object = pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
5.检测单个细胞中高度可变基因
目的:检测高度可变基因并用于下游分析。
用FindVariableGenes计算平均表达量和离差,将这些基因放入bins,计算bin的z-score。可以用于控制变异和平均表达量之间的关系。
官方建议使用者设置这些参数标记dispersion plot上可视化的异常点,但具体的参数设置必须基于数据类型、样本异质性和标准化策略。
这里的参数可以鉴定~2000可变基因,是标准化到1e4 分子数的UMI数据的经典参数设置。
pbmc <- FindVariableGenes(object = pbmc, mean.function = ExpMean, dispersion.function = LogVMR, x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5)
length(x = pbmc@var.genes)
6.数据缩放和除去不要变异来源
常见混杂的变异来源有:技术噪音,批次效应,生物学来源的变异(如细胞周期)。
Buettneret al, NBT, 2015, 认为剔除这些信号可以改善下游降维分析和聚类分析效果(https://www.nature.com/articles/nbt.3102)
Seurat基于使用者给定的变量建立线性模型从而预测基因表达。这些模型缩放的z-score残差存储在scale.data slot,用于降维和聚类分析。可以回归批次效应、cell alignment rate(针对Drop-seq data)检测分子数以及线粒体基因表达带来的细胞-细胞变异。
对于细胞周期细胞,可以学习“cell-cycle”评分,进行回归化。
v2版本将regression作为scale的一部分,所以原来的RegressOut弃用,代为ScaleData的vars.to.regress。
pbmc <- ScaleData(object = pbmc, vars.to.regress = c("nUMI", "percent.mito"))