01_pbmc3k_sctf

clp

03 June, 2020

请注意，本教程借用自Seurat website，完成于2020年6月3日。

Seurat Object

在本教程中，我们将分析10X基因组公司(10X Genomics)免费提供的外周血单核细胞(PBMC)数据集。在Illumina NextSeq 500平台上对2700个单细胞进行了测序。

任务：在工作目录找到zip文件filtered_gene_bc_matrices.zip并解压(“解压到这里”)。例如，解压缩文件的位置应该如下所示， 原始数据下载地址：https://satijalab.org/seurat/vignettes.html

image.png

我们从读入数据开始。Read10X函数从10X读入cellranger pipeline的输出文件，返回唯一的分子标识符计数矩阵((UMI) count matrix)。该矩阵中的值表示在每个细胞(列)中检测到的每个特征值(即基因；行)的分子数量。

接下来，我们使用count matrix创建一个Seurat对象。该对象充当一个容器，包含单个细胞数据集的数据(如计数矩阵)和分析(如PCA或聚类结果)。有关Seurat对象结构的技术讨论，请查看我们的GitHub Wiki。例如，计数矩阵存储在pbmc[["RNA"]]@counts中。

首先加载必要的R包,

library(dplyr)
library(Seurat)
library(ggplot2)

读取我们上面下载的Cellranger Count输出文件，并创建一个Seurat对象，

# 该目录由10x Genomics cellranger count生成。目录应该包含3个文件：`barcodes.tsv`, `genes.tsv`和`matrix.mtx`。输出将是sparse矩阵(`dgCMatrix`)。
pbmc.data <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19")
# 使用原始(非标准化数据)初始化Seurat对象。
# 考虑出现在至少3个细胞中的基因，以及表达至少100个独特基因的细胞
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 100)
#> Warning: Feature names cannot have underscores ('_'), replacing with dashes
#> ('-')
dim(pbmc)
#> [1] 13714  2700

筛选指标1：考虑保留出现在至少3个细胞中的基因，以及表达至少100个独特基因的细胞

任务：探索PBMC刚刚创建的Seurat对象，该对象看起来是什么样子？

class(pbmc)
#> [1] "Seurat"
#> attr(,"package")
#> [1] "Seurat"
getSlots("Seurat")
#>            assays         meta.data      active.assay      active.ident 
#>            "list"      "data.frame"       "character"          "factor" 
#>            graphs         neighbors        reductions      project.name 
#>            "list"            "list"            "list"       "character" 
#>              misc           version          commands             tools 
#>            "list" "package_version"            "list"            "list"
str(pbmc)
#> Formal class 'Seurat' [package "Seurat"] with 12 slots
#>   ..@ assays      :List of 1
#>   .. ..$ RNA:Formal class 'Assay' [package "Seurat"] with 8 slots
#>   .. .. .. ..@ counts       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#>   .. .. .. .. .. ..@ i       : int [1:2282976] 29 73 80 148 163 184 186 227 229 230 ...
#>   .. .. .. .. .. ..@ p       : int [1:2701] 0 779 2131 3260 4220 4741 5522 6304 7094 7626 ...
#>   .. .. .. .. .. ..@ Dim     : int [1:2] 13714 2700
#>   .. .. .. .. .. ..@ Dimnames:List of 2
#>   .. .. .. .. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ...
#>   .. .. .. .. .. .. ..$ : chr [1:2700] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
#>   .. .. .. .. .. ..@ x       : num [1:2282976] 1 1 2 1 1 1 1 41 1 1 ...
#>   .. .. .. .. .. ..@ factors : list()
#>   .. .. .. ..@ data         :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#>   .. .. .. .. .. ..@ i       : int [1:2282976] 29 73 80 148 163 184 186 227 229 230 ...
#>   .. .. .. .. .. ..@ p       : int [1:2701] 0 779 2131 3260 4220 4741 5522 6304 7094 7626 ...
#>   .. .. .. .. .. ..@ Dim     : int [1:2] 13714 2700
#>   .. .. .. .. .. ..@ Dimnames:List of 2
#>   .. .. .. .. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ...
#>   .. .. .. .. .. .. ..$ : chr [1:2700] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
#>   .. .. .. .. .. ..@ x       : num [1:2282976] 1 1 2 1 1 1 1 41 1 1 ...
#>   .. .. .. .. .. ..@ factors : list()
#>   .. .. .. ..@ scale.data   : num[0 , 0 ] 
#>   .. .. .. ..@ key          : chr "rna_"
#>   .. .. .. ..@ assay.orig   : NULL
#>   .. .. .. ..@ var.features : logi(0) 
#>   .. .. .. ..@ meta.features:'data.frame':   13714 obs. of  0 variables
#>   .. .. .. ..@ misc         : NULL
#>   ..@ meta.data   :'data.frame': 2700 obs. of  3 variables:
#>   .. ..$ orig.ident  : Factor w/ 1 level "pbmc3k": 1 1 1 1 1 1 1 1 1 1 ...
#>   .. ..$ nCount_RNA  : num [1:2700] 2419 4903 3147 2639 980 ...
#>   .. ..$ nFeature_RNA: int [1:2700] 779 1352 1129 960 521 781 782 790 532 550 ...
#>   ..@ active.assay: chr "RNA"
#>   ..@ active.ident: Factor w/ 1 level "pbmc3k": 1 1 1 1 1 1 1 1 1 1 ...
#>   .. ..- attr(*, "names")= chr [1:2700] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
#>   ..@ graphs      : list()
#>   ..@ neighbors   : list()
#>   ..@ reductions  : list()
#>   ..@ project.name: chr "pbmc3k"
#>   ..@ misc        : list()
#>   ..@ version     :Classes 'package_version', 'numeric_version'  hidden list of 1
#>   .. ..$ : int [1:3] 3 1 5
#>   ..@ commands    : list()
#>   ..@ tools       : list()

