10x Genomics PBMC(二):10x Genomics test1

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']]@datapbmc[['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
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 211,561评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,218评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,162评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,470评论 1 283
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,550评论 6 385
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,806评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,951评论 3 407
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,712评论 0 266
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,166评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,510评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,643评论 1 340
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,306评论 4 330
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,930评论 3 313
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,745评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,983评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,351评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,509评论 2 348