File Coverage

blib/lib/CracTools/Annotator.pm

Criterion	Covered	Total	%
statement	171	221	77.3
branch	33	86	38.3
condition	3	30	10.0
subroutine	23	24	95.8
pod	12	12	100.0
total	242	373	64.8

package CracTools::Annotator;
{
  $CracTools::Annotator::DIST = 'CracTools';
}
# ABSTRACT: Generic annotation base on CracTools::GFF::Query::File
$CracTools::Annotator::VERSION = '1.25';
use strict;
use warnings;

use Carp;
use List::Util qw[min max];

use CracTools::Const;
use CracTools::GFF::Annotation;
use CracTools::Interval::Query;
use CracTools::Interval::Query::File;


sub new {
  my $class = shift;
  my $gff_file = shift;
  my $mode = shift;

  $mode = 'light' if !defined $mode;

  if(!defined $gff_file) {
    croak "Missing GFF file argument in CracTools::Annotator constructor";
  }

  my $self = bless {
    gff_file => $gff_file,
    mode => $mode,
  }, $class;

  $self->_init();

  return $self;
}


sub mode {
  my $self = shift;
  return $self->{mode};
}


sub foundAnnotation {
  my $self = shift;
  #my ($chr,$pos_start,$pos_end,$strand) = @_;
  my @candidates = @{ $self->getAnnotationCandidates(@_)};
  return (scalar @candidates > 0);
}


sub foundGene {
  my $self = shift;
  my ($chr,$pos_start,$pos_end,$strand) = @_;
  my @candidates = @{ $self->getAnnotationCandidates($chr,$pos_start,$pos_end,$strand)};
  foreach my $candidate (@candidates) {
    return 1 if defined $candidate->{gene};
  }
  return 0;
}


sub foundSameGene {
  my $self = shift;
  my ($chr,$pos_start1,$pos_end1,$pos_start2,$pos_end2,$strand) = @_;
  my @candidates1 = @{ $self->getAnnotationCandidates($chr,$pos_start1,$pos_end1,$strand)};
  my @candidates2 = @{ $self->getAnnotationCandidates($chr,$pos_start2,$pos_end2,$strand)};
  my $found_same_gene = 0;
  my @genes1;
  my @genes2;
  foreach my $candi1 (@candidates1) {
    if(defined $candi1->{gene}) {
      push @genes1,$candi1->{gene}->attribute('ID');
    }
  }
  foreach my $candi2 (@candidates2) {
    if(defined $candi2->{gene}) {
      push @genes2,$candi2->{gene}->attribute('ID');
    }
  }
  foreach my $gene_id (@genes1) {
    foreach (@genes2) {
      if($gene_id eq $_) {
        $found_same_gene = 1;
        last;
      }
    }
    last if $found_same_gene == 1;
  }
  return $found_same_gene;
}


sub getBestAnnotationCandidate {
  my $self = shift;
  my ($best_candidates,$best_priority,$best_type) = $self->getBestAnnotationCandidates(@_);
  if(@{$best_candidates}) {
    return $best_candidates->[0],$best_priority,$best_type;
  } else {
    return undef,undef,undef;
  }
}


sub getBestAnnotationCandidates {
  my $self = shift;
  my ($chr,$pos_start,$pos_end,$strand,$prioritySub,$compareSub) = @_;

  if(!defined $prioritySub && !defined $compareSub) {
    $prioritySub = \&getCandidatePriorityDefault unless defined $prioritySub;
    $compareSub = \&compareTwoCandidatesDefault unless defined $compareSub;
  }

  my @candidates = @{ $self->getAnnotationCandidates($chr,$pos_start,$pos_end,$strand)};
  my @best_candidates;
  my ($best_priority,$best_type);
  foreach my $candi (@candidates) {
    my ($priority,$type);
    ($priority,$type) = $prioritySub->($pos_start,$pos_end,$candi) if defined $prioritySub;
    if(defined $priority && $priority != -1) {
      if(!defined $best_priority) {
        $best_priority = $priority;
        push @best_candidates, $candi;
        $best_type = $type;
      } elsif($priority < $best_priority) {
        @best_candidates = ($candi);
        $best_priority = $priority;
        $best_type = $type;
      }
      #we should compare two candidates with equal priority to always choose the one
      elsif (!defined $priority || $priority == $best_priority){
        my $candidate_chosen;
        my $found_better_candidate = 0;
        foreach my $best_candidate (@best_candidates) {
          $candidate_chosen = $compareSub->($best_candidate,$candi,$pos_start,$pos_end) if defined $compareSub;
          # They are both equal
          if (!defined $candidate_chosen) {
            # We cannnot say if this candidate is better
            next;
          } elsif ($candidate_chosen == $candi) {
            # We have found a better candidate that previously register ones
            # we save it and remove the others
            @best_candidates = ($candi);
            $found_better_candidate = 1;
            last;
          } else {
            # The better candidate is not "candi", so this candidates
            # does not belong the the best_candidate array.
            # We can stop looping
            $found_better_candidate = 1;
            last;
          }
        }
        push @best_candidates, $candi if !$found_better_candidate;
      }
    }
  }
  # TODO We should not return variable in that order,
  # it is not easy to only retrieve the best candidatse...
  return \@best_candidates,$best_priority,$best_type;
}


