[发明专利]一种批量提取基因组基因信息并翻译比对分析序列的方法

权利要求书

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文件夹中，关闭所有相关文件。