單細(xì)胞 RNA-seq：質(zhì)量控制設(shè)置

image

基因表達(dá)量化后，我們需要將此數(shù)據(jù)讀入 R 中，執(zhí)行 QC 的指標(biāo)。我們將討論可以預(yù)期的數(shù)據(jù)格式，以及如何將其讀入 R 以便我們可以繼續(xù)流程中的 QC 步驟。我們還將討論我們將使用的數(shù)據(jù)集和相關(guān)的metadata
。

探索示例數(shù)據(jù)集

我們將使用單細(xì)胞 RNA-seq 數(shù)據(jù)集，該數(shù)據(jù)集是Kang 等人于2017年進(jìn)行的一項(xiàng)較大研究的一部分。在本文中，作者提出了一種計算算法，該算法利用遺傳變異 (eQTL) 來確定每個包含單個細(xì)胞 (singlet) 的液滴的遺傳身份，并識別包含來自不同個體的兩個細(xì)胞 (doublet) 的液滴。

用于測試其算法的數(shù)據(jù)由來自八名狼瘡患者的外周血單核細(xì)胞 (PBMC) 組成，分為對照和干擾素 β 治療（刺激）條件。

image

圖片來源：Kang 等人，2017 年(https://www.nature.com/articles/nbt.4042)

原始數(shù)據(jù)

該數(shù)據(jù)集在 GEO ( GSE96583 ) 上可用，但是可用的計數(shù)矩陣缺少線粒體讀數(shù)，因此我們從 SRA ( SRP102802 )下載了 BAM 文件。這些 BAM 文件被轉(zhuǎn)換回 FASTQ 文件，然后通過Cell Ranger運(yùn)行以獲得計數(shù)數(shù)據(jù)。

注意： 此數(shù)據(jù)集的計數(shù)數(shù)據(jù)也可從 10X Genomics 免費(fèi)獲得，并在Seurat 教程中使用。

metadata(元數(shù)據(jù))

除了原始數(shù)據(jù)，我們還需要收集有關(guān)數(shù)據(jù)的信息；這稱為元數(shù)據(jù)。人們往往拿到數(shù)據(jù)就會直接來開始探索這些數(shù)據(jù)，但如果我們對這些數(shù)據(jù)的來源樣本一無所知，那就沒有太大的意義了。

下面提供了我們數(shù)據(jù)集的一些相關(guān)元數(shù)據(jù)：

• 使用 10X Genomics V2 版本化學(xué)試劑盒制備文庫

• 樣品在 Illumina NextSeq 500 上測序

• 將來自 8 名狼瘡患者的 PBMC 樣品各分成兩等份。

  • 用100 U/mL 重組 IFN-β 激活一份PBMC，持續(xù) 6 小時。
  • 第二等份不處理。
  • 6 小時后，將每種條件下的 8 個樣品混合在兩個最終池（受刺激細(xì)胞和對照細(xì)胞）中。我們將使用這兩個合并樣本。（我們沒有對樣本進(jìn)行多路分解，因?yàn)檎撐闹惺褂昧?SNP 基因型信息進(jìn)行多路分解，而這些數(shù)據(jù)的條形碼/樣本 ID 不容易獲得。通常，您會對每個單獨(dú)的樣本進(jìn)行多路分解和 QC，而不是合并樣本。 )

• 用于對照和刺激的合并樣品分別鑒定了 12138 和 12167 個細(xì)胞（去除雙聯(lián)體后）。

• 由于是 PBMC樣本，預(yù)計有以下免疫細(xì)胞，例如：

  • B細(xì)胞
  • T細(xì)胞
  • NK細(xì)胞
  • 單核細(xì)胞
  • 巨噬細(xì)胞
  • 可能有巨核細(xì)胞

建議您在執(zhí)行 QC 之前對您希望在數(shù)據(jù)集中看到的細(xì)胞類型有一些期望。這會告訴您是否有任何低復(fù)雜度的細(xì)胞類型（來自少數(shù)基因的大量轉(zhuǎn)錄本）或線粒體表達(dá)水平較高的細(xì)胞。我們將在分析工作流程中考慮這些生物因素。

預(yù)計上述細(xì)胞類型均不具有低復(fù)雜性或預(yù)計具有高線粒體含量。

設(shè)置R運(yùn)行環(huán)境

我們將在 RStudio 項(xiàng)目中工作。為了繼續(xù)學(xué)習(xí)，您應(yīng)該已經(jīng)下載了 R 項(xiàng)目。