sub getAnnotationCandidates {
  my $self = shift;
  my ($chr,$pos_start,$pos_end,$strand) = @_;
  # TODO if no strand is provided we should return annotations from both strands

  # get GFF annotations that overlap the region to annotate
  my $annotations = $self->{gff_query}->fetchByRegion($chr,$pos_start,$pos_end,$strand);
  # get a ref of an array of hash of candidates
  my $candidatates = $self->_constructCandidatesFromAnnotation($annotations);
  return $candidatates;
}


sub getAnnotationNearestDownCandidates {
  my $self = shift;
  my ($chr,$pos_start,$strand) = @_;

  # get GFF annotations that overlap the pos_start to annotate
  my $annotations_overlap = $self->{gff_query}->fetchByLocation($chr,$pos_start,$strand);
  # get GFF annotations of nearest down intervals that not overlaped [pos_start,pos_end] pos
  my @annotations_down;

  push @annotations_down, @{$self->{gff_query}->fetchAllNearestDown($chr,$pos_start,$strand)};

  # get a ref of an array of hash of candidates
  my @annotations = (@$annotations_overlap,@annotations_down);
  my $candidatates = $self->_constructCandidatesFromAnnotation(\@annotations);
  return $candidatates;
}


sub getAnnotationNearestUpCandidates {
  my $self = shift;
  my ($chr,$pos_end,$strand) = @_;

  # get GFF annotations that overlap the pos_end to annotate
  my $annotations_overlap = $self->{gff_query}->fetchByLocation($chr,$pos_end,$strand);
  # get GFF annotations of nearest up intervals that not overlaped [pos_start,pos_end] pos
  my @annotations_up;

  push @annotations_up, @{$self->{gff_query}->fetchAllNearestUp($chr,$pos_end,$strand)};

  # get a ref of an array of hash of candidates
  my @annotations = (@$annotations_overlap,@annotations_up);
  my $candidatates = $self->_constructCandidatesFromAnnotation(\@annotations);
  return $candidatates;
}


sub getCandidatePriorityDefault {
  my ($pos_start,$pos_end,$candidate) = @_;
  my ($priority,$type) = (-1,'');
  my ($mRNA,$exon) = ($candidate->{mRNA},$candidate->{exon});
  if(defined $mRNA) {
    if(defined $mRNA->attribute('type') && $mRNA->attribute('type') =~ /protein_coding/i) {
      if(defined $exon) {
        if(($exon->start <= $pos_start) && ($exon->end >= $pos_end)) {
          $priority = 1;
          if(defined $candidate->{three}) {
            $type = '3PRIM_UTR';
          } elsif(defined $candidate->{five}) {
            $type = '5PRIM_UTR';
          # } elsif(defined $candidate->{cds}) {
          #   $type = 'CDS';
          } else {
            $type = 'EXON';
          }
        } else {
          $priority = 2;
          $type = 'INXON';
        }
      } else {
        $priority = 4;
        $type = 'INTRON';
      }
    } else {
      if(defined $exon) {
        if(($exon->start <= $pos_start) && ($exon->end >= $pos_end)) {
          $priority = 3;
          $type = 'NON_CODING';
        }
      }
    }
  }
  return ($priority,$type);
}

sub compareTwoCandidatesDefault{
  my ($candidate1,$candidate2,$pos_start) = @_;
  # If both candidates are exons we try to find wich one is closer to the pos_start of the region to annotate
  if ($candidate1->{exon} && $candidate2->{exon}) {
    my $dist1= min(abs($candidate1->{exon}->end - $pos_start),abs($candidate1->{exon}->start - $pos_start));
    my $dist2= min(abs($candidate2->{exon}->end - $pos_start),abs($candidate2->{exon}->start - $pos_start));
    if ($dist1 > $dist2) {
      return $candidate2;
    } elsif ($dist1 < $dist2) {
      return $candidate1;
    }
  }
  # If we have not found a better candidate, we use the lexicographic order of the mRNA ID
  my ($mRNA1,$mRNA2) = ($candidate1->{mRNA},$candidate2->{mRNA});
  if(defined $mRNA1 && defined $mRNA1->attribute('ID') && defined $mRNA2 && defined $mRNA2->attribute('ID')) {
    if($mRNA1->attribute('ID') lt $mRNA2->attribute('ID')) {
      return $candidate1;
    } else {
      return $candidate2;
    }
  }
  # If nothing has worked we return "undef"
  return undef;
}


