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