perl计算基因在基因组上的距离

介绍

根据输入注释文件及目的基因文件，找出目的基因在基因组上的距离并输出到文件内

一、文件格式

基因注释文件：为细菌全基因组利用prokka注释文件，格式如下
gene和CDS前为基因的起始和结束位置，locus_tag为基因编号，product为注释结果，以制表符分隔。

目的基因文件xxx.txt：

以制表符分隔，第一行为描述，第二行为全基因组长度bp，第三行及之后为所要计算距离的基因。

二、代码

这里输入注释文件为ZJSH63.tbl，目的基因文件为genes.txt
代码distance.pl如下：

#! /usr/bin/perl -w
use strict;
my @list = @ARGV;
#运行脚本示例：perl distance.pl ZJSH63.tbl  genes.txt
#ZJSH63.tbl文件为prokka注释文件
#genes.txt为所要计算距离的基因编号文件
open(Genes,"$ARGV[1]") or die ("No genes file finded: $!");
my @gene;
my %location;
my @temp1;
my @temp2;
my @temp3;
my $num =0;
my $temp;
while(<Genes>){chomp;#geen.txt第0列为信息，第2列为菌株、基因组长度和基因编号@temp1 = split(/\t/,$_);#@gene 0位为菌名，1位为基因组长度bp，后续为基因编号$temp = $temp1[1];chomp($temp);if($temp =~/\r$/){chop($temp);}#print "$temp\t";$gene[$num] = $temp;#print "$gene[$num]\t";$num++;
}
close (Genes) or die ("can not close genes file: $!");print "number of genes: $num\n";
my $name;
my $left = "left";
my $right = "right";
my $l;
my $r;
foreach $name (@gene[2..$num-1]){chomp($name);#print "$name\t";$l = $left.$name;$r = $right.$name;$location{$l} = 0;$location{$r} = 0;#print "$l\t";#print "$location{$l}\t","$r\t","$location{$r}\n";
}
#my @key = keys(%location);
#print "@key\n";
my $judge =0;
my $word1 = "locus_tag";
my $word2 = "CDS";
open(Annotate,"$ARGV[0]") or die ("No annotate file finded: $!");
open(OUT,">result.txt") or die ("can not write result: $!");
print OUT "Genes\tleft\tright\n";
while(<Annotate>){chomp;if(/$word1/){if($judge == 2){$judge =0; next;}@temp3 = split(/\t/,$_);$temp = $temp3[-1];chomp($temp);if($temp=~/\r$/){chop($temp);}if(grep /$temp/, @gene){$l = $left.$temp;$r = $right.$temp;print "find $temp\t";print OUT "$temp\t";#print "$l\t";$judge =1;#$temp_gene = $temp3[-1];next;}}if($judge ==1){if(/$word2/){@temp2 = split(/\t/, $_);my $te1 = $temp2[0];chomp($te1);if($te1=~/\r$/){chop($te1);}my $te2 = $temp2[1];chomp($te2);if($te2=~/\r$/){chop($te2);}if($te1 > $ te2){$location{$l}=$te2;$location{$r}=$te1;}else{$location{$l}=$te1;$location{$r}=$te2;}#print "$l\t";print "$location{$l}\t";print OUT "$location{$l}\t";print "$location{$r}\n";print OUT "$location{$r}\n";$judge =2;#$temp_gene ="";next;}}
}
close (Annotate) or die ("can not close annotate file: $!");#my %distance;
#print "$gene[1]\n";
my $row1 = join("\t",@gene[2..$num-1]);
#print "$row1\n";
print OUT "\n";
print OUT "Distance\t$row1\n";
my $dis;
my $dis1;
my $dis2;
my $i;
my $j;
#print $num."\n";
foreach $name (@gene[2..$num-1]){print OUT "$name\t";#基因1的起始位置key值$l1my $l1 = $left.$name;#基因1的终止位置key值$r1my $r1 = $right.$name;foreach my $name1 (@gene[2..$num-1]){#基因2的起始位置key值$l2my $l2 = $left.$name1;#基因2的起始位置key值$r2my $r2 = $right.$name1;if($name1 eq $gene[$num-1]){if($location{$l1} >= $location{$r2}){$dis1 = $location{$l1} - $location{$r2};$dis2 = $gene[1] - ($location{$r1} - $location{$l2});if($dis1 < $dis2){$dis = $dis1}else{$dis=$dis2}#print $dis."\n";print OUT "$dis\n";}elsif($location{$l2} >= $location{$r1}){$dis1 = $location{$l2} - $location{$r1};$dis2 = $gene[1] - ($location{$r2} - $location{$l1});if($dis1 < $dis2){$dis = $dis1}else{$dis=$dis2}#print $dis."\n";print OUT "$dis\n";}else{$dis = $location{$l1}-$location{$l2};if($dis < 0){$dis = -$dis}#print $dis."\n";print OUT "$dis\n";}next;}if($location{$l1} >= $location{$r2}){$dis1 = $location{$l1} - $location{$r2};$dis2 = $gene[1] - ($location{$r1} - $location{$l2});if($dis1 < $dis2){$dis = $dis1}else{$dis=$dis2}#print $dis."\t";print OUT "$dis\t";}elsif($location{$l2} >= $location{$r1}){$dis1 = $location{$l2} - $location{$r1};$dis2 = $gene[1] - ($location{$r2} - $location{$l1});if($dis1 < $dis2){$dis = $dis1}else{$dis=$dis2}#print $dis."\t";print OUT "$dis\t";}else{$dis = $location{$l1}-$location{$l2};if($dis < 0){$dis = -$dis}#print $dis."\t";print OUT "$dis\t";}}
}
close (OUT) or die ("can not close result file: $!");

原理是利用正则表达式以locus_tag和CDS识别基因编号和起始至终止位置；利用hash数组储存不同目的基因的起始和终止位置；再判断两两基因位置，当两基因A、B不重叠或包含，则计算A终止至B起始碱基数、B终止至A起始之间碱基数，取两者最小值；当两基因A、B重叠或包含时取A起始至B起始之间碱基数。

三、运行代码

以linux为例：

perl distance.pl ZJSH63.tbl genes.txt

三、输出文件

输出文件为result.txt
内容如下：

会打印目的基因的起始和终止位置，再打印两基因距离。

不足

对于NCBI的注释文件，其结构与prokka注释文件结构不一致，如GenBank注释文件xxx.gbff：

该文件识仍可以采用locus_tag和CDS，但是locus_tag后又引号，需要去除，分隔方式也不为制表符；利用CDS识别起始和终止位置也需要特殊处理。
需要直接表征目的基因在基因组上前后关系时，需要更改代码内基因位置判断关系的语句。
（懒得做了…第一次发布，请见谅）

【perl计算基因在基因组上的距离】相关推荐

DSSM、CNN-DSSM、LSTM-DSSM等深度学习模型在计算语义相似度上的应用+距离运算
在NLP领域,语义相似度的计算一直是个难题:搜索场景下query和Doc的语义相似度.feeds场景下Doc和Doc的语义相似度.机器翻译场景下A句子和B句子的语义相似度等等.本文通过介绍DSSM.C ...
计算地球上两点距离（震中距）的Matlab函数（兼容度数和度分秒）及另外三种方法
目录写在前面方法1: taup 方法2: ObsPy 方法3: Mapping Toolbox的distance函数方法4: 自己写的Matlab函数参数公式函数写在前面最近要计算震中 ...
js 计算在AB两点连线上，距离A点一定距离的点的坐标
/*** @description 获取在AB两点连线上,以AB为方向,距离A点,L处的点的坐标* @param A:{x,z} 点A* @param B:{x,z} 点B* @param L 距A点 ...
MPB：微生物所东秀珠组-基于16S rRNA基因和基因组序列对细菌物种的初步鉴定...
为进一步提高<微生物组实验手册>稿件质量,本项目新增大众评审环节.文章在通过同行评审后,采用公众号推送方式分享全文,任何人均可在线提交修改意见.公众号格式显示略有问题,建议电脑端点击文末阅 ...
RNAmmer:预测基因组上的核糖体RNA
核糖体RNA, 缩写为rRNA, 是细胞内含量最多的一类RNA, 能够与蛋白质结合形成核糖体,完成氨基酸的合成.rRNA分子量较大,通常利用沉降系数来区分不同类别的rRNA,沉降系数越大,分子量越大. ...
CB：中国农大胡永飞组构建整合的鸡肠道微生物组的参考基因和基因组集
鸡肠道微生物组的宏基因组组装基因组和基因集助力破译耐药基因组 Metagenome-assembled genomes and gene catalog from the chicken gut mi ...
分子生物学第三章基因、基因组及基因组学
文章目录第三章基因.基因组及基因组学第一节基因 1 基因认识的三个阶段 2 基因的特征 (1)跳跃基因 (2)断裂基因 3 基因的分类 4 基因的结构 5 基因的大小 6 基因的数目第二节 ...
html计算平均分,Calculate phastCon Score for a gene —- 计算基因的phastCon平均分，判断基因保守型...
Calculate phastCon Score for a gene -- 计算基因的phastCon平均分,判断基因保守型 PhastCon socre is the score from 0 t ...
html设置图片与边框的距离,css图片如何设置上边框距离
在css中,可以使用padding-top属性来设置图片的上边框距离,只需要给图片元素添加"padding-top:距离值;"样式即可.padding-top属性可以设置元素的上内 ...

【perl计算基因在基因组上的距离】