經(jīng)典賞析

NGS系列文章包括NGS基礎(chǔ)客扎、轉(zhuǎn)錄組分析?（Nature重磅綜述|關(guān)于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）霎褐、單細(xì)胞測(cè)序分析?(重磅綜述：三萬(wàn)字長(zhǎng)文讀懂單細(xì)胞RNA測(cè)序分析的最佳實(shí)踐教程 （原理、代碼和評(píng)述）)亚亲、DNA甲基化分析救拉、重測(cè)序分析难审、GEO數(shù)據(jù)挖掘（典型醫(yī)學(xué)設(shè)計(jì)實(shí)驗(yàn)GEO數(shù)據(jù)分析 (step-by-step) - Limma差異分析、火山圖亿絮、功能富集）等內(nèi)容告喊。

單細(xì)胞轉(zhuǎn)錄組測(cè)序進(jìn)展迅速，伴隨而來(lái)的是許許多多的內(nèi)容與講義派昧，很多課程都老長(zhǎng)黔姜。。蒂萎。嗨像我這樣莫有耐心的小孩子秆吵，能完整看完的只有一個(gè) —— Single-cell RNA-seq analysis workshop(https://github.com/hbctraining/scRNA-seq)。我將把教程內(nèi)容盡量匯總岖是，你需要做的就是點(diǎn)擊收藏慢慢看帮毁！

（一）Why single-cell RNA-seq

在整個(gè)人體組織中，細(xì)胞類型豺撑、狀態(tài)和相互作用非常多樣。為了更好地了解這些組織存在的細(xì)胞類型黔牵，單細(xì)胞RNA-seq（scRNA-seq）提供了在單細(xì)胞水平上表達(dá)基因的信息聪轿。

單細(xì)胞轉(zhuǎn)錄組測(cè)序可用于 ：

常見(jiàn)方法包括：

Challenges of scRNA-seq analysis

在scRNA-seq之前灯抛，我們通常使用bulk RNA-seq進(jìn)行轉(zhuǎn)錄組分析。bulk RNA-seq是一種直接比較細(xì)胞的平均表達(dá)值的方法音瓷，在尋找疾病生物標(biāo)志物对嚼，或者不關(guān)心樣品中大量細(xì)胞異質(zhì)性的情況下，這可能是最佳方法绳慎。

盡管bulk RNA-seq可以探索不同條件（例如治療或疾沧菔）之間基因表達(dá)的差異，但無(wú)法充分捕獲細(xì)胞水平的差異杏愤。例如靡砌，在下面的圖像中，如果進(jìn)行bulk分析（左）珊楼，我們將無(wú)法檢測(cè)到基因A和基因B的表達(dá)之間的正確關(guān)聯(lián)通殃。但是，如果我們按細(xì)胞類型或細(xì)胞狀態(tài)正確地對(duì)細(xì)胞進(jìn)行分組厕宗，我們可以看到基因之間的正確相關(guān)性画舌。

scRNA-seq也有一定的局限性，除了制樣和建庫(kù)價(jià)格高昂外已慢，它在數(shù)據(jù)分析中也具有一定的復(fù)雜性曲聂，包括：

數(shù)據(jù)量大

scRNA-seq實(shí)驗(yàn)的數(shù)據(jù)來(lái)自捕獲的成千上萬(wàn)個(gè)細(xì)胞，也就是數(shù)十萬(wàn)條reads蛇受，需要更多的內(nèi)存和存儲(chǔ)空間句葵。

細(xì)胞的測(cè)序深度低

基于液滴的scRNA-seq方法的測(cè)序深度較淺，通常每個(gè)細(xì)胞只能檢測(cè)到轉(zhuǎn)錄組的10-50％兢仰。這導(dǎo)致細(xì)胞中許多基因的計(jì)數(shù)為零乍丈。但是，在特定的細(xì)胞中把将，基因的零計(jì)數(shù)可能意味著該基因未表達(dá)或僅表示該轉(zhuǎn)錄本未被檢測(cè)到轻专。在整個(gè)細(xì)胞中，具有較高表達(dá)水平的基因測(cè)到0值的機(jī)率較低察蹲。由于此特征请垛，在任何細(xì)胞中都不會(huì)檢測(cè)到許多基因，并且細(xì)胞之間的基因表達(dá)差異很大洽议。

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）.

跨細(xì)胞/樣品的生物變異性

我們不感興趣的某些生物變異可能導(dǎo)致細(xì)胞之間的基因表達(dá)比實(shí)際生物細(xì)胞的類型/狀態(tài)更為相似/不同宗收，并掩蓋細(xì)胞類型。這些變異（除非實(shí)驗(yàn)研究的一部分）包括（以下為原文）：

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)

