学习目标
知道如何导入和读取数据，并了解数据的质控，能够对数据进行质控和分析。
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：
-
文库是使用
10X Genomics第 2 版制备的 -
样本在
Illumina NextSeq 500上进行测序 -
来自八名狼疮患者的 PBMC 样本被分成两个等分试样
- 一份 PBMC 被 100 U/mL 重组 IFN-β 激活 6 小时。
- 第二个等分样未处理。
- 6 小时后，将每种条件的 8 个样品汇集到两个池中。
-
分别鉴定了 12,138 和 12,167 个细胞，用于对照和刺激的合并样本。
-
由于样本是 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 序列数据的技术或管道如何，定量后表达数据的输出通常是相同的。也就是说，对于每个单独的样本，将拥有以下三个文件：
-
具有细胞
ID的文件，代表所有定量的细胞 -
具有基因
ID的文件，代表所有定量的基因 -
每个细胞的每个基因的计数矩阵
以上数据存放在data/ctrl_raw_feature_bc_matrix文件夹内。
barcodes.tsv
这是一个文本文件，其中包含该样本的所有细胞条形码。条形码按矩阵文件中显示的数据顺序列出
barcodes.tsv
features.tsv
这是一个包含定量基因标识符的文本文件。标识符的来源可能是 Ensembl、NCBI、UCSC，但大多数情况下这些是官方基因符号。这些基因的顺序对应于矩阵文件中的行顺序。
features.tsv
matrix.mtx
这是一个包含计数值矩阵的文本文件。行与上面的基因 ID 相关联，列对应于细胞条形码。请注意，此矩阵中有许多零值。
matrix.mtx
将此数据加载到 R 中，需要将这三个数据整合为一个计数矩阵，并且考虑到减少计算的原因，此计数矩阵是一个稀疏矩阵。
- 不同的读取数据方法：
-
readMM(): 这个函数来自Matrix包，它将标准矩阵转换为稀疏矩阵。features.tsv文件和barcodes.tsv必须先单独加载到R中，然后才能将它们组合起来。 -
Read10X(): 此函数来自Seurat包，将直接使用Cell Ranger输出目录作为输入。使用这种方法，不需要加载单个文件，而是该函数将加载并将它们组合成一个稀疏矩阵。本文将采取这个办法。
使用 Cell Ranger 处理 10X数据后，将拥有一个 outs目录。在此目录中，有下列文件：
-
web_summary.html: 报告不同的QC指标，包括映射指标、过滤阈值、过滤后估计的细胞数，以及过滤后每个细胞的读数和基因数量的信息。 -
BAM alignment files: 用于可视化映射读取和重新创建
FASTQ文件的文件（如果需要） -
filtered_feature_bc_matrix:包含使用Cell Ranger过滤的数据构建计数矩阵所需的所有文件的文件夹 -
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)
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 添加一个特定于样本的前缀。
本文由mdnice多平台发布
网友评论