少量细胞输入的自制转录组测序文库构建试剂评测

丁 蕾 ，高彩霞 ，刘兆远 ，陈 磊

1.上海交通大学基础医学院免疫学与微生物学系，上海 200025；2.上海交通大学医学院，上海市免疫学研究所，上海 200025

伴随着高通量测序技术的迅猛发展和测序成本的降低，转录组测序借助其灵敏度高、检测范围广、提供转录组信息全面等独特的优势在转录组学研究领域逐渐占据主导地位[1]，是目前深入研究转录组复杂性的强大工具。

针对转录组测序的生物信息学分析主要包括以下几个步骤：数据质量控制、序列比对和基因表达分析。其中数据质量控制主要包括评估建库试剂在实验中捕获转录本的能力和所测序列在基因覆盖度的均匀性。基因表达分析主要包括差异表达基因的鉴定和差异表达基因富集分析，如KEGG（Kyoto encyclopedia of genes and genomes）通路富集分析，这是转录组测序最重要的分析结果。基因表达分析结果可用于更深入的研究，包括转录本结构研究（如基因的可变剪切）、转录本变异研究（如基因融合、单核苷酸突变）、差异基因表达水平的比较，甚至包括全新转录本或稀有转录本的发现[2]。尤其对于如肿瘤干细胞等来源极为有限的生物样品，转录组测序技术的优势明显[3]，要求的起始样品量要比芯片技术少得多且技术重复性好。SMART（switching mechanism at 5′ end of RNA template）技术经优化后，选择使用锁核酸（locked nucleic acid，LNA）、更高浓度的氯化镁以及甜菜碱的新方案，使得起始仅使用1～1 000个细胞或10 pg～10 ng的总RNA就可获得高质量的测序文库，而且还使得序列在基因的覆盖度上获得很大改善[4]，能够更好地实现单核苷酸变异等的检测。基于该技术的商业化试剂盒有多种，均价格昂贵，其中普遍使用的是TaKaRa Bio公司名为SMART-Seq v4 Ultra Low Input RNA Kit的试剂盒（下文简称为TaKaRa试剂）。本研究利用SMART技术自制相对价格更低的建库试剂（下文简称为DIY试剂）用于少量细胞输入的转录组测序文库构建实验，并通过生物信息学分析，从多个方面比较DIY试剂和TaKaRa试剂对转录组测序结果的影响，从而验证用DIY试剂替代TaKaRa试剂的可行性，以降低转录组测序文库构建的成本。

1 材料与方法

1.1 材料

1.1.1 实验细胞和动物 小鼠腹腔巨噬细胞，分别来自于4只8周龄雌性SPF级C57BL/6小鼠。小鼠饲养于上海交通大学医学院实验动物科学部屏障环境内。生产许可证号为SCXK（沪）2018-0007，使用许可证号为SYXK（沪）2018-0027。

1.1.2 主要试剂 TaKaRa试剂：SMART-Seq v4 Ultra Low Input RNA Kit，购于TaKaRa Bio公司。DIY试剂：不含钙、镁的磷酸缓冲液DPBS（Dulbecco′s phosphate buffered saline），1×TrypLE Express 酶（无酚红），购于Gibco公司；RNaseZap，购于Ambion公司；DNA去除试剂（DNAOFF），购于TaKaRa Bio公司；2′-脱氧核苷酸-5′-三磷酸混合物（10 mmol/L dNTP mix），购 于Fermentas公 司；5×First-strand Buffer（250 mmol/L Tris-HCl，pH 8.3，室温；375 mmol/L 氯化钾；15 mmol/L 氯化镁；二硫苏糖醇）、Superscript Ⅱ Reverse Transcriptase，购于Invitrogen公司；重组RNase 抑制剂，购于Clonetech公司；甜菜碱、氯化镁（无水），购于Sigma-Aldrich公司。其他试剂：蒸馏水，购于Gibco公司；2×KAPA HiFi HotStart ReadyMix，购于KAPA Biosystems公司；Agencourt Ampure XP beads，购于Beckman Coulter公司；溴化乙锭（ethidium bromide，EB，10 mmol/L Tris-HCl，pH 8.5），购于Qiagen公司；双端测序TruSeq Dual-index Sequencing Primer试剂盒、Nextera XT DNA Sample Preparation试剂盒、Nextera XT 24-index试剂盒、Adapter oligos，购于Illumina公司；99.5%乙醇，购于Sigma-Aldrich公司；TRIzol 试剂、Glycoblue，购于Invitrogen公司；氯仿、异丙醇，购于国药试剂公司；4%巯基乙酸盐肉汤，购于BD公司。

