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

在第三節(jié)作者展示了數(shù)據(jù)載入和合并，作者代碼如下，我略微做了更改。

學(xué)習(xí)目標(biāo)：

• 了解如何從單細(xì)胞RNA序列實(shí)驗(yàn)中引入數(shù)據(jù)

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

image

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

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

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

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

image

圖片來源：Kang等，2017

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

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

注意：此數(shù)據(jù)集的計(jì)數(shù)也可以從10X Genomics免費(fèi)獲得，并用作Seurat教程的一部分。

元數(shù)據(jù)

除了原始數(shù)據(jù)外，我們還需要收集有關(guān)數(shù)據(jù)的信息。這就是所謂的元數(shù)據(jù)。只是探索數(shù)據(jù)，但是如果我們不知道該數(shù)據(jù)所源自的樣本，那沒有意義。

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

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

• 樣品在Illumina NextSeq 500上測(cè)序

• 將來自八名狼瘡患者的PBMC樣本分別分成兩等份。

• 用100 U / mL重組IFN-β激活PBMC一份，持續(xù)6小時(shí)。

• 第二等分試樣未經(jīng)處理。

• 6小時(shí)后，將每種條件的8個(gè)樣品一起收集到兩個(gè)最終的收集池中（受激細(xì)胞和對(duì)照細(xì)胞）。我們將處理這兩個(gè)合并的樣本。

• 分別鑒定了12138和12167個(gè)細(xì)胞（去除雙峰后）作為對(duì)照樣品和刺激后的合并樣品。

• 由于樣本是PBMC，因此我們期望有免疫細(xì)胞，例如：

• B細(xì)胞

• T細(xì)胞

• NK細(xì)胞

• 單核細(xì)胞

• 巨噬細(xì)胞

• 可能是巨核細(xì)胞

建議您對(duì)執(zhí)行QC之前希望在數(shù)據(jù)集中看到的像元類型有一些期望。這將告知您是否具有低復(fù)雜度的細(xì)胞類型（許多基因的轉(zhuǎn)錄本）或線粒體表達(dá)水平較高的細(xì)胞。這將使我們能夠在分析工作流程中考慮這些生物學(xué)因素。

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

設(shè)置R環(huán)境

研究涉及大量數(shù)據(jù)的最重要部分之一就是如何最好地對(duì)其進(jìn)行管理。我們傾向于對(duì)分析進(jìn)行優(yōu)先級(jí)排序，但是數(shù)據(jù)管理中還有許多其他重要方面經(jīng)常被人們忽略，因此對(duì)新數(shù)據(jù)的初學(xué)者一見傾心。該HMS數(shù)據(jù)管理工作組，討論深入一些事情要考慮超出數(shù)據(jù)創(chuàng)建和分析。

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

打開RStudio并創(chuàng)建一個(gè)名為的新R項(xiàng)目single_cell_rnaseq。然后，創(chuàng)建以下目錄：

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

下載資料

右鍵單擊下面的鏈接，以將每個(gè)樣本的輸出文件夾從Cell Ranger下載到data文件夾中：

• 對(duì)照樣品

• 受激樣品

現(xiàn)在，讓我們通過雙擊解壓縮剛剛下載的兩個(gè)壓縮文件夾，以便我們可以在RStudio中查看它們的內(nèi)容。

新項(xiàng)目

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

＃ 2020年2月

＃ HBC單細(xì)胞RNA-seq的流程

＃單細(xì)胞RNA序列分析-QC

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

image

加載庫

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

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

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

無論處理單細(xì)胞RNA-seq序列數(shù)據(jù)的技術(shù)或流程如何，其輸出通常都是相同的。也就是說，對(duì)于每個(gè)單獨(dú)的樣本，您將擁有以下三個(gè)文件：

1. 具有細(xì)胞ID的文件，代表所有量化的細(xì)胞

2. 一個(gè)帶有基因ID的文件，代表所有量化的基因

3. 一個(gè)計(jì)數(shù)的矩陣每個(gè)基因的每一個(gè)細(xì)胞

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

1 barcodes.tsv

這是一個(gè)文本文件，其中包含該樣品存在的所有細(xì)胞條形碼。條形碼以矩陣文件中顯示的數(shù)據(jù)順序列出（即，這些是列名）。

cell_ids_new.png (384×788)

2 features.tsv

這是一個(gè)文本文件，其中包含定量基因的標(biāo)識(shí)符。標(biāo)識(shí)符的來源可能會(huì)有所不同，具體取決于您在定量方法中使用的參考文獻(xiàn)（即Ensembl，NCBI，UCSC），但最常見的是官方基因符號(hào)。這些基因的順序與矩陣文件中行的順序相對(duì)應(yīng)（即，這些是行名）。

image

3. matrix.mtx

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

image

