[Bioperl-l] What's wrong with the script about module Bio::SeqFeature::Gene::GeneStructure

[Hans-Rudolf Hotz]

Hi Tao

I don't fully understand your script, but I do see a major problem:

why do you select for "primary_tag eq 'mRNA'" first?


this simple loop will probably do what you want:

while( my $seq_obj = $catch_seq -> next_seq) {

	my @all_exon_features = grep {$_->primary_tag eq 'exon'}
$seq_obj ->get_SeqFeatures;
	$exon_number = scalar(@all_exon_features);
	print "$exon_number\n";
}


Alternatively, I recommend to read the HOWTO page:

http://www.bioperl.org/wiki/HOWTO:Feature-Annotation

which has a nice example to print out all 'primary_tag' (ie the 'feature 
keys' in a GenBank formated file)

I also recommend to read up on the DDBJ/EMBL/GenBank Feature Table 
definition:
http://www.ncbi.nlm.nih.gov/projects/collab/FT/index.html


Regards, Hans



On 03/18/2011 05:33 AM, Tao Zhu wrote:
> I wrote my script like this,
>
> #!/usr/bin/perl -w
> use Bio::SeqIO;
> my $catch_seq = Bio::SeqIO ->  new(-file =>  'test.gbk',-format=>
> 'genbank');
> while( my $seq_obj = $catch_seq ->  next_seq)
> {
> 	my @all_mRNA_features = grep {$_->primary_tag eq 'mRNA'} $seq_obj ->
> get_SeqFeatures;
> 	for my $mRNA_feature (@all_mRNA_features)
> 	{
> 			if($mRNA_feature->isa('Bio::SeqFeature::Gene::GeneStructure'))
> 		{
> 			@exons=$mRNA_feature->exons;
> 			$exon_number = scalar(@exons);
> 			print "$exon_number\n";
> 		}
> 	}
> }
>
> I hope to count exon number in every mRNA. But it print nothing(You can
> arbitrarily get a genbank file to test it). What's wrong?
>
>