如果您還沒有這樣做，可以使用此鏈接訪問該項(xiàng)目。

下載后，您應(yīng)該會在您的計算機(jī)上看到一個名為single_cell_rnaseq.zip的文件（可能在您的Downloads文件夾中）。

1. 解壓此文件。這將產(chǎn)生一個同名的文件夾。
2. 將文件夾移動到計算機(jī)分析的位置。
3. 打開文件夾。內(nèi)容將類似于下面的屏幕截圖。
4. 找到.Rproj file并雙擊它。這將打開帶有“single_cell_rnaseq”項(xiàng)目加載的 RStudio。

image

項(xiàng)目管理

研究中的大量數(shù)據(jù)最重要的部分之一是如何最好地管理它。我們傾向于優(yōu)先考慮分析，但數(shù)據(jù)管理的許多其他重要方面往往在第一次看到新數(shù)據(jù)中被忽視。該HMS數(shù)據(jù)管理工作組，深入討論一些超出數(shù)據(jù)創(chuàng)建和分析。

數(shù)據(jù)管理的另一個重點(diǎn)是組織。對于您處理和分析每個實(shí)驗(yàn)的數(shù)據(jù)，通過創(chuàng)建計劃的存儲空間（目錄結(jié)構(gòu)）來進(jìn)行組織被認(rèn)為是最佳實(shí)踐。我們將為我們的單細(xì)胞分析做到這一點(diǎn)。

查看您的項(xiàng)目空間，您會發(fā)現(xiàn)已經(jīng)為您設(shè)置了一個目錄結(jié)構(gòu)：

single_cell_rnaseq/
├── data
├── results
└── figures

WINDOWS OS 用戶注意事項(xiàng)： 當(dāng)您解壓后打開項(xiàng)目文件夾時，請檢查您是否有一個data文件夾，其子文件夾也稱為data. 如果是這種情況，請將子文件夾中的所有文件移動到父data文件夾中。

創(chuàng)建新的腳本

接下來，打開一個新的 Rscript 文件，開始寫一些注釋，以指示該文件將包含的內(nèi)容：

# July/August 2021
# HBC single-cell RNA-seq workshop

# Single-cell RNA-seq analysis - QC

將 Rscript 保存為quality_control.R. 您的工作目錄應(yīng)如下所示：

image

庫加載

現(xiàn)在，我們可以加載必要的庫：

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

加載單細(xì)胞 RNA-seq 計數(shù)數(shù)據(jù)

無論用于處理原始單細(xì)胞 RNA-seq 序列數(shù)據(jù)的技術(shù)或pipeline如何，輸出的表達(dá)的矩陣通常是相同的。也就是說，對于每個單獨(dú)的樣本，都擁有以下三個文件：

1. 帶有cell ID的文件，表示所有量化的細(xì)胞
2. 帶有gene ID的文件，代表所有量化的基因
3. 一個計數(shù)的矩陣每個基因?qū)?yīng)每一個細(xì)胞

我們可以通過單擊data/ctrl_raw_feature_bc_matrix文件夾來瀏覽這些文件：

1. barcodes.tsv

這是一個文本文件，其中包含該樣本的所有細(xì)胞barcode。條形碼按照矩陣文件中數(shù)據(jù)的順序列出（即這些是列名稱）。

image

2. features.tsv

這是一個文本文件，其中包含量化基因的標(biāo)識符。標(biāo)識符可能因您在量化方法中使用的參考的來源（即 Ensembl、NCBI、UCSC）而異，但大多數(shù)情況下，這些是官方基因符號。這些基因的順序?qū)?yīng)于矩陣文件中行的順序（即這些是行名稱）。

image

3. matrix.mtx

這是一個包含計數(shù)值矩陣的文本文件。行與上面的基因 ID 相關(guān)聯(lián)，列對應(yīng)于細(xì)胞條形碼。請注意，此矩陣中有許多零值。

image

將此數(shù)據(jù)加載到 R 中需要我們使用將這三個文件組合成單個計數(shù)矩陣的函數(shù)。然而，我們將使用的函數(shù)不是創(chuàng)建常規(guī)矩陣數(shù)據(jù)結(jié)構(gòu)，而是創(chuàng)建一個稀疏矩陣，以減少使用我們龐大的計數(shù)矩陣所需的內(nèi)存 (RAM)、處理容量 (CPU) 和存儲量。

