10x Genomics PBMC（二）：10x Genomic
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
