一种基于单细胞多组学测序分析T淋巴细胞TCR特点的方法

摘要

本发明公开了一种基于单细胞多组学测序分析T淋巴细胞TCR特点的方法。该方法对特定病毒感染前后、治疗前后的不同患者样本进行单细胞测序，基于单细胞转录组测序数据构建单细胞基因表达矩阵并进行细胞类型的注释和分群，基于单细胞免疫组测序数据组装TCR序列并识别TCR链组成基因，然后评估和识别各细胞亚群稳定TCR的表达率，评估不同疾病状态的各细胞亚群克隆状态和各样本TCR组装与使用偏好性，以及分析同一患者治疗前后TCR的变化特点。本发明对于特定病毒感染机制的解析以及针对宿主免疫反应研发特定疫苗和药物疗法有着重要的意义。

说明书

技术领域

本发明涉及单细胞转录组测序、单细胞免疫组测序、病毒感染与宿主免疫反应和生物信息学分析领域，具体涉及一种基于单细胞多组学数据全面分析病毒感染后T淋巴细胞TCR特点的方法。

背景技术

高通量测序技术的出现使得生物学的各个领域的研究进入高通量的时代。传统基于大块组织的转录组测序技术(RNA-seq)获得的是一个组织所有细胞类型的均值，然而生物样本和组织往往是一个细胞类型多样的十分复杂的异质性系统。基于此背景，单细胞的转录组测序技术(single-cell RNA-sequencing,scRNA-seq)应运而生，首次被报道于2009年，从此生物学研究领域进入了单细胞的时代。随着scRNA-seq技术的方便与成熟，以及花费的降低，scRNA-seq技术在生物学研究的各个领域如火如荼的开展，单细胞测序研究的论文数量与日俱增，研究中细胞数量也达到了井喷式增长，从开始的几百几千个细胞，发展到现在的几十万甚至上百万个细胞数量级。

随着scRNA-seq技术的迅猛发展，单细胞研究也迎来了多组学的时代，众多单细胞多组学测序技术相应产生，被有效地利用到生命科学领域的各种问题中，例如：单细胞转录组测序技术，单细胞免疫组测序技术，单细胞基因组测序技术，单细胞染色质可及性测序技术。特别的，随着原位杂交和特定RNA探针的出现，单细胞空间组学技术今年发展迅速，而且被多个期刊评为年度技术，因而被寄予厚望。目前相对成熟且应用最多的是单细胞转录组测序技术，代表性的为10x Genomics公司的基于3’或5’端的转录本测序的Chromium技术。因为该技术只测序转录本的一端，因此价格相对便宜，细胞通量也很高。

病毒感染性疾病是目前威胁人类健康的主要杀手之一。例如人类免疫缺陷病毒(Human immunodeficiency virus，HIV)感染可以直接摧毁人类免疫系统，引起艾滋病(Acquired immune deficiency syndrome，AIDS)；狂犬病毒(Rabies virus，RABV)是一种神经感染病毒，发病率接近100％，发病后至今没有有效的药物和治疗手段；特别地，2020年年初，新冠肺炎(The coronavirus disease 2019，COVID-19)疫情在世界范围内快速暴发。COVID-19由严重急性呼吸综合征冠状病毒2(Severe acute respiratory syndromecoronavirus 2，SARS-CoV-2)感染引起，是一种急性感染性肺炎，患者以肺部病变为主。

