一种单细胞RNA测序的制作方法

一种单细胞rna测序
技术领域
1.本发明属于分子生物学和药物开发领域。具体来说，本发明涉及哺乳动物细胞培养、单细胞rna测序(scrna-seq)和药物筛选。


背景技术：

2.药物筛选可以从转录组(rna)、蛋白质组以及代谢组的数据中获得。在用药物治疗前后，通过比较细胞或组织的基因表达谱获得转录组特征。用先前的批量方法开发的药物遇到一些问题，其中对某些病人有效的药物可能对其他病人效果很小，甚至没有效果，或者药物可能在某些病人的治疗过程中逐渐失去疗效，尽管它在治疗初期似乎是有效的。以上问题可部分归因于不同水平的异质性，也就是说，在不同个体或一个个体中的不同肿瘤以及一个肿瘤中的不同细胞之中，其表达谱可能不同，从而对药物的反应也不同。因此，药物可能只对一些而非所有肿瘤细胞(在此表示为亚群)起作用，而当停止治疗时，其他亚群仍存活并恢复繁荣。还有一种可能性是，在用某些药物治疗期间，新的肿瘤细胞亚群甚至可能进化。这造成大量通过批量方法发现的类似药物的化学物质在临床试验中(甚至之前)终止，导致药物开发中时间和经济成本的巨大浪费。
3.肿瘤细胞的异质性是导致肿瘤复发和/或对治疗产生耐药的原因之一。scrna-seq能够揭示耐药亚群的存在和表达特征，以解释和支持由细胞异质性导致的药物失效。scrna-seq已经广泛用于肿瘤异质性的转化研究。然而，scrna-seq分析的能力还未被用于高通量药物筛选。
4.培养的原代肿瘤细胞或肿瘤源性细胞系多年来一直被用作药物筛选的模型，因为与基因工程小鼠或异种移植等其他模型相比，其具有便捷性、相对较低的成本及物理相关性。
5.基因表达谱(例如微阵列或rna测序)已被应用于药物开发中，并推动药物开发从常用的基于靶点的方法进入到基于分子表型的阶段，后者不需要已知的分子机制。
6.利用转录组的变化作为对细胞系中候选药物反应的度量，积累了大量的数据，这些数据构成了连接库，将药物与疾病关联起来用于药物开发，并比较疾病与药物反应特征以用于药物再利用。
7.用于癌症治疗的药物是为杀死癌细胞而设计的。由于药物浓度迅速下降，细胞易处于非致死剂量的环境中。当它们处于非致死范围内时，细胞可以选择保持增殖或变得衰老。然而，批量rna转录组不能区分细胞的发育方向。
8.scrna-seq能够揭示耐药亚群的存在和表达特征，以解释和支持由细胞异质性导致的药物失效。并且可以筛选出针对不同细胞亚型的有效药物，或者发现针对不同细胞亚型的两种或多种药物的组合。
9.在癌症治疗期间，一些癌细胞不能被药物杀死并变得衰老。一方面，scrna-seq能够筛选出诱导衰老的药物，作为一种有可能降低细胞毒性并使免疫系统参与其中的治疗策略。另一方面，衰老表型的缓慢发展为改进联合用药等药物治疗策略赢得了挑战。
10.单细胞水平的测序能够揭示细胞的异质性，这种异质性在批量rna测序中，是被检测到的全部细胞的平均表达所掩盖的。单细胞基因表达分析的优点包括：(1)更准确地表征单个细胞的基因表达谱，和(2)识别细胞类型，包括罕见的新细胞类型或亚型，(3)获得细胞从一种状态过渡到另一种状态期间的基因表达快照(snapshot)。
11.本发明涉及以更通用和系统的方式将scrna-seq应用于细胞系中的药物筛选和再利用，作为解决由肿瘤异质性导致的药效丧失问题的平台。


技术实现要素：

