CN108491689B - 基于转录组的肿瘤新抗原鉴定方法 - Google Patents

Abstract

本发明公开了基于转录组的肿瘤抗原鉴定方法，包括获取患者肿瘤组织的RNA样本，对RNA样本建库和扩增，获得肿瘤组织的RNA样本测序结果；将RNA样本测序结果的短读段比对到人类参考基因组，获得RNA比对结果；根据RNA比对结果计算基因表达量，根据RNA比对结果进行突变检测和预测融合基因事件；根据比对结果预测转录组HLA分型；计算基因表达量、突变检测和预测融合基因事件；将转录组样本的基因表达量，转录组突变位点在全外测序样本中的深度，新生短肽与患者HLA分型的结合力作为分析结果交给下游分析人员四个步骤。本发明提供了一种能够从肿瘤患者转录组NGS数据中鉴定出个体样本的肿瘤特异性抗原的方法。

Description

基于转录组的肿瘤新抗原鉴定方法

技术领域

本发明涉及一种肿瘤新生抗原的生物信息鉴定方法。

背景技术

目前常规的基因变异检测技术包括Sanger测序，焦磷酸测序，扩增阻滞辨析系统(ARMS)、荧光原为杂交技术及等位基因特异PCR等。这些基因检测手段普遍通量较低，且费用高，耗时长。

随着新一代测序技术(Next-Generation Sequencing,NGS)的出现，大规模平行测序成为可能。相较于第一代测序技术，NGS技术检测速度快，准确率高，成本低，覆盖度广。利用NGS技术可以进行肿瘤基因组序列的重测序，继而对肿瘤相关基因的不同类型变异(点突变、插入、缺失等)、拷贝数变异以及基因相关性和多态性进行鉴定。目前肿瘤NGS检测可应用于遗传性肿瘤筛查和突变鉴定。

在肿瘤免疫治疗中，鉴定肿瘤组织细胞产生的新抗原(neoantigen)是决定下游临床治疗的关键步骤。正常细胞在癌变过程中，由于遗传物质发生改变，导致其DNA序列与其他正常细胞有差异，因此会产生肿瘤相关抗原(肿瘤细胞高表达)或肿瘤特异性抗原(只在肿瘤细胞中表达)。这些抗原由于具有标识肿瘤细胞的特异性抗原决定簇，理论上可以被抗原呈递细胞上的人类白细胞抗原(HLA)识别，然后与TCR结合，进而激活T细胞，启动免疫反应，是肿瘤免疫治疗潜在的靶点。通过对病人肿瘤组织和正常组织NGS数据的分析，可以鉴定出该组织的突变，例如点突变、插入缺失突变、结构变异、基因融合等。从这些突变事件可以预测肿瘤细胞在下游转录与表达过程发生的变化，进而推测出可能存在的新抗原，转录组测序数据相对于基因组而言，可以更为直接地反映遗传物质的变化，不但可以鉴定基因组水平突变在中心法则上的呈现情况，还可以印证基因组测序数据分析结果的可靠性，为肿瘤新生抗原鉴定提供更多依据和证据，进而指导临床治疗。

发明内容

本发明的目的在于提供一种能够从肿瘤患者转录组NGS数据中鉴定出个体样本的肿瘤特异性抗原的方法。

基于转录组的肿瘤抗原鉴定方法，包括以下步骤：

S1：获取患者肿瘤组织的RNA样本，对RNA样本建库和扩增，获得肿瘤组织的RNA样本测序结果；

S2：将RNA样本测序结果的短读段比对到人类参考基因组，获得RNA比对结果；

S3、根据RNA比对结果计算基因表达量，根据RNA比对结果进行突变检测，根据RNA比对结果预测融合基因事件；根据比对结果预测转录组HLA分型；计算基因表达量、突变检测和预测融合基因事件按指定顺序进行，或者同步进行；

S3A：计算基因表达量包含：去除RNA比对结果中的内含子短读段，形成仅由外显子短读段组成的的RNA比对结果，计算基因表达量；

S3B：突变检测包含：

S3B1：从RNA比对结果中检测RNA样本中的突变，去掉正常组织的全外显子样本中具有的突变，剩余的突变作为RNA水平的体细胞突变；对RNA的突变进行功能注释；

S3B2：计算突变位点的深度：针对每个突变位点，计算当前突变位点在全外显子测序样本中的深度，全外显子样本包括肿瘤组织的全外显子样本和正常组织的全外显子样本；

