CN111415707B

CN111415707B - 临床个体化肿瘤新抗原的预测方法

Info

Publication number: CN111415707B
Application number: CN202010162494.9A
Authority: CN
Inventors: 许恒; 舒洋; 杨莉; 丁振宇; 魏于全
Original assignee: Sichuan University
Current assignee: Sichuan University
Priority date: 2020-03-10
Filing date: 2020-03-10
Publication date: 2023-04-25
Anticipated expiration: 2040-03-10
Also published as: CN111415707A

Abstract

本发明涉及预测技术，解决了现有肿瘤新抗原分析无法采用一步法获得最后结果及根据临床结果的反馈提升预测的准确率的缺点，提供了一种临床个体化肿瘤新抗原的预测方法，其技术方案可概括为：将患者肿瘤细胞的测序数据与对应参考基因组进行比对得到DNA及RNA比对结果，进行预处理后分析肿瘤细胞中的体细胞变异、克隆类型、肿瘤纯度、肿瘤细胞中的基因表达量及体细胞变异等位基因表达丰度，得到分析结果，再分析潜在新抗原的亲和性及计算多肽递呈效力，根据分析结果、潜在新抗原的亲和性及多肽递呈效力对各对应的新抗原进行评分及排序后呈现。其有益效果是，标准化的过程可提高准确性，且具有更便捷的操作性及可重复性，适用于临床个体化肿瘤新抗原的预测。

Description

临床个体化肿瘤新抗原的预测方法

技术领域

本发明涉及预测技术，特别涉及个体化肿瘤新抗原的预测技术。

背景技术

肿瘤的免疫治疗被认为是继手术、放疗、化疗和小分子靶向治疗之后的新一代肿瘤治疗方法，它与之前传统治疗方法的不同之处在于：免疫治疗将治疗手段从直接药物杀伤肿瘤细胞的思路转变为激活免疫细胞，通过增强患者自身的抗肿瘤免疫能力来***，与传统治疗相比具有杀伤精准、副作用小、疗效持久及个性化程度高等优势。此外，机体免疫***具有免疫记忆的特性，因此免疫治疗可以帮助患者形成记忆型免疫，在防止肿瘤复发和转移上具有显著优势。

目前免疫治疗主要分为两类：基于免疫检查点抑制剂的免疫治疗和细胞免疫治疗。免疫检查点抑制剂已经有多种上市的药物(包括施贵宝和默克公司的PD-1抑制剂等)，在晚期肿瘤的治疗中已经取得了显著的成效，特别是很多对传统疗法耐受的患者也能通过这种疗法获益。对部分肿瘤的治疗中该疗法已经作为一线治疗手段。但是该疗法的总体获益人群占总患者人群的比例约为10～30％，代表着仍然有大量的患者不能获益。

细胞免疫治疗又称细胞过继免疫治疗，主要将病人的自身T细胞进行外界修饰或者刺激，让普通T细胞能够特异识别肿瘤细胞，从而对其进行攻击和杀伤。而作为细胞免疫治疗的代表，CAR-T已在血液肿瘤及部分实体瘤中取得了巨大的成功。然而对于实体瘤而言，免疫细胞治疗仍面临着很大的挑战，需要研发新的更有效的免疫治疗关键技术。T细胞特异识别肿瘤细胞的机制在于：肿瘤细胞中存在区别于正常细胞的基因突变，这些突变被翻译成蛋白质后有异于正常的蛋白，其中一些突变蛋白经过泛素化修饰和降解后形成的多肽能被免疫***识别，称为多肽的抗原性，这些具有抗原性的多肽被肿瘤细胞自身的主要组织相容性复合体(MHC)递呈到细胞表面，继而被免疫***识别为外来物，使得肿瘤细胞被免疫细胞攻击并杀伤。这类具有抗原性的突变蛋白也被称为叫肿瘤新抗原(neoantigen)。细胞过继免疫治疗根据其发展历程主要分为以下几类：自体淋巴因子激活的杀伤细胞(LAK)、自体肿瘤浸润性淋巴细胞(Tumor Infiltrating Lymphocytes,TILs)、树突状细胞结合细胞因子诱导的杀伤细胞(DC-CIK)、以及经基因修饰改造的T细胞(CAR-T)。CAR-T疗法通过修饰设计了特异针对某种抗原，目前在治疗B淋巴细胞白血病和淋巴瘤有较好的效果，但对实体瘤效果较差，主要是因为其能利用的肿瘤抗原较少且靶点单一。除CAR-T外，其他的细胞治疗均没有特殊设计的抗原，需要免疫细胞通过自然反应识别肿瘤的特殊抗原，因此这些疗法的受益人群并不多(约10％)，主要是因为肿瘤自身的抗原不足以激活免疫细胞。近期有临床研究发现将肿瘤细胞裂解物与免疫细胞进行共培养可以提高治疗效果，同时在患者体内检测到了特异识别肿瘤抗原的特殊T细胞，证明免疫***通过识别肿瘤抗原来攻击肿瘤细胞。与此类似的是，美国国家癌症研究所外科分部的主任Steven Rosenberg博士从切除的肿瘤或活检样本中提取出肿瘤浸润T淋巴细胞(这类细胞来源于肿瘤，是已知的能够识别和攻击癌症的免疫细胞)，挑选出抗肿瘤活性最强的细胞，进行体外增殖并输送回患者体内，该治疗方法取得了很好的疗效，在少数患者中甚至能使肿瘤完全消退。然而，这些方法都比较被动，因为很多患者的肿瘤细胞还不足以激活自身免疫***，或者淋巴浸润T细胞数量太少达不到最终治疗要求，因此有效提高细胞治疗疗效的一个策略应当是主动寻找肿瘤新抗原，通过体外培养对免疫细胞进行高强度的抗原刺激，再将能特异识别肿瘤新抗原的细胞回输到患者体内用于攻击具有肿瘤新抗原的肿瘤细胞。在临床操作中，研究人员通常从肿瘤患者中抽取血细胞培养出成熟的DC细胞，然后用“突变抗原”刺激DC细胞，制备成DC疫苗。美国华盛顿大学的研究表明，使用该类DC疫苗可在晚期黑色素瘤患者的体内激发出有效的免疫反应，产生精准识别癌细胞的具有高杀伤活性的T细胞。

