scRNA-seq/03_SC_quality_control-setup.md at master · hbctraining/scRNA-seq · GitHub
https://github.com/hbctraining/scRNA-seq/blob/master/lessons/03_SC_quality_control-setup.md
在第三节作者展示了数据载入和合并，作者代码如下，我略微做了更改。
学习目标：
- 了解如何从单细胞RNA序列实验中引入数据
单细胞RNA序列：质量控制设置
定量基因表达后，我们需要将此数据带入R中以生成用于执行QC的指标。在本课程中，我们将讨论可以预期的计数数据格式，以及如何将其读入R，以便我们继续进行工作流程中的QC步骤。我们还将讨论将要使用的数据集和相关的元数据。
探索示例数据集
我们将使用单细胞RNA-seq数据集，该数据集是Kang等人在2017年进行的一项较大研究的一部分。在本文中，作者提出了一种计算算法，该算法利用遗传变异（eQTL）来确定每个包含单个细胞的液滴的遗传同一性（单一），并鉴定包含来自不同个体的两个细胞的液滴（双重）。
用于测试其算法的数据由八名狼疮患者的合并外周血单个核细胞（PBMC）组成，分为对照和干扰素β治疗（刺激）条件。
图片来源：Kang等，2017
原始数据
该数据集可在GEO（GSE96583）上获得，但是可用的计数矩阵缺少线粒体读数，因此我们从SRA（SRP102802）下载了BAM文件。这些BAM文件被转换回FASTQ文件，然后通过Cell Ranger运行以获得我们将要使用的计数数据。
注意：此数据集的计数也可以从10X Genomics免费获得，并用作Seurat教程的一部分。
元数据
除了原始数据外，我们还需要收集有关数据的信息。这就是所谓的元数据。只是探索数据，但是如果我们不知道该数据所源自的样本，那没有意义。
下面提供了与我们的数据集相关的一些元数据：
使用10X Genomics第2版化学试剂盒制备文库
样品在Illumina NextSeq 500上测序
将来自八名狼疮患者的PBMC样本分别分成两等份。
用100 U / mL重组IFN-β激活PBMC一份，持续6小时。
第二等分试样未经处理。
6小时后，将每种条件的8个样品一起收集到两个最终的收集池中（受激细胞和对照细胞）。我们将处理这两个合并的样本。
分别鉴定了12138和12167个细胞（去除双峰后）作为对照样品和刺激后的合并样品。
由于样本是PBMC，因此我们期望有免疫细胞，例如：
B细胞
T细胞
NK细胞
单核细胞
巨噬细胞
可能是巨核细胞
建议您对执行QC之前希望在数据集中看到的像元类型有一些期望。这将告知您是否具有低复杂度的细胞类型（许多基因的转录本）或线粒体表达水平较高的细胞。这将使我们能够在分析工作流程中考虑这些生物学因素。
预期上述细胞类型都不具有低复杂性或预期具有高线粒体含量。
设置R环境
研究涉及大量数据的最重要部分之一就是如何最好地对其进行管理。我们倾向于对分析进行优先级排序，但是数据管理中还有许多其他重要方面经常被人们忽略，因此对新数据的初学者一见倾心。该HMS数据管理工作组，讨论深入一些事情要考虑超出数据创建和分析。
数据管理的一个重要方面是组织。对于您进行和分析数据的每个实验，通过创建计划的存储空间（目录结构）来组织起来被认为是最佳实践。我们将对单细胞分析进行此操作。
打开RStudio并创建一个名为的新R项目single_cell_rnaseq。然后，创建以下目录：
single_cell_rnaseq/
├── data
├── results
└── figures
下载资料
右键单击下面的链接，以将每个样本的输出文件夹从Cell Ranger下载到data文件夹中：
现在，让我们通过双击解压缩刚刚下载的两个压缩文件夹，以便我们可以在RStudio中查看它们的内容。
新项目
接下来，打开一个新的Rscript文件，并从一些注释开始以指示该文件将包含的内容：
＃ 2020年2月
＃ HBC单细胞RNA-seq的流程
＃单细胞RNA序列分析-QC
将Rscript另存为quality_control.R。您的工作目录应如下所示：
加载库
现在，我们可以加载必要的R包：
# Load libraries
library(SingleCellExperiment)
library(Seurat)
library(tidyverse)
library(Matrix)
library(scales)
library(cowplot)
library(RCurl)
加载单细胞RNA-seq计数数据
无论处理单细胞RNA-seq序列数据的技术或流程如何，其输出通常都是相同的。也就是说，对于每个单独的样本，您将拥有以下三个文件：
具有细胞ID的文件，代表所有量化的细胞
一个带有基因ID的文件，代表所有量化的基因
一个计数的矩阵每个基因的每一个细胞
我们可以通过单击data/ctrl_raw_feature_bc_matrix文件夹浏览这些文件：
1 barcodes.tsv
这是一个文本文件，其中包含该样品存在的所有细胞条形码。条形码以矩阵文件中显示的数据顺序列出（即，这些是列名）。
2 features.tsv
这是一个文本文件，其中包含定量基因的标识符。标识符的来源可能会有所不同，具体取决于您在定量方法中使用的参考文献（即Ensembl，NCBI，UCSC），但最常见的是官方基因符号。这些基因的顺序与矩阵文件中行的顺序相对应（即，这些是行名）。
3. matrix.mtx
这是一个文本文件，其中包含计数值的矩阵。行与上面的基因ID相关联，列与细胞条形码相对应。请注意，此矩阵中有许多零值。
将此数据加载到R中要求我们使用函数，这些函数使我们能够将这三个文件有效地组合为一个计数矩阵。但是，我们将使用创建函数来创建稀疏矩阵，而不是创建规则的矩阵数据结构，以改善处理庞大计数矩阵所需的空间，内存和CPU。
读取数据的不同方法包括：
readMM()：此函数来自Matrix包，它将把我们的标准矩阵变成稀疏矩阵。该features.tsv文件barcodes.tsv必须首先单独加载到R中，然后再合并。有关如何执行此操作的特定代码和说明，请参阅我们的其他材料。Read10X()：此功能来自Seurat软件包，并将使用Cell Ranger的输出目录作为输入。这样，不需要加载单个文件，而是该函数将加载并将它们合并为一个稀疏矩阵。我们将使用此功能加载数据！
读取一个样本（read10X()）
当使用10X数据及其专有软件Cell Ranger时，您将始终有一个outs目录。在此目录中，您可以找到许多不同的文件，包括：
web_summary.html：该报告探讨了不同的QC指标，包括映射指标，过滤阈值，过滤后的估计细胞数以及过滤后每个细胞的读取数和基因数的信息。BAM对齐文件：用于可视化映射的读取和重新创建FASTQ文件的文件（如果需要）
filtered_feature_bc_matrix：包含使用Cell Ranger过滤的数据构造计数矩阵所需的所有文件的文件夹raw_feature_bc_matrix：包含使用原始未过滤数据构建计数矩阵所需的所有文件的文件夹
我们主要对感兴趣，raw_feature_bc_matrix因为我们希望执行自己的质量控制和过滤，同时考虑到我们的实验/生物学系统的生物学特性。
如果我们只有一个样本，我们可以生成计数矩阵，然后创建一个Seurat对象：
ctrl_counts <- Read10X(data.dir = "data/ctrl_raw_feature_bc_matrix")
＃将计数矩阵转换为Seurat对象（输出是Seurat对象）
ctrl <- CreateSeuratObject(counts = ctrl_counts,min.features = 100)
注意：该
min.features参数指定每个细胞需要检测的最小基因数量。此参数将过滤掉质量较差的单元，这些单元可能只是封装了随机条形码而没有任何单元。通常，检测不到少于100个基因的细胞不考虑进行分析。
当您使用该Read10X()功能读取数据时，Seurat会为每个细胞自动创建一些元数据。此信息存储在Seurat对象的meta.data插槽中（请参阅下面的注释中的更多信息）。
Seurat对象是一个自定义的类似列表的对象，具有定义明确的空间来存储特定的信息/数据。您可以在此链接中找到有关Seurat对象中插槽的更多信息。
＃探索元数据
head(ctrl@meta.data)
元数据列是什么意思？
orig.ident：通常包含示例身份（如果已知），但默认为“ SeuratProject”nCount_RNA：每个单元的UMI数量nFeature_RNA：每个细胞检测到的基因数量
读取多个样本for loop
实际上，您可能需要读取几个样本，才能使用我们前面讨论的2个功能之一（Read10X()或readMM()）读取数据。因此，为了使数据更有效地导入R中，我们可以使用一个for循环，该循环将为给定的每个输入插入一系列命令。
在R中，它具有以下结构/语法：
for (variable in input){
command1
command2
command3
}
for我们今天将要使用的循环将迭代两个样本“文件”，并对每个样本执行两个命令-（1）读入计数数据（Read10X()）和（2）从读入数据（CreateSeuratObject()）创建Seurat对象：
＃创建每个单独的Seurat对象对于每个样本
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_matrix
stim_raw_feature_bc_matrix我们可以使用在输入部分中将这些样本指定为
for loop向量的元素c()。我们将这些变量分配给变量，然后我们可以将其命名为任意变量（尝试为其赋予一个有意义的名称）。在此示例中，我们将变量称为file。在执行上述循环期间，
file将首先包含值“ ctrl_raw_feature_bc_matrix”，从头至尾一直执行命令assign()。接下来，它将包含值“ stim_raw_feature_bc_matrix”，并再次遍历所有命令。如果您有15个文件夹作为输入，而不是2个，则对于每个数据文件夹，以上代码将运行15次。
＃创建每个单独的Seurat对象
# 创建Seurat对象
for (file in c("ctrl_raw_feature_bc_matrix", "stim_raw_feature_bc_matrix")){
步骤2：读入数据作为输入
我们可以
for loop通过添加一行以读取数据来继续操作Read10X()：
- 在这里，我们需要指定文件的路径，因此我们将
data/使用paste0()函数将目录添加到示例文件夹名称的前面。
seurat_data <- Read10X(data.dir = paste0("data/", file))
步骤3：根据10倍计数数据创建Seurat对象
现在，我们可以使用
CreateSeuratObject()函数创建Seurat对象，并添加参数project，在其中可以添加样品名称。
seurat_obj <- CreateSeuratObject(counts =seurat_data, min.features = 100, project = file)
步骤4：根据样本将Seurat对象分配给新变量
最后一个命令
assign是将创建的Seurat对象（seurat_obj）添加到新变量。这样，当我们迭代并移至我们的下一个样本时，input我们将不会覆盖在先前迭代中创建的Seurat对象：
assign(file, seurat_obj)
}
现在我们已经创建了这两个对象，让我们快速看一下元数据以了解其外观：
#检查新的Seurat对象中的元数据
head(ctrl_raw_feature_bc_matrix@meta.data)
head(stim_raw_feature_bc_matrix@meta.data)
接下来，我们需要将这些对象合并到单个Seurat对象中。这将使两个样品组的质量控制步骤一起运行变得更加容易，并使我们能够轻松比较所有样品的数据质量。
我们可以使用merge()Seurat包中的函数来执行此操作：
# 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。如果我们查看合并对象的元数据，我们应该能够在行名中看到前缀：
#检查合并的对象是否具有适当的特定于样本的前缀head(merged_seurat@meta.data)tail(merged_seurat@meta.data)
全部代码：
1. 载入包
rm(list = ls())
gc()
library(SingleCellExperiment)
library(Seurat)
library(tidyverse)
library(Matrix)
library(scales)
library(cowplot)
library(RCurl)
library(ggsci)
2. 载入并合并数据
# How to read in 10X data for a single sample (output is a sparse matrix)
ctrl_counts <- Read10X(data.dir = "./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)
#ctrl <-Read10X(data.dir = "./ctrl_raw_feature_bc_matrix")%>%
# CreateSeuratObject(min.features = 100)
stim<- Read10X(data.dir = "./stim_raw_feature_bc_matrix/")%>%
CreateSeuratObject(min.features = 100)
scRNAlist1<-list(ctrl,stim)
#scRNAlist1<-c(ctrl,stim)
#给列表命名并保存数据
names(scRNAlist) <- samples_name
# Check the metadata in the new Seurat objects
head(scRNAlist$ctrl@meta.data)
tail(scRNAlist$stim@meta.data)
####根据原文，个人更改后的代码
file<-list.files("./data/")
samples_name = c('ctrl', 'stim')
for (file in sample){
seurat_data <- Read10X(data.dir = paste0("data/", file))
seurat_obj <- CreateSeuratObject(counts = seurat_data,
min.features = 100,
project = samples_name)
assign(samples_name, seurat_obj)#assign是将创建的Seurat对象（seurat_obj）添加到新变量
}
其实我在这里更推荐吴晓琪老师的代码，简洁高效：
dir <- dir("data/")
dir <- paste0("data/", dir)
samples_name = c('ctrl', 'stim')
##使用循环命令批量创建seurat对象
scRNAlist <- list()
for(i in 1:length(dir)){
#Insufficient data values to produce 24 bins.
counts <- Read10X(data.dir = dir[i])
scRNAlist[[i]] <- Read10X(data.dir = dir[i]) %>%
CreateSeuratObject(project=samples_name[i],
#min.cells=3,#可以指定每个样本最少的细胞数目
min.features = 100)
#给细胞barcode加个前缀，防止合并后barcode重名
scRNAlist[[i]] <- RenameCells(scRNAlist[[i]], add.cell.id = samples_name[i])
}
#给列表命名并保存数据
names(scRNAlist) <- samples_name
＃检查合并的对象是否具有适当的特定于样本的前缀
head(scRNAlist$ctrl@meta.data)
head(scRNAlist$stim@meta.data)
数据合并
merged_seurat<- merge(scRNAlist[[1]], scRNAlist[2:length(scRNAlist)])
＃检查合并的对象
head(merged_seurat@meta.data)
tail(merged_seurat@meta.data)
