本发明涉及转录组分析领域，更特别地，涉及一种基于转录组测序数据检测融合基因的方法。
背景技术：
基因重排是生物体中遗传物质之间时有发生的现象，由于基因重排常常导致原本不在一个顺反子下的一个或多个基因或基因片段形成融合基因，并作为一个顺反子转录，这将导致某些基因的活化、失活或产生新功能。许多疾病的发生都伴随有融合基因现象，例如，白血病常伴随有bcr/abl、AML1/ETO、CBFβ/MYH11、PML/RARα等融合基因，多种实体瘤中也发现融合基因，非小细胞肺癌中有EML4-ALK，前列腺癌中有SLC45A3-ELK4，横纹肌肉瘤中有PAX3-FOXO1等等。科学研究发现，一些融合基因参与了相关疾病的致病过程，因此，这些融合基因的检测可作为诊断标准之一，甚至可将其作为治疗靶点。
目前，对融合基因的检测一般基于对二代转录组测序数据进行分析来得到。其主要通过比对两种类型的读序(reads)来进行。一种为非一致成对读序(discordant paired-end reads，即，成对的reads分别比对到参与融合的5’伴侣基因和3’伴侣基因)，另一种为结合读序(junction reads，即，reads的比对跨越了融合位点)。通过鉴定这两类读序的支持情况，可检测出转录产物中的融合基因。目前基于这类方法已经开发出来的软件很多，包括SOAPfusion、Defuse、SOAPfuse、FusionCatcher、FusionMap、Tophat-fusion、ChimeraScan、Star-fusion等。但是由于转录组的复杂性和二代RNA-seq测序读序读长的限制，短读序比对本身就面临非常大的挑战。参考基因组上的重复序列以及参考基因组本身的不完整，很容易导致对读序比对位置和比对唯一性做出错误判断。为了控制比对错误导致的假阳性，需要设置比较严格的过滤参数来对候选融合基因结果进行过滤，但这样做往往导致很多真阳性结果也被过滤掉。另一方面，二代测序在文库构建过程中会随机产生一些来自不同基因的片段随机连接在一起形成的嵌合读列，采用现有的融合基因检测方法，我们无法区分这些嵌合序列和真正的融合基因序列，这些随机产生的嵌合读列也会被当做融合基因检测出来，造成结果的假阳性。由于上述原因，单纯基于二代RNA-seq测序的融合基因检测方法很难在准确性和假阴性率两个方面做到平衡。
三代转录组测序又称全长转录组测序，其每个读序的长度都远远大于二代转录组测序读长。相比二代测序的短读长读序，三代测序的长读长读序在比对过程中可以更有效地降低比对错误造成的假阳性，但是也无法完全避免这种类型的假阳性。同时，三代转录组测序在文库构建过程中也会随机产生一些来自不同基因的片段随机连接形成的嵌合读序，导致假阳性融合。
因此，需要一种新的检测转录产物中融合基因的方法。
技术实现要素：
为解决以上问题，本发明提供了一种基于转录组测序数据检测融合基因的方法，其包括以下步骤：
S1：对样本进行二代转录组测序和三代转录组测序，分别得到二代转录组测序数据和三代转录组测序数据；
S2：将所述三代转录组测序数据与参考基因组进行比对，鉴定可能发生了基因融合的FLNC读序以及可能参与融合的基因对，提取所述可能发生了基因融合的FLNC读序的序列，并判断融合位置；
S3：将所述二代转录组测序数据比对到S2中得到的可能的融合基因FLNC读序，根据比对结果中非一致性成对读序以及结合读序的个数，鉴定确实发生了融合的基因对。
进一步地，S2包括以下步骤：
S2.1：将所述三代转录组测序数据与参考基因组进行比对，得到比对到参考基因组上两个不同位置的FLNC读序；
S2.2：判断所述参考基因组上两个不同位置以及所述FLNC读序上与所述两个不同位置对应的片段是否满足融合基因判定条件，当满足所有所述融合基因判定条件时，则将所述FLNC判定为可能发生了基因融合的FLNC读序并得到可能参与融合的基因对，提取所述可能发生了基因融合的FLNC读序的序列，并判断融合位置。
进一步地，所述融合基因判定条件为：
1)所述参考基因组上的两个不同位置分别对应所述FLNC读序的5’片段和3’片段；
2)所述5’片段和3’片段在所述FLNC读序上的位置满足不超过最大重叠长度和最大间隔长度，并且不小于最小总长度；
3)所述5’片段和3’片段在所述参考基因组上比对满足最低比对一致性；
4)所述参考基因组上的两个不同位置满足以下条件之一：a、处于不同染色体上；b、处于同一染色体上但方向相反；c、处于同一染色体上且方向相同，但是距离超过基因组注释的最大内含子长度；以及
5)所述参考基因组上的两个不同位置处有基因注释信息，并且根据所述基因注释信息可确定所述参考基因组上的两个不同位置的基因注释结构分别与其对应的所述5’片段和3’片段的基因结构一致。
进一步地，所述最大重叠长度和最大间隔长度为5-20bp，所述最小总长度为所述FLNC读序长度的10-20％，所述最低比对一致性为80-95％，所述最大内含子长度为50kb。
进一步地，S3包括以下步骤：
S3.1：将所述二代转录组测序数据与所述可能发生了基因融合的FLNC读序进行比对，针对每个所述可能发生了基因融合的FLNC读序鉴定非一致性成对读序和结合读序；
S3.2：当支持S2中鉴定的可能参与融合的基因对的所述可能发生了基因融合的FLNC读序的个数和所述结合读序的个数满足判定最小个数并且所述非一致性成对读序的对数满足最小判定对数时，判定所述可能参与融合的基因对发生了融合。
进一步地，所述可能发生了基因融合的FLNC读序的判定最小个数为1个。
进一步地，所述结合读序的判定最小个数为1个。
进一步地，所述非一致性读序的判定最小对数为1对。
三代全长转录组测序能够覆盖绝大多数的转录本序列，因此也能完整测序融合基因的融合转录本序列。二代测序数据测序深度高，对于低丰度的融合基因也能提供足够的reads支持。另外，二代测序和三代测序作为两个不同的技术手段，能够有效避免单一技术系统误差导致的假阳性或假阴性问题。例如针对两种测序技术文库构建过程中各自产生的嵌合读序，由于这些嵌合读序是随机产生的，我们通过要求在两种文库中都检测到融合的支持证据，即可有效避免嵌合读序造成的假阳性。FLNC本身作为全长转录本序列，也给二代RNA-seq数据比对提供了一个准确的候选融合转录本的参考序列，极大的提高了二代数据比对和检测融合基因的效率和准确性。
本发明通过结合三代转录组测序与二代转录组测序来检测融合基因，避免了单纯使用二代转录组测序因为读序长度短而带来的假阳性，以及单一测序技术因为嵌合读序导致的假阳性，同时也避免了原来二代数据做融合基因检测时融合转录本难以重构的问题，使得结合二代和三代测序支持证据的融合基因检测结果更为可靠。
具体实施方式
以下结合实例对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。
我们将本发明所述的方法的一个实施例应用到一个大豆转录组的项目当中。该项目对大豆转录组用三代ISO-seq测序方法和二代RNA-seq测序方法分别进行了测序。测序样品为大豆不同组织和发育时期的混合样品。其中三代转录组建了两个文库，文库大小分别为0.6-2.5kb和>1.5kb，两个文库用PacBio RSII测序平台各测序产生了16个cell和7个cell。二代测序建了一个文库，文库大小为200bp，共得到6G的RNA-seq数据。后续主要分析过程如下：
1)三代测序原始数据经过SMRT analysis软件中的RS_subreads、RS_ReadsofInsert和IRS_Isoseq三个pipeline对测序数据进行数据预处理和质控，并得到全长FLNC读序。二代数据经过FastQC软件进行预处理和质控，得到clean RNA-seq读序；
2)FLNC读序通过GMAP比对软件比对到参考基因组，得到比对结果文件；
3)比对结果文件进行筛选，找到分段比对到基因组上两个位置并满足以下条件的FLNC比对结果：(1)两个比对位置分别对应FLNC的5’和3’(2)两个片段间的重叠长度和间隔长度都小于10bp(3)每个片段的区域比对一致性都大于90％，且FLNC读序总比对长度大于90％(4)如果是同染色体同方向，则两个比对位置至少相距50k bp；
4)根据大豆基因组注释文件，挑选两个比对位置都有基因注释，且FLNC比对的exon/intron结构跟注释基因的结构相同的FLNC读序作为候选融合转录本，并记录参与融合的基因位置信息和融合位点信息；
5)提取支持候选融合基因的FLNC读序作为全长融合转录本序列，并把二代RNA-seq读序比对到融合转录本序列。根据比对结果检测支持融合基因的非一致性成对读序和结合读序；
6)筛选至少有一条FLNC读序支持，且至少有一对非一致性成对读序和一个结合读序支持的候选融合基因作为最终的融合基因检测结果。
通过上述分析，我们一共检测到了225个融合基因，其中染色体间的融合有209个，染色体内的融合16个。我们也通过JBrowser可视化软件对一些融合基因结果做了可视化作图，进一步确认了结果的可靠性。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。