將此數(shù)據(jù)加載到R中要求我們使用函數(shù)，這些函數(shù)使我們能夠?qū)⑦@三個(gè)文件有效地組合為一個(gè)計(jì)數(shù)矩陣。但是，我們將使用創(chuàng)建函數(shù)來創(chuàng)建稀疏矩陣，而不是創(chuàng)建規(guī)則的矩陣數(shù)據(jù)結(jié)構(gòu)，以改善處理龐大計(jì)數(shù)矩陣所需的空間，內(nèi)存和CPU。

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

1. readMM()：此函數(shù)來自Matrix包，它將把我們的標(biāo)準(zhǔn)矩陣變成稀疏矩陣。該features.tsv文件barcodes.tsv必須首先單獨(dú)加載到R中，然后再合并。有關(guān)如何執(zhí)行此操作的特定代碼和說明，請(qǐng)參閱我們的其他材料。

2. Read10X()：此功能來自Seurat軟件包，并將使用Cell Ranger的輸出目錄作為輸入。這樣，不需要加載單個(gè)文件，而是該函數(shù)將加載并將它們合并為一個(gè)稀疏矩陣。我們將使用此功能加載數(shù)據(jù)！

讀取一個(gè)樣本（read10X()）

當(dāng)使用10X數(shù)據(jù)及其專有軟件Cell Ranger時(shí)，您將始終有一個(gè)outs目錄。在此目錄中，您可以找到許多不同的文件，包括：

• web_summary.html：該報(bào)告探討了不同的QC指標(biāo)，包括映射指標(biāo)，過濾閾值，過濾后的估計(jì)細(xì)胞數(shù)以及過濾后每個(gè)細(xì)胞的讀取數(shù)和基因數(shù)的信息。

• BAM對(duì)齊文件：用于可視化映射的讀取和重新創(chuàng)建FASTQ文件的文件（如果需要）

• filtered_feature_bc_matrix：包含使用Cell Ranger過濾的數(shù)據(jù)構(gòu)造計(jì)數(shù)矩陣所需的所有文件的文件夾

• raw_feature_bc_matrix：包含使用原始未過濾數(shù)據(jù)構(gòu)建計(jì)數(shù)矩陣所需的所有文件的文件夾

我們主要對(duì)感興趣，raw_feature_bc_matrix因?yàn)槲覀兿Ｍ麍?zhí)行自己的質(zhì)量控制和過濾，同時(shí)考慮到我們的實(shí)驗(yàn)/生物學(xué)系統(tǒng)的生物學(xué)特性。

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

ctrl_counts <- Read10X(data.dir = "data/ctrl_raw_feature_bc_matrix")
＃將計(jì)數(shù)矩陣轉(zhuǎn)換為Seurat對(duì)象（輸出是Seurat對(duì)象）
ctrl <- CreateSeuratObject(counts = ctrl_counts,min.features = 100)                           

注意：該min.features參數(shù)指定每個(gè)細(xì)胞需要檢測(cè)的最小基因數(shù)量。此參數(shù)將過濾掉質(zhì)量較差的單元，這些單元可能只是封裝了隨機(jī)條形碼而沒有任何單元。通常，檢測(cè)不到少于100個(gè)基因的細(xì)胞不考慮進(jìn)行分析。

當(dāng)您使用該Read10X()功能讀取數(shù)據(jù)時(shí)，Seurat會(huì)為每個(gè)細(xì)胞自動(dòng)創(chuàng)建一些元數(shù)據(jù)。此信息存儲(chǔ)在Seurat對(duì)象的meta.data插槽中（請(qǐng)參閱下面的注釋中的更多信息）。

Seurat對(duì)象是一個(gè)自定義的類似列表的對(duì)象，具有定義明確的空間來存儲(chǔ)特定的信息/數(shù)據(jù)。您可以在此鏈接中找到有關(guān)Seurat對(duì)象中插槽的更多信息。

＃探索元數(shù)據(jù)
head(ctrl@meta.data)

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

• orig.ident：通常包含示例身份（如果已知），但默認(rèn)為“ SeuratProject”

• nCount_RNA：每個(gè)單元的UMI數(shù)量

• nFeature_RNA：每個(gè)細(xì)胞檢測(cè)到的基因數(shù)量

讀取多個(gè)樣本for loop

實(shí)際上，您可能需要讀取幾個(gè)樣本，才能使用我們前面討論的2個(gè)功能之一（Read10X()或readMM()）讀取數(shù)據(jù)。因此，為了使數(shù)據(jù)更有效地導(dǎo)入R中，我們可以使用一個(gè)for循環(huán)，該循環(huán)將為給定的每個(gè)輸入插入一系列命令。

在R中，它具有以下結(jié)構(gòu)/語法：

for (variable in input){
    command1
    command2
    command3
}

for我們今天將要使用的循環(huán)將迭代兩個(gè)樣本“文件”，并對(duì)每個(gè)樣本執(zhí)行兩個(gè)命令-（1）讀入計(jì)數(shù)數(shù)據(jù)（Read10X()）和（2）從讀入數(shù)據(jù)（CreateSeuratObject()）創(chuàng)建Seurat對(duì)象：

