CN111312331B - 一种利用二代和三代转录组测序数据的基因组注释方法 - Google Patents

Abstract

本发明涉及基因组注释技术领域，提供一种利用二代和三代转录组测序数据的基因组注释方法。本发明的方法包括下述步骤：步骤1：将三代全长转录组测序序列比对到目标基因组，得到各编码基因的初始结构信息；步骤2：将二代转录组测序短序列比对到目标基因组，从比对文件中提取内含子剪切位点信息；步骤3：结合各编码基因的初始结构信息和内含子剪切位点信息，获得各编码基因的最终结构信息。本发明既能充分利用三代转录组无须拼接即可得到全长转录本序列从而准确性高的优点，又能充分考虑二代转录组测序数据能够提供大量的内含子剪切位点证据的优点，大大提高了基因组注释的准确性和效率。

Description

一种利用二代和三代转录组测序数据的基因组注释方法

技术领域

本发明涉及基因组注释技术领域，涉及一种真核生物全基因组编码基因结构的注释方法，特别是涉及一种利用二代和三代转录组测序数据的基因组注释方法。

背景技术

一般全基因组序列组装完毕之后，需要对其蛋白编码基因的结构进行预测。预测通常综合采用三种策略：从头预测、基于近缘物种的序列同源性进行预测以及基于转录组数据的预测。由于转录组是一个物种所表达基因的最直接的证据，因此基于转录组数据的预测结果被认为是最可信的，在整合所有策略的预测结果得到最终基因集的时候该策略所给的权重一般也是最高的。

基于转录组数据预测编码基因的方法中，比较有代表性的有PASA、tophat+cufflinks、hisat2+stringtie等。二代和三代转录组测序数据都可以作为PASA的输入。三代全长转录组因为无须组装即可得到转录本的全长序列，具有准确可靠的优点。通常三代数据可以直接输入给PASA，二代数据则需先用组装软件如Trinity组装后才能输入。而tophat+cufflinks和hisat2+stringtie则比较相似，都是基于二代测序的预测策略，前者可以将测序的短序列比对到基因组上，而后者可以将比对结果信息提取出来从而达到预测编码基因的目的。这两种策略的优点是能充分利用二代测序数据提供的信息，但不能直接利用三代测序数据。PASA的缺点则是不能充分利用二代测序数据所能提供的大量内含子位置信息，且速度较慢。

可见，现有的基因组注释方法各有优缺点，但都无法同时利用二代和三代测序的优点，对基因组预测的精度和效率都不够高。而且随着测序技术的进步，使得每个物种可以利用的二代和三代转录组数据越来越多，如何充分利用这些数据得到更为可靠的预测结果成为迫切需要解决的难题。

发明内容

针对现有技术存在的问题，本发明提供一种利用二代和三代转录组测序数据的基因组注释方法，能够同时利用二代和三代转录组测序数据的优点，大大提高基因组注释的准确性和效率。

本发明的技术方案为：

一种利用二代和三代转录组测序数据的基因组注释方法，其特征在于，包括下述步骤：

步骤1：将三代全长转录组测序序列比对到目标基因组，得到各编码基因的初始结构信息；

步骤2：将二代转录组测序短序列比对到目标基因组，从比对文件中提取内含子剪切位点信息；

步骤3：结合各编码基因的初始结构信息和内含子剪切位点信息，获得各编码基因的最终结构信息。

进一步的，所述步骤1包括下述步骤：

步骤1.1：预测三代全长转录组测序序列的蛋白编码序列：对三代全长转录组测序序列进行分析得到全长非嵌合序列，预测全长非嵌合序列可能的蛋白编码序列；

步骤1.2：将蛋白编码序列比对到目标基因组：将预测得到的所有蛋白编码序列比对到目标基因组，获得比对结果文件；

步骤1.3：提取外显子序列并重新比对回目标基因组：根据比对结果文件及基因组文件提取各外显子序列，将各外显子序列比对回目标基因组，得到各编码基因的初始结构信息。

本发明的有益效果为：

(1)本发明综合利用二代和三代转录组测序数据对基因组进行注释，既能充分利用三代转录组无须拼接即可得到全长转录本序列从而准确性高的优点，又能充分考虑二代转录组测序数据碱基水平的错误率很低同时通量很高从而比对回基因组能够提供大量的内含子剪切位点证据的优点，大大提高了基因组注释的准确性。

(2)本发明的大部分步骤都能并行进行，不需使用mysql数据库，可操作性强，大大提高了基因组注释的效率。

附图说明

图1为本发明的利用二代和三代转录组测序数据的基因组注释方法的原理图。

具体实施方式

下面将结合附图和具体实施方式，对本发明作进一步描述。

本实施例中，使用本发明的基因组注释方法对高原鱼类全基因组基因结构进行注释。如图1所示，本发明的利用二代和三代转录组测序数据的基因组注释方法，包括下述步骤：

步骤1：将三代全长转录组测序序列比对到目标基因组，得到各编码基因的初始结构信息，具体包括下述步骤：