12.本发明涉及一种用于基于细胞的药物筛选的分析方法，包括以下步骤：
[0013]-在培养环境中培养细胞群；
[0014]-将细胞群暴露于药物物质中，其中设置了明确的时间监测点和药物物质的浓度监测点；
[0015]-暴露于药物物质中后，形成单细胞文库；
[0016]-对细胞文库中的细胞的多个基因进行单细胞文库的单细胞rna测序；以及
[0017]-进行生物信息学分析，以确定差异表达基因、富集通路和发育轨迹。
[0018]
本发明涉及一种用于基于细胞的药物筛选的分析方法的应用，其中差异表达基因是潜在的标记基因，并用于筛选癌症药物。
具体实施方式
[0019]
本文所用的术语“亚群”是指通过生物信息学算法对单细胞rna测序进行聚类，每个聚类(cluster)都是一个亚群。
[0020]
本文所用的术语“聚类”是指根据单细胞rna测序数据之间的相似性，将数据归类为相应的类别。同类别之间有高度相似性，不同类别之间的差异程度最大。
[0021]
本文所用的术语“聚类分析”是指将单细胞rna测序数据的集合分组为由相似基因表达谱构成的多个类别的过程。
[0022]
本文所用的术语“生物信息学分析”是指确定差异表达基因、聚类分析、富集通路和发育轨迹，并指本专利中的所有分析。
[0023]
本文所用的术语“差异表达基因(deg)”是指其表达在细胞亚群之间有显著差异的基因。
[0024]
本文所用的术语“再利用”是指在另一种疾病中发现药物的新用途。
[0025]
本文所用的术语“促调亡基因”是指促进细胞调亡的基因。
[0026]
本文所用的术语“耐药性癌症”是指对药物治疗有耐受性的癌症。
[0027]
本文所用的术语“发育轨迹”是指根据不同的细胞状态来模拟细胞分化过程。
[0028]
本文所用的术语“富集通路”是指通过差异表达基因来富集kegg通路。
[0029]
本文所用的术语“癌症”是指但不限于：胆道癌；膀胱癌；脑癌，包括胶质母细胞瘤和髓母细胞瘤；乳腺癌；宫颈癌；绒毛膜癌；结肠癌；子宫内膜癌；食道癌；胃癌；血液肿瘤，包括急性淋巴细胞性和骨髓性白血病；多发性骨髓瘤；艾滋病相关白血病和成人t细胞白血病淋巴瘤；上皮内肿瘤，包括鲍恩病(bowen’s disease)和佩吉特氏病(paget’s disease)；肝癌；肺癌；淋巴瘤，包括霍奇金病(hodgkin’s disease)和淋巴细胞性淋巴瘤；神经母细胞
瘤；口腔癌，包括鳞状细胞癌；卵巢癌，包括源自上皮细胞、基质细胞、生殖细胞和间质细胞的卵巢癌；胰腺癌；前列腺癌；直肠癌；肉瘤，包括平滑肌肉瘤、横纹肌肉瘤、脂肪肉瘤、纤维肉瘤和骨肉瘤；皮肤癌，包括黑色素瘤、卡波西肉瘤(kaposi’s sarcoma)、基底细胞癌和鳞状细胞癌；睾丸癌，包括生殖肿瘤，如精原细胞瘤、非精原细胞瘤(畸胎瘤、绒毛膜癌)、间质瘤和生殖细胞瘤；甲状腺癌，包括甲状腺腺癌和髓样癌；以及包括腺癌和肾母细胞瘤在内的肾癌、乳腺癌、前列腺癌、肺癌、卵巢癌、结直肠癌和脑癌。
[0030]
本文所表述的术语“癌症治疗”可包括但不限于：化疗、放疗、辅助治疗或上述方法的任意组合。可以变化的治疗方面包括但不限于：剂量、给药时间或持续时间；并且可以结合或者不结合其他治疗，这些治疗也可以在剂量、时间或持续时间上有所不同。
[0031]
用于药物筛选的单细胞rna测序技术。首先对药物处理过的单细胞进行测序，然后对测序数据进行生物信息学分析。scrna-seq能够揭示耐药亚群的存在和表达特征，以解释和支持由细胞异质性所导致的药物失效。并且可以筛选出针对不同细胞亚型的有效药物，或者能够发现针对不同细胞亚型的两种或多种药物的组合。
[0032]
发明人利用单细胞rna测序技术，对肿瘤细胞系k562与奥沙利铂(oxaliplatin)和lsd1等药物反应后的细胞群组成和基因表达的变化进行分析，以选择是否将这些药物用于肿瘤治疗。本发明的区别和技术特征是：(1)虽然已有文献通过rt-qpcr、ngs和流式细胞术分析k562与药物奥沙利铂或lsd1反应的基因表达，但没有通过单细胞rna测序来进行的；(2)细胞系与药物反应的时间和浓度是本专利的发明；(3)在本专利中，我们发现奥沙利铂或lsd1可用于白血病治疗；以及(4)本发明涉及将scrna-seq应用在细胞系中进行药物筛选和再利用。
[0033]
发明人发现，细胞系与药物相互作用一定的时间和浓度后可用于白血病的治疗，特别是奥沙利铂或lsd1可用于治疗。本发明涉及将scrna-seq应用在细胞系中进行药物筛选和再利用。
[0034]
本发明涉及一种用于基于细胞的药物筛选的分析方法，包括以下步骤：
[0035]-在培养环境中培养细胞群；
[0036]-将细胞群暴露于药物物质中，其中设置明确的时间监测点和药物物质的浓度监测点；
[0037]-暴露于药物物质后，形成单细胞文库；
[0038]-对细胞文库中的细胞的多个基因进行单细胞文库的单细胞rna测序；以及
[0039]-进行生物信息学分析，以确定差异表达基因、富集通路和发育轨迹。
[0040]
在一个实施方案中，细胞群选自原代细胞、细胞系、癌细胞。
[0041]
在一个实施方案中，细胞群是k562细胞系。培养环境是指dmem细胞培养基、10％fbs(胎牛血清)和1％ps(青霉素和链霉素)。
[0042]
在一个实施方案中，对细胞文库中单独分选的细胞的多个基因进行单细胞文库的单细胞rna测序。
[0043]
在一个实施方案中，药物物质可以是任何其他抗肿瘤药物或具有抗肿瘤潜力的药物或以前未用于抗肿瘤目的的药物。
[0044]
在一个实施方案中，药物物质是奥沙利铂或赖氨酸特异性组蛋白去甲基化酶1a(lsd1)抑制剂。
[0045]
在一个实施方案中，时间监测点包括0至6小时之间的时间点。在优选的实施方案中，时间监测点包括0至24小时之间的时间点。
[0046]
在一个实施方案中，药物物质的浓度监测点包括50nm和5μm。
[0047]
在一个实施方案中，单细胞rna测序包括分析polya+mrna。
[0048]
在一个实施方案中，单细胞rna测序包括分析mrna。
[0049]
在一个实施方案中，所确定的差异表达基因包括参与药物代谢、细胞生长、增殖和分化通路的细胞亚型的标记基因。
[0050]
在一个实施方案中，所确定的差异表达基因包括bnip3、bnip3l、hemgn和kdm1a。
[0051]
在一个实施方案中，所确定的差异表达基因包括任何在细胞亚群之间差异表达的基因。数值为avg_logfc》0，该亚群中表达该基因的细胞比例显著高于其他亚群。
[0052]
在一个实施方案中，用奥沙利铂和/或任何抗癌药物处理的细胞群中，所确定的差异表达基因是与rna转运和剪接体有关的基因。
[0053]
在一个实施方案中，对细胞文库中单独分选的细胞的多个基因进行单细胞文库的单细胞rna测序。
[0054]
在一个实施方案中，用组蛋白赖氨酸特异性去甲基化酶1a(lsd1)抑制剂处理的细胞群中，所确定的差异表达基因是促调亡基因或对组蛋白赖氨酸特异性去甲基化酶1a(lsd1)抑制剂有抗性的基因。
[0055]
在一个实施方案中，用组蛋白赖氨酸特异性去甲基化酶1a(lsd1)抑制剂处理的细胞群中，所确定的差异表达基因包括bnip3、bnip3l、hemgn和kdm1a基因。
[0056]
本发明涉及一种用于基于细胞的药物筛选的分析方法的应用，其中差异表达基因是潜在的标记基因，并用于筛选癌症药物。
[0057]
所述用于基于细胞的药物筛选的分析方法的应用，其中差异表达基因是潜在的标记基因，并用于筛选针对耐药性癌症的药物，包括以下步骤：
[0058]-在培养环境中培养细胞群；
[0059]-将细胞群暴露于药物物质，其中设置了明确的时间监测点和药物物质的浓度监测点；
[0060]-暴露于药物物质后，形成单细胞文库；
[0061]-对细胞文库中的细胞的多个基因进行单细胞文库的单细胞rna测序；以及
[0062]-进行生物信息学分析，以确定差异表达基因、富集通路和发育轨迹。
[0063]
在一个实施方案中，细胞群选自原代细胞、细胞系、癌细胞。
[0064]
在一个实施方案中，细胞群是k562细胞系。培养环境是指dmem细胞培养基、10％fbs(胎牛血清)和1％ps(青霉素和链霉素)。
[0065]
在一个实施方案中，药物物质可以是任何其他抗肿瘤药物或具有抗肿瘤潜力的药物或以前未用于抗肿瘤目的的药物。
[0066]
在一个实施方案中，药物物质是奥沙利铂或赖氨酸特异性组蛋白去甲基化酶1a(lsd1)抑制剂。
[0067]
在一个实施方案中，时间监测点包括0至6小时的时间点。在优选的实施方案中，时间监测点包括0至24小时的时间点。
[0068]
在一个实施方案中，单细胞rna测序包括分析polya+mrna。
[0069]
在一个实施方案中，单细胞rna测序包括分析mrna。
[0070]
在一个实施方案中，所确定的差异表达基因包括细胞亚型的标记基因，参与药物代谢、细胞生长、增殖和分化通路。
[0071]
在一个实施方案中，所确定的差异表达基因包括bnip3、bnip3l、hemgn和kdm1a。
[0072]
在一个实施方案中，所确定的差异表达基因包括任何在细胞亚群之间差异表达的基因。在一个实施方案中，用奥沙利铂和/或任何抗癌药物处理的细胞群中，所确定的差异表达的基因是与rna转运和剪接体有关的基因。
[0073]
在一个实施方案中，用组蛋白赖氨酸特异性去甲基化酶1a(lsd1)抑制剂处理的细胞群中，所确定的差异表达基因是促凋亡基因或对组蛋白赖氨酸特异性去甲基化酶1a(lsd1)抑制剂有抗性的基因。
[0074]
在一个实施方案中，用组蛋白赖氨酸特异性去甲基化酶1a(lsd1)抑制剂处理的细胞群中，所确定的差异表达基因包括bnip3、bnip3l、hemgn和kdm1a基因。
[0075]
实施例
[0076]
实施例1
[0077]
药物处理、scrna-seq以及测序数据分析
[0078]
在培养基中培养k562细胞，然后将其暴露于50nm的奥沙利铂中。在暴露0、2、6、24小时后，清洗细胞并将其重悬在pbs中，样品表示为k00_0、k02_50、k06_50、k24_50。通过单细胞3’凝胶珠、微流控芯片以及文库试剂盒处理单细胞。在illumina测序仪上进行文库测序。
[0079]
通过umap图显示细胞系群的组成变化。对细胞进行基因集富集分析(gsea)。
[0080]
使用allprep dna/rna mini kit(qiagen，货号：80204)分离总rna，并使用qubit2.0荧光仪(life-technologies，ca，usa)和nanophotometer分光光度计(implen，ca，usa)进行测量。
[0081]
制备浓度为1
×
105个细胞/ml的pbs(hyclone)单细胞悬液。然后将单细胞悬液上样到微流控设备上，根据新格元标准操作规程，通过单细胞rna文库试剂盒(新格元生物科技)构建scrna-seq文库。将单个文库稀释到4nm并混合(pool)以待测序。在illumina hiseq x上以150bp成对末端测序读段(reads)对混合样本(pool)进行测序。
[0082]
使用内部流程对原始测序读段进行处理以形成基因表达谱。具有相同细胞条形码、umi和基因的测序读段被分为一组，以计算每个细胞的每个基因的umi数。各细胞条形码的umi计数表被用于进一步分析。使用seurat程序(http://satjalab.org/seurat/，r软件包，v.3.0.1)进行细胞类型识别和聚类分析。选择符合统计显著性原则的成分用于seurat中的umap和findneighbors分析。
[0083]
通过clusterprofiler软件对基因集进行kegg功能富集分析，以寻找与特异性表达基因显著相关的生物功能或通路。使用gseapy(v0.9.12)进行差异表达基因gsea分析。
[0084]
利用scrna-seq来揭示抗性亚群
[0085]
用seurat findclusters进行无监督聚类，细胞被聚类为7个不同的亚群(图1)。我们分析每个亚群的差异表达基因(deg)作为该组的特征，并比较各时间点亚群的变化。
[0086]
根据结果可知，亚群3的比例增加，而亚群1的比例减少，且在药物治疗期间出现了
几个新的亚群(亚群5-7)(图2)。
[0087]
亚群3的deg富集在与rna转运和剪接体有关的通路中，表明转录活跃，可能支持连续细胞分裂和快速生长(图3)。作为推论，我们用这些deg作为亚群3、5、6的特征来修正我们的数据库，这可能有助于针对含有类似耐药癌细胞亚群的患者进行药物筛选。
[0088]
实施例2
[0089]
药物处理、scrna-seq和测序数据分析
[0090]
在培养基中培养k562细胞，然后将其暴露在不同浓度(0μm，5μm，20μm)的赖氨酸特异性组蛋白去甲基化酶1a(lsd1)抑制剂中0、6、24小时。清洗细胞并使其重悬在pbs中，样品表示为lsd1_0_0、lsd1_6_0、lsd1_6_5、lsd1_6_20、lsd1_24_0、lsd1_24_5和lsd1_24_20。通过单细胞3’凝胶珠、微流控芯片以及文库试剂盒处理单细胞。在illumina测序仪上进行文库测序
[0091]
通过t-sne图显示细胞系群体的组成变化。通过cluster profiler对考虑不同数据集(基因本体数据库和kegg信号通路数据库)的细胞进行gsea。通过使用默认参数下的monocle2识别每个单一样本的伪时间轨迹。对全部细胞系的伪时间轨迹进行了分析，并对每个聚类的细胞在轨迹上的位置进行了分析。使用scrna-seq揭示耐药亚群
[0092]
评估细胞的集体异质性，我们观察到基于mrna表达的5个不同的聚类。分析每个亚群的deg作为该组的特征，并比较各个药物浓度和时间点的亚群的变化。
[0093]
根据结果可知，随着处理时间逐渐延长，亚群3首先增加，然后减少，而亚群1和2不断减少，直到最终可以忽略不计。
[0094]
在药物处理的24小时期间出现了新的亚群(聚类5)(图4)。亚群5是进入调亡通路的细胞群，在细胞刺激的反应、分泌和细胞表面相关过程中存在显著的富集(图5)(es》0，pval《0.05，其中es为富集分数，pval为p值)。纵坐标是富集分数es和基因表达与表型之间的相关性。绝对值越大，相关性越强。大于0的值表示正相关，小于0的值表示负相关。
[0095]
bnip3、bnip3l、hemgn和kdm1a存在差异性表达(图6)。hemgn是抗调亡的基因，在亚群4中有较高的表达水平，表明该聚类对该药物有一定的耐受性。从k562细胞系的细胞亚群来看，其中一些细胞在药物处理前就表现出了耐受性特征(图4)，这突出表明了即使在细胞系的群体中也存在着显著的异质性。
[0096]
聚类1和聚类2都是药物敏感细胞群。然而，可以看出聚类2中hemgn基因的表达较低，更容易被lsd1抑制剂杀死或更容易进入调亡通路。
[0097]
在lsd1抑制剂的作用下，我们观察到细胞轨迹有明显的分叉结构，聚类1、2、3可以分化为两条通路：调亡(聚类5)和耐受(聚类4)(图8)。
[0098]
bnip3和bnip3l基因是参与线粒体自噬通路的促调亡基因，表明聚类3是一类转变向凋亡通路的细胞群(图7)。
[0099]
这些deg可以作为潜在的标记基因，帮助我们在具有耐药特征的癌症中筛选治疗的药物。
附图说明
[0100]
图1.来自多个细胞系样本的基因表达谱的umap图。
[0101]
根据k562细胞系在暴露于50nm奥沙利铂中0、2、6和24小时的基因表达差异，umap
图以颜色划分的聚类显示了7个不同的聚类。
[0102]
图2.药物处理期间不同时间点的细胞系群体组成的变化。在暴露于50nm奥沙利铂0、2、6和24小时的k562细胞系中，从图1得到的7个聚类的比例。
[0103]
图3.kegg富集分析。kegg通路分析显示，deg参与rna转运和剪接体，表明转录活跃，可能支持连续细胞分裂和快速生长。
[0104]
图4.多个细胞系样本的基因表达谱的t-sne图(a)。基于暴露于0μm、5μm和20μm的组蛋白赖氨酸特异性去甲基化酶1a(lsd1)抑制剂中0、6和24小时的k562细胞系的基因表达差异，t-sne图用颜色划分的聚类显示了5个不同聚类(a)。以及暴露于0μm、5μm和20μm的lsd1中0、6和24小时的k562细胞系中，5个聚类的组成变化(b)。
[0105]
图5.通过gseapy(v0.9.12)对差异表达基因进行的gsea分析。进入调亡通路的细胞群在细胞刺激的反应、分泌和细胞表面相关过程中具有显著的正富集(es》0，pval《0.05)。
[0106]
图6.细胞系群体的潜在新标志物。特征图和小提琴图表明潜在的新标记基因(bnip3、bnip3l、hemgn和kdm1a)在每个聚类中的表达。
[0107]
图7.kegg富集分析。kegg通路分析显示，deg参与了线粒体自噬通路，表明聚类3是一类转变向凋亡通路(a)和蛋白质相互作用网络(b)的细胞群。
[0108]
图8.由单细胞转录组定义的细胞轨迹。通过monocle 2的无监督方式来构建细胞系的潜在发育轨迹。根据伪时间分析的结果可知，从分支的右端到左端的分化呈现从前到后的时间(a)。细胞系的分支轨迹显示出明显的分叉结构(b)。