S3B3：HLA分型亲和力预测：根据突变的注释结果，针对每一个突变事件编辑相应转录本序列，获得一条包含所有突变的多肽序列，去掉在野生型蛋白序列中可以匹配到的肽段形成转录组特有的新生多肽序列，将转录组特有的多肽序列按照HLA分型结合需要的长度截取成短肽，将短肽与步骤S1的RNA样本测序结果中鉴定的患者HLA分型做亲和力预测；

其中，计算突变位点的深度和HLA分型亲和力预测按指定顺序进行或同步进行；

S3C、预测融合基因事件：根据步骤S2的RNA比对结果预测转录组的融合基因事件，从中筛选出可信的融合基因后，翻译成蛋白序列，将蛋白序列按照HLA分型结合需要的长度截取成短肽；将短肽与步骤S1的RNA样本测序结果中鉴定的患者HLA分型做亲和力预测；

S4：将转录组样本的基因表达量，转录组突变位点在全外测序样本中的深度，新生短肽与患者HLA分型的结合力作为分析结果交给下游分析人员。

进一步，步骤S3A中，计算基因表达量包含以下步骤：

S3A1、去除建库过程中由PCR扩增产生的重复的短读段序列；

S3A2、将比对到内含子的短读段去除,获得由外显子短读段组成的RNA比对结果；

S3A3、计算每个外显子短读段数，用以下公式换算成基因的表达量：

其中total exon reads表示比对到某个外显子区域的总短读段数，mapped reads(millions)表示每1000000条短读段中比对到该区域的短读段数，exon length(KB)表示这段外显子的长度。

进一步，步骤S2获得基因表达量后，对RNA比对结果进行质控，质控的内容包含：整个捕获区域的覆盖度，平均测序深度，比对率，重复序列的比重和唯一比对reads的比重；若比对结果符合质控要求，则进入S3；若比对结果不符合质控要求，则重新获取样本。

进一步，步骤S3B2包含所有突变的多肽序列的获取方法为：根据转录本编号从ensembl数据库中找出相应的CDS区域核苷酸序列，拼接下游3’端序列(3’UTR)，按照突变的功能注释将核苷酸序列上相应位置的核苷酸做出修改，得到一条包含所有突变的核苷酸序列；如功能注释给出的是G100A，则将拼接核苷酸序列的第100位的碱基G修改为A；将包含所有突变的转录本核苷酸序列翻译成多肽序列。

进一步，步骤S3B和S3C中截取短肽的方法为：

SI、顺序获取多肽序列上的突变为当前突变，以当前突变为中心、向前截取n个氨基酸和向后截取m个氨基酸获得多肽序列，n为HLA分型所能呈递的最大长度，m为HLA分型所能呈递的最大长度、或者m为从当前突变到第一个终止密码子的长度；

SII，在多肽序列上依次截取长度为N的短肽，N为HLA分型结合需要的长度；获得N条包含当前突变的短肽。例如：HLA分型结合需要的长度为8，N＝8；从突变位置开始向前截取7个氨基酸、与该突变位置形成一条长度为8个氨基酸的短肽；从突变位置开始向前截取6个氨基酸，向后截取1个氨基酸，这7个氨基酸与突变位置形成第2条长度为8个氨基酸的短肽，如此，一共能够获得8条含有该突变的短肽。

进一步，步骤SI截取含有突变位点的多肽序列的规则为：

A、对于点突变：以发生突变的位置为中心，向前截取n个氨基酸，向后截取m个氨基酸，n＝m＝HLA分型所能呈递的最长肽段，若前段或后段长度不足时，则有多少截多少；如果点突变属于stop loss(终止密码子丢失)，则m为从当前突变到第一个终止密码子的长度；

B、对于非移码突变：非移码插入要从插入序列的第一个氨基酸向前截取n个氨基酸；从插入序列的最后一个氨基酸向后截取m氨基酸，n＝m＝HLA分型所能呈递的最长肽段；以插入序列为中心，插入序列、插入序列之前的n个氨基酸和插入序列之后的m个氨基酸和共同组成新生多肽序列；

非移码缺失则以缺失位点为中心，分别向前截取n个氨基酸，向后截取m个氨基酸，n＝m＝HLA分型所能呈递的最长肽段；

C、对于移码突变：从开始移码突变的第一个氨基酸为中心，向前截取n个氨基酸，n＝HLA分型所能呈递的最长肽段；向后截取m个氨基酸，m为从当前突变到第一个终止密码子的长度，即向后截取至第一个终止密码子。

