scRNA-seq—读入数据详解
学习目标
了解如何导入单细胞rna-seq实验的数据。
质量控制
流程
在量化基因表达之后，我们需要将该数据导入R，以生成用于执行QC的矩阵。在本课中，我们将讨论盘点数据可以采用的格式，以及如何将其读入R，以便我们可以继续工作流程中的QC步骤。我们还将讨论我们将使用的数据集和相关的元数据
探索示例数据集
这次，我们将使用单细胞RNA-seq数据集，该数据集是Kang等人在2017年进行的一项较大研究的一部分。在本文中，作者提出了一种利用遗传变异(eQTL)的计算算法，以确定包含单个细胞(单胞体)的每个液滴的遗传同一性，并识别包含来自不同个体(双胞体)的两个细胞的液滴。
用于测试其算法的数据由八名狼疮患者的外周血单核细胞（PBMC）组成，分为对照和干扰素β治疗（刺激）条件。
图片来源：Kang等，2017
Raw data
该数据集在GEO（GSE96583），但是可用的计数矩阵缺少线粒体读数，因此我们从SRA（SRP102802）下载了BAM文件。这些BAM文件被转换回FASTQ文件，然后通过Cell Ranger运行以获得我们将要使用的计数矩阵。
注意： 此数据集的计数也可以从10X Genomics免费获得，并用作Seurat教程的的示例数据。
Metadata
除了原始数据之外，我们还需要收集有关数据的信息；这称为元数据。人们往往拿到数据就会直接来开始探索这些数据，但如果我们对这些数据的来源样本一无所知，那就没有太大的意义了。
下面提供了我们数据集的一些相关元数据：
- 使用10X Genomics V2版化学试剂盒制备文库
- 样品在Illumina NextSeq 500上测序
- 将来自八名狼疮患者的PBMC样本各分成两等份
- 用100 U/mL重组IFN-β激活一份PBMC，持续6小时。
- 第二等分试样未经处理
- 6小时后，将每个条件下的8个样本混合在两个最终池(刺激细胞和对照细胞)中。我们将使用这两个混合样本
- 对照和刺激混合样本分别鉴定了12138和12167个细胞(去除双峰后)。
- 由于样本是PBMC，我们预计有以下免疫细胞：
- B细胞
- T细胞
- NK细胞
- 单核细胞
- 巨噬细胞
- 可能有巨核细胞
建议您在执行QC之前，对希望在数据集中看到的细胞类型有一定的期望值。这告知您是否具有低复杂度的细胞类型(几个基因的大量转录本)或线粒体表达水平较高的细胞。这将使我们能够在分析工作流程中考虑这些生物因素。
上述细胞类型都不是低复杂性的，也不是线粒体含量高的。
设置R环境
涉及大量数据的研究中最重要的部分之一是如何最好地管理这些数据。我们倾向于确定分析的优先顺序，但在第一眼看到新数据的兴奋中，数据管理的许多其他重要方面经常被忽略。HMS数据管理工作组深入讨论了数据创建和分析之外需要考虑的一些问题。
数据管理的一个重要方面是组织。对于您进行和分析数据的每个实验，最佳实践是通过创建计划的存储空间(目录结构)来组织。我们将在单细胞分析中这样做。
打开RStudio并创建一个名为single_cell_rnaseq的新R项目。然后，创建以下目录：
single_cell_rnaseq/
├── data
├── results
└── figures
下载资料
将每个样本的输出文件夹从Cell Ranger下载到data文件夹
- Control sample
- Stimulated sample
解压下载好的文件夹，并在Rstudio中查看
新建脚本 新建R脚本，并注释（这个注释可忽略，无关紧要）
# February 2020
# HBC single-cell RNA-seq workshop
# Single-cell RNA-seq analysis - QC
将R脚本保存为quality_control.R。此时的工作目录如下：
工作目录
加载R包 没有安装的要提前安装。至于如何安装，可以看这个教程“【紧急通知】下载R包却联网失败？初学者的痛”
# Load libraries
library(SingleCellExperiment)
library(Seurat)
library(tidyverse)
library(Matrix)
library(scales)
library(cowplot)
library(RCurl)
导入scRNA-seq数据
无论使用哪种技术或流程来处理您的单细胞RNA-seq序列数据，输出通常都是相同的。也就是说，对于每个单独的样本，您将拥有以下三个文件：
- 包含细胞ID的文件，表示量化的所有细胞
- 包含基因ID的文件，表示量化的所有基因
- 每个细胞的每个基因的表达矩阵
可以通过单击data/ctrl_raw_feature_bc_matrix文件夹浏览这些文件：
- barcodes.tsv 这是一个文本文件，其中包含该样品的所有细胞条形码。条形码按照矩阵文件中显示的数据顺序列出(即这些是列名)。
cell_id
- features.tsv 这是一个文本文件，其中包含量化基因的标识符。根据您在量化方法中使用的参考(即Ensembl、NCBI、UCSC)的不同，标识符的来源可能会有所不同，但大多数情况下，这些都是官方的基因符号。这些基因的顺序对应于矩阵文件中行的顺序(即，这些是行名)。
gene_id
- matrix.mtx 这是一个包含计数值矩阵的文本文件。行与上面的基因ID相关联，列与细胞条形码相对应。需注意的是，此矩阵中有许多零值。
matrix
将这些数据加载到R中需要使用允许我们有效地将这三个文件组合成单个计数矩阵的函数。但是，我们将使用的函数不是创建常规矩阵数据结构，而是创建稀疏矩阵，以改进处理庞大计数矩阵所需的空间量、内存和CPU。
读取数据的不同方法：
readMM()：此函数来自Matrix包，它将把我们的标准矩阵转换为稀疏矩阵。首先必须先将features.tsv文件和barcodes.tsv分别加载到R中，然后再将它们合并。有关如何执行此操作的具体代码和说明，请参阅其他的材料。 2.Read10X()：此功能来自Seurat软件包，并将使用Cell Ranger输出目录作为输入。这样，不需要加载单个文件，而是该函数将加载并将它们合并为一个稀疏矩阵。我们将使用此功能加载数据！
读取一个样本（`Read10x()`）
当使用10X数据及其专有软件Cell Ranger时，您将始终拥有outs目录。在此目录中，您将发现许多不同的文件，包括：
web_summary.html：该报告探讨了不同的QC指标，包括映射指标，过滤阈值，过滤后估计的细胞数以及过滤后每个细胞的读取数和基因数的信息。BAM alignment files：用于可视化映射的读取和重新创建FASTQ文件的文件（如果需要）。filtered_feature_bc_matrix：该文件夹包含使用Cell Ranger 过滤数据构造计数矩阵所需的所有文件。raw_feature_bc_matrix：该文件夹包含使用原始数据(未经过滤)构造计数矩阵所需的所有文件。
我们主要对raw_feature_bc_matrix感兴趣，因为我们希望执行QC和过滤标准的同时考虑上我们自己的实验/生物学问题。
如果我们只有一个样本，我们可以生成计数矩阵，然后创建一个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 参数指定每个细胞需要检测的最小基因数量。此参数将过滤掉质量较差的细胞，这些细胞可能只是封装了随机barcodes，而没有任何真实的细胞。通常，检测到的基因少于100个的细胞不会被考虑进行分析。
当您使用Read10X()函数读入数据时，Seurat会自动为每个细胞创建一些元数据。此信息存储在seurat对象的meta.data槽中(更多内容请参阅下面的注释)。
Seurat对象是一个自定义的类列表对象，具有定义明确的空间来存储特定的信息/数据。您可以在此链接中找到有关Seurat对象插槽的更多信息。
# Explore the metadata
head(ctrl@meta.data)
元数据列是什么？
orig.ident：这通常包含样本标识(如果已知)，但默认为“SeuratProject”nCount_RNA：每个细胞的UMI号nFeature_RNA：每个细胞检测到的基因数量
读取多个样本`for loop`
在实践中，一般可能需要读取几个样本，同样使用我们前面讨论的两个函数(read10X()或readMM())中的一个来读入数据。为了更有效地将数据导入到R中，我们可以使用for循环，该循环将对给定的每个输入执行一系列命令。
在R中，语法如下：
## DO NOT RUN
for (variable in input){
command1
command2
command3
}
我们今天将使用的for循环将遍历两个样本“file”，并为每个样本执行两个命令
- (1)读入计数数据(
Read10X()) - (2)从读入数据创建Seurat对象(
CreateSeuratObject())：
# Create each individual Seurat object for every 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循环: Step 1: 指定输入
对于此数据集，我们有两个样本想要为以下对象创建Seurat对象：
- ctrl_raw_feature_bc_matrix
- stim_raw_feature_bc_matrix
我们可以使用c()在for循环的输入部分中将这些样本指定为向量的元素。我们将这些赋值给一个变量，我们可以随心所欲地给该变量命名(尽量给它起一个有意义的名称)。在本例中，我们将变量命名为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")){
step2：读入数据作为输入
通过对for loop添加一行以读取数据来继续操作Read10X()：
- 在这里，我们需要指定文件的路径，因此我们将data/使用
paste0()函数将目录添加到样本文件夹名称的前面。
## DO NOT RUN
seurat_data <- Read10X(data.dir = paste0("data/", file))
根据10x计数数据创建seurat对象
使用CreateSeuratObject()函数创建Seurat对象，并添加参数project，在其中可以添加样品名称。
## DO NOT RUN
seurat_obj <- CreateSeuratObject(counts = seurat_data,
min.features = 100,
project = file)
根据样本将Seurat对象分配给新变量
最后一个命令assign是将创建的Seurat对象（seurat_obj）添加到新变量。这样，当我们迭代并移动到输入中的下一个样本时，我们将不会覆盖在前一个迭代中创建的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"))
因为相同的细胞IDs可以用于不同的样本，所以我们使用add.cell.id参数为每个细胞ID添加一个特定于样本的前缀。如果我们查看合并对象的元数据，我们应该能够看到行名中的前缀：
# Check that the merged object has the appropriate sample-specific prefixes
head(merged_seurat@meta.data)
tail(merged_seurat@meta.data)
References
[1] Kang等人在2017年: https://www.nature.com/articles/nbt.4042
[2] GSE96583: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE96583
[3] SRP102802: https://www-ncbi-nlm-nih-gov.ezp-prod1.hul.harvard.edu/sra?term=SRP102802
[4] Cell Ranger: https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger
[5] Seurat教程的: https://satijalab.org/seurat/v3.0/immune_alignment.html
[6] HMS数据管理工作组: https://datamanagement.hms.harvard.edu/hms-data-management-working-group
[7] Control sample: https://www.dropbox.com/sh/73drh0ipmzfcrb3/AADMlKXCr5QGoaQN13-GbeKSa?dl=1
[8] Stimulated sample: https://www.dropbox.com/sh/cii4j356moc08w5/AAC2c3jfvh2hHWPmEaVsZKRva?dl=1
[9] 其他的材料: https://hbctraining.github.io/scRNA-seq/lessons/readMM_loadData.html
[10] 一个Seurat对象: https://github.com/satijalab/seurat/wiki/Seurat
[11] 在此链接中: https://github.com/satijalab/seurat/wiki/Seurat
注：以上内容来自哈佛大学生物信息中心（HBC）的教学团队的生物信息学培训课程。原文链接：https://hbctraining.github.io/scRNA-seq/schedule/