1.1.3 主要仪器 Agilent 2100 生物分析仪系统，购于安捷伦科技有限公司；Illumina NextSeq 500仪器，购于Illumina公司；移液枪、低温超速离心机、高速台式离心机，购于Eppendorf公司；EasyCycler 96 PCR仪，购于Analytik Jena AG公司。

1.1.4 主要软件 使用基于JAVA 8 的FastQC 0.11.9对高通量测序数据进行质量评估；使用基于Python 3.7 的RSeQC 3.0.1 评估高通量测序尤其是转录组测序数据的覆盖均匀性；使用基于C++ 的fastp 0.20.0对原始测序数据进行随机抽样等操作；使用HISAT2 （hierarchical indexing for spliced alignment of transcripts 2） 2.1.0进行有参考基因组的序列比对；使用基于C 语言的SAMtools 1.4处理经HISAT2 等软件比对后生成的sam/bam 文件；基于Python 3.7 的HTSeq 0.11.3处理高通量测序数据；使用基于R 3.6.1的DESeq2、edgeR、ggplot2软件包进行基因差异分析以及数据可视化分析等。

1.2 实验方法

1.2.1 腹腔巨噬细胞的诱导 将4只小鼠随机均分为2组，再分别编号为M1、M2、M3、M4，其中M1、M2号小鼠作为对照组不进行任何处理，M3、M4号小鼠作为实验组向其腹腔内注射1 mL浓度为4%的巯基乙酸盐肉汤以诱导巨噬细胞[5]。经72 h后，利用流式细胞分选技术分离腹腔巨噬细胞。

1.2.2 腹腔巨噬细胞RNA的提取 腹腔巨噬细胞离心5 min，4 ℃，500×g；向每个样品中加入500 μL TRIzol，混合均匀，室温孵育5 min；置于-80 ℃至少1 h；向每个样品中添加适量的氯仿（TRIzol与氯仿的体积比为5:1），剧烈摇动约15 s，室温孵育2～3 min；离心15 min，4 ℃，12 000×g；将水相转移至新管中（TRIzol与氯仿的体积比为2:1），此后在冰上工作；向每个样品中添加1 μL共沉淀剂GlycoBlue；在每个样品中加入适量的异丙醇（TRIzol与异丙醇的体积比为2:1），混合均匀，-80 ℃孵育1 h以上；离心20 min，4 ℃，12 000×g；用适量的75%预冷过的乙醇洗涤（TRIzol与异丙醇的体积比为1:1）；离心15 min，4 ℃，7 400×g，风干以尽可能地除去乙醇；在12.5～15.0 μL无核酸酶的水中将每个RNA沉淀重悬；室温孵育2～3 min，充分混合并置于冰上；使用Qubit仪器定量RNA；保存RNA于-80℃条件下。

1.2.3 cDNA的构建及二代测序 将编号为M2和M3的小鼠的腹腔巨噬细胞提取的RNA各均分成2份，用DIY试剂和TaKaRa试剂分别进行cDNA的构建，所得4个样品的编号为C_M2_DIY、C_M2_TaKaRa、T_M3_DIY和T_M3_TaKaRa（C，对照组；T，实验组）。另外，从编号为M1和M4的小鼠的腹腔巨噬细胞提取的RNA皆用DIY试剂进行cDNA的构建，所得2个样品的编号为C_M1_DIY、T_M4_DIY。使用Agilent 2100 生物分析仪系统对6个样品的cDNA进行质检。在完成转录组测序文库构建后，对上述6个样品利用Illumina NextSeq 500平台进行双端测序，测序模式为2×75 bp。

1.2.4 测序数据质量控制和序列比对 使用FastQC软件对测序数据进行初步质量控制，输出结果会给每个碱基一个相应的质量评分，用于衡量测序精确度。一个给定碱基的测序质量评分Q定义为：Q=-10×lge，其中e为预计碱基检出不正确的概率。故碱基的质量评分与碱基检出精确度存在一定关系，即Q较高表示出错的概率较小。其中，质量控制结果中的重要指标Q30（%）代表碱基质量评分≥30的碱基的数量占全部碱基数的百分比。一般而言，Q30＞85%时，测序数据质量合格。使用fastp软件对原始的序列文件进行随机抽样，使用HISAT2软件[6]进行有参考基因组的序列比对，并使用SAMtools软件[7]处理比对结果，再用HTSeq软件[8]统计比对到参考基因组上基因区间内的读序数目，计算不同测序数据量下可检测的转录本数，比较DIY试剂和TaKaRa试剂捕获转录本的能力。使用RSeQC软件[9]对测序数据进行质量控制，获得测序数据在基因上的覆盖均匀性。从数据质量以及序列比对的层面探究2种试剂对转录组测序结果的影响。

