CN110534157B - 一种批量提取基因组基因信息并翻译比对分析序列的方法 - Google Patents

一种批量提取基因组基因信息并翻译比对分析序列的方法 Download PDF

Info

Publication number
CN110534157B
CN110534157B CN201910684539.6A CN201910684539A CN110534157B CN 110534157 B CN110534157 B CN 110534157B CN 201910684539 A CN201910684539 A CN 201910684539A CN 110534157 B CN110534157 B CN 110534157B
Authority
CN
China
Prior art keywords
file
gene
pro
res
sequence
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
Application number
CN201910684539.6A
Other languages
English (en)
Other versions
CN110534157A (zh
Inventor
郭月
刘静
胡茂龙
浦惠明
张洁夫
龙卫华
张维
周晓婴
孙程明
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Jiangsu Academy of Agricultural Sciences
Original Assignee
Jiangsu Academy of Agricultural Sciences
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Jiangsu Academy of Agricultural Sciences filed Critical Jiangsu Academy of Agricultural Sciences
Priority to CN201910684539.6A priority Critical patent/CN110534157B/zh
Publication of CN110534157A publication Critical patent/CN110534157A/zh
Application granted granted Critical
Publication of CN110534157B publication Critical patent/CN110534157B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/30—Detection of binding sites or motifs
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • G16B25/10—Gene or protein expression profiling; Expression-ratio estimation or normalisation
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/10—Sequence alignment; Homology search
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Genetics & Genomics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Medical Informatics (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Molecular Biology (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
  • Apparatus Associated With Microorganisms And Enzymes (AREA)

Abstract

本发明公开了一种批量提取基因组基因信息并翻译比对分析序列的方法。本发明所提供的批量提取基因组基因信息并翻译比对分析序列的方法综合运用了基于多序列比对分析的MUSCLE程序，并结合多个Perl脚本语言编程的方法。实验证明，本发明所提供的批量提取基因组基因信息并翻译比对分析序列的方法比较系统，能够完成目标基因序列和转录本序列的提取、目标基因或转录本的基因组关键信息获取、DNA序列翻译以及翻译后序列的多重比较，获取各相关结果文件的重复效果好，速度快，易实现批量化、自动化和流程化。

Description

一种批量提取基因组基因信息并翻译比对分析序列的方法
技术领域
本发明属于生物技术领域，涉及一种批量提取基因组基因信息并翻译比对分析序列的方法。
背景技术
人类通过数千年的驯化和近代以来有计划性的育种，形成了当今多样化的动物和作物品种，从而提供丰富的动植物源性蛋白满足人类需求。在过去的100年里，数量遗传学应用于动植物育种领域引发了相关育种技术的革命，但动物和作物机体遗传发育体系相当复杂，一些性状仍然难以通过基于系谱的育种值进行高效选育，遗传潜能尚未充分发掘。人类基因组计划带来的理念和技术极大促进了动植物基因组学的发展，使得人们可以从全基因组水平精准定位功能变异，挖掘功能元件的生物学意义，为动物和作物分子设计育种提供重要的理论基础。
然而基因组学的飞速发展同时带来的却是基因组数据爆炸式的增长，进而导致数据孤岛与数据海洋等问题日益严重。育种家如何从海量的基因组数据中提取对他们有帮助的数据信息难上加难。本专利的目的在于帮助生物学科研人员尤其是育种家们，便捷高效提取利用辅助其科研工作的基因组信息。着力打造批量化、流程化和自动化的简单便捷、通俗易懂、容易上手的一站式操作方法。
发明内容
本发明的目的是提供一种批量提取基因组基因信息并翻译比对分析序列的方法。根据某一物种的转录本ID或者基因ID，依据供试基因组cds文件、蛋白质文件、gff文件和染色体fasta文件等信息，通过6个perl脚本程序，实现目标转录本或基因在基因组上的位置、长度、正反义链等结构信息的提取，并在染色体fasta文件上提取该转录本或基因的cds或基因序列，在基因组蛋白文件上提取该转录本的蛋白质序列。最后对所需cds序列进行翻译，或者直接用所获的蛋白质序列，调用Linux系统程序完成蛋白质的多序列比对工作。
上述方法具体包括如下步骤：
(1)建立工作文件夹work_dir，将待测某一物种的转录本ID文件记为A数据集A，所述A数据集A的文件名为“XXX1”，运行“perl script1.pl XXX1”命令，在当前工作文件夹work_dir下得到“res_Gene_ID”文件；所述“XXX1”在运行“perl script1.pl XXX1”程序时已置于包含脚本“script1.pl”的当前工作文件夹work_dir内，关闭所有相关文件；所述“res_Gene_ID”文件为该物种转录本ID对应的基因ID文件，记为数据集B，命名为“XXX3”；
如果上述步骤直接提供的是某一物种基因ID，则将该基因ID文件记为数据集B，命名为“XXX3”。
(2)将该物种基因组gff文件记为C数据集，所述C数据集的文件名为“XXX2”，运行“perl script2.pl XXX2XXX3”命令，在当前工作文件夹work_dir下得到“res_Geneinfo”文件；
所述“res_Geneinfo”文件为根据该物种基因ID文件提取的基因组信息文件，记为数据集D；所述“XXX2”、“XXX3”在运行“perl script2.pl XXX2XXX3”程序时已置于包含脚本“script2.pl”的当前工作文件夹work_dir内，关闭所有相关文件。
(3)为Strawberry Perl软件安装Bioperl模块，将该物种基因组cds的fasta格式文件记为数据集E，所述数据集E的文件名为“XXX4”，运行“perl script3.pl XXX1”命令，在当前工作文件夹work_dir下得到“res_CDS_seq”文件；
所述“res_CDS_seq”文件为根据该物种转录本ID文件提取的基因cds序列fasta文件，记为数据集G；所述“XXX4”在运行“perl script3.pl XXX1”程序时已置于包含脚本“script3.pl”的当前工作文件夹work_dir内，关闭所有相关文件。
(4)将该物种基因组染色体的fasta格式文件记为数据集F，所述数据集F的文件名为“XXX5”，运行“perl script4.pl res_Geneinfo”命令，在当前工作文件夹work_dir下得到“res_Gene_seq”文件；
所述“res_Gene_seq”文件为根据该物种基因ID文件从该物种基因组染色体文件中提取的基因序列fasta文件，记为数据集H；所述“XXX5”在运行“perl script4.pl res_Geneinfo”程序时已置于包含脚本“script4.pl”的当前工作文件夹work_dir内，关闭所有相关文件。
(5)在当前工作文件夹work_dir内运行“perl script 5.pl”命令，得到“PRO_1st.fa”、“PRO_2nd.fa”、“PRO_3rd.fa”、“PRO_RC_1st.fa”、“PRO_RC_2nd.fa”、“PRO_RC_3rd.fa”和“PRO_last.fa”7个文件；
所述“PRO_1st.fa”、“PRO_2nd.fa”、“PRO_3rd.fa”、“PRO_RC_1st.fa”、“PRO_RC_2nd.fa”和“PRO_RC_3rd.fa”6个文件为根据该物种基因ID文件提取的基因序列或者转录本cds序列翻译后的蛋白质fasta文件，分别记为数据集I、J、K、L、M和N；所述“PRO_last.fa”文件为筛选出的用于后续多序列比对计算的蛋白质序列文件，记为数据集O；所述“res_CDS_seq”文件在运行“perl script 5.pl”程序时已置于包含脚本“perl script 5.pl”的当前工作文件夹work_dir内，关闭所有相关文件。
(6)如果通过下载获得该物种基因组蛋白质的fasta格式文件，则将其记为P数据集，所述P数据集的文件名为“XXX6”，运行“perl script6.pl XXX1”命令，在当前工作文件夹work_dir下得到“res_PRO_seq”文件；
所述“res_PRO_seq”文件为根据该物种转录本ID文件提取的基因蛋白质序列fasta文件，记为数据集Q；所述“XXX6”在运行“perl script6.pl XXX1”程序时已置于包含脚本“script6.pl”的当前工作文件夹work_dir内，关闭所有相关文件。
(7)在当前工作文件夹work_dir内运行“muscle-inPRO_last.fa–outPRO_last.out”命令，如果存在上述步骤(6)，则运行“muscle-inres_PRO_seq–outres_PRO_seq.out”命令，在当前工作文件夹中得到多重序列比对的结果文件；
所述“PRO_last.out”和“res_PRO_seq.out”文件为MUSCLE软件计算后的输出文件，记为数据集R；且在运行“muscle-inPRO_last.fa–out PRO_last.out”命令或者“muscle-in res_PRO_seq–out res_PRO_seq.out”命令后所产生的结果文件在当前工作文件夹work_dir中，关闭所有相关文件。
在上述方法步骤(1)中，所述脚本“script1.pl”中关于获得“res_Gene_ID”文件的内容，是基于如下原理进行编程的：While循环对“XXX1”文件进行逐行处理，对每行进行模式匹配，将Bn开头到“.”符号前的基因ID进行提取并存入变量$gene_id中，将结果打印至同一文件中，文件名即为“res_Gene_ID”，同时把该文件置于当前工作目录work_dir文件夹中，关闭所有相关文件。
在上述方法步骤(2)中，所述脚本“script2.pl”中关于获得“res_Geneinfo”文件的内容，是基于如下原理进行编程的：将res_Gene_ID文件读入数组@name_can中，打开供试基因组gff文件“XXX2”，while循环逐条处理并分割“XXX2”文件。模式匹配识别“mRNA”标识的行并提取该行的基因ID到变量$id_tmp，for循环遍历数组@name_can的每一行，当变量$id_tmp与数组某行基因ID相同时，计算该基因的长度并存入到变量$genelen中，把基因ID、所在染色体号、基因的起始位点、终止位点、基因长度以及正反义链等信息，逐行打印至同一文件，文件名为“res_Geneinfo”，同时把该文件置于当前工作目录work_dir文件夹中，关闭所有相关文件。
在上述方法步骤(3)中，所述脚本“script3.pl”中关于获得“res_CDS_seq”文件的内容，是基于如下原理进行编程的：利用Bio::SeqIO模块和while循环将供试基因组cds文件“XXX4”逐条读入哈希％hash中，打开供试转录本ID文件“XXX1”，While循环对“XXX1”文件进行逐行处理，if判别如果存在以“XXX1”文件中某行的转录本ID为key值的哈希value值$hash{$line}，则去掉$hash{$line}后的最后一个“*”号，并将转录本ID以及对应的哈希value值即cds序列，以fasta的格式逐条打印至同一结果文件中，文件名为“res_CDS_seq”，else条件如果不存在上述哈希value值，则在屏幕输出该转录本ID未找到。把该结果文件“res_CDS_seq”置于当前工作目录work_dir文件夹中，关闭所有相关文件。
在上述方法步骤(4)中，所述脚本“script4.pl”中关于获得“res_Gene_seq”文件的内容，是基于如下原理进行编程的：利用Bio::SeqIO模块和while循环将供试基因组染色体文件“XXX5”逐条读入哈希％hash中，打开文件“res_Geneinfo”，While循环对其进行逐行处理，next if语句去掉以字母“G”开头的行后逐行分割文件，通过substr函数依靠文件中基因的起始、终止位置变量和基因长度变量$row[1]、$row[2]和$row[4]，提取位于染色体$hash{$row[1]}上的基因序列，并存入变量$seq_tmp中。If判别如果该基因的方向是反义链“-”，则将该序列的反向互补序列求出，存于变量$seq_tmp中。最后将所有结果以基因ID对应序列的fasta文件格式打印至同一文件中，文件名为“res_Gene_seq”，同时把该文件置于当前工作目录work_dir文件夹中，关闭所有相关文件。
在上述方法步骤(5)中，所述脚本“script5.pl”中关于获得“PRO_last.fa”文件的内容，是基于如下原理进行编程的：首先把20种氨基酸的64种密码子在程序里面写完整并存入哈希％genetic_code中。打开待翻译DNA序列的fasta文件，依靠Bio::SeqIO模块接收供试DNA序列。
然后利用while循环逐条读取输入文件DNA序列，利用uc函数将序列字母转换为大写，利用reverse函数及正则表达式tr///获取读取DNA序列的反向互补序列，length函数计算序列长度。利用存储密码子简表的哈希％genetic_code，分别从读取DNA序列起始位置的第一、二和三位密码子开始向后进行翻译(以三个相连密码子为翻译单位)，将翻译后的蛋白质序列及其ID以fasta文件格式写入结果文件PRO_1st.fa、PRO_2nd.fa和PRO_3rd.fa中；同时分别从计算所得DNA序列反向互补序列起始位置的第一、二和三位密码子开始向后进行翻译(以三个相连密码子为翻译单位)，将翻译后的蛋白质序列及其ID以fasta文件格式写入结果文件PRO_RC_1st.fa、PRO_RC_2nd.fa和PRO_RC_3rd.fa中，同时将所得6个结果文件置于当前工作目录work_dir文件夹中，关闭所有相关文件。
第三，stat函数分别取6个结果文件的文件大小，存入数组@array_size中，并分别以文件大小和文件名为哈希％hash_size的key值和value值；将数组中元素按着从大到小的顺序排序后存入新数组@array_sort中，然后筛选出@array_sort中最大的元素$array_sort[0]，以及以该最大元素为key值对应的哈希value值$hash_size{$array_sort[0]}存入变量$file_biggest中，最后将$file_biggest文件打开，并利用Bio::SeqIO模块嵌套while循环，将该文件内容逐行打印至结果文件“PRO_last.fa”中，同时把该文件置于当前工作目录work_dir文件夹中，关闭所有相关文件。
在上述方法步骤(6)中，所述脚本“script6.pl”中关于获得“res_PRO_seq”文件的内容，是基于如下原理进行编程的：
利用Bio::SeqIO模块和while循环将供试基因组蛋白质文件“XXX6”逐条读入哈希％hash中，打开供试转录本ID文件“XXX1”，While循环对“XXX1”文件进行逐行处理，if判别如果存在以“XXX1”文件中某行的转录本ID为key值的哈希value值$hash{$line}，则去掉$hash{$line}后的最后一个“*”号，并将转录本ID以及对应的哈希value值即蛋白质序列，以fasta的格式逐条打印至同一结果文件中，文件名为“res_PRO_seq”，else条件如果不存在上述哈希value值，则在屏幕输出该转录本ID未找到。把该结果文件“res_PRO_seq”置于当前工作目录work_dir文件夹中，关闭所有相关文件。
进一步地，本发明中所述脚本“script1.pl”具体为：
进一步地，本发明中所述脚本“script2.pl”具体为：
进一步地，本发明中所述脚本“script3.pl”具体为：
进一步地，本发明中所述脚本“script4.pl”具体为：
进一步地，本发明中所述脚本“script5.pl”具体为：
/>
/>
进一步地，本发明中所述脚本“script6.pl”具体为：
/>
在本发明中，步骤(1)中的所述待测物种为任意物种。
在所述方法中，所述待测基因组染色体、cds编码区、蛋白质序列以及基因组信息gff文件可以通过下载已公开的全基因组序列获得或通过全基因组测序得到。具体地，本发明所述待测基因组具体为油菜(Brassica napus)中双11基因组。所述油菜基因组记录于中国农业科学院油料作物研究所的油料作物基因组数据库(http://ocri-genomics.org/ Brassia_napus_genome_ZS11/)
本发明具有以下优点：
一是为计算批量提取基因组基因信息并翻译比对分析序列提供了一种简单便捷、批量高效的方法；提取基因组基因信息较精确，对基因cds序列的翻译和比对较为全面和准确，效果好，速度快；
二是能够高效整合各基因组相关信息并进行全面系统的解析，能一次性获得批量的有用有效结果，易实现流程化、批量化、自动化；
三是本发明将高效常用的多序列比对软件与多个Perl脚本语言编程完美流畅的结合起来，实现了软件之间的良好衔接，在很大程度上弥补了提取基因组基因信息以及进行翻译比对分析序列过程中出现的耗时费力、流程化批量化欠缺等不足。
本方法在提取基因组基因信息以及翻译比对分析序列过程中可以发挥重要的作用。
附图说明
图1为本发明批量高效提取基因组基因信息并翻译比对分析序列的方法流程图。
图2为实施例2中步骤1)中Perl脚本“script1.pl”分析后所获“res_Gene_ID”文件格式的图示。
图3为实施例2中步骤2)中Perl脚本“script2.pl”分析后所获“res_Geneinfo”文件格式的图示。
图4为实施例2中步骤3)中Perl脚本“script3.pl”分析后所获“res_CDS_seq”文件格式的图示。
图5为实施例2中步骤4)中Perl脚本“script4.pl”分析后所获“res_Gene_seq”文件格式的图示。
图6为实施例2中步骤5)中Perl脚本“script5.pl”分析后所获“PRO_1st.fa”的图示。
图7为实施例2中步骤5)中Perl脚本“script5.pl”分析后所获“PRO_2nd.fa”的图示。
图8为实施例2中步骤5)中Perl脚本“script5.pl”分析后所获“PRO_3rd.fa”的图示。
图9为实施例2中步骤5)中Perl脚本“script5.pl”分析后所获“PRO_RC_1st.fa”的图示。
图10为实施例2中步骤5)中Perl脚本“script5.pl”分析后所获“PRO_RC_2nd.fa”的图示。
图11为实施例2中步骤5)中Perl脚本“script5.pl”分析后所获“PRO_RC_3rd.fa”的图示。
图12为实施例2中步骤5)中Perl脚本“script5.pl”分析后所获“PRO_last.fa”的图示。
图13为实施例2中步骤6)中Perl脚本“script6.pl”分析后所获“res_PRO_seq”文件格式的图示。
图14为实施例2中步骤7)中MUSCLE软件运行之后获得的“PRO_last.out”文件格式的图示。
图15为实施例2中步骤7)中MUSCLE软件运行之后获得的“res_PRO_seq.out”文件格式的图示。
具体实施方式
下面将通过实施例更详细地说明本发明，而这些实施例并不试图限制本发明的保护范围。
下述实施例中所使用的实验方法如无特殊说明，均为常规方法。
下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。
实施例1、批量高效提取基因组基因信息并翻译比对分析序列的方法建立
本发明所提供的批量高效提取基因组基因信息并翻译比对分析序列的方法流程图见图1，具体包括如下步骤：
(1)在Linux或Window系统下操作，建立工作文件夹work_dir，根据某一物种的供试转录本ID(A数据集)，采用Perl脚本“script1.pl”，按照如下步骤从供试转录本ID中提取每行前面的基因ID部分后逐条打印至“res_Gene_ID”文件(B数据集)。
获取“res_Gene_ID”文件的步骤：建立工作目录work_dir，把所述脚本“script1.pl”和供试转录本ID文件XXX1放在工作目录work_dir文件夹下，运行“perlscript1.pl XXX1”命令，得到“res_Gene_ID”文件，记为B数据集。
所述“XXX1”代表A数据集的文件名。所述“res_Gene_ID”文件中放置的是从供试物种转录本ID文件(A数据集)中通过模式匹配批量筛选出的对应的基因ID。该基因ID文件有如下特点，文件每行都包含一个基因ID，且文件中的基因ID与“XXX1”文件中的转录本ID顺序相同。
其中，所述脚本“script1.pl”中关于获得“res_Gene_ID”文件的内容具有如下特点：根据“XXX1”文件中提供的供试转录本ID，利用perl语言的模式匹配语法将对应的基因ID逐个提取出来并打印至结果文件，文件名为“res_Gene_ID”，将该文件放置于work_dir文件夹下。
此外，所述脚本“script1.pl”关于获取“res_Gene_ID”文件的部分是基于如下原理进行编程的：While循环对“XXX1”文件进行逐行处理，对每行进行模式匹配，将Bn开头到“.”符号前的基因ID进行提取并存入变量$gene_id中，将结果写入“res_Gene_ID”文件中，以备使用。
script1.pl(脚本1)
(2)根据所获得的待测物种的供试基因ID文件(B数据集，或者事先准备的基因ID文件)和供试物种基因组信息文件(C数据集)，采用Perl脚本“script2.pl”，按照如下步骤从供试基因ID中提取每个基因的基因组信息后逐条打印至“res_Geneinfo”文件(D数据集)。
获取“res_Geneinfo”文件的步骤：把所述脚本“script2.pl”、获取的供试基因ID文件“res_Gene_ID”以及供试物种基因组信息文件“XXX2”，都放在工作目录work_dir文件夹下，运行perl script2.pl XXX2res_Gene_ID”命令，得到“res_Geneinfo”文件，记为D数据集。
所述“XXX2”代表C数据集的文件名。所述“res_Geneinfo”文件中放置的是从基因组信息文件(C数据集)中通过与供试物种基因ID文件(B数据集)中的基因ID进行基因ID匹配提取批量筛选出的对应的基因组信息内容。该基因组信息文件有如下特点，文件每行都包含基因ID、基因所在染色体号、基因的起始位置、基因的终止位置、基因的长度以及基因所在链的方向等内容。
其中，所述脚本“script2.pl”中关于获得“res_Geneinfo”文件的内容具有如下特点：根据“XXX2”文件中提供的供试基因组信息，利用perl语言模式匹配及两字符串相同的语法将所需基因在基因组上信息逐条提取出来并打印至结果文件，文件名为“res_Geneinfo”，将该文件放置于work_dir文件夹下。
此外，所述脚本“script2.pl”关于获取“res_Geneinfo”文件的部分是基于如下原理进行编程的：将res_Gene_ID文件读入数组@name_can中，打开供试基因组gff文件“XXX2”，while循环逐条处理并分割“XXX2”文件。模式匹配识别“mRNA”标识的行并提取该行的基因ID到变量$id_tmp，for循环遍历数组@name_can的每一行，当变量$id_tmp与数组某行基因ID相同时，计算该基因的长度并存入到变量$genelen中，把基因ID、所在染色体号、基因的起始位点、终止位点、基因长度以及正反义链等信息，逐行打印至文件“res_Geneinfo”中，以备使用。
script2.pl(脚本2)
(3)为StrawberryPerl软件安装Bioperl模块，根据待测物种的供试转录本ID(A数据集)和cds文件(E数据集)，采用Perl脚本“script3.pl”，按照如下步骤利用供试转录本ID匹配cds文件中的ID，进而从cds文件中提取对应的序列后逐条打印至“res_CDS_seq”文件(G数据集)。
获取“res_CDS_seq”文件的步骤：把所述脚本“script1.pl”、供试转录本ID文件XXX1以及供试物种基因组cds文件XXX4都放在工作目录work_dir文件夹下，运行“perlscript3.pl XXX1”命令，得到“res_CDS_seq”文件，记为G数据集。
所述“XXX1”代表A数据集的文件名，“XXX4”代表E数据集的文件名。所述“res_CDS_seq”文件中放置的是从物种基因组cds文件(E数据集)中通过对转录本ID进行匹配后批量筛选出的序列文件。该基因ID文件有如下特点，文件中包含转录本ID及其对应的cds序列，该文件为fasta格式，且文件中的转录本ID顺序与“XXX1”文件中的转录本ID顺序相同。
其中，所述脚本“script3.pl”中关于获得“res_CDS_seq”文件的内容具有如下特点：根据“XXX1”文件中提供的供试转录本ID，利用perl语言的字符串匹配语法去对应的供试物种cds序列文件“XXX4”中将对应的cds序列逐个提取出来并打印至结果文件，文件名为“res_CDS_seq”，将该文件放置于work_dir文件夹下。
此外，所述脚本“script3.pl”关于获取“res_CDS_seq”文件的部分是基于如下原理进行编程的：利用Bio::SeqIO模块和while循环将供试基因组cds文件“XXX4”逐条读入哈希％hash中，打开供试转录本ID文件“XXX1”，While循环对“XXX1”文件进行逐行处理，if判别如果存在以“XXX1”文件中某行的转录本ID为key值的哈希value值$hash{$line}，则去掉$hash{$line}后的最后一个“*”号，并将转录本ID以及对应的哈希value值即cds序列，以fasta的格式逐条打印至同一结果文件中，文件名为“res_CDS_seq”，else条件如果不存在上述哈希value值，则在屏幕输出该转录本ID未找到的语句。把该结果文件“res_CDS_seq”置于当前工作目录work_dir文件夹中，以备使用。
script3.pl(脚本3)
(4)根据所获得的待测物种的供试基因信息文件(D数据集)和供试物种基因组染色体文件(F数据集)，采用Perl脚本“script4.pl”，按照如下步骤从供试基因组染色体序列中提取每个基因的碱基序列后逐条打印至“res_Gene_seq”文件(H数据集)。
获取“res_Gene_seq”文件的步骤：把所述脚本“script4.pl”、获取的供试基因信息文件“res_Geneinfo”以及供试物种基因组染色体文件“XXX5”，都放在工作目录work_dir文件夹下，运行“perl script 4.pl res_Geneinfo”命令，得到“res_Gene_seq”文件，记为H数据集。
所述“XXX5”代表F数据集的文件名。所述“res_Gene_seq”文件中放置的是从供试物种基因组染色体文件(F数据集)中通过所获基因组信息文件(D数据集)中的染色体信息、基因ID、位置以及长度信息提取批量筛选出的对应基因的碱基序列。该基因序列文件有如下特点，文件中包含基因ID及其对应的碱基序列，该文件为fasta格式，且文件中的基因ID顺序与“res_Geneinfo”文件中的基因ID顺序相同。
其中，所述脚本“script4.pl”中关于获得“res_Gene_seq”文件的内容具有如下特点：根据获取的“res_Geneinfo”文件中提供的基因组信息，利用perl语言语法去对应的供试物种染色体文件“XXX5”中将对应的基因序列逐个提取出来并打印至结果文件，文件名为“res_Gene_seq”，将该文件放置于work_dir文件夹下。
此外，所述脚本“script4.pl”关于获取“res_Gene_seq”文件的部分是基于如下原理进行编程的：利用Bio::SeqIO模块和while循环将供试基因组染色体文件“XXX5”逐条读入哈希％hash中，打开文件“res_Geneinfo”，While循环对其进行逐行处理，next if语句去掉以字母“G”开头的行后逐行分割文件，通过substr函数依靠文件中基因的起始、终止位置变量和基因长度变量$row[1]、$row[2]和$row[4]，提取位于染色体$hash{$row[1]}上的基因序列，并存入变量$seq_tmp中。If判别如果该基因的方向是反义链“-”，则将该序列的反向互补序列求出，存于变量$seq_tmp中。最后将所有结果以基因ID对应序列的fasta文件格式打印至“res_Gene_seq”文件中，以备使用。
script4.pl(脚本4)
(5)根据获取的待测物种cds序列(G数据集)，采用Perl脚本“script5.pl”，按照如下步骤对供试待测物种cds序列进行逐行翻译后将最终选取的的蛋白质文件名写入“PRO_last.fa”文件(O数据集)。
获取“PRO_last.fa”文件的步骤：把所述脚本“script5.pl”和供试物种cds序列文件“res_CDS_seq”放在工作目录work_dir文件夹下，运行“perl script 5.pl”命令，得到“PRO_1st.fa”、“PRO_2nd.fa”、“PRO_3rd.fa”、“PRO_RC_1st.fa”、“PRO_RC_2nd.fa”、“PRO_RC_3rd.fa”和“PRO_last.fa”7个文件，分别记为记为I、J、K、L、M、N和O数据集。
所述“res_CDS_seq”是获取的供试基因的cds序列文件。所述“PRO_1st.fa”、“PRO_2nd.fa”、“PRO_3rd.fa”、“PRO_RC_1st.fa”、“PRO_RC_2nd.fa”和“PRO_RC_3rd.fa”文件中放置的是分别按着不同的起始位置以及正反义链方向对供试基因cds序列文件(G数据集)进行逐条批量筛翻译后的蛋白质文件。该蛋白序列文件有如下特点，文件中包含基因ID及其对应的蛋白质序列，该文件为fasta格式。“PRO_last.fa”文件中放置的是最后筛选出的用于后续多序列比对分析的蛋白质序列文件，该蛋白序列文件有如下特点，文件中包含基因ID及其对应的蛋白质序列，该文件为fasta格式。
其中，所述脚本“script5.pl”中关于获得“PRO_1st.fa”、“PRO_2nd.fa”、“PRO_3rd.fa”、“PRO_RC_1st.fa”、“PRO_RC_2nd.fa”和“PRO_RC_3rd.fa”等文件的内容具有如下特点：根据“res_CDS_seq”文件中提供的供试基因cds序列，利用perl语言事先创建包含翻译密码子简表的哈希，后根据起始位置不同以及链的正反方向等因素进行翻译，进而得到六种不同的蛋白质文件，文件名分别为“PRO_1st.fa”、“PRO_2nd.fa”、“PRO_3rd.fa”、“PRO_RC_1st.fa”、“PRO_RC_2nd.fa”和“PRO_RC_3rd.fa”；关于获得“PRO_last.fa”文件的内容具有如下特点：根据得出的六个不同的蛋白质文件，利用perl语言选出其中最大的文件，即为翻译最全面正确的文件，并将该文件信息逐行提取出来并打印至结果文件，文件名为“PRO_last.fa”，将该文件放置于work_dir文件夹下。
此外，所述脚本“script5.pl”关于获取“PRO_1st.fa”、“PRO_2nd.fa”、“PRO_3rd.fa”、“PRO_RC_1st.fa”、“PRO_RC_2nd.fa”、“PRO_RC_3rd.fa”和“PRO_last.fa”文件的部分是基于如下原理进行编程的：首先把20种氨基酸的64种密码子在程序里面写完整并存入哈希％genetic_code中。打开待翻译DNA序列的fasta文件，依靠Bio::SeqIO模块接收供试DNA序列。
然后利用while循环逐条读取输入文件DNA序列，利用uc函数将序列字母转换为大写，利用reverse函数及正则表达式tr///获取读取DNA序列的反向互补序列，length函数计算序列长度。利用存储密码子简表的哈希％genetic_code，分别从读取DNA序列起始位置的第一、二和三位密码子开始向后进行翻译(以三个相连密码子为翻译单位)，将翻译后的蛋白质序列及其ID以fasta文件格式写入结果文件PRO_1st.fa、PRO_2nd.fa和PRO_3rd.fa中；同时分别从计算所得DNA序列反向互补序列起始位置的第一、二和三位密码子开始向后进行翻译(以三个相连密码子为翻译单位)，将翻译后的蛋白质序列及其ID以fasta文件格式写入结果文件PRO_RC_1st.fa、PRO_RC_2nd.fa和PRO_RC_3rd.fa中，同时将所得6个结果文件置于当前工作目录work_dir文件夹中，关闭所有相关文件。
第三，stat函数分别取6个结果文件的文件大小，存入数组@array_size中，并分别以文件大小和文件名为哈希％hash_size的key值和value值；将数组中元素按着从大到小的顺序排序后存入新数组@array_sort中，然后筛选出@array_sort中最大的元素$array_sort[0]，以及以该最大元素为key值对应的哈希value值$hash_size{$array_sort[0]}存入变量$file_biggest中，最后将$file_biggest文件打开，并利用Bio::SeqIO模块嵌套while循环，将该文件内容逐行打印至结果文件“PRO_last.fa”中，以备使用。
script5.pl(脚本5)
/>
/>
/>
/>
/>
(6)根据待测物种的供试转录本ID(A数据集)和蛋白质文件(P数据集)，采用Perl脚本“script6.pl”，按照如下步骤利用供试转录本ID匹配蛋白质文件中的ID，进而从蛋白质文件中提取对应的序列后逐条打印至“res_PRO_seq”文件(Q数据集)。
获取“res_PRO_seq”文件的步骤：把所述脚本“script6.pl”、供试转录本ID文件XXX1以及供试物种基因组蛋白质文件XXX6都放在工作目录work_dir文件夹下，运行“perlscript6.pl XXX1”命令，得到“res_PRO_seq”文件，记为Q数据集。
所述“XXX1”代表A数据集的文件名，“XXX6”代表P数据集的文件名。所述“res_PRO_seq”文件中放置的是从物种基因组蛋白质文件(P数据集)中通过对转录本ID进行匹配后批量筛选出的序列文件。该基因蛋白文件有如下特点，文件中包含转录本ID(蛋白质ID)及其对应的蛋白质序列，该文件为fasta格式，且文件中的转录本ID顺序与“XXX1”文件中的转录本ID顺序相同。
其中，所述脚本“script6.pl”中关于获得“res_PRO_seq”文件的内容具有如下特点：根据“XXX1”文件中提供的供试转录本ID，利用perl语言的字符串匹配语法去对应的供试物种蛋白质序列文件“XXX6”中将对应的蛋白质序列逐个提取出来并打印至结果文件，文件名为“res_PRO_seq”，将该文件放置于work_dir文件夹下。
此外，所述脚本“script6.pl”关于获取“res_PRO_seq”文件的部分是基于如下原理进行编程的：利用Bio::SeqIO模块和while循环将供试基因组蛋白质文件“XXX6”逐条读入哈希％hash中，打开供试转录本ID文件“XXX1”，While循环对“XXX1”文件进行逐行处理，if判别如果存在以“XXX1”文件中某行的转录本ID为key值的哈希value值$hash{$line}，则去掉$hash{$line}后的最后一个“*”号，并将转录本ID以及对应的哈希value值即蛋白质序列，以fasta的格式逐条打印至同一结果文件中，文件名为“res_PRO_seq”，else条件如果不存在上述哈希value值，则在屏幕输出该转录本ID未找到的语句。把该结果文件“res_PRO_seq”置于当前工作目录work_dir文件夹中，以备使用。
script6.pl(脚本6)
(7)采用蛋白质多序列比对的MUSCLE程序，对供试物种中的蛋白质进行多序列比对分析，在Linux系统下操作，采用默认参数设置进行分析；按照如下步骤获取完整的多序列比对结果(R数据集)；
获取完整的蛋白质多序列比对文件步骤：打开待测蛋白质文件所在文件夹，运行“muscle-in PRO_last.fa–out PRO_last.out”命令(或者如果存在上述(6)，则运行“muscle-in res_PRO_seq–outres_PRO_seq.out”命令)，其中，“PRO_last.fa”为根据cds序列翻译后的蛋白质序列(如果能够直接下载到基因组蛋白质文件，则“res_PRO_seq”为依据供试转录本ID提取的蛋白质序列文件)。此步骤将得到“PRO_last.out”文件(或者“res_PRO_seq.out”文件)。
所述“PRO_last.out”文件(或者“res_PRO_seq.out”文件)，文件中的数据为供试蛋白质多序列比对后的结果文件，记为R数据集。
实施例2、利用实施例1建立的方法批量高效提取中双11油菜基因组基因信息并完成翻译比对序列分析
进入中国农业科学院油料作物研究所的油料作物基因组数据库(http://ocri- genomics.org/Brassia_napus_genome_ZS11/)下载中双11油菜(Brassica napus)的基因组序列(19条染色体，976Mb)。在Windows系统或本地Linux运算服务器，进行高效提取中双11油菜基因组基因信息并完成翻译序列比对计算。计算过程中，所涉及的常用程序名称、运行环境及地址如表1所示。计算方法具体操作步骤如下：
1)参照实施例1的步骤(1)进行。
利用所获中双11油菜271个供试转录本ID文件(A数据集)，在Linux或Window系统下建立工作文件夹work_dir，采用Perl脚本“script1.pl”，根据Perl语言模式匹配的语法对271个转录本ID文件进行逐行提取基因ID，并将271个基因ID结果依次写入如图2所示结构格式的“res_Gene_ID”文件(B数据集)，且A和B数据集皆放置于work_dir文件夹下。
2)参照实施例1的步骤(2)进行。
利用所获中双11油菜271个供试基因ID文件(B数据集)和670060行的中双11油菜基因组信息gff文件(C数据集)，采用Perl脚本“script2.pl”，通过271个供试基因的ID匹配从中双11油菜基因组信息文件(B数据集)批量提取筛选出的对应271个基因的基因组信息，并将该271条基因组信息依次写入如图3所示结构格式(基因ID、基因所在染色体号、基因的起始位置、基因的终止位置、基因的长度以及基因所在链的方向)的“res_Geneinfo”文件(D数据集)，且B、C、D数据集皆放置于work_dir文件夹下。
3)参照实施例1的步骤(3)进行。
利用所获中双11油菜271个供试转录本ID文件(A数据集)和101942条双11油菜基因组cds序列文件(E数据集)，为StrawberryPerl软件安装Bioperl模块，采用Perl脚本“script3.pl”，通过271个供试转录本ID匹配中双11油菜基因组cds序列文件(E数据集)中的转录本ID信息，批量提取筛选出对应的271条cds序列，并将其依次写入如图4所示结构格式(转录本ID对应cds序列的fasta序列文件格式)的“res_CDS_seq”文件(G数据集)，且A、E、G数据集皆放置于work_dir文件夹下。
4)参照实施例1的步骤(4)进行。
利用中双11油菜基因组19条染色体文件(F数据集)和所获271条基因组信息文件(D数据集)，采用Perl脚本“script4.pl”，通过所获271条基因组信息文件(D数据集)中的染色体信息、基因ID、位置以及长度信息，在基因组染色体文件(F数据集)中提取批量筛选出对应基因的碱基序列，并将其依次写入如图5所示结构格式(基因ID对应基因碱基序列的fasta序列文件格式)的“res_Gene_seq”文件(H数据集)，且D、F、H数据集皆放置于work_dir文件夹下。
5)参照实施例1的步骤(5)进行。
利用所获中双11油菜271条基因cds序列文件(G数据集)，采用Perl脚本“script5.pl”，按着不同的起始位置以及正反义链方向对271条基因cds序列文件(G数据集)中碱基序列进行逐条批量翻译蛋白质。并将271条翻译结果分别依次写入如图6-11所示结构格式(转录本ID对应蛋白质序列的fasta序列文件格式)的“PRO_1st.fa”、“PRO_2nd.fa”、“PRO_3rd.fa”、“PRO_RC_1st.fa”、“PRO_RC_2nd.fa”和“PRO_RC_3rd.fa”六个文件(I、J、K、L、M和N数据集)。然后从中筛选出唯一正确可用于进行后续分析的结果文件将其写入如图12所示结构格式(转录本ID对应蛋白质序列的fasta序列文件格式)的“PRO_last.fa”文件(G数据集)，且G、I、J、K、L、M、N和O数据集皆放置于work_dir文件夹下。
6)参照实施例1的步骤(6)进行。
利用所获中双11油菜271个供试转录本ID文件(A数据集)和101942条双11油菜基因组蛋白质序列文件(P数据集)，采用Perl脚本“script6.pl”，通过271个供试转录本ID匹配中双11油菜基因组蛋白质序列文件(P数据集)中的转录本ID信息，批量提取筛选出对应的271条蛋白质序列，并将其依次写入如图13所示结构格式(转录本ID对应蛋白质序列的fasta序列文件格式)的“res_PRO_seq”文件(Q数据集)，且A、P、Q数据集皆放置于work_dir文件夹下。
7)参照实施例1的步骤(7)进行。
在Linux系统下操作，在将步骤5)所获取的“PRO_last.fa”文件(O数据集)所在的work_dir文件夹内，运行“muscle-inPRO_last.fa–outPRO_last.out”命令，得到如图14所示格式的“PRO_last.out”文件。所述“PRO_last.out”文件为MUSCLE软件计算后的结果文件，记为R数据集；或者如果能够直接得到供试基因组蛋白序列的话，在将步骤6)所获取的“res_PRO_seq”文件(Q数据集)所在的work_dir文件夹内，运行“muscle-in res_PRO_seq–outres_PRO_seq.out”命令，得到如图15所示格式的“res_PRO_seq.out”文件。所述“res_PRO_seq.out”文件为MUSCLE软件计算后的结果文件，记为R数据集。
表1说明书中的常用软件
本发明的发明人从最终所得的271条基因序列(数据集H)中随机选择10条序列，手动找出它们在染色体文件(数据集F)上的位置、长度、正反义链然后与基因组信息文件结果(数据集D)进行人工校对，发现位置、长度以及链方向等信息都准确无误；发明人从最终所得的271条翻译的蛋白序列(数据集O)文件中随机选择10条序列，手动找出它们在271条cds文件(数据集G)中的ID，并对相同ID下的蛋白质序列以及cds序列进行人工校对，发现蛋白质和cds序列信息都准确无误；发明人从最终所获取的271条蛋白质文件(数据集Q)中随机选择10条序列，手动找出它们在供试基因组蛋白质文件(数据集P)中的ID，并对相同ID下的蛋白质序列进行人工校对，发现蛋白质序列信息准确无误，从而证实了以上本发明方法的准确性。

