lib/Bio/Tools/Est2Genome.pm

   1 #
   2 # BioPerl module for Bio::Tools::Est2Genome
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Jason Stajich <jason-at-bioperl.org>
   7 #
   8 # Copyright Jason Stajich
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11
  12 # POD documentation - main docs before the code
  13
  14 =head1 NAME
  15
  16 Bio::Tools::Est2Genome - Parse est2genome output, makes simple Bio::SeqFeature::Generic objects
  17
  18 =head1 SYNOPSIS
  19
  20   use Bio::Tools::Est2Genome;
  21
  22   my $featureiter = Bio::Tools::Est2Genome->new(-file => 'output.est2genome');
  23
  24   # This is going to be fixed to use the SeqAnalysisI next_feature
  25   # Method eventually when we have the objects to put the data in
  26   # properly
  27   while( my $f = $featureiter->parse_next_gene ) {
  28    # process Bio::SeqFeature::Generic objects here
  29   }
  30
  31 =head1 DESCRIPTION
  32
  33 This module is a parser for C<est2genome> [EMBOSS] alignments of est/cdna
  34 sequence to genomic DNA.  This is generally accepted as the best
  35 program for predicting splice sites based on est/dnas (as far as I know).
  36
  37 This module currently does not try pull out the ungapped alignments
  38 (Segment) but may in the future.
  39
  40 =head1 FEEDBACK
  41
  42 =head2 Mailing Lists
  43
  44 User feedback is an integral part of the evolution of this and other
  45 Bioperl modules. Send your comments and suggestions preferably to
  46 the Bioperl mailing list.  Your participation is much appreciated.
  47
  48   bioperl-l@bioperl.org                  - General discussion
  49   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  50
  51 =head2 Support
  52
  53 Please direct usage questions or support issues to the mailing list:
  54
  55 I<bioperl-l@bioperl.org>
  56
  57 rather than to the module maintainer directly. Many experienced and
  58 reponsive experts will be able look at the problem and quickly
  59 address it. Please include a thorough description of the problem
  60 with code and data examples if at all possible.
  61
  62 =head2 Reporting Bugs
  63
  64 Report bugs to the Bioperl bug tracking system to help us keep track
  65 of the bugs and their resolution. Bug reports can be submitted the
  66 web:
  67
  68   https://github.com/bioperl/bioperl-live/issues
  69
  70 =head1 AUTHOR - Jason Stajich
  71
  72 Email jason-at-bioperl.org
  73
  74 =head1 APPENDIX
  75
  76 The rest of the documentation details each of the object methods.
  77 Internal methods are usually preceded with a _
  78
  79 =cut
  80
  81
  82 # Let the code begin...
  83
  84
  85 package Bio::Tools::Est2Genome;
  86
  87 use strict;
  88
  89 # Object preamble - inherits from Bio::Root::Root
  90
  91 use Bio::Root::Root;
  92 use Bio::SeqFeature::Gene::Exon;
  93 use Bio::SeqFeature::Gene::Transcript;
  94 use Bio::SeqFeature::Gene::Intron;
  95 use Bio::SeqFeature::Gene::GeneStructure;
  96 use Bio::SeqFeature::SimilarityPair;
  97
  98 use base qw(Bio::Tools::AnalysisResult);
  99
 100 =head2 new
 101
 102  Title   : new
 103  Usage   : my $obj = Bio::Tools::Est2Genome->new();
 104  Function: Builds a new Bio::Tools::Est2Genome object
 105  Returns : an instance of Bio::Tools::Est2Genome
 106  Args    : -file => 'output.est2genome' or
 107            -fh   => \*EST2GENOMEOUTPUT
 108            -genomefirst => 1  # genome was the first input (not standard)
 109
 110 =cut
 111
 112 sub _initialize_state {
 113     my($self,@args) = @_;
 114
 115     # call the inherited method first
 116     my $make = $self->SUPER::_initialize_state(@args);
 117
 118     my ($genome_is_first) = $self->_rearrange([qw(GENOMEFIRST)], @args);
 119
 120     delete($self->{'_genome_is_first'});
 121     $self->{'_genome_is_first'} = $genome_is_first if(defined($genome_is_first));
 122     $self->analysis_method("est2genome");
 123 }
 124
 125 =head2 analysis_method
 126
 127  Usage     : $sim4->analysis_method();
 128  Purpose   : Inherited method. Overridden to ensure that the name matches
 129              /est2genome/i.
 130  Returns   : String
 131  Argument  : n/a
 132
 133 =cut
 134
 135 #-------------
 136 sub analysis_method {
 137 #-------------
 138     my ($self, $method) = @_;
 139     if($method && ($method !~ /est2genome/i)) {
 140         $self->throw("method $method not supported in " . ref($self));
 141     }
 142     return $self->SUPER::analysis_method($method);
 143 }
 144
 145 =head2 parse_next_gene
 146
 147  Title   : parse_next_gene
 148  Usage   : @gene = $est2genome_result->parse_next_gene;
 149            foreach $exon (@exons) {
 150                # do something
 151            }
 152
 153  Function: Parses the next alignments of the est2genome result file and
 154            returns the found exons as an array of
 155            Bio::SeqFeature::SimilarityPair objects. Call
 156            this method repeatedly until an empty array is returned to get the
 157            results for all alignments.
 158
 159            The $exon->seq_id() attribute will be set to the identifier of the
 160            respective sequence for both sequences.
 161            The length is accessible via the seqlength()
 162            attribute of $exon->query() and
 163            $exon->est_hit().
 164  Returns : An array (or array reference) of Bio::SeqFeature::SimilarityPair and Bio::SeqFeature::Generic objects
 165            or Bio::SeqFeature::Gene::GeneStructure
 166  Args    : flag(1/0) indicating to return Bio::SeqFeature::Gene::GeneStructure or Bio::SeqFeature::SimilarityPair
 167            defaults to 0
 168
 169 =cut
 170
 171 sub parse_next_gene {
 172    my ($self,$return_gene) = @_;
 173    return $self->_parse_gene_struct if $return_gene;
 174    my $seensegment = 0;
 175    my @features;
 176    my ($qstrand,$hstrand) = (1,1);
 177    my $lasthseqname;
 178    while( defined($_ = $self->_readline) ) {
 179        if( /Note Best alignment is between (reversed|forward) est and (reversed|forward) genome, (but|and) splice\s+sites imply\s+(forward gene|REVERSED GENE)/) {
 180            if( $seensegment ) {
 181                $self->_pushback($_);
 182                return wantarray ? @features : \@features;
 183            }
 184            $hstrand = -1 if $1 eq 'reversed';
 185            $qstrand = -1 if $4 eq 'REVERSED GENE';
 186            #$self->debug( "1=$1, 2=$2, 4=$4\n");
 187        }
 188        elsif( /^Exon/ ) {
 189            my ($name,$score,$perc_ident,$qstart,$qend,$qseqname,
 190                $hstart,$hend, $hseqname) = split;
 191            $lasthseqname = $hseqname;
 192            my $query = Bio::SeqFeature::Similarity->new(-primary => $name,
 193                                                        -source  => $self->analysis_method,
 194                                                        -seq_id => $qseqname, # FIXME WHEN WE REDO THE GENERIC NAME CHANGE
 195                                                        -start   => $qstart,
 196                                                        -end     => $qend,
 197                                                        -strand  => $qstrand,
 198                                                        -score   => $score,
 199                                                        -tag => {
 200 #                                                          'Location' => "$hstart..$hend",
 201                                                            'Sequence' => "$hseqname",
 202                                                            'identity' => $perc_ident,
 203                                                            }
 204                                                        );
 205            my $hit = Bio::SeqFeature::Similarity->new(-primary => 'exon_hit',
 206                                                      -source  => $self->analysis_method,
 207                                                      -seq_id => $hseqname,
 208                                                      -start   => $hstart,
 209                                                      -end     => $hend,
 210                                                      -strand  => $hstrand,
 211                                                      -score   => $score,
 212                                                      -tag => {
 213 #                                                        'Location' => "$qstart..$qend",
 214                                                          'Sequence' => "$qseqname",
 215                                                          'identity' => $perc_ident,
 216                                                      }
 217                                                      );
 218            push @features, Bio::SeqFeature::SimilarityPair->new
 219                (-query => $query,
 220                 -hit   => $hit,
 221                 -source => $self->analysis_method);
 222        } elsif( /^([\-\+\?])(Intron)/) {
 223            my ($name,$score,$perc_ident,$qstart,$qend,$qseqname) = split;
 224            push @features, Bio::SeqFeature::Generic->new(-primary => $2,
 225                                                         -source => $self->analysis_method,
 226                                                         -start => $qstart,
 227                                                         -end   => $qend,
 228                                                         -strand => $qstrand,
 229                                                         -score  => $score,
 230                                                         -seq_id => $qseqname,
 231                                                         -tag => {
 232                                                          'identity' => $perc_ident,
 233                                                             'Sequence' => $lasthseqname});
 234        } elsif( /^Span/ ) {
 235        } elsif( /^Segment/ ) {
 236            $seensegment = 1;
 237        } elsif( /^\s+$/ ) { # do nothing
 238        } else {
 239            $self->warn( "unknown line $_\n");
 240        }
 241    }
 242    return unless( @features );
 243    return wantarray ? @features : \@features;
 244 }
 245
 246 sub _parse_gene_struct {
 247    my ($self) = @_;
 248    my $seensegment = 0;
 249    my @features;
 250    my ($qstrand,$hstrand) = (1,1);
 251    my $lasthseqname;
 252    my $gene = Bio::SeqFeature::Gene::GeneStructure->new(-source => $self->analysis_method);
 253    my $transcript = Bio::SeqFeature::Gene::Transcript->new(-source => $self->analysis_method);
 254    my @suppf;
 255    my @exon;
 256    while( defined($_ = $self->_readline) ) {
 257        if( /Note Best alignment is between (reversed|forward) est and (reversed|forward) genome, (but|and) splice\s+sites imply\s+(forward gene|REVERSED GENE)/) {
 258               if( $seensegment ) {
 259                $self->_pushback($_);
 260                return $gene;
 261               }
 262            $hstrand = -1 if $1 eq 'reversed';
 263            $qstrand = -1 if $4 eq 'REVERSED GENE';
 264        }
 265        elsif( /^Exon/ ) {
 266            my ($name,$score,$perc_ident,$qstart,$qend,$qseqname,$hstart,$hend, $hseqname) = split;
 267            $lasthseqname = $hseqname;
 268            my $exon = Bio::SeqFeature::Gene::Exon->new(-primary => $name,
 269                                                        -source  => $self->analysis_method,
 270                                                        -seq_id => $qseqname, # FIXME WHEN WE REDO THE GENERIC NAME CHANGE
 271                                                        -start   => $qstart,
 272                                                        -end     => $qend,
 273                                                        -strand  => $qstrand,
 274                                                        -score   => $score,
 275                                                        -tag => {
 276                             #'Location' => "$hstart..$hend",
 277                                                          'identity' => $perc_ident,
 278                                                                    'Sequence' => "$hseqname",
 279                                                                       }
 280                                                        );
 281           $transcript->seq_id($qseqname) unless $transcript->seq_id;
 282           $exon->add_tag_value('phase',0);
 283           push @exon, $exon;
 284
 285        } elsif( /^([\-\+\?])(Intron)/) {
 286          next; #intron auto matically built from exons..hope thats ok..
 287        } elsif( /^Span/ ) {
 288        } elsif( /^Segment/ ) {
 289             my ($name,$score,$perc_ident,$qstart,$qend,$qseqname,$hstart,$hend, $hseqname) = split;
 290                  my $query = Bio::SeqFeature::Similarity->new(-primary => $name,
 291                                                        -source  => $self->analysis_method,
 292                                                        -seq_id => $qseqname, # FIXME WHEN WE REDO THE GENERIC NAME CHANGE
 293                                                        -start   => $qstart,
 294                                                        -end     => $qend,
 295                                                        -strand  => $qstrand,
 296                                                        -score   => $score,
 297                                                        -tag => {
 298 #                                                          'Location' => "$hstart..$hend",
 299                                                            'Sequence' => "$hseqname",
 300                                                          'identity' => $perc_ident,
 301                                                            }
 302                                                      );
 303            my $hit = Bio::SeqFeature::Similarity->new(-primary => 'exon_hit',
 304                                                                                              -source  => $self->analysis_method,
 305                                                                                              -seq_id => $hseqname,
 306                                                                                              -start   => $hstart,
 307                                                                                              -end     => $hend,
 308                                                                                              -strand  => $hstrand,
 309                                                                                              -score   => $score,
 310                                                                                              -tag => {
 311                                                             #   'Location' => "$qstart..$qend",
 312                                                                                                  'Sequence' => "$qseqname",
 313                                                          'identity' => $perc_ident,
 314                                                                                                 }
 315                                                      );
 316                    my $support =  Bio::SeqFeature::SimilarityPair->new(-query => $query,
 317                                                                         -hit   => $hit,
 318                                                                         -source => $self->analysis_method);
 319              push @suppf, $support;
 320        } elsif( /^\s+$/ ) { # do nothing
 321        } else {
 322            $self->warn( "unknown line $_\n");
 323        }
 324    }
 325    return unless $#exon >=0;
 326    foreach my $e(@exon){
 327     my @add;
 328     foreach my $sf(@suppf){
 329       if($sf->overlaps($e)){
 330           push @add,$sf;
 331       }
 332     }
 333     $e->add_tag_value('supporting_feature',@add);
 334     $transcript->add_exon($e);
 335   }
 336
 337    $gene->add_transcript($transcript);
 338    $gene->seq_id($transcript->seq_id);
 339    return $gene;
 340 }
 341
 342 =head2 next_feature
 343
 344  Title   : next_feature
 345  Usage   : $seqfeature = $obj->next_feature();
 346  Function: Returns the next feature available in the analysis result, or
 347            undef if there are no more features.
 348  Example :
 349  Returns : A Bio::SeqFeatureI implementing object, or undef if there are no
 350            more features.
 351  Args    : none
 352
 353 =cut
 354
 355 sub next_feature {
 356     my ($self) = shift;
 357     $self->throw("We haven't really done this right, yet, use parse_next_gene");
 358 }
 359
 360
 361 1;