細(xì)胞/樣品之間的技術(shù)差異

Image credit: Hicks SC, et al., bioRxiv (2015)(https://www.biorxiv.org/content/early/2015/08/25/025528)

批次不良所產(chǎn)生的問(wèn)題在這篇文章中有很好的介紹：https://f1000research.com/articles/4-121/v1

如何知道你的實(shí)驗(yàn)中具有批次呢汰翠？

如果以上任一問(wèn)題的答案為“否”苟穆，說(shuō)明你的實(shí)驗(yàn)中具有批次。

有關(guān)批次的最佳做法：

• 在實(shí)驗(yàn)設(shè)計(jì)文件中加上批次信息雳旅，這樣可以在分析過(guò)程中退還批次引起的差異。

建議：

（二）Single-cell RNA-seq data - raw data to count matrix

根據(jù)所用文庫(kù)制備方法的不同尚蝌，RNA序列（也被稱為reads或tag）將從轉(zhuǎn)錄本（(10X Genomics, CEL-seq2, Drop-seq, inDrops）的3'端（或5'端）或全長(zhǎng)轉(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)

不同測(cè)序方式的優(yōu)點(diǎn)：

3’（或5’）末端測(cè)序：

全長(zhǎng)測(cè)序：

我們將主要介紹3’端測(cè)序，重點(diǎn)是基于液滴的方法(inDrops, Drop-seq, 10X Genomics)热某。

3’-end reads (includes all droplet-based methods)

在3’端測(cè)序中鲜锚，同一轉(zhuǎn)錄本的不同reads片段僅會(huì)源自轉(zhuǎn)錄本的3’端突诬，相同序列的可能性很高苫拍，同時(shí)在建庫(kù)過(guò)程中的PCR步驟可能導(dǎo)致reads的重復(fù)芜繁，因此為了區(qū)分是生物學(xué)還是技術(shù)上的重復(fù)，我們使用唯一標(biāo)識(shí)符（unique molecular identifiers绒极，UMI）進(jìn)行標(biāo)注骏令。

我們以下圖為例审丘，下圖中分子ACTB的UMI均相同吏够，因此只能記為1個(gè)molecule，而ARL1的UMI不同所以可以記為2個(gè)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)

在細(xì)胞水平進(jìn)行正確定量都需要以下條件：

• Added during library preparation - needs to be documented

• Each library preparation method has a stock of cellular barcodes used during the library preparation

• The UMI will be used to collapse PCR duplicates

例如锅知，使用inDrops v3庫(kù)準(zhǔn)備方法時(shí)，以下內(nèi)容是reads的所有信息：

Image credit: Sarah Boswell(https://scholar.harvard.edu/saboswell), Director of the Single Cell Sequencing Core at HMS

對(duì)于不同的基于液滴的scRNA-seq方法脓钾，scRNA-seq的分析工作流程相似售睹，但是UMI、細(xì)胞ID和樣品索引的解析會(huì)有所不同可训。例如昌妹，以下是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方法能通過(guò)測(cè)序的reads解析barcodes和UMI，它們?cè)谔囟ú襟E里會(huì)輕微地不同钩蚊，但除了方法外川蒙，大致流程都是一致的蚜厉，常規(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)

工作流程的步驟是：

無(wú)論進(jìn)行那種分析，生物學(xué)重復(fù)都是必要的！

Generation of count matrix

我們聚焦于基于液滴型的3’端測(cè)序（比如inDrops脊另、10X Genomics和Drop-seq）导狡，將原始測(cè)序數(shù)據(jù)轉(zhuǎn)換為count矩陣。

測(cè)序工具將以BCL或FASTQ格式輸出原始測(cè)序數(shù)據(jù)偎痛，或生成count矩陣旱捧。如果reads是BCL格式，我們將需要轉(zhuǎn)換為FASTQ格式踩麦。有一個(gè)有用的命令行工具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.

對(duì)于許多scRNA-seq方法谓谦，從原始測(cè)序數(shù)據(jù)中生成count矩陣都將經(jīng)歷相似的步驟贫橙。

umis（https://github.com/vals/umis）和zUMIs（https://github.com/vals/umis）是命令行工具，可用于估計(jì)測(cè)轉(zhuǎn)錄本3'端的scRNA-seq數(shù)據(jù)的表達(dá)反粥。此過(guò)程中的步驟包括：

當(dāng)然，如果使用10X Genomics建庫(kù)方法娜膘，Cell Ranger pipeline(https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger)將負(fù)責(zé)執(zhí)行以上的所有步驟逊脯。

格式化reads并過(guò)濾嘈雜的細(xì)胞barcodes：

FASTQ文件能解析得到細(xì)胞barcodes、UMIs和樣本barcodes竣贪。對(duì)于基于液滴型的方法军洼，一些細(xì)胞barcodes會(huì)對(duì)應(yīng)的低的reads數(shù)(< 1000 reads) ，原因是：

在比對(duì)reads之前演怎，需要從序列數(shù)據(jù)中過(guò)濾掉多余的條形碼匕争。為了進(jìn)行這種過(guò)濾，提取并保存每個(gè)細(xì)胞的“細(xì)胞條形碼”和“分子條形碼”爷耀。例如甘桑，如果使用“umis”工具，則信息將以以下格式添加到每次reads的標(biāo)題行中：

@HWI-ST808:130:H0B8YADXX:1:1101:2088:2222:CELL_GGTCCA:UMI_CCCT
AGGAAGATGGAGGAGAGAAGGCGGTGAAAGAGACCTGTAAAAAGCCACCGN
+@@@DDBD>=AFCF+<CAFHDECII:DGGGHGIGGIIIEHGIIIGIIDHII#

建庫(kù)中使用的細(xì)胞條形碼應(yīng)該是已知的歹叮，未知的條形碼會(huì)被丟棄跑杭，同時(shí)對(duì)于已知的細(xì)胞條形碼可以有一定量的錯(cuò)配。

Demultiplexing the samples：

如果測(cè)序多于一個(gè)樣品執(zhí)行此步驟咆耿，這是一步不由“umis”工具處理德谅，而由“zUMIs”完成的步驟，這步會(huì)解析reads以確定與每個(gè)與細(xì)胞相關(guān)的樣本條形碼萨螺。

比對(duì)/偽比對(duì)到轉(zhuǎn)錄：

通過(guò)傳統(tǒng)（STAR）或輕量型（Kallisto/RapMap）方法窄做，將reads比對(duì)回基因愧驱。

折疊UMI和定量reads：

使用Kallisto或featureCounts之類的工具僅對(duì)唯一的UMI進(jìn)行量化，得到

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)

