CN116779035B - 多倍体转录组亚基因组分型方法及计算机可读存储介质 - Google Patents
多倍体转录组亚基因组分型方法及计算机可读存储介质 Download PDFInfo
- Publication number
- CN116779035B CN116779035B CN202310605118.6A CN202310605118A CN116779035B CN 116779035 B CN116779035 B CN 116779035B CN 202310605118 A CN202310605118 A CN 202310605118A CN 116779035 B CN116779035 B CN 116779035B
- Authority
- CN
- China
- Prior art keywords
- comparison
- sequence
- transcriptome
- sample
- total
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 108091027544 Subgenomic mRNA Proteins 0.000 title claims abstract description 65
- 208000020584 Polyploidy Diseases 0.000 title claims abstract description 50
- 238000000034 method Methods 0.000 title claims abstract description 30
- 108090000623 proteins and genes Proteins 0.000 claims abstract description 87
- 230000010354 integration Effects 0.000 claims abstract description 22
- 102000004169 proteins and genes Human genes 0.000 claims abstract description 20
- 238000012216 screening Methods 0.000 claims abstract description 20
- 238000001914 filtration Methods 0.000 claims abstract description 19
- 241000894007 species Species 0.000 claims description 55
- 230000008030 elimination Effects 0.000 claims description 16
- 238000003379 elimination reaction Methods 0.000 claims description 16
- 238000012360 testing method Methods 0.000 claims description 7
- 210000001808 exosome Anatomy 0.000 claims description 4
- 241000938605 Crocodylia Species 0.000 claims description 3
- 241000124008 Mammalia Species 0.000 claims description 3
- 238000004590 computer program Methods 0.000 claims description 2
- 108091008053 gene clusters Proteins 0.000 abstract description 2
- 230000014509 gene expression Effects 0.000 description 15
- 238000004364 calculation method Methods 0.000 description 12
- 241000209140 Triticum Species 0.000 description 11
- 235000021307 Triticum Nutrition 0.000 description 11
- 238000012163 sequencing technique Methods 0.000 description 8
- 238000011002 quantification Methods 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 4
- 210000000349 chromosome Anatomy 0.000 description 4
- 150000007523 nucleic acids Chemical group 0.000 description 4
- 238000011144 upstream manufacturing Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 241000209758 Aegilops Species 0.000 description 2
- 241000209148 Aegilops speltoides Species 0.000 description 2
- 241000196324 Embryophyta Species 0.000 description 2
- 108091028043 Nucleic acid sequence Proteins 0.000 description 2
- 240000002805 Triticum turgidum Species 0.000 description 2
- 241000209147 Triticum urartu Species 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 239000012634 fragment Substances 0.000 description 2
- 230000002068 genetic effect Effects 0.000 description 2
- 238000003205 genotyping method Methods 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000005945 translocation Effects 0.000 description 2
- 241001522110 Aegilops tauschii Species 0.000 description 1
- 235000003826 Artemisia Nutrition 0.000 description 1
- 235000003261 Artemisia vulgaris Nutrition 0.000 description 1
- 240000006891 Artemisia vulgaris Species 0.000 description 1
- 241000972773 Aulopiformes Species 0.000 description 1
- 108020004414 DNA Proteins 0.000 description 1
- 241000206602 Eukaryota Species 0.000 description 1
- 208000035199 Tetraploidy Diseases 0.000 description 1
- 244000098338 Triticum aestivum Species 0.000 description 1
- 235000007247 Triticum turgidum Nutrition 0.000 description 1
- 235000014814 Triticum turgidum subsp turgidum Nutrition 0.000 description 1
- 235000009052 artemisia Nutrition 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 210000004027 cell Anatomy 0.000 description 1
- 230000002759 chromosomal effect Effects 0.000 description 1
- 238000012217 deletion Methods 0.000 description 1
- 230000037430 deletion Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 102000054766 genetic haplotypes Human genes 0.000 description 1
- 238000012268 genome sequencing Methods 0.000 description 1
- 238000003780 insertion Methods 0.000 description 1
- 230000037431 insertion Effects 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 230000008774 maternal effect Effects 0.000 description 1
- 108020004999 messenger RNA Proteins 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000001273 protein sequence alignment Methods 0.000 description 1
- 230000010076 replication Effects 0.000 description 1
- 235000019515 salmon Nutrition 0.000 description 1
- 210000001082 somatic cell Anatomy 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000013518 transcription Methods 0.000 description 1
- 230000035897 transcription Effects 0.000 description 1
- 230000002103 transcriptional effect Effects 0.000 description 1
- 238000011222 transcriptome analysis Methods 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- 238000013519 translation Methods 0.000 description 1
Abstract
本发明涉及多倍体转录组亚基因组分型方法及计算机可读存储介质，涉及基因组学领域，包括如下步骤：(1)获取参考物种外群基因序列和待测样本转录组序列，将序列翻译成蛋白序列，进行比对，得到比对结果；(2)设定筛选条件，筛选所述比对结果，获得可信度高的比对结果；(3)将可信度大于预设值的所述比对结果进行连接、整合，得到整合结果；(4)设定过滤参数的阈值，筛选所述整合结果，获得大于等于所述阈值的整合结果并将所述阈值的整合结果按照总得分大小逆序排列；(5)保留总得分最高的一条记录。本发明通过将待测样本转录组序列的CDS序列或者转录本序列比对到参考物种外群基因，来对多倍体转录本进行同源基因簇的划分。
Description
技术领域
本发明涉及基因组学领域，具体涉及多倍体转录组亚基因组分型方法及计算机可读存储介质。
背景技术
多倍体(polyploid)是指体细胞中含有三个或三个以上染色体组的个体，多倍体生物在真核生物中很常见，根据染色体组来源可分为同源多倍体(染色体复制后细胞不分裂所致)和异源多倍体(由不同物种杂交形成)。至今，多倍体起源问题仍是许多物种的研究热点。亚基因组是基因组的子集，在多倍体物种中它是指从一个母本/父本中遗传的一组基因组。如异源多倍体六倍体小麦(AABBDD)，其亚基因组就有三个，亚基因组A，亚基因组B，亚基因组D。多倍体生物在经过测序后各亚基因组信息混杂在一起，需要利用一些方法将各亚基因组信息分开，这就是亚基因组分型。通过亚基因组分型可以确定亚基因组的组成，实现对多倍体基因组的溯源，进而对下游分析奠定基础，比如进化、功能基因分析、多倍化时间鉴定等。
目前亚基因组分型主流方法是基于基因组测序进行分型。比如，利用近源物种作为参考，通过相似性将各亚基因组区分出来；还有利用自身全基因组数据，通过推断各亚基因组的特异性k-mer，将各亚基因组区分开来。
另外，在多倍体转录组分析中，转录组的组装和表达定量都受到了多倍体的影响。在组装上尚无成熟的专门针对多倍体的转录组装方案，为了组装的准确性和完整性，基本都依赖于参考基因组/转录组，而很多非模式多倍体物种缺乏高质量的参考物种基因组/转录组，因此组装质量很低，很难应用于下游分析。
在表达上，多倍体物种存在多种多样的表达，包括有偏向的同源表达(同源基因之间的表达不平衡)、特定条件的差异同源表达、侵袭表达(所有同源基因表达水平或高/低于父母本的表达水平)和优势表达(一个基因的所有同源表达水平大致等于父/母本的表达水平)。如果能够明晰各个亚基因组，则可以实现对单个同源性表达水平的量化，并且可以通过每个同源性集合的总和得出聚合表达，进而分析各种不同的表达模式。
多倍体转录组定量主要面临两个问题：(1)是否有对应多倍体物种的参考基因组/转录组，(2)是同源单倍型间测序片段(reads)的相似程度。现有研究仅讨论了存在参考多倍体物种基因组/转录组信息的情况，总结出一套多倍体转录表达定量的方法。当同源单倍型间reads不存在相似性或相似性不大，则选择传统的比对模式，对应到各自参考基因组/转录组的亚基因组上即可；当reads存在较高的相似性时可选择三种方法：假比对、概率评估、亚基因组分类。假比对的方法通过Salmon、Kallisto等软件，是不基于比对计数而是利用kmer特征推断亚基因组，进而对基因进行定量。亚基因组分类则是对reads映射到参考亚基因组的数目进行统计，并联合亚基因组变异信息，进行映射假设，判断出reads到底映射到哪个亚基因组。
研究发现，虽然已经有多种方法可用于多倍体的表达定量，但是现有方法的策略基本是在准确度和舍弃相似度高的reads之间做取舍，导致组装结果不尽如人意。鉴于此，本发明提供多倍体转录组亚基因组分型方法及计算机可读存储介质。
发明内容
针对现有的将亚基因组分型和表达定量结合到了一起，存在计算资源耗费过大、同源基因分辨率低、定量数据缺省等问题，以及利用基因组测序的方法进行分型，成本又会很高的问题。本发明提供多倍体转录组亚基因组分型方法及计算机可读存储介质。目的是通过将待测样本转录组序列的CDS序列或者转录本序列比对到参考物种外群基因序列，来对多倍体转录组亚基因组进行同源基因簇的划分。
本发明为了解决上述技术问题，第一个方面是提供多倍体转录组亚基因组分型方法，包括如下步骤：
(1)获取参考物种外群基因序列和待测样本转录组序列，将所述参考物种外群基因序列和所述待测样本转录组序列翻译成蛋白序列，进行所述蛋白序列的比对，得到比对结果；
(2)设定筛选条件，根据所述筛选条件筛选所述比对结果，获得可信度大于预设值的比对结果；
上述的筛选条件可以为一致性大于或等于10％且覆盖长度大于或等于10％；
(3)将可信大于预设值的所述比对结果进行连接、整合，得到整合结果；
(4)设定过滤参数的阈值，根据所述阈值筛选所述整合结果，获得大于或等于所述阈值的整合结果，并将所述阈值的整合结果按照得分大小逆序排列；
(5)将大于或等于所述阈值的同一比对的所述整合结果保留总得分最高的一条记录，即为多倍体转录组亚基因组分型的结果。
上述将所述参考物种外群基因序列和所述待测样本转录组序列翻译成蛋白序列，进行所述蛋白序列的比对可以最大程度的保证更多转录本比对到参考物种外群基因序列上，避免遗漏。
本发明的有益效果是：
(1)与传统的基因家族聚类相比，本发明所采用的参考物种外群基因包含了上下游延伸序列，结合了基因的多种可变剪接形式，最大程度上保证了测序序列片段(reads)能够比对上更准确的基因座；由于参考物种外群基因和待测样本转录组序列均为核酸序列，在比对时均被翻译为了多种蛋白序列，这在最大程度上保证了更多reads能够被确定到某个基因上而不被漏掉，再通过后续计算的多项指标来抉择最优比对，由此确定出来的各亚基因组基因代表序列将更为可靠；比如IsoSeq(isoform-sequencing，PacBio三代转录本测序方法)测序数据，因为内含子保留(intron-retention，一种可变剪接形式)或者测序错误导致产生小的indel(insertion and deletion，插入缺失)都可能使转录本无法通过软件翻译为蛋白，而通过tblastx比对则能够更好地保留这部分数据；
(2)本发明在运行上，支持多任务并行，运算效率得到了极大的提升。
综上，本发明基于多倍体物种亚基因组的多种特征，在参考序列的选择、数据比对、分型等方面做了精细的设计，做到了在保证高运算效率的基础上，实现了结果更高的准确性、完整性。
在上述技术方案的基础上，本发明还可以做如下改进。
进一步，步骤(1)中所述参考物种外群基因序列为待测样本近缘物种的基因信息，所述参考物种为已知基因组亚型的二倍体。
上述近缘物种的基因信息可以通过基因库中筛选得到。例如近缘物种的同源基因序列可以通过设定待测样本基因序列相似度阀值，在基因库中筛选超过一定阀值，便可以得到同源基因序列。
上述中参考物种外群基因序列包括基因的延伸序列，上下游2000bp；参考物种应为所有已知基因组亚型的二倍体，因为多倍体参考物种的各亚基因组之间存在基因的渐渗，会影响多倍体样本分型的结果。选择基因延伸序列的原因是大多数的基因注释文件仅给一个基因的一个可变剪接形式，比如，一个基因的第1,3,5个外显子连接的转录本，如果待测样本的这个基因是第1,2,3,4个外显子连接的转录本，结果比对到这个外群，覆盖长度会不高，可能会被过滤掉。
上述中的渐渗是指基因或遗传物质通过群体中的杂种个体与其亲本个体之间的不断回交而导致基因在群体或个体之间转移和传递的过程。发生转移和传递的基因区段称为渐渗区段。而渐渗与易位是因果关系，渐渗就是易位的结果。
进一步，步骤(1)中将所述参考物种外群基因序列和所述待测样本转录组序列翻译成蛋白序列，进行所述蛋白序列的比对具体采用tblastx。
上述的tblastx为BLAST套件的子工具，其将DNA被检索的序列(待测样本转录组序列)和核算序列数据库中的序列(参考物种外群基因序列)按不同的阅读框全部翻译成蛋白质序列，然后进行蛋白质序列比对。
本发明中是在具体比对的应用中将参考序列作为tblastx的subject，而非与reads混合，可避免旁系同源基因与直系同源基因的混淆，也保证了各基因家族内部的直系同源关系，同时，也方便了增加样本数据(同样的流程再跑就可以)。
进一步，步骤(3)包括如下的具体的步骤：
(3.1)将所述参考物种外群基因序列的方向固定为正向，将可信度高的所述比对结果分为正向比对记录和反向比对记录，按照所述待测样本转录组序列的比对位置，连接同向的比对记录，得到所述待测样本转录组序列的正向连接结果序列和所述待测样本序列的反向连接结果序列；
(3.2)设置两个同向的所述比对记录相交的重叠区域的长度占比的阈值，去除冗余，获得去冗余后的所述正向连接结果序列和去冗余后的所述反向连接结果序列；计算每一条去冗余后的所述正向连接结果序列和去冗余后的所述反向连接结果序列的对比长度，同一比对的总一致性，同一比对的总覆盖度，及同一比对的所有比对块的总得分，得到连接结果序列的比对记录。
本发明综合了基因序列的正向和反向链特征，设置两个同向的所述比对记录相交的重叠区域的长度占比的阈值去除冗余(mRNA的可变剪接)，同一比对的总一致性、总覆盖度、及比对块的总得分(亚基因组间基因的相似性)以及分析样本的追加等因素，进行准确的待测样本转录组序列的基因家族聚类，从而达到亚基因组分型的目的。
对上述步骤(2)筛选获得可信度高的比对结果进行连接、整合。但是因为具体的比对过程中的核酸序列参考物种外群基因序列和待测样本转录组序列被翻译成蛋白质，根据核酸序列翻译的相位、方向，比对后的结果可能出现冗余、方向顺序杂乱的情况。此外，同一待测样本转录组序列(query)与同一subject比对时会根据比对到的区块分为很多条比对记录。因此，为了将同一subject与query的上述步骤(2)筛选获得可信度高的比对结果进行连接、整合为一条记录，需要进行上述的步骤(3.1)-(3.2)。
针对步骤(3.1)，由于tblastx是从正负链各三个相位，总共六个相位来翻译的，所以在相同query，subject，相同位置的记录会有多个。因此，首先固定subject的方向为正(subject的起始位置小于subject的终止位置)，然后分别考虑query的两个方向(起始位置小于终止位置；起始位置大于终止位置)，以这样的基础来连接同向的比对记录，并按query比对位置排序。
进一步，步骤(3.2)中，相交的所述重叠区域的长度占比的阈值为大于或等于0.5。
进一步，步骤(3.2)中，所述对比长度为每一条去冗余后的所述正向连接结果序列或去冗余后的所述反向连接结果序列中所有比对块的长度的加和；所述总一致性为每一条去冗余后的所述正向连接结果序列或去冗余后的所述反向连接结果序列中各比对块的长度与所述对比长度的比值，并将各所述比值与相应的比对块的一致性相乘，后求和；所述总覆盖度为所述对比长度与对应去冗余后的所述正向连接结果序列或去冗余后的所述反向连接结果序列的长度的比值；所述总得分为去冗余后的所述正向连接结果序列或去冗余后的所述反向连接结果序列中所有比对块的得分的加和。
若去冗余后的所述正向连接结果序列包括三个比对块，具体的计算过程：例如所述对比长度(AL，alignment length)的计算公式如公式(1)，所述总一致性(I，identity)的计算公式如公式(2)，所述总覆盖度(C，coverage)的计算公式如公式(3)，所述总得分(C，score)(S)的计算公式如公式(4)；
AL＝L1+L2+L3 (1)；
I＝I1*(L1/AL)+I2*(L2/AL)+I3*(L3/AL) (2)；
C＝AL/L (3)；
S＝SL1+SL2+SL3 (4)；
其中，L1、L2、L3分别为同一比对的比对块的长度，I1、I2、I3分别为同一比对的比对块的一致性，L为整条所述待测样本转录组序列的长度，SL1、SL2、SL3分别为同一比对的比对块的得分。
进一步，步骤(4)中，所述过滤参数为总一致性和总覆盖度，所述过滤参数的总一致性的阈值根据不同的待测样本类别和与所述参考物种的亲缘关系远近确定；所述过滤参数的总覆盖度的阈值根据所述待测样本转录组序列的类型确定。
由于步骤(3)得到的整合结果存在多对一的情况，为了保证结果的准确度，需要从总一致性和总覆盖度上进一步进行过滤。
进一步，若为哺乳动物间的比对，总一致性>＝70％；若为两栖爬行动物间的比对，总一致性>＝50％；若为植物间的比对，总一致性>＝50％；若所述参考物种与所述待测样本亲缘关系较远，总一致性>＝30％；若所述待测样本转录组序列的类型为CDS序列，总覆盖度>＝50％；若所述待测样本转录组序列的类型为转录本序列，总覆盖度>＝30％。
进一步，步骤(5)具体为：保留总得分最高的一条记录，即为多倍体转录组亚基因组分型的结果。
第二个方面是提供一种计算机可读存储介质，所述计算机可读存储介上存储有计算机程序，所述计算机程序被处理器执行时实现如上述任一项所述的多倍体转录组亚基因组分型方法。
上述计算机可读介质包括永久性和非永久性、可移动和非可移动，媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存(PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁带磁磁盘存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体(transitory media)，如调制的数据信号和载波。
附图说明
图1为本发明多倍体转录组亚基因组分型方法步骤图；
图2为本发明多倍体转录组亚基因组分型方法具体实施例步骤图；
图3为本发明整合结果图。
具体实施方式
以下对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。
本申请中的亚基因组，是指组成多倍体生物体的基因组中的来自不同祖先的染色体系统，例如多倍体作物的单倍体基因组，如六倍体小麦(AABBDD)，其亚基因组就有三个，亚基因组A，亚基因组B，亚基因组D。
实施例1
本实施例涉及多倍体转录组亚基因组分型方法(如图1)，包括如下步骤：
(1)获取参考物种外群基因序列和待测样本转录组序列，将所述参考物种外群基因序列和所述待测样本转录组序列翻译成蛋白序列，进行所述蛋白序列的比对，得到比对结果；
优选的，步骤(1)中所述参考物种外群基因序列为待测样本近缘物种的基因信息，所述参考物种为已知基因组亚型的二倍体。
上述中参考物种外群基因序列包括基因的延伸序列，上下游2000bp；参考物种应为所有已知基因组亚型的二倍体，因为多倍体参考物种的各亚基因组之间存在基因的渐渗，会影响多倍体样本分型的结果。选择基因延伸序列的原因是大多数的基因注释文件仅给一个基因的一个可变剪接形式(比如，一个基因的第1,3,5个外显子连接的转录本)，如果待测样本的这个基因是第1,2,3,4个外显子连接的转录本，结果比对到这个外群，覆盖长度(coverage)会不高，可能会被过滤掉。
优选的，步骤(1)中将所述参考物种外群基因序列和所述待测样本转录组序列翻译成蛋白序列，进行所述蛋白序列的比对具体可以采用tblastx。
具体的(如图2)，a)tblastx比对：参考物种：乌拉尔图小麦(Triticum urartu)，二倍体，A亚基因组+拟斯皮尔脱山羊草(Aegilops speltoides Tausch)，二倍体，B亚基因组+节节麦(Aegilops tauschii)，二倍体，D亚基因组。
样本：中国春小麦(Triticum aestivum)，六倍体，AABBDD。
研究表明中国春小麦的A、B、D亚基因组分别来自于乌拉尔图小麦、拟斯皮尔脱山羊草和节节麦的多次多种杂交，因此选取这三个物种作为参考物种。
根据各参考物种的基因组注释信息，提取了其所有基因及其上下游2000bp的核酸序列，作为tblastx的subject序列。选取中国春小麦的全基因组序列作为tblastx的query序列。
(2)设定一致性大于或等于10％且覆盖长度大于或等于10％的筛选条件，根据所述筛选条件筛选所述比对结果，获得可信的比对结果；
具体的(如图2)，b)初步处理tblastx结果：
过滤掉一致性低于10％和覆盖度低于10％的比对结果，获得可信的比对结果；
(3)将可信度高的所述比对结果进行连接、整合，得到整合结果；
优选地，步骤(3)包括如下的具体的步骤：
(3.1)将所述参考物种外群基因序列的方向固定为正向，将可信度高的所述比对结果分为正向比对记录和反向比对记录，按照所述待测样本转录组序列的比对位置，连接同向的比对记录，得到所述待测样本转录组序列的正向连接结果序列和所述待测样本序列的反向连接结果序列；
(3.2)设置两个同向的所述比对记录相交的重叠区域的长度占比的阈值，去除冗余，获得去冗余后的所述正向连接结果序列和去冗余后的所述反向连接结果序列；计算每一条去冗余后的所述正向连接结果序列和去冗余后的所述反向连接结果序列的对比长度，同一比对的总一致性，同一比对的总覆盖度，及同一比对的所有比对块的总得分，得到连接结果序列的比对记录；
优选的，步骤(3.2)中，相交所述重叠区域的长度占比的阈值为大于或等于0.5。
优选的，步骤(3.2)中，所述对比长度为每一条去冗余后的所述正向连接结果序列或去冗余后的所述反向连接结果序列中所有比对块的长度的加和；所述总一致性为每一条去冗余后的所述正向连接结果序列或去冗余后的所述反向连接结果序列中各比对块的长度与所述对比长度的比值，并将各所述比值与相应的比对块的一致性相乘，后求和；所述总覆盖度为所述对比长度与对应去冗余后的所述正向连接结果序列或去冗余后的所述反向连接结果序列的长度的比值；所述总得分为去冗余后的所述正向连接结果序列或去冗余后的所述反向连接结果序列中所有比对块的得分的加和。
若去冗余后的所述正向连接结果序列包括三个比对块(如图3)，具体的计算过程：例如所述对比长度(AL，alignment length)的计算公式如公式(1)，所述总一致性(I，identity)的计算公式如公式(2)，所述总覆盖度(C，coverage)的计算公式如公式(3)，所述总得分(C，score)(S)的计算公式如公式(4)；
AL＝L1+L2+L3 (1)；
I＝I1(L1/AL)+I2(L2/AL)+I3*(L3/AL) (2)；
C＝AL/L (3)；
S＝SL1+SL2+SL3 (4)；
其中，L1、L2、L3分别为同一比对的比对块的长度，I1、I2、I3分别为同一比对的比对块的一致性，L为整条所述待测样本转录组序列的长度，SL1、SL2、SL3分别为同一比对的比对块的得分。
具体的(如图2)，c)连接分段比对记录：
i.固定subject的方向为正向，将上述结果分为正向比对记录(query为正向)和反向(query为反向)比对记录，并按query的坐标顺序排序；
ii.设置重叠区段阈值为0.5，过滤掉低于0.5的记录，并保留得分最高的记录；
iii.通过计算获得query的每条序列的长度，通过计算同一query和同一subject比对(后面称为该比对)的所有区块的长度百分比，并与各自的比对一致性相乘，求和后将其作为该比对的总一致性；另外将总的比对块的长度求和后计算占该query序列的百分比，并作为该比对的总覆盖度；后将该比对的所有比对块的得分加起来，作为该比对的总得分。
(4)设定过滤参数的阈值，根据所述阈值筛选所述整合结果，获得大于或等于所述阈值的整合结果，并将所述阈值的整合结果按照总得分大小逆序排列；
优选的，步骤(4)中，所述过滤参数为总一致性和总覆盖度，所述过滤参数的总一致性的阈值根据不同的待测样本类别和与所述参考物种的亲缘关系远近确定；所述过滤参数的总覆盖度的阈值根据所述待测样本转录组序列的类型确定。
由于步骤(3)得到的整合结果存在多对一的情况，为了保证结果的准确度，需要从总一致性和总覆盖度上进一步进行过滤。
优选的，若为哺乳动物间的比对，总一致性>＝70％；若为两栖爬行动物间的比对，总一致性>＝50％；若为植物间的比对，总一致性>＝50％；若所述参考物种与所述待测样本亲缘关系较远，总一致性>＝30％；若所述待测样本转录组序列的类型为CDS序列，总覆盖度>＝50％；若所述待测样本转录组序列的类型为转录本序列，总覆盖度>＝30％。
具体的(如图2)，d)过滤并将结果按总得分逆序排序：
将前面处理好的正向和反向比对记录合并，过滤掉一致性低于30％和覆盖度低于30％的比对记录，并以总得分逆序排列；
(5)将大于或等于所述阈值的同一比对的所述整合结果保留总得分最高的一条记录，即为多倍体转录组亚基因组分型的结果。
优选的，步骤(5)具体为：保留总得分最高的一条记录，即为多倍体转录组亚基因组分型的结果。
具体的(如图2)，e)同源基因聚类：
通过计算保留每同一比对(同一query和同一subject比对)中得分最高的记录，去除多个相同的最高得分比对记录，最后经整理得到中国春小麦基因的聚类情况，从而获得中国春小麦的亚基因组情况。
本实例共得到51072条聚类结果，其中真阳性率(正确聚类到对应亚基因组的比率，即(queryA-subjectA+queryB-subjectB+queryD-subjectD)/subjectABD对应到的总聚类条数)为45.32％；假阳性率为54.68％。由于A、B和D同源性较高，本身区分难度较大；此外，多倍体化发生后，染色体组间存在结构变异，导致亚基因组染色体上存在外源亚基因组基因。因此假阳性率较高属正常现象。
实施例2
参考物种：乌拉尔图小麦(Triticum urartu)，二倍体，A亚基因组+拟斯皮尔脱山羊草(Aegilops speltoides Tausch)，二倍体，B亚基因组。
样本：圆锥小麦(Triticum turgidum L.ssp)，四倍体，AABB。
步骤同上实施例1的记载。
本实施例共得到41311条聚类结果，其中真阳性率(正确聚类到对应亚基因组的比率)为65.58％，假阳性率为34.42％。
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。
Claims (5)
1.多倍体转录组亚基因组分型方法，其特征在于，包括如下步骤：
（1）获取参考物种外群基因序列和待测样本转录组序列，将所述参考物种外群基因序列和所述待测样本转录组序列翻译成蛋白序列，进行所述蛋白序列的比对，得到比对结果；
（2）设定筛选条件，根据所述筛选条件筛选所述比对结果，获得可信度大于预设值的比对结果；
（3）将可信度大于预设值的所述比对结果进行连接、整合，得到整合结果；
（4）设定过滤参数的阈值，根据所述阈值筛选所述整合结果，获得大于或等于所述阈值的整合结果，并将所述阈值的整合结果按照总得分大小降序排列；
（5）将大于或等于所述阈值的同一比对的所述整合结果保留总得分最高的一条记录，即为多倍体转录组亚基因组分型的结果；
步骤（1）中所述参考物种外群基因序列为待测样本近缘物种的基因信息，所述参考物种为已知基因组亚型的二倍体；
步骤（3）包括如下的具体的步骤：
（3.1）将所述参考物种外群基因序列的方向固定为正向，将可信度高的所述比对结果分为正向比对记录和反向比对记录，按照所述待测样本转录组序列的比对位置，连接同向的比对记录，得到所述待测样本转录组序列的正向连接结果序列和所述待测样本转录组序列的反向连接结果序列；
（3.2）设置两个同向的所述比对记录相交的重叠区域的长度占比的阈值，去除冗余，获得去冗余后的所述正向连接结果序列和去冗余后的所述反向连接结果序列；计算每一条去冗余后的所述正向连接结果序列和去冗余后的所述反向连接结果序列的对比长度，同一比对的总一致性，同一比对的总覆盖度，及同一比对的所有比对块的总得分，得到连接结果序列的比对记录；
步骤（3.2）中，相交的所述重叠区域的长度占比的阈值为大于或等于0.5；
步骤（3.2）中，所述对比长度为每一条去冗余后的所述正向连接结果序列或去冗余后的所述反向连接结果序列中所有比对块的长度的加和；所述总一致性为每一条去冗余后的所述正向连接结果序列或去冗余后的所述反向连接结果序列中各比对块的长度与所述对比长度的比值，并将各所述比值与相应的比对块的一致性相乘，后求和；所述总覆盖度为所述对比长度与对应去冗余后的所述正向连接结果序列或去冗余后的所述反向连接结果序列的长度的比值；所述总得分为去冗余后的所述正向连接结果序列或去冗余后的所述反向连接结果序列中所有比对块的得分的加和；
步骤（4）中，所述过滤参数为总一致性和总覆盖度，所述过滤参数的总一致性的阈值根据不同的待测样本类别和与所述参考物种的亲缘关系远近确定；所述过滤参数的总覆盖度的阈值根据所述待测样本转录组序列的类型确定。
2.根据权利要求1所述多倍体转录组亚基因组分型方法，其特征在于，步骤（1）中将所述参考物种外群基因序列和所述待测样本转录组序列翻译成蛋白序列，进行所述蛋白序列的比对具体采用tblastx。
3.根据权利要求1所述多倍体转录组亚基因组分型方法，其特征在于，若为哺乳动物间的比对，总一致性>=70%；若为两栖爬行动物间的比对，总一致性>=50%；若为植物间的比对，总一致性>=50%；若所述参考物种与所述待测样本亲缘关系较远，总一致性>=30%；若所述待测样本转录组序列的类型为CDS序列，总覆盖度>=50%；若所述待测样本转录组序列的类型为转录本序列，总覆盖度>=30%。
4.根据权利要求1所述多倍体转录组亚基因组分型方法，其特征在于，步骤（5）具体为：保留总得分最高的一条记录，即为多倍体转录组亚基因组分型的结果。
5.一种计算机可读存储介质，其特征在于，所述计算机可读存储介上存储有计算机程序，所述计算机程序被处理器执行时实现如权利要求1至4中任一项所述的多倍体转录组亚基因组分型方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310605118.6A CN116779035B (zh) | 2023-05-26 | 2023-05-26 | 多倍体转录组亚基因组分型方法及计算机可读存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310605118.6A CN116779035B (zh) | 2023-05-26 | 2023-05-26 | 多倍体转录组亚基因组分型方法及计算机可读存储介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116779035A CN116779035A (zh) | 2023-09-19 |
| CN116779035B true CN116779035B (zh) | 2024-03-15 |
Family
ID=87990505
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310605118.6A Active CN116779035B (zh) | 2023-05-26 | 2023-05-26 | 多倍体转录组亚基因组分型方法及计算机可读存储介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116779035B (zh) |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009121091A1 (en) * | 2008-04-04 | 2009-10-08 | Molecular Plant Breeding Nominees Ltd | Mapping method for polyploid subjects |
| CN101974629A (zh) * | 2010-10-26 | 2011-02-16 | 西南大学 | 一种虚拟合成物种调查异源多倍体物种起源的方法 |
| WO2013103759A2 (en) * | 2012-01-04 | 2013-07-11 | Dow Agrosciences Llc | Haplotype based pipeline for snp discovery and/or classification |
| WO2018232580A1 (zh) * | 2017-06-20 | 2018-12-27 | 深圳华大基因研究院 | 基于三代捕获测序对二倍体基因组单倍体分型的方法和装置 |
| CN111445953A (zh) * | 2020-03-27 | 2020-07-24 | 武汉古奥基因科技有限公司 | 一种利用全基因组比对拆分四倍体鱼类亚基因组的方法 |
| WO2020257719A1 (en) * | 2019-06-21 | 2020-12-24 | Coopersurgical, Inc. | Systems and methods for determining genome ploidy |
| CN113496760A (zh) * | 2020-04-01 | 2021-10-12 | 深圳华大基因科技服务有限公司 | 基于第三代测序的多倍体基因组组装方法和装置 |
-
2023
- 2023-05-26 CN CN202310605118.6A patent/CN116779035B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009121091A1 (en) * | 2008-04-04 | 2009-10-08 | Molecular Plant Breeding Nominees Ltd | Mapping method for polyploid subjects |
| CN101974629A (zh) * | 2010-10-26 | 2011-02-16 | 西南大学 | 一种虚拟合成物种调查异源多倍体物种起源的方法 |
| WO2013103759A2 (en) * | 2012-01-04 | 2013-07-11 | Dow Agrosciences Llc | Haplotype based pipeline for snp discovery and/or classification |
| WO2018232580A1 (zh) * | 2017-06-20 | 2018-12-27 | 深圳华大基因研究院 | 基于三代捕获测序对二倍体基因组单倍体分型的方法和装置 |
| WO2020257719A1 (en) * | 2019-06-21 | 2020-12-24 | Coopersurgical, Inc. | Systems and methods for determining genome ploidy |
| CN111445953A (zh) * | 2020-03-27 | 2020-07-24 | 武汉古奥基因科技有限公司 | 一种利用全基因组比对拆分四倍体鱼类亚基因组的方法 |
| CN113496760A (zh) * | 2020-04-01 | 2021-10-12 | 深圳华大基因科技服务有限公司 | 基于第三代测序的多倍体基因组组装方法和装置 |
Non-Patent Citations (6)
| Title |
|---|
| Cytochrome P450 monooxygenase genes in the wild silkworm.Bombyx mandarina.2021,全文. * |
| Tracking the ancestry of known and ‘ghost’ homeologous subgenomes in model grass Brachypodium polyploids;Ruben Sancho等;The Plant Journal;全文 * |
| 同源区段靶向测序数据的基因型鉴定与分析流程开发;武建楠;中国硕士学位论文全文库 基础科学辑;全文 * |
| 天然异源四倍体甘蓝型油菜中亚基因组的分离及互作;朱斌;中国硕士学位论文全文库 农业科技辑;全文 * |
| 植物多倍化与多倍体基因组进化研究进展;李霖锋;刘宝;;中国科学:生命科学(第04期);全文 * |
| 鸡lncRNA-MSTRG.15568.9及其预测靶基因的表达;黄子妍;中国畜牧兽医;第46卷(第7期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116779035A (zh) | 2023-09-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Zhang et al. | Assembly of allele-aware, chromosomal-scale autopolyploid genomes based on Hi-C data | |
| Marchant et al. | The C-Fern (Ceratopteris richardii) genome: insights into plant genome evolution with the first partial homosporous fern genome assembly | |
| Kellis et al. | Methods in comparative genomics: genome correspondence, gene identification and regulatory motif discovery | |
| CN107423578B (zh) | 检测体细胞突变的装置 | |
| CN110010193B (zh) | 一种基于混合策略的复杂结构变异检测方法 | |
| CN104302781B (zh) | 一种检测染色体结构异常的方法及装置 | |
| CN104204221A (zh) | 一种检验融合基因的方法及系统 | |
| WO2017143585A1 (zh) | 对分隔长片段序列进行组装的方法和装置 | |
| CN110739027A (zh) | 一种基于染色质区域覆盖深度的癌症组织定位方法及系统 | |
| CN107480470A (zh) | 基于贝叶斯与泊松分布检验的已知变异检出方法和装置 | |
| Forsythe et al. | Biased gene retention in the face of introgression obscures species relationships | |
| CN113362889A (zh) | 基因组结构变异注释方法 | |
| CN116779035B (zh) | 多倍体转录组亚基因组分型方法及计算机可读存储介质 | |
| CN114530200B (zh) | 基于计算snp熵值的混合样本鉴定方法 | |
| Sun et al. | Biased mutations and gene losses underlying diploidization of the tetraploid broomcorn millet genome | |
| CN115148289B (zh) | 同源四倍体基因组分型组装的方法和装置、构建染色体的装置 | |
| CN115992261A (zh) | 一种构建牡蛎单倍型基因组的方法 | |
| Molina et al. | ddRADseq-mediated detection of genetic variants in sugarcane | |
| CN113035279A (zh) | 基于miRNA测序数据的帕金森疾病演化关键模块识别方法 | |
| CN114913918A (zh) | 一种针对孤独症的高通量测序数据分析方法及装置 | |
| CN109935274A (zh) | 一种基于k-mer分布特征的长读数重叠区检测方法 | |
| CN116072222B (zh) | 病毒基因组鉴定和拼接的方法及应用 | |
| Howe et al. | Illumina sequencing artifacts revealed by connectivity analysis of metagenomic datasets | |
| CN114067909B (zh) | 一种矫正同源重组缺陷评分的方法、装置和存储介质 | |
| CN111554347B (zh) | 构建用于手足口样本归类的模型的方法及其应用 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |