CN110322925A

CN110322925A - 一种预测融合基因产生新生抗原的方法

Info

Publication number: CN110322925A
Application number: CN201910649862.XA
Authority: CN
Inventors: 罗凯; 莫凡; 林志伟
Original assignee: Hangzhou New Ann Tianjin Biological Technology Co Ltd
Current assignee: Hangzhou New Ann Tianjin Biological Technology Co Ltd
Priority date: 2019-07-18
Filing date: 2019-07-18
Publication date: 2019-10-11
Anticipated expiration: 2039-07-18
Also published as: CN110322925B

Abstract

本发明公开了一种预测融合基因产生新生抗原的方法，包括如下内容：将测序得到的原始数据进行质控；将质控后的数据比对到参考基因组，并对重复序列进行标记；使用软件预测HLA一类、二类分型；使用软件预测融合基因事件；对预测的融合基因事件进行注释、打分、排序和筛选；对筛选后的融合基因断点上下游的核苷酸序列进行提取和翻译得到新生多肽；对新生多肽进行MHC亲和力预测，使用多个HLA分型亲和力预测软件进行亲和力预测，再在预测结果中筛选出亲和力较强的多肽；根据融合基因新生多肽的表位结合能力和多肽特性进行筛选和设计新生多肽；本发明能够在优化检测速度的同时，大幅提升检测目的基因融合新生抗原的灵敏度和特异性。

Description

一种预测融合基因产生新生抗原的方法

技术领域

本发明涉及基因检测技术领域，特别是一种利用生物信息学方法预测融合基因产生新生抗原的方法。

背景技术

恶性肿瘤的治疗面临诸多困境，亟需要新的治疗策略。近年来，肿瘤免疫治疗受到越来越多的关注并取得了突破性进展。肿瘤的免疫治疗是指利用机体自身的免疫***来清除肿瘤细胞，包括抗体治疗、细胞治疗、肿瘤疫苗等。

肿瘤的发生经常伴有多个基因的融合基因断点，新生抗原指的是由肿瘤细胞突变所产生的表位特异性抗原，只在肿瘤细胞上表达，从而不会导致机体的免疫耐受。已有较多的研究显示以新生抗原为靶标的免疫治疗，在一些癌症患者上已经取得了不错的临床效果。因此筛选鉴定出肿瘤特异性新生抗原是改善肿瘤免疫治疗效果的关键技术之一，也是实现个体化免疫治疗的基础。

融合基因(Fusion Gene)是由两个基因的部分序列拼接在一起形成的，其表达产物为融合转录本。融合基因作为染色体结构变异的一种产物，已被证明和某些癌症的产生有关，是研究某些癌症需要重点攻克的对象。融合基因由染色体变异导致，发生融合的基因往往是致癌基因。基因融合事件会产生新的融合基因，得到新的编码融合蛋白产物，从而引起了癌症的发生。当调控细胞增殖、分化和凋亡的基因发生融合，会直接影响下游信号传递途径，导致细胞增殖能力增强、凋亡障碍，分化障碍等，影响正常形状的表达，引发癌症。此外，基因融合会导致某些原癌基因与一个强启动子结合，受上游融合强启动子的刺激表达致癌活性。这会激活原癌基因或使抑癌基因失活，导致癌症的发生。致癌融合基因转录本也可以通过反式剪接或读段跨越事件引起。融合基因最早是在血液***恶性肿瘤中被发现。1960年，Peter Nowell和David Hungerford在慢性粒细胞白血病(CML)中发现一种异常染色体，即费城染色体。大约95％的该病患者染色体有异常：9号染色体上的原癌基因ABL易位到22号染色体上的断裂点形成了BCR-ABL融合基因。其表达产物为融合蛋白P-20，具有络氨酸激酶活性，会刺激白细胞增殖，导致白血病。

测序技术的飞速发展，尤其是二代测序技术的应用降低了检测基因融合的难度。其中末端配对(PE)测序的应用对成功检测融合基因尤为关键。基因融合的类型包括重复、倒位、覆盖和单碱基***缺失。全基因组测序及末端配对序列是目前能够最准确、最全面的检测所有基因融合类型的工具。全基因组测序在发现全新的基因融合断裂点方面更胜一筹，测序的深度和更长的测序序列实现了在融合连接单元单碱基的分辨率，为研究产生融合的机制提供了线索。

转录组测序是一种非常高效的方法，能为大量样品的筛查提供一种相对经济高效的途径。转录组研究基于高表达融合基因具有最大生物学影响的假设。尽管它对检测高表达癌基因特别有效，但很难获得基因间隔区域和非编码区域上的信息。目前大多数转录组分析比对的算法是假设转录组序列中的所有相关序列片段均来自同一个染色体，而基因融合的两个基因往往来自两个不同染色体，例如断裂-融合-桥接循环所产生的基因融合。对于这类情况，传统算法往往无法检测或正确拼接。

市场需要一种能够提高融合基因事件检测灵敏度和准确性的预测方法，本发明解决这样的问题。

发明内容

为解决现有技术的不足，本发明的目的在于提供一种预测融合基因产生新生抗原的方法，能够在优化检测速度的同时，大幅提升检测目的基因融合新生抗原的灵敏度和特异性。

为了实现上述目标，本发明采用如下的技术方案：

一种预测融合基因产生新生抗原的方法，包括如下内容：

一，将测序得到的原始数据进行质控；

二，将质控后的数据CleanReads比对到参考基因组，并对重复序列DuplicatedReads进行标记；

三，基于肿瘤HLA质谱数据库的机器学习训练HLA预测模型，并计算预测结果各HLA分型支持Reads数和HLA基因的表达，筛选得到准确的HLA一类分型和HLA二类分型；

四，将Duplicated Reads与参考基因组注释文件比对，根据split reads预测融合基因事件；

五，建立融合基因Meta分析预测模型，对预测的融合基因事件进行结构注释、打分、排序和筛选；

六，对筛选后的融合基因断点上下游的核苷酸序列进行提取和翻译得到新生多肽，所述新生多肽包括：具有融合基因断点的短肽，具有融合基因断点和突变组合的短肽；

七，对新生多肽进行MHC亲和力预测，使用多个HLA分型亲和力预测软件进行亲和力预测，预测结果会给出每一个HLA分型和每一条新生多肽的亲和力，筛选出其中亲和力较强的多肽，亲和力较强的标准是亲和力结合值≤500；