矩陣中的每個(gè)值代表源自相應(yīng)基因的細(xì)胞中的reads數(shù)椭盏。

（三）Single-cell RNA-seq: Quality control set-up

在生成count矩陣后组砚，我們需要對(duì)其進(jìn)行QC分析，并將其導(dǎo)入R進(jìn)行后續(xù)步驟：

探索示例集

該數(shù)據(jù)集來(lái)自Kang et al, 2017(https://www.nature.com/articles/nbt.4042)的文章庸汗，是八名狼瘡患者的PBMC（Peripheral Blood Mononuclear Cells）數(shù)據(jù)惫确，將其分為對(duì)照組和干擾素beta處理（刺激）組。

Image credit: Kang et al, 2017(https://www.nature.com/articles/nbt.4042)

Raw data

研究團(tuán)隊(duì)發(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).（我說(shuō)我咋這么眼熟。揍鸟。兄裂。）

Metadata

相關(guān)信息也叫metadata。關(guān)于以上數(shù)據(jù)的metadata如下：

• 一份用100 U/mL重組IFN-β激活PBMC，處理時(shí)間為6小時(shí)腥泥。

• 另一份樣品未經(jīng)處理匾南。

• 6小時(shí)后，將不同條件的8個(gè)樣品一起收集起來(lái)（受激細(xì)胞和對(duì)照細(xì)胞）蛔外。

• 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.

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

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

為了更好的管理數(shù)據(jù)，使得整個(gè)項(xiàng)目具有組織性矛纹，先創(chuàng)建一個(gè)項(xiàng)目叫“single_cell_rnaseq”臂聋，然后構(gòu)建以下目錄：

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

下載數(shù)據(jù)

解壓數(shù)據(jù)后新建一個(gè)R腳本，命名為“quality_control.R”或南，并且保證整個(gè)的工作目錄看起來(lái)像這樣：

加載R包

library(SingleCellExperiment)
library(Seurat)
library(tidyverse)
library(Matrix)
library(scales)
library(cowplot)library(RCurl)

加載count矩陣

一共有3個(gè)文件孩等，分別為：分別為cell IDs、gene IDs和matrix of counts（每個(gè)細(xì)胞的每個(gè)基因）迎献。

有兩種方法可以進(jìn)行矩陣讀入瞎访，分別為readMM()和Read10X()。

1. 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）.

2. 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()吁恍。原因就要從頭說(shuō)起了扒秸，一般在軟件Cell Ranger處理10X數(shù)據(jù)后會(huì)生成一個(gè)outs目錄播演，在outs文件夾中有這樣幾個(gè)文件：

web_summary.html: 包括許多QC指標(biāo)，預(yù)估細(xì)胞數(shù)伴奥，比對(duì)率等写烤；
BAM alignment files：用于可視化比對(duì)的reads和重新創(chuàng)建FASTQ文件；
filtered_feature_bc_matrix:經(jīng)過(guò)Cell Ranger過(guò)濾后構(gòu)建矩陣所需要的所有文件拾徙；
raw_feature_bc_matrix: 未過(guò)濾的可以用于構(gòu)建矩陣的文件洲炊；

如果你有一個(gè)單個(gè)樣本的話，可以直接構(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ù)時(shí)尼啡，Seurat會(huì)為每個(gè)細(xì)胞自動(dòng)創(chuàng)建metadata暂衡，存儲(chǔ)在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 metadatahead(ctrl@meta.data)

metadata 的每一列代表什么呢崖瞭？

使用for循環(huán)對(duì)多個(gè)數(shù)據(jù)進(jìn)行讀入

# 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兩個(gè)seurat對(duì)象，使用merge函數(shù)對(duì)不同數(shù)據(jù)進(jì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"))

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

# 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步驟中驯杜，我們的目標(biāo)和挑戰(zhàn)主要包括：

目標(biāo)：

挑戰(zhàn)：

建議：在執(zhí)行質(zhì)量控制之前核畴，需要對(duì)細(xì)胞類型進(jìn)行一個(gè)預(yù)估。例如冲九，您是否需要樣品中的復(fù)雜性較低的細(xì)胞或線粒體表達(dá)水平較高的細(xì)胞谤草？如果是，那么在評(píng)估我們的數(shù)據(jù)質(zhì)量時(shí)莺奸，我們需要考慮這種生物學(xué)特征丑孩。

生成質(zhì)控指標(biāo)：
除了我們上面提到過(guò)的orig.ident、nCount_RNA和nFeature_RNA外灭贷，我們還可以計(jì)算其他QC指標(biāo)如number of genes detected per UMI（能反映數(shù)據(jù)的復(fù)雜度温学，UMI數(shù)越大數(shù)據(jù)復(fù)雜度越高）和mitochondrial ratio（可以得知來(lái)源于線粒體基因的reads）。

#將每個(gè)細(xì)胞的UMI數(shù)量進(jìn)行l(wèi)og10轉(zhuǎn)換并加入到metadatamerged_seurat$log10GenesPerUMI <- log10(merged_seurat$nFeature_RNA)/log10(merged_seurat$nCount_RNA)

#計(jì)算線粒體相關(guān)基因比例甚疟，注意仗岖！("^MT-")只限于人哦逃延。。轧拄。merged_seurat$mitoRatio <- PercentageFeatureSet(object = merged_seurat, pattern = "^MT-")merged_seurat$mitoRatio <- merged_seurat@meta.data$mitoRatio / 100

當(dāng)然揽祥，我們也可以使用 $ 添加其他質(zhì)控指標(biāo)到metadata中（R變量索引 - 什么時(shí)候使用 @或$），這樣不會(huì)影響merged_seurat對(duì)象檩电。我們可以直接從Seurat對(duì)象中提取meta.dataslot 創(chuàng)建metadata：

# Create metadata dataframemetadata <- merged_seurat@meta.data

你應(yīng)該看到每個(gè)細(xì)胞ID都有一個(gè)ctrl_或stim_前綴拄丰，這是我們合并Seurat對(duì)象時(shí)指定的，它們應(yīng)與orig.ident中列出的樣本匹配俐末。我們可以先添加帶有細(xì)胞ID的列料按，并更改當(dāng)前列名稱以使其更加直觀：

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

基于細(xì)胞前綴得到每個(gè)細(xì)胞的樣本名。給sample這一列進(jìn)行命名：

# 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鹅搪，包含了用來(lái)質(zhì)控的指標(biāo):

# 將新的metadata重新添加回seurat
merged_seurat@meta.data <- metadata

#保存save(merged_seurat, file="data/merged_filtered_seurat.RData")

評(píng)估質(zhì)量指標(biāo)

我們將對(duì)多個(gè)質(zhì)量指標(biāo)進(jìn)行控制并去除不合格的細(xì)胞站绪，這些指標(biāo)包括：

什么是doublets？簡(jiǎn)單的說(shuō)就是兩個(gè)細(xì)胞混在一起丽柿，可能發(fā)生在細(xì)胞捕獲過(guò)程中，并且可能會(huì)誤導(dǎo)認(rèn)為是兩種細(xì)胞類型的過(guò)渡態(tài)（transitory states）魂挂，所以應(yīng)該去除（單細(xì)胞預(yù)測(cè)Doublets軟件包匯總-過(guò)渡態(tài)細(xì)胞是真的嗎甫题？）。

我們?yōu)槭裁床粰z查doublets呢涂召？許多的工作流程都是通過(guò)設(shè)置UMI或genes的最大閾值進(jìn)行的坠非，其原理為大量的reads或基因表明存在著多個(gè)細(xì)胞，盡管此理由似乎很直觀果正，但并不準(zhǔn)確炎码。Scrublet(https://github.com/AllonKleinLab/scrublet)是檢測(cè)doublets的重要工具，但是我們還沒(méi)有對(duì)它進(jìn)行基準(zhǔn)測(cè)試秋泳。我們不建議使用閾值進(jìn)行篩選潦闲，當(dāng)我們確定了每個(gè)cluster的marker后，我們建議探索這些標(biāo)記迫皱，以確定這些marker是否適用于一種以上的細(xì)胞類型歉闰。

細(xì)胞計(jì)數(shù)

細(xì)胞計(jì)數(shù)由檢測(cè)到的barcode的數(shù)量確定。對(duì)于此實(shí)驗(yàn)卓起，預(yù)計(jì)將有12,000-13,000個(gè)細(xì)胞和敬。inDrops的細(xì)胞捕獲效率略高（70-80％），而10X則為50-60％戏阅。注意：如果用于文庫(kù)制備的細(xì)胞濃度不準(zhǔn)確昼弟，捕獲效率可能會(huì)低得多。我們不應(yīng)通FACS或Bioanalyzer確定細(xì)胞濃度（這些工具無(wú)法準(zhǔn)確確定濃度）奕筐，而應(yīng)使用血細(xì)胞計(jì)數(shù)器或自動(dòng)細(xì)胞計(jì)數(shù)器來(lái)計(jì)算細(xì)胞濃度舱痘。

同時(shí)在inDrops和10X中变骡，也會(huì)發(fā)生得到的細(xì)胞數(shù)目（即細(xì)胞barcodes數(shù)）大于細(xì)胞數(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")

從上圖可以看出每個(gè)sample的細(xì)胞數(shù)均在15000以上衰粹。

每個(gè)細(xì)胞的UMI

每個(gè)細(xì)胞的UMI計(jì)數(shù)通常應(yīng)高于500锣光，500是預(yù)期的下限。如果UMI在500-1000計(jì)數(shù)之間铝耻，雖然可以使用誊爹，但可能應(yīng)該對(duì)細(xì)胞進(jìn)行更深的測(cè)序。

# 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ù)細(xì)胞的UMI數(shù)均在1000以上瓢捉。

每個(gè)細(xì)胞觀察到的基因

正常情況應(yīng)該和UMI相同频丘，也是有一個(gè)大大的峰，如果我們?cè)谥鞣宓挠疫呉部吹搅艘粋€(gè)小峰泡态，出現(xiàn)的原因也是多種多樣的搂漠，可能是由于捕獲單個(gè)細(xì)胞失敗，或者具有某些生物學(xué)意義某弦，需要仔細(xì)評(píng)估桐汤。

對(duì)于高質(zhì)量數(shù)據(jù)，比例直方圖應(yīng)包含一個(gè)大峰靶壮，該峰代表封裝的細(xì)胞怔毛。如果我們?cè)谥饕宓挠覀?cè)看到一個(gè)小峰（我們使用的數(shù)據(jù)中不存在這種情況），或者細(xì)胞的雙峰分布腾降，則可能表明有兩件事拣度。一個(gè)可能是某些原因?qū)е乱唤M細(xì)胞失敗。一個(gè)可能是存在生物學(xué)上不同類型的細(xì)胞（例如螃壤，靜止細(xì)胞群體抗果，不太復(fù)雜的目標(biāo)細(xì)胞），和/或一種類型比另一種類型小得多（即奸晴，具有高計(jì)數(shù)的細(xì)胞可能是尺寸較大的細(xì)胞） ）的情況冤馏。

