单细胞系列教程：质控（四）

学习目标

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

1. 质控准备

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

2. 数据来源

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

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

• Raw data

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

注意：此数据集的计数数据也可从 10X Genomics 获得，并在 Seurat 教程中使用。

• 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. 数据准备

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

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

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

1. features.tsv

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

1. matrix.mtx

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

将此数据加载到 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对象

# 如何读取单个样本的 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)

元数据的列:

• 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 添加一个特定于样本的前缀。