病毒侵染人体后会激发机体的抗病毒免疫反应，涉及天然免疫反应和适应性免疫反应两种。天然免疫反应相比适应性免疫反应响应时间更早，响应范围更大。适应性免疫反应又分为细胞免疫和体液免疫两大类，其中T淋巴细胞主要介导了细胞免疫反应。T淋巴细胞来源于骨髓中的淋巴干细胞，在胸腺中分化发育成熟，通过血液循环和淋巴循环分布到全身各处，从而介导细胞免疫反应发挥作用。初始状态T细胞经过抗原递呈细胞刺激激活后，分化为效应T细胞和记忆T细胞，前者发挥杀伤靶细胞的功能，后者当相同的抗原再次出现时，可以快速分化为效应T细胞发挥作用。T细胞的激活依赖T细胞受体(T cellreceptor,TCR)和共刺激信号的双信号刺激。激活的效应T细胞可以与靶细胞结合，同时释放干扰素、穿孔素和颗粒酶等物质杀死被病毒感染的细胞或肿瘤细胞。T淋巴细胞主要由α和β链组成TCR的T细胞为主，在分化成熟的过程中，α和β链基因在特异性重组酶的作用下将DNA上分隔的、原本没有转录活性的基因片段连接成一个完整的、有转录功能的活性基因。大致过程为：β链依次进行D、J连接，再进行V、DJ连接，后再与C区基因相连，形成一个完整的β链功能基因。β基因的重排会诱导α基因的重排。α链无D片段，直接进行V、J连接，再与C区相连，形成一个完整的α链功能基因。正常的α和β链TCR基因重组是相对随机的过程，病毒感染后，T淋巴细胞发生克隆扩增。TCR克隆扩增的比例、状态、多样性和偏好性等对病毒的清除和理解病毒感染后机体清除病毒机制进程有着重要作用。

当今，scRNA-seq已经成功应用到许多前沿的生物医学领域，其中就包括病毒感染和宿主免疫反应的研究。scRNA-seq成功研究了多种病毒感染的过程机制，并揭示了宿主天然免疫和适应性免疫的反应特征，为疫苗的研究和临床治疗提供了宝贵的思路。然而，目前单一利用scRNA-seq转录组测序技术或免疫组测序技术对病毒感染过程的研究主要基于对宿主特定细胞和免疫系统组成和转录组状态的定量分析，或特异性序列克隆程度的分析。同时结合单细胞转录组和免疫组技术的分析也并未全面系统地实现不同疾病状态(感染前后、康复治疗前后)，不同患者，以及同一患者不同疾病状态的特定亚群的TCR克隆特点的细致深入分析。因此无法全面深入地分析病毒感染的疾病进程中，不同疾病状态下的各个细胞亚群的稳定TCR表达率和克隆状态，各样本在病毒感染前后以及治疗前后TCR组装与使用偏好性，以及同一患者治疗前后的TCR变化的状态、特异性和多样性等。

发明内容

本发明的目的是针对上述现有分析技术的不足，提供一种基于单细胞多组学测序数据全面分析病毒感染后T淋巴细胞TCR特点的方法。该方法能够实现评估和识别各细胞亚群稳定TCR的表达率，评估不同疾病状态的各细胞亚群克隆状态和各样本TCR组装与使用偏好性，以及分析同一患者治疗前后TCR的变化特点，对病毒感染机制的研究以及针对宿主免疫反应研发特定疫苗和药物疗法有着重要的意义。

本发明的技术方案如下：

一种基于单细胞多组学测序数据全面分析病毒感染后T淋巴细胞TCR特点的方法，采集特定病毒感染前后、治疗前后的不同患者的样本进行单细胞测序，然后通过以下步骤进行分析：

1)基于单细胞转录组测序数据构建单细胞基因表达矩阵并进行细胞类型的注释和分群：在相应的参考基因组序列网站下载宿主参考基因组序列，同时下载对应版本的参考基因组的基因/功能亚基注释文件(GTF文件)。通过数据比对、基因表达定量获得单细胞基因表达矩阵。通过单细胞转录组测序分析软件进行数据降维以及细胞类型的分群与注释。

2)基于单细胞免疫组测序数据组装TCR序列并识别TCR链组成基因：在已经下载好的宿主参考基因组序列和基因/功能亚基注释文件中，特异性提取T淋巴细胞TCR相关的基因序列和注释信息，构建TCR参考基因组信息。通过数据比对和基因表达定量组装TCR序列并识别TCR链组成基因。

