view release on metacpan or  search on metacpan

Bio/Assembly/IO/tigr.pm  view on Meta::CPAN

                          'cds_end'    => $$contiginfo{'cds_end'},
                          'ed_pn'      => $$contiginfo{'ed_pn'},
                          'ed_date'    => $$contiginfo{'ed_date'},
                          'comment'    => $$contiginfo{'comment'},
                          'frameshift' => $$contiginfo{'frameshift'} }
   );
   $singletobj->add_features([ $contigtags ], 1);

   # Add read location and sequence to singlet features (in 'gapped consensus' coordinates)
   $$readinfo{'aln_start'} = $$readinfo{'offset'} + 1; # seq offset is in gapped coordinates
   $$readinfo{'aln_end'} = $$readinfo{'aln_start'} + length($$readinfo{'lsequence'}) - 1; # lsequence is aligned seq

   my $alncoord = Bio::SeqFeature::Generic->new(
       -primary     => '_aligned_coord',
       -source      => $readid,
       -start       => $$readinfo{'aln_start'},
       -end         => $$readinfo{'aln_end'},
       -strand      => $$readinfo{'strand'},
       -tag         => { 'contig' => $contigid }
   );
   $alncoord->attach_seq($singletobj->seqref);
   $singletobj->add_features([ $alncoord ], 0);

   # Add quality clipping read information in singlet features
   # (from 'aligned read' to 'gapped consensus' coordinates)
   $$readinfo{'clip_start'} = $$readinfo{'seq_lend'};
   $$readinfo{'clip_end'}   = $$readinfo{'seq_rend'};
   my $clipcoord = Bio::SeqFeature::Generic->new(
       -primary     => '_quality_clipping',
       -source      => $readid,
       -start       => $$readinfo{'clip_start'},
       -end         => $$readinfo{'clip_end'},
       -strand      => $$readinfo{'strand'},
       -tag         => { 'contig' => $contigid }
   );
   $clipcoord->attach_seq($singletobj->seqref);
   $singletobj->add_features([ $clipcoord ], 0);

   # Add other misc read information as subsequence feature
   my $readtags = Bio::SeqFeature::Generic->new(
       -primary     => '_main_read_feature',
       -source      => $readid,
       -start       => $$readinfo{'aln_start'},
       -end         => $$readinfo{'aln_end'},
       -strand      => $$readinfo{'strand'},
       -tag         => { 'best'    => $$readinfo{'best'},
                         'comment' => $$readinfo{'comment'} }
   );
   $singletobj->get_features_collection->add_features([$readtags]);
   $singletobj->get_features_collection->add_SeqFeature($alncoord, $readtags);

   return $singletobj;
}


=head2 write_assembly

    Title   : write_assembly
    Usage   : $asmio->write_assembly($assembly)
    Function: Write the assembly object in TIGR Assembler compatible format. The
              contig IDs are sorted naturally if the Sort::Naturally module is
              present, or lexically otherwise. Internally, write_assembly use
              the write_contig, write_footer and write_header methods. Use these
              methods if you want more control on the writing process.
    Returns : 1 on success, 0 for error
    Args    : A Bio::Assembly::Scaffold object
              1 to write singlets in the assembly file, 0 otherwise

=cut


=head2 write_contig

    Title   : write_contig
    Usage   : $asmio->write_contig($contig)
    Function: Write a contig or singlet object in TIGR compatible format. Quality
              scores are automatically generated if the contig does not contain
              any
    Returns : 1 on success, 0 for error
    Args    : A Bio::Assembly::Contig or Singlet object

=cut

sub write_contig {
    my ($self, @args) = @_;
    my ($contigobj) = $self->_rearrange([qw(CONTIG)], @args);

    # Sanity check
    if ( !$contigobj || !$contigobj->isa('Bio::Assembly::Contig') ) {
        $self->throw("Must provide a Bio::Assembly::Contig or Singlet object when calling write_contig");
    }

    my $decimal_format = '%.2f';
    my $contigid = $contigobj->id;
    my $numseqs = $contigobj->num_sequences;

    if ( $contigobj->isa('Bio::Assembly::Singlet') ) {
        # This is a singlet
        my $readid     = $contigobj->seqref->id;      
        my $singletobj = $contigobj;

        # Get contig information
        my ($contanno) = $singletobj->get_features_collection->get_features_by_type("_main_contig_feature:$contigid");

        my %contiginfo;
        $contiginfo{'sequence'}   = $singletobj->seqref->seq;
        $contiginfo{'lsequence'}  = $contiginfo{'sequence'};
        $contiginfo{'quality'}    = $self->_qual_dec2hex(
            join ' ', @{$singletobj->get_consensus_quality->qual} );
        $contiginfo{'asmbl_id'}   = $contigid;
        $contiginfo{'seq_id'}     = ($contanno->get_tag_values('seq_id'))[0];   
        $contiginfo{'com_name'}   = ($contanno->get_tag_values('com_name'))[0];
        $contiginfo{'type'}       = ($contanno->get_tag_values('type'))[0];
        $contiginfo{'method'}     = ($contanno->get_tag_values('method'))[0];
        $contiginfo{'ed_status'}  = ($contanno->get_tag_values('ed_status'))[0];
        $contiginfo{'redundancy'} = sprintf($decimal_format, 1);
        $contiginfo{'perc_N'}     = sprintf(
            $decimal_format, $self->_perc_N($contiginfo{'sequence'}));
        $contiginfo{'seqnum'}     = 1;
        $contiginfo{'full_cds'}   = ($contanno->get_tag_values('full_cds'))[0];
        $contiginfo{'cds_start'}  = ($contanno->get_tag_values('cds_start'))[0];