八，根据融合基因新生多肽的表位结合能力和多肽特性进行筛选和设计新生多肽，所述多肽特性包括：多肽的亲水性、疏水性、毒性以及二级结构。

前述的一种预测融合基因产生新生抗原的方法，

一，将测序得到的原始数据进行质控；

具体步骤包括：

1)去除接头污染的Reads，所述接头污染的Reads的碱基数大于5bp；

2)去除低质量的Reads，所述低质量的Reads为：质量值Q≤19的碱基占总碱基的50％以上的Reads；

3)去除含N比例大于5％的Reads；

4)舍弃去adapter及质量修剪后长度小于20bp的序列。

前述的一种预测融合基因产生新生抗原的方法，四，将Duplicated Reads与参考基因组注释文件比对，根据split reads预测融合基因事件；

即采用FusionCatacher、Star_Fusion和SOAPFuse 3个软件预测融合基因事件。

前述的一种预测融合基因产生新生抗原的方法，

对预测的融合基因事件进行注释，数据库为本地化的Chimerdb、Pancreases、TCGA和osmic注释信息；

注释是将预测得到的融合基因结果，基于融合基因断点附近Reads情况和数据库Chimerdb、Pancreases、TCGA和osmic进行注释；

注释内容包括：已知和新的融合基因，融合基因上下游基因的编码开放阅读框ORF是否发生改变；以及融合基因断点上下游两个基因的距离。

前述的一种预测融合基因产生新生抗原的方法，建立融合基因Meta分析预测模型的具体内容包括：采用多元线性回归模型，根据数据库注释融合基因事件类型、融合基因事件预测软件个数、融合基因上下游基因的编码开放阅读框ORF是否发生改变和融合基因断点上下游两个基因的距离四个因素进行打分；

具体打分公式如下：

Ei＝β1N1i+β2N2i+β3N1i+β4N2i+εi；

其中，Ei表示融合基因事件i的得分，N1i是融合基因i的数据库注释融合基因事件类型，β1是N1i对应的权重值，N2i是融合基因i的预测软件个数，β2是N2i对应的权重值，N3i是融合基因i对应的上下游基因的编码开放阅读框ORF是否发生改变状态，β3是N3i对应的权重值，N4i是融合基因i的融合基因断点上下游两个基因的距离，β4是N4i对应的权重值，εi是融合基因事件i的残差值。

前述的一种预测融合基因产生新生抗原的方法，

对预测的融合基因事件进行打分；

打分是用朴素贝叶斯算法建立一个融合基因测机器学习模型对预测基因融合结果进行打分；打分Score的具体计算公式如下：

N是本次分析预测融合基因事件总数，E是Meta分析预测模型对融合基因进行打分的得分，gene_Exp是基因表达，j和i分别代表融合基因断点上下游基因，m是支持融合基因事件的SpanReads数，n是支持融合基因事件的JunctionReads数，W是该因素分配的权限，Support SpanReads是SpanReads支持数；All SpanReads是可以比对到两个基因的转录组测序片段的数目；Support JunctionReads是JunctionReads支持数；All junctionReads是可以比对到推断的融合基因断点的转录组测序片段的数目；是融合基因事件发生支持SpanReads的频率；是融合基因事件发生支持JunctionReads的频率；

前述的一种预测融合基因产生新生抗原的方法，

对预测的融合基因事件进行排序和筛选；

根据融合基因测机器学习模型的打分进行排序，过滤reads覆盖度小于2的融合基因事件；过滤破坏外显子结果的融合基因事件；过滤不能翻译氨基酸序列的融合基因事件，筛选得到可信度较高的融合基因事件作为下游分析的输入；

融合基因筛选条件包括：SpanReads支持数≥2；JunctionReads支持数≥2；分隔成假双端读段的片段长度不小于25bp，并且不小于总长度的40％；split reads中融合位点坐标必须落在人工双端读段的中间部分，误差不超过4bp。

前述的一种预测融合基因产生新生抗原的方法，

具有融合基因断点的短肽的获取过程包括：以融合基因断点为中心、向前截取N个氨基酸和向后截取M个氨基酸获得多肽序列，N为HLA分型所能呈递的最大长度，M为HLA分型所能呈递的最大长度、或者M为从当前融合基因断点到第一个终止密码子的长度；在多肽序列依次截取长度为n的短肽，n为HLA分型结合需要的长度；得到包含当前融合基因断点的短肽；

具有融合基因断点和突变组合的短肽的获取过程包括：以融合基因断点为中心、向前截取N个氨基酸和向后截取M个氨基酸获得多肽序列，N为HLA分型所能呈递的最大长度，M为HLA分型所能呈递的最大长度、或者M为从当前融合基因断点到第一个终止密码子的长度；在多肽序列依次截取长度为n的短肽，n为HLA分型结合需要的长度；若多肽序列中除了当前融合基因断点以外，还包含其他突变，并且在n个短肽中具有同时包含融合基因断点和其他突变的短肽，则成为具有融合基因断点和突变组合的短肽。

前述的一种预测融合基因产生新生抗原的方法，

八，根据融合基因新生多肽的表位结合能力和多肽特性进行筛选和设计新生多肽，所述多肽特性包括：多肽的亲水性、疏水性、毒性以及二级结构；

筛选的标准包括：过滤高疏水性的多肽，过滤有毒性的多肽序列，选择高亲和力的融合基因新生抗原，选择高SpanReads支持数和JunctionReads支持数的融合基因新生抗原。

本发明的有益之处在于：

本发明先将测序得到的原始数据进行质控，保证了后续信息分析的质量；

本发明采用Meta分析方法，结合基于不同算法的融合基因预测软件，对预测结果进行汇总、打分，再根据最少reads支持准则以及校验工具进行进一步的矫正分析的结果，有效提高了融合基因事件预测的准确性；

本发明对预测的融合基因事件进行注释、打分、排序和筛选；通过SpanReads支持数、JunctionReads支持数、分隔成假双端读段的片段长度等因素，以及机器学习的手段，通过对预测结果的多因素分析，对融合基因事件进行打分；结果显示，这样的方法在不降低敏感度的同时能很好的控制假阳性率；