3)结合表达数据识别各细胞亚群的稳定TCR表达率：设定合理阈值对步骤1)获得的单细胞表达数据进行数据过滤和去噪，只保留高质量的单细胞表达数据。TCR序列定量和识别后，只保留同时具有α和β链数据的TCR。最终只保留同时具有高质量单细胞表达数据和TCR数据的细胞，从而获得各细胞亚群的稳定TCR表达率，即各细胞亚群中检测到稳定TCR的细胞数占总细胞数的百分比。

4)评估不同疾病状态的各细胞亚群克隆状态：将数据按不同的疾病状态分类，再将同一种疾病状态下的数据根据不同的细胞类型分组，定义克隆和非克隆状态，评估不同疾病状态和细胞类型的整体克隆特点。

5)分析不同疾病状态下各样本TCR组装与使用偏好性：通过分析不同患者、不同疾病状态下各样本的TCRα和β链V-J基因配对使用的特点，解析不同疾病状态下各样本TCR组装与使用偏好性。

6)分析同一患者治疗前后整体和特定细胞亚群的TCR变化的状态、特异性和多样性：针对所有细胞类型或者表达量和细胞比例具有明显差异的感兴趣的特定细胞类型，对同一患者治疗前后整体TCR克隆变化，所有细胞亚群的TCR克隆变化，以及特定细胞亚群的TCR克隆变化进行深入分析。

上述步骤1)中对宿主参考基因组序列和注释以及细胞类型注释的要求如下：

宿主参考基因组序列和对应的基因注释信息可以来自但不限于UCSC、NCBI、Ensemble、Genecode以及10x Genomics等数据库整理和维护的人类参考基因组hg18(GRCh36)、hg19(GRCh37)、hg20(GRCh38)或者染色体组的子集。宿主参考基因组序列文件应为标准的FASTA格式的文件，基因注释文件应该是标准的GTF格式的文件。细胞类型的分群可以使用但不限于Seurat、Scanpy、Cellranger Loupe等软件。细胞类型的注释基于文献报道的经典细胞表面标志基因或特异性表达的标志性基因。

上述步骤2)组装TCR序列并识别TCR链组成基因，其中：

宿主TCR序列和对应的基因注释提取自宿主参考基因组序列和注释信息。提取标签主要是但不限于：--attribute＝gene_biotype:TR_V_gene/TR_V_pseudogene/TR_D_gene/TR_J_gene/TR_J_pseudogene/TR_C_gene。TCR参考基因组序列文件应为标准的FASTA格式的文件，基因注释文件应该是标准的GTF格式的文件。组装TCR序列并识别TCR链组成基因可以使用但不限于Cellranger vdj软件。

上述步骤3)结合表达数据过滤不稳定的TCR从而识别稳定的表达率，其中：

对单细胞表达数据进行数据过滤和去噪，只保留高质量的单细胞表达数据，过滤细胞的指标包括但不限于：基因数检测过低的细胞，UMI数检测过低的细胞，线粒体比例过高的细胞，具有两个或两个以上经典细胞标志基因的双细胞，没有任何明显标志性基因表达的细胞碎片。TCR序列组装和定量识别后，仅保留具有至少一条可表达的TCRα链(TRA)和一条可表达的TCRβ链(TRB)的细胞用于进一步分析。最终将两个集合的细胞取交集，只保留同时具有高质量表达数据和TCR数据的细胞进行下游分析。

上述步骤4)评估不同疾病状态的各细胞亚群克隆状态，其中：

在获得不同细胞TCR序列、基因组成和表达量后，每个独特的TRA(s)-TRB(s)对被定义为一个克隆型。如果一种克隆型存在于至少两个细胞中，则携带该克隆型的细胞被认为是克隆的。在此，可以同时评估某个疾病状态下某个细胞类型中没检测到TCR、TCR不克隆以及TCR克隆的情况。

上述步骤5)分析不同疾病状态下各样本TCR组装与使用偏好性，其中：

正常的α和β链TCR基因重组是相对随机的过程，抗原刺激后发生克隆扩增，β链依次进行D、J连接，再进行V、DJ连接，后再与C区基因相连，形成一个完整的β链功能基因。α链无D片段，直接进行V、J连接，再与C区相连，形成一个完整的α链功能基因。因此在同时分析TRA与TRB时，特异性分析不同样本中V基因与J基因配对使用的偏好性。

