單細胞轉(zhuǎn)錄組測序進展迅速曹质，伴隨而來的是許許多多的內(nèi)容與講義，很多課程都老長涩禀。料滥。。嗨像我這樣莫有耐心的小孩子艾船，能完整看完的只有一個 —— Single-cell RNA-seq analysis workshop(https://github.com/hbctraining/scRNA-seq)葵腹。我將把教程內(nèi)容盡量匯總，你需要做的就是點擊收藏慢慢看屿岂！
（一）Why single-cell RNA-seq
在整個人體組織中践宴，細胞類型、狀態(tài)和相互作用非常多樣爷怀。為了更好地了解這些組織存在的細胞類型阻肩，單細胞RNA-seq（scRNA-seq）提供了在單細胞水平上表達基因的信息。
單細胞轉(zhuǎn)錄組測序可用于 ：
- 探索組織中存在哪些細胞類型
- 識別未知/稀有的細胞類型或狀態(tài)
- 闡明分化過程中或跨時間或狀態(tài)的基因表達變化
- 識別在不同條件（例如治療或疾苍耸凇）下特定細胞類型中差異表達的基因
- 探索細胞類型之間的表達變化烤惊，同時納入空間、調(diào)控和/或蛋白質(zhì)信息
常見方法包括：
Challenges of scRNA-seq analysis
在scRNA-seq之前吁朦，我們通常使用bulk RNA-seq進行轉(zhuǎn)錄組分析柒室。bulk RNA-seq是一種直接比較細胞的平均表達值的方法，在尋找疾病生物標志物逗宜，或者不關(guān)心樣品中大量細胞異質(zhì)性的情況下雄右，這可能是最佳方法。
盡管bulk RNA-seq可以探索不同條件（例如治療或疾卜慕病）之間基因表達的差異擂仍，但無法充分捕獲細胞水平的差異。例如熬甚，在下面的圖像中逢渔，如果進行bulk分析（左），我們將無法檢測到基因A和基因B的表達之間的正確關(guān)聯(lián)乡括。但是肃廓，如果我們按細胞類型或細胞狀態(tài)正確地對細胞進行分組冲簿，我們可以看到基因之間的正確相關(guān)性。
scRNA-seq也有一定的局限性亿昏，除了制樣和建庫價格高昂外峦剔，它在數(shù)據(jù)分析中也具有一定的復雜性，包括：
- 數(shù)據(jù)量大
- 細胞的測序深度低
- 細胞/樣品之間的技術(shù)差異
- 跨細胞/樣品的生物變異性
數(shù)據(jù)量大
scRNA-seq實驗的數(shù)據(jù)來自捕獲的成千上萬個細胞角钩，也就是數(shù)十萬條reads吝沫，需要更多的內(nèi)存和存儲空間。
細胞的測序深度低
基于液滴的scRNA-seq方法的測序深度較淺递礼，通常每個細胞只能檢測到轉(zhuǎn)錄組的10-50％惨险。這導致細胞中許多基因的計數(shù)為零。但是脊髓，在特定的細胞中辫愉，基因的零計數(shù)可能意味著該基因未表達或僅表示該轉(zhuǎn)錄本未被檢測到。在整個細胞中将硝，具有較高表達水平的基因測到0值的機率較低恭朗。由于此特征，在任何細胞中都不會檢測到許多基因依疼，并且細胞之間的基因表達差異很大痰腮。
Zero-inflated?：scRNA-seq data is often referred to as zero-inflated; however, recent analyses suggest that it does not contain more zeros than what would be expected given the sequencing depth （Valentine Svensson’s blog post：http://www.nxn.se/valent/2017/11/16/droplet-scrna-seq-is-not-zero-inflated）.
跨細胞/樣品的生物變異性
我們不感興趣的某些生物變異可能導致細胞之間的基因表達比實際生物細胞的類型/狀態(tài)更為相似/不同，并掩蓋細胞類型律罢。這些變異（除非實驗研究的一部分）包括（以下為原文）：
- Transcriptional bursting: Gene transcription is not turned on all of the time for all genes. Time of harvest will determine whether gene is on or off in each cell.
- Varying rates of RNA processing: Different RNAs are processed at different rates.
- Continuous or discrete cell identities (e.g. the pro-inflammatory potential of each individual T cell): Continuous phenotypes are by definitition variable in gene expression, and separating the continuous from the discrete can sometimes be difficult.
- Environmental stimuli: The local environment of the cell can influence the gene expression depending on spatial position, signaling molecules, etc.
- Temporal changes: Fundamental fluxuating cellular processes, such as cell cycle, can affect the gene expression profiles of individual cells.
Image credit: Wagner, A, et al. Revealing the vectors of cellular identity with single-cell genomics, Nat Biotechnol. 2016 (doi:https://dx.doi.org/10.1038%2Fnbt.3711)
細胞/樣品之間的技術(shù)差異
- 細胞特異性捕獲效率：不同細胞捕獲的轉(zhuǎn)錄本不同膀值，導致測序深度不同（例如，轉(zhuǎn)錄組的10-50％）误辑。
- 文庫質(zhì)量：降解的RNA沧踏、低存活力/瀕死細胞、大量自由漂浮的RNA以及細胞定量不準確會導致質(zhì)量指標降低巾钉。
- 擴增偏差：在文庫制備的擴增步驟中翘狱，并非所有轉(zhuǎn)錄本都擴增到相同水平。
- 批次效應：對于scRNA-Seq分析睛琳，批次效應是一個重要的問題盒蟆，因為看到的顯著差異表達可能只是因為批次效應引起的踏烙。
Image credit: Hicks SC, et al., bioRxiv (2015)(https://www.biorxiv.org/content/early/2015/08/25/025528)
批次不良所產(chǎn)生的問題在這篇文章中有很好的介紹：https://f1000research.com/articles/4-121/v1
如何知道你的實驗中具有批次呢师骗？
- 是否在同一天進行了所有RNA的分離？
- 是否在同一天進行了所有建庫工作讨惩？
- 是否由同一個人對所有樣品進行RNA分離/文庫制備辟癌？
- 是否對所有樣品使用相同的試劑？
- 是否在同一位置進行RNA分離/文庫制備荐捻？
如果以上任一問題的答案為“否”黍少，說明你的實驗中具有批次寡夹。
有關(guān)批次的最佳做法：
- 盡可能避免以批次的方式設計實驗。
- 如果無法避免：
- Do NOT confound your experiment by batch:
- 將不同樣品組的重復樣品分成多個批次：
- 在實驗設計文件中加上批次信息厂置，這樣可以在分析過程中退還批次引起的差異菩掏。
建議：
- 在實驗開始之前與專家討論實驗設計。
- 同時從樣品中分離RNA昵济。
- 同時準備樣品庫或備用樣品組智绸，以避免批次混淆。
- 不要混淆性別访忿、年齡或批次的樣本組瞧栗。
（二）Single-cell RNA-seq data - raw data to count matrix
根據(jù)所用文庫制備方法的不同，RNA序列（也被稱為reads或tag）將從轉(zhuǎn)錄本（(10X Genomics, CEL-seq2, Drop-seq, inDrops）的3'端（或5'端）或全長轉(zhuǎn)錄本（Smart-seq）中獲得海铆。
Image credit: Papalexi E and Satija R. Single-cell RNA sequencing to explore immune cell heterogeneity, Nature Reviews Immunology 2018 (https://doi.org/10.1038/nri.2017.76)
不同測序方式的優(yōu)點：
3’（或5’）末端測序：
- 通過使用UMI進行更準確的定量迹恐，從而將生物學重復與擴增重復（PCR）區(qū)別開來；
- 測序的細胞數(shù)量更多卧斟，可以更好地鑒定細胞類型群殴边；
- 每個細胞成本更低；
- 大于10,000個細胞的結(jié)果最佳
全長測序：
- 檢測亞型水平（
isoform-level）表達差異珍语； - 鑒定等位基因特異性差異表達找都；
- 對較少數(shù)量的細胞進行更深的測序；
- 最適用于細胞數(shù)少的樣品廊酣。
我們將主要介紹3’端測序能耻，重點是基于液滴的方法(inDrops, Drop-seq, 10X Genomics)。
3’-end reads (includes all droplet-based methods)
在3’端測序中亡驰，同一轉(zhuǎn)錄本的不同reads片段僅會源自轉(zhuǎn)錄本的3’端晓猛，相同序列的可能性很高，同時在建庫過程中的PCR步驟可能導致reads的重復凡辱，因此為了區(qū)分是生物學還是技術(shù)上的重復戒职，我們使用唯一標識符（unique molecular identifiers，UMI）進行標注透乾。
- 比對到相同的轉(zhuǎn)錄本洪燥、UMI不同的reads來源于不同的分子，為生物學重復乳乌，每個read都被計數(shù)捧韵。
- UMI相同的reads來自同一分子，為技術(shù)重復汉操，計為1個read再来。
我們以下圖為例，下圖中分子ACTB的UMI均相同，因此只能記為1個molecule芒篷，而ARL1的UMI不同所以可以記為2個molecule搜变。
Image credit: modified from Macosko EZ et al. Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets, Cell 2015 (https://doi.org/10.1016/j.cell.2015.05.002)
在細胞水平進行正確定量都需要以下條件：
- Sample index: 樣本來源
- Added during library preparation - needs to be documented
- Cellular barcode: 細胞來源
- Each library preparation method has a stock of cellular barcodes used during the library preparation
- Unique molecular identifier (UMI): 轉(zhuǎn)錄本來源
- The UMI will be used to collapse PCR duplicates
- Sequencing read1: the Read1 sequence
- Sequencing read2: the Read2 sequence
例如，使用inDrops v3庫準備方法時针炉，以下內(nèi)容是reads的所有信息：
Image credit: Sarah Boswell(https://scholar.harvard.edu/saboswell), Director of the Single Cell Sequencing Core at HMS
- R1 (61 bp Read 1): sequence of the read (Red top arrow)
- R2 (8 bp Index Read 1 (i7)): cellular barcode - which cell read originated from (Purple top arrow)
- R3 (8 bp Index Read 2 (i5)): sample/library index - which sample read originated from (Red bottom arrow)
- R4 (14 bp Read 2): read 2 and remaining cellular barcode and UMI - which transcript read originated from (Purple bottom arrow)
對于不同的基于液滴的scRNA-seq方法挠他，scRNA-seq的分析工作流程相似，但是UMI篡帕、細胞ID和樣品索引的解析會有所不同绩社。例如，以下是10X序列reads的示意圖赂苗，其中index愉耙，UMI和barcode的位置不同：
Image credit: Sarah Boswell(https://scholar.harvard.edu/saboswell), Director of the Single Cell Sequencing Core at HMS
Single-cell RNA-seq workflow
scRNA-seq方法能通過測序的reads解析barcodes和UMI，它們在特定步驟里會輕微地不同拌滋，但除了方法外朴沿，大致流程都是一致的，常規(guī)工作流程如下所示：
Image credit: Luecken, MD and Theis, FJ. Current best practices in single‐cell RNA‐seq analysis: a tutorial, Mol Syst Biol 2019 (doi: https://doi.org/10.15252/msb.20188746)
工作流程的步驟是：
- 生成count矩陣（
method-specific steps）：reads格式化赌渣，對樣本進行多路分解（demultiplexing，即通過barcodes確定reads的來源）昌犹，比對和定量坚芜。 - 原始count的質(zhì)量控制：過濾質(zhì)量較差的細胞。
- 過濾count的聚類：基于轉(zhuǎn)錄活性的相似性對細胞進行聚類（細胞類型數(shù)=簇數(shù)）斜姥。
-
marker識別：識別每個cluster的標記基因鸿竖。 - 可選的下游步驟。
無論進行那種分析铸敏，生物學重復都是必要的缚忧！
Generation of count matrix
我們聚焦于基于液滴型的3’端測序（比如inDrops、10X Genomics和Drop-seq）杈笔，將原始測序數(shù)據(jù)轉(zhuǎn)換為count矩陣闪水。
測序工具將以BCL或FASTQ格式輸出原始測序數(shù)據(jù)，或生成count矩陣蒙具。如果reads是BCL格式球榆，我們將需要轉(zhuǎn)換為FASTQ格式。有一個有用的命令行工具bcl2fastq禁筏，可以輕松執(zhí)行此轉(zhuǎn)換持钉。
NOTE: We do not demultiplex at this step in the workflow. You may have sequenced 6 samples, but the reads for all samples may be present all in the same BCL or FASTQ file.
對于許多scRNA-seq方法，從原始測序數(shù)據(jù)中生成count矩陣都將經(jīng)歷相似的步驟融师。
umis（https://github.com/vals/umis）和zUMIs（https://github.com/vals/umis）是命令行工具右钾，可用于估計測轉(zhuǎn)錄本3'端的scRNA-seq數(shù)據(jù)的表達蚁吝。此過程中的步驟包括：
- 格式化reads并過濾嘈雜的細胞
barcodes旱爆； -
Demultiplexing the samples（通過barcodes確定reads的來源）舀射； - 比對/偽比對到轉(zhuǎn)錄；
- 折疊UMI和定量reads怀伦。
當然脆烟，如果使用10X Genomics建庫方法，Cell Ranger pipeline(https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger)將負責執(zhí)行以上的所有步驟房待。
格式化reads并過濾嘈雜的細胞barcodes：
FASTQ文件能解析得到細胞barcodes邢羔、UMIs和樣本barcodes。對于基于液滴型的方法桑孩，一些細胞barcodes會對應的低的reads數(shù)(< 1000 reads) 拜鹤，原因是：
- encapsulation of free floating RNA from dying cells
- simple cells (RBCs, etc.) expressing few genes
- cells that failed for some reason
在比對reads之前，需要從序列數(shù)據(jù)中過濾掉多余的條形碼流椒。為了進行這種過濾敏簿，提取并保存每個細胞的“細胞條形碼”和“分子條形碼”钥飞。例如募狂，如果使用“umis”工具，則信息將以以下格式添加到每次reads的標題行中：
@HWI-ST808:130:H0B8YADXX:1:1101:2088:2222:CELL_GGTCCA:UMI_CCCT
AGGAAGATGGAGGAGAGAAGGCGGTGAAAGAGACCTGTAAAAAGCCACCGN
+
@@@DDBD>=AFCF+<CAFHDECII:DGGGHGIGGIIIEHGIIIGIIDHII#
建庫中使用的細胞條形碼應該是已知的撞牢，未知的條形碼會被丟棄绣硝，同時對于已知的細胞條形碼可以有一定量的錯配蜻势。
Demultiplexing the samples：
如果測序多于一個樣品執(zhí)行此步驟，這是一步不由“umis”工具處理鹉胖，而由“zUMIs”完成的步驟握玛，這步會解析reads以確定與每個與細胞相關(guān)的樣本條形碼。
比對/偽比對到轉(zhuǎn)錄：
通過傳統(tǒng)（STAR）或輕量型（Kallisto/RapMap）方法甫菠，將reads比對回基因败许。
折疊UMI和定量reads：
使用Kallisto或featureCounts之類的工具僅對唯一的UMI進行量化，得到
Image credit: extracted from Lafzi et al. Tutorial: guidelines for the experimental design of single-cell RNA sequencing studies, Nature Protocols 2018 (https://doi.org/10.1038/s41596-018-0073-y)
矩陣中的每個值代表源自相應基因的細胞中的reads數(shù)淑蔚。
（三）Single-cell RNA-seq: Quality control set-up
在生成count矩陣后市殷，我們需要對其進行QC分析，并將其導入R進行后續(xù)步驟：
探索示例集
該數(shù)據(jù)集來自Kang et al, 2017(https://www.nature.com/articles/nbt.4042)的文章刹衫，是八名狼瘡患者的PBMC（Peripheral Blood Mononuclear Cells）數(shù)據(jù)醋寝，將其分為對照組和干擾素beta處理（刺激）組。
Image credit: Kang et al, 2017(https://www.nature.com/articles/nbt.4042)
Raw data
研究團隊發(fā)現(xiàn)GEO (GSE96583：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE96583)提供的矩陣缺少線粒體的reads带迟，因此從SRA (SRP102802：https://www-ncbi-nlm-nih-gov.ezp-prod1.hul.harvard.edu/sra?term=SRP102802)下載BAM文件音羞，然后轉(zhuǎn)化為FASTQ文件，并使用Cell Ranger(https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger)獲得count矩陣仓犬。
NOTE: The counts for this dataset is also freely available from 10X Genomics and is used as part of the
Seurat tutorial(https://satijalab.org/seurat/v3.0/immune_alignment.html).（我說我咋這么眼熟嗅绰。。。）
Metadata
相關(guān)信息也叫metadata窘面。關(guān)于以上數(shù)據(jù)的metadata如下：
- 使用10X Genomics版本2化學試劑制備文庫翠语；
- 樣品在Illumina NextSeq 500上測序；
- 將來自八名狼瘡患者的PBMC樣本分別分成兩等份：
- 一份用100 U/mL重組IFN-β激活PBMC财边，處理時間為6小時肌括。
- 另一份樣品未經(jīng)處理。
- 6小時后酣难，將不同條件的8個樣品一起收集起來（受激細胞和對照細胞）谍夭。
- 分別鑒定了
12138和12167個細胞（去除doublets后）作為對照樣品和刺激后的合并樣品。 - 由于樣品是PBMC憨募，因此我們期望有免疫細胞紧索，例如：
- B cells
- T cells
- NK cells
- monocytes
- macrophages
- possibly megakaryocytes
It is recommended that you have some expectation regarding the cell types you expect to see in a dataset prior to performing the QC. This will inform you if you have any cell types with low complexity (lots of transcripts from a few genes) or cells with higher levels of mitochondrial expression. This will enable us to account for these biological factors during the analysis workflow.
上述細胞類型都不預估有低復雜度或高線粒體含量。
設置R環(huán)境
為了更好的管理數(shù)據(jù)菜谣，使得整個項目具有組織性齐板，先創(chuàng)建一個項目叫“single_cell_rnaseq”，然后構(gòu)建以下目錄：
single_cell_rnaseq/
├── data
├── results
└── figures
下載數(shù)據(jù)
- Control sample(https://www.dropbox.com/sh/73drh0ipmzfcrb3/AADMlKXCr5QGoaQN13-GeKSa?dl=1)
- Stimulated sample(https://www.dropbox.com/sh/cii4j356moc08w5/AAC2c3jfvh2hHWPmEaVsZKRva?dl=1)
解壓數(shù)據(jù)后新建一個R腳本葛菇，命名為“quality_control.R”甘磨，并且保證整個的工作目錄看起來像這樣：
加載R包
library(SingleCellExperiment)
library(Seurat)
library(tidyverse)
library(Matrix)
library(scales)
library(cowplot)
library(RCurl)
加載count矩陣
一共有3個文件，分別為：分別為cell IDs眯停、gene IDs和matrix of counts（每個細胞的每個基因）济舆。
barcodes.tsv
這是一個文本文件，其中包含該樣品存在的所有細胞條形碼莺债。條形碼以矩陣文件中顯示的數(shù)據(jù)順序列出（這些是列名）滋觉。
features.tsv
這是一個文本文件，其中包含定量基因的標識符齐邦。標識符的來源可能會有所不同椎侠，具體取決于定量過程中使用的參考（即Ensembl，NCBI措拇，UCSC）我纪，大多數(shù)情況下這些都是官方gene symbol。這些基因的順序與矩陣文件中各行的順序相對應（這些是行名）丐吓。
-
matrix.mtx
注意該矩陣中有許多zero values浅悉。
有兩種方法可以進行矩陣讀入，分別為readMM()和Read10X()券犁。
- readMM(): This function is from the Matrix package and will turn our standard matrix into a sparse matrix. The features.tsv file and barcodes.tsv must first be individually loaded into R and then they are combined. For specific code and instructions on how to do this please see our additional material（https://github.com/hbctraining/scRNA-seq/blob/master/lessons/readMM_loadData.md）.
- Read10X(): This function is from the Seurat package and will use the Cell Ranger output directory as input. In this way individual files do not need to be loaded in, instead the function will load and combine them into a sparse matrix for you. We will be using this function to load in our data!
Reading in a single sample (read10X())
我們通常使用Read10X()术健。原因就要從頭說起了，一般在軟件Cell Ranger處理10X數(shù)據(jù)后會生成一個outs目錄粘衬，在outs文件夾中有這樣幾個文件：
web_summary.html: 包括許多QC指標荞估，預估細胞數(shù)咳促，比對率等；
BAM alignment files：用于可視化比對的reads和重新創(chuàng)建FASTQ文件勘伺；
filtered_feature_bc_matrix:經(jīng)過Cell Ranger過濾后構(gòu)建矩陣所需要的所有文件跪腹；
raw_feature_bc_matrix: 未過濾的可以用于構(gòu)建矩陣的文件；
如果你有一個單個樣本的話娇昙，可以直接構(gòu)建a Seurat object(https://github.com/satijalab/seurat/wiki/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)
NOTE: The min.features argument specifies the minimum number of genes that need to be detected per cell. This argument will filter out poor quality cells that likely just have random barcodes encapsulated without any cell present. Usually, cells with less than 100 genes detected are not considered for analysis.
在使用Read10X（）函數(shù)讀取數(shù)據(jù)時尺迂，Seurat會為每個細胞自動創(chuàng)建metadata笤妙，存儲在meta.dataslot冒掌。
The Seurat object is a custom list-like object that has well-defined spaces to store specific information/data. You can find more information about the slots in the Seurat object at this link:https://github.com/satijalab/seurat/wiki/Seurat
# Explore the metadata
head(ctrl@meta.data)
metadata 的每一列代表什么呢？
-
orig.ident:細胞聚類的cluster蹲盘，含有樣本的身份信息（已知的情況下）股毫； -
nCount_RNA:每個細胞的UMI數(shù)目； -
nFeature_RNA:檢測到的每個細胞中的genes數(shù)目召衔。
使用for循環(huán)對多個數(shù)據(jù)進行讀入
# Create each individual Seurat object for every 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)
}
于是分別形成了ctrl_raw_feature_bc_matrix和stim_raw_feature_bc_matrix兩個seurat對象铃诬，使用merge函數(shù)對不同數(shù)據(jù)進行合并：
# 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可以用于不同的樣本，所以我們使用add.cell.id參數(shù)為每個細胞ID添加特定于樣本的前綴苍凛。如果我們查看合并對象的metadata趣席，我們應該能夠在行名中看到前綴：
# Check that the merged object has the appropriate sample-specific prefixes
head(merged_seurat@meta.data)
tail(merged_seurat@meta.data)
（四）Single-cell RNA-seq: Quality control
在QC步驟中，我們的目標和挑戰(zhàn)主要包括：
目標：
- 篩選數(shù)據(jù)以僅包括高質(zhì)量的真實細胞醇蝴，以便在對細胞進行聚類時宣肚，更容易識別不同的細胞類型.
- 識別任何失敗的樣本并嘗試挽救數(shù)據(jù)或從分析中刪除，此外還試圖了解樣本失敗的原因.
挑戰(zhàn)：
- 從不太復雜的細胞中劃定質(zhì)量較差的細胞悠栓。
- 選擇合適的閾值進行過濾霉涨，保留高質(zhì)量的細胞而不去除生物學相關(guān)的細胞類型。
建議：在執(zhí)行質(zhì)量控制之前惭适，需要對細胞類型進行一個預估笙瑟。例如，您是否需要樣品中的復雜性較低的細胞或線粒體表達水平較高的細胞癞志？如果是往枷，那么在評估我們的數(shù)據(jù)質(zhì)量時，我們需要考慮這種生物學特征凄杯。
生成質(zhì)控指標：
除了我們上面提到過的orig.ident师溅、nCount_RNA和nFeature_RNA外，我們還可以計算其他QC指標如number of genes detected per UMI（能反映數(shù)據(jù)的復雜度盾舌，UMI數(shù)越大數(shù)據(jù)復雜度越高）和mitochondrial ratio（可以得知來源于線粒體基因的reads）墓臭。
#將每個細胞的UMI數(shù)量進行l(wèi)og10轉(zhuǎn)換并加入到metadata
merged_seurat$log10GenesPerUMI <- log10(merged_seurat$nFeature_RNA)/log10(merged_seurat$nCount_RNA)
#計算線粒體相關(guān)基因比例，注意妖谴！("^MT-")只限于人哦窿锉。酌摇。。merged_seurat$mitoRatio <- PercentageFeatureSet(object = merged_seurat, pattern = "^MT-")
merged_seurat$mitoRatio <- merged_seurat@meta.data$mitoRatio / 100
當然嗡载，我們也可以使用 $ 添加其他質(zhì)控指標到metadata中（R變量索引 - 什么時候使用 @或$）窑多，這樣不會影響merged_seurat對象。我們可以直接從Seurat對象中提取meta.dataslot 創(chuàng)建metadata：
# Create metadata dataframe
metadata <- merged_seurat@meta.data
你應該看到每個細胞ID都有一個ctrl_或stim_前綴洼滚，這是我們合并Seurat對象時指定的埂息，它們應與orig.ident中列出的樣本匹配。我們可以先添加帶有細胞ID的列遥巴，并更改當前列名稱以使其更加直觀：
# Add cell IDs to metadata
metadata$cells <- rownames(metadata)
# Rename columns
metadata <- metadata %>%
dplyr::rename(seq_folder = orig.ident,
nUMI = nCount_RNA,
nGene = nFeature_RNA)
基于細胞前綴得到每個細胞的樣本名千康。給sample這一列進行命名：
# Create sample column
metadata$sample <- NA
metadata$sample[which(str_detect(metadata$cells, "^ctrl_"))] <- "ctrl"
metadata$sample[which(str_detect(metadata$cells, "^stim_"))] <- "stim"
最終的metadata，包含了用來質(zhì)控的指標:
# 將新的metadata重新添加回seurat
merged_seurat@meta.data <- metadata
#保存
save(merged_seurat, file="data/merged_filtered_seurat.RData")
評估質(zhì)量指標
我們將對多個質(zhì)量指標進行控制并去除不合格的細胞铲掐，這些指標包括：
- 細胞計數(shù)（
Cell counts） - 每個細胞的UMI（
UMI counts per cell） - 每個細胞觀察到的基因（
Genes detected per cell） UMIs vs. genes detected- 線粒體counts比例（
Mitochondrial counts ratio） Novelty
什么是doublets拾弃？簡單的說就是兩個細胞混在一起，可能發(fā)生在細胞捕獲過程中摆霉，并且可能會誤導認為是兩種細胞類型的過渡態(tài)（transitory states）豪椿，所以應該去除（單細胞預測Doublets軟件包匯總-過渡態(tài)細胞是真的嗎？）携栋。
我們?yōu)槭裁床粰z查doublets呢搭盾？許多的工作流程都是通過設置UMI或genes的最大閾值進行的，其原理為大量的reads或基因表明存在著多個細胞婉支，盡管此理由似乎很直觀鸯隅，但并不準確。Scrublet(https://github.com/AllonKleinLab/scrublet)是檢測doublets的重要工具磅摹，但是我們還沒有對它進行基準測試滋迈。我們不建議使用閾值進行篩選，當我們確定了每個cluster的marker后户誓，我們建議探索這些標記饼灿，以確定這些marker是否適用于一種以上的細胞類型。
細胞計數(shù)
細胞計數(shù)由檢測到的barcode的數(shù)量確定帝美。對于此實驗碍彭，預計將有12,000-13,000個細胞。inDrops的細胞捕獲效率略高（70-80％）悼潭，而10X則為50-60％庇忌。注意：如果用于文庫制備的細胞濃度不準確，捕獲效率可能會低得多舰褪。我們不應通FACS或Bioanalyzer確定細胞濃度（這些工具無法準確確定濃度）皆疹，而應使用血細胞計數(shù)器或自動細胞計數(shù)器來計算細胞濃度。
同時在inDrops和10X中占拍，也會發(fā)生得到的細胞數(shù)目（即細胞barcodes數(shù)）大于細胞數(shù)的情況略就。
# Visualize the number of cell counts per sample
metadata %>%
ggplot(aes(x=sample, fill=sample)) +
geom_bar() +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
theme(plot.title = element_text(hjust=0.5, face="bold")) +
ggtitle("NCells")
從上圖可以看出每個sample的細胞數(shù)均在15000以上捎迫。
每個細胞的UMI
每個細胞的UMI計數(shù)通常應高于500，500是預期的下限表牢。如果UMI在500-1000計數(shù)之間窄绒，雖然可以使用，但可能應該對細胞進行更深的測序崔兴。
# Visualize the number UMIs/transcripts per cell
metadata %>%
ggplot(aes(color=sample, x=nUMI, fill= sample)) +
geom_density(alpha = 0.2) +
scale_x_log10() +
theme_classic() +
ylab("Cell density") +
geom_vline(xintercept = 500)
我們可以看到大多數(shù)細胞的UMI數(shù)均在1000以上彰导。
每個細胞觀察到的基因
正常情況應該和UMI相同，也是有一個大大的峰敲茄，如果我們在主峰的右邊也看到了一個小峰位谋，出現(xiàn)的原因也是多種多樣的，可能是由于捕獲單個細胞失敗折汞，或者具有某些生物學意義倔幼，需要仔細評估盖腿。
對于高質(zhì)量數(shù)據(jù)爽待，比例直方圖應包含一個大峰，該峰代表封裝的細胞翩腐。如果我們在主要峰的右側(cè)看到一個小峰（我們使用的數(shù)據(jù)中不存在這種情況）鸟款，或者細胞的雙峰分布，則可能表明有兩件事茂卦。一個可能是某些原因?qū)е乱唤M細胞失敗何什。一個可能是存在生物學上不同類型的細胞（例如，靜止細胞群體等龙，不太復雜的目標細胞）处渣，和/或一種類型比另一種類型小得多（即，具有高計數(shù)的細胞可能是尺寸較大的細胞） ）的情況蛛砰。
因此應該使用本課中介紹的其他指標來評估此閾值罐栈。
# Visualize the distribution of genes detected per cell via histogram
metadata %>%
ggplot(aes(color=sample, x=nGene, fill= sample)) +
geom_density(alpha = 0.2) +
theme_classic() +
scale_x_log10() +
geom_vline(xintercept = 300)
# Visualize the distribution of genes detected per cell via boxplot
metadata %>%
ggplot(aes(x=sample, y=log10(nGene), fill=sample)) +
geom_boxplot() +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
theme(plot.title = element_text(hjust=0.5, face="bold")) +
ggtitle("NCells vs NGenes")
UMIs vs. genes
通常我們會一起評估UMI的數(shù)量和每個細胞檢測到的基因數(shù)量。下圖是通過線粒體reads的分數(shù)繪制的基因數(shù)和UMI數(shù)之間的關(guān)系泥畅，顯示出線粒體reads分數(shù)在帶有低基因數(shù)的低計數(shù)細胞中很高荠诬。
為什么我們需要去除線粒體基因占比較高的細胞呢，主要是由于破損細胞的胞質(zhì)中的mRNA會由于細胞破損游離出來位仁，而線粒體中的被保留了下來柑贞，線粒體基因比例高。這樣的細胞可以通過count數(shù)和基因數(shù)聯(lián)合去除聂抢。
質(zhì)量較差的細胞的基因和UMI可能很低钧嘶，對應于該圖左下象限中的數(shù)據(jù)點。好的細胞通常會同時顯示每個細胞更多的基因和更高數(shù)量的UMI琳疏。
通過此圖有决，我們還可以評估線的斜率以及該圖右下象限中數(shù)據(jù)點的分布摄欲。具有大量的UMI，但只有少數(shù)基因的細胞可能是垂死的細胞疮薇，也可能是低復雜度的細胞類型（即紅細胞）胸墙。
# Visualize the correlation between genes detected and number of UMIs and determine whether strong presence of cells with low numbers of genes/UMIs
metadata %>%
ggplot(aes(x=nUMI, y=nGene, color=mitoRatio)) +
geom_point() +
scale_colour_gradient(low = "gray90", high = "black") +
stat_smooth(method=lm) +
scale_x_log10() +
scale_y_log10() +
theme_classic() +
geom_vline(xintercept = 500) +
geom_hline(yintercept = 250) +
facet_wrap(~sample)
線粒體counts比例
我們可以識別死細胞或垂死細胞是否存在大量線粒體污染。除非必要按咒，我們將線粒體比率超過0.2的細胞定義為線粒體計數(shù)質(zhì)量差的樣品迟隅。
# Visualize the distribution of mitochondrial gene expression detected per cell
metadata %>%
ggplot(aes(color=sample, x=mitoRatio, fill=sample)) +
geom_density(alpha = 0.2) +
scale_x_log10() +
theme_classic() +
geom_vline(xintercept = 0.2)
復雜度
其實進行少量測序的每個細胞具有高度的復雜性，這主要是由于我們沒有對每個樣品進行基因深度測序励七。樣品中的異常細胞可能其RNA復雜度較低智袭，我們可以通過該指標檢測篩選低復雜度（如紅細胞）的污染。一般情況下掠抬，我們把評分（novelty score）設置為0.8吼野。
# Visualize the overall complexity of the gene expression by visualizing the genes detected per UMI
metadata %>%
ggplot(aes(x=log10GenesPerUMI, color = sample, fill=sample)) +
geom_density(alpha = 0.2) +
theme_classic() +
geom_vline(xintercept = 0.8)
NOTE: Reads per cell is another metric that can be useful to explore; however, the workflow used would need to save this information to assess. Generally, with this metric you hope to see all of the samples with peaks in relatively the same location between 10,000 and 100,000 reads per cell.
過濾
線粒體數(shù)相對高的細胞也可能參與了呼吸過程，是某研究中正需要的两波，因此作者認為在設置條件時應綜合考慮瞳步，盡量寬松，防止有意義的生物學變化被刪除腰奋。
細胞層面過濾
- nUMI > 500
- nGene > 250
- log10GenesPerUMI > 0.8
- mitoRatio < 0.2
用subset()函數(shù)對Seurat對象進行設置：
# Filter out low quality reads using selected thresholds - these will change with experiment
filtered_seurat <- subset(x = merged_seurat,
subset= (nUMI >= 500) &
(nGene >= 250) &
(log10GenesPerUMI > 0.80) &
(mitoRatio < 0.20))
基因?qū)用孢^濾
刪除0表達值的基因和在少于10個細胞中表達的基因单起。
# Output a logical vector for every gene on whether the more than zero counts per cell
# Extract counts
counts <- GetAssayData(object = filtered_seurat, slot = "counts")
# Output a logical vector for every gene on whether the more than zero counts per cell
nonzero <- counts > 0
# Sums all TRUE values and returns TRUE if more than 10 TRUE values per gene
keep_genes <- Matrix::rowSums(nonzero) >= 10
# Only keeping those genes expressed in more than 10 cells
filtered_counts <- counts[keep_genes, ]
# Reassign to filtered Seurat object
filtered_seurat <- CreateSeuratObject(filtered_counts, meta.data = filtered_seurat@meta.data)
Re-assess QC metrics
過濾之后可以回顧一下過濾前后的效果，以評判過濾指標是否合適劣坊。
Saving filtered cells
# Create .RData object to load at any time
save(filtered_seurat, file="data/seurat_filtered.RData")
（五）Count Normalization and Principal Component Analysis
獲得高質(zhì)量的單細胞后嘀倒，單細胞RNA-seq（scRNA-seq）分析工作流程的下一步就是執(zhí)行聚類。聚類的目標是將不同的細胞類型分成獨特的細胞亞群局冰。為了進行聚類测蘑，我們確定了在細胞之間表達差異最大的基因。然后康二，我們使用這些基因來確定哪些相關(guān)基因集是造成細胞之間表達差異最大的原因碳胳。
數(shù)值歸一化
歸一化最重要的目的就是使表達水平在細胞之間和/或細胞內(nèi)更具有可比性。那么在歸一化中主要需要處理的因素包括：
測序深度：考慮測序深度是比較細胞之間基因表達的必要條件赠摇。在下面的示例中固逗，每個基因在細胞2中的表達似乎都增加了一倍，但這是細胞2具有兩倍測序深度的結(jié)果藕帜。
因此烫罩，要準確比較細胞之間的表達，有必要對測序深度進行歸一化洽故。
基因長度：需要基因長度來比較同一細胞內(nèi)不同基因之間的表達贝攒。比對到較長基因的reads數(shù)似乎與較高表達的較短基因具有相同的計數(shù)/表達。
- 如果進行的5’末端或3’末端的測序时甚，則不需要考慮基因長度的影響隘弊；
- 如果使用全長測序則需要考慮哈踱。
主成分分析（PCA）
PCA是一個用來強調(diào)差異和相似性的技術(shù)（一文看懂PCA主成分分析），這里推薦一個學習視頻：StatQuest's video（https://www.youtube.com/watch?v=_UVHneBUBW0）
下面是PCA的示例模擬過程梨熙，幫助理解：
如果你已經(jīng)定量了兩個樣本（或細胞）中四個基因的表達开镣，則可以繪制這些基因的表達值，其中一個樣本在x軸上表示咽扇，另一個樣本在y軸上表示邪财，如下所示：
我們可以沿代表最大變化的方向在數(shù)據(jù)上畫一條線，在此示例中為對角線质欲。數(shù)據(jù)集中的最大變異是在組成兩個端點的基因树埠。我們還看到基因在該線的上方和下方有些不同。我們可以在數(shù)據(jù)上繪制另一條線嘶伟，代表數(shù)據(jù)中變化第二大的變量怎憋。
末端附近的基因?qū)⑹亲儺愖畲蟮幕颉＿@些基因在數(shù)學上對線的方向影響最大九昧。
例如绊袋，基因C值的微小變化將極大地改變較長線的方向，而基因A或基因D的微小變化對其幾乎沒有影響耽装。
我們還可以旋轉(zhuǎn)整個圖愤炸，保證線條方向是從左到右和從上到下∑诰荆現(xiàn)在掉奄，可以將這些線視為代表變化的軸。這些軸本質(zhì)上是“主要組分”凤薛，其中PC1代表數(shù)據(jù)的最大差異姓建，PC2代表數(shù)據(jù)的第二大差異。
如果有N個細胞缤苫，以此類推速兔。。活玲。
確定PCs后涣狗，則需要對每個PC進行評分，按照以下步驟對所有樣本PC對（sample-PC pairs）計算分數(shù)：
（1）首先舒憾，根據(jù)基因?qū)γ總€PC的影響程度镀钓，為其分配“影響力”評分。對給定PC沒有任何影響的基因得分接近零镀迂，而具有更大影響力的基因得分更高丁溅。PC線末端的基因?qū)a(chǎn)生更大的影響，因此它們將獲得更大的分數(shù)探遵，但兩端的符號相反窟赏。
（2）確定影響分數(shù)后妓柜，使用以下公式計算每個樣本的分數(shù)：
Sample1 PC1 score = (read count * influence) + ... for all genes
以我們的2個樣本示例，以下是分數(shù)的計算方式：
## Sample1
PC1 score = (4 * -2) + (1 * -10) + (8 * 8) + (5 * 1) = 51
PC2 score = (4 * 0.5) + (1 * 1) + (8 * -5) + (5 * 6) = -7
## Sample2
PC1 score = (5 * -2) + (4 * -10) + (8 * 8) + (7 * 1) = 21
PC2 score = (5 * 0.5) + (4 * 1) + (8 * -5) + (7 * 6) = 8.5
（3）一旦為所有PC計算了這些分數(shù)涯穷，就可以將其繪制在簡單的散點圖上棍掐。下面是示例圖：
對于具有大量樣本或細胞的數(shù)據(jù)集，通常會繪制每個樣本/細胞的PC1和PC2分數(shù)拷况。由于這些PC解釋了數(shù)據(jù)集中最大的變化塌衰，因此更相似的樣本/細胞將在PC1和PC2聚在一起。請參見下面的示例：
Image credit: https://github.com/AshwiniRS/Medium_Notebooks/blob/master/PCA/PCA_Iris_DataSet.ipynb
對于我們的單細胞數(shù)據(jù)蝠嘉，我們最終會選擇10-100 PC去對細胞進行聚類分析最疆，而不是全部基因。
未完待續(xù)蚤告。