因此應(yīng)該使用本課中介紹的其他指標(biāo)來(lái)評(píng)估此閾值。

# 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

通常我們會(huì)一起評(píng)估UMI的數(shù)量和每個(gè)細(xì)胞檢測(cè)到的基因數(shù)量蚁滋。下圖是通過(guò)線粒體reads的分?jǐn)?shù)繪制的基因數(shù)和UMI數(shù)之間的關(guān)系宿接，顯示出線粒體reads分?jǐn)?shù)在帶有低基因數(shù)的低計(jì)數(shù)細(xì)胞中很高。

為什么我們需要去除線粒體基因占比較高的細(xì)胞呢辕录，主要是由于破損細(xì)胞的胞質(zhì)中的mRNA會(huì)由于細(xì)胞破損游離出來(lái)睦霎，而線粒體中的被保留了下來(lái)，線粒體基因比例高走诞。這樣的細(xì)胞可以通過(guò)count數(shù)和基因數(shù)聯(lián)合去除副女。

質(zhì)量較差的細(xì)胞的基因和UMI可能很低，對(duì)應(yīng)于該圖左下象限中的數(shù)據(jù)點(diǎn)蚣旱。好的細(xì)胞通常會(huì)同時(shí)顯示每個(gè)細(xì)胞更多的基因和更高數(shù)量的UMI碑幅。