本发明发明设计了从二代测序数据获得融合基因新生多肽的流程，通过原始数据质控保证，HLA预测模型提高预测的灵敏度和准确性，提高融合基因预测灵敏度，打分模型进一步地过滤掉假阳性融合基因集合，从而设计出可精准的融合基因新生多肽，确保优选出具有临床价值的融合基因事件。

附图说明

图1是本发明的一种实施例的流程图；

图2是本发明的融合基因事件DCC-DYM进行一代测序PCR验证的结果图；

图3是本发明的融合基因事件MED13-BRIP1进行一代测序PCR验证的结果图；

图4是本发明的融合基因事件MAP2K4-SPOCK1进行一代测序PCR验证的结果图；

图5是本发明实验三中DCC-DYM的融合基因新生抗原的实验结果图；

图6是本发明实验三中MED13-BRIP1的融合基因新生抗原的实验结果图；

图7是本发明实验三中MAP2K4-SPOCK1的融合基因新生抗原的实验结果图；

图8是本发明实验三中NBEA-PDS5B的未入选融合基因位点突变型多肽的阴性对照组的实验结果图；

图9是本发明实验四中AGS(WILD-3)诱导的细胞毒性实验图；

图10是本发明实验四中AGS(MUT-3)诱导的细胞毒性杀伤实验图；

图11是本发明实验五的实验结果图；

图12是本发明实验六的实验结果图。

具体实施方式

以下结合附图和具体实施例对本发明作具体的介绍。

本发明中涉及的术语具体含义说明如下。

“参考基因组”，是指物种参考的标准基因组序列，一般是从官方网站下载(如：NCBI数据库等)，例如人参考基因组(GRCh38)等。

“Sanger测序”，是指双脱氧测序法，是通过一系列实验，从而获得可见的基因组碱基序列的方法。

“二代测序”，二代测序(Next-generation sequencing)是Sanger测序基础上通过技术创新产生的一种测序手段，核心思想是边合成边测序。

“PE测序”，是指双端测序，是一种测序方法。

“adaptor”，是指接头序列，为一段已知的短核苷酸序列，用于链接未知的目标测序片段。

“PCR”，是指多聚酶链式反应)，是体外酶促合成特异DNA片段的一种方法。

“Reads”，是指测序所得的基因组序列片段。

“CleanReads”，是指原始数据质控以后得到的序列片段。

“融合基因断点”，是指融合基因中的两个基因连接的点。

“HLA”，是指人类的主要组织相容性复合体(MHC)的表达产物。

“MHC”，是指一组编码动物主要组织相容性抗原的基因群的统称。

“Duplicated Reads”，是指重复序列，由于PCR扩增引入。

“STAR”，是指一种转录组主流的基因组比对软件，生产BAM格式文件。

“BAM”，是指基因组比对后产生的标准二进制格式文件。

“PICARD”，是指一种序列处理工具。

“net-MHCpan”，是指MHC亲和力预测分析的一种软件。

“SoapFuse”,是指融合基因事件预测分析的一种软件。

“TCGA”,是指目前最大的癌症基因信息的数据库(Cancer Genome Atlas)。

“Cosmic”,是指一个突变、融合基因和临床信息综合数据库。

“Chimerdb”,是指一个融合基因数据库。

如图1所示，一种预测融合基因产生新生抗原的方法，包括如下步骤：

步骤一，对转录组测序的原始数据进行质控；

具体内容包括：去除接头污染的Reads、去除低质量的Reads、去除含N比率大于5％的Reads、舍弃adapter及质量修剪后长度小于20bp的序列。

测序得到的某些原始下机序列，会含有测序接头序列以及低质量序列，为保证后续信息分析的质量，我们对原始序列进行过滤得到高质量的CleanReads，再进行后续分析。

作为一种实施例，具体如下：

1.去除接头污染的Reads(接头污染的Reads的碱基数大于5bp。对于双端测序，若一端受到接头污染，则去掉两端的Reads)；

2.去除低质量的Reads(Reads中质量值Q≤19的碱基占总碱基的50％以上，对于双端测序，若一端为低质量Reads，则会去掉两端Reads)；

3.去除含N比例大于5％的Reads(对于双端测序，若一端含N比例大于5％，则会去掉两端Reads)；

4.舍弃去adapter及质量修剪后长度小于20bp的序列。

步骤二，将质控后的数据比对到参考基因组；

具体内容包括：将CleanReads比对到参考基因组并进行排序和对DuplicatedReads进行标记。

质控后得到的CleanReads，通过软件STAR比对到参考基因组(例如人参考基因组)，使用两次比对模型，筛选出split reads(含有两个基因融合基因断点的reads)和discordant pair(比对到不同基因上的reads)作为候选的融合基因序列。PICARD标记Duplicated Reads，为后续的分析做准备。

步骤三，预测转录组数据的HLA(人类白细胞抗原)分型预测；

具体内容包括：使用软件预测一类和二类HLA分型。

基于肿瘤HLA质谱数据库的机器学习训练HLA预测模型，并计算预测结果各HLA分型支持Reads数和HLA基因的表达，筛选得到准确的HLA一类分型和HLA二类分型，HLA一类分型包括：HLA-A，HLA-B，HLA-C；HLA二类分型包括：HLA-DR，HLA-DQ，HLA-DP，为后续MHC亲和力预测做准备。

步骤四，在进行步骤三的同时，将候选融合基因与参考基因组注释文件比对，根据split reads预测融合基因事件；即采用FusionCatacher、Star_Fusion和SOAPFuse 3个软件预测融合基因事件。

步骤五，过滤融合基因事件；

具体内容包括：建立基因Meta分析预测模型，对预测的融合基因事件进行注释、打分、排序和筛选。

建立基因Meta分析预测模型的具体内容包括：采用多元线性回归模型，根据数据库注释融合基因事件类型、融合基因事件预测软件个数、融合基因上下游基因的编码开放阅读框ORF是否发生改变和融合基因断点上下游两个基因的距离四个因素进行打分，数据库为本地化的Chimerdb、Pancreases、TCGA和Cosmic注释信息。

具体打分公式如下：

Ei＝β1N1i+β2N2i+β3N1i+β4N2i+εi；

Meta分析采用R语言撰写程序实现。

注释是将预测得到的融合基因结果，基于融合基因断点附近Reads情况和数据库Cosmic和Chimerdb进行注释；