Claims (3)

1.一种批量提取基因组基因信息并翻译比对分析序列的方法，其特征在于，将某一物种的转录本ID或者基因ID，依据供试基因组cds文件、蛋白质文件、gff文件和染色体fasta文件信息，通过script1.pl、script2.pl、script3.pl、script 4.pl、script 5.pl、script6.pl 6个perl脚本程序，实现目标转录本或基因在基因组上的位置、长度、正反义链结构信息的提取，并在染色体fasta文件上提取该转录本或基因的cds或基因序列，在基因组蛋白文件上提取该转关闭所有相关文录本的蛋白质序列；最后对所需cds序列进行翻译，或者直接用所获的蛋白质序列，调用Linux系统程序完成蛋白质的多序列比对工作；
包括如下步骤：
(1)建立工作文件夹work_dir，将某一物种的转录本ID文件记为数据集A，所述数据集A的文件名为“XXX1”，运行“perl script1.pl XXX1”命令，在当前工作文件夹work_dir下得到“res_Gene_ID”文件；所述“XXX1”在运行“perl script1.pl XXX1”程序时已置于包含脚本“script1.pl”的当前工作文件夹work_dir内，关闭所有相关文件；所述“res_Gene_ID”文件为该物种转录本ID对应的基因ID文件，记为数据集B，命名为“XXX3”；
如果上述步骤直接提供的是某一物种基因ID，则将该基因ID文件记为数据集B，命名为“XXX3”；
(2)将该物种基因组gff文件记为C数据集，所述C数据集的文件名为“XXX2”，运行“perlscript2.pl XXX2 XXX3”命令，在当前工作文件夹work_dir下得到“res_Geneinfo”文件；
所述“res_Geneinfo”文件为根据该物种基因ID文件提取的基因组信息文件，记为数据集D；所述“XXX2”、“XXX3”在运行“perl script2.pl XXX2 XXX3”程序时已置于包含脚本“script2.pl”的当前工作文件夹work_dir内，关闭所有相关文件；
(3)为Strawberry Perl软件安装Bioperl模块，将该物种基因组cds的fasta格式文件记为数据集E，所述数据集E的文件名为“XXX4”，运行“perl script3.pl XXX1”命令，在当前工作文件夹work_dir下得到“res_CDS_seq”文件；
所述“res_CDS_seq”文件为根据该物种转录本ID文件提取的基因cds序列fasta文件，记为数据集G；所述“XXX4”在运行“perl script3.pl XXX1”程序时已置于包含脚本“script3.pl”的当前工作文件夹work_dir内，关闭所有相关文件；
(4)将该物种基因组染色体的fasta格式文件记为数据集F，所述数据集F的文件名为“XXX5”，运行“perl script 4.pl res_Geneinfo”命令，在当前工作文件夹work_dir下得到“res_Gene_seq”文件；
所述“res_Gene_seq”文件为根据该物种基因ID文件从该物种基因组染色体文件中提取的基因序列fasta文件，记为数据集H；所述“XXX5”在运行“perl script 4.pl res_Geneinfo”程序时已置于包含脚本“script 4.pl”的当前工作文件夹work_dir内，关闭所有相关文件；
(5)在当前工作文件夹work_dir内运行“perl script 5.pl”命令，得到“PRO_1st.fa”、“PRO_2nd.fa”、“PRO_3rd.fa”、“PRO_RC_1st.fa”、“PRO_RC_2nd.fa”、“PRO_RC_3rd.fa”和“PRO_last.fa”7个文件；
所述“PRO_1st.fa”、“PRO_2nd.fa”、“PRO_3rd.fa”、“PRO_RC_1st.fa”、“PRO_RC_2nd.fa”和“PRO_RC_3rd.fa”6个文件为根据该物种基因ID文件提取的基因序列或者转录本cds序列翻译后的蛋白质fasta文件，分别记为数据集I、J、K、L、M和N；所述“PRO_last.fa”文件为筛选出的用于后续多序列比对计算的蛋白质序列文件，记为数据集O；所述“res_CDS_seq”文件在运行“perl script 5.pl”程序时已置于包含脚本“perl script 5.pl”的当前工作文件夹work_dir内，关闭所有相关文件；
(6)如果通过下载获得该物种基因组蛋白质的fasta格式文件，则将其记为P数据集，所述P数据集的文件名为“XXX6”，运行“perl script6.pl XXX1”命令，在当前工作文件夹work_dir下得到“res_PRO_seq”文件；
所述“res_PRO_seq”文件为根据该物种转录本ID文件提取的基因蛋白质序列fasta文件，记为数据集Q；所述“XXX6”在运行“perl script6.pl XXX1”程序时已置于包含脚本“script6.pl”的当前工作文件夹work_dir内，关闭所有相关文件；
(7)在当前工作文件夹work_dir内运行“muscle-in PRO_last.fa–out PRO_last.out”命令，如果存在上述步骤(6)，则运行“muscle-inres_PRO_seq–out res_PRO_seq.out”命令，在当前工作文件夹中得到多重序列比对的结果文件；
所述“PRO_last.out”和“res_PRO_seq.out”文件为MUSCLE软件计算后的输出文件，记为数据集R；且在运行“muscle-in PRO_last.fa–out PRO_last.out”命令或者“muscle-inres_PRO_seq–out res_PRO_seq.out”命令后所产生的结果文件在当前工作文件夹work_dir中，关闭所有相关文件；
步骤(1)中，所述脚本“script1.pl”中关于获得“res_Gene_ID”文件是基于如下方法进行编程的：
While循环对“XXX1”文件进行逐行处理，对每行进行模式匹配，将Bn开头到“.”符号前的基因ID进行提取并存入变量$gene_id中，将结果打印至同一文件中，文件名即为“res_Gene_ID”，同时把该文件置于当前工作目录work_dir文件夹中，关闭所有相关文件;
步骤(2)中，所述脚本“script2.pl”中关于获得“res_Geneinfo”文件是基于如下方法进行编程的：
将res_Gene_ID文件读入数组@name_can中，打开该物种基因组gff文件“XXX2”，while循环逐条处理并分割“XXX2”文件；模式匹配识别“mRNA”标识的行并提取该行的基因ID到变量$id_tmp，for循环遍历数组@name_can的每一行，当变量$id_tmp与数组某行基因ID相同时，计算该基因的长度并存入到变量$genelen中，把基因ID、所在染色体号、基因的起始位点、终止位点、基因长度以及正反义链信息，逐行打印至同一文件，文件名为“res_Geneinfo”，同时把该文件置于当前工作目录work_dir文件夹中，关闭所有相关文件；
步骤(3)中，所述脚本“script3.pl”中关于获得“res_CDS_seq”文件是基于如下方法进行编程的：
利用Bio::SeqIO模块和while循环将供试基因组cds文件“XXX4”逐条读入哈希％hash中，打开供试转录本ID文件“XXX1”，While循环对“XXX1”文件进行逐行处理，if判别如果存在以“XXX1”文件中某行的转录本ID为key值的哈希value值$hash{$line}，则去掉$hash{$line}后的最后一个“*”号，并将转录本ID以及对应的哈希value值即cds序列，以fasta的格式逐条打印至同一结果文件中，文件名为“res_CDS_seq”，else条件如果不存在上述哈希value值，则在屏幕输出转录本ID未找到，把该结果文件“res_CDS_seq”置于当前工作目录work_dir文件夹中，关闭所有相关文件；
步骤(4)中，所述脚本“script4.pl”中关于获得“res_Gene_seq”文件是基于如下方法进行编程的：
利用Bio::SeqIO模块和while循环将待测物种基因组染色体文件“XXX5”逐条读入哈希％hash中，打开文件“res_Geneinfo”，While循环对其进行逐行处理，next if语句去掉以字母“G”开头的行后逐行分割文件，通过substr函数依靠文件中基因的起始、终止位置变量和基因长度变量$row[1]、$row[2]和$row[4]，提取位于染色体$hash{$row[1]}上的基因序列，并存入变量$seq_tmp中；If判别如果该基因的方向是反义链“-”，则将该序列的反向互补序列求出，存于变量$seq_tmp中；最后将所有结果以基因ID对应序列的fasta文件格式打印至同一文件中，文件名为“res_Gene_seq”，同时把该文件置于当前工作目录work_dir文件夹中，关闭所有相关文件；
步骤(5)中，所述脚本“script5.pl”中关于获得“PRO_last.fa”文件是基于如下方法进行编程的：
首先把20种氨基酸的64种密码子在程序里面写完整并存入哈希％genetic_code中，打开待翻译DNA序列的fasta文件，依靠Bio::SeqIO模块接收供试DNA序列；
然后利用while循环逐条读取输入文件DNA序列，利用uc函数将序列字母转换为大写，利用reverse函数及正则表达式tr///获取读取DNA序列的反向互补序列，length函数计算序列长度；利用存储密码子简表的哈希％genetic_code，分别从读取DNA序列起始位置的第一、二和三位密码子开始向后进行翻译，以三个相连密码子为翻译单位，将翻译后的蛋白质序列及其ID以fasta文件格式写入结果文件PRO_1st.fa、PRO_2nd.fa和PRO_3rd.fa中；同时分别从计算所得DNA序列反向互补序列起始位置的第一、二和三位密码子开始向后进行翻译，以三个相连密码子为翻译单位，将翻译后的蛋白质序列及其ID以fasta文件格式写入结果文件PRO_RC_1st.fa、PRO_RC_2nd.fa和PRO_RC_3rd.fa中，同时将所得6个结果文件置于当前工作目录work_dir文件夹中，关闭所有相关文件；
第三，stat函数分别获取6个结果文件的文件大小，存入数组@array_size中，并分别以文件大小和文件名为哈希％hash_size的key值和value值；将数组中元素按着从大到小的顺序排序后存入新数组@array_sort中，然后筛选出@array_sort中最大的元素$array_sort[0]，以及以该最大元素为key值对应的哈希value值$hash_size{$array_sort[0]}存入变量$file_biggest中，最后将$file_biggest文件打开，并利用Bio::SeqIO模块嵌套while循环，将该文件内容逐行打印至结果文件“PRO_last.fa”中，同时把该文件置于当前工作目录work_dir文件夹中，关闭所有相关文件；
步骤(6)中，所述脚本“script6.pl”中关于获得“res_PRO_seq”文件是基于如下方法进行编程的：
利用Bio::SeqIO模块和while循环将供试基因组蛋白质文件“XXX6”逐条读入哈希％hash中，打开供试转录本ID文件“XXX1”，While循环对“XXX1”文件进行逐行处理，if判别如果存在以“XXX1”文件中某行的转录本ID为key值的哈希value值$hash{$line}，则去掉$hash{$line}后的最后一个“*”号，并将转录本ID以及对应的哈希value值即蛋白质序列，以fasta的格式逐条打印至同一结果文件中，文件名为“res_PRO_seq”，else条件如果不存在上述哈希value值，则在屏幕输出该转录本ID未找到，把该结果文件“res_PRO_seq”置于当前工作目录work_dir文件夹中，关闭所有相关文件。
2.根据权利要求1中任一所述的方法，其特征在于：步骤(1)中，所述某一物种为任意完成全基因组测序的物种。
3.根据权利要求1中任一所述的方法，其特征在于：所述某一物种基因组序列通过下载已公开的全基因组注释gff文件、cds序列、染色体序列和蛋白质序列获得，或通过全基因组测序获得相关文件。
CN201910684539.6A 2019-07-26 2019-07-26 一种批量提取基因组基因信息并翻译比对分析序列的方法 Active CN110534157B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910684539.6A CN110534157B (zh) 2019-07-26 2019-07-26 一种批量提取基因组基因信息并翻译比对分析序列的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910684539.6A CN110534157B (zh) 2019-07-26 2019-07-26 一种批量提取基因组基因信息并翻译比对分析序列的方法

