冷冻工厂

V1

2022/11/03阅读:39主题:科技蓝

单细胞分析之质控(四)

学习目标

知道如何导入和读取数据,并了解数据的质控,能够对数据进行质控和分析。

1. 质控准备

在基因表达定量后,需要将这些数据导入到 R 中,以生成用于执行 QC(质控)。下面将讨论定量数据的格式,以及如何将其导入 R,以便可以继续工作流程中的 QC 步骤。

2. 数据来源

在本教程中,将使用scRNA-seq 数据集,该数据集是 Kang 等人 2017[1] 年一项大规模研究的一部分。在本文中,作者提出了一种算法,该算法利用遗传变异 (eQTL) 来确定每个包含单个细胞的液滴 (singlet) 的遗传身份,并识别包含来自不同个体的两个细胞的液滴 (doublet)。

用于测试他们算法的数据取自八名狼疮患者的外周血单核细胞 (PBMC) 组成,分为对照和干扰素 β 治疗(刺激)条件。

  • Raw data

该数据集在 GEO (GSE96583) 上可下载,但是可用的计数矩阵缺少线粒体读数,因此从 SRA (SRP102802) 下载了 BAM 文件。这些 BAM 文件被转换回 FASTQ 文件,然后通过 Cell Ranger 运行以获得将使用的计数数据。

注意:此数据集的计数数据也可从 10X Genomics 获得,并在 Seurat[2] 教程中使用。

  • Metadata

除了原始数据,还需要收集有关数据的信息;这称为Metadata。常常有一种直接放手去做的冲动,但如果对这些数据的来源样本一无所知,这并不是一个好的习惯。

下面提供了数据集的一些相关Metadata

  1. 文库是使用 10X Genomics 第 2 版制备的

  2. 样本在 Illumina NextSeq 500 上进行测序

  3. 来自八名狼疮患者的 PBMC 样本被分成两个等分试样

    • 一份 PBMC 被 100 U/mL 重组 IFN-β 激活 6 小时。
    • 第二个等分样未处理。
    • 6 小时后,将每种条件的 8 个样品汇集到两个池中。
  4. 分别鉴定了 12,138 和 12,167 个细胞,用于对照和刺激的合并样本。

  5. 由于样本是 PBMC,预计包含免疫细胞,例如:

    • B细胞
    • T细胞
    • NK细胞
    • 单核细胞
    • 巨噬细胞
    • 巨核细胞(可能)

推荐在质控或分析前,对自己的样本有充分的了解,这对于后续的分析十分有帮助。

3. 数据准备

  • 环境准备:> 参考文末往期推荐
  • 数据下载:> 数据地址[3]

4. 项目结构

涉及大量数据的研究中,最重要的部分之一是如何管理它。倾向于优先分析,但数据管理的许多其他重要方面,往往在第一次看到新数据中被忽视。哈佛大学的生物医学数据管理[4] 很好的讲述了这一过程。

数据管理的一个重要方面是组织。对于处理和分析数据的每个实验,通过创建计划的存储空间(目录结构)来组织被认为是最佳实践。

  • 创建目录结构
single_cell_rnaseq/
├── data
├── results
└── figures

5. 数据处理

  • 新建Rscript
touch quality_control.R
  • 加载包
# 在前面创建的脚本中,用R打开
library(SingleCellExperiment)
library(Seurat)
library(tidyverse)
library(Matrix)
library(scales)
library(cowplot)
library(RCurl)
  • 加载scRNA-seq count数据

无论用于处理原始scRNA-seq 序列数据的技术或管道如何,定量后表达数据的输出通常是相同的。也就是说,对于每个单独的样本,将拥有以下三个文件:

  1. 具有细胞ID的文件,代表所有定量的细胞

  2. 具有基因ID的文件,代表所有定量的基因

  3. 每个细胞的每个基因的计数矩阵

以上数据存放在data/ctrl_raw_feature_bc_matrix文件夹内。

  1. barcodes.tsv

这是一个文本文件,其中包含该样本的所有细胞条形码。条形码按矩阵文件中显示的数据顺序列出

barcodes.tsv
barcodes.tsv
  1. features.tsv

这是一个包含定量基因标识符的文本文件。标识符的来源可能是 Ensembl、NCBI、UCSC,但大多数情况下这些是官方基因符号。这些基因的顺序对应于矩阵文件中的行顺序。

features.tsv
features.tsv
  1. matrix.mtx

这是一个包含计数值矩阵的文本文件。行与上面的基因 ID 相关联,列对应于细胞条形码。请注意,此矩阵中有许多零值。

matrix.mtx
matrix.mtx