注释内容包括：已知和新的融合基因，融合基因上下游基因的编码开放阅读框ORF是否发生改变(in-Frame ORF未改变；Frame-shift ORF改变)；以及融合基因断点上下游两个基因的距离等；

打分是用朴素贝叶斯算法建立一个融合基因测机器学习模型对预测基因融合结果进行打分；

打分Score的具体计算公式如下：

融合基因筛选条件包括：SpanReads支持数≥2；JunctionReads支持数≥2；分隔成假双端读段的片段长度不小于25bp，并且不小于总长度的40％，保证足够长度的序列唯一的能比对到基因组上同时确保融合位点的发现；split reads中融合位点坐标必须落在人工双端读段的中间部分，误差不超过4bp。

筛选后再根据模型打分进行排序，过滤reads覆盖度小于2的融合基因事件；过滤破坏外显子结果的融合基因事件；过滤不能翻译氨基酸序列的融合基因事件，得到筛选出的可信度较高的融合基因事件作为下游分析的输入。本发明通过过对预测结果的多因素分析，对融合基因事件进行打分，有效的提高了预测结果的真阳性率。

步骤六，对筛选后融合基因断点上下游的核苷酸序列进行提取和翻译。

以融合基因断点为中心、向前截取N个氨基酸和向后截取M个氨基酸获得多肽序列，N为HLA分型所能呈递的最大长度，M为HLA分型所能呈递的最大长度、或者M为从当前融合基因断点到第一个终止密码子的长度；在多肽序列依次截取长度为n的短肽，n为HLA分型结合需要的长度；将包含当前融合基因断点的短肽用于步骤七；

以n＝11为例，将获得11个包含当前融合基因断点的短肽。如果多肽序列中除了当前融合基因断点以外，还包含其他突变，并且在n个短肽中具有某些短肽同时包含融合基因断点和其他突变，则成为具有融合基因断点和突变组合的短肽。

步骤七，新生多肽的MHC亲和力预测的步骤，该步骤包括对融合基因新生多肽的表位结合能力进行排序和打分。

分别使用多个HLA分型亲和力预测软件进行亲和力预测，多个HLA分型亲和力预测软件包括：针对HLA I型的netMHCpan，针对HLA II型的netMHC；软件的所有参数均为默认。预测结果会给出每一个HLA分型和每一条新生多肽的亲和力，筛选出其中亲和力较强的多肽，亲和力较强的标准是亲和力结合值≤500。

步骤八，对融合基因新生多肽的表位结合能力强弱和多肽特性进行筛选、设计。

筛选出有较强MHC亲和力的融合基因新生多肽，进一步的根据多肽的亲水性、疏水性、毒性以及二级结构综合考虑。

筛选的标准是：过滤高疏水性的多肽，避免多肽难合成；过滤有毒性的多肽序列；优选高亲和力的融合基因新生抗原；优选高SpanReads支持数和JunctionReads支持数的融合基因新生抗原；

多肽设计要求：限定长度在20-30bp；多肽合成委托有GMP资质的生产商。

以下通过实验验证本发明的效果；

制备实验材料；

实施例1

按照图一所示的流程，利用免疫组化法(IHC)和荧光原位杂交法(FISH)都验证3个融合基因事件(DCC-DYM、MED13-BRIP1和MAP2K4-SPOCK1)阳性的样本，进行高通量测序得到原始数据，然后作为本发明的方法的输入数据进行融合基因新生抗原的预测。具体的步骤和检测参数设定如下：

测序原始数据质控，去除接头序列污染碱基数大于5bp。对于双端测序，若一端受到接头污染，则去掉两端的Reads。去除低质量的Reads，Reads中质量值Q≤19的碱基占总碱基的50％以上，对于双端测序，若一端为低质量Reads，则会去掉两端Reads。去除含N比率大于5％的Reads，对于双端测序，若一端含N比例大于5％，则会去掉两端Reads。舍弃adapter及质量修剪后长度小于20bp的序列；质控后的数据比对到参考基因组，比对软件STAR将CleanReads比对到参考基因组，软件PICARD和SAMTOOLS进行排序和对Duplicated Reads进行标记；预测转录组数据的HLA(人类白细胞抗原)分型，使用软件预测一类和二类HLA分型，利用HLA丰度排除假阳性的HLA分型；同时进行融合基因事件预测，使用3个软件(FusionCatacher、Star_Fusion和SOAPFuse)预测融合基因事件；融合基因事件过滤，对预测的融合基因事件进行注释、打分、排序和筛选，支持融合基因事件的最低Junction Reads为2，支持融合基因事件的SpanReads为1；融合基因新生多肽序列提取，对筛选后融合基因断点上下游核苷酸序列的提取和翻译，筛除破坏外显子结构的融合基因事件，对于未改变开放阅读框的融合基因事件截取融合基因断点前后长度为48bp的核苷酸序列，对于改变开放阅读框的融合基因事件截取融合基因断点前长度为48bp的核苷酸序列和融合基因断点后的全部序列；融合基因新生多肽的MHC亲和力预测，对融合基因新生多肽的表位结合能力进行排序和打分，计算rank值；融合基因新生多肽设计，对融合基因新生多肽的表位结合能力强弱和多肽特性进行筛选、设计，最终得到融合基因产生的个体化新生抗原多肽。

实验一，用实施例1得到的个体化新生抗原多肽验证预测流程灵敏度和准确性；

检测所用的计算机资源峰值为4G，运行时长为15小时，质控结果如表1所示；得到质控结果的方法为：融合基因新生抗原生信预测使用上文的方法(不再赘述)和已有流程INTEGRATE Neo进行比较，使用服务器基本的资源管理命令即可得到资源和运行时间。

筛选后最终得到的融合基因检测结果如表2所示，测序是委托测序公司，生信分析步骤为上文的流程步骤，不再赘述。

评估效果如表3所示。

评估效果的方法为软件评估的一种基本方法，通过阳性数据对不同软件预测结果进行统计：

(1)真阳性(True Positive，TP)：检测有基因融合，且实际有基因融合；正确肯定的匹配数目；

(2)假阳性(False Positive，FP)：检测有基因融合，但实际无基因融合；误报，给出的匹配是不正确的；

(3)真阴性(True Negative，TN)：检测无基因融合，且实际无基因融合；正确拒绝的非匹配数目；