Publications (2)

Publication Number Publication Date
CN110534157A CN110534157A (zh) 2019-12-03
CN110534157B true CN110534157B (zh) 2023-07-25

Family

ID=68661941

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910684539.6A Active CN110534157B (zh) 2019-07-26 2019-07-26 一种批量提取基因组基因信息并翻译比对分析序列的方法

Country Status (1)

Country Link
CN (1) CN110534157B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111816254A (zh) * 2020-06-01 2020-10-23 上海派森诺生物科技股份有限公司 一种基于perl语言快速批量去除载体序列的方法
CN112712850A (zh) * 2020-12-29 2021-04-27 中南大学 一种可应用于传染病病原体测序读段映射的种子序列定位方法
CN113066530B (zh) * 2021-03-31 2024-05-10 江苏省农业科学院 一种批量合并eQTL分析结果中存在连锁不平衡SNP的方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103824000A (zh) * 2014-02-24 2014-05-28 江苏省农业科学院 一种批量检测植物基因组ltr-反转座子的方法
CN105274092A (zh) * 2015-11-30 2016-01-27 中国人民解放军军事医学科学院卫生学环境医学研究所 一种特异性等温寡核苷酸探针的批量获取方法
CN105426700A (zh) * 2015-12-18 2016-03-23 江苏省农业科学院 一种批量计算基因组直系同源基因进化速率的方法
CN107091929A (zh) * 2016-02-25 2017-08-25 安徽省农业科学院水稻研究所 一种启动子批量捕获方法
CN107122624A (zh) * 2017-05-01 2017-09-01 杨永臣 人类基因突变的hgvs名称生成及分析系统的实现方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040002818A1 (en) * 2001-12-21 2004-01-01 Affymetrix, Inc. Method, system and computer software for providing microarray probe data

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103824000A (zh) * 2014-02-24 2014-05-28 江苏省农业科学院 一种批量检测植物基因组ltr-反转座子的方法
CN105274092A (zh) * 2015-11-30 2016-01-27 中国人民解放军军事医学科学院卫生学环境医学研究所 一种特异性等温寡核苷酸探针的批量获取方法
CN105426700A (zh) * 2015-12-18 2016-03-23 江苏省农业科学院 一种批量计算基因组直系同源基因进化速率的方法
CN107091929A (zh) * 2016-02-25 2017-08-25 安徽省农业科学院水稻研究所 一种启动子批量捕获方法
CN107122624A (zh) * 2017-05-01 2017-09-01 杨永臣 人类基因突变的hgvs名称生成及分析系统的实现方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
张大勇等.《基于Perl 脚本的大豆核苷酸序列高通量提取》.《南京农业大学学报》.2012,第第35卷卷(第第35卷期),第7-13页. *
郭景康等.《基因和蛋白质的批量注释系统UBROAD》.《上海大学学报（ 自然科学版）》.2007,第第13卷卷(第第13卷期),第99-104,110页. *

