单细胞测序的设计与分析

Posted 2023-04-16

参考技术A 缩略词表：
fluorescence-activated cell sorting (FACS) —— 流式细胞荧光分选技术
whole-genome amplification (WGA)—— 全基因组扩增技术
multiple displacement amplification (MDA) —— 多重置换扩增
multiple annealing and looping-based amplification cycles（MALBAC）—— 多次退火环状循环扩增技术
micro-well displacement amplification system (MIDAS) —— 微孔置换扩增系统
UMI （Unique molecularidentifier)——特异性分子标签（UMI）

由于哺乳动物单细胞DNA含量少于10pg，所以全基因组扩增技术对于单细胞测序和微阵列分析至关重要。
目前WGA有以下几种技术：
PCR、MDA、MALBAC

但是所有方法都会引入基因组覆盖度不均造成的技术伪像。尤其是GC偏差、碱基错配、DNA嵌合体。所以我们应该根据所需的结果选择合适的方法。
例如，基于随机引物PCR的方法可实现高度均匀的扩增，但产量高
仅覆盖基因组的稀疏区域，因此非常适合长度少于长度60 kb的低分辨率拷贝数变异检测。
MDA具有更好的基因组覆盖度，适用于SNP的检测，但是由于MDA复制的高度不均匀性，做CNV检测则有很强的不确定性。
PCR和MDA都会产生嵌合DNA分子会被认为是插入或结构重排。

MDA中嵌合的发生机制
Lasken, R.S., Stockwell, T.B. Mechanism of chimera formation during the Multiple Displacement Amplification reaction. BMC Biotechnol 7, 19 (2007).

MALBAC复制均匀且覆盖度高