(4)假阴性(False Negative，FN)：检测无基因融合，但实际有基因融合；漏报，没有正确找到的匹配的数目。

正确率(Precision)：

真阳性率(True Positive Rate，TPR)，灵敏度(Sensitivity)，召回率(Recall)：

表1

表2

表3

对比实施例1，

使用INTEGRATE-Neo是实施例1中的原始数据进行默认参数的检测，检测所用的计算资源峰值为3.6G，运行时长18个小时，得到的评估效果如表4所示。

实验结果分析：从以上结果中可以看出，与对比例1相比，本发明实施例1的资源需求差不多，检测时长较现有方法更低，同时检测的敏感性和特异性也具有明显的优势。

实验二，代测序PCR验证融合基因事件真实性测试：

1.融合基因事件DCC-DYM进行一代测序PCR验证，使用的引物序列为：

F:5'TAGTGGGAGGAAGCAACTTACGGAT 3'

R:5'TGTGTGAAGGTACTTGTCTCGTGTC 3'

验证结果如图2所示。

2.融合基因事件MED13-BRIP1进行一代测序PCR验证，设计的反向引物序列，使用的引物序列为：

F:5'TGGAGAACTTTTAGAAACTTGTATCATT 3'

R:5'GCTGGAAAGAATAGTCATTGGAACC 3'

验证结果如图3所示。

3.融合基因事件MAP2K4-SPOCK1进行一代测序PCR验证，使用的引物序列为：

F:5'TATCAAACCTTCCAATATTCTTCTGGA 3'

R:5'CTCCTTGTCTGTGCAGGTGCCAT 3'

验证结果如图4所示。

结果分析：从以上PCR验证结果进一步的证明，本发明的融合基因检测和筛选的准确性高。

实验三，多肽的免疫原性测定实验；

实验目的：通过ELISpot实验，验证本发明中3条多肽在人源化小鼠体内均能引起免疫反应；并且，对比例1中的未入选融合基因位点NBEA-PDS5B几乎不能引起免疫反应。

实验方法：

为了检测多肽的免疫反应，实施IFN-γ酶联免疫吸附(ELISPOT)测定法。详细的实验过程如下：选用8周龄人源化小鼠B-NSG(CD34+)15只，随机分为5组，每组3只。适应一周后，分为阴性对照组编号1，3个单条融合基因新生多肽组编号3-5(共3组)、未入选融合基因位点NBEA-PDS5B对应多肽组编号6。采用CpG为佐剂(0.2μg/只)，多肽50μg每只,再与弗氏不完全佐剂Freund’s adjuvant(FIA，Sigma-Aldrich)1：1混匀，乳化30分钟，PBS与弗氏不完全佐剂1：1混合乳化30分钟作为阴性对照，四次于颈背部右胸皮下免疫，总剂量0.5mL/只，1周一次，共三周，第三次免疫后10天，取的小鼠脾脏，制备小鼠淋巴细胞悬液，用于ELISPOT检测。留取编号3-6部分小鼠淋巴细胞悬液用于实验3。

ELISPOT检测结果中，IFN-γ呈阳性结果的多肽，即判定为阳性候选多肽。将小鼠淋巴细胞稀释成浓度为1-2*106/mL，铺24孔板，每孔1mL，分为对照组(与多肽相同浓度的DMSO)、3个单条融合基因新生多肽组、对比例1中NBEA-PDS5B对应的多肽组，PHA阳性对照组编号2(PHA组淋巴细胞来源于阴性对照组)，共计6组，分别加入相应的多肽(10μg/mL)，预孵育72h后离心分离细胞，调整细胞浓度为2*106/mL，上IFN-γElispot板，按照IFN-γELISPOT试剂盒的说明书方法进行显色，运用CTL-ImmunoSpotS5系列酶联斑点分析仪读取产生的斑点数。IFN-γ阳性结果表明有抗原特异性T细胞产生，视为多肽能引起机体的免疫反应，斑点数的多少反映其免疫的强弱。

实验结果：

实验结果如图5-8所示；DCC-DYM、MED13-BRIP1和MAP2K4-SPOCK1分别代表每组3个单条融合基因新生多肽组注射小鼠，NBEA-PDS5B是未入选融合基因位点突变型多肽的阴性对照组注射小鼠，每只小鼠3个重复。每百万个细胞中，3条野生型多肽基本无斑点产生，3个单条融合基因新生多肽组都能产生约112-709个斑点，说明3条融合基因新生多肽单独都能引起免疫反应。对比例1中NBEA-PDS5B对应多肽组基本无斑点产生，说明该多肽几乎不能引起免疫反应。

实验四，融合基因新生抗原多肽疫苗有效性实验验证；

△实验4-1构建含特定突变位点的稳转细胞系；

实验目的：为了验证本发明中胃癌多肽疫苗的治疗及预防效果，因此需要构建一套含本发明中特定突变位点的稳转细胞系。

a.构建突变位点真核表达质粒

采用人工合成的方法获得能够表达本发明中突变的全部3条疫苗多肽的mini-gene，由以下几部分组成：信号肽部分(lysosome-associated membrane glycoprotein 1，LAMP1)，3条融合基因新生多肽和MHC class I trafficking domain(MITD)，3条融合基因新生多肽之间用柔性连接肽GGSGGGGSGG连接，基因进行密码子优化后，上游引入GATATC(EcoR V)，下游引入CTCGAG(Xho I)，将其克隆到真核表达载体pcDNA3.1-hygro(+)中，命名为MUT-3-pcDNA3.1(+)。同时合成相应野生型多肽作为对照，命名为WILD-3-pcDNA3.1(+)。为了验证靶的特异性，另设计了含有本发明中的1个未通过筛选标准的融合基因位点—NBEA-PDS5B的融合基因新生多肽，作为突变型多肽的阴性对照，命名为MUT-1-pcDNA3.1(+)。所有基因片段由南京金斯瑞生物科技有限公司代合成和构建，氨基酸序列如下所示：

其中第1-28个氨基酸为信号肽区域，用加粗斜体表示；下划线标注部分为MITD序列；

WILD-3氨基酸序列

MUT-3氨基酸序列

MUT-1氨基酸序列

b.构建能够稳定表达融合基因新生多肽的细胞系