Also Published As

Publication number Publication date
CN110534157A (zh) 2019-12-03

Similar Documents

Publication Publication Date Title
CN110534157B (zh) 一种批量提取基因组基因信息并翻译比对分析序列的方法
Gremme et al. Engineering a software tool for gene structure prediction in higher organisms
CN106599614B (zh) 一种高通量测序数据处理及分析流程控制方法及系统
CN111192630B (zh) 一种宏基因组数据挖掘方法
CN114420212B (zh) 一种大肠杆菌菌株鉴定方法和系统
CN113299344A (zh) 基因测序分析方法、装置、存储介质和计算机设备
CN112786102A (zh) 一种基于宏基因组学分析精准识别水体中未知微生物群落的方法
CN105426700B (zh) 一种批量计算基因组直系同源基因进化速率的方法
CN111180013B (zh) 检测血液病融合基因的装置
CN109243531B (zh) 一种批量计算近缘物种间基因组编码区snp位点的方法
Contreras-Moreira et al. RSAT:: Plants: motif discovery within clusters of upstream sequences in plant genomes
CN110570901B (zh) 一种基于测序数据进行ssr分型的方法及系统
CN108182348A (zh) 基于种子序列信息的dna甲基化数据检测方法及其装置
Sun et al. Designing patterns for profile HMM search
CN115295084A (zh) 一种肿瘤新抗原免疫组库数据可视化分析方法和系统
EP3693970A1 (en) Biological sequence information handling
CN111492436A (zh) 使用k聚体在没有比对的情况下进行测序数据的快速质量控制
CN111696629A (zh) 一种rna测序数据的基因表达量计算方法
Gao et al. AUSPP: a universal short-read pre-processing package
CN109493918A (zh) 一种生物数据管理及系统发育分析流程化方法
CN109741788A (zh) 一种snp位点分析方法及系统
Thallinger Comparison of ddRAD Analysis Pipelines
Conklin et al. Mining of assembled expressed sequence tag (EST) data for protein families: Application to the G protein-coupled receptor superfamily
Royaux et al. Hands-on: Hands-on: Metabarcoding/eDNA through Obitools
Looso et al. Data mining in newt-omics, the repository for omics data from the newt

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