从读取数据开始。

Read10X ()函数 读取 cellranger pipeline输出的10X单细胞数据，返回一个(UMI)count矩阵。此矩阵中行（row）代表基因，列（col）代表细胞

raw data

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

# 加载 PBMC dataset
pbmc.data <- Read10X(data.dir = "pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19/")

加载完得到pbmc.data

接下来我们使用 count 矩阵来创建一个 Seurat 对象。

该对象充当一个容器，它包含单细胞数据(如count矩阵)和分析结果(如 PCA 或聚类)。
比如,count matrix储存在pbmc[["RNA"]]@counts.

#使用raw data(non-normalized data)创建Seurat 对象
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc

pbmc

Seurat 的 scRNA-seq 数据的标准预处理流程。

包括 QC 、数据标准化、缩放和高度可变基因的检测。

1. 质量控制（QC）及选择细胞进行下游分析

Seurat允许你简单的过滤一下你的细胞，质控的一般指标包括：

(1)在每个细胞里检测到的基因

-->低质量细胞或者空的droplets只有非常少量的基因
-->细胞doublets 或者 multiplets 有非常高的基因count数

(2)在一个细胞内检测到的分子总数

(3)读取到线粒体基因组的比例

-->低质量/死细胞通常有很高的线粒体污染
-->我们使用PercentageFeatureSet()函数计算线粒体QC指标 （计算来自一系列特征的计数百分比）
-->使用所有以MT-开头的基因作为线粒体基因

# The "[[ ]]"operator can add columns to object metadata. This is a great place to stash QC stats
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
head(pbmc@meta.data, 5)

image.png

我们还可以可视化QC指标，并用它们来过滤细胞：

（1）将unique基因count数超过2500，或者小于200的细胞过滤掉。
（2）把线粒体count数占5%以上的细胞过滤掉

# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

image.png

（画小提琴一直报错，暂未解决）

我们还可以用FeatureScatter函数来可视化特征-特征之间的关系，可以使用Seurat对象里的任何东西，如对象中的列、PC分数等。

plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2

2. 数据标准化

在去除掉不想要的细胞后，就可以标准化数据了。在默认情况下，我们使用global-scaling标准化方法，称为“LogNormalize”，这种方法是利用总的表达量对每个细胞里的基因表达值进行标准化，乘以一个scale factor（默认值是10000），再用log转换一下。标准化后的数据存放在pbmc[["RNA"]]@data里。

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

image.png

3. 鉴定高度变化的基因

接下来计算数据集中表现出细胞间高度差异的特征子集(即，它们在某些细胞中高表达，而在其他细胞中低表达)。在下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

#识别出10个差异最大的基因
top10 <- head(VariableFeatures(pbmc), 10)

# 带标签或不带标签作图展示差异基因
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

image.png