Seurat -- variable features select

brief

• seurat提供的教学里面包含了Standard pre-processing workflow,workflow包括QC，normalization，scale data ，detection of highly variable features。
• 其中 normalization就有蛮多方法的，seurat自己就提供了两种，一种是"LogNormalization"，另一种是目前正在大力推荐的"SCTransform",而且每一种方法参数众多，我觉得有必要对其进行仔细的探究，所以需要分开学习和记录。
• 还有 detection of highly variable features以及scale data，也是方法和参数众多，我觉得有必要对其进行仔细的探究，所以需要分开学习和记录
• 概要图以及系列博文可以参见链接。

这里主要记录了 FindVariableFeatures的学习过程。

实例

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

rm(list=ls())

# 使用read10X读取output of the cellranger pipeline from 10X，包括barcodes，genes，matrix.mtx三个文件
pbmc.data <- Read10X(data.dir = "D:/djs/pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19")
# 使用 CreateSeuratObject函数构造seurat对象
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k",
                           min.cells = 3, min.features = 200,
                           names.delim = "-",names.field = 1)

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

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

上面完成了Seurat对象的构建与QC.

==============================================================================================

下面开始进行normalization，然后开始 FindVariableFeatures

• normalization

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- SCTransform(pbmc, vst.flavor = "v2", verbose = FALSE,variable.features.n = 2000)

str(pbmc)
# 此处我们直接把两种normalization方法都进运行了
# LogNormalize把数据放在了 assays$RNA@data
# SCTransform把数据放在了 assays$SCT@data，同时assays$SCT@counts与assays$RNA@counts也有区别
# SCTransform运行过程中已经做好了 FindVariableFeatures --> ScaleData ，分别放在了assays$SCT@scale.data /assays$SCT@var.features 

• FindVariableFeatures

# 因为上面把两种normalization方法都进运行了，有两个data，所以需要指定assay = "RNA" / assay = "SCT"
pbmc <- FindVariableFeatures(pbmc,assay = "RNA" ,selection.method = "vst", nfeatures = 2000)
# Identification of highly variable features (feature selection)
# 这里得到的feature会在下面的pca降维种默认使用

# Identify the 10 most highly variable genes
head(VariableFeatures(pbmc,assay = "RNA"), 10)
head(VariableFeatures(pbmc,assay = "SCT"), 10)

这里你可以看到前十的高变基因就有差异了。
看到下面一个数据就更值得思考了。

• seurat提供的可视化方法

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc,assay = "RNA")
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc,assay = "RNA")
plot2 <- VariableFeaturePlot(pbmc,assay = "SCT")
plot1 + plot2

然后variable 选择有几种方法，我试了一下不同方法在默认设置下的结果：

x <- intersect(VariableFeatures(pbmc_vst),VariableFeatures(pbmc_disp))
length(x)
length(VariableFeatures(pbmc_mvp))
x <- intersect(VariableFeatures(pbmc_vst),VariableFeatures(pbmc_mvp))
length(x)
x <- intersect(VariableFeatures(pbmc_disp),VariableFeatures(pbmc_mvp))
length(x)

后续很多分析是基于variable features进行的，不同的normalization 流程和 feature selectiion对后续分析的具体影响暂时我还没搞清楚。

但是你可以看到normalization 流程的不同会影响 feature selection ⇒ feature selection 的方法不同会导致highly variable features不同 ⇒ highly variable features不同你的RunPCA and cluster以及annotation 都会有些区别。

seurat官方比较推荐走SCTransform normalization流程，上一篇博客也简单的比较一下两种normalization方法，最后我个人也比较推荐使用 SCTransform normalization流程。SCTransform normalization流程会默认完成feature selection的选择，后续直接RunPCA --> FindNeighbors -->FindClusters一套流程直接走完。