通過(guò)此圖戴陡，我們還可以評(píng)估線的斜率以及該圖右下象限中數(shù)據(jù)點(diǎn)的分布。具有大量的UMI沟涨，但只有少數(shù)基因的細(xì)胞可能是垂死的細(xì)胞恤批，也可能是低復(fù)雜度的細(xì)胞類型（即紅細(xì)胞）。

# 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比例

我們可以識(shí)別死細(xì)胞或垂死細(xì)胞是否存在大量線粒體污染裹赴。除非必要喜庞，我們將線粒體比率超過(guò)0.2的細(xì)胞定義為線粒體計(jì)數(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)

復(fù)雜度

其實(shí)進(jìn)行少量測(cè)序的每個(gè)細(xì)胞具有高度的復(fù)雜性棋返，這主要是由于我們沒(méi)有對(duì)每個(gè)樣品進(jìn)行基因深度測(cè)序延都。樣品中的異常細(xì)胞可能其RNA復(fù)雜度較低，我們可以通過(guò)該指標(biāo)檢測(cè)篩選低復(fù)雜度（如紅細(xì)胞）的污染睛竣。一般情況下晰房，我們把評(píng)分（novelty score）設(shè)置為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.

過(guò)濾

線粒體數(shù)相對(duì)高的細(xì)胞也可能參與了呼吸過(guò)程射沟，是某研究中正需要的殊者，因此作者認(rèn)為在設(shè)置條件時(shí)應(yīng)綜合考慮，盡量寬松验夯，防止有意義的生物學(xué)變化被刪除幽污。