胃癌细胞株系AGS购于ATCC。将人胃癌细胞株AGS以2*105/孔种于6孔板，待细胞覆盖70-80％时开始转染。分别将2.5μg MUT-3-pcDNA3.1(+)质粒、2.5μg WILD-3-pcDNA3.1(+)质粒和2.5μg MUT-1-pcDNA3.1(+)质粒稀释于3份100μl无血清RPMI-1640培养基中，再分别加入2.5μl PLUS^TM Reagents，室温孵育5min后，分别与含有5μl Lipofectamine^TM LTX的无血清RPMI-1640 100μl体系混合，室温孵育30min。将脂质体质粒络合物分别滴加于3份含有1000μl无血清RPMI-1640的待转染细胞中，前后轻轻摇匀，静置6h后更换为含10％血清的RPMI-1640培养基，继续培养48h后，换为含700μg/mL G418及400μg/mL潮霉素B的10％血清的培养基进行细胞筛选。抗性筛选10-14天，待对照组细胞全部死亡，转染组细胞大量生长，将细胞消化，采用有限稀释法种入96孔板中。显微镜下挑选单克隆细胞，用含700μg/mLG418，400μg/mL潮霉素B的培养基继续培养，隔天换液。继续培养约10天后，单克隆细胞长成较大一团，消化，转入24孔板培养。同时转染pcDNA3.1(+)-Hygro空载体质粒作同法抗性筛选，作为对照细胞，命名为AGS(含有pcDNA3.1(+))。上述方法中，用MUT-3-pcDNA3.1(+)质粒，WILD-3-pcDNA3.1(+)质粒和MUT-1-pcDNA3.1(+)质粒构建的细胞系分别命名为MUT-3-AGS(含有MUT-3-pcDNA3.1(+))、WILD-3-AGS(含有MUT-3-pcDNA3.1(+))，MUT-1-AGS(含有MUT-1-pcDNA3.1(+))。

△实验4-2：体外细胞杀伤实验；

实验目的：验证3条多肽在体外细胞水平能够引起肿瘤细胞的杀伤效果。

实验方法：

(1)5-(6)-Carboxy-fluorescein succinimidyl ester(CFSE)染料购自Invitrogen公司。操作步骤按照试剂盒说明书进行。无菌条件下用CFSE标记MUT-3-AGS(含有MUT-3-pcDNA3.1(+))和WILD-3-AGS(含有MUT-3-pcDNA3.1(+))靶细胞，分别作为实验组和对照组用的靶细胞。

(2)杀伤实验

a.准备效应细胞：

将实验2中留取的4组小鼠淋巴细胞悬液，用RPMI 1640培养基重悬，台盼蓝染色计数。

b.效应细胞CTL的诱导培养：

分别将4组淋巴细胞稀释浓度为(1-2)*106/mL，铺6孔板，每孔3mL，用RPMI1640+10％FBS+1×penicillin(100μg/mL)+streptomycin(100μg/mL)+1×MEM non-essentialamino acids+1mM sodium pyruvate+10mM HEPES buffer培养基进行培养，补加50U/mL的rhIL2。往每孔中加入10μg/mL对应的抗原肽，培养7天，每3天半量换液，并补加相应的抗原肽和rhIL2；一周之后细胞重悬，用PBS洗涤2次，制备成效应细胞CTL。

将MUT-3-AGS(含有MUT-3-pcDNA3.1(+))靶细胞分别与7种效应细胞CTL按照1:5、1:10和1:20的细胞数目比进行混合，加入到U型96孔板内，每孔体积200μL，作为实验组，每个实验组设置三个平行对照孔。

将WILD-3-AGS(含有MUT-3-pcDNA3.1(+))靶细胞分别与7种效应细胞CTL按照1:5、1:10和1:20的细胞数目比进行混合，加入到U型96孔板内，每孔体积200μL，作为对照组，每个对照组设置三个平行对照孔。

将96孔板放在37℃培养箱培养4h。

将96孔板离心去上清，用200μL预冷的PBS将细胞沉淀重悬，转到流式上样管中，用碘化丙啶(Prodium Iodide，PI)染色标记，浓度为1μg/m L，染色3min，马上进行流式上机检测。

实验结果和分析：

如图9所示，AGS(WILD-3)诱导的细胞毒性实验图；如图10所示为AGS(MUT-3)诱导的细胞毒性杀伤实验图；E:T代表效应T细胞个数:靶细胞个数，纵坐标表示细胞裂解率。DCC-DYM、MED13-BRIP1和MAP2K4-SPOCK1分别代表本发明中的3条多肽，Control代表阴性对照，NBEA-PDS5B是未入选的突变肽阴性对照组多肽。实验组诱导的效应T细胞其杀伤效率在(18.39±3.51)％-(67.71±2.77)％不等，杀伤效率明显高于对照组，说明融合基因新生多肽组对靶细胞起到了杀伤作用。融合基因新生多肽组中，随着效靶比的升高，T细胞的杀伤作用越来越强。该实验说明3条融合基因新生多肽能有效地激活特异性T细胞免疫应答。

实验五，人源化小鼠模型中本发明疫苗对肿瘤治疗效力的评价；

实验目的：验证3条多肽在体内能够激发T细胞对肿瘤细胞的杀伤。

实验过程：

选用8周龄人源化小鼠B-NSG(CD34+)12只，适应一周后，收集对数生长期的AGS(MUT-3)细胞，制备成5*10⁶/mL细胞悬液，以0.2mL种于小鼠左前肢腋下。以皮下肿瘤基因融合直径达5mm左右为成瘤标准，9-12天内成瘤，选择无出血、无坏死、无感染的小鼠，随机分为4组，每组3只，分3个单条融合基因新生多肽组和空白组(PBS)，分组之后当天免疫，采用CpG为佐剂(0.2μg/只)，多肽50μg每只，再与弗氏不完全佐剂Freund’s adjuvant(FIA，Sigma-Aldrich)1：1混匀，乳化30分钟，PBS与弗氏不完全佐剂1：1混合乳化30分钟作为阴性对照，四次于颈背部右胸皮下免疫，总剂量0.5mL/只，1周一次，免疫三次。第28天处死小鼠。小鼠肿瘤体积计算方法，用游标卡尺日测量肿瘤长径和短径，肿瘤体积计算公式如下；