本发明的优点在于：

1.整个鉴定过程从fastq文件开始，使用者无需准备其他输入文件。

2.所有鉴定步骤均有相应质控步骤，提高了结果的准确性。

3.对于突变位点的鉴定更加全面，不仅考虑单个细胞的情况，也考虑整块组织的突变分布。

4.有多组学鉴定结果的印证。

5.关键步骤都有多种算法同时进行，既可相互验证，也降低了结果的假阴性。

6.关键步骤都优化了算法，并引入并行计算，缩短了单个样本的分析时间。

具体实施方式

基于转录组的肿瘤抗原鉴定方法，包括以下步骤：

S1：获取患者肿瘤组织的RNA样本，对RNA样本建库和扩增，获得肿瘤组织的RNA样本测序结果。

S2：将RNA样本测序结果的短读段比对到人类参考基因组，获得RNA比对结果。

S3、根据RNA比对结果并计算基因表达量，根据RNA比对结果进行突变检测，和根据RNA比对结果预测融合基因事件；根据比对结果预测HLA分型；计算基因表达量、突变检测和预测融合基因事件按指定顺序进行，或者同步进行。

S4：将转录组样本的基因表达量；转录组突变位点在全外测序样本中的深度；新生短肽与患者HLA分型的结合力作为分析结果交给下游分析人员。

计算基因表达量

去除RNA比对结果中的内含子短读段，形成仅由外显子短读段的RNA比对结果，计算基因表达量。

计算基因表达量包含以下步骤：

S3A1、去除建库过程中由PCR扩增产生的重复的短读段序列；在测序建库阶段，随机打断的核苷酸序列在flowcell里会扩增成簇，一簇的序列都是一样的，称为重复序列，主要是为了增加碱基被测到的概率和可信度，在本步骤中，将重复序列去除，捕获测序的时候，如果目标区域的重复序列不去除，就无法获知该区域真实的覆盖度和深度，对突变鉴定和表达量计算都会产生干扰。

S3A2、使用GATK工具包中的split_N_cigar功能将将转录到RNA上的内含子短读段片段切除获得仅由外显子短读段组成的RNA比对结果；虽然我们测序得到的是成熟RNA，但有些短读段会比对到内含子区域，因为短读段只有150bp的长度，有很大可能会比对到实际上不属于它的地方，因此需要在转录组鉴定中内含子区域的短读段去掉。

S3A3、使用HTseq计算每个基因区域比对到的短读段数，用以下公式换算成基因的表达量：

其中total exon reads表示比对到某个外显子区域的总短读段数，mapped reads(millions)表示每1000000条短读段中比对到该区域的短读段数，exon length(KB)表示这段外显子的长度。

计算突变位点深度

计算突变点在全外显子样本中的深度包含以下步骤：

S3B.1：从RNA比对结果中检测RNA样本中的突变，去掉正常组织的全外显子样本中具有的突变，剩余的突变作为RNA水平的体细胞突变；对RNA的突变进行功能注释；

S3B.2：计算突变位点的深度：针对每个突变位点，计算当前突变位点在全外显子测序样本中的深度，全外显子样本包括肿瘤组织的全外显子样本和正常组织的全外显子样本，可使用Bam-recount等工具计算转录组突变位点在全外显子样本中的深度，Bam-recount只能计算SNV的深度，javakit可以计算indel的深度，两者结果合并作为最终的结果。全外显子样本包括肿瘤组织的全外显子样本和正常组织的全外显子样本。

HLA分型亲和力预测S3C.1：从RNA比对结果中检测RNA样本中的突变，去掉正常组织的全外显子样本中已鉴定的突变，剩余的突变作为RNA水平的突变；对RNA的突变进行功能注释；

S3C.2：根据突变的注释结果，针对每一个突变事件编辑相应转录本序列，获得一条包含所有突变的核苷酸序列，然后翻译成多肽序列，去掉多肽序列中在野生型蛋白可以匹配到的肽段形成转录组特有的新生多肽序列。将转录组特有的多肽序列按照HLA分型结合需要的长度截取成短肽，将短肽与步骤S1的RNA样本测序结果中鉴定的患者HLA分型做亲和力预测。

转录本核苷酸序列的获得方法为：根据转录本编号从ensembl数据库中找出相应的CDS区域核苷酸序列，拼接下游3’端序列，按照突变的功能注释将核苷酸序列上相应位置的核苷酸做出修改，得到一条包含所有突变的核苷酸序列。如功能注释给出的是G100A，则将拼接核苷酸序列的第100位的碱基G修改为A；将包含所有突变的转录本核苷酸序列翻译成多肽序列。

