欢迎光临
我们一直在努力

perl如何提取GFF中所有转录本的位置信息

本篇内容主要讲解“perl如何提取GFF中所有转录本的位置信息”,感兴趣的朋友不妨来看看。本文介绍的方法操作简单快捷,实用性强。下面就让小编来带大家学习“perl如何提取GFF中所有转录本的位置信息”吧!

提取基因组注释文件GFF中所有基因转录本的位置信息,以及转录本对应的基因的ID:

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=();
open IN,"$gff" or die "$!";
open OUT ,">$ARGV[1]" or die "$!";
print OUT "#mRNA_ID\tgene_ID\tchr\tstart\tend\tstrand\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}){
print OUT "$id\t$pid\t$tmp[0]\t$tmp[3]\t$tmp[4]\t$tmp[6]\n";
}
#print "mRNA:$id\n";
}
}
close(IN);
close(OUT);

到此,相信大家对“perl如何提取GFF中所有转录本的位置信息”有了更深的了解,不妨来实际操作一番吧!这里是云搜网网站,更多相关内容可以进入相关频道进行查询,关注我们,继续学习!

赞(0)
【声明】:本博客不参与任何交易,也非中介,仅记录个人感兴趣的主机测评结果和优惠活动,内容均不作直接、间接、法定、约定的保证。访问本博客请务必遵守有关互联网的相关法律、规定与规则。一旦您访问本博客,即表示您已经知晓并接受了此声明通告。