不同的讀入數(shù)據(jù)方法包括：

1. readMM()：這個函數(shù)來自Matrix包，它將把我們的標(biāo)準(zhǔn)矩陣轉(zhuǎn)換成一個稀疏矩陣。features.tsv和barcodes.tsv文件必須先被單獨(dú)地加載到R，然后再將它們組合。有關(guān)如何執(zhí)行此操作的特定代碼和說明，請參閱這些附加材料。
2. Read10X()：此函數(shù)來自Seurat包，將直接使用 Cell Ranger 輸出目錄作為輸入。使用這種方法不需要加載單個文件，而是該函數(shù)將加載它們并將它們組合成一個稀疏矩陣。我們將使用這個函數(shù)來加載我們的數(shù)據(jù)！

讀取單個樣本

使用 Cell Ranger 軟件處理 10X 數(shù)據(jù)后，您將擁有一個outs目錄（始終）。在此目錄中，您將找到許多不同的文件，包括下面列出的文件：

• web_summary.html：不同 QC 指標(biāo)的報告，包括映射指標(biāo)、過濾閾值、過濾后的估計細(xì)胞數(shù)以及有關(guān)過濾后每個細(xì)胞的read數(shù)和基因數(shù)的信息。
• BAM 對齊文件：用于可視化映射read和重新創(chuàng)建 FASTQ 文件的文件（如果需要）
• filtered_feature_bc_matrix:包含使用 Cell Ranger 過濾的數(shù)據(jù)構(gòu)建計數(shù)矩陣所需的所有文件的文件夾
• raw_feature_bc_matrix：包含使用原始未過濾數(shù)據(jù)構(gòu)建計數(shù)矩陣所需的所有文件的文件夾

雖然 Cell Ranger 對表達(dá)計數(shù)執(zhí)行過濾（請參閱下面的注釋），但我們希望執(zhí)行我們自己的 QC 和過濾，因?yàn)槲覀冃枰紤]自己的實(shí)驗(yàn)/生物學(xué)問題。鑒于此，我們只對raw_feature_bc_matrixCell Ranger 輸出中的****文件夾****感興趣。

注意：為什么我們不使用filtered_feature_bc_matrix文件夾？該filtered_feature_bc_matrix應(yīng)用通過細(xì)胞內(nèi)部過濾條件，我們沒有什么細(xì)胞的控制，以保持或放棄。

Cell Ranger 在生成filtered_feature_bc_matrix時執(zhí)行的過濾通常很好；然而，有時數(shù)據(jù)的質(zhì)量可能非常高，而 Seurat 過濾過程可以去除高質(zhì)量的細(xì)胞。此外，通常最好探索您自己的數(shù)據(jù)，同時考慮在過濾過程中應(yīng)用閾值的實(shí)驗(yàn)的生物學(xué)特性。例如，如果您希望數(shù)據(jù)集中的特定細(xì)胞類型更小和/或轉(zhuǎn)錄活性不如數(shù)據(jù)集中的其他細(xì)胞類型，則這些細(xì)胞有可能被過濾掉。然而，在 Cell Ranger v3 中，將考慮不同大小的細(xì)胞（例如，腫瘤與浸潤淋巴細(xì)胞），現(xiàn)在可能無法根據(jù)需要過濾盡可能多的低質(zhì)量細(xì)胞。

如果我們只有一個樣本，我們可以生成計數(shù)矩陣，然后創(chuàng)建一個 Seurat 對象：

Seurat 對象是一個自定義的類似列表的對象，它具有明確定義的空間來存儲特定的信息/數(shù)據(jù)。您可以在此鏈接中找到有關(guān) 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參數(shù)指定每個細(xì)胞需要檢測的最小基因數(shù)。這將過濾掉質(zhì)量較差的細(xì)胞，這些細(xì)胞可能只封裝了隨機(jī)barcode，而沒有任何細(xì)胞存在。通常，檢測到的基因少于 100 個的細(xì)胞不用于分析。

當(dāng)您使用該Read10X()函數(shù)讀入數(shù)據(jù)時，Seurat 會自動為每個細(xì)胞創(chuàng)建一些元數(shù)據(jù)。此信息存儲在Seurat 對象內(nèi)的meta.data中。

# Explore the metadata
head(ctrl@meta.data)

元數(shù)據(jù)的列是什么意思？

