#!/usr/bin/env perl
use warnings;
use strict;
use Bio::SeqIO;
die "perl $0 <fasta> > <outfile>\n" if(@ARGV != 1);
my @len = `fastalength $ARGV[0]`;
my %xsh = map { chomp; my($length,$name)=split /\s+/; $name=>$length } @len;
my %hash;
open FA, $ARGV[0] or die $!;
while(<FA>)
{
chomp;
next if($_ !~ /^>/);
if(/gene:([^ ]+) transcript:([^ ]+)/)
{
if(exists $xsh{$2})
{
push @{$hash{$1}}, "$2 $xsh{$2}";
}
}
}
my %trans;
open OUT, ">longest_transcript.list" or die $!;
foreach my $g(keys %hash)
{
my $tmp = 0;
my $lt;
foreach my $t(@{$hash{$g}})
{
my @x = split / /, $t;
($x[1] > $tmp) ? $lt = $x[0] : next;
}
$trans{$lt} = 0;
print OUT "$lt\n";
}
$/ = "\n>";
open FB, $ARGV[0] or die $!;
while(<FB>)
{
chomp;
s/^>//;
my @tmp = split /\n/;
my @t = s
取转录本fasta最长的当作基因fasta
最新推荐文章于 2024-08-17 18:17:45 发布
这段Perl脚本从多个转录本fasta文件中找出最长的转录本,将其作为对应基因的fasta序列。脚本首先计算每个转录本的长度,然后保存最长的转录本信息,并最终输出最长的转录本序列。
摘要由CSDN通过智能技术生成