使用多个HLA分型鉴定工具鉴定转录组样本的HLA分型，如SOAP-HLA等，将各个算法获得的结果综合考量作为最终结果。

分别使用多个HLA分型亲和力预测软件进行亲和力预测，如针对HLA I型(netMHC4.0等)，针对HLA II型(netMHCII 2.2等)；预测结果中，保留至少有一个软件判定为有亲和力的结果。

预测融合基因事件

根据步骤S2的RNA比对结果预测转录组的融合基因事件，通过SOAPfuse(v1.27)、STAR-Fusion等工具进行预测。从中筛选出可信的融合基因后，翻译成蛋白序列，将蛋白序列按照HLA分型结合需要的长度截取成短肽；将短肽与步骤S1的RNA样本测序结果中鉴定的患者HLA分型做亲和力预测；

步骤S2获得基因表达量后，对RNA比对结果进行质控，质控的内容包含：整个捕获区域的覆盖度，平均测序深度，比对率，重复序列的比重和唯一比对reads的比重；若比对结果符合质控要求，则进入S4；若比对结果不符合质控要求，则重新获取样本。

截取短肽

截取短肽的方法为：

SI、顺序获取多肽序列上的突变为当前突变，以当前突变为中心、向前截取n个氨基酸和向后截取m个氨基酸获得多肽序列，n为HLA分型所能呈递的最大长度，m为HLA分型所能呈递的最大长度、或者m为从当前突变到第一个终止密码子的长度；

SII，在多肽序列依次截取长度为N的短肽，N为HLA分型结合需要的长度；获得N条包含当前突变的短肽。例如：HLA分型结合需要的长度为8，N＝8；从突变位置开始向前截取7个氨基酸、与该突变位置形成一条长度为8个氨基酸的短肽；从突变位置开始向前截取6个氨基酸，向后截取1个氨基酸，这7个氨基酸与突变位置形成第2条长度为8个氨基酸的短肽，如此，一共能够获得8条含有该突变的短肽。

针对不同类型的突变、截取短肽的规则为：

A、对于点突变：以发生突变的位置为中心，向前截取n个氨基酸，向后截取m个氨基酸，n＝m＝HLA分型所能呈递的最长肽段，若前段或后段长度不足时，则有多少截多少；如果点突变属于stop loss(终止密码子丢失)，则m为从当前突变到第一个终止密码子的长度；

B、对于非移码突变：非移码插入要从插入序列的第一个氨基酸向前截取n个氨基酸；从插入序列的最后一个氨基酸向后截取m氨基酸，n＝m＝HLA分型所能呈递的最长肽段；以插入序列为中心，插入序列、插入序列之前的n个氨基酸和插入序列之后的m个氨基酸和共同组成多肽序列；

非移码缺失则以缺失位点为中心，分别向前截取n个氨基酸，向后截取m个氨基酸，n＝m＝HLA分型所能呈递的最长肽段；

C、对于移码突变：从开始移码突变的第一个氨基酸为中心，向前截取n个氨基酸，n＝HLA分型所能呈递的最长肽段；向后截取m个氨基酸，m为从当前突变到第一个终止密码子的长度，即向后截取至第一个终止密码子。

具体实施例是对本发明进行进一步说明，但本发明的并不局限于这些具体实施方式。

Claims (4)

1.基于转录组的肿瘤抗原鉴定方法，其特征在于，包括以下步骤：

S1：获取患者肿瘤组织的RNA样本，对RNA样本建库和扩增，获得肿瘤组织的RNA样本测序结果，根据测序结果鉴定患者的HLA分型；

S2：将RNA样本测序结果的短读段比对到人类参考基因组，获得RNA比对结果；

S3、根据RNA比对结果计算基因表达量，根据RNA比对结果进行突变检测，根据RNA比对结果预测融合基因事件；根据比对结果预测转录组HLA分型；计算基因表达量、突变检测和预测融合基因事件按指定顺序进行，或者同步进行；

S3A：计算基因表达量包含：去除RNA比对结果中的内含子短读段，形成仅由外显子短读段组成的RNA比对结果，计算基因表达量，

计算基因表达量包含以下步骤：

S3A1、去除建库过程中由PCR扩增产生的重复的短读段序列；

S3A2、将比对到内含子的短读段去除,获得由外显子短读段组成的RNA比对结果；

