03 载入并合并数据

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)组成,分为对照和干扰素β治疗(刺激)条件。

image

图片来源: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。您的工作目录应如下所示:

image

加载库

现在,我们可以加载必要的R包:

# Load libraries
library(SingleCellExperiment)
library(Seurat)
library(tidyverse)
library(Matrix)
library(scales)
library(cowplot)
library(RCurl)

加载单细胞RNA-seq计数数据

无论处理单细胞RNA-seq序列数据的技术或流程如何,其输出通常都是相同的。也就是说,对于每个单独的样本,您将拥有以下三个文件

  1. 具有细胞ID的文件,代表所有量化的细胞

  2. 一个带有基因ID的文件,代表所有量化的基因

  3. 一个计数的矩阵每个基因的每一个细胞

我们可以通过单击data/ctrl_raw_feature_bc_matrix文件夹浏览这些文件:

1 barcodes.tsv

这是一个文本文件,其中包含该样品存在的所有细胞条形码。条形码以矩阵文件中显示的数据顺序列出(即,这些是列名)。

2 features.tsv

这是一个文本文件,其中包含定量基因的标识符。标识符的来源可能会有所不同,具体取决于您在定量方法中使用的参考文献(即Ensembl,NCBI,UCSC),但最常见的是官方基因符号。这些基因的顺序与矩阵文件中行的顺序相对应(即,这些是行名)。

image

3. matrix.mtx

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

将此数据加载到R中要求我们使用函数,这些函数使我们能够将这三个文件有效地组合为一个计数矩阵。但是,我们将使用创建函数来创建稀疏矩阵,而不是创建规则的矩阵数据结构,以改善处理庞大计数矩阵所需的空间,内存和CPU。

读取数据的不同方法包括:

  1. readMM():此函数来自Matrix包,它将把我们的标准矩阵变成稀疏矩阵。该features.tsv文件barcodes.tsv必须首先单独加载到R中,然后再合并。有关如何执行此操作的特定代码和说明,请参阅我们的其他材料

  2. 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)
image.png
image.png
image.png
image.png
image.png

推荐阅读更多精彩内容