細(xì)胞層面過(guò)濾

用subset()函數(shù)對(duì)Seurat對(duì)象進(jìn)行設(shè)置：

# 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ū)用孢^(guò)濾

刪除0表達(dá)值的基因和在少于10個(gè)細(xì)胞中表達(dá)的基因。

# 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 objectfiltered_seurat <- CreateSeuratObject(filtered_counts, meta.data = filtered_seurat@meta.data)

Re-assess QC metrics

過(guò)濾之后可以回顧一下過(guò)濾前后的效果簿姨，以評(píng)判過(guò)濾指標(biāo)是否合適。

Saving filtered cells

# Create .RData object to load at any timesave(filtered_seurat, file="data/seurat_filtered.RData")

（五）Count Normalization and Principal Component Analysis

獲得高質(zhì)量的單細(xì)胞后簸搞，單細(xì)胞RNA-seq（scRNA-seq）分析工作流程的下一步就是執(zhí)行聚類扁位。聚類的目標(biāo)是將不同的細(xì)胞類型分成獨(dú)特的細(xì)胞亞群。為了進(jìn)行聚類趁俊，我們確定了在細(xì)胞之間表達(dá)差異最大的基因域仇。然后，我們使用這些基因來(lái)確定哪些相關(guān)基因集是造成細(xì)胞之間表達(dá)差異最大的原因寺擂。