TV＝1/2×a×b²，其中a、b分别表示长径和短径。

根据测量结果计算出相对肿瘤体积(relative tumor volume,RTV)，RTV＝Vt/V0。其中V0为分笼时(即d0)测量所得肿瘤体积，Vt为每次测量时的肿瘤体积，绘制相对肿瘤体积曲线，并记录各组动物体重。

实验结果和分析：

分组处理之前各组移植瘤的生长情况基本一致，体积大小无明显差异(P>0.05)。如图11所示，横坐标表示疫苗注射后的天数，纵坐标表示相对肿瘤体积。DCC-DYM、MED13-BRIP1和MAP2K4-SPOCK1分别代表本发明中的3条多肽，Control代表阴性对照，NBEA-PDS5B是未入选的突变肽阴性对照组多肽。免疫后1周，各组移植瘤生长与对照组相比，其体积大小无明显差异(P>0.05)，随着免疫次数增多和时间的延长，实验组与对照组相比，瘤体大小受抑程度越来越明显，说明3条多肽均能在小鼠体内起到抑制肿瘤生长的作用。免疫后2周，其肿瘤已出现缩小的情况，说明融合基因新生多肽疫苗激发的免疫响应对肿瘤细胞产生了杀伤作用。

实验六，人源化小鼠模型中本发明疫苗对肿瘤预防效力的评价；

实验目的：验证3条多肽在人源化小鼠上安全有效且能有效预防肿瘤的发生。

实验过程：

选用7周龄人源化小鼠B-NSG(CD34+)12只，随机分为4组，每组3只，适应一周后，采用CpG为佐剂(0.2μg/只)，多肽50μg每只，再与弗氏不完全佐剂Freund’s adjuvant(FIA，Sigma-Aldrich)1：1混匀，乳化30分钟，PBS与弗氏不完全佐剂1：1混合乳化30分钟作为阴性对照，分3个单条融合基因新生多肽组和空白组(PBS)，四次于颈背部右胸皮下免疫，总剂量0.5mL/只，3-5天一次，免疫三次。第三次免疫后1周，收集对数生长期的AGS(MUT-3)细胞，制备成5*106个/mL细胞悬液，以0.2mL种于小鼠左前肢腋下。

每天记录各组动物的肿瘤发生情况，肿瘤长出来后，用游标卡尺日测量肿瘤长径和短径，计算肿瘤体积：TV＝1/2×a×b²，其中a、b分别表示长径和短径，并记录各组动物体重。评价：接种肿瘤四周后，处死所有小鼠，终止实验，计算抑瘤率并进行评价。

抑瘤率(％)＝(对照组肿瘤平均体积-实验组肿瘤平均体积)/对照组肿瘤平均体积×100％。

数据处理：采用R语言进行统计学分析，采用单因素方差分析进行比较，P<0.05为差异有统计学意义。

实验结果：

如图12所示，横坐标表示疫苗注射后的天数，纵坐标表示肿瘤体积。DCC-DYM、MED13-BRIP1和MAP2K4-SPOCK1分别代表本发明中的3条多肽，Control代表阴性对照，NBEA-PDS5B是未入选的突变肽阴性对照组多肽。接种肿瘤10天之后，对照组肿瘤开始缓慢增大，各实验组肿瘤体积抑制明显。与对照组相比，第16天时，各实验组仍无肿瘤形成。随着时间的延长，各实验组瘤体大小缓慢增大，但增长速度明显低于对照组。

结果分析：本试验证实了3条融合基因新生多肽能预防特异肿瘤发生。

本发明通过原始数据质控保证，HLA预测模型提高预测的灵敏度和准确性，提高融合基因预测灵敏度，打分模型进一步地过滤掉假阳性融合基因集合，从而设计出可精准的融合基因新生多肽。

以上显示和描述了本发明的基本原理、主要特征和优点。本行业的技术人员应该了解，上述实施例不以任何形式限制本发明，凡采用等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。

序列表

<110> 杭州纽安津生物科技有限公司

<120> 一种预测融合基因产生新生抗原的方法

<141> 2019-07-17

<160> 9

<170> SIPOSequenceListing 1.0

<210> 1

<211> 25

<212> DNA

<213> Artificial Sequence

<400> 1

tagtgggagg aagcaactta cggat 25

<210> 2

<211> 25

<212> DNA

<213> Artificial Sequence

<400> 2

tgtgtgaagg tacttgtctc gtgtc 25

<210> 3

<211> 28

<212> DNA

<213> Artificial Sequence

<400> 3

tggagaactt ttagaaactt gtatcatt 28

<210> 4

<211> 25

<212> DNA

<213> Artificial Sequence

<400> 4

gctggaaaga atagtcattg gaacc 25

<210> 5

<211> 27

<212> DNA

<213> Artificial Sequence

<400> 5

tatcaaacct tccaatattc ttctgga 27

<210> 6

<211> 23

<212> DNA

<213> Artificial Sequence

<400> 6

ctccttgtct gtgcaggtgc cat 23

<210> 7

<211> 209

<212> PRT

<213> Aminoacid sequence

<400> 7

Met Ala Ala Pro Gly Ser Ala Arg Arg Pro Leu Leu Leu Leu Leu Leu

1 5 10 15

Leu Leu Leu Leu Gly Leu Met His Cys Ala Ser Ala Glu Ala Gly Asn

20 25 30

Ala Gln Thr Ser Ala Gln Leu Ile Val Pro Lys Pro Ala Ile Pro Ser

35 40 45

Ser Ser Val Leu Pro Ser Ala Pro Arg Asp Val Val Gly Gly Ser Gly

50 55 60

Gly Gly Gly Ser Gly Gly Tyr Gly Glu Leu Leu Glu Thr Cys Ile Ile

65 70 75 80

Asn Ile Asp Val Pro Asn Arg Ala Arg Arg Lys Lys Ser Ser Ala Arg

85 90 95

Lys Phe Gly Leu Gln Lys Gly Gly Ser Gly Gly Gly Gly Ser Gly Gly

100 105 110

Asp Ser Ile Ala Lys Thr Arg Asp Ala Gly Cys Arg Pro Tyr Met Ala

115 120 125

Pro Glu Arg Ile Asp Pro Ser Ala Ser Arg Gln Gly Tyr Asp Val Arg

130 135 140