• orig.ident：如果已知，這通常包含樣本身份，但將默認(rèn)為“SeuratProject”
• nCount_RNA: 每個細(xì)胞的 UMI 數(shù)量
• nFeature_RNA: 每個細(xì)胞檢測到的基因數(shù)量

for loop讀取多個樣本

在現(xiàn)實(shí)中，您可能有多個樣本需要讀取數(shù)據(jù)，如果一次讀取一個樣本，可能會變得乏味且容易出錯。因此，為了使數(shù)據(jù)導(dǎo)入 R 更有效，我們可以使用for循環(huán)，該循環(huán)將為每個輸入給定一系列命令，并為我們的每個樣本創(chuàng)建 seurat 對象。

在 R 中，for循環(huán)具有以下結(jié)構(gòu)/語法：

## DO NOT RUN

for (variable in input){
    command1
    command2
    command3
}

今天我們將使用它來遍歷兩個示例文件夾并為每個示例執(zhí)行兩個命令，就像我們上面對單個示例所做的一樣 。
(1) 讀取計數(shù)數(shù)據(jù) ( Read10X())
(2) 從讀取的數(shù)據(jù)中創(chuàng)建 Seurat 對象數(shù)據(jù) ( 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：指定輸入

對于這個數(shù)據(jù)集，我們有兩個樣本和兩個關(guān)聯(lián)的文件夾，我們將創(chuàng)建兩個 Seurat 對象的輸入：

• ctrl_raw_feature_bc_matrix
• stim_raw_feature_bc_matrix

我們可以在輸入部分指定這些示例文件夾for loop作為向量的元素使用c(). 我們將這些分配給一個變量，我們可以隨意調(diào)用該變量（嘗試給它一個有意義的名稱）。在這個例子中，我們調(diào)用了變量 file。

在上述循環(huán)的執(zhí)行過程中，file將首先包含值“ctrl_raw_feature_bc_matrix”，通過命令一直運(yùn)行到assign()。接下來，它將包含值“stim_raw_feature_bc_matrix”并再次運(yùn)行所有命令。如果您有 15 個文件夾作為輸入，而不是 2 個，則上述代碼將針對您的每個數(shù)據(jù)文件夾運(yùn)行 15 次。

## DO NOT RUN

# Create each individual Seurat object
for (file in c("ctrl_raw_feature_bc_matrix", "stim_raw_feature_bc_matrix")){

步驟 2：輸入讀入數(shù)據(jù)

我們可以for loop通過添加一行來繼續(xù)我們的讀取數(shù)據(jù)Read10X()：

• 在這里，我們需要指定文件的路徑，因此我們將使用paste0()函數(shù)將data/目錄添加到我們的示例文件夾名稱中。

## DO NOT RUN

        seurat_data <- Read10X(data.dir = paste0("data/", file))

步驟 3：從 10X 計數(shù)數(shù)據(jù)創(chuàng)建 Seurat 對象

現(xiàn)在，我們可以使用該CreateSeuratObject()函數(shù)創(chuàng)建 Seurat 對象，添加參數(shù)project，我們可以在其中添加樣本名稱。

## DO NOT RUN

        seurat_obj <- CreateSeuratObject(counts = seurat_data, 
                                         min.features = 100, 
                                         project = file)        

第 4 步：將 Seurat 對象分配給基于樣本的新變量

最后一個命令assign是將 Seurat 對象創(chuàng)建 ( seurat_obj) 到一個新變量。這樣，當(dāng)我們迭代并移動到我們的下一個樣本時，我們input不會覆蓋在前一次迭代中創(chuàng)建的 Seurat 對象：

## DO NOT RUN
  
        assign(file, seurat_obj)
}

現(xiàn)在我們已經(jīng)創(chuàng)建了這兩個對象，讓我們快速瀏覽一下元數(shù)據(jù)：

# 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 對象。同時運(yùn)行兩個樣品組的 QC 步驟，我們能夠輕松比較所有樣品的數(shù)據(jù)質(zhì)量。

我們可以使用Seurat 包中的merge()函數(shù)來做到這一點(diǎn)：

# 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"))

由于相同的細(xì)胞 ID 可用于不同的樣本，因此我們使用參數(shù)為每個細(xì)胞 ID添加特定的樣本的前綴add.cell.id。如果我們查看合并對象的元數(shù)據(jù)，我們應(yīng)該能夠看到行名中的前綴：

# Check that the merged object has the appropriate sample-specific prefixes
head(merged_seurat@meta.data)
tail(merged_seurat@meta.data)