數(shù)值歸一化

歸一化最重要的目的就是使表達(dá)水平在細(xì)胞之間和/或細(xì)胞內(nèi)更具有可比性暇务。那么在歸一化中主要需要處理的因素包括：

測(cè)序深度：考慮測(cè)序深度是比較細(xì)胞之間基因表達(dá)的必要條件。在下面的示例中怔软，每個(gè)基因在細(xì)胞2中的表達(dá)似乎都增加了一倍垦细，但這是細(xì)胞2具有兩倍測(cè)序深度的結(jié)果。

因此挡逼，要準(zhǔn)確比較細(xì)胞之間的表達(dá)括改，有必要對(duì)測(cè)序深度進(jìn)行歸一化。

基因長(zhǎng)度：需要基因長(zhǎng)度來(lái)比較同一細(xì)胞內(nèi)不同基因之間的表達(dá)家坎。比對(duì)到較長(zhǎng)基因的reads數(shù)似乎與較高表達(dá)的較短基因具有相同的計(jì)數(shù)/表達(dá)嘱能。

主成分分析（PCA）

PCA是一個(gè)用來(lái)強(qiáng)調(diào)差異和相似性的技術(shù)（（https://www.youtube.com/watch?v=_UVHneBUBW0）

下面是PCA的示例模擬過(guò)程苏携，幫助理解：