先进行预扩增，MALBAC引物随机退火至DNA模板。在高温下具有置换活性的聚合酶会放大模板，生成“半扩增子”。随着扩增和退火过程的重复，半扩增子被扩增为完整的扩增子，其3'端与5'端互补。结果，全扩增子末端杂交形成环状结构，抑制了环状扩增子的进一步扩增，而仅半扩增子和基因组DNA经历了扩增。经过5次的预扩增之后进入PCR流程。最终可获得93%的基因组覆盖度和平均25×的测序深度。
与MDA相比提高了复制的均一性，等位基因缺失率大大降低。 (1% for MALBAC versus 31%–65% for MDA )。MALBAC的假阳性率(4 ×10-5）这是由于聚合酶的忠实度较低，可以多用几个细胞做重复以降低假阳性率。MALBAC适用于同样表征的细胞的SNPs和CNVs检测。

MIDAS相比于MDA所需样本量减少1000倍，增加对模板的忠实度，减少污染。纳米孔反应也有这样的作用，假阳性率降低至4×10-9.

WGA之后的定量可由测序或者基因芯片完成。

首先要比对到参考基因组上，参考基因组可以从UCSC genome browser或者ensembl获得。
再比对之前需要检查reads的质量，切除低质量序列以及接头。长度过短的reads也应该舍弃以避免多重比对。之后比对到基因组上保留只比对到一个位置的reads，对于多重比对的reads有两种处理方法。一是舍弃，二是计算权重（每个reads的权重总和为1）。
对于拷贝数变异的检测，将肿瘤细胞和肺肿瘤细胞的拷贝数用归一化因子标准化之后进行比较，通常使用genome analysis toolkit （GATK）。为增加可信度一般会设置重复，另外细胞周期也会对CNV产生影响，应使用G1或G2/M期细胞，避免使用S期细胞。

单细胞测序面临的主要问题就是获得的遗传物质的量很少，上面我们已经介绍了扩增方法，但是这些方法都存在扩增偏差，这会使不同基因的mRNA的丰度受到影响。
在过去，扩增的单细胞RNA用微阵列芯片检测（2002）。不过目前已经发展出一些灵敏度较高的单细胞测序技术，第一个单细胞测序方案发布于2009 Surani的实验室。最初的扩增方法是利用带有特殊锚定序列的poly(T)引物捕获以及反转录poly（A）RNA，获得的单链cDNA经过多聚腺苷酸化再结合带有特殊锚定序列的poly(T)引物，得到双链cDNA。cDNA通过针对锚定序列的引物进行PCR扩增，在建立文库之前将产物片段化。

A. CEL-seq
多聚腺苷酸mRNA被oligo dT引物反转录，该引物含有Illumina P1 接头，细胞条形码，T7启动子，通常还会含有一个UMI。紧接着进行第二条链合成，从T7启动子开始，之后双链cDNA进行片段化，连接上含有Illumina P2接头。对reads的测序开始于mRNA的3‘端。
B. STRT-seq
使用Oligo-dT引物对多聚腺苷酸RNA逆转录，该引物还包含Illumina P1 接头和Pvul 限制性酶切位点。用一段带有Illumina P1 接头、UMI和template switch oligo (TSO)的引物接在转录本的5’端，然后合成双链cDNA。双链cDNA通过与Illumina P1 接头互补的引物扩增，片段化产物，用Tn5转座酶在片段上链接Illumina P2接头以及细胞条形码。3‘端被Pvul限制酶消化，仅保留5’端进行测序。
C. Smart-seq2
使用含有Oligo-dT的PCR引物对多聚腺苷酸RNA逆转录，同样的引物作为TSO的一部分被添加在模板链5‘末端。PCR扩增后，产物片段化，由Tn5转座酶在3’和5‘末端分别加上不同的引物。新一轮的扩增使用Nextera sequencing primers ，这样可以实现全长覆盖，但是没有UMI定量。

此外还有：液滴测序（Drop-seq），可以快速低成本的测多个细胞，并且多个细胞同时出现在一个液滴中也避免了上样量过低对测序造成的影响。
http://www.merrybio.com.cn/blog/Drop-seq.html
https://www.jianshu.com/p/0800a07cfa37

为了量化敏感度，我们通常会采用已知浓度的外源spike-in RNA
https://www.youtube.com/watch?v=YVlrzKMJ2uc
加入spike-in的浓度通常为mRNA总数的1%~5%，通常会使用ERCC的产品，这些涉及的RNA比哺乳动物的RNA短，有较短的poly(A)尾，缺乏5’ 帽。

分析的第一步就是进行质控（fastqc）并修剪(bwa)reads，对于人和鼠来说最终应保留长度>35bp的reads。
在比对到参考基因组之前，应确保barcode\UMI等primer来源的序列都被除去。不过对于1对reads来说，其中一条read保留索引信息，另一条read比对到参考基因组上[见前文图：三种测序方法]。通常，可以将读段映射到基因组，然后通过将基因组的读段与基因模型注释相交来进行表达定量。建议仅保留单一比对的reads。

由于单细胞测序对基因的覆盖度低，不同转录本的鉴定（Cufflinks）成为一个难题。如果异构体的信息对于你的研究不是必须的，你可以把这些异构体合并到同一个基因位点。
除了依照参考文献，更重要的是考虑实验策略。如果我们的测序方法回富集3‘或5’端的序列，那么基因注释的质量就会对实验的里灵敏度产生很大的影响。因为基因模型在转录本的两端可信度较低，改善3‘或5’端注释可能会更好，尤其是对于那些非标准的模式生物。例如，Junker等人运用一种修正的CEL-seq进行长读段低深度测序以精确检测斑马鱼胚胎的3’poly(A)位点。

一旦细胞中的所有的reads或者转录本被计数，我们建议滤除reads含量低的细胞。这可能是样品准备过程造成的问题，比如细胞凋亡、应激、不当裂解、RNA降解或者扩增测序的效率较低。每个细胞中reads的总数或者UMIs代表的转录本的数量，最先预示着样本的质量。应该设置阈值以去除read counts分布左尾的细胞，防范由低质量细胞产生的伪像。
spike-in RNA 的表达可以用来鉴别和剔除测序效率不高的样本。由于所有样品的spike-in RNA数量应相同，因此鉴定低产量样品非常简单。

单细胞测序数据的差异表达分析方法总结

无论是传统的多细胞转录组测序（bulk RNA-seq）还是单细胞转录组测序（scRNA-seq），差异表达分析（differential expression analysis）是比较两组不同样本基因表达异同的基本方法，可获得一组样本相对于另一组样本表达显著上调（up-regulated）和下调的基因（down-regulated），从而可进一步研究这些差异表达基因的功能，包括富集的通路（pathway）或生物学过程（biological process）。

 

由于单细胞测序技术的局限性，单细胞测序数据通常具有高噪音，有较高的dropout问题，即很多低表达或中度表达的基因无法有效检测到。所以，以前针对传统多细胞转录组测序数据开发的差异表达检测方法或软件不一定完全适用于单细胞测序数据。若想比较不同细胞亚型或不同条件下的细胞表达差异时，为了能得到可靠的结果，需要选定一个好的差异表达分析方法（微信公众号：AIPuFuBio）。

 

近年来，有不少专门针对单细胞转录组测序数据的差异表达分析方法相继被开发出来，如MAST (Finak et al., 2015)、SCDE (Kharchenko et al., 2014)、 DEsingle (Miao et al., 2018)、 Census (Qiu et al., 2017)、 BCseq (Chen and Zheng, 2018)等。具体可以见下表所示：

红线上方是专门针对单细胞测序数据开发的差异表达分析软件或R包，红色下方是针对bulk转录组数据开发的软件或R包

 

 

图1、一些比较流行的差异表达分析软件（Chen et al. Frontiers in Genetics, 2019） 

这里要值得提一下SCDE（全名：Single Cell Differential Expression）软件，其属于最早一批专门针对单细胞测序数据开发的差异表达分析软件，地址为：https://hms-dbmi.github.io/scde/。下图是原文章中SCDE与其他传统差异表达分析软件的性能比较，显示SCDE具有不错的性能。

 

图2、SCDE与其他软件在单细胞测序数据集上鉴定差异表达基因的性能比较（Kharchenko et al. Nature Methods, 2014）


最近，Wang et al.等人比较了11款经典的软件在单细胞测序测序数据集上的差异表达分析性能，这些软件具体如下表所示：

图3、不同差异表达软件的相关信息（Wang et al. BMC Bioinformatics, 2019）

图4、不同差异表达软件ROC曲线比较（ Wang et al. BMC Bioinformatics, 2019）

图5、不同差异表达软件各主要指标的比较（ Wang et al. BMC Bioinformatics, 2019）

 

图6、不同差异表达软件之间在真实数据集上检测到的差异表达基因比较（ Wang et al. BMC Bioinformatics, 2019）。差异表达基因的定义为：adjusted p-value< 0.05

 

图7、样本数量对不同差异表达软件各方面性能的影响比较（ Wang et al. BMC Bioinformatics, 2019）

 

图8、不同差异表达软件鉴定到的top 300个差异表达基因富集的显著KEGG通路数和GO条目数比较（ Wang et al. BMC Bioinformatics, 2019） 。（条件：FDR<0.05）


总的来说，不同的差异表达软件有不同的优缺点。有些软件具有高灵敏性，但检测精度却比较低，有些则刚好相反。这11款软件中，DEsingle 和SigEMD这两个方法较好的平衡了差异表达基因检测灵敏性和准确性。值得注意的是，Wang et al. 的比较发现，现在这些专门针对单细胞测序数据开发的差异表达分析软件和传统的方法相比，并没有显示出太多的优势（ Wang et al. BMC Bioinformatics, 2019）。这也进一步说明，还需不断开发新的单细胞测序差异表达分析方法，以更好的检测单细胞测序数据的差异表达基因。（更多经典，可见大型免费综合生物信息学资源和工具平台AIPuFu：www.aipufu.com）。笔者建议，做单细胞测序数据的差异表达分析，最好还是选择专门针对单细胞测序数据开发的软件，如SCDE、DEsingle 和SigEMD等。

希望今天的内容对大家有用哦，会持续更新的，欢迎留言~~


参考文献

1. Chen et al. Single-Cell RNA-Seq Technologies and Related Computational Data Analysis，Frontiers in Genetics, 2019

2. Wang et al. Comparative analysis of differential gene expression analysis tools for single-cell RNA sequencing data, BMC Bioinformatics, 2019

3. Kharchenko et al. Bayesian approach to single-cell differential expression analysis, Nature Methods, 2014