＃創(chuàng)建每個(gè)單獨(dú)的Seurat對(duì)象對(duì)于每個(gè)樣本
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：指定輸入

對(duì)于此數(shù)據(jù)集，我們有兩個(gè)樣本想要為以下對(duì)象創(chuàng)建Seurat對(duì)象：

• ctrl_raw_feature_bc_matrix

• stim_raw_feature_bc_matrix

我們可以使用在輸入部分中將這些樣本指定為for loop向量的元素c()。我們將這些變量分配給變量，然后我們可以將其命名為任意變量（嘗試為其賦予一個(gè)有意義的名稱）。在此示例中，我們將變量稱為file。

在執(zhí)行上述循環(huán)期間，file將首先包含值“ ctrl_raw_feature_bc_matrix”，從頭至尾一直執(zhí)行命令assign()。接下來，它將包含值“ stim_raw_feature_bc_matrix”，并再次遍歷所有命令。如果您有15個(gè)文件夾作為輸入，而不是2個(gè)，則對(duì)于每個(gè)數(shù)據(jù)文件夾，以上代碼將運(yùn)行15次。

＃創(chuàng)建每個(gè)單獨(dú)的Seurat對(duì)象

# 創(chuàng)建Seurat對(duì)象
for (file in c("ctrl_raw_feature_bc_matrix", "stim_raw_feature_bc_matrix")){

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

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

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

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

步驟3：根據(jù)10倍計(jì)數(shù)數(shù)據(jù)創(chuàng)建Seurat對(duì)象

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

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

步驟4：根據(jù)樣本將Seurat對(duì)象分配給新變量

最后一個(gè)命令assign是將創(chuàng)建的Seurat對(duì)象（seurat_obj）添加到新變量。這樣，當(dāng)我們迭代并移至我們的下一個(gè)樣本時(shí)，input我們將不會(huì)覆蓋在先前迭代中創(chuàng)建的Seurat對(duì)象：

        assign(file, seurat_obj)
}

現(xiàn)在我們已經(jīng)創(chuàng)建了這兩個(gè)對(duì)象，讓我們快速看一下元數(shù)據(jù)以了解其外觀：

#檢查新的Seurat對(duì)象中的元數(shù)據(jù)
head(ctrl_raw_feature_bc_matrix@meta.data)
head(stim_raw_feature_bc_matrix@meta.data)

接下來，我們需要將這些對(duì)象合并到單個(gè)Seurat對(duì)象中。這將使兩個(gè)樣品組的質(zhì)量控制步驟一起運(yùn)行變得更加容易，并使我們能夠輕松比較所有樣品的數(shù)據(jù)質(zhì)量。

我們可以使用merge()Seurat包中的函數(shù)來執(zhí)行此操作：

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

因?yàn)橄嗤募?xì)胞ID可以用于不同的樣本，所以我們使用參數(shù)為每個(gè)細(xì)胞ID添加一個(gè)特定于樣本的前綴add.cell.id。如果我們查看合并對(duì)象的元數(shù)據(jù)，我們應(yīng)該能夠在行名中看到前綴：

#檢查合并的對(duì)象是否具有適當(dāng)?shù)奶囟ㄓ跇颖镜那熬Yhead(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. 載入并合并數(shù)據(jù)

# 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)
#給列表命名并保存數(shù)據(jù)
names(scRNAlist) <- samples_name
# Check the metadata in the new Seurat objects
head(scRNAlist$ctrl@meta.data)
tail(scRNAlist$stim@meta.data)

####根據(jù)原文，個(gè)人更改后的代碼
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是將創(chuàng)建的Seurat對(duì)象（seurat_obj）添加到新變量
}

其實(shí)我在這里更推薦吳曉琪老師的代碼，簡(jiǎn)潔高效：

dir <- dir("data/")
dir <- paste0("data/", dir)
samples_name = c('ctrl', 'stim')

##使用循環(huán)命令批量創(chuàng)建seurat對(duì)象
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,#可以指定每個(gè)樣本最少的細(xì)胞數(shù)目
                   min.features = 100)
#給細(xì)胞barcode加個(gè)前綴，防止合并后barcode重名
scRNAlist[[i]] <- RenameCells(scRNAlist[[i]], add.cell.id = samples_name[i])
}
#給列表命名并保存數(shù)據(jù)
names(scRNAlist) <- samples_name

＃檢查合并的對(duì)象是否具有適當(dāng)?shù)奶囟ㄓ跇颖镜那熬Y
head(scRNAlist$ctrl@meta.data)
head(scRNAlist$stim@meta.data)

數(shù)據(jù)合并

merged_seurat<- merge(scRNAlist[[1]], scRNAlist[2:length(scRNAlist)])
＃檢查合并的對(duì)象
head(merged_seurat@meta.data)
tail(merged_seurat@meta.data)

image.png

image.png

image.png

image.png

image.png