一种过滤全长转录本的方法和系统与流程

1.本发明涉及生物信息技术领域，具体涉及一种过滤全长转录本的方法和系统。


背景技术：

2.传统的基因表达研究，是在个体或者组织水平上展开，即将个体或者组织样本所包含的数百万或更多个细胞作为一个整体。但是，明显的，同一样本中不同类型的细胞、相同类型不同时期的细胞，其基因表达模式是不同的，主要体现在具体表达什么基因、同一个基因怎么表达(转录本的差异)以及具体基因表达了多少(即表达量差异)几个方面。平均化的传统基因表达研究明显无法解决这个技术问题，因此单细胞技术孕育而生。
3.2009年第一个单细胞转录组测序技术被成功开发(tang f,barbacioru c,wang y,et al.mrna-seq whole-transcriptome analysis of a single cell[j].nature methods,2009,6(5):377-382.)，随后smart-seq技术(daniel,luo,shujun,wang,yu-chieh,et al.full-length mrna-seq from single-cell levels of rna and individual circulating tumor cells[j].nature biotechnology,2012,30(8):777-782.)、smart-seq2技术(serra,lorrayne,chang,dennis,macchietto,marissa,et al.adapting the smart-seq2 protocol for robust single worm rna-seq[j].bio protocol,2018,8(4).)的开发，10
×
genomics chromium(https://www.10xgenomics.com/)、illumina bio-rad single-cell sequencing solution(https://www.bio-rad.com/)等商业化平台的推出，使得单细胞研究迅速发展。
[0004]
现阶段，较成熟单细胞研究都是与第二代测序技术相结合。第二代测序技术(也称下一代测序技术)，以短读长、高准确度、高通量、低价格为显著特征。根据二代测序数据(短读长读段)可以在基因层面进行比对，即确定每一条短读长读段与具体基因的关系。但是，由于每个基因在转录时存在多种剪切模式，产生不同转录本，受限于短读长，无法直接确定具体每一条短读长读段与转录本的关系。因此，第二代测序技术只能利用生物信息学算法进行拼接以获得转录本信息，其数据在可变剪切模式、转录本分类上的适用性非常有限。
[0005]
第三代测序技术平台，以pacific biosciences和oxford nanopore technologies为代表，具有长读长的显著特点，为测序领域带来了重大的变革。作为新兴的应用方向，三代测序技术与单细胞技术的结合还有诸多技术问题。长读长测序比短读长测序技术具有更高的错误率和更低的测序深度。在数据分析时，如果使用传统的“一刀切”的处理方式，势必无法获得低深度但反应生物体真实转录情况的数据，而这些数据对于细胞水平的研究而言非常重要，它们很可能与单细胞的异质性相关。


技术实现要素：

[0006]
有鉴于此，本发明提供了一种过滤全长转录本的方法和系统，实现了对全长转录本数据的精确过滤，具备较高的准确性和灵敏性。
[0007]
为了实现上述目的，本发明的技术方案具体如下：
[0008]
一种过滤全长转录本的方法，包括以下步骤：
[0009]
s1、将测序获得的转录本读段与参考基因组进行比对，根据转录本剪切模式，进行聚类去冗余，得到第一转录本；所述第一转录本包括，与参考基因组剪切模式相同的转录本集合a，以及与参考基因组剪切模式不同的转录本集合b；
[0010]
s2、根据平均测序深度、外显子支持数以及转录本读段支持数，过滤转录本集合b，得到转录本集合c；
[0011]
s3、取转录本集合a与转录本集合c并集，与参考基因组重比对，再次聚类去冗余，得到第二转录本，即过滤后的全长转录本。
[0012]
在上述技术方案中，可以理解的是，所述第一转录本为转录本集合a和转录本集合b的并集；所述参考基因组中包括序列信息和结构注释信息；转录本读段与参考基因组比对后，能够获得各条读段的转录本剪切模式。
[0013]
优选的，步骤s1中所述测序为三代测序，所述转录本读段的样本类型为单细胞。
[0014]
优选的，步骤s1与步骤s3所述聚类去冗余的方法相同，具体为：根据3'端的外显子和剪切位点的完全匹配，将相同剪切模式下的转录本读段聚为一类，且一种剪切模式下只保留唯一最长的转录本。在上述技术方案中，由于一个基因模型下面有很多相似的转录本，故根据3'端的外显子和剪切位点的完全匹配对相似转录本进行聚类；而且由于rna易存在5'端的降解，所以仅5'端不同的转录本在本发明中被归于一类。
[0015]
优选的，步骤s2所述外显子支持数的计数方法为：根据集合b内不同转录本中具有重叠关系的外显子之间的相似性，判断是否为相同外显子，并进行计数。
[0016]
更加优选的，所述外显子之间的相似性大于等于70％时判定为相同外显子。
[0017]
所述外显子之间的相似性的计算方法优选为：设置窗口区域，并对其进行切分得到切分窗口，对切分窗口进行赋值并均一化，将各外显子所落入的切分窗口赋值之和计为该外显子得分，通过比较各外显子得分获得外显子之间的相似性。
[0018]
更加优选的，所述窗口区域的起点为存在重叠比对关系的外显子的最上游5’端，所述窗口区域的终点为存在重叠比对关系的外显子的最下游3’端；并取所述切分窗口赋值的对数值进行均一化。其中，所述窗口区域、切分窗口、外显子均按本领域惯用5’至3’方向表示。
[0019]
优选的，步骤s2所述转录本集合b过滤方法具体为：
[0020]
1)对于参考基因组中已知基因但剪切模式不同的转录本，转录本中的每一个外显子均满足以下条件，则保留：
[0021]
若外显子与参考基因组中的外显子有重叠，外显子支持数＜第一阈值且转录本读段支持数≥第二阈值，或外显子支持数≥第一阈值；
[0022]
若外显子与参考基因组中的外显子无重叠，外显子长度＜第三阈值且外显子支持数≥第四阈值，或外显子长度≥第三阈值且外显子支持数≥第五阈值；
[0023]
2)对于参考基因组中未知基因的转录本，保留平均测序深度≥第六阈值的转录本。
[0024]
更加优选的，所述第一阈值为2，所述第二阈值为2，所述第三阈值为200bp，所述第四阈值为10，所述第五阈值为5，所述第六阈值为1000。
[0025]
在上述技术方案中，所述转录本集合b中平均测序深度为转录本所含每个碱基测
序深度总和除以该转录本碱基总数，所述外显子支持数为外显子在转录本集合b所有的转录本中出现的总次数，所述转录本读段支持数为转录本集合b中每条转录本下的转录本读段数。
[0026]
优选的，所述方法还包括步骤s4，即对第二转录本基于剪切位点进行转录本分类。
[0027]
更加优选的，分类类型具体为：完全匹配的转录本，部分匹配的转录本，包含已知剪切位点重新组合的转录本，包含新剪切位点的转录本，反义转录本，融合转录本，内含子区的转录本，基因间区的转录本以及跨内含子和外显子区域的转录本。其中，所述完全匹配的转录本是指剪切位点与参考基因组中的已知注释信息完全一致的转录本；所述部分匹配的转录本是指剪切位点与参考基因组中的已知注释信息只有部分序列信息相同的转录本；所述包含已知剪切位点重新组合的转录本是指各剪切位点均是已知的，因剪切位点间形成了新的排列组合，生成了不同的转录本；所述包含新剪切位点的转录本是指出现新的剪切位点的转录本；所述反义转录本是指序列方向和参考基因组的序列方向相反的转录本；所述融合转录本是指跨越两个不同的基因的转录本；所述内含子区的转录本是指完全比对到了基因组的内含子区域的转录本；所述基因间区的转录本是指位于不同基因间的区域的转录本；所述跨内含子和外显子区域的转录本是指一半位于内含子区域一半位于外显子区域的转录本。
[0028]
本发明还提供了一种过滤全长转录本的系统，包括：
[0029]
比对模块：用于转录本读段与参考基因组的比对，获得各条读段的转录本剪切模式；
[0030]
聚类去冗余模块：根据转录本剪切模式，将相同剪切模式下的读段聚为一类且仅保留该类中唯一最长的转录本；
[0031]
转录本过滤模块：以平均测序深度、外显子支持数以及转录本读段支持数为参数，对新剪切模式的转录本进行过滤；所述新剪切模式的转录本为聚类去冗余后，转录本的剪切方式在参考基因组中没有报道，所述没有报道具体为：至少有一处的剪接方式与参考基因组中已知的剪切模式不同；
[0032]
过滤转录本输出模块：对过滤后的转录本进行输出。
[0033]
本发明的有益效果为：本发明提供的全长转录本过滤方法能够针对性的处理参考基因组中未知剪切模式的转录本，有效保留低丰度可信转录本；以平均测序深度、外显子支持数以及转录本读段支持数为参数，实现新剪切模式的转录本的过滤，且过滤结果相较于现有技术均具备较高的灵敏性和准确性；在实际应用中，可通过把控外显子之间的相似性，实现过滤结果的优化。
附图说明
[0034]
图1为实施例1中切分窗口的结构示意图(图中，双竖线标识省略部分)；
[0035]
图2为本发明提供的过滤方法与现有技术的技术效果对比图。
具体实施方式
[0036]
为了更好的理解本发明，下面结合实施例对本发明做进一步的详细说明。
[0037]
需要说明的是，本发明中的术语“a”、“b”、“第一”、“第二”等是用于区别不同的对
象，而不是用于描述特定的顺序；本发明中的“集合a”、“集合b”包含至少一个全长转录本。
[0038]
实施例1
[0039]
使用ont测序平台，对人系细胞系样品进行单细胞cdna测序，测序数据量如下表所示：
[0040][0041]
采用本发明的方法和系统，对上述测序数据进行过滤，具体步骤如下：
[0042]
(1)将上述测序数据与参考基因组进行比对，获得各条转录本测序读段的转录本剪切模式。
[0043]
在具体的实施过程中，比对软件可采用minimap2、gmap等，选择比对结果中得分最高的比对为最终比对结果，同时过滤掉包含嵌合体序列的比对。嵌合体序列是指，测序过程中，两段来源于基因组不同位置的dna分子由于外界因素(例如实验建库)连接在一起，形成的非自然存在的dna分子；在比对结果中体现为，一条读段比对至基因组的两个或多个位置。
[0044]
(2)根据转录本的剪切模式，对测序数据进行聚类去冗余，得到第一转录本。聚类去冗余的具体过程为：根据3'端的外显子和剪切位点的完全匹配，将相同剪切模式下的转录本读段聚为一类，且一种剪切模式下只保留唯一最长的转录本。
[0045]
根据转录本剪切模式是否已知，将第一转录本分为转录本集合a和转录本集合b。其中，转录本集合a中的转录本的剪切模式与参考基因组中已报道的剪切模式完全一样，转录本集合b中的转录本的剪切模式与参考基因组中已报道的剪切模式不同。
[0046]
(3)根据平均测序深度、外显子支持数以及转录本读段支持数，过滤转录本集合b，得到转录本集合c。
[0047]
首先对转录本集合b中存在的具有重叠关系的外显子进行判定，具体过程如下：
[0048]
①
设置窗口区域，其中，窗口区域的起点为存在重叠比对关系的外显子的最上游5’端，窗口区域的终点为存在重叠比对关系的外显子的最下游3’端；
[0049]
②
对窗口区域进行切分得到切分窗口；本实施例中，将窗口区域切分成200个切分窗口，如图1中所示；
[0050]
③
对切分窗口进行赋值并均一化，将各外显子所落入的切分窗口赋值之和计为该外显子得分；本实施例中，从2开始对切分窗口赋值，取所述切分窗口赋值的对数值进行均一化处理，以图1中存在重叠关系的a、b、c、d这4个转录本为例，各外显子得分依次为：
[0051]
score(a)a＝log23 log24 log25
…
log2201；
[0052]
score(b)＝log22 log23 log24
…
log2157；
[0053]
score(c)＝log23 log24 log25
…
log279；
[0054]
score(d)＝log22 log23 log24
…
log278；
[0055]
对上述各分值进行两两比较，将差值小于30％(即相似性大于等于70％)的记为同一个外显子。
[0056]
其次，记录每一个外显子在转录本集合b中出现的次数，记为该外显子的支持数。外显子支持数越大，表明转录本中的这个外显子可信度越高。
[0057]
转录本集合b过滤的过程具体为：
[0058]
①
对于参考基因组中已知基因但剪切模式不同的转录本，转录本中的每一个外显子均满足以下条件，则保留：
[0059]
若外显子与参考基因组中的外显子有重叠，外显子支持数＜2且转录本读段支持数≥2，或外显子支持数≥2；
[0060]
若外显子与参考基因组中的外显子无重叠，外显子长度＜200bp且外显子支持数≥10，或外显子长度≥200bp且外显子支持数≥5；
[0061]
②
对于参考基因组中未知基因的转录本，保留平均测序深度≥1000的转录本。
[0062]
过滤后得到的转录本集合记为转录本集合c。
[0063]
(4)将集合a与集合c合并，对其并集采用与步骤(1)、(2)同样的方法，得到第二转录本，即过滤后的全长转录本。
[0064]
需要说明的是，本实施例中的窗口区域的切分数目、切分窗口的赋值方式、均一化的具体方法、取对数值的底数、过滤条件中各指标的具体数值等均为具体示例，并不限制本发明的保护范围。
[0065]
实施例2
[0066]
对实施例1得到的第二转录本进行分类，具体采用基于剪切位点进行转录本分类的方法。本实施例具体实施过程可参考专利201910388054.2。
[0067]
通过该分类方法，可将实施例1中的第二转录本分为以下九种类型，具体为：
[0068]
完全匹配的转录本：剪切位点与参考基因组中的已知注释信息完全一致的转录本；
[0069]
部分匹配的转录本：剪切位点与参考基因组中的已知注释信息只有部分序列信息相同；
[0070]
包含已知剪切位点重新组合的转录本：各剪切位点均是已知的，因剪切位点间形成了新的排列组合，生成了不同的转录本；
[0071]
包含新剪切位点的转录本：出现新的剪切位点的转录本；
[0072]
反义转录本：序列方向和参考基因组的序列方向相反的转录本；
[0073]
融合转录本：跨越两个不同的基因的转录本；
[0074]
内含子区的转录本：完全比对到了基因组的内含子区域的转录本；
[0075]
基因间区的转录本：位于不同基因间的区域的转录本；
[0076]
跨内含子和外显子区域的转录本：一半位于内含子区域一半位于外显子区域的转录本。
[0077]
把分类结果中的“包含已知剪切位点重新组合形成的转录本”和“包含新的剪切位点的转录本”归于新转录本。
[0078]
对比例1
[0079]
采用与实施例1相同的数据，先使用软件flair(https://github.com/brookslabucsc/flair)获得过滤后的全长转录本数据，然后使用实施例2的方式进行转录本的分类。
[0080]
对比例2
[0081]
采用与实施例1相同的数据，先使用ont官方软件pinfish(https://github.com/
nanoporetech/pipeline-pinfish-analysis)获得过滤后的全长转录本数据，然后使用实施例2的方式进行转录本的分类。
[0082]
对比例3
[0083]
与实施例1不同的是，步骤(3)中进行外显子分值比较时，将差值小于40％(即相似性大于等于60％)的记为同一个外显子；其他过程同实施例1完全一致；然后，使用实施例2的方式进行转录本的分类。
[0084]
通过对实施例2、对比例1、对比例2和对比例3的分类结果进行准确性和灵敏性评估，以验证全长转录本数据的过滤效果。其中，灵敏性是指判定的已报道转录本占参考基因组报道转录本的比例，反应了过滤方法在发现转录本方面的能力；准确性是指判定的已报道转录本占判定的总转录本的比例，反应了过滤方法在正确判定转录本方面的能力。
[0085]
实施例2与对比例1、对比例2的分类结果对比如下表所示：
[0086][0087]
由上表可知，实施例2的准确性和灵敏性均优于对比例1，准确性略低于对比例2，但是灵敏性大大高于对比例2。截取某基因分析结果图进行直观展示，如图2所示(其中填充方块表示外显子、线条表示内含子、箭头表示转录方向)，图中1为人类参考基因组hg38中已报道的转录本，2为对比例2获取的转录本，只获取了一条转录本，说明灵敏性较差；3是对比例1获取的转录本，虽然转录本数目较多，但是可以明显看到存在很多错误的转录本(与1相比剪接模式不同)；4是实施例2获取的转录本，不仅转录本数目较多，并且不存在错误的转录本，说明在准确性和灵敏性上表现优异。由此可见，使用本发明所述过滤方法可以获得更准确、可信、全面的全长转录本数据。
[0088]
实施例2与对比例3分类结果对比如下表所示：
[0089][0090][0091]
由此可见，将具有重叠关系的外显子相似性由70％调节成60％后，灵敏性和准确性都会下降。
[0092]
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。