S3A3、计算每个外显子短读段数，用以下公式换算成基因的表达量：

其中total exon reads表示比对到某个外显子区域的总短读段数，mapped reads(millions)表示每1000000条短读段中比对到该区域的短读段数，exon length(KB)表示这段外显子的长度；

S3B：突变检测包含：

S3B1：从RNA比对结果中检测RNA样本中的突变，去掉正常组织的全外显子样本中具有的突变，剩余的突变作为RNA水平的体细胞突变；对RNA的突变进行功能注释；

S3B2：计算突变位点的深度：针对每个突变位点，计算当前突变位点在全外显子测序样本中的深度，全外显子样本包括肿瘤组织的全外显子样本和正常组织的全外显子样本；

S3B3：HLA分型亲和力预测：根据突变的注释结果，针对每一个突变事件编辑相应转录本序列，获得一条包含所有突变的多肽序列，去掉在野生型蛋白序列中可以匹配到的肽段形成转录组特有的新生多肽序列，将转录组特有的多肽序列按照HLA分型结合需要的长度截取成短肽，将短肽与步骤S1的RNA样本测序结果中鉴定的患者HLA分型做亲和力预测；

其中，计算突变位点的深度和HLA分型亲和力预测按指定顺序进行或同步进行；

S3C、预测融合基因事件：根据步骤S2的RNA比对结果预测转录组的融合基因事件，从中筛选出可信的融合基因后，翻译成蛋白序列，将蛋白序列按照HLA分型结合需要的长度截取成短肽；将短肽与步骤S1的RNA样本测序结果中鉴定的患者HLA分型做亲和力预测；

步骤S3B3和S3C中截取短肽的方法为：

SI、顺序获取多肽序列上的突变为当前突变，以当前突变为中心、向前截取n个氨基酸和向后截取m个氨基酸获得多肽序列，n为HLA分型所能呈递的最大长度，m为HLA分型所能呈递的最大长度、或者m为从当前突变到第一个终止密码子的长度；

SII，在多肽序列上依次截取长度为N的短肽，N为HLA分型结合需要的长度；获得N条包含当前突变的短肽；

S4：将转录组样本的基因表达量，转录组突变位点在全外测序样本中的深度，新生短肽与患者HLA分型的亲和力作为分析结果交给下游分析人员。

2.根据权利要求1所述的基于转录组的肿瘤抗原鉴定方法，其特征在于：步骤S2获得基因表达量后，对RNA比对结果进行质控，质控的内容包含：整个捕获区域的覆盖度，平均测序深度，比对率，重复序列的比重和唯一比对reads的比重；若比对结果符合质控要求，则进入S3；若比对结果不符合质控要求，则重新获取样本。

3.根据权利要求2所述的基于转录组的肿瘤抗原鉴定方法，其特征在于：步骤S3B3包含所有突变的多肽序列的获取方法为：根据转录本编号从ensembl数据库中找出相应的CDS区域核苷酸序列，拼接下游3’端序列，按照突变的功能注释将核苷酸序列上相应位置的核苷酸做出修改，得到一条包含所有突变的核苷酸序列；将包含所有突变的转录本核苷酸序列翻译成多肽序列。

4.根据权利要求3所述的基于转录组的肿瘤抗原鉴定方法，其特征在于：步骤SI截取含有突变位点的多肽序列的规则为：

A、对于点突变：以发生突变的位置为中心，向前截取n个氨基酸，向后截取m个氨基酸，n＝m＝HLA分型所能呈递的最长肽段，若前段或后段长度不足时，则有多少截多少；如果点突变属于stop loss，则m为从当前突变到第一个终止密码子的长度；

B、对于非移码突变：非移码插入要从插入序列的第一个氨基酸向前截取n个氨基酸；从插入序列的最后一个氨基酸向后截取m氨基酸，n＝m＝HLA分型所能呈递的最长肽段；以插入序列为中心，插入序列、插入序列之前的n个氨基酸和插入序列之后的m个氨基酸共同组成新生多肽序列；

非移码缺失则以缺失位点为中心，分别向前截取n个氨基酸，向后截取m个氨基酸，n＝m＝HLA分型所能呈递的最长肽段；

C、对于移码突变：从开始移码突变的第一个氨基酸为中心，向前截取n个氨基酸，n＝HLA分型所能呈递的最长肽段；向后截取m个氨基酸，m为从当前突变到第一个终止密码子的长度，即向后截取至第一个终止密码子。