国际大型的肿瘤基因组计划(TCGA)已经找到了大量不同肿瘤的体细胞突变并绘制成突变图谱，目前已经有研究开展了针对突变图谱中新抗原的分析和验证。不过值得特别注意的是，由于东西方人种之间的差异，以及肿瘤本身巨大的异质性的特点，基于欧美患者获得的肿瘤突变谱在我国病人中有可能不突变或是突变频率较低，而在我国属于高频突变的位点在西方人群中则可能未被报道；此外，决定突变是否有免疫原性的关键决定因素HLA位点变异在不同人种中也存在着巨大差异，提示我们需要对中国人群进行人种特异性的研究和分析。

要主动找到肿瘤新抗原，必须具备两个条件：第一是能对每个患者个体肿瘤基因组进行高通量的DNA序列检测，以寻找所有的肿瘤特异性突变；另一方面则是对这些突变表达出来的蛋白和多肽进行抗原性预测，将潜在的新抗原肽段进行体外合成(这些肽段也被称为肿瘤疫苗)，最后使用高浓度和高强度的肿瘤疫苗激活免疫细胞。目前，针对肿瘤的新抗原免疫治疗已经有较好的理论和实践基础，主要体现在以下几方面：

1、随着二代测序技术和生物信息学的发展，短时间大规模的测序和分析已经成为可能，目前对人基因组的全扫描只需要一周的时间，特别是还可以仅针对基因编码区进行全外显子测序，大大减少了测序量和分析的难度，使鉴定肿瘤内部所有突变的过程变得越来越简单。另一方面，针对肿瘤体细胞检测的软件已经非常成熟，能高效快速地辨别出肿瘤组织特有的体细胞突变。此外，肿瘤新抗原预测程序和软件的成功开发和应用也为我们寻找新抗原提供有力的工具。这些成熟的生物信息学分析方法和软件使得我们可以在一天内即完成病人的体细胞变异检测、HLA基因型预测、突变所在RNA的丰度计算、肿瘤克隆状态分析以及新抗原亲和性预测等工作，为临床治疗节约宝贵的时间。

2、科学家长期以来立足于免疫学的基础研究，其中关于T细胞抗原识别的研究已经取得了长足的进展。研究发现在肿瘤中的某些变异蛋白在经过降解后形成的多肽(也被称为抗原表位)可以被MHC分子特异识别并形成MHC-抗原肽复合体，继而复合体被呈递到抗原递呈细胞的表面并被T细胞识别为非己成分，使得T细胞活化成效应T细胞而对肿瘤细胞进行攻击和清除，其中最为关键的步骤就是MHC分子与多肽分子的结合。人体内的MHC分子又称为人类白细胞抗原(HLA)，主要分为MHC I类分子和MHC II类分子，前者由大多数的细胞表达，主要参与内源抗原的递呈；而后者主要由抗原递呈细胞表达，主要参与外源抗原的递呈。其中人的基因组中存在三个编码MHC I(HLA-A、HLA-B和HLA-C)和三个MHC II的基因(HLA-DR、HLA-DQ、HLA-DP)，这些基因的变异非常大，目前已经发现有数万种不一样的基因型，不同的基因型一般对同一多肽有不同的亲和力，也就是识别其作为抗原的能力，亲和力越强则识别能力越强。

3、截止目前为止，相关领域的科学家已经构建了多个抗原表位数据库并开发了多个抗原表位预测相关软件。

IEDB数据库是由Peter等人开发能为用户提供预测以及相关注释的服务，包括抗原表位结合预测，基于细胞内存在的肽对抗原表位进行预测，以及对多肽/MHC复合物引发免疫应答的相对能力进行预测等等。目前抗原预测的主流方法都在预测MHC分子结合的抗原肽，根据预测时使用的特征可将预测方法分为三类：基于序列信息的预测方法，基于结构信息的预测方法，以及结合序列信息和结构信息的预测方法。MHC I、II类分子结合抗原肽的结构位于该分子远膜端的抗原结合槽。其中，MHC I类分子凹槽封闭，与之结合的抗原肽长度为8～10个氨基酸残基，而MHC II类分子凹槽开放，与之结合的抗原长度比MHC I类分子更长，为13～17个氨基酸残基。因而与MHC I类分子结合的抗原肽比与MHC II类分子结合的抗原肽更具有特定的锚定基。由于目前MHC I类分子和MHC II类分子结合肽的预测方法没有很大差别，所以预测MHC I类分子结合肽的结果会更加准确。基于序列信息的预测方法是因为相似的序列局有相似的结构和相同的功能。最早的预测方法就是经过统计序列信息来预测抗原表位，通过锚定基的氨基酸的出现频率构建序列模体。后来发现除了锚定基的氨基酸之外，其他位置的氨基酸也对MHC分子结合的抗原肽有影响。因此考虑结合肽的所有位置上出现20种氨基酸的频率，构建出L X 20的矩阵模型。研究发现氨基酸的物化性质和相邻氨基酸的相关系数都能为预测提供有效特征，此外随着计算科学的发展，机器学习(ANN,HMM)和线性规划(SMM)等预测方法被引入到抗原预测。

4、最新的四个独立的研究表明，用结合序列信息和结构信息的预测方法预测出来的新抗原可以作为肿瘤疫苗对患者进行治疗，并取得了很好的效果。

这些成功的案例提示利用新抗原进行的免疫治疗具有很高的应用前景。然而目前该方法的个体化非常强，从测序到最终新抗原选择需要对每个个体都进行步骤繁琐和复杂的分析流程，特别是一定要有遗传学和生物信息学专家的参与。在报道的临床研究中，由于纳入的患者人数较少，这样的操作是可以接受的。但是如果这样成功的治疗方案一旦进行常规临床治疗，绝大多数医院都缺乏负责新抗原选择的专业人员，不但因为新抗原的选择延误治疗时间，而且在不同的医院很难进行标准化，可能在最终治疗效果上造成极大的差异。因此，基于新抗原的免疫治疗的大规模推广需要一个可标准化的流程，如采用一步法的软件和分析流程。另一方面，目前绝大多数的临床研究都是基于欧美人群的，由于前述的人种差异的存在，还需要对新抗原预测过程和分析进行一个人种特异性的优化过程。

发明内容

本发明的目的是要克服目前肿瘤新抗原分析无法采用一步法获得最后结果及根据临床结果的反馈提升预测的准确率的缺陷，提供一种临床个体化肿瘤新抗原的预测方法。

本发明解决上述技术问题，采用的技术方案是，临床个体化肿瘤新抗原的预测方法，包括以下步骤：

步骤1、将患者肿瘤细胞的测序数据与预设的对应参考基因组进行比对，得到DNA比对结果及RNA比对结果；