标准预处理流程

以下步骤包含Seurat中scRNA-seq数据的标准预处理流程。它们代表了基于QC度量的单元格选择和过滤、数据标准化和缩放，以及高度可变特征的检测。

QC质控和筛选进一步分析的细胞数据

使用Seurat，您可以轻松浏览QC指标，并根据任何用户定义的标准过滤细胞。几个QC指标commonly used包括：

• 在每个细胞中检测到的唯一基因的数量。
  • 低质量的细胞或空滴通常只有很少的基因。
  • 细胞二倍体或多倍体可能表现出异常高的基因计数。
• 同样，细胞内检测到的分子总数(与唯一基因密切相关)。
• 映射到线粒体基因组的读数百分比。
  • 低质量/濒临死亡的细胞通常表现出广泛的线粒体污染。
  • 我们使用PercentageFeatureSet函数计算线粒体QC度量，该函数计算源自一组特征的计数的百分比。
  • 我们使用以MT-开头的所有基因的集合作为一组线粒体基因。

# `[[`运算符可以向对象元数据添加列。这是存放QC统计数据的好地方
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

从特征计数矩阵中，对象收集一些有用的度量并将它们存储到slot(槽)meta.data中。

head(pbmc@meta.data)
#>                  orig.ident nCount_RNA nFeature_RNA percent.mt
#> AAACATACAACCAC-1     pbmc3k       2419          779  3.0177759
#> AAACATTGAGCTAC-1     pbmc3k       4903         1352  3.7935958
#> AAACATTGATCAGC-1     pbmc3k       3147         1129  0.8897363
#> AAACCGTGCTTCCG-1     pbmc3k       2639          960  1.7430845
#> AAACCGTGTATGCG-1     pbmc3k        980          521  1.2244898
#> AAACGCACTGGTAC-1     pbmc3k       2163          781  1.6643551

在下面的示例中，我们将QC指标可视化，并使用这些指标来过滤细胞。绘制特征数、读取计数和线粒体百分比的分布。

# 在QC之前将QC指标可视化为小提琴曲线图
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size=0)

image.png

通过检查小提琴图，我们进行质量控制。我们过滤具有独特特征数超过2500或小于200的单元格。我们过滤线粒体数量大于5%的细胞。

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
# 将QC指标可视化为QC后的小提琴曲线图
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size=0)

image.png

# FeatureScatter通常用于可视化要素-要素关系，但也可以使用对象计算的任何内容，即对象元数据中的列、PC分数等。

plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt") + theme(legend.position = 'none')
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") + theme(legend.position = 'none')
CombinePlots(plots = list(plot1, plot2))
#> Warning: CombinePlots is being deprecated. Plots should now be combined
#> using the patchwork system.

image.png

归一化和缩放数据

由于技术变异或生物变异，每个细胞都有不同的测序深度。测序深度与UMI计数高度相关，但它也随每个基因的不同而不同。通常，我们使用全局缩放归一化方法LogNormalize，该方法用总表达式归一化每个细胞的特征表达式测量值，将其乘以比例因子(默认情况下为10,000)，然后对结果进行对数转换。

识别高度可变的特征(特征选择)

基因(特征)的数量超过2万个。在归一化之后，细胞间的许多基因丰度值几乎是相似的。我们想要选择细胞间高度可变的基因子集。通常，归一化读取计数是z变换的，并选择具有top-K最高方差的基因。

SCTransform: 三步并为一步

上述三个步骤由SCTransform执行。SCTransform使用正则化负二项回归计算scRNA-seq数据中的技术噪声模型。该模型的残差是regularized negative binomial regression规格化值，可以是正值，也可以是负值。给定细胞中给定基因的正残差表明，考虑到该基因在群体中的平均表达和细胞测序深度，我们观察到的UMI比预期的要多，而负残差则表明相反。

# run sctransform
pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = FALSE)

sctransfrom的结果存储在SCT分析中。检测有3个槽：``pbmc[['SCT']]@counts,pbmc[['SCT']]@data和pbmc[['SCT']]@scale.data`。

• pbmc[['SCT']]@scale.data 包含残差(归一化值)，直接作为<ins style="box-sizing: border-box;">PCA(即降维)</ins>的输入。请注意，此矩阵是非稀疏的，因此，如果存储所有基因，可能会占用大量内存。为了节省内存，通过在SCTransform函数调用中默认设置return.only.var.genes = TRUE，我们只存储变量基因的这些值。

• 校正的UMI计数存储在pbmc[['SCT']]@counts中，它由Pearson残差调整。此数据通常用于<ins style="box-sizing: border-box;">差异基因表达检测</ins>。

• 我们将这些校正后的计数的对数归一化版本存储在pbmc[['SCT']]@data中，这对于<ins style="box-sizing: border-box;">可视化差异表达式</ins>非常有帮助。

此时，您可以保存对象，以便可以轻松地将其重新载入，而不必重新运行上面执行的计算密集型步骤，或者可以轻松地与协作者共享。

save(pbmc,file='filtered_gene_bc_matrices/hg19/01_pbmc3k_sctf.rd',compress=TRUE)

需要掌握的要点

• Seurat object structure
• QC
• Normalization
• Save R environment variables