perl如何輸出基因的位置信息

這篇文章主要介紹了perl如何輸出基因的位置信息的相關知識，內容詳細易懂，操作簡單快捷，具有一定借鑒價值，相信大家閱讀完這篇perl如何輸出基因的位置信息文章都會有所收獲，下面我們一起來看看吧。

perl輸出基因的位置信息按照基因所在染色體，和位置信息排序

我們在整理基因組的gff文件，想輸出基因的位置信息，以及基因所對應的多個轉錄本信息，需要對基因按照染色體排序，這里使用到了perl里面hash按照值來排序，而且還用了兩個值基因型排序。示例代碼如下，以下代碼可以從gff文件當中提取所有基因的位置信息以及對應的多個轉錄本信息：

perl代碼如下：

#!/usr/bin/perl -w
use strict;
use Cwd qw(abs_path getcwd);
use Getopt::Long;
use Data::Dumper;

die "perl $0 <gff> <outfile>" unless(@ARGV==2);


my$gff=$ARGV[0];
my%gene=();
my%gene_region=();
my%mRNA2Gene=();
my%Gene2mRNA=();
open IN,"$gff" or die "$!";
open OUT ,">$ARGV[1]" or die "$!";
print OUT "#gene_ID\tchr\tstart\tend\tstrand\ttranscript_id\n";
while(<IN>){
chomp;
next if (/^#/);
my@tmp=split(/\t/);


if($tmp[2] =~/^gene/){
my($id)=($tmp[8]=~/ID=([^;]+)/);
$gene{$id}=1;
$gene_region{$id}=[$tmp[0],$tmp[3],$tmp[4],$tmp[6]];


#print "gene:$id\n";
#my$gene_chr->{$id}=$tmp[0];
}
if($tmp[2] =~/mRNA|transcript/i){
my($id)=($tmp[8]=~/ID=([^;]+)/);
my($pid)=($tmp[8]=~/Parent=([^;]+)/);


if(exists $gene{$pid}){

push @{$Gene2mRNA{$pid}},$id;
}
#print "mRNA:$id\n";
}
}
close(IN);
#多層排序，先按染色體排序，再按基因位置排序
for my$id(sort {$gene_region{$a}->[0] cmp $gene_region{$b}->[0] 
or $gene_region{$a}->[1] <=> $gene_region{$b}->[1] } keys %gene_region){
print OUT "$id\t".join("\t",(@{$gene_region{$id}},sort @{$Gene2mRNA{$id}})    )."\n";
}
close(OUT);

關于“perl如何輸出基因的位置信息”這篇文章的內容就介紹到這里，感謝各位的閱讀！相信大家對“perl如何輸出基因的位置信息”知識都有一定的了解，大家如果還想學習更多知識，歡迎關注億速云行業資訊頻道。