1.2.5 基因表达分析 经过HTSeq软件完成基因表达量的计算后，使用R语言软件包DESeq2[10]和edgeR[11]进行样品的相似性分析以及实验组与对照组之间的差异表达基因分析[12]，选择符合|log2fold change|＞1、P＜0.005的基因为差异表达基因，探究2种试剂处理样品的差异表达基因鉴定结果是否一致。log2fold change表示以2为底、基因表达量差异的对数，代表实验组中的基因表达量与对照组中的基因表达量的差异倍数。若log2fold change的数值为正，则反映该基因在实验组中的表达量比在对照组中的表达量高；相反若数值为负，则反映该基因在实验组中的表达量比在对照组中的表达量低。

1.2.6 差异表达基因通路富集分析 使用在线富集分析工具DAVID（database for annotation, visualization, and integrated discovery）6.8[13]对差异表达基因进行通路富集分析，并使用R语言软件包ggplot2[14]进行可视化分析。通路富集分析方法以通路为单位，以物种已知的全部基因为背景，通过Fisher精确检验来分析计算各个通路基因富集度的显著性水平，从而确定受到显著影响的代谢和信号转导途径，探究2种试剂处理情况下差异表达基因通路富集分析的结果是否一致。以错误发现率（false discovery rate，FDR）表示富集的显著性水平。FDR表示以10为底、P值的对数，FDR数值越大，代表在统计学上该通路的影响或变化越显著。

2 结果

2.1 2种试剂处理样品的数据质量比较

6个样品的cDNA片段均分布在200～3 000 bp范围内，主峰均在1 000 bp左右，cDNA质量良好。使用FastQC软件对测序原始数据进行质量控制，结果如表1所示。Q30均＞91%，表明数据质量良好。经HISAT2软件比对后，发现比对到基因组的序列占比均＞90%，捕获的转录本数也基本相近。

表1 样品基本信息及数据质量控制Tab 1 Sample basic information and data quality control

2.2 2种试剂捕获转录本的能力比较

在一定范围内，转录组测序结果显示捕获到的转录本数会随着测序数据量的增加而增加，并最终趋于饱和。为了比较DIY试剂和TaKaRa试剂捕获转录本的能力，使用fastp软件对原始的序列文件进行随机抽样，获得具有0.5、1、1.5、2、2.5、3、4、5、6、7 M的序列子集文件，再用HTSeq软件计算不同测序数据量下可检测的转录本数。如图1所示：可检测的转录本数随着测序数据量的增加而增加；测序数据量一定时，对照组（C_M1_DIY、C_M2_DIY、C_M2_TaKaRa）和 实 验 组（T_M3_DIY、T_M4_DIY、T_M3_TaKaRa）可检测的转录本数差异明显，而不同试剂处理的样品各自在对照组（C_M2_DIYvsC_M2_TaKaRa）和实验组（T_M3_DIYvsT_M3_TaKaRa）内无明显差异。

图1 2种试剂捕获转录本的能力比较Fig 1 Comparison of the capability in transcripts capture of DIY reagent and TaKaRa reagent

2.3 2种试剂处理的样品基因覆盖度比较

使用RSeQC软件计算所测序列在基因上的覆盖度，结果如图2所示。DIY试剂和TaKaRa试剂处理的样品基因覆盖度都较均匀、无偏差，且两者的一致性较高。

图2 2种试剂处理的样品基因覆盖度Fig 2 Gene coverage plot of the samples treated with the two reagents

2.4 主成分分析

分别使用R软件的DESeq2包和ggplot2包对样品进行主成分分析（principal component analysis，PCA）以及可视化分析。图3清晰地体现了样品的聚类情况，横坐标PC1反映实验处理导致的样品间差异，即巯基乙酸盐肉汤对腹腔巨噬细胞持续72 h的诱导作用，存在91%的差异性；纵坐标PC2反映2种试剂导致的样本间差异，该差异性仅为8%，说明DIY试剂和TaKaRa试剂对样品的影响远小于实验处理导致的基因表达谱差异。而且，DIY试剂处理的对照组样品（C_M1_DIYvsC_M2_DIY）和实验组样品（T_M3_DIYvsT_M4_DIY）的组内差异均较小，表示DIY试剂的实验重复性好且结果稳定。

图3 PCA体现的样品聚类情况Fig 3 Sample clustering reflected by PCA

2.5 2种试剂处理样品的差异表达基因鉴定结果

