单细胞 RNA-seq：质量控制设置
基因表达量化后，我们需要将此数据读入 R 中，执行 QC 的指标。我们将讨论可以预期的数据格式，以及如何将其读入 R 以便我们可以继续流程中的 QC 步骤。我们还将讨论我们将使用的数据集和相关的metadata
。
探索示例数据集
我们将使用单细胞 RNA-seq 数据集，该数据集是Kang 等人于2017年进行的一项较大研究的一部分。在本文中，作者提出了一种计算算法，该算法利用遗传变异 (eQTL) 来确定每个包含单个细胞 (singlet) 的液滴的遗传身份，并识别包含来自不同个体的两个细胞 (doublet) 的液滴。
用于测试其算法的数据由来自八名狼疮患者的外周血单核细胞 (PBMC) 组成，分为对照和干扰素 β 治疗（刺激）条件。
图片来源：Kang 等人，2017 年(https://www.nature.com/articles/nbt.4042)
原始数据
该数据集在 GEO ( GSE96583 ) 上可用，但是可用的计数矩阵缺少线粒体读数，因此我们从 SRA ( SRP102802 )下载了 BAM 文件。这些 BAM 文件被转换回 FASTQ 文件，然后通过Cell Ranger运行以获得计数数据。
注意： 此数据集的计数数据也可从 10X Genomics 免费获得，并在Seurat 教程中使用。
metadata(元数据)
除了原始数据，我们还需要收集有关数据的信息；这称为元数据。人们往往拿到数据就会直接来开始探索这些数据，但如果我们对这些数据的来源样本一无所知，那就没有太大的意义了。
下面提供了我们数据集的一些相关元数据：
使用 10X Genomics V2 版本化学试剂盒制备文库
样品在 Illumina NextSeq 500 上测序
-
将来自 8 名狼疮患者的 PBMC 样品各分成两等份。
- 用100 U/mL 重组 IFN-β 激活一份PBMC，持续 6 小时。
- 第二等份不处理。
- 6 小时后，将每种条件下的 8 个样品混合在两个最终池（受刺激细胞和对照细胞）中。我们将使用这两个合并样本。（我们没有对样本进行多路分解，因为论文中使用了 SNP 基因型信息进行多路分解，而这些数据的条形码/样本 ID 不容易获得。通常，您会对每个单独的样本进行多路分解和 QC，而不是合并样本。 )
用于对照和刺激的合并样品分别鉴定了 12138 和 12167 个细胞（去除双联体后）。
-
由于是 PBMC样本，预计有以下免疫细胞，例如：
- B细胞
- T细胞
- NK细胞
- 单核细胞
- 巨噬细胞
- 可能有巨核细胞
建议您在执行 QC 之前对您希望在数据集中看到的细胞类型有一些期望。这会告诉您是否有任何低复杂度的细胞类型（来自少数基因的大量转录本）或线粒体表达水平较高的细胞。我们将在分析工作流程中考虑这些生物因素。
预计上述细胞类型均不具有低复杂性或预计具有高线粒体含量。
设置R运行环境
我们将在 RStudio 项目中工作。为了继续学习，您应该已经下载了 R 项目。
如果您还没有这样做，可以使用此链接访问该项目。
下载后，您应该会在您的计算机上看到一个名为single_cell_rnaseq.zip的文件（可能在您的Downloads文件夹中）。
- 解压此文件。这将产生一个同名的文件夹。
- 将文件夹移动到计算机分析的位置。
- 打开文件夹。内容将类似于下面的屏幕截图。
-
找到
.Rproj file并双击它。这将打开带有“single_cell_rnaseq”项目加载的 RStudio。
项目管理
研究中的大量数据最重要的部分之一是如何最好地管理它。我们倾向于优先考虑分析，但数据管理的许多其他重要方面往往在第一次看到新数据中被忽视。该HMS数据管理工作组，深入讨论一些超出数据创建和分析。
数据管理的另一个重点是组织。对于您处理和分析每个实验的数据，通过创建计划的存储空间（目录结构）来进行组织被认为是最佳实践。我们将为我们的单细胞分析做到这一点。
查看您的项目空间，您会发现已经为您设置了一个目录结构：
single_cell_rnaseq/
├── data
├── results
└── figures
WINDOWS OS 用户注意事项： 当您解压后打开项目文件夹时，请检查您是否有一个
data文件夹，其子文件夹也称为data. 如果是这种情况，请将子文件夹中的所有文件移动到父data文件夹中。
创建新的脚本
接下来，打开一个新的 Rscript 文件，开始写一些注释，以指示该文件将包含的内容：
# July/August 2021
# HBC single-cell RNA-seq workshop
# Single-cell RNA-seq analysis - QC
将 Rscript 保存为quality_control.R. 您的工作目录应如下所示：
库加载
现在，我们可以加载必要的库：
# Load libraries
library(SingleCellExperiment)
library(Seurat)
library(tidyverse)
library(Matrix)
library(scales)
library(cowplot)
library(RCurl)
加载单细胞 RNA-seq 计数数据
无论用于处理原始单细胞 RNA-seq 序列数据的技术或pipeline如何，输出的表达的矩阵通常是相同的。也就是说，对于每个单独的样本，都拥有以下三个文件：
- 带有cell ID的文件，表示所有量化的细胞
- 带有gene ID的文件，代表所有量化的基因
- 一个计数的矩阵每个基因对应每一个细胞
我们可以通过单击data/ctrl_raw_feature_bc_matrix文件夹来浏览这些文件：
1. barcodes.tsv
这是一个文本文件，其中包含该样本的所有细胞barcode。条形码按照矩阵文件中数据的顺序列出（即这些是列名称）。
2. features.tsv
这是一个文本文件，其中包含量化基因的标识符。标识符可能因您在量化方法中使用的参考的来源（即 Ensembl、NCBI、UCSC）而异，但大多数情况下，这些是官方基因符号。这些基因的顺序对应于矩阵文件中行的顺序（即这些是行名称）。
3. matrix.mtx
这是一个包含计数值矩阵的文本文件。行与上面的基因 ID 相关联，列对应于细胞条形码。请注意，此矩阵中有许多零值。
将此数据加载到 R 中需要我们使用将这三个文件组合成单个计数矩阵的函数。然而，我们将使用的函数不是创建常规矩阵数据结构，而是创建一个稀疏矩阵，以减少使用我们庞大的计数矩阵所需的内存 (RAM)、处理容量 (CPU) 和存储量。
不同的读入数据方法包括：
-
readMM()：这个函数来自Matrix包，它将把我们的标准矩阵转换成一个稀疏矩阵。features.tsv和barcodes.tsv文件必须先被单独地加载到R，然后再将它们组合。有关如何执行此操作的特定代码和说明，请参阅这些附加材料。 -
Read10X()：此函数来自Seurat包，将直接使用 Cell Ranger 输出目录作为输入。使用这种方法不需要加载单个文件，而是该函数将加载它们并将它们组合成一个稀疏矩阵。我们将使用这个函数来加载我们的数据！
读取单个样本
使用 Cell Ranger 软件处理 10X 数据后，您将拥有一个outs目录（始终）。在此目录中，您将找到许多不同的文件，包括下面列出的文件：
-
web_summary.html：不同 QC 指标的报告，包括映射指标、过滤阈值、过滤后的估计细胞数以及有关过滤后每个细胞的read数和基因数的信息。 - BAM 对齐文件：用于可视化映射read和重新创建 FASTQ 文件的文件（如果需要）
-
filtered_feature_bc_matrix:包含使用 Cell Ranger 过滤的数据构建计数矩阵所需的所有文件的文件夹 -
raw_feature_bc_matrix：包含使用原始未过滤数据构建计数矩阵所需的所有文件的文件夹
虽然 Cell Ranger 对表达计数执行过滤（请参阅下面的注释），但我们希望执行我们自己的 QC 和过滤，因为我们需要考虑自己的实验/生物学问题。鉴于此，我们只对raw_feature_bc_matrixCell Ranger 输出中的****文件夹****感兴趣。
注意：为什么我们不使用
filtered_feature_bc_matrix文件夹？该filtered_feature_bc_matrix应用通过细胞内部过滤条件，我们没有什么细胞的控制，以保持或放弃。Cell Ranger 在生成
filtered_feature_bc_matrix时执行的过滤通常很好；然而，有时数据的质量可能非常高，而 Seurat 过滤过程可以去除高质量的细胞。此外，通常最好探索您自己的数据，同时考虑在过滤过程中应用阈值的实验的生物学特性。例如，如果您希望数据集中的特定细胞类型更小和/或转录活性不如数据集中的其他细胞类型，则这些细胞有可能被过滤掉。然而，在 Cell Ranger v3 中，将考虑不同大小的细胞（例如，肿瘤与浸润淋巴细胞），现在可能无法根据需要过滤尽可能多的低质量细胞。
如果我们只有一个样本，我们可以生成计数矩阵，然后创建一个 Seurat 对象：
Seurat 对象是一个自定义的类似列表的对象，它具有明确定义的空间来存储特定的信息/数据。您可以在此链接中找到有关 Seurat 对象中的更多信息。
# How to read in 10X data for a single sample (output is a sparse matrix)
ctrl_counts <- Read10X(data.dir = "data/ctrl_raw_feature_bc_matrix")
# Turn count matrix into a Seurat object (output is a Seurat object)
ctrl <- CreateSeuratObject(counts = ctrl_counts,
min.features = 100)
注意：该
min.features参数指定每个细胞需要检测的最小基因数。这将过滤掉质量较差的细胞，这些细胞可能只封装了随机barcode，而没有任何细胞存在。通常，检测到的基因少于 100 个的细胞不用于分析。
当您使用该Read10X()函数读入数据时，Seurat 会自动为每个细胞创建一些元数据。此信息存储在Seurat 对象内的meta.data中。
# Explore the metadata
head(ctrl@meta.data)
元数据的列是什么意思？
-
orig.ident：如果已知，这通常包含样本身份，但将默认为“SeuratProject” -
nCount_RNA: 每个细胞的 UMI 数量 -
nFeature_RNA: 每个细胞检测到的基因数量
for loop读取多个样本
在现实中，您可能有多个样本需要读取数据，如果一次读取一个样本，可能会变得乏味且容易出错。因此，为了使数据导入 R 更有效，我们可以使用for循环，该循环将为每个输入给定一系列命令，并为我们的每个样本创建 seurat 对象。
在 R 中，for循环具有以下结构/语法：
## DO NOT RUN
for (variable in input){
command1
command2
command3
}
今天我们将使用它来遍历两个示例文件夹并为每个示例执行两个命令，就像我们上面对单个示例所做的一样 。
(1) 读取计数数据 ( Read10X())
(2) 从读取的数据中创建 Seurat 对象数据 ( CreateSeuratObject()):
# Create a Seurat object for each sample
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)
}
让我们分解
for loop并查看不同的代码行：步骤 1：指定输入
对于这个数据集，我们有两个样本和两个关联的文件夹，我们将创建两个 Seurat 对象的输入：
ctrl_raw_feature_bc_matrixstim_raw_feature_bc_matrix我们可以在输入部分指定这些示例文件夹
for loop作为向量的元素使用c(). 我们将这些分配给一个变量，我们可以随意调用该变量（尝试给它一个有意义的名称）。在这个例子中，我们调用了变量file。在上述循环的执行过程中，
file将首先包含值“ctrl_raw_feature_bc_matrix”，通过命令一直运行到assign()。接下来，它将包含值“stim_raw_feature_bc_matrix”并再次运行所有命令。如果您有 15 个文件夹作为输入，而不是 2 个，则上述代码将针对您的每个数据文件夹运行 15 次。## DO NOT RUN # Create each individual Seurat object for (file in c("ctrl_raw_feature_bc_matrix", "stim_raw_feature_bc_matrix")){步骤 2：输入读入数据
我们可以
for loop通过添加一行来继续我们的读取数据Read10X()：
- 在这里，我们需要指定文件的路径，因此我们将使用
paste0()函数将data/目录添加到我们的示例文件夹名称中。## DO NOT RUN seurat_data <- Read10X(data.dir = paste0("data/", file))步骤 3：从 10X 计数数据创建 Seurat 对象
现在，我们可以使用该
CreateSeuratObject()函数创建 Seurat 对象，添加参数project，我们可以在其中添加样本名称。## DO NOT RUN seurat_obj <- CreateSeuratObject(counts = seurat_data, min.features = 100, project = file)第 4 步：将 Seurat 对象分配给基于样本的新变量
最后一个命令
assign是将 Seurat 对象创建 (seurat_obj) 到一个新变量。这样，当我们迭代并移动到我们的下一个样本时，我们input不会覆盖在前一次迭代中创建的 Seurat 对象：## DO NOT RUN assign(file, seurat_obj) }
现在我们已经创建了这两个对象，让我们快速浏览一下元数据：
# Check the metadata in the new Seurat objects
head(ctrl_raw_feature_bc_matrix@meta.data)
head(stim_raw_feature_bc_matrix@meta.data)
接下来，我们需要将这些对象合并为一个 Seurat 对象。同时运行两个样品组的 QC 步骤，我们能够轻松比较所有样品的数据质量。
我们可以使用Seurat 包中的merge()函数来做到这一点：
# Create a merged Seurat object
merged_seurat <- merge(x = ctrl_raw_feature_bc_matrix,
y = stim_raw_feature_bc_matrix,
add.cell.id = c("ctrl", "stim"))
由于相同的细胞 ID 可用于不同的样本，因此我们使用参数为每个细胞 ID添加特定的样本的前缀add.cell.id。如果我们查看合并对象的元数据，我们应该能够看到行名中的前缀：
# Check that the merged object has the appropriate sample-specific prefixes
head(merged_seurat@meta.data)
tail(merged_seurat@meta.data)