上述步骤6)分析同一患者治疗前后整体和特定亚群的TCR变化的状态、特异性和多样性，其中：

对于每一个克隆型，具有该克隆型的所有细胞的个数称为该克隆型的克隆大小。如果一个克隆型在治疗前后都存在，则该克隆型为治疗前后共享克隆，称为维持性克隆；如果一个克隆型在治疗前或后特异性存在，则该克隆型为特异性克隆。通过序列分析鉴定不同的克隆型，通过分析克隆大小和特异性来解析同一患者治疗前后整体和特定细胞亚群的TCR克隆变化特点。

本发明针对目前利用scRNA-seq技术对病毒感染过程中TCR在整体细胞类型、特定细胞亚群、疾病状态、单一患者和同一患者治疗前后的研究的不足，提供一种基于单细胞多组学数据全面分析病毒感染后T淋巴细胞TCR特点的方法。该方法能够结合表达数据识别各细胞亚群的稳定TCR表达率、评估不同疾病状态的各细胞亚群克隆状态、分析不同疾病状态下各样本TCR组装与使用偏好性，以及分析同一患者治疗前后的TCR变化的状态、特异性和多样性，对特定病毒感染宿主的机制解析以及特异性疫苗和药物的研发有着重要的意义。

附图说明

图1.本发明实施例中分析病毒感染后T淋巴细胞TCR特点的整体分析流程。

图2.本发明实施例中结合表达数据识别各细胞亚群的稳定TCR表达率整体情况，其中：左图为细胞分类图，数字代表不同的细胞类型；中间图显示了细胞是否检测到稳定TCR的情况；右图为不同细胞类型的细胞亚群中检测到稳定TCR信号的百分比。

图3.本发明实施例中不同疾病状态的各细胞亚群克隆状态评估。

图4.本发明实施例中分析不同疾病状态下各样本TCR组装与使用偏好性。

图5.本发明实施例中对同一患者治疗前后整体TCR克隆变化的分析。

图6.本发明实施例中对同一患者治疗前后所有细胞亚群的TCR克隆变化的分析。

图7.本发明实施例中对同一患者治疗前后特定亚群的TCR克隆变化的分析。

具体实施方式

以下所述是本发明的更加细致的实施描述，其中的参数以及具体实施细节用以解释本发明的可行性及实施效果，并不构成对本发明的限定。

本实施例以14例HIV感染的患者(包含9例治疗前患者和8例抗病毒治疗后患者，其中3例患者具有配对的抗病毒治疗前后的样本)和4例健康人的外周血单个核细胞(PBMC)中磁珠分离的CD4

1.基于Single Cell 5'Library&Gel Bead Kit(10x Genomics)和ChromiumSingle Cell A Chip Kit(10x Genomics)的样本建库与高通量单细胞测序

使用CD4(130-045-101,Miltenyi Biotech)分离试剂盒立即从新鲜PBMC中纯化CD4

2.构建单细胞转录组测序数据基因表达矩阵及细胞类型注释

使用10x Genomics公司配套开发的Cell Ranger(v.3.0.2)软件进行测序数据的比对、计数以及单细胞基因表达矩阵的构建等，主要分为两步：(1)cellranger mkfastq；(2)cellranger count。如此，对每个样本都可以获得一个行代表不同基因、列代表不同细胞的单细胞基因表达矩阵。单细胞转录组测序数据的过滤、标准化、降维以及聚类等步骤均使用R包Seurat(版本3.5.3)完成，步骤如下：(1)“CreateSeuratObject”函数用来构建Seurat对象，只保留至少在0.1％细胞中表达的基因和至少有200个基因表达的细胞来构建对象；(2)“subset”函数用来进行细胞的二次过滤。过滤掉基因数小于500、UMI数小于1000和线粒体含量大于10％的细胞；(3)“NormalizeData”函数完成数据的标准化；(4)“FindVariableFeatures”构建细胞间高度差异的基因集合；(5)“ScaleData”函数进行数据的均一化，去除无意变量的影响；(6)“RunPCA”进行主成分分析；(7)“FindNeighbors”和“FindClusters”函数进行细胞间距离的计算和细胞的聚类；(8)利用“RunTSNE”函数完成高维数据的非线性空间降维，将细胞映射到二维空间进行结果的可视化；(9)最后使用“FindAllMarkers”进行不同类群标志性基因的鉴定。细胞类型的鉴定主要依据之前的先验知识。以上步骤详细的操作参考于Seurat官方教程(https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html)。