将6个样品分为2组：第1组为DIY试剂处理的样品，包括C_M1_DIY、C_M2_DIY、T_M3_DIY、T_M4_DIY；第2组为TaKaRa试剂处理的样品，包括C_M2_TaKaRa和T_M3_TaKaRa。使用edgeR包对2组数据就对照组和实验组之间的差异表达基因进行鉴定，最终DIY试剂处理组样品共有3 877个差异表达的基因，TaKaRa试剂处理组样品共有3 855个差异表达的基因，重合的基因数为2 737个。结果如图4所示，从2组结果中各取显著性差异最大的前1 000个基因（按照P值排序），重合的基因数有748个，重合率达75%，并且相同基因在2组数据中的|log2fold change|值比较接近。

图4 2种试剂处理样品的差异表达基因鉴定结果Fig 4 Differential gene expression analysis results of the samples treated with the two reagents

2.6 2种试剂处理样品的基因通路富集分析

根据上述按照所用试剂将样品分为2组的策略，各自选取差异表达分析结果中最显著的1 000个基因（按照P值排序），使用DAVID在线分析工具进行KEGG通路富集分析，结果显示2组差异基因中富集程度显著的通路一致性较高，重合率达80%。使用R语言的ggplot2软件包对各自重合的前10个通路进行可视化，结果如图5所示。横坐标数值为正数代表在实验组中上调基因富集出来的结果，为负值则代表在实验组中下调基因富集出来的结果；数值的大小则体现通路富集的显著性水平。结果表明，DIY试剂处理组和TaKaRa试剂处理组的基因富集通路的显著差异情况较相近。

图5 2种试剂处理样品的基因通路富集分析结果Fig 5 Results of KEGG pathway analysis of the samples treated with the two reagents

3 讨论

转录组测序的研究对象为特定细胞在某一功能状态下能转录出来的所有RNA的总和，可以用于基因功能以及基因结构的研究[15]。转录组测序结果可以揭示特定生物学过程以及疾病发生过程中的分子机制，因此该技术被广泛应用于生物学、基础医学、临床诊断和药物研发等多个领域[16]，并成为一个识别细胞类型、鉴定标记基因、识别信号通路和研究调控机制的重要工具。转录组测序技术可为不同个体和疾病的组织带来更好的生物靶标识别和药物标靶[17]，从而有助于进一步发展精准医学。

本研究使用DIY试剂和TaKaRa试剂分别进行转录组测序文库构建，并通过生物信息学分析，从不同的方面来验证DIY试剂代替TaKaRa试剂的可行性。结果发现，2种试剂处理样品得到的数据质量皆良好，比对到基因组的序列占比也基本相等；在相同的测序深度下，2种试剂捕获的转录本数较相近；2种试剂处理的样品所测序列在基因上的覆盖度均匀，且一致性较高。此外，分析结果显示DIY试剂在仅有少量细胞输入的情况下，能够捕获相对较多的转录本数量，表明该试剂使得来源极为有限的生物样品分析成为可能。所以从数据质量、捕获转录本能力、基因覆盖度均匀性等方面来看，DIY试剂可以很好地替代TaKaRa试剂进行少量细胞输入的转录组测序文库构建。差异表达基因的鉴定和差异基因通路富集分析是转录组测序较重要的分析结果[18]，所以有必要比较这2个分析结果在2种不同试剂处理的情况下是否一致。结果显示，2组分析结果相近，无明显差别。从构建文库的成本方面来看，DIY试剂的优势更加明显：若用于建库的起始细胞量为100个细胞，TaKaRa试剂建库实验所花费用高达6万元，DIY试剂则仅需约5 700元。然而在进行一项科学研究时，动辄需要对上千万乃至更大数量级的细胞进行文库构建及测序，测序的费用更是昂贵，所以使用DIY试剂代替此类商业化试剂不仅能在保证数据质量良好、差异基因分析等结果可靠的前提下降低约90%的建库成本，还能根据具体实验进行酶反应条件的调整以及试剂使用的优化等。因此综合多方面的评测，DIY试剂可以替代昂贵的商业化试剂用于少量细胞输入的转录组测序文库构建。此外，SMART技术不仅可用于少量细胞输入的转录组测序文库构建，也可对单个细胞进行建库并测序，用以在免疫、神经等复杂系统中研究细胞高异质性这一类问题；该技术不仅为测序技术提供了新的发展方向，还有望在单细胞水平上进行遗传变异的检测、肿瘤的诊断及免疫治疗等研究[19]。因此，如何优化DIY试剂并将其用于单细胞的转录组测序文库构建还亟待解决。

综上所述，本研究比较了DIY试剂与TaKaRa试剂对少量细胞的转录组测序文库构建以及转录组测序结果，表明DIY试剂可以替代昂贵的商业化试剂进行转录组测序文库构建。该研究将使得基于SMART技术的转录组测序得到真正的广泛应用，将为生物学以及基础医学等各个研究领域提供良好的技术支持。