本文介绍了单细胞RNA-seq技术及其相对于传统RNA-seq的优势，探讨了不同单细胞测序方法如SMART-seq2、Drop-seq等，强调了Seurat在数据处理和分析中的关键作用。Seurat用于质控、分析和数据探索，能完成从原始数据到表达矩阵的构建。文章还讨论了实验方法的挑战，如细胞选择、文库构建的差异以及如何选择合适的平台，并提供了数据处理的初步步骤，包括FastQC质量控制、STAR比对和Kallisto的伪比对。
单细胞RNA-seq简介
混合RNA-seq2000年末的重大技术突破，取代微阵列表达芯片被广泛使用
通过混合大量细胞获取足够RNA用于建库测序，来定量每个基因的平均表达水平
用于比较转录组，例如比较不同物种的同一组织样本
量化整体表达特征，如疾病研究中的表达模式
研究异质系统方面还有力所不及之处，例如对早期发育的研究,复杂组织(大脑)的研究
在基因表达随机性研究方面心有余而力不足
scRNA-seq是一项由汤富酬等人在2009年首次发表的新技术。文章发表于Nature Method，测序了7个单细胞，两个卵裂球，两个野生型卵子，两个Dicer敲除的卵 子，一个Ago2敲除的卵子。
这项技术在2013年被Nature评为年度技术,更简便的操作流程和较低的测序成本促成单细胞技术的广泛流行。2018年底，单细胞技术应用于胚胎发育追踪评为Science年度突破。
检测每个基因在大量细胞中的表达水平分布。
可以研究细胞类型特异性转录调控的新型生物问题，例如细胞类型鉴定，细胞应答的异质性，细胞表达的随机性，细胞间基因调控网络的推断等
研究中细胞数目范围从100个变到10^6个且每年递增。
目前有许多不同的单细胞Protocol，例如 SMART-seq2 , CELL-seq 和 Drop-seq 。
还有商业平台，包括 Fluidigm C1, Wafergen ICELL8和the 10X Genomics Chromium。
Bulk RNA-seq技术中一些计算分析方法可应用于单细胞分析。
多数情况下单细胞计算分析需要调整现有方法或者开发新方法
工作流程
总体而言，scRNA-seq的实验方案和bulk RNA-seq的相似。我们将在下一节一起讨论一些最通用的方法。
计算分析
本课程内容是scRNA-seq实验中得到的数据进行计算分析。总体流程如下图所示，前面三步(黄色)对于任何高通量测序数据是通用的，紧随其后的四步(橙色)是要将传统RNA-Seq分析中已有的方法和新开发的方法结合起来解决scRNA-seq的技术差异问题，最后的部分(蓝色)是使用专门为scRNA-seq开发的方法来进行生物分析解读。
scRNA-seq分析的综述有几篇，包括 Computational and Analytical Challenges in Single-Cell Transcriptomics.” Nat Rev Genet 16 (3) 。
目前还有其他平台可以执行上述流程图中的一步或多步操作：Falco：是一个单细胞RNA-seq的云处理平台，更像是一个流程部署和管理工具，一年多未更新了，一般也用不上。能部署的应该都有自己 的一套部署工具，初学者不需要学这么复杂的。有精力，可以学习下其部署理念应用于自己的流程。
SCONE(Single-Cell Overview of Normalized Expression)：单细胞RNA-seq质量控制和标准化的R包 (一年多没更新了, Yosef研究 组2018年在Nature method发表一个单细胞分型的深度学习平台，scVI，效果不错，值得尝试)Seurat ：单细胞质控，分析和数据探索而设计的R包，可以完成获得定量数据后的几乎所有分析。不少文章的几个主图都是来自这个软件包 。这个软件包可以作为学习的入门，官网的教程示例写的很详细。
ASAP(Automated Single-cell Analysis Pipeline) ：是一款单细胞分析的交互式网络平台。从基因表达矩阵开始到后期分析。功能相对比较全，定制化弱一些。学完这份教程，里面的功能都可以自己实现。
挑战
Bulk RNA-seq和scRNA-seq的主要差别是每个测序文库代表一个单细胞还是一群细胞。比较不同细胞(不同测序文库)的结果需要格外注意。文库之间差异的主要来源是：扩增效率和扩增偏好性(部分文库可扩增多达100万倍)
基因 ‘dropouts’: 基因在一个细胞中呈现中等表达水平,但在另一个细胞中未检测到表达，这可能来源于scRNA-seq中RNA总量低导致的扩增建库丢失或RNA表达的随机性。
取自于单独一个细胞的低转录本总量是这两个文库差异的一个主要原因。提高转录本捕获效率和降低扩增偏好可以降低差异，是目前活跃的研究方向。从后续课程学习中也可以看 到，合适的标准化和校正方法也可以抵消一部分文库构建引入的噪音。
单细胞测序实验方法
开发scRNA-seq的新实验方法和操作手册是目前很火的研究领域，而且最近这些年已经发表了一些改进方法。上图可以看出检测的细胞数目以指数形式增加，以下是一份不完全的列表：STRT-seq
Smart-seq
CEL-seq
MARS-seq
SCRB-seq
Smart-seq2
Drop-seq
InDrop-seq
CEL-seq2
Seq-well
SMARTer
Microwell-seq 【2018，郭国冀团队】
总体流程如下图所示：
这些方法可以用不同的方式分类，但两个最主要的优化是：定量和捕获。
定量有两种类型—-全长型和标签型 (tag-based)。全长型力图捕获并均匀测序整条转录本，标签型只捕获转录本的5'或3'端。不同定量方式需要自己对应的计算分析方法。全长方案理论上可以对整个转录本进行均匀测序，但实际上总会有测序覆盖偏好性的存在。标签型的主要优点是可以与唯一的分子标识符(UMIs)结合进行更精确定量 (后面详细描述)。其缺点是，测序限制在转录本的5'或3'端，可能会降低比对率，并且难以区分不同剪接体的表达。
捕获的策略决定了实验通量、细胞如何被选择和除测序外的哪些额外信息可获得。最常用的三种捕获方式是基于微孔- (microwell-)，微流- (microfluidic-)，液滴- (droplet-)，细分如下图所示。
对于基于微孔的捕获平台，先用移液管或者激光切割的方式分离细胞并放到微流孔中。它的一个优点是可以结合流式细胞荧光分选(FACS, fluorescent activated cell sorting)根据表面Marker分选细胞。因此特别适合分选细胞子集用于测序。另一个优点是可以获得细胞形态全览图，提供多一个维度的信息，可用于鉴定微孔中是否有损伤的细胞或双份细胞，主要缺点是通量低且每个细胞所需的工作量相当大。
微流型平台，比如Fluidigm’s C1，提供了一个更加整合的系统，同时可以捕获细胞和完成文库构建的准备过程。因此它们比微孔型平台通量更高。但是微流型平台大约只能捕获10%的细胞，不适合处理稀有细胞或者细胞量很少的情况。此外，微流控芯片相对昂贵，不过更小的反应体积可以节省试剂的费用。
液滴型方法是将单独的细胞和一个包含建库所学酶的珠粒 (bead)包裹在一个纳米级液滴里面。特殊地，每个珠粒(bead)包含一段独特的条形码序列 (barcode)，会加到所有来自于液滴里面这个细胞的序列上，用于区分不同细胞的转录本。因此所有的液滴可以混合在一起测序，然后再根据barcode序列确定其是否归属于同一细胞。液滴型平台通常有最高的通量，因为文库的准备成本很低，约为0.05美元/每个细胞。随之而来的，测序成本往往是其限制因素，通常测序深度比较低，只检测几千个转录本的表达。
Microwell-seq是浙江大学郭国冀老师研究组综合以上技术的优势开发出的新的大规模低成本单细胞捕获测序技术，单个细胞制备成本可以到0.02美元。
采用光刻技术制作微孔矩阵硅片(微孔直径28 um，深度35 um，100,000个微孔)，以此为模具制作PDMS微柱模具。这两个模具可以反复使用。最终用于富集的微孔板是通过倾到5%的琼脂糖凝胶到PDMS微柱模具上生成的。细胞悬液加到凝胶微孔模具上，利用重力使细胞落入微孔，通常一个微孔只能容纳一个细胞，一块板子可以同时捕获约10000个单细胞。每一步操作都可视、可控制，doublets可以通过镜检洗除。随后每个空加入包含10^7-10^8特定探针集的与孔径大小匹配的磁珠，标记每个细胞中的mRNA(每个磁珠的寡核苷酸序列中都有一段特异的序列用于标记细胞来源)，然后使用Smart-seq2方法进行后续的反转录、扩增。扩增后的cDNA片段使用转座酶片段化(这步倒有些类似ATAC-seq)，富集3'末端转录本序列测序。
如何选择合适的平台？
平台选择取决于手上的生物问题，举个例子，如果一个人对描述组织的细胞类型构成感兴趣，能够捕获大量细胞的液滴型平台很可能是最合适的，而如果感兴趣的是有已知的表面Marker的稀有细胞群体，那最好用流式细胞FACS分选法富集细胞并对少量细胞进行测序。
如果对研究不同剪接体的表达感兴趣，全长转录本定量会更适合。相反，UMIs只适用于tagged protocol，更有利于稳定定量基因水平的表达。
Enard团队和Teichmann团队的最近两项研究对几种不同单细胞分选和建库方案进行了比较，Ziegenhain等用同一种老鼠胚胎干细胞 (mESCs)比较了五种不同方案。通过控制细胞数量和测序深度，作者能直接比较不同方法检测敏感性，噪音水平和费用差异。他们的一个结论如下图，不同方法检测到的基因数量 (检测阈值固定)差别挺大。如图所示，drop-seq检测到的基因数目最少，和Smart-seq2检测到的基因数差了近乎两倍，意味着不同方案的选择对研究影响很大。
Svensson等人则采用了一种不同的方法，即通过使用已知浓度的合成转录本 (spike-ins, 后面详细介绍)来评估不同方案的准确性和敏感性。通过广泛比较同样发现不同的细胞分离建库方式差别较大。Bulk测序的准确性和敏感性相对最好，其它方法准确性高的，敏感性就会差一些。
随着实验操作技术的发展和计算方法的改进，后续研究可以帮助我们进一步了解不同方法的适用优缺点。这些比较研究不仅有助于研究人员决定使用哪种方案，而且可以通过基准测试确定哪个方法组合最有效来研发新的单细胞实验和计算方法。
scRNA-seq原始数据加工
FastQC
得到单细胞RNA-seq测序数据后，首先检查测序reads的质量。为了完成这个任务，我们使用的工具是FastQC。FastQC是一款质控工具，能对bulk RNA-seq和单细胞RNA-seq的原始数据进行质量控制 (其他类型的高通量测序结果也适用)。FastQC以原始测序reads为输入(fastq格式)，输出序列质量报告。复制粘贴下面的链接到你的浏览器进入FastQC官网：
https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
这个网址包含下载和安装FastQC及示例报告文件的链接。向下滚动页面到Example Reports并点击Good Illumina Data，会看到高质量Illumina测序Reads的理想质控结果。如果使用Docker镜像，则FastQC已经安装好。如果是自己的服务器，FastQC下载下来即可使用(依赖于Java)。生信宝典的推文NGS基础 - FASTQ格式解释和质量评估对FASTQ原始数据和FastQC的使用和结果描述有比较详细的介绍，如果不熟悉，建议阅读。
测序文库拆分 (Demultiplexing)
文库拆分因使用的前期Protocol不同或构建的流程不同需要有对应的处理方式。我们认为最灵活可用的文库拆分工具是zUMIs (https://github.com/sdparekh/zUMIs/wiki/Usage)，可以用来拆分和比对大部分基于UMI的建库方式。对于Smartseq2或其他双端全长转录本方案，数据通常已经拆分好了。例如GEO或ArrayExpress之类的公共数据存储库会要求小规模或plate-based scRNASeq数据拆分好再上传，并且很多测序服务商提供的数据都是自动拆分好的。如果使用的分析流程依赖于拆分好的数据但测序服务商提供的数据没有拆分时就需要自己拆分。因为不同的建库方案引入的barcode序列的长度和位置不同，通常都需要自己写脚本解决。
对于所有数据类型，”文库拆分”涉及从一端或双端短序列中识别和移除细胞条形码序列 (cell barcode)。如果提前知道加入的细胞条形码，比如数据来自基于PCR板的方案，只需要找到条形码并与条形码库作比对，归类于与之最相似的那个就可以 (根据条形码的设计，一般允许最多1-2个错配)。这些数据通常在比对之前先做拆分，从而可以并行比对，提高效率。
鉴定含有细胞的液滴/微孔
液滴型scRNA-seq方法中只有一小部分的液滴包含珠状物和一个完整细胞。然而生物实验不会那么理想，有些RNA会从死细胞或破损细胞中漏出来。所以没有完整细胞的液滴有可能捕获周围环境游离出的少了RNA并且走完测序环节出现在最终测序结果中。液滴大小、扩增效率和测序环节中的波动会导致”背景”和真实细胞最终获得的文库大小变化很大，使得区分哪些文库来源于背景哪些来源于真实细胞变得复杂。
大多数方法使用每个barcode对应的总分子数(如果是UMI)或总reads数的分布来寻找一个”break point”区分来自于真实细胞的较大的文库和来自于背景的较小的文库。下面加载一些包含大细胞和细胞的模拟数据实际操作下：
umi_per_barcode
truth
使用STAR比对read
现在我们已经对测序原始数据进行了质控，获得了高质量的Clean data，下一步就是把它们比对到参考基因组。如果我们想定量基因表达或鉴定样本之间差异表达的基因，则通常需要某种形式的比对。
用于短序列比对的工具已经开发了很多(转录组分析工具哪家强？)，但今天我们主要涉及两个。第一个工具是STAR。对于测序数据中的每条reads，STAR尝试找到能与参考基因组中一个或多个位置匹配的最长可能序列。例如，在下图中，有一个短序列(蓝色)，它跨越两个外显子和一个选择性剪接点(紫色)。STAR能够发现短序列的第一部分与第一外显子的序列匹配，而第二部分与第二外显子中的序列匹配。因为STAR能够用这种方式识别剪接事件，所以它被称为splice aware的比对工具 (一般转录组分析的比对工具都需要有这个功能)。
通常STAR把短序列比对到参考基因组时允许检测新的剪接事件或染色体重排事件。然而，STAR的一个问题是它需要大量的内存，尤其是参考基因组很大(例如老鼠和人类，大约需要30 G内存)的时候。为了加速今天的分析，我们将使用STAR把reads比对到只包含2000个转录本的参考转录组上。请注意，这不是常用或推荐的做法，只是考虑时间问题才这样做。我们通常建议比对到参考基因组 (但过程与此类似)。
执行STAR比对需要两个步骤。在第一步中，用户向STAR提供参考基因组序列(FASTA)和注释(GTF)(见NGS基础 - 参考基因组和基因注释文件)，来创建基因组索引 (生信宝典注：建索引的目的是为了提高比对速度)。第二步，STAR将用户的短序列数据比对到基因组索引。
Kallisto和Pseudo-Alignment
STAR是序列比对工具，而Kallisto是伪比对工具 [@bray_2016]。它们的区别是：比对工具是把reads比对回参考基因组或转录组，而伪比对工具是把k-mers比对到参考转录组。
什么是k-mer？
k-mer是来源于测序短序列中的长度为k的子序列。例如，假设有短序列ATCCCGGGTTAT，想从中获得7-mer。为此，我们将提取前七个碱基作为第一个7-mer，然后向下移动一个碱基获得第二个7-mer，以此类推。下图显示了从序列中可以得到的所有7-mers：
为什么比对k-mers而不是reads？
有两个主要原因：伪比对工具利用算法技巧使得比对k-mers比比对reads速度快很多，具体可以参阅 [@bray_2016]了解详情。
在某些情况下，伪对齐工具可能比传统比对工具更好的处理测序错误问题。例如，假设序列上第一个碱基中存在测序错误，本来是T却测序成了A。对伪比对工具来说，只会影响第一个7-mer而不会影响后续7-mer的比对。
Kallisto的伪比对模式
Kallisto有一个为单细胞转录组特别设计的伪比对模式。和STAR不同，Kallisto比对到的是参考转录组而不是参考基因组，意味着Kallisto是将序列比对到剪接异构体而不是基因上，对单细胞来讲，这是有挑战性的：单细胞RNA-seq的覆盖率低于bulk RNA-seq，这意味着可以从序列中获得的信息总量减少了。
许多单细胞RNA-seq方案具有3'覆盖偏好性，如果两种剪接异构体只在5’末端不同，则很难确定序列来源于哪个剪接异构体。
一些单细胞RNA-seq方案测序读长短，也难以区分来源于哪个剪接异构体。
Kallisto的伪模式采用了略微不同的伪比对方法。它不与剪接异构体比对，而是与等价类 (equivalence classes)比对。所以如果read比对到多个异构体，Kallisto会记录read比对到包含有所有异构体的等价类，因此可以使用等价类计数而非基因或转录本计数用于下游的聚类等分析。具体见下图：
构建表达矩阵
scRNA-seq数据的许多分析以表达矩阵为起点。一般来讲，表达矩阵的每一行代表一个基因，每一列代表一个细胞(但是一些作者会做个转置)。每个条目代表特定基因在给定细胞中的表达水平。而表达值的测量单位取决于建库方案和所用的标准化方法。
reads质控
见前面章节FastQC部分。
另外，使用Integrative Genomics Browser(IGV)或SeqMonk通常对数据可视化很有帮助，具体见下。测序数据可视化 (一)
IGV基因组浏览器可视化高通量测序数据
高通量数据分析必备-基因组浏览器使用介绍 - 1
高通量数据分析必备-基因组浏览器使用介绍 - 2
高通量数据分析必备-基因组浏览器使用介绍 - 3
Reads比对
见前面章节STAR部分和Kallisto部分。
注释的很好的模式生物(例如小鼠，人)有着大量全长转录本数据集，伪比对方法(例如Kallisto，Salmon)可能优于常规比对方法。drop-seq方法获得的数据集有数以千万条reads，伪比对工具的运行时间比传统比对工具会少几个数量级，更有时间优势。从39个转录组分析工具，120种组合评估(转录组分析工具哪家强-导读版)一文中可以看出，伪比对工具的准确性和稳定性也相对比较高。
用STAR比对的操作示例 (前面章节部分更详细)
STAR --runThreadN 1 --runMode alignReads
--readFilesIn reads1.fq.gz reads2.fq.gz --readFilesCommand zcat --genomeDir
--parametersFiles FileOfMoreParameters.txt --outFileNamePrefix /output
注意，如果用了spike-ins(已知浓度的外源RNA分子)，在比对前应该将参考基因组和spike-in分子的DNA序列合并作为共同”参考基因组”。具体见之前文件格式部分。
注意，使用UMI时，应从read序列中删除其条形码。常见的是将条形码加到read名称上。
一旦reads完成了到基因组的比对，我们需要检查比对率和确保有足够多的reads比对回了参考基因组。根据我们的经验，小鼠或人类细胞中read的比对率为60-70％。但是这个结果可能会因建库方案、read长度和比对工具参数设置而有所不同。常规上，我们希望所有细胞都具有相似的read比对率。如果有样品比对率异常低或比对回去的reads异常低，则需要多加注意甚至从后续分析中移除。较低的read比对率通常表示存在污染。生信宝典建议取出几十条未必对回去的reads做个blast，看下能否比对到其它物种来确定是污染还是测序错误还是比对参数设置的问题。
一个用Salmon量化表达操作的例子
salmon quant -i salmon_transcript_index -1 reads1.fq.gz -2 reads2.fq.gz -p #threads -l A -g genome.gtf --seqBias --gcBias --posBias
注意 Salmon操作会得到一个估计的read counts和transcripts per million(tpm)。根据我们的经验，TPM对单细胞测序中长基因的表达做了过度校正，因此我们建议使用read counts。
比对示例
下面直方图 (http://www.ehbio.com/ImageGP)显示了scRNA-seq实验中每个细胞不同比对状态的reads的数目。每个柱子代表一个细胞，按细胞的总read数升序排列。三个红色箭头标记的是比对到基因组的reads较低的异常样本，应该在后续分析中移除。两个黄色箭头指的是unmapped reads数目十分大的细胞。该例中，在比对质控期间这两个细胞会保留下来，但后期细胞质控时这两个细胞会因为核糖体RNA reads比例过高而移除。
Mapping QC
在把原始序列比对到基因组后，需要评估比对质量。这可以从多个角度进行评估，包括：rRNA/tRNAs的reads的占比或总量，reads在基因组上唯一比对位置的比例，比对到splice junction的reads比例，reads在转录本的覆盖均一性或深度。而为Bulk RNA-seq开发的方法如RSeQC，也适用于单细胞数据：
#pip install RSeQC
geneBody_coverage.py -i input.bam -r genome.bed -o output.txt
bam_stat.py -i input.bam -r genome.bed -o output.txt
split_bam.py -i input.bam -r rRNAmask.bed -o output.txt
然而，预期结果的评估取决于采用的建库方案，例如许多scRNA-seq用poly-A selection捕获转录本。这个方法可以排除核糖体RNA污染，但会导致3'区域更容易测到。下图展示了测序reads分布的3'偏好性，和去除的三个异常细胞的结果 (应该是最下面3条，推测是降解严重)。
Reads量化
scRNA-seq基因定量计算可以用bulk RNA-seq一样的工具，比如HT-seq or FeatureCounts。
# include multimapping
featureCounts -O -M -Q 30 -p -a genome.gtf -o outputfile input.bam
# exclude multimapping
featureCounts -Q 30 -p -a genome.gtf -o outputfile input.bam
唯一分子标识符UMI让计算转录本的绝对量成为可能，并且在scRNA-seq中很受欢迎。我们将在下一章讨论如何处理UMI。
唯一分子标识符(Unique Molecular Identifiers, UMI)
感谢EMBL Monterotondo的 Andreas Buness 在本节的合作。
UMI添加到每个转录本上
唯一分子标识符 (UMI)是在反转录过程中添加到转录本上的短的(4-10 bp)随机条形码序列。它们使得测序reads可以对应到单个转录本，从而去除扩增噪声和偏好性。
当测序含UMI的文库时，仅对包含UMI的转录本末端 (通常为3'末端)进行性测序。
比对UMI条形码
由于UMI数量(, N是UMIs的长度值)比每个细胞中的RNA分子数(~)少得多，每个UMI条形码可能会连接到多个转录本，因此需要借助条形码序列和reads比对位置两个条件鉴定起始的转录本分子。第一步是比对UMI reads，推荐用STAR来处理，因为它处理速度快且输出高质量的BAM比对。此外，比对位置的准确性对识别新的3'UTR区域很有意义。
UMI测序通常由双端reads组成，其中一端read是捕获细胞和UMI的条形码，而另一端read包含转录本的外显子序列。注意，推荐去除reads中的poly-A序列部分，以免这些reads比对到转录本内部poly-A或poly-T序列而产生错误。
处理UMI实验中的reads，通常有以下惯例：UMI被添加到另一个配对read的序列名称中。
reads按细胞条形码分类到单独的文件中 (见前面的文章)。但对于细胞量极大的低深度测序数据集 (drop-seq)，可以将细胞条形码添加到read名称中而不是拆分为单独文件以减少文件数量。
Counting 条形码
理论上，每个唯一的UMI-转录本对应该对应来源于一个RNA分子的所有reads。但是现实往往并非如此，最常见的原因是：不同的UMI序列不一定表示它们是不同的UMI分子由于PCR或测序错误，碱基替换可能导致新的UMI序列。较长的UMI出现碱基替换的机会更多。根据细胞条码测序错误估计，7-10％的10 bp长度的UMI中至少有一个碱基替换。如果错误没有纠正，将会过高估计转录本的数目。
不同的转录本不一定是不同的分子比对错误或多个比对位置可能导致某些UMI对应到错误的基因/转录本。这种类型的错误也会导致过高估计转录本的数目。
相同的UMI不一定意味着相同的分子UMI频次的不同和短UMI可导致同一UMI和相同基因的不同mRNA分子相连，进而可能低估转录本数量。
错误纠正
如何最好的校正UMIs中的这些问题仍然是一个活跃的研究领域。我们自己认为的最好的解决上述问题的方法是：UMI-tools’，设计了directional-adjacency算法，同时考虑错配数目和相似UMI的相对频率来识别PCR和测序错误。(alevin, cellranger都是不错的选择，后面详细介绍)
问题还无法完全解决，但通过删除只有很少read支持的UMIs-转录本对，或者移除所有多比对位置的reads，可能会减轻该问题。
Simple saturation (也称为”collision probability”)方法来估计分子的数量
其中N=唯一UMI条形码的总数，n=观察到的条形码数。
这个方法的一个重要缺陷是它假设所有UMI出现频率相同。但因为序列GC含量不同引入的偏差使得这一假设在大多数情况下这是不正确的。
如何最好地处理和使用UMI在目前生物信息学界是一个活跃的研究领域。而我们了解到的几种最近开发的方法有：UMI-tools
PoissonUMIs
zUMIs
dropEst
下游分析
当前的UMI平台(DropSeq，InDrop，ICell8)展现出从低到高变化很大的捕获效率，如下图所示。
这一高可变性可能会引入很强的偏差，需要在下游分析时考虑到。现在的分析通常根据细胞类型或生物通路把细胞/gene混合一起增加检测能力。更合适的统计分析方法亟待研究以便更好地调整这些偏差，使得结果更能反映真实现象。
练习1 数据是三个不同来源的诱导多功能干细胞的UMI counts和read counts (有关此数据集的详细信息请参阅后续文章)。
umi_counts
read_counts
使用此数据：绘制捕获效率的变化
确定扩增率：每个UMI对应的平均reads数。
# Exercise 1
# Part 1
plot(colSums(umi_counts), colSums(umi_counts > 0), xlab="Total Molecules Detected", ylab="Total Genes Detected")
# Part 2
amp_rate
amp_rate
Tabula Muris
Tabula Muris是测序小鼠20个器官和组织的单细胞转录组图谱的国际合作项目 (Transcriptomic characterization of 20 organs and tissues from mouse at single cell resolution creates a Tabula Muris)。
简介
我们使用 Tabula Muris最开始释放的数据做为测试数据来完成完整的单细胞数据分析。The Tabula Muris是一个国际合作组织，目的是采用标准方法生成小鼠每个细胞的图谱。建库测序方法包括通量高覆盖率低的10X数据和通量低覆盖率高的FACS筛选+Smartseq2建库技术。
起始数据于2017年12月20日释放，包含20个组织/器官的100,000细胞的转录组图谱。
下载数据
与其它sc-RNASeq数据上传到GEO或ArrayExpress不同，Tabula Muris通过figshare平台释放数据。可以通过搜索文章的doi号10.6084/m9.figshare.5715040下载FACS/Smartseq2数据和10.6084/m9.figshare.5715025下载10X数据。数据可以直接点击链接下载，或使用下面的wget命令：
终端下载FACS数据:
wget https://ndownloader.figshare.com/files/10038307
unzip 10038307
wget https://ndownloader.figshare.com/files/10038310
mv 10038310 FACS_metadata.csv
wget https://ndownloader.figshare.com/files/10039267
mv 10039267 FACS_annotations.csv
终端下载10X数据:
wget https://ndownloader.figshare.com/files/10038325
unzip 10038325
wget https://ndownloader.figshare.com/files/10038328
mv 10038328 droplet_metadata.csv
wget https://ndownloader.figshare.com/files/10039264
mv 10039264 droplet_annotation.csv
如果数据是手动下载的，也需要像上面一样解压和重命名。
现在应该有两个文件夹: FACS和droplet，每个对应一个annotation和metadata文件。使用head命令查看前10行：
head -n 10 droplet_metadata.csv
使用wc -l查看文件的行数:
wc -l droplet_annotation.csv
练习：FACS和10X数据中各有多少细胞有注释信息？
答案: FACS : 54,838 cells; Droplet : 42,193 cells
读入数据 (Smartseq2)
读入逗号分隔的count matrix，存储为数据框:
dat = read.delim("FACS/Kidney-counts.csv", sep=",", header=TRUE)
dat[1:5,1:5]
数据库第一列是基因名字，把它移除作为列名字：
dim(dat)
rownames(dat)
dat
这是Smartseq2数据集，可能含有spike-ins:
rownames(dat)[grep("^ERCC-", rownames(dat))]
从列名字中提取metadata信息:
cellIDs
cell_info
Well
Well
Plate
Mouse
检测每种metadata类型的数据分布：
summary(factor(Mouse))
查看有没有技术因子是cofounded，实验批次与供体小鼠批次一致:
table(Mouse, Plate)
最后读入计算预测的细胞类型注释，并与表达矩阵中的细胞注释做比较:
ann
ann
celltype
table(celltype)
构建scater对象
为了构建SingleCellExperiment对象，先把所有的细胞注释放到一个数据框中。因为实验批次(pcr plate)和供体小鼠完全重合，所以只保留一个信息。
suppressMessages(require("SingleCellExperiment"))
suppressMessages(require("scater"))
cell_anns
rownames(cell_anns)
sceset
str(sceset)
如果数据集包含spike-ins，我们在SingleCellExperiment对象中定义一个变量记录它们。
isSpike(sceset, "ERCC")
str(sceset)
读入10X的数据
因为10X技术细胞通量高但测序覆盖度低，所以其count matrix是一个大的稀疏矩阵(矩阵中高达90%的数据的数值为0)。CellRanger默认的输出格式是.mtx文件用于存储这个稀疏矩阵，第一列是基因的坐标(0-based)，第二列是细胞的坐标(0-based)，第三列是大于0的表达值 (长表格形式)。 打开.mtx文件会看到两行标题行后面是包含总行数 (基因数)、列数 (样本数)和稀疏矩阵总行数 (生信宝典注：所有细胞中表达不为0的基因的总和)的一行数据。
%%MatrixMarket matrix coordinate integer general
%
23433 610 1392643
5 1 1
28 1 1
40 1 2
鉴于.mtx文件中只存储了基因和样品名字的坐标，而实际的基因和样品的名字必须单独存储到文件genes.tsv和barcodes.tsv。
下面使用Matrix包读入稀疏矩阵：
suppressMessages(require("Matrix"))
cellbarcodes
genenames
molecules
下一步增加合适的行或列的名字。首先查看read的cellbarcode信息会发现这个文件只有barcode序列。考虑到10X数据每一批的cellbarcode是有重叠的，所以在合并数据前，需要把批次信息与barcode信息合并一起。
head(cellbarcodes)rownames(molecules)
colnames(molecules)
读入计算注释的细胞类型信息:
meta
head(meta)
我们需要用10X_P4_5获得这批数据对应的metadata信息。这时需要注意metadata表格中mouse ID与前面plate-based (FACS SmartSeq2)数据集的mouse ID不同，这里用-而非_作为分隔符，并且性别在中间。通过查阅文献中的描述得知droplet (10X)和plate-based (FACS SmartSeq2)的技术用了同样的8只老鼠。所以对数据做下修正，使得10X与FACS的数据一致。
meta[meta$channel == "10X_P4_5",]
mouseID
注意：有些组织的10X数据可能来源于多个小鼠的样品，如mouse id = 3-M-5/6。也需要格式化这些信息，但可能这些与FACS数据的mouse id会不一致，进而影响下游分析。如果小鼠不是纯系，可能需要通过exonic-SNP把细胞和对应的小鼠联系起来 (本课程不会涉及)。
ann
head(ann)
注释中的cellID和cellbarcodes也存在细微差别，少了最后的-1，在匹配前需要做下校正。(生信宝典注：这种数据不一致是经常要处理的问题，每一步检查结果。如果与预期不符，考虑有没有未考虑到的数据不一致的地方。)
ann[,1]
ann_subset
celltype
构建cell-metadata数据框:
cell_anns
rownames(cell_anns)
head(cell_anns)
表达矩阵质控
UMI表达定量 (UMI)
UMI表达定量简介
基因定量后会整理成一个行为基因(或转录本)列为细胞的表达矩阵。虽然前面做了原始数据质控和测序数据质控移除了一部分从reads数层面就不合格的细胞，还需要进一步根据表达矩阵移除其它类型低质量细胞。如果未能识别并移除低质量细胞会混淆下游分析中的有意义的生物信息。
鉴于scRNASeq还没有标准方法，后续用到的质控质量值的标准会因为实验方法不同而有很大变化。因此，执行质控时，我们是通过数据集内部比较找到异常细胞，而不是依赖于其它独立的质量标准。因此比较不同的建库方法获得的不同数据集时需要格外注意。
Tung数据集
我们使用芝加哥大学Yoav Gilad实验室的3个不同来源的诱导多能性干细胞 (iPSC)的数据集 (http://jdblischak.github.io/singleCellSeq/analysis/) in Yoav Gilad。细胞分选采用Fluidigm C1微流控台，同时使用UMIs和ERCC spike in进行质控为了保证可重复性，数据是2016年3月15生成的原始数据的拷贝，存储于tung文件夹下。
library(SingleCellExperiment)
library(scater)
options(stringsAsFactors = FALSE)
读入数据和注释：
molecules
anno
查看读入的数据：
head(molecules[ , 1:3])
结果
## NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03
## ENSG00000237683 0 0 0
## ENSG00000187634 0 0 0
## ENSG00000188976 3 6 1
## ENSG00000187961 0 0 0
## ENSG00000187583 0 0 0
## ENSG00000187642 0 0 0
查看注释数据
head(anno)
结果
## individual replicate well batch sample_id
## 1 NA19098 r1 A01 NA19098.r1 NA19098.r1.A01
## 2 NA19098 r1 A02 NA19098.r1 NA19098.r1.A02
## 3 NA19098 r1 A03 NA19098.r1 NA19098.r1.A03
## 4 NA19098 r1 A04 NA19098.r1 NA19098.r1.A04
## 5 NA19098 r1 A05 NA19098.r1 NA19098.r1.A05
## 6 NA19098 r1 A06 NA19098.r1 NA19098.r1.A06
数据包含 3个个体，3次重复和9个批次。
通过使用SingleCellExperiment (SCE) 和scater包标准化分析过程。首先构建SCE对象:
umi
assays = list(counts = as.matrix(molecules)),
colData = anno
)
移除不在任何细胞表达的基因:
keep_feature 0) > 0
umi
定义对照特征基因集：ERCC spike-ins和线粒体基因 (作者提供，见http://jdblischak.github.io/singleCellSeq/analysis/qc-filter-ipsc.html):
isSpike(umi, "ERCC")
isSpike(umi, "MT")
c("ENSG00000198899", "ENSG00000198727", "ENSG00000198888",
"ENSG00000198886", "ENSG00000212907", "ENSG00000198786",
"ENSG00000198695", "ENSG00000198712", "ENSG00000198804",
"ENSG00000198763", "ENSG00000228253", "ENSG00000198938",
"ENSG00000198840")
计算质量矩阵:
umi
umi,
feature_controls = list(
ERCC = isSpike(umi, "ERCC"),
MT = isSpike(umi, "MT")
)
)
str(umi)
SingleCellExperiment对象结构。str (structure)是一个很好的工具，可以用来查看数据的结构组成。(RStudio中的View对于较大的对象会给出更好的展现方式。)
如下面的展示SingleCellExperiment有10个数据槽 (slots)，使用@索引。比如想获得降维后的结果，使用umi@reduceDims，使用umi@colData可以获取质控信息，都有哪些质控变量，进一步的使用umi@colData@listData$total_features_by_counts可以获得每个细胞的基因数。umi@rowRanges@elementMetadata@listData基因的属性信息。
(在Rstudio中有个便利的地方，输入变量名再输@, $可以弹出对应的子属性，方便快速输入。)
细胞质控
文库大小
查看每个样品(细胞)检测到的总分子数 (UMI count)或总reads数 (reads count)，拥有很少的reads或分子数的样品可能是细胞破损或捕获失败，应该移除。
hist(
umi$total_counts,
breaks = 100
)
abline(v = 25000, col = "red")
基因分析
基因表达
除了移除低质量的细胞，也会排除受技术操作影响较大的一部分基因。而且查看基因表达结果，可以帮助改进实验操作。
通常会看top 50表达的基因占据了多少reads。
plotHighestExprs(umi, exprs_values = "counts")
基因过滤
通常建议移除那些表达水平极低以至于可以视为”未检测出”的基因。这里针对UMI数据，“检出”定义为至少有2个细胞检测到某个基因存在多于一个转录本。如果是reads counts数据, “检出”可以定义为至少有2个细胞检测到某个基因有至少5个reads count支持。请注意，对两种表达量计算方式，阈值的选择都与测序深度有关。自己的数据可以做相应的修改。另外一个需要注意的点是基因的过滤必须在细胞过滤后面，因为部分基因可能只在低质量细胞中能检测的到 (注意下面的colData(umi)$use过滤).
数据可视化
简介
我们继续使用上一步生成的Tung数据集。我们将探索不同的数据可视化方式来让你意识到质控后发生了什么。scater包提供了几个有用的简化可视化的函数。
scRNA-seq分析的一个重要方面是控制批次效应。批次效应是实验过程中引入的技术偏差。例如两个不同实验室准备的样品或同一实验室不同时间准备的样品，同时或同地准备的样品相似度更高。最差的情况，批次效应引入的差异可能大于生物本身的偏差，进而获得错误结果 (http://f1000research.com/articles/4-121/v1)。因为样品处理过程记录详细，所以Tung数据集将允许我们探索这一现象。理想上，如果同一个体来源的细胞聚在一起，分组的数目与个体是对应的，就说明存在批次效应。
PCA分析
查看数据整体分布的最简单方式是PCA主成分分析,在其前两个主成分轴查看数据的分布。
主成分分析 (PCA)是一种数学降维方法, 利用正交变换 (orthogonal transformation)把一系列可能线性相关的变量转换为一组线性不相关的新变量，也称为主成分，从而利用新变量在更小的维度下展示数据的特征。
PCA分析不是简单地选取2个或3个变化最大的基因，而是先把原始的变量做一个评估，计算各个变量各自的变异度(方差)和两两变量的相关度(协方差)，得到一个协方差矩阵。在这个协方差矩阵中，对角线的值为每一个变量的方差，其它值为每两个变量的协方差。随后对原变量的协方差矩阵对角化处理，即求解其特征值和特征向量。原变量与特征向量的乘积(对原始变量的线性组合)即为新变量(回顾下线性代数中的矩阵乘法)；新变量的协方差矩阵为对角协方差矩阵且对角线上的方差由大到小排列；然后从新变量中选择信息最丰富也就是方差最大的的前2个或前3个新变量也就是主成分用以可视化。
tSNE聚类
scRNA-seq中另一个常用的可视化方法是tSNE。tSNE (t-Distributed Stochastic Neighbor Embedding) 通过组合降维(如PCA)和最近邻网络随机行走算法在保留细胞的局部距离的基础上实现高维数据(14000维基因表达)到二维空间的映射。与PCA不同，tSNE算法有随机性，每次运行结果都会不同。因为它的非线性和随机性特征，tSNE结果难以直观解释。为了保证重复性，我们固定一个随机数种子使得每次结果都一致。
PCA线性降维算法在集中将不相似的数据点放置在较低维度区域时，尽可能多的保留原始数据的差异，因此导致小部分数据点相距甚远，大部分数据重叠放置。tSNE把点的高维空间的距离转换成点的相似度的概率，维持高维空间和低维空间中一对点之间的条件概率差值总和最小。同时使用t-分布的长尾性解决高维数据映射到低维时的重叠问题。t-SNE算法定义了数据的局部和全局结构之间的软边界，既可以使点在局部分散，又在全局聚集，同时照顾近距离和远距离的点。其性能优于任何非线性降维算法。具体见推文还在用PCA降维？快学学大牛最爱的t-SNE算法吧(附Python/R代码)和B站视频https://www.bilibili.com/video/av19592239?from=search&seid=1152699498338102899。
扫一扫
398
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