3.TCR序列组装、克隆型过滤与鉴定

使用Cell Ranger(v.3.0.2)V(D)J流程程序和GRCh38作为参考进行TCR组装与克隆型鉴定。在获得了包含克隆型组成和频率列表后，仅保留具有至少一条可表达的TCRα链(TRA)和一条可表达的TCRβ链(TRB)的细胞用于进一步分析。每个独特的TRA(s)-TRB(s)对被定义为一个克隆型。如果一种克隆型存在于至少两个细胞中，则携带该克隆型的细胞被认为是克隆的。某个克隆型的细胞数量为该克隆型的克隆大小，表明克隆型的克隆程度。

4.分析不同细胞亚群TCR的捕获效率和不同疾病状态下的克隆状态

使用细胞条形码信息(Barcode)，将具有TCR克隆型的数据投影到表达数据的t-SNE图上。由此进行再过滤，如果一群细胞表达两个或两个以上不同细胞类型的经典标志基因，则暗示该群细胞大概率是建库过程中产生的双细胞。如果一个类群标志性基因多为线粒体相关基因、核糖体基因或者没有显著性基因，往往这些细胞的TCR捕获效率也非常低，则提示该亚群高度可能是低质量细胞。所有的双细胞亚群和低质量细胞亚群都从总的细胞矩阵中剔除，不进行下游分析。在获得了高质量表达数据和TCR数据后，计算不同的细胞类型中TCR检测的稳定比例，为该细胞类型中检测到稳定TCR的细胞数占该细胞类型总细胞数的百分比，结果如图2所示。进一步，将不同的细胞类型按照不同的疾病状态分开，计算某一个特定细胞类型在某一个特定疾病状态下的克隆状态，包括未检测到TCR的比例、检测到TCR为非克隆状态的比例，以及检测到TCR状态为克隆状态的比例，结果如图3所示。

5.不同疾病状态下各样本TCR组装与使用偏好性的分析

对于每个样本，我们获得该样本每个细胞TCR的克隆型组成，分别提取该细胞TRA和TRB中V-J基因的配对信息，进而统计该样本所有TCR中TRA和TRB的V-J基因配对使用的频率。统计好所有样本的TCR V-J配对使用情况后，将样本按照不同的疾病状态和不同的患者分组，选取每个样本中高频使用V-J基因对(Top5)绘制成图，有效地展示不同患者之间特异性的和共享的V基因、J基因和V-J基因对，以及不同疾病状态之间特异性的和共享的V基因、J基因和V-J基因对。由此解析病毒感染后相比健康状态，以及抗病毒治疗后差异性的高频V-J基因对使用的偏好性。分析结果如图4所示。

6.分析同一患者治疗前后整体和特定细胞亚群的TCR变化

对于同一患者具有抗病毒治疗前后配对的数据，我们全面分析治疗前后克隆的变化情况。首先分析克隆种类的变化，将克隆分为治疗前特异性克隆、治疗后特异性克隆，以及治疗前后共享的克隆。通过散点图比较HIV感染后，3个配对患者克隆变化的规律是否一致，主要发生的是克隆的新生扩张，还是克隆的维持，结果如图5所示。接着我们统计特异性克隆和维持性克隆在不同疾病状态下的不同细胞类型中所占的相对比例，以揭示特定疾病状态克隆变化的多样性以及特定细胞亚群在抗病毒治疗前后克隆发生变化的特异性，例如我们发现表达GNLY基因的效应CD4

以上所述仅作为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之类的所做的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。