如果你已經(jīng)定量了兩個(gè)樣本（或細(xì)胞）中四個(gè)基因的表達(dá)，則可以繪制這些基因的表達(dá)值对粪，其中一個(gè)樣本在x軸上表示右冻，另一個(gè)樣本在y軸上表示，如下所示：

我們可以沿代表最大變化的方向在數(shù)據(jù)上畫(huà)一條線衩侥，在此示例中為對(duì)角線国旷。數(shù)據(jù)集中的最大變異是在組成兩個(gè)端點(diǎn)的基因。我們還看到基因在該線的上方和下方有些不同茫死。我們可以在數(shù)據(jù)上繪制另一條線跪但，代表數(shù)據(jù)中變化第二大的變量。

末端附近的基因?qū)⑹亲儺愖畲蟮幕蚵臀＿@些基因在數(shù)學(xué)上對(duì)線的方向影響最大屡久。

例如，基因C值的微小變化將極大地改變較長(zhǎng)線的方向爱榔，而基因A或基因D的微小變化對(duì)其幾乎沒(méi)有影響被环。

我們還可以旋轉(zhuǎn)整個(gè)圖，保證線條方向是從左到右和從上到下∠暧模現(xiàn)在筛欢，可以將這些線視為代表變化的軸。這些軸本質(zhì)上是“主要組分”唇聘，其中PC1代表數(shù)據(jù)的最大差異版姑，PC2代表數(shù)據(jù)的第二大差異。

如果有N個(gè)細(xì)胞迟郎，以此類推剥险。。宪肖。

確定PCs后表制，則需要對(duì)每個(gè)PC進(jìn)行評(píng)分，按照以下步驟對(duì)所有樣本PC對(duì)（sample-PC pairs）計(jì)算分?jǐn)?shù)：

（1）首先控乾，根據(jù)基因?qū)γ總€(gè)PC的影響程度么介，為其分配“影響力”評(píng)分。對(duì)給定PC沒(méi)有任何影響的基因得分接近零阱持，而具有更大影響力的基因得分更高夭拌。PC線末端的基因?qū)a(chǎn)生更大的影響，因此它們將獲得更大的分?jǐn)?shù)，但兩端的符號(hào)相反鸽扁。

（2）確定影響分?jǐn)?shù)后蒜绽，使用以下公式計(jì)算每個(gè)樣本的分?jǐn)?shù)：

Sample1 PC1 score = (read count * influence) + ... for all genes

以我們的2個(gè)樣本示例，以下是分?jǐn)?shù)的計(jì)算方式：

## 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) = 21PC2 score = (5 * 0.5) + (4 * 1) + (8 * -5) + (7 * 6) = 8.5

（3）一旦為所有PC計(jì)算了這些分?jǐn)?shù)桶现，就可以將其繪制在簡(jiǎn)單的散點(diǎn)圖上躲雅。下面是示例圖：

對(duì)于具有大量樣本或細(xì)胞的數(shù)據(jù)集，通常會(huì)繪制每個(gè)樣本/細(xì)胞的PC1和PC2分?jǐn)?shù)骡和。由于這些PC解釋了數(shù)據(jù)集中最大的變化相赁，因此更相似的樣本/細(xì)胞將在PC1和PC2聚在一起。請(qǐng)參見(jiàn)下面的示例：

Image credit: https://github.com/AshwiniRS/Medium_Notebooks/blob/master/PCA/PCA_Iris_DataSet.ipynb

對(duì)于我們的單細(xì)胞數(shù)據(jù)慰于，我們最終會(huì)選擇10-100 PC去對(duì)細(xì)胞進(jìn)行聚類分析钮科，而不是全部基因。

未完待續(xù)婆赠。

你可能還想看