从读取数据开始。
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
网友评论