步骤2、对DNA比对结果及RNA比对结果进行预处理；

步骤3、利用DNA比对结果及RNA比对结果分析肿瘤细胞中的体细胞变异、克隆类型、肿瘤纯度、肿瘤细胞中的基因表达量及体细胞变异等位基因表达丰度，得到分析结果；

步骤4、预测患者肿瘤细胞的MHC I类分子的基因型，并根据DNA分析结果及患者肿瘤细胞的MHC I类分子的基因型分析潜在新抗原的亲和性及计算多肽递呈效力；

步骤5、根据分析结果、潜在新抗原的亲和性及多肽递呈效力对各对应的新抗原进行评分及排序后呈现。

具体的，步骤1中，所述将患者肿瘤细胞的测序数据与预设的对应参考基因组进行比对中，采用Burrows-Wheeler转化算法得到DNA比对结果；采用STAR算法得到RNA比对结果。

进一步的，步骤1之前，还包括以下步骤：

步骤0、对患者肿瘤细胞的原始测序数据进行预处理，得到患者肿瘤细胞的测序数据。

具体的，步骤0中，所述预处理包括修剪测序接头序列、修剪测序标签序列及去除测序质量较差的序列，其可采用Trimmomatic软件进行。

再进一步的，步骤2中，所述预处理包括去重、添加组信息及碱基质量重计算，其可采用GATK软件进行。

具体的，步骤3包括以下步骤：

步骤3A、利用DNA对比结果分析肿瘤细胞中的体细胞突变位点；

步骤3B、根据体细胞突变位点结合RNA比对结果分析得到体细胞变异等位基因表达丰度及肿瘤细胞中的基因表达量；

步骤3C、通过体细胞变异等位基因表达丰度过滤体细胞突变位点中的假阳性位点，得到实际体细胞突变位点；

步骤3D、对实际体细胞突变位点进行注释及进一步过滤，得到体细胞变异；

步骤3E、计算肿瘤的拷贝数变异文件；

步骤3F、根据体细胞变异及肿瘤的拷贝数变异文件计算肿瘤纯度及克隆结构；

步骤3G、根据RNA比对结果分析肿瘤细胞中的基因表达量。

再进一步的，步骤3D中，所述体细胞变异是指：针对注释后的实际体细胞突变位点进行进一步过滤，仅保留外显子区域中非同义体细胞突变位点。

具体的，步骤3A采用Mutect2软件进行分析；步骤3B采用GATK的HaplotyperCaller软件分析；步骤3C采用FilterMutectCalls软件进行分析；步骤3D采用ANNOVAR软件进行注释；步骤3E采用cnvkit软件进行计算；步骤3F中，所述计算肿瘤纯度采用ABSOULTE软件进行计算，所述克隆结构采用pyClone软件进行计算；所述克隆结构是指突变基因克隆比，即每个体细胞变异位点的克隆比例。

再进一步的，步骤3G中，基因表达量的分析是指采用featureCounts软件进行分析得到counts数值，并将所得到的counts数值转换成每百万测序片段中转录本测序片段数(TPM数值)。

具体的，所述将所得到的counts数值转换成每百万测序片段中转录本测序片段数中，其转换公式为:

其中，R为基因的count数目，L为基因的长度，

R_k为基因k的count数目，L_k为基因k的长度，n为基因的总数量。

再进一步的，步骤4包括以下步骤：

步骤4A、预测患者肿瘤细胞的MHC I类分子的基因型；

步骤4B、根据注释后的实际体细胞突变位点获取非同义突变蛋白数据；

步骤4C、根据非同义突变蛋白数据预测MHC I类分子结合野生型多肽和变异型多肽的能力，得到亲和力预测；

步骤4D、分析多肽序列，得到多肽解离效率、转运效力及MHC结合效力，并根据所得到的多肽解离效率、转运效力及MHC结合效力计算多肽递呈效力；

步骤4E、根据亲和力预测及非同义突变蛋白数据得到潜在新抗原的亲和性。

具体的，步骤4A采用polsover软件进行分析；步骤4C采用IEDB提供的MHC class Iepitope predictors软件进行分析；步骤4D采用netCTLpan软件进行分析计算。

再进一步的，步骤4B具体为：下载UCSC网站提供的人基因组蛋白序列fasta文件，根据注释后的实际体细胞突变位点提取出变异氨基酸前后各15个氨基酸的序列，分别获得野生型蛋白质和突变型蛋白质的侧翼氨基酸序列，将所有的序列组合成多肽序列的fasta文件，即非同义突变蛋白数据。

具体的，步骤5中，所述根据分析结果、潜在新抗原的亲和性及多肽递呈效力对各对应的新抗原进行评分的方法为：

score＝abundance*dissimilarity*clonality

其中，

abundance＝L(IC_m)*Freq*tanh(TPM)；

dissimilarity＝(1–L(IC_w/50000))/2；

clonality＝NCTL*MC；

L(x)＝1/(1+e^x),tanh()为双曲正切函数计算，score是新抗原计算获得的分数，abundance是突变多肽的丰度,TPM为基因表达丰度，Dissimilarity是不相似度(dissimilarity)衡量突变肽段与对应正常肽段亲和力的差异，克隆性(clonality)衡量新抗原由突变至短肽的效率NCTL及新抗原在所有肿瘤细胞中分布的比例MC，Freq为体细胞变异等位基因表达丰度比例，所述体细胞变异等位基因表达丰度比例是指体细胞变异等位基因表达丰度与肿瘤细胞中的基因表达量的比值，IC_m为亲和力预测时得到的突变多肽亲和力，IC_w为亲和力预测时得到的野生多肽结合亲和力,NCTL为递呈效力，MC为突变基因克隆比。

具体的，还包括以下步骤：

步骤6、基于实验结果，利用机器学习方法优化步骤4的潜在新抗原的亲和性的分析方法。

本发明的有益效果是：在本发明方案中，采用上述临床个体化肿瘤新抗原的预测方法，可见，其利用已有的较成熟的高通量分析软件建立了一个自动化、规范化的新抗原预测平台。在该方法中，有机整合了从原始数据比对到潜在新抗原排序的多个步骤，可以通过简单的一个计算机命令，分析得到病人所有潜在的新抗原信息，并提供临床上可用于新抗原免疫治疗的多肽序列信息；虽然都是基于现有开源软件的方法进行数据的处理和分析，但相对于软件包相对独立，需要各自运行后整合而且软件之间存在复杂的依赖关系，本发明的规范化数据分析流程具有更便捷的操作性及可重复性。

