lib/Bio/SeqIO/tigrxml.pm

   1 #
   2 # BioPerl module for Bio::SeqIO::tigrxml
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Jason Stajich <jason-at-bioperl-dot-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::SeqIO::tigrxml - Parse TIGR (new) XML
  17
  18 =head1 SYNOPSIS
  19
  20   use Bio::SeqIO;
  21   my $in = Bio::SeqIO->new(-format => 'tigrcoordset',
  22                           -file   => 'file.xml');
  23
  24   while( my $seq = $in->next_seq ) {
  25      # do something...
  26   }
  27
  28 =head1 DESCRIPTION
  29
  30 This is a parser for TIGR Coordset XML for their in-progress
  31 annotation dbs.
  32
  33 =head1 FEEDBACK
  34
  35 =head2 Mailing Lists
  36
  37 User feedback is an integral part of the evolution of this and other
  38 Bioperl modules. Send your comments and suggestions preferably to
  39 the Bioperl mailing list.  Your participation is much appreciated.
  40
  41   bioperl-l@bioperl.org                  - General discussion
  42   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  43
  44 =head2 Support
  45
  46 Please direct usage questions or support issues to the mailing list:
  47
  48 I<bioperl-l@bioperl.org>
  49
  50 rather than to the module maintainer directly. Many experienced and
  51 reponsive experts will be able look at the problem and quickly
  52 address it. Please include a thorough description of the problem
  53 with code and data examples if at all possible.
  54
  55 =head2 Reporting Bugs
  56
  57 Report bugs to the Bioperl bug tracking system to help us keep track
  58 of the bugs and their resolution. Bug reports can be submitted via the
  59 web:
  60
  61   https://github.com/bioperl/bioperl-live/issues
  62
  63 =head1 AUTHOR - Jason Stajich
  64
  65 Email jason-at-bioperl-dot-org
  66
  67 =head1 APPENDIX
  68
  69 The rest of the documentation details each of the object methods.
  70 Internal methods are usually preceded with a _
  71
  72 =cut
  73
  74
  75 # Let the code begin...
  76
  77
  78 package Bio::SeqIO::tigrxml;
  79 use vars qw($Default_Source);
  80 use strict;
  81 use XML::SAX;
  82 use XML::SAX::Writer;
  83 use Data::Dumper;
  84 use Bio::Seq::SeqFactory;
  85 use Bio::Species;
  86 use Bio::SeqFeature::Generic;
  87 use Bio::Annotation::Reference;
  88 use Bio::Annotation::Comment;
  89 use Bio::Annotation::DBLink;
  90 use List::Util qw(min max);
  91
  92 use base qw(Bio::SeqIO XML::SAX::Base);
  93
  94
  95 $Default_Source = 'TIGR';
  96
  97 sub _initialize {
  98     my ($self) = shift;
  99     $self->SUPER::_initialize(@_);
 100     $self->{'_parser'} = XML::SAX::ParserFactory->parser('Handler' => $self);
 101     if( ! defined $self->sequence_factory ) {
 102         $self->sequence_factory(Bio::Seq::SeqFactory->new
 103                                 (-verbose => $self->verbose(),
 104                                  -type => 'Bio::Seq::RichSeq'));
 105     }
 106     return;
 107 }
 108
 109 sub next_seq {
 110     my $self = shift;
 111     if( @{$self->{'_seendata'}->{'_seqs'} || []} ||
 112         eof($self->_fh)) {
 113         return shift @{$self->{'_seendata'}->{'_seqs'}};
 114     }
 115     $self->{'_parser'}->parse_file($self->_fh);
 116     return shift @{$self->{'_seendata'}->{'_seqs'}};
 117 }
 118
 119 # XML::SAX::Base methods
 120
 121 sub start_document {
 122     my ($self,$doc) = @_;
 123     $self->{'_seendata'} = {'_seqs'    => [],
 124                             '_authors' => [],
 125                             '_feats'   => [] };
 126     $self->SUPER::start_document($doc);
 127 }
 128
 129 sub end_document {
 130     my ($self,$doc) = @_;
 131     $self->SUPER::end_document($doc);
 132 }
 133
 134 sub start_element {
 135     my ($self,$ele) = @_;
 136     # attributes
 137     my $name = uc $ele->{'LocalName'};
 138     my $attr = $ele->{'Attributes'};
 139     my $seqid = defined $self->{'_seendata'}->{'_seqs'}->[-1] ?
 140         $self->{'_seendata'}->{'_seqs'}->[-1]->display_id : undef;
 141
 142     # we're going to try and be SO-nice here
 143     if( $name eq 'ASSEMBLY' ) { # New sequence
 144         my ($len) = $attr->{'{}COORDS'}->{'Value'} =~ /\d+\-(\d+)/;
 145         push @{$self->{'_seendata'}->{'_seqs'}},
 146         $self->sequence_factory->create
 147             (
 148              -display_id => $attr->{'{}ASMBL_ID'}->{'Value'},
 149              -length     => $len,
 150              );
 151     } elsif( $name eq 'HEADER' ) {
 152     } elsif( $name eq 'CLONE_NAME' ) {
 153     } elsif( $name eq 'ORGANISM' ) {
 154     } elsif( $name eq 'AUTHOR_LIST' ) {
 155         $self->{'_seendata'}->{'_authors'} = [];
 156     } elsif( $name eq 'TU' ) { # gene feature
 157         my ($s,$e) = ($attr->{'{}COORDS'}->{'Value'} =~ /(\d+)\-(\d+)/);
 158         my $strand = 1;
 159         if( $s > $e) {
 160             ($s,$e,$strand) = ( $e,$s,-1);
 161         }
 162         my $fname = $attr->{'{}FEAT_NAME'}->{'Value'};
 163         my $f = Bio::SeqFeature::Generic->new
 164             (-seq_id      => $seqid,
 165              -start       => $s,
 166              -end         => $e,
 167              -strand      => $strand,
 168              -primary_tag => 'gene', # what does this really map to?
 169              -source_tag  => $Default_Source,
 170              -tag         => {
 171                  'Note'         => $attr->{'{}COM_NAME'}->{'Value'},
 172                  'ID'           => $fname,
 173                  'locus'        => $attr->{'{}LOCUS'}->{'Value'},
 174                  'pub_locus'    => $attr->{'{}PUB_LOCUS'}->{'Value'},
 175                  'alt_locus'    => $attr->{'{}ALT_LOCUS'}->{'Value'},
 176                  'pub_comment'  => $attr->{'{}PUB_COMMENT'}->{'Value'},
 177              }
 178              );
 179         push @{$self->{'_seendata'}->{'_feats'}}, $f;
 180         # add this feature to the current sequence
 181         $self->{'_seendata'}->{'_seqs'}->[-1]->add_SeqFeature($f);
 182     } elsif( $name eq 'MODEL' ) { # mRNA/transcript
 183         # reset the UTRs
 184         $self->{'_seendata'}->{"five_prime_UTR"}= undef;
 185         $self->{'_seendata'}->{"three_prime_UTR"} = undef;
 186         my ($s,$e) = ($attr->{'{}COORDS'}->{'Value'} =~ /(\d+)\-(\d+)/);
 187         my $strand = 1;
 188         if( $s > $e) {
 189             ($s,$e,$strand) = ( $e,$s,-1);
 190         }
 191         my $parent = $self->{'_seendata'}->{'_feats'}->[-1];
 192         my ($parentid) = $parent->get_tag_values('ID');
 193         my $f = Bio::SeqFeature::Generic->new
 194             (-primary_tag => 'transcript',
 195              -source_tag  => $Default_Source,
 196              -start       => $s,             # we use parent start/stop because 'MODEL' means CDS start/stop
 197              -end         => $e,             # but we want to reflect
 198              -strand      => $strand,
 199              -seq_id      => $seqid,
 200              -tag         => {
 201                  'ID'     => $attr->{'{}FEAT_NAME'}->{'Value'},
 202                  'Parent' => $parentid,
 203                  'Note'   => $attr->{'{}COMMENT'}->{'Value'},
 204              });
 205         $parent->add_SeqFeature($f);
 206         push @{$self->{'_seendata'}->{'_feats'}}, $f;
 207         $self->{'_seendata'}->{'_seqs'}->[-1]->add_SeqFeature($f);
 208     } elsif( $name eq 'EXON' ) { # exon feature
 209         my ($s,$e) = ($attr->{'{}COORDS'}->{'Value'} =~ /(\d+)\-(\d+)/);
 210         my $strand = 1;
 211         if( $s > $e) {
 212             ($s,$e,$strand) = ( $e,$s,-1);
 213         }
 214         my $parent = $self->{'_seendata'}->{'_feats'}->[-1];
 215
 216         my ($parentid) = $parent->get_tag_values('ID');
 217
 218         my $f = Bio::SeqFeature::Generic->new
 219             (-primary_tag => 'exon',
 220              -source_tag  => $Default_Source,
 221              -seq_id      => $seqid,
 222              -start       => $s,
 223              -end         => $e,
 224              -strand      => $strand,
 225              -tag         => {
 226                  'ID'     => $attr->{'{}FEAT_NAME'}->{'Value'},
 227                  'Parent' => $parentid,
 228              });
 229         $parent->add_SeqFeature($f,'EXPAND');
 230         $self->{'_seendata'}->{'_seqs'}->[-1]->add_SeqFeature($f);
 231         # we'll still just add exons to the transcript
 232     } elsif( $name eq 'PROTEIN_SEQ' ) {
 233
 234     } elsif( $name eq 'CDS' ) {
 235         # CDS will be the translation of the transcript
 236         my ($s,$e) = ($attr->{'{}COORDS'}->{'Value'} =~ /(\d+)\-(\d+)/);
 237         my $strand = 1;
 238         if( $s > $e) {
 239             ($s,$e,$strand) = ( $e,$s,-1);
 240         }
 241         my $parent = $self->{'_seendata'}->{'_feats'}->[-1];
 242         my ($parentid) = $parent->get_tag_values('ID');
 243         $self->assert($parent->primary_tag eq 'transcript', 'Testing for primary tag equivalent to mRNA');
 244         $self->assert($parent->strand == $strand || abs($s-$e) == 0, 'Testing that parent feature and current feature strand are equal '. $parentid. ' '.$attr->{'{}FEAT_NAME'}->{'Value'});
 245         my $f = Bio::SeqFeature::Generic->new
 246             (-primary_tag => 'CDS',
 247              -source_tag  => $Default_Source,
 248              -seq_id      => $seqid,
 249              -start       => $s,
 250              -end         => $e,
 251              -strand      => $parent->strand,
 252              -tag         => {
 253                  'ID'     => $attr->{'{}FEAT_NAME'}->{'Value'},
 254                  'Parent' => $parentid, # should be the mRNA
 255              });
 256         $parent->add_SeqFeature($f);
 257         $self->{'_seendata'}->{'_seqs'}->[-1]->add_SeqFeature($f);
 258     } elsif( $name eq 'RNA-EXON' ) {
 259
 260         my ($s,$e) = ($attr->{'{}COORDS'}->{'Value'} =~ /(\d+)\-(\d+)/);
 261         my $strand = 1;
 262         if( $s > $e) {
 263             ($s,$e,$strand) = ( $e,$s,-1);
 264         }
 265         my $parent = $self->{'_seendata'}->{'_feats'}->[-1];
 266         my ($parentid) = $parent->get_tag_values('ID');
 267         my $f = Bio::SeqFeature::Generic->new
 268             (-primary_tag => 'tRNA_exon', # tRNA_exon?
 269              -source_tag  => $Default_Source,
 270              -seq_id      => $seqid,
 271              -start       => $s,
 272              -end         => $e,
 273              -strand      => $strand,
 274              -tag         => {
 275                  'ID'     => $attr->{'{}FEAT_NAME'}->{'Value'},
 276                  'Parent' => $parentid,
 277              }
 278              );
 279         $parent->add_SeqFeature($f);
 280         $self->{'_seendata'}->{'_seqs'}->[-1]->add_SeqFeature($f);
 281     } elsif( $name eq 'PRE-TRNA' ) { # tRNA gene
 282         my ($s,$e) = ( $attr->{'{}COORDS'}->{'Value'} =~/(\d+)\-(\d+)/);
 283         my $strand = 1;
 284         if( $s > $e) {
 285             ($s,$e,$strand) = ( $e,$s,-1);
 286         }
 287         my $f = Bio::SeqFeature::Generic->new
 288             ( -primary_tag => 'tRNA_coding_gene',
 289               -source_tag  => $Default_Source,
 290               -seq_id      => $seqid,
 291               -start       => $s,
 292               -end         => $e,
 293               -strand      => $strand,
 294               -tag         => {'ID' => $attr->{'{}FEAT_NAME'}->{'Value'},
 295                            }
 296               );
 297         push  @{$self->{'_seendata'}->{'_feats'}}, $f;
 298         $self->{'_seendata'}->{'_seqs'}->[-1]->add_SeqFeature($f);
 299     } elsif( $name eq 'TRNA' ) { # tRNA transcript
 300         my ($s,$e) = ($attr->{'{}COORDS'}->{'Value'} =~ /(\d+)\-(\d+)/);
 301         my $strand = 1;
 302         if( $s > $e) {
 303             ($s,$e,$strand) = ( $e,$s,-1);
 304         }
 305         my $parent = $self->{'_seendata'}->{'_feats'}->[-1];
 306         my ($parentid) = $parent->get_tag_values('ID');
 307         my $f = Bio::SeqFeature::Generic->new
 308             (-primary_tag => 'tRNA_primary_transcript',
 309              -source_tag  => $Default_Source,
 310              -start       => $s,
 311              -end         => $e,
 312              -strand      => $strand,
 313              -seq_id      => $seqid,
 314              -tag         => {
 315                  'ID'     => $attr->{'{}FEAT_NAME'}->{'Value'},
 316                  'Parent' => $parentid,
 317                  'Note'   => $attr->{'{}COM_NAME'}->{'Value'},
 318                  'anticodon' => $attr->{'{}ANTICODON'}->{'Value'},
 319                  'pub_locus' => $attr->{'{}PUB_LOCUS'}->{'Value'},
 320
 321              });
 322         $parent->add_SeqFeature($f);
 323         push  @{$self->{'_seendata'}->{'_feats'}}, $f;
 324         $self->{'_seendata'}->{'_seqs'}->[-1]->add_SeqFeature($f);
 325     } elsif( $name eq 'REPEAT_LIST' ) {
 326     } elsif( $name eq 'REPEAT' ) {
 327         my ($s,$e) = ($attr->{'{}COORDS'}->{'Value'} =~ /(\d+)\-(\d+)/);
 328         my $strand = 1;
 329         if( $s > $e) {
 330             ($s,$e,$strand) = ( $e,$s,-1);
 331         }
 332         my $f = Bio::SeqFeature::Generic->new
 333             (-primary_tag => 'simple_repeat',
 334              -source_tag  => $Default_Source,
 335              -seq_id      => $seqid,
 336              -start       => $s,
 337              -end         => $e,
 338              -stand       => $strand,
 339              -tag         => {
 340                  'ID'     => $attr->{'{}FEAT_NAME'}->{'Value'},
 341              });
 342
 343         push @{$self->{'_seendata'}->{'_feats'}}, $f;
 344         $self->{'_seendata'}->{'_seqs'}->[-1]->add_SeqFeature($f);
 345     } elsif ( $name  eq 'AUTHOR' ) {
 346     } elsif( $name eq 'GB_DESCRIPTION' ) {
 347
 348     } elsif( $name eq 'GB_COMMENT' ) {
 349     } elsif( $name eq 'LINEAGE' ) {
 350
 351     } else {
 352         $self->warn("Unknown element $name, ignored\n");
 353     }
 354     push @{$self->{'_state'}}, $name;
 355     $self->SUPER::start_element($ele);
 356 }
 357
 358 sub end_element {
 359     my ($self,$ele) = @_;
 360     pop @{$self->{'_state'}};
 361     my $name = $ele->{'LocalName'};
 362     my $curseq = $self->{'_seendata'}->{'_seqs'}->[-1];
 363     if( $name eq 'AUTHOR_LIST' ) {
 364         if( $curseq->can('annotation') ) {
 365             $curseq->annotation->add_Annotation
 366                 ('reference',Bio::Annotation::Reference->new
 367                  (-authors => join(',',@{$self->{'_seendata'}->{'_authors'}}))
 368                  );
 369         }
 370         $self->{'_seendata'}->{'_authors'} = [];
 371     } elsif( $name eq 'ASSEMBLY' ) {
 372         if( @{$self->{'_seendata'}->{'_feats'} || []} ) {
 373             $self->warn("Leftover features which were not finished!");
 374         }
 375         $self->debug("end element for ASSEMBLY ". $curseq->display_id. "\n");
 376     } elsif( $name eq 'TU' ||
 377              $name eq 'TRNA'  || $name eq 'PRE-TRNA' ||
 378              $name eq 'REPEAT' ) {
 379         pop @{$self->{'_seendata'}->{'_feats'}};
 380     } elsif( $name eq 'MODEL' ) {
 381         # This is all to for adding UTRs
 382
 383         my $model = pop @{$self->{'_seendata'}->{'_feats'}};
 384         my $curseq = $self->{'_seendata'}->{'_seqs'}->[-1];
 385         # sort smallest to largest, don't forget about
 386         # strandedness
 387         my ($parentid) = $model->get_tag_values('Parent');
 388
 389         my @features = $model->get_SeqFeatures();
 390         my @exons = sort { $a->start <=> $b->start }
 391                     grep { $_->primary_tag eq 'exon' } @features;
 392
 393         my @cdsexons = sort { $a->start <=> $b->start }
 394                        grep { $_->primary_tag eq 'CDS' } @features;
 395
 396         # look at the exons, find those which come after the model start
 397         my $cdsexon = shift @cdsexons;
 398         my $exon = shift @exons; # first exon
 399         if( ! defined $cdsexon ) {
 400             $self->warn( "no CDS exons $parentid!");
 401             return;
 402         } elsif( ! defined $exon ) {
 403             $self->warn("no exons $parentid!" );
 404             return;
 405         }
 406         my $utrct = 1;
 407         while( defined $exon && $exon->start < $cdsexon->start ) {
 408             my ($pid) = $exon->get_tag_values('Parent');
 409             $self->debug("LeftPhase: tu-id $parentid mrna-id $pid exon is ".
 410                          $exon->location->to_FTstring.
 411                          " CDSexon is ".$cdsexon->location->to_FTstring."\n");
 412
 413             my $utr = Bio::SeqFeature::Generic->new
 414                 (-seq_id      => $exon->seq_id,
 415                  -strand      => $exon->strand,
 416                  -primary_tag => $exon->strand > 0 ? "five_prime_UTR" : "three_prime_UTR",
 417                  -source_tag  => $Default_Source,
 418                  -tag         => {
 419                      'ID'     => "$pid.UTR".$utrct++,
 420                      'Parent' => $pid },
 421                  );
 422             my ($ns,$ne);
 423             if( $utr->primary_tag eq 'five_prime_UTR' ) {
 424                 $ns = $exon->start;
 425                 $ne = min ( $exon->end, $cdsexon->start - 1);
 426             } else {
 427                 $ne = min( $exon->end, $cdsexon->start - 1);
 428                 $ns = $exon->start;
 429             }
 430             $utr->start($ns); $utr->end($ne);
 431             $model->add_SeqFeature($utr);
 432             $curseq->add_SeqFeature($utr);
 433             $exon = shift @exons;
 434         }
 435         @exons = sort { $a->start <=> $b->start }
 436                  grep {$_->primary_tag eq 'exon' } @features;
 437         @cdsexons = sort { $a->start <=> $b->start }
 438                     grep { $_->primary_tag eq 'CDS' } @features;
 439
 440         $cdsexon = pop @cdsexons;
 441         $exon = pop @exons;
 442         if( ! defined $cdsexon ) {
 443             $self->warn( "no CDS exons $parentid!");
 444             return;
 445         } elsif( ! defined $exon ) {
 446             $self->warn("no exons $parentid!" );
 447             return;
 448         }
 449         $utrct = 1;
 450         while( defined $exon &&$exon->end > $cdsexon->end ) {
 451             my ($pid) = $exon->get_tag_values('Parent');
 452             $self->debug("RightPhase: tu-id $parentid mrna-id $pid exon is ".
 453                          $exon->location->to_FTstring.
 454                          " CDSexon is ".$cdsexon->location->to_FTstring."\n");
 455
 456             my $utr = Bio::SeqFeature::Generic->new
 457                (-seq_id      => $exon->seq_id,
 458                 -strand      => $exon->strand,
 459                 -primary_tag => $exon->strand < 0 ? "five_prime_UTR" : "three_prime_UTR",
 460                 -source_tag  => $Default_Source,
 461                 -tag         => {
 462                     'Parent' => $pid,
 463                     'ID'     => "$pid.UTR".$utrct++,
 464                 }
 465                 );
 466             my ($ns,$ne);
 467             if( $utr->primary_tag eq 'three_prime_UTR' ) {
 468                 $ns = max ( $exon->start, $cdsexon->end + 1);
 469                 $ne = $exon->end;
 470             } else {
 471                 $ns = $cdsexon->end+1;
 472                 $ne = max ( $exon->end, $cdsexon->start + 1);
 473             }
 474             $utr->start($ns); $utr->end($ne);
 475
 476             $model->add_SeqFeature($utr);
 477             $curseq->add_SeqFeature($utr);
 478             $exon = pop @exons;
 479         }
 480     }
 481     $self->SUPER::end_element($ele);
 482 }
 483
 484 sub characters {
 485     my ($self,$data) = @_;
 486     if( ! @{$self->{'_state'}} ) {
 487         $self->warn("Calling characters with no previous start_element call. Ignoring data");
 488     } else {
 489         my $curseq = $self->{'_seendata'}->{'_seqs'}->[-1];
 490         my $curfeat = $self->{'_seendata'}->{'_feats'}->[-1];
 491         my $name = $self->{'_state'}->[-1];
 492         if( defined $curseq ) {
 493             if( $name eq 'CLONE_NAME' ) {
 494                 $self->debug("Clone name is ",$data->{'Data'}, "\n");
 495                 $curseq->display_id($data->{'Data'});
 496             } elsif( $name eq 'ORGANISM' ) {
 497                 my ($genus,$species,$subspec) = split(/\s+/,$data->{Data},3);
 498                 $curseq->species(Bio::Species->new(
 499                                                    -classification =>
 500                                                    [$species,$genus],
 501                                                    -sub_species => $species));
 502             } elsif( $name eq 'LINEAGE' ) {
 503                 $curseq->species->classification(
 504                                   [
 505                                     $curseq->species->species,
 506                                     $curseq->species->genus,
 507                                     reverse  (map { s/^\s+//;
 508                                                     s/\s+$//; $_; }
 509                                               split /[;\.]+/,$data->{'Data'} ),
 510                                     ]
 511                                   );
 512             } elsif( $name eq 'AUTHOR' ) {
 513                 push @{$self->{'_seendata'}->{'_authors'}}, $data->{'Data'};
 514             }
 515         }
 516         if( defined $curfeat ) {
 517             if( $name eq 'EXON' ) { # exon feature
 518             } elsif( $name eq 'RNA-EXON' ) {
 519
 520             } elsif( $name eq 'PROTEIN_SEQ' ) {
 521                 $curfeat->add_tag_value('translation',$data->{'Data'});
 522             } elsif( $name eq 'CDS' ) {
 523             } elsif( $name eq 'PRE-TRNA' ) { # tRNA gene
 524             } elsif( $name eq 'TRNA' ) { # tRNA transcript
 525             } elsif( $name eq 'REPEAT_LIST' ) {
 526             } elsif( $name eq 'REPEAT' ) {
 527                 $curfeat->add_tag_value('Note',$data->{'Data'});
 528             } elsif( $name eq 'GB_COMMENT' ) {
 529                 $curseq->annotation->add_Annotation
 530                     ('comment',
 531                      Bio::Annotation::Comment->new(-text => $data->{'Data'}));
 532             } elsif( $name eq 'GB_DESCRIPTION' ) {
 533                 $curseq->description($data->{'Data'});
 534             }
 535         }
 536     }
 537     $self->SUPER::characters($data);
 538 }
 539
 540
 541 sub assert {
 542     my ($self,$test,$msg) = @_;
 543     $self->throw($msg) unless $test;
 544 }
 545 1;