将此数据加载到 R 中,需要将这三个数据整合为一个计数矩阵,并且考虑到减少计算的原因,此计数矩阵是一个稀疏矩阵。

  • 不同的读取数据方法:
  1. readMM(): 这个函数来自 Matrix 包,它将标准矩阵转换为稀疏矩阵。 features.tsv 文件和barcodes.tsv 必须先单独加载到R 中,然后才能将它们组合起来。
  2. Read10X(): 此函数来自 Seurat 包,将直接使用 Cell Ranger 输出目录作为输入。使用这种方法,不需要加载单个文件,而是该函数将加载并将它们组合成一个稀疏矩阵。本文将采取这个办法。

使用 Cell Ranger 处理 10X 数据后,将拥有一个 outs 目录。在此目录中,有下列文件:

  1. web_summary.html: 报告不同的 QC 指标,包括映射指标、过滤阈值、过滤后估计的细胞数,以及过滤后每个细胞的读数和基因数量的信息。
  2. BAM alignment files: 用于可视化映射读取和重新创建FASTQ文件的文件(如果需要)
  3. **filtered_feature_bc_matrix:**包含使用 Cell Ranger 过滤的数据构建计数矩阵所需的所有文件的文件夹
  4. raw_feature_bc_matrix: 包含使用原始未过滤数据构建计数矩阵所需的所有文件的文件夹

虽然 Cell Ranger 对表达计数执行过滤,但希望执行自己的 QC 和过滤。鉴于此,只对 Cell Ranger 输出中的 raw_feature_bc_matrix 文件夹感兴趣。

如果有一个样本,可以生成计数矩阵,然后创建一个 Seurat 对象:

关于Seurat[5]对象

# 如何读取单个样本的 10X 数据(输出为稀疏矩阵)
ctrl_counts <- Read10X(data.dir = "data/ctrl_raw_feature_bc_matrix")

# 将计数矩阵转为 Seurat 对象
ctrl <- CreateSeuratObject(counts = ctrl_counts, min.features = 100)
# min.features 参数指定每个细胞需要检测的最小基因数。

min.features 参数将过滤掉质量差的细胞,这些细胞可能只是封装了随机条形码而没有任何细胞存在。通常,检测到的基因少于 100 个的细胞不被考虑用于分析。

当使用 Read10X() 函数读入数据时,Seurat 会自动为每个单元格创建一些元数据。此信息存储在Seurat对象内的 meta.data 中。

# 探索元数据
head(ctrl@meta.data)
meta.data
meta.data

元数据的列:

  • orig.ident: 如果已知,这通常包含样本标识,但默认为“SeuratProject

  • nCount_RNA: 每个单元格的 UMI

  • nFeature_RNA: 每个细胞检测到的基因数量

  • 使用 for 循环读取多个样本

在实践中,可能有几个样本需要读取数据,如果一次只读取一个,可能会变得乏味且容易出错。因此,为了使数据导入R更有效,可以使用 for循环,它将为给定的每个输入迭代一系列命令,并为每个样本创建 seurat 对象。

# 仅测试,无法运行。
for (variable in input){
 command1
 command2
 command3
}


# 利用循环读取数据
for (file in c("ctrl_raw_feature_bc_matrix""stim_raw_feature_bc_matrix")){
        seurat_data <- Read10X(data.dir = paste0("data/", file))
        seurat_obj <- CreateSeuratObject(counts = seurat_data, 
                                         min.features = 100
                                         project = file)
        assign(file, seurat_obj)  # 将对象赋值给变量
}

接下来,将这些对象合并到一个单独的 Seurat 对象中。这将使两个样品组一起运行 QC 步骤变得更加容易,并能够轻松地比较所有样品的数据质量。

# 创建一个合并的 Seurat 对象
merged_seurat <- merge(x = ctrl_raw_feature_bc_matrix, 
                       y = stim_raw_feature_bc_matrix, 
                       add.cell.id = c("ctrl""stim"))

# 合并两个以上的样本
merged_seurat <- merge(x = ctrl_raw_feature_bc_matrix, 
                       y = c(stim1_raw_feature_bc_matrix, stim2_raw_feature_bc_matrix,
                       stim3_raw_feature_bc_matrix), 
                       add.cell.id = c("ctrl""stim1""stim2""stim3"))

# 查看对象的元数据
head(merged_seurat@meta.data)
tail(merged_seurat@meta.data)

因为相同的单元格 ID 可用于不同的样本,所以使用add.cell.id参数为每个单元格 ID 添加一个特定于样本的前缀。

参考资料

[1]

Kang: https://www.nature.com/articles/nbt.4042

[2]

Seurat: https://satijalab.org/seurat/v3.0/immune_alignment.html

[3]

数据地址: https://www.dropbox.com/s/we1gmyb9c8jej2u/single_cell_rnaseq.zip?dl=1

[4]

数据管理: https://datamanagement.hms.harvard.edu/

[5]

Seurat: https://github.com/satijalab/seurat/wiki/Seurat

分类:

其他

标签:

医学

作者介绍

冷冻工厂
V1

欢迎关注同名公众号“冷冻工厂”