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ù)據(jù)
單細(xì)胞RNA序列：質(zhì)量控制設(shè)置
定量基因表達后，我們需要將此數(shù)據(jù)帶入R中以生成用于執(zhí)行QC的指標(biāo)篡诽。在本課程中反浓，我們將討論可以預(yù)期的計數(shù)數(shù)據(jù)格式间涵，以及如何將其讀入R，以便我們繼續(xù)進行工作流程中的QC步驟耙厚。我們還將討論將要使用的數(shù)據(jù)集和相關(guān)的元數(shù)據(jù)。
探索示例數(shù)據(jù)集
我們將使用單細(xì)胞RNA-seq數(shù)據(jù)集岔霸，該數(shù)據(jù)集是Kang等人在2017年進行的一項較大研究的一部分薛躬。在本文中，作者提出了一種計算算法呆细，該算法利用遺傳變異（eQTL）來確定每個包含單個細(xì)胞的液滴的遺傳同一性（單一）型宝，并鑒定包含來自不同個體的兩個細(xì)胞的液滴（雙重）。
用于測試其算法的數(shù)據(jù)由八名狼瘡患者的合并外周血單個核細(xì)胞（PBMC）組成絮爷，分為對照和干擾素β治療（刺激）條件趴酣。
圖片來源：Kang等，2017
原始數(shù)據(jù)
該數(shù)據(jù)集可在GEO（GSE96583）上獲得坑夯，但是可用的計數(shù)矩陣缺少線粒體讀數(shù)岖寞，因此我們從SRA（SRP102802）下載了BAM文件。這些BAM文件被轉(zhuǎn)換回FASTQ文件柜蜈，然后通過Cell Ranger運行以獲得我們將要使用的計數(shù)數(shù)據(jù)仗谆。
注意：此數(shù)據(jù)集的計數(shù)也可以從10X Genomics免費獲得，并用作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上測序
將來自八名狼瘡患者的PBMC樣本分別分成兩等份捷绒。
用100 U / mL重組IFN-β激活PBMC一份瑰排，持續(xù)6小時。
第二等分試樣未經(jīng)處理暖侨。
6小時后椭住，將每種條件的8個樣品一起收集到兩個最終的收集池中（受激細(xì)胞和對照細(xì)胞）。我們將處理這兩個合并的樣本字逗。
分別鑒定了12138和12167個細(xì)胞（去除雙峰后）作為對照樣品和刺激后的合并樣品京郑。
由于樣本是PBMC，因此我們期望有免疫細(xì)胞葫掉，例如：
B細(xì)胞
T細(xì)胞
NK細(xì)胞
單核細(xì)胞
巨噬細(xì)胞
可能是巨核細(xì)胞
建議您對執(zhí)行QC之前希望在數(shù)據(jù)集中看到的像元類型有一些期望些举。這將告知您是否具有低復(fù)雜度的細(xì)胞類型（許多基因的轉(zhuǎn)錄本）或線粒體表達水平較高的細(xì)胞。這將使我們能夠在分析工作流程中考慮這些生物學(xué)因素俭厚。
預(yù)期上述細(xì)胞類型都不具有低復(fù)雜性或預(yù)期具有高線粒體含量户魏。
設(shè)置R環(huán)境
研究涉及大量數(shù)據(jù)的最重要部分之一就是如何最好地對其進行管理。我們傾向于對分析進行優(yōu)先級排序挪挤，但是數(shù)據(jù)管理中還有許多其他重要方面經(jīng)常被人們忽略叼丑，因此對新數(shù)據(jù)的初學(xué)者一見傾心。該HMS數(shù)據(jù)管理工作組扛门，討論深入一些事情要考慮超出數(shù)據(jù)創(chuàng)建和分析鸠信。
數(shù)據(jù)管理的一個重要方面是組織。對于您進行和分析數(shù)據(jù)的每個實驗论寨，通過創(chuàng)建計劃的存儲空間（目錄結(jié)構(gòu)）來組織起來被認(rèn)為是最佳實踐星立。我們將對單細(xì)胞分析進行此操作。
打開RStudio并創(chuàng)建一個名為的新R項目single_cell_rnaseq葬凳。然后绰垂，創(chuàng)建以下目錄：
single_cell_rnaseq/
├── data
├── results
└── figures
下載資料
右鍵單擊下面的鏈接，以將每個樣本的輸出文件夾從Cell Ranger下載到data文件夾中：
現(xiàn)在火焰，讓我們通過雙擊解壓縮剛剛下載的兩個壓縮文件夾辕坝，以便我們可以在RStudio中查看它們的內(nèi)容。
新項目
接下來荐健，打開一個新的Rscript文件酱畅，并從一些注釋開始以指示該文件將包含的內(nèi)容：
＃ 2020年2月
＃ HBC單細(xì)胞RNA-seq的流程
＃單細(xì)胞RNA序列分析-QC
將Rscript另存為quality_control.R。您的工作目錄應(yīng)如下所示：
加載庫
現(xiàn)在江场，我們可以加載必要的R包：
# 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ù)或流程如何纺酸，其輸出通常都是相同的。也就是說址否，對于每個單獨的樣本餐蔬，您將擁有以下三個文件：
具有細(xì)胞ID的文件碎紊，代表所有量化的細(xì)胞
一個帶有基因ID的文件，代表所有量化的基因
一個計數(shù)的矩陣每個基因的每一個細(xì)胞
我們可以通過單擊data/ctrl_raw_feature_bc_matrix文件夾瀏覽這些文件：
1 barcodes.tsv
這是一個文本文件樊诺，其中包含該樣品存在的所有細(xì)胞條形碼仗考。條形碼以矩陣文件中顯示的數(shù)據(jù)順序列出（即，這些是列名）词爬。
2 features.tsv
這是一個文本文件秃嗜，其中包含定量基因的標(biāo)識符。標(biāo)識符的來源可能會有所不同顿膨，具體取決于您在定量方法中使用的參考文獻（即Ensembl锅锨，NCBI，UCSC）恋沃，但最常見的是官方基因符號必搞。這些基因的順序與矩陣文件中行的順序相對應(yīng)（即，這些是行名）囊咏。
3. matrix.mtx
這是一個文本文件恕洲，其中包含計數(shù)值的矩陣。行與上面的基因ID相關(guān)聯(lián)梅割，列與細(xì)胞條形碼相對應(yīng)霜第。請注意，此矩陣中有許多零值炮捧。
將此數(shù)據(jù)加載到R中要求我們使用函數(shù)，這些函數(shù)使我們能夠?qū)⑦@三個文件有效地組合為一個計數(shù)矩陣惦银。但是咆课，我們將使用創(chuàng)建函數(shù)來創(chuàng)建稀疏矩陣，而不是創(chuàng)建規(guī)則的矩陣數(shù)據(jù)結(jié)構(gòu)扯俱，以改善處理龐大計數(shù)矩陣所需的空間书蚪，內(nèi)存和CPU。
讀取數(shù)據(jù)的不同方法包括：
readMM()：此函數(shù)來自Matrix包迅栅，它將把我們的標(biāo)準(zhǔn)矩陣變成稀疏矩陣殊校。該features.tsv文件barcodes.tsv必須首先單獨加載到R中，然后再合并读存。有關(guān)如何執(zhí)行此操作的特定代碼和說明为流，請參閱我們的其他材料。Read10X()：此功能來自Seurat軟件包让簿，并將使用Cell Ranger的輸出目錄作為輸入敬察。這樣，不需要加載單個文件尔当，而是該函數(shù)將加載并將它們合并為一個稀疏矩陣莲祸。我們將使用此功能加載數(shù)據(jù)！
讀取一個樣本（read10X()）
當(dāng)使用10X數(shù)據(jù)及其專有軟件Cell Ranger時，您將始終有一個outs目錄锐帜。在此目錄中田盈，您可以找到許多不同的文件，包括：
web_summary.html：該報告探討了不同的QC指標(biāo)缴阎，包括映射指標(biāo)允瞧，過濾閾值，過濾后的估計細(xì)胞數(shù)以及過濾后每個細(xì)胞的讀取數(shù)和基因數(shù)的信息药蜻。BAM對齊文件：用于可視化映射的讀取和重新創(chuàng)建FASTQ文件的文件（如果需要）
filtered_feature_bc_matrix：包含使用Cell Ranger過濾的數(shù)據(jù)構(gòu)造計數(shù)矩陣所需的所有文件的文件夾raw_feature_bc_matrix：包含使用原始未過濾數(shù)據(jù)構(gòu)建計數(shù)矩陣所需的所有文件的文件夾
我們主要對感興趣瓷式，raw_feature_bc_matrix因為我們希望執(zhí)行自己的質(zhì)量控制和過濾，同時考慮到我們的實驗/生物學(xué)系統(tǒng)的生物學(xué)特性语泽。
如果我們只有一個樣本贸典，我們可以生成計數(shù)矩陣，然后創(chuàng)建一個Seurat對象：
ctrl_counts <- Read10X(data.dir = "data/ctrl_raw_feature_bc_matrix")
＃將計數(shù)矩陣轉(zhuǎn)換為Seurat對象（輸出是Seurat對象）
ctrl <- CreateSeuratObject(counts = ctrl_counts,min.features = 100)
注意：該
min.features參數(shù)指定每個細(xì)胞需要檢測的最小基因數(shù)量踱卵。此參數(shù)將過濾掉質(zhì)量較差的單元廊驼，這些單元可能只是封裝了隨機條形碼而沒有任何單元。通常惋砂，檢測不到少于100個基因的細(xì)胞不考慮進行分析妒挎。
當(dāng)您使用該Read10X()功能讀取數(shù)據(jù)時，Seurat會為每個細(xì)胞自動創(chuàng)建一些元數(shù)據(jù)西饵。此信息存儲在Seurat對象的meta.data插槽中（請參閱下面的注釋中的更多信息）酝掩。
Seurat對象是一個自定義的類似列表的對象，具有定義明確的空間來存儲特定的信息/數(shù)據(jù)眷柔。您可以在此鏈接中找到有關(guān)Seurat對象中插槽的更多信息期虾。
＃探索元數(shù)據(jù)
head(ctrl@meta.data)
元數(shù)據(jù)列是什么意思？
orig.ident：通常包含示例身份（如果已知）驯嘱，但默認(rèn)為“ SeuratProject”nCount_RNA：每個單元的UMI數(shù)量nFeature_RNA：每個細(xì)胞檢測到的基因數(shù)量
讀取多個樣本for loop
實際上镶苞，您可能需要讀取幾個樣本，才能使用我們前面討論的2個功能之一（Read10X()或readMM()）讀取數(shù)據(jù)鞠评。因此茂蚓，為了使數(shù)據(jù)更有效地導(dǎo)入R中，我們可以使用一個for循環(huán)剃幌，該循環(huán)將為給定的每個輸入插入一系列命令聋涨。
在R中，它具有以下結(jié)構(gòu)/語法：
for (variable in input){
command1
command2
command3
}
for我們今天將要使用的循環(huán)將迭代兩個樣本“文件”负乡，并對每個樣本執(zhí)行兩個命令-（1）讀入計數(shù)數(shù)據(jù)（Read10X()）和（2）從讀入數(shù)據(jù)（CreateSeuratObject()）創(chuàng)建Seurat對象：
＃創(chuàng)建每個單獨的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：指定輸入
對于此數(shù)據(jù)集牛郑，我們有兩個樣本想要為以下對象創(chuàng)建Seurat對象：
ctrl_raw_feature_bc_matrix
stim_raw_feature_bc_matrix我們可以使用在輸入部分中將這些樣本指定為
for loop向量的元素c()。我們將這些變量分配給變量敬鬓，然后我們可以將其命名為任意變量（嘗試為其賦予一個有意義的名稱）淹朋。在此示例中笙各，我們將變量稱為file。在執(zhí)行上述循環(huán)期間础芍，
file將首先包含值“ ctrl_raw_feature_bc_matrix”杈抢，從頭至尾一直執(zhí)行命令assign()。接下來仑性，它將包含值“ stim_raw_feature_bc_matrix”惶楼，并再次遍歷所有命令栏笆。如果您有15個文件夾作為輸入劳曹，而不是2個，則對于每個數(shù)據(jù)文件夾你画，以上代碼將運行15次晨汹。
＃創(chuàng)建每個單獨的Seurat對象
# 創(chuàng)建Seurat對象
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倍計數(shù)數(shù)據(jù)創(chuàng)建Seurat對象
現(xiàn)在剥扣，我們可以使用
CreateSeuratObject()函數(shù)創(chuàng)建Seurat對象，并添加參數(shù)project铝穷，在其中可以添加樣品名稱钠怯。
seurat_obj <- CreateSeuratObject(counts =seurat_data, min.features = 100, project = file)
步驟4：根據(jù)樣本將Seurat對象分配給新變量
最后一個命令
assign是將創(chuàng)建的Seurat對象（seurat_obj）添加到新變量。這樣曙聂，當(dāng)我們迭代并移至我們的下一個樣本時晦炊，input我們將不會覆蓋在先前迭代中創(chuàng)建的Seurat對象：
assign(file, seurat_obj)
}
現(xiàn)在我們已經(jīng)創(chuàng)建了這兩個對象，讓我們快速看一下元數(shù)據(jù)以了解其外觀：
#檢查新的Seurat對象中的元數(shù)據(jù)
head(ctrl_raw_feature_bc_matrix@meta.data)
head(stim_raw_feature_bc_matrix@meta.data)
接下來宁脊，我們需要將這些對象合并到單個Seurat對象中断国。這將使兩個樣品組的質(zhì)量控制步驟一起運行變得更加容易，并使我們能夠輕松比較所有樣品的數(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"))
因為相同的細(xì)胞ID可以用于不同的樣本并思，所以我們使用參數(shù)為每個細(xì)胞ID添加一個特定于樣本的前綴add.cell.id庐氮。如果我們查看合并對象的元數(shù)據(jù)语稠，我們應(yīng)該能夠在行名中看到前綴：
#檢查合并的對象是否具有適當(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ù)原文，個人更改后的代碼
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對象（seurat_obj）添加到新變量
}
其實我在這里更推薦吳曉琪老師的代碼弄砍，簡潔高效：
dir <- dir("data/")
dir <- paste0("data/", dir)
samples_name = c('ctrl', 'stim')
##使用循環(huán)命令批量創(chuàng)建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,#可以指定每個樣本最少的細(xì)胞數(shù)目
min.features = 100)
#給細(xì)胞barcode加個前綴仙畦，防止合并后barcode重名
scRNAlist[[i]] <- RenameCells(scRNAlist[[i]], add.cell.id = samples_name[i])
}
#給列表命名并保存數(shù)據(jù)
names(scRNAlist) <- samples_name
＃檢查合并的對象是否具有適當(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)])
＃檢查合并的對象
head(merged_seurat@meta.data)
tail(merged_seurat@meta.data)