目前已发布的亲和力预测数据库更新速度较慢，并且对于某些特定的HLA基因型没有相关的数据收录，我们根据临床上不断累积的数据对亲和力预测数据库进行扩充，可以取得更准确的预测结果，这也是本发明方法的一个重要突破。

本发明方法通过对多个高通量测序数据分析软件及抗原亲和力分析软件进行整合，得到一个标准化的分析流程。通过多种条件的运行测试，该预测流程具有较好的代码稳健性。

而利用本发明挑选出来的排名靠前的多肽，在临床应用后通过ELISPOT检测其特异性，结果与本发明新抗原分析流程的预测结果具有一致性，预测准确率达到50％以上(已检测的两个患者阳性率为10/13和8/12，其中肿瘤特异性强抗原阳性率为6/13和5/12)，显著高于其他相关报道，我们还可在临床反馈的基础上进一步改进和完善算法进一步提高预测的特异性和准确性。

具体实施方式

下面结合实施例，详细描述本发明的技术方案。

本发明的临床个体化肿瘤新抗原测序的预测方法，包括以下步骤：

步骤2、对DNA比对结果及RNA比对结果进行预处理；

步骤1中，所述将患者肿瘤细胞的测序数据与预设的对应参考基因组进行比对中，采用Burrows-Wheeler转化算法得到DNA比对结果；采用STAR算法得到RNA比对结果。

步骤1之前，还可包括以下步骤：

步骤0中，所述预处理包括修剪测序接头序列、修剪测序标签序列及去除测序质量较差的序列，其可采用Trimmomatic软件进行。

步骤2中，所述预处理可包括去重、添加组信息及碱基质量重计算等，其可采用GATK软件进行。

步骤3可包括以下步骤：

步骤3E、计算肿瘤的拷贝数变异文件；

步骤3G、根据RNA比对结果分析肿瘤细胞中的基因表达量。

步骤3D中，所述体细胞变异可为：针对注释后的实际体细胞突变位点进行进一步过滤，仅保留外显子区域中非同义体细胞突变位点。

步骤3A优选采用Mutect2软件进行分析；步骤3B优选采用GATK的HaplotyperCaller软件分析；步骤3C优选采用FilterMutectCalls软件进行分析；步骤3D优选采用ANNOVAR软件进行注释；步骤3E优选采用cnvkit软件进行计算；步骤3F中，所述计算肿瘤纯度优选采用ABSOULTE软件进行计算，所述克隆结构优选采用pyClone软件进行计算；所述克隆结构一般是指突变基因克隆比，即每个体细胞变异位点的克隆比例。

步骤3G中，基因表达量的分析可以为采用featureCounts软件进行分析得到counts数值，并将所得到的counts数值转换成每百万测序片段中转录本测序片段数(TPM数值)。

所述将所得到的counts数值转换成每百万测序片段中转录本测序片段数中，其转换公式可为:

其中，R为基因的count数目，L为基因的长度，

步骤4可包括以下步骤：

步骤4A、预测患者肿瘤细胞的MHC I类分子的基因型；

步骤4A采用polsover软件进行分析；步骤4C采用IEDB提供的MHC class Iepitopepredictors软件进行分析；步骤4D采用netCTLpan软件进行分析计算。

步骤4B具体可以为：下载UCSC网站提供的人基因组蛋白序列fasta文件，根据注释后的实际体细胞突变位点提取出变异氨基酸前后各15个氨基酸的序列，分别获得野生型蛋白质和突变型蛋白质的侧翼氨基酸序列，将所有的序列组合成多肽序列的fasta文件，即非同义突变蛋白数据。

步骤5中，所述根据分析结果、潜在新抗原的亲和性及多肽递呈效力对各对应的新抗原进行评分的方法优选为：

score＝abundance*dissimilarity*clonality

其中，

abundance＝L(IC_m)*Freq*tanh(TPM)；

dissimilarity＝(1–L(IC_w/50000))/2；

clonality＝NCTL*MC；

L(x)＝1/(1+e^x),tanh()为双曲正切函数计算，score是新抗原计算获得的分数，abundance是突变多肽的丰度,TPM为基因表达丰度，Dissimilarity是不相似度(dissimilarity)衡量突变肽段与对应正常肽段亲和力的差异，克隆性(clonality)衡量新抗原由突变至短肽的效率NCTL及新抗原在所有肿瘤细胞中分布的比例MC，Freq为体细胞变异等位基因表达丰度比例所述体细胞变异等位基因表达丰度比例是指体细胞变异等位基因表达丰度与肿瘤细胞中的基因表达量的比值，IC_m为亲和力预测时得到的突变多肽亲和力，IC_w为亲和力预测时得到的野生多肽结合亲和力,NCTL为递呈效力，MC为突变基因克隆比。

还可包括以下步骤：

实施例

本实施例中，具体的临床个体化肿瘤新抗原测序的预测方法如下：

具体为：获取患者原始的肿瘤细胞DNA/RNA以及患者血液的DNA材料的高通量测序文件(fastq文件)后，使用Trimmomatic软件对测序fastq文件进行修剪测序接头序列、修剪测序标签序列、去除测序质量较差的。

其中对于参数的选择为：1)选择PE模式也即双端测序(pair-end)模式进行数据分析；2)ILLUMINACLIP参数设置为TruSeq3-PE.fa:2:30:10，即使用ILLUMINA公司的TruSeq3-PE建库方案中的接头序列信息进行接头的去除，该信息存在在TruSeq3-PE.fa文件中，且目标序列与源序列不能超过2个碱基错配，回文序列比对分阈值设置为30(如果文库***片段比测序读长短，则正反向测序片段中一段碱基可以完全反向互补，将两个接头序列与测序片段进行比对，同时两条配对的测序片段之间也互相比对，可以将测序片段末端的接头序列准确去除，该参数即控制两条配对的测序片段之间回文结构的比对质量值超过30才用此种方法进行接头的切除)，最后设置切除接头序列的最低比对分为10；3)LEADING参数设定为3，即将测序片段最前端测序质量值低于3的连续若干个碱基切除掉(测序片段两端测序质量相对较差)；4)TRAILING参数设定为3，即将测序片段最后端测序质量值低于3的连续若干个碱基切除掉，理由同上；5)SLIDINGWINDOW参数设置为4:15，即软件滑动窗口设置为4碱基的大小，如果这4个碱基的质量平均值低于15则将该窗口中的碱基切除；6)MINLEN参数设置为36，即修剪后的测序片段长度低于36碱基数，则将该测序片段丢弃不用。

