File Coverage

blib/lib/Bio/ViennaNGS/AnnoC.pm

Criterion	Covered	Total	%
statement	30	96	31.2
branch	0	26	0.0
condition	0	8	0.0
subroutine	10	13	76.9
pod	3	3	100.0
total	43	146	29.4

# -*-CPerl-*-
# Last changed Time-stamp: <2014-10-02 19:17:12 mtw>

package Bio::ViennaNGS::AnnoC;

use Exporter;
use version; our $VERSION = qv('0.06');
use strict;
use strict;
use warnings;
use Data::Dumper;
use Bio::Tools::GFF;
use Bio::DB::Fasta;
use Path::Class;
use Carp;

our @ISA       = qw(Exporter);
our @EXPORT_OK = qw(parse_gff feature_summary get_fasta_ids
		    $feat $fstat $fastadb
		    @fastaids
		    %features %featsta);
our @EXPORT    = ();

#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#
#^^^^^^^^^^ Variables ^^^^^^^^^^^#
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#

our ($fastadb);
our %features  = ();
our %featstat  = ();
our $feat      = \%features;
our $fstat     = \%featstat;
our @fastaids  = ();

#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#
#^^^^^^^^^^^ Subroutines ^^^^^^^^^^#
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#

sub parse_gff {
  my $in_file = shift;
  my ($i,$gffio,$feature,$gbkey);
  my $this_function = (caller(0))[3];

  $gffio = Bio::Tools::GFF->new(-file        => $in_file,
				-gff_version  => 3,
			       );
  $gffio->ignore_sequence(1);
  if (my $header = $gffio->next_segment() ){
    $featstat{accession}= $header->display_id();
  }
  else{
    carp "[$this_function]: Cannot parse GFF header\n";
  }

  while($feature = $gffio->next_feature()) {
    my ($uid,$feat_name);
    my @name = my @id = my @gbkeys = ();

    next if ($feature->primary_tag() eq "exon");

    # 1) determine gbkey of the current feature
    @gbkeys = $feature->get_tag_values("gbkey");
    $gbkey  = $gbkeys[0];

    # 2) get a unique ID for each feature
    if ($feature->has_tag('ID')){
      @id = $feature->get_tag_values('ID');
      $uid = $id[0]; # ID=id101
    }
    else {
      croak "ERROR [$this_function] Feature '$gbkey' at pos.\
             $feature->start does not have \'ID\' attribute\n";
    }

    # 3) assign parent's unique ID in case a parent record exists
    if ($feature->has_tag('Parent')){
      @id = $feature->get_tag_values('Parent');
      $uid = $id[0]; # ID=id101
    }

    # 4) find a name for the current feature, use 'Name' or 'ID' attribute
    if ($feature->has_tag('Name')){
      @name = $feature->get_tag_values('Name');
      $feat_name = $name[0];
    }
    elsif ($feature->has_tag('ID')){
      @id = $feature->get_tag_values('ID');
      $feat_name = $id[0]; # ID=id101, use ID as feature name
    }
    else {
      croak "ERROR [$this_function] Cannot set name for feature \
              $feature->gbkey at pos. $feature->start\n";
    }

    unless (exists $features{$uid}) { # gene / ribosome_entry_site / etc.
      $features{$uid}->{start}     = $feature->start;
      $features{$uid}->{end}       = $feature->end;
      $features{$uid}->{strand}    = $feature->strand;
      $features{$uid}->{length}    = $feature->length;
      $features{$uid}->{seqid}     = $feature->seq_id;
      $features{$uid}->{score}     = $feature->score || 0;
      $features{$uid}->{gbkey}     = $gbkey;
      $features{$uid}->{name}      = $feat_name;
      $features{$uid}->{uid}       = $uid;
    }
    else { # CDS / tRNA / rRNA / etc
      $features{$uid}->{gbkey} = $gbkey;  # gbkey for tRNA/ rRNA/ CDS etc
    }
  }

  # finally generate some statistics on features present in this annotation
  $featstat{total} = 0;
  $featstat{origin} = "$this_function ".$VERSION;
  foreach my $ft (keys %features){
    $featstat{total}++;
    my $key = $features{$ft}->{gbkey};
    unless (exists $featstat{$key}){
      $featstat{$key} = 0;
    }
    $featstat{$key} += 1;
  }
  $gffio->close();
 # print Dumper (\%features);
}

# feature_summary($fsR,$dest)
# Print summary of %featstat hash
#
# ARG1: reference to %featstat hash
# ARG2: path for output file
sub feature_summary {
  my ($fsR, $dest) = @_;
  my ($fn,$fh);
  my $this_function = (caller(0))[3];

  croak "ERROR [$this_function] $dest does not exist\n"
    unless (-d $dest);

  $fn = dir($dest,$$fsR{accession}.".summary.txt");
  open $fh, ">", $fn or croak $!;

  print $fh "Accession\t $$fsR{accession}\n";
  print $fh "Origin   \t $$fsR{origin}\n";
  foreach my $ft (sort keys %$fsR){
    next if ($ft =~ /total/ || $ft =~ /accession/ || $ft =~ /origin/);
    print $fh "$ft\t$$fsR{$ft}\n";
  }
  print $fh "Total\t$$fsR{total}\n";
  close $fh;
}

# get_fasta_ids($fa_in)
# Returns an array with all fasta ids from a (multi)fasta file
#
# ARG1: fa_in
sub get_fasta_ids {
  my $fa_in = shift;
  unless (defined $fa_in){
    confess "Fasta input not available";
  }
  $fastadb = Bio::DB::Fasta->new($fa_in) or die $!;
  @fastaids = $fastadb->ids;
  return @fastaids;
}

1;

__END__