Gly Gly Ser Gly Gly Gly Gly Ser Gly Gly Ile Val Gly Ile Val Ala

145 150 155 160

Gly Leu Ala Val Leu Ala Val Val Val Ile Gly Ala Val Val Ala Thr

165 170 175

Val Met Cys Arg Arg Lys Ser Ser Gly Gly Lys Gly Gly Ser Tyr Ser

180 185 190

Gln Ala Ala Ser Ser Asp Ser Ala Gln Gly Ser Asp Val Ser Leu Thr

195 200 205

Ala

<210> 8

<211> 243

<212> PRT

<213> Aminoacid sequence

<400> 8

Met Ala Ala Pro Gly Ser Ala Arg Arg Pro Leu Leu Leu Leu Leu Leu

1 5 10 15

Leu Leu Leu Leu Gly Leu Met His Cys Ala Ser Ala Glu Ala Gly Asn

20 25 30

Ala Gln Thr Ser Ala Gln Leu Ile Val Pro Lys Pro Asp Thr Lys Lys

35 40 45

Tyr Tyr Leu Val Phe Arg Thr Ser Phe Asn Gly Gly Ser Gly Gly Gly

50 55 60

Gly Ser Gly Gly Tyr Gly Glu Leu Leu Glu Thr Cys Ile Ile Asn Ile

65 70 75 80

Asp Val Pro Asn Ser Pro Leu Ala Thr Val Leu Gly Ala Val Val Leu

85 90 95

Leu Asn Lys Glu Thr Val Lys Ser His Arg Ile Pro Leu Arg Arg Ile

100 105 110

Ile Gln Gly Asn Pro Arg Tyr Pro Lys Tyr Ile Leu Gly His Ala His

115 120 125

Thr Ser Arg Leu Leu Arg Leu Leu Glu Ser Ser Gly Gly Arg His Ile

130 135 140

Gln Gly Phe Gln Gly Gly Ser Gly Gly Gly Gly Ser Gly Gly Asp Ser

145 150 155 160

Ile Ala Lys Thr Arg Asp Ala Gly Cys Arg Pro Tyr Met Ala Pro Ala

165 170 175

Gln Thr Arg Ser Gly Gly Ser Gly Gly Gly Gly Ser Gly Gly Ile Val

180 185 190

Gly Ile Val Ala Gly Leu Ala Val Leu Ala Val Val Val Ile Gly Ala

195 200 205

Val Val Ala Thr Val Met Arg Arg Lys Ser Ser Gly Gly Lys Gly Gly

210 215 220

Tyr Ser Gln Ala Ala Ser Ser Asp Ser Ala Gln Gly Ser Asp Val Ser

225 230 235 240

Leu Thr Ala

<210> 9

<211> 123

<212> PRT

<213> Aminoacid sequence

<400> 9

Met Ala Ala Pro Gly Ser Ala Arg Arg Pro Leu Leu Leu Leu Leu Leu

1 5 10 15

Leu Leu Leu Leu Gly Leu Met His Cys Ala Ser Ala Gly Lys Gly Phe

20 25 30

Gln His Cys Val Lys Tyr Asp Phe Gln Pro Arg Lys Asp Ile Phe Met

35 40 45

Phe Ile Thr Arg Gln Leu Lys Gly Leu Glu Asp Thr Gly Gly Ser Gly

50 55 60

Gly Gly Gly Ser Gly Gly Ile Val Gly Ile Val Ala Gly Leu Ala Val

65 70 75 80

Leu Ala Val Val Val Ile Gly Ala Val Val Ala Thr Val Met Arg Arg

85 90 95

Lys Ser Ser Gly Gly Lys Gly Gly Tyr Ser Gln Ala Ala Ser Ser Asp

100 105 110

Ser Ala Gln Gly Ser Asp Val Ser Leu Thr Ala

115 120

Claims

1.一种预测融合基因产生新生抗原的方法，其特征在于，包括如下内容：

一，将测序得到的原始数据进行质控；

二，将质控后的数据CleanReads比对到参考基因组，并对重复序列Duplicated Reads进行标记；

2.根据权利要求1所述的一种预测融合基因产生新生抗原的方法，其特征在于，

一，将测序得到的原始数据进行质控；

具体步骤包括：

3)去除含N比例大于5％的Reads；

4)舍弃去adapter及质量修剪后长度小于20bp的序列。

3.根据权利要求1所述的一种预测融合基因产生新生抗原的方法，其特征在于，四，将Duplicated Reads与参考基因组注释文件比对，根据split reads预测融合基因事件；

即采用FusionCatacher、Star_Fusion和SOAPFuse3个软件预测融合基因事件。

4.根据权利要求1所述的一种预测融合基因产生新生抗原的方法，其特征在于，

5.根据权利要求4所述的一种预测融合基因产生新生抗原的方法，其特征在于，所述建立融合基因Meta分析预测模型的具体内容包括：采用多元线性回归模型，根据数据库注释融合基因事件类型、融合基因事件预测软件个数、融合基因上下游基因的编码开放阅读框ORF是否发生改变和融合基因断点上下游两个基因的距离四个因素进行打分；

具体打分公式如下：

Ei＝β1N1i+β2N2i+β3N1i+β4N2i+εi；

6.根据权利要求5所述的一种预测融合基因产生新生抗原的方法，其特征在于，

对预测的融合基因事件进行打分；

N是本次分析预测融合基因事件总数，E是Meta分析预测模型对融合基因进行打分的得分，gene_Exp是基因表达，j和i分别代表融合基因断点上下游基因，m是支持融合基因事件的SpanReads数，n是支持融合基因事件的JunctionReads数，W是该因素分配的权限，Support SpanReads是SpanReads支持数；AllSpanReads是可以比对到两个基因的转录组测序片段的数目；Support JunctionReads是JunctionReads支持数；All junctionReads是可以比对到推断的融合基因断点的转录组测序片段的数目；是融合基因事件发生支持SpanReads的频率；是融合基因事件发生支持JunctionReads的频率。

7.根据权利要求6所述的一种预测融合基因产生新生抗原的方法，其特征在于，

对预测的融合基因事件进行排序和筛选；

8.根据权利要求1所述的一种预测融合基因产生新生抗原的方法，其特征在于，

9.根据权利要求1所述的一种预测融合基因产生新生抗原的方法，其特征在于，