sub _init {
  my $self = shift;
  my $gff_query;

  # Create a GFF file to query exons
  if($self->mode eq "fast") {
    $gff_query = CracTools::Interval::Query->new();
    my $gff_it = CracTools::Utils::getFileIterator(file => $self->{gff_file},
      parsing_method => sub { CracTools::GFF::Annotation->new(@_) },
      header_regex => "^#",
    );
    while(my $gff_annot = $gff_it->()) {
      $gff_query->addInterval($gff_annot->chr,
        $gff_annot->start+1,
        $gff_annot->end+1,
        $gff_annot->strand,
        $gff_annot,
      );
    }
  } else {
    $gff_query = CracTools::Interval::Query::File->new(file => $self->{gff_file}, type => 'gff');
  }

  $self->{gff_query} = $gff_query;
}


sub _constructCandidates {
  my ($annot_id,$candidate,$annot_hash) = @_;

  # We init the "leaf_feature" value if this is the first recursion step
  $candidate->{leaf_feature} = $annot_hash->{$annot_id}->feature if !defined $candidate->{leaf_feature};

  my @candidates;
  if (!defined $annot_hash->{$annot_id}){
      carp("Missing feature for $annot_id in the gff file");
  }
  $candidate->{$annot_hash->{$annot_id}->feature} = $annot_hash->{$annot_id};
  my $parents = $annot_hash->{$annot_id}->parents;
  if(@$parents) {
    foreach my $parent (@{$parents}) {

      #Test to avoid a deep recursion
      if($parent eq $annot_id) {
  carp("Parent could not be the candidat itself, please check your gff file for $annot_id");
  next;
      # If there is already a parent with this feature type we duplicated
      # the candidate since we are branching in the annotation tree
      }elsif(!defined $annot_hash->{$parent}) {
        carp("Parent not found, please check your gff file for $annot_id (Parent: $parent)");

      }elsif(defined $candidate->{$annot_hash->{$parent}->feature}) {
        my %copy_candidate = %{$candidate};
        my %copy_parent_feature = %{$candidate->{parent_feature}};
        $copy_candidate{parent_feature} = \%copy_parent_feature;
        # We register in parent_feature links
        $copy_candidate{parent_feature}->{$annot_hash->{$annot_id}->feature} = $annot_hash->{$parent}->feature;
        my $copy_ref = \%copy_candidate;
        push(@candidates,@{_constructCandidates($parent,$copy_ref,$annot_hash)});
      # If not we only go up to the parent node in order to continue candidate
      # construction
      } else {
        # We register in parent_feature links
        $candidate->{parent_feature}->{$annot_hash->{$annot_id}->feature} = $annot_hash->{$parent}->feature;
        push(@candidates,@{_constructCandidates($parent,$candidate,$annot_hash)});
      }
    }
    return \@candidates;
  } else {
    return [$candidate];
  }
}


sub _constructCandidatesFromAnnotation {
  my $self = shift;
  my $annotations = shift;
  my %annot_hash = ();
  my @candidates = ();

  # Construct annotation hash with annot ID as key
  foreach my $annot_line (@{$annotations}) {
    if($self->mode eq "fast") {
      $annot_hash{$annot_line->attribute('ID')} = $annot_line;
    } else {
      my $annot = CracTools::GFF::Annotation->new($annot_line,'gff3');
      $annot_hash{$annot->attribute('ID')} = $annot;
    }
  }

  # Find leaves in annotation tree
  my %hash_leaves;
  foreach my $annot_id (keys %annot_hash) {
    #my @parents = $annot_hash{$annot_id}->parents;
    foreach my $parent (@{$annot_hash{$annot_id}->parents}){
      $hash_leaves{$parent} = 1 unless (defined $hash_leaves{$parent});
    }
  }
  foreach my $annot_id (keys %annot_hash) {
    # check if annot_id is a leaf
    if (!defined $hash_leaves{$annot_id}){
      # Get all possible path from this leaf to the root
      push @candidates, @{_constructCandidates($annot_id,my $new_candidate,\%annot_hash)};
    }
  }

  return \@candidates;
}

1;

__END__