步骤1.1：预测三代全长转录组测序序列的蛋白编码序列：对三代全长转录组测序序列进行分析得到全长非嵌合序列(FLNC)，预测全长非嵌合序列可能的蛋白编码序列。

其中，使用TransDecoder预测全长FLNC的蛋白编码区域，具体步骤如下：

1)提取FLNC序列中最长的开放阅读框(ORF)

TransDecoder.LongOrfs-t target_transcripts.fasta

可以通过-m参数设置ORF的最小长度，一般ORF长度越小，假阳性率会越高，默认值是100。

2)通过Blast对Swissprot蛋白库搜索已知蛋白

blastp-query transdecoder_dir/longest_orfs.pep\

-db uniprot_sprot.fasta-max_target_seqs 1\

-outfmt 6-evalue 1e-5-num_threads 10>blastp.outfmt6

3)预测潜在的蛋白编码区域

TransDecoder.Predict-t target_transcripts.fasta--retain_blastp_hitsblastp.outfmt6

步骤1.2：将蛋白编码序列比对到目标基因组：将预测得到的所有蛋白编码序列比对到目标基因组，获得比对结果文件，具体为：

1)使用blat软件将预测的CDS序列比对到基因组，获得psl格式的比对结果文件

blat genome_seq.fasta CDS_seq.fasta test.psl

2)使用blat2gff(http://eugenes.org/gmod/tandy/blat2gff.pl)将psl格式的比对结果文件转换成gff格式

blat2gff test.psl test.gff

步骤1.3：提取外显子序列并重新比对回目标基因组：根据比对结果文件及基因组文件提取各外显子序列，将各外显子序列比对回目标基因组，得到各编码基因的初始结构信息，具体为：

1)使用GeMoMa根据gff格式的比对结果文件和基因组文件提取各外显子序列

java–jar GeMoMa-1.5.3.jar CLI Extractor a＝test.gffg＝genome_seq.fasta

生成cds-parts.fasta文件

2)对基因组序列建库

makeblastdb-dbtype nucl-in genome_seq.fasta-logfile formatdb.log

3)使用tblastn将各外显子序列比对到基因组

tblastn-outfmt"6std sallseqid score nident positive gaps ppos qframesframe qseq sseq qlen slen salltitles"–db genome_seq.fasta–out test.bla–evalue 1

步骤2：将二代转录组测序短序列比对到目标基因组，从比对文件中提取内含子剪切位点信息，具体为：

1)构建hisat2索引

hisat2-build genome_seq.fasta genome

2)利用hisat2索引将二代转录组测序短序列比对到基因组

Hisat2–x test-1reads1.fq-2reads2.fq|samtools view–bS>hisat2.bam

3)从比对结果中提取内含子剪切位点信息

java-jar GeMoMa-1.5.3.jar CLI ERE c＝true s＝FR_UNSTRANDED m＝hisat2.bam

步骤3：结合各编码基因的初始结构信息和内含子剪切位点信息，获得各编码基因的最终结构信息，具体为：

java-jar GeMoMa-1.5.3.jar CLI GeMoMa t＝test.bla tg＝genome_seq.fastac＝cds-parts.fasta i＝introns.gff得到最终的结果文件predicted_annotation.gff

本发明充分利用了三代转录组无须拼接即可得到全长转录本序列从而准确性高的优点，又充分考虑了二代转录组测序数据提供的大量内含子剪切位点证据，大大提高了基因组注释的准确性。而且本发明的大部分步骤都能并行进行，不需使用mysql数据库，可操作性强，大大提高了基因组注释的效率。

显然，上述实施例仅仅是本发明的一部分实施例，而不是全部的实施例。上述实施例仅用于解释本发明，并不构成对本发明保护范围的限定。基于上述实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，也即凡在本申请的精神和原理之内所作的所有修改、等同替换和改进等，均落在本发明要求的保护范围内。

Claims (1)

1.一种利用二代和三代转录组测序数据的基因组注释方法，其特征在于，包括下述步骤：

步骤1：将三代全长转录组测序序列比对到目标基因组，得到各编码基因的初始结构信息；

步骤2：将二代转录组测序短序列比对到目标基因组，从比对文件中提取内含子剪切位点信息；

步骤3：结合各编码基因的初始结构信息和内含子剪切位点信息，获得各编码基因的最终结构信息；

所述步骤1包括下述步骤：

步骤1.1：预测三代全长转录组测序序列的蛋白编码序列：对三代全长转录组测序序列进行分析得到全长非嵌合序列，预测全长非嵌合序列的蛋白编码序列；

步骤1.2：将蛋白编码序列比对到目标基因组：将预测得到的所有蛋白编码序列比对到目标基因组，获得比对结果文件；

步骤1.3：提取外显子序列并重新比对回目标基因组：根据比对结果文件及基因组文件提取各外显子序列，将各外显子序列比对回目标基因组，得到各编码基因的初始结构信息。

Images