Trimmomatic程序将原始的fastq文件作为输入文件，运行得到处理后的fastq文件，称为clean fastq文件，用于后继步骤1的测序片段的比对操作。

步骤1、将患者肿瘤细胞的测序数据与预设的对应参考基因组进行比对，得到DNA比对结果及RNA比对结果。

具体为：步骤0已得到纯化过的fastq文件(即clean fastq文件)，接下来将这些测序片段比对到对应的参考基因组上。考虑到全外显子测序和全转录组建库方式的差异，采用BWA软件进行全外显子数据的基因组比对，采用STAR软件进行全转录组数据的基因组比对：

使用BWA的重要参数为：1)k参数设置为19，即种子序列的最短长度为19；2)w参数设置为100，如果比对缺口长度大于100则忽略掉；3)；4)r参数设置为1.5，如果种子长度超过k*r，也即19*1.5＝28.5bp，则在该种子序列内部寻找新的种子序列。该参数是启发式算法调节算法性能的关键参数，如果该值设置过大可以导致更快的比对速度，但是会牺牲部分的准确性；5)c参数设置为500，如果种子序列的重复次数超过500则忽略掉该种子序列；6)chimSegmentMin设置为0，只输出正确比对的测序片段，不输出嵌合片段；

使用STAR的参数为：1)将参数genomeDir设置为到存储参考基因组文件的目录路径；2)将参数readFilesIn设置为输入片段文件的路径；3)将参数twopassMode设置为基本(Basic)模式，通道映射模式，即采用双通道比对模式；4)将参数outReadUnmapp None设置为None，即不输出比对失败或部分比对成功的片段；5)将参数chimSegmentMin设置为12，即令嵌片段的最小长度为12bp；6)将参数chimJunctionOverhangMin设置为12，即嵌合结合处的最低允许突出长度为12bp；7)将参数alignSJDBoverhangMin设置为10，即在有注释的结合区域，最低运行突出10bp的单链序列；8)将参数alignMatesGapMax设置为200000，即两个配合之间的最大间隙设为200000bp；9)将参数alignIntronMax设置为200000，即最大内含子大小设置为200000bp；10)将参数chimSegmentReadGapMax设置为3，即嵌合片段之间的最大间隙为3bp；11)将参数alignSJstitchMismatchNmax设置为5-1 5 5，即能用于拼接缝合片段的最长不匹配数限制为：(1)非经典模式不超过5bp，(2)GT/AG和CT/AC模式不限制，(3)GC/AG和CT/GC模式不超过5bp,(4)AT/AC和GT/AT模式不超过5bp；12)将参数runThreadN设置成需要的并行线程数；13)将参数outFilterMultimapNmax设置为2，即允许单个片段最多能匹配到两个区域；14)将参数outFilterMismatchNoverLmax设置为0.04，即只有错配长度与片段长度比值不大于0.04时才输出该片段比对结果；15)将参数outSAMtype设置为BAMSortedByCoordinate，即设定输出为按照坐标排序的BAM文件；16)将参数outfileNameprefix设置为输出文件名的前缀，程序能自动根据该前缀生成相应的多个结果文件；17)将参数readFilesCommand设置为处理输入文件格式的shell命令，使其输出结果为FASTA或FASTQ文件。

经过该步骤，得到了已经比对到选定参考基因组上的测序序列文件，即BAM文件，其即为DNA比对结果及RNA比对结果。

步骤2、对DNA比对结果及RNA比对结果进行预处理。

具体为：在使用BWA或者STAR将测序片段比对到参考基因组后，需要对比对结果(包括DNA比对结果及RNA比对结果)进行去重、添加组信息、碱基质量重计算等工作。这些步骤可使用GATK进行操作。

使用GATK进行测序比对片段的去重时的参数为：1)将参数PG设置为MarkDuplicates，以启用PG记录创建的程序记录ID。此字符串可能附加后缀以避免与其他程序记录ID冲突；2)将参数PG_NAME设置为MarkDuplicates，即创建PG记录的PN标记的值；3)将MAX_FILE_HANDLES设为8000，即将片段结尾溢出到磁盘时保持打开的最大文件句柄数设为8000；4)将参数SORTING_COLLECTION_SIZE_RATIO设置为0.25，即确定某些排序集合使用的内存占用量为内存总容量的0.25倍；5)将参数OPTICAL_DUPLICATE_PIXEL_DISTANCE设置为100，即将两个重复簇之间的最大偏移量设置为100。

使用GATK进行测序比对片段组信息的添加时的参数为：1)参数ID设置为1，即读取组ID默认值设为1；2)参数version设置为false即不显示软件的版本号。

使用GATK进行测序片段中碱基的质量重计算时的参数为：1)将参数R设置为hg19参考基因组fasta文件；2)将known-sites参数设置为1000G_phase1.indels.hg19.sites.vcf.gz，Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz和dbsnp150，将这些已知的变异作为参考，重新评估BAM文件中碱基的质量。

经过该步骤，得到了质控后的BAM文件，该文件用于后继的步骤3A、步骤3E、步骤3G、步骤4A。

步骤3、利用DNA比对结果及RNA比对结果分析肿瘤细胞中的体细胞变异、克隆类型、肿瘤纯度、肿瘤细胞中的基因表达量及体细胞变异等位基因表达丰度，得到分析结果。

本步骤中，包括以下具体步骤：

步骤3A、利用DNA对比结果分析肿瘤细胞中的体细胞突变位点。

经过步骤2对BAM文件处理，得到了肿瘤样本与正常对照样本在基因组上的比对结果BAM文件，接下来通过mutect2软件对比分析正常--肿瘤成对样本的BAM文件，找出所有潜在的突变位点。

使用mutect2的参数为：1)将参数base-quality-score-threshold设置为18，即将碱基质量分数阈值设为18，低于此阈值的碱基将降至最低值6；2)将参数heterozygosity设置为0.001，即设任何基因座杂合度的先验可能性值为0.001；3)将indel-heterozygosity设置为1.25E-4，即设任何基因座杂合度的先验可能性值为1.25E-4；4)将参数mbq设置为10，即在检索突变时，允许的合格碱基的最低质量为10；5)将参native-pair-hmm-threads设为4，即本机pairHMM实现使用的线程数为4；6)将参数minimum-mapping-quality设为4，即比对质量低于20的片段都舍去；7)将参min-read-length设为30，即仅保留长度大于等于30都片段。

经过该步骤，得到了包含体细胞变异位点信息的vcf文件，该文件用于后继的步骤3B、步骤3C。

步骤3B、根据体细胞突变位点结合RNA比对结果分析得到体细胞变异等位基因表达丰度及肿瘤细胞中的基因表达量。

具体为：使用vcf2bed软件将存储体细胞变异位点信息的vcf文件转换成位点坐标的bed文件，使用GATK的HaplotyperCaller分析RNA测序比对文件中这些位点的信息。

使用HaplotyperCaller的参数为：1)-L提供vcf2bed生成的bed文件；2)--reference使用GRC37/hg19参考基因组；3)-ERC选择GVCF格式，输出gvcf文件；4)--output-mode选择EMIT_ALL_SITES，无论是否判别为变异，将目标区域中所有的位点输出。

接下来，使用GATK的GenotypeGVCFs命令对gvcf进行基因型的鉴定，并获得体细胞突变位点在RNA水平上野生型碱基数量和突变型碱基数量。

经过该步骤，得到了体细胞变异等位基因表达丰度及肿瘤细胞中的基因表达量，用于步骤3C和步骤5。

步骤3C、通过体细胞变异等位基因表达丰度过滤体细胞突变位点中的假阳性位点，得到实际体细胞突变位点。

具体为：实施完步骤3A后我们可以找出肿瘤细胞区别于正常细胞的体细胞突变位点，但是其中还掺杂了许多由于测序技术等原因导致的假阳性位点，接下来在本步骤中，使用FilterMutectCalls软件过滤其中的假阳性位点，结合体细胞变异等位基因表达丰度进一步筛选出有RNA表达的变异位点。

使用FilterMutectCalls软件的参数为：1)将参数base-quality-score-threshold设置为18，即低于阈值18的碱基质量将降至最低值6；2)将参数heterozygosity设置为0.001，即设任何基因座杂合度的先验可能性值为0.001；3)将参数heterozygosity-STDEV设置为0.01，即SNP和indel的杂合性标准偏差值设为0.01；4)将参数indel-heterozygosity设置为1.25E-4，即设任何基因座杂合度的先验可能性值为1.25E-4；5)将参数max-alt-allele-count设置为1，即删除具有多个等位基因的突变位点；6)将参数max-germline-posterior设置为0.1，即等位基因是种系变异的最大后验概率值为0.1；7)将参数max-strand-artifact-probability设置为0.99，如果链伪影的概率超过0.99，则过滤该变异位点；8)将参数mbq设置为10，即过滤掉质量10以下的碱基；9)将参数min-median-base-quality设置为20，即过滤掉等位片段碱基质量中位值小于20的突变位点；10)将参数min-median-mapping-quality设置为30，即过滤掉等位片段比对质量中位值小于30的突变位点；11)将参数min-median-read-position设置为5，即过滤掉靠近片段末端5个碱基以内的突变位点；12)将参数native-pair-hmm-threads设置为4，即本机pairHMM实现使用的线程数为4；13)将参数normal-p-value-threshold设置为0.001，即将正常伪影过滤器的P值阈值设为0.001；14)将参数output-mode设置为EMIT_VARIANTS_ONLY，即只输出突变位点；15)将参数stand-call-conf设为10.0，即将搜索变异位点的最小phred-scaled置信度阈值设为10.0。

经过本步骤，得到了过滤后的体细胞突变vcf文件，用于后继的步骤3D。

步骤3D、对实际体细胞突变位点进行注释及进一步过滤，得到体细胞变异。

具体为：使用ANNOVAR软件注释实际体细胞突变位点的功能，ANNOVAR软件的使用参数为：1)将参数buildver设为hg19，即使用hg19的基因构建的数据库作为参考；2)将参数maxgenethread设置为6。指定基因注释中允许的最大线程数为6；3)使用refGene,avsnp150,clinvar,cosmic,1000gEAS,EXAC,gnomad,icgc,dbnsfp33作为参考进行注释。根据注释结果进一步过滤突变位点，只保留外显子区域的非同义突变位点。

经过本步骤，得到了包含外显子区域的非同义突变位点注释信息的vcf文件，即体细胞变异，用于后继的步骤3F及步骤4B。

步骤3E、计算肿瘤的拷贝数变异文件。

具体为：为了估计肿瘤的纯度，需要准备肿瘤的拷贝数变异文件(seg文件)。为了计算肿瘤的体细胞拷贝数变异，使用cnvkit软件对肿瘤和癌旁的bam文件进行分析，分析的具体参数为：1)将参数method设置为hybrid，即测序方案设为杂交捕获；2)将参数process设置为process each BAM in serial，即并行处理每个BAM。

经过本步骤，得到了包含肿瘤拷贝数变异的seg文件，用于后继的步骤3F。

步骤3F、根据体细胞变异及肿瘤的拷贝数变异文件计算肿瘤纯度及克隆结构。

本步骤具体可包括以下3个部分：

I、得到MAF文件

为了估计肿瘤的纯度，不仅需要准备肿瘤的拷贝数变异文件(seg文件)，还需要包含外显子区域非同义突变位点的MAF文件，在获得包含外显子区域的非同义突变位点注释信息的vcf文件(即步骤3D获得的vcf文件)后，利用该vcf文件用于肿瘤纯度的估计，首先需要使用vcf2maf软件将vcf文件转换成MAF文件，vcf2maf软件的具体使用参数为：1)将参数vep-forks设置为4，即运行VEP时使用的并行进程数为4；2)将参数ncbi-build设置为GRCh37，即使用NCBI的GRCh37为参考，组装突变位点的MAF。

经过该部分，得到了包含外显子区域非同义突变位点的MAF文件，用于后继的II、III部分。

II、计算肿瘤纯度

在上述I部分及步骤3E中分别得到了包含外显子区域非同义突变位点的MAF文件和肿瘤的拷贝数变异文件(seg文件)，在该部分中，使用这两个文件作为输入文件，利用ABSOLUTE软件进行肿瘤纯度的计算。ABSOLUTE软件的相关参数为：1)将参数sigma p设为0，即用于模式搜索的超出样本水平方差的临时值设置为0；2)将参数max sigma h设置为0.015。超出样本水平方差的最大值设为0.015；3)将参数min ploidy设置为0.95N，指定算法要考虑的最小倍性值N为0.95N；4)将参数max ploidy设为10N，即指定要考虑的最大倍性值N为10N，并丢弃可能更大倍性值的模型；5)将参数primary disease设置为对应的样本名。输入'NA'以使用泛癌核型为参考。如果提供的输入与列表项中类型不匹配，则ABSOLUTE默认为泛癌为参考；6)将参数max as seg count设为1500，即将等位基因区段数大于1500的样本标记为“失败”；7)将参数max sigma h设为0.005。有时，由于数据中的噪声，ABSOLUTE建模可以将归因于肿瘤亚克隆的基因组部分小于零。此参数指定基因组的最大允许组分为0.005，可以在不放弃给定解决方案的情况下建模为小于零；8)将参数max nonclonal设为0.05，即可以被建模为非克隆的最大基因组分值设为0.05。对于预期具有较高比例的异质的样品，可以适当增加该参数。

经过该部分，得到了肿瘤组织的纯度值，用于后继的III部分。

III、分析肿瘤克隆结构

在计算肿瘤纯度的同时，还可以根据包含外显子区域非同义突变位点的MAF文件和肿瘤的拷贝数变异文件(seg文件)分析肿瘤的克隆状态。利用包含外显子区域非同义突变位点的MAF文件、肿瘤的拷贝数变异文件(seg文件)及肿瘤组织的纯度值作为输入文件，采用pyClone软件进行克隆数的分析。pyClone软件使用的具体参数为：1)使用run_analysis_pipeline子命令；2)in_files参数使用根据maf整理得到的文件；3)prior参数使用total_copy_number；4)tumour_contents使用ABSOLUTE计算得到的肿瘤纯度值。

经过该部分，得到所有突变位点的克隆比例文件，即肿瘤克隆结构，用于后继的步骤5。

步骤3G、根据RNA比对结果分析肿瘤细胞中的基因表达量。

具体为：使用featureCounts软件分析RNA-seq比对结果(包含在步骤2中得到的质控后的BAM文件中)，得到基因表达量值FPKM(即counts数值)，并进一步对结果进行转换为TPM值。

使用featureCounts软件的参数为：1)将参数minOverlap设为1。片段分配使用片段的最小重叠碱基数为1；2)将参数T设为8，即使用8个线程数并行运行。

其转换公式为:

其中，R为基因的count数目，L为基因的长度，

经过该步骤得到包含所有基因表达量的文件，即基因表达量，用于步骤5。

步骤4、预测患者肿瘤细胞的MHC I类分子的基因型，并根据DNA分析结果及患者肿瘤细胞的MHC I类分子的基因型分析潜在新抗原的亲和性及计算多肽递呈效力。

本步骤包括以下具体步骤：

步骤4A、预测患者肿瘤细胞的MHC I类分子的基因型。

具体为：使用polysover软件对HLA等位基因进行推测，其相关参数如下:1)bam参数设置为步骤2得到的BAM文件；2)将参数race设置为Asia，即设置病人的种族为亚洲人(当然根据情况需要，也可设置为欧洲人或非洲人等)；3)将参数includeFreq设置为1，即群体级别等位基因频率用作先验条件；4)将参数build设置为hg19，即设置测序数据比对基因组为GRCh37/hg19；5)将参数format设置为STDFQ，即测序数据采用的是Hiseq，返回的测序数据是标准fastq文件；6)将insertCalc设置为1，即在模型中使用经验***大小分布模型。

经过该步骤，得到了肿瘤中HLA的分型信息，即患者肿瘤细胞的MHC I类分子的基因型，用于后继的步骤4C。

步骤4B、根据注释后的实际体细胞突变位点获取非同义突变蛋白数据。

具体为：将步骤3D中得到的体细胞变异进行筛选，保留其中的非同义突变。同时，下载UCSC网站提供的人基因组蛋白序列fasta文件，根据体细胞变异的注释文件中的变异位置提取出变异氨基酸前后各15个氨基酸的序列，分别获得野生型和突变型的侧翼氨基酸序列。将所有的序列组合成多肽序列的fasta文件。

经过该步骤，得到了包含体细胞突变位点多肽信息的蛋白质fasta文件，即非同义突变蛋白数据，用于后继的步骤4C和步骤4E。

步骤4C、根据非同义突变蛋白数据预测MHC I类分子结合野生型多肽和变异型多肽的能力，得到亲和力预测。

具体为：使用IEDB官网上提供的mhc_i软件对MHC I类分子结合多肽的能力进行预测，其相关参数如下：1)method参数设置为IEDB_recommended，即使用IEDB推荐的设置(其中含有ann，netmhcpan，pickpocket等方法)；2)mhc参数设置为polysolver得到的HLA分型结果；3)peptide_length长度设置为潜在的新抗原多肽长度，我们通常设置为8，9，10，11，12,13,14；4)input_file参数为步骤13整理得到的fastq文件。

经过该步骤，得到了所有潜在多肽与每种HLA分子的亲和力，即亲和力预测，用于后继的步骤4E。

步骤4D、分析多肽序列，得到多肽解离效率、转运效力及MHC结合效力，并根据所得到的多肽解离效率、转运效力及MHC结合效力计算多肽递呈效力。

具体为：使用netCTLpan软件分析多肽序列，得到多肽解离效率(Proteasomecleavage prediction)，转运效力(transport efficiency)，MHC结合效力，并根据得到结果计算递呈效力。使用netCTLpan软件的参数为:1)-a HLA等位基因,设置与多肽序列结合的HLA等位其中含有ann，netmhcpan，pickpocket等方法基因，可设置多个。2)-s按招某列排序输出，设置为0即按照综合结果输出(递呈效力)。3)-l肽段长度,设置需要计算太短的长度。

经过本步骤，得到了每个多肽被HLA分子递呈的效力值，即多肽递呈效力，用于后继的步骤5。

具体为：根据步骤4C得到的亲和力预测和步骤4B得到的非同义突变蛋白数据，使用脚本对结果文件进行处理。主要准则如下：1)变异氨基酸位于预测的多肽中；2)如果野生型序列和非同义突变序列不同时存在于预测结果则去除；3)输出野生型和非同义突变的相关信息，包括HLA类型，在原始的起始终止位置，长度，预测多肽的序列，预测结合多肽亲和力值(IC50)并将野生型和非同义突变信息整理至一行。

经过本步骤，得到了野生型/突变型有效多肽(包含由于体细胞突变被改变的那个氨基酸)和HLA结合的亲和力，即潜在新抗原的亲和性，用于后继的步骤5。

具体为：将结合多肽亲和力预测结果、突变频率、肿瘤纯度、多肽表达量等进行整合，并按照相关规则对预测结果进行综合打分，主要指标为在野生型和突变型多肽序列的亲和力值(IC50)、多肽表达水平、多肽长度，并对潜在新抗原进行排序。评分规则如下：

score＝abundance*dissimilarity*clonality

其中，

abundance＝L(IC_m)*Freq*tanh(TPM)；

dissimilarity＝(1–L(IC_w/50000))/2；

clonality＝NCTL*MC；

本例中，还具有如下步骤：

步骤6、基于实验结果，利用机器学习方法优化步骤5的潜在新抗原的亲和性的分析方法。

可具体包括以下三个具体步骤：

步骤6A、数据收集

从IEDB收集(http://www.iedb.org/database_export_v3.php)最新的MHC结合肽酶的综合数据集，并添加自己实验的MHC结合肽数据。

步骤6B、数据整理

整理已收集的数据，结合肽段的亲和力值(IC50)大于50000nM的测量值取为50000nM。将得到的结合肽段的亲和力值进行转换，转换公式为t＝(1-log(aff)/log(50000))。如若有重复序列的数据存在，取其结合肽段亲和力的均值作为最后的数据合并结果，并将之添加到训练数据。

步骤6C、数据集训练

将步骤6B的收集整理的结果分别使用不同的方法进行训练，包括pickpocket、netMHCpan、netMHCcons和/或netMHC等。用新训练获得的模型文件更新mhc_i软件原有的数据库文件，用于提高预测的准确性。

Claims

1.临床个体化肿瘤新抗原的预测方法，包括以下步骤：

步骤2、对DNA比对结果及RNA比对结果进行预处理；

步骤3、利用DNA比对结果及RNA比对结果分析肿瘤细胞中的体细胞变异、克隆类型、肿瘤纯度、肿瘤细胞中的基因表达量及体细胞变异等位基因表达丰度，得到分析结果；步骤3包括以下步骤：

步骤3E、计算肿瘤的拷贝数变异文件；

步骤3G、根据RNA比对结果分析肿瘤细胞中的基因表达量；

步骤4、预测患者肿瘤细胞的MHC I类分子的基因型，并根据DNA分析结果及患者肿瘤细胞的MHC I类分子的基因型分析潜在新抗原的亲和性及计算多肽递呈效力；步骤4包括以下步骤：

步骤4A、预测患者肿瘤细胞的MHC I类分子的基因型；

步骤4E、根据亲和力预测及非同义突变蛋白数据得到潜在新抗原的亲和性；

步骤5、根据分析结果、潜在新抗原的亲和性及多肽递呈效力对各对应的新抗原进行评分及排序后呈现；

所述的根据分析结果、潜在新抗原的亲和性及多肽递呈效力对各对应的新抗原进行评分的方法为：

score＝abundance*dissimilarity*clonality

其中，

abundance＝L(IC_m)*Freq*tanh(TPM)；

dissimilarity＝(1–L(IC_w/50000))/2；

clonality＝NCTL*MC；

L(x)＝1/(1+e^x),tanh()为双曲正切函数计算，score是新抗原计算获得的分数，abundance是突变多肽的丰度,TPM为基因表达丰度，Dissimilarity是不相似度衡量突变肽段与对应正常肽段亲和力的差异，克隆性衡量新抗原由突变至短肽的效率NCTL及新抗原在所有肿瘤细胞中分布的比例MC，Freq为体细胞变异等位基因表达丰度比例，所述体细胞变异等位基因表达丰度比例是指体细胞变异等位基因表达丰度与肿瘤细胞中的基因表达量的比值，IC_m为亲和力预测时得到的突变多肽亲和力，IC_w为亲和力预测时得到的野生多肽结合亲和力,NCTL为递呈效力，MC为突变基因克隆比

2.如权利要求1所述的临床个体化肿瘤新抗原的预测方法，其特征在于，步骤1中，所述将患者肿瘤细胞的测序数据与预设的对应参考基因组进行比对中，采用Burrows-Wheeler转化算法得到DNA比对结果；采用STAR算法得到RNA比对结果。

3.如权利要求1所述的临床个体化肿瘤新抗原的预测方法，其特征在于，步骤1之前，还包括以下步骤：

4.如权利要求3所述的临床个体化肿瘤新抗原的预测方法，其特征在于，步骤0中，所述预处理包括修剪测序接头序列、修剪测序标签序列及去除测序质量较差的序列，采用Trimmomatic软件进行。

5.如权利要求1所述的临床个体化肿瘤新抗原的预测方法，其特征在于，步骤2中，所述预处理包括去重、添加组信息及碱基质量重计算，采用GATK软件进行。

6.如权利要求1所述的临床个体化肿瘤新抗原的预测方法，其特征在于，步骤3D中，所述体细胞变异是指：针对注释后的实际体细胞突变位点进行进一步过滤，仅保留外显子区域中非同义体细胞突变位点。

7.如权利要求1所述的临床个体化肿瘤新抗原的预测方法，其特征在于，步骤3A采用Mutect2软件进行分析；步骤3B采用GATK的HaplotyperCaller软件分析；步骤3C采用FilterMutectCalls软件进行分析；步骤3D采用ANNOVAR软件进行注释；步骤3E采用cnvkit软件进行计算；步骤3F中，所述计算肿瘤纯度采用ABSOULTE软件进行计算，所述克隆结构采用pyClone软件进行计算；所述克隆结构是指突变基因克隆比，即每个体细胞变异位点的克隆比例。

8.如权利要求1所述的临床个体化肿瘤新抗原的预测方法，其特征在于，步骤3G中，基因表达丰度的分析是指采用featureCounts软件进行分析得到counts数值，并将所得到的counts数值转换成每百万测序片段中转录本测序片段数。

9.如权利要求8所述的临床个体化肿瘤新抗原的预测方法，其特征在于，所述将所得到的counts数值转换成每百万测序片段中转录本测序片段数中，其转换公式为:

其中，R为基因的count数目，L为基因的长度，

10.如权利要求1所述的临床个体化肿瘤新抗原的预测方法，其特征在于，步骤4A采用polsover软件进行分析；步骤4C采用IEDB提供的MHC class I epitope predictors软件进行分析；步骤4D采用netCTLpan软件进行分析计算。

11.如权利要求1所述的临床个体化肿瘤新抗原的预测方法，其特征在于，步骤4B具体为：下载UCSC网站提供的人基因组蛋白序列fasta文件，根据注释后的实际体细胞突变位点提取出变异氨基酸前后各15个氨基酸的序列，分别获得野生型蛋白质和突变型蛋白质的侧翼氨基酸序列，将所有的序列组合成多肽序列的fasta文件，即非同义突变蛋白数据。