Bio::Assembly::IO maq
SummaryIncluded librariesPackage variablesSynopsisDescriptionGeneral documentationMethods
Toolbar
WebCvs
Summary
Bio::Assembly::IO::maq - Driver to read assembly files in maq format *BETA*
Package variables
Privates (from "my" definitions)
$progname = 'maq'
Included modules
Bio::Assembly::Contig
Bio::Assembly::IO
Bio::Assembly::Scaffold
Bio::Assembly::Singlet
Bio::LocatableSeq
Bio::Seq::PrimaryQual
Bio::Seq::Quality
Bio::SeqIO
File::Basename
File::Spec
Inherit
Bio::Assembly::IO
Synopsis
    # convert the native maq map format to plain text
$ maq mapview all.map > all.maq
# Building an input stream use Bio::Assembly::IO; # Assembly loading methods my $asmio = Bio::Assembly::IO->new( -file => 'all.maq', -format => 'maq' ); my $scaffold = $asmio->next_assembly;
Description
This package loads and writes map information in/from maq map files converted by the maq mapview utility. This module is a driver module for
Bio::Assembly::IO input/output.
Parsing is based on Heng Li's description of maq mapview output, found
at the maq manpage: http://maq.sourceforge.net/maq-manpage.shtml.
The basic maq workflow is: map reads to a reference sequence (with
maq map), then create a consensus from the map (with maq
assemble). To read a complete assembly with this module, the
following files need to be available:
 [basename].maq : created by maq mapview [basename].map > [basename].maq
[basename].cns.fastq : created as follows
$ maq assemble [basename].cns [refseq].bfa [basename].map
$ maq cns2fq [basename].cns > [basename].cns.fastq
maq produces only one "contig"; all reads map to the reference
sequence, which covers everything. This module breaks the reads into
contigs by dividing the maq consensus into pieces for which there
are contiguous non-zero quality values.
The module Bio::Tools::Run::Maq will help in this process (eventually).
This module has no write capability. Assemblies are loaded into Bio::Assembly::Scaffold objects composed of
Bio::Assembly::Contig and Bio::Assembly::Singlet objects. Contigs are
not explicitly specified in map files; the division of the map into
contigs is calculated in this module.
Additional assembly information is stored as features. Contig objects have
SeqFeature information associated with the primary_tag:
    _main_contig_feature:$contig_id -> misc contig information
Read objects have sub_seqFeature information associated with the
primary_tag:
    _main_read_feature:$read_id     -> misc read information
Singlets are contigs of a single sequence, as calculated within this module.
They are cataloged separately, as specified in Bio::Assembly::Scaffold.
Methods
next_assemblyDescriptionCode
next_contigDescriptionCode
_init_contigDescriptionCode
_store_contigDescriptionCode
_parse_cns_fileDescriptionCode
_consDescriptionCode
_next_consDescriptionCode
_store_readDescriptionCode
_store_singletDescriptionCode
write_assemblyDescriptionCode
_basenameDescriptionCode
Methods description
next_assemblycode    nextTop
 Title   : next_assembly
Usage : $scaffold = $stream->next_assembly()
Function: return the assembly defined by the map and cns files
Returns : Bio::Assembly::Scaffold object
Args : none
next_contigcodeprevnextTop
 Title   : next_contig
Usage : $scaffold = $stream->next_contig()
Function: Returns the next contig or singlet in the ACE stream.
Returns : a Bio::Assembly::Contig or Bio::Assembly::Single object
Args : none
_init_contig()codeprevnextTop
    Title   : _init_contig
Usage : my $contigobj; $contigobj = $self->_init_contig(
\%contiginfo, $scaffoldobj);
Function: store information of a contig belonging to a scaffold in the
appropriate object
Returns : Bio::Assembly::Contig object
Args : hash, Bio::Assembly::Scaffold
_store_contig()codeprevnextTop
    Title   : _store_contig
Usage : my $contigobj; $contigobj = $self->_store_contig(
\%contiginfo, $contigobj);
Function: store information of a contig belonging to a scaffold
in the appropriate object
Returns : Bio::Assembly::Contig object
Args : hash, Bio::Assembly::Contig
_parse_cns_file()codeprevnextTop
 Title   : _parse_cns_file
Usage : $self->_parse_cns_file
Function: parse the .cns.fastq (consensus) file
associated with the present map;
set the objects cns attribute
Returns : true on success; nil if file dne
Args : none
_cons()codeprevnextTop
 Title   : _cons
Usage : @cons = $self->_cons
Function: get the array of consensus fastq Bio::Seq::Quality objects
Returns : array of Bio::Seq::Quality objects
Args : none
_next_cons()codeprevnextTop
_store_read()codeprevnextTop
    Title   : _store_read
Usage : my $readobj = $self->_store_read(\%readinfo, $contigobj);
Function: store information of a read belonging to a contig
in the appropriate object
Returns : a Bio::LocatableSeq object
Args : hash, Bio::Assembly::Contig
_store_singlet()codeprevnextTop
    Title   : _store_singlet
Usage : my $singletobj = $self->_store_read(\%readinfo, \%contiginfo);
Function: store information of a singlet belonging to a scaffold in a singlet object
Returns : Bio::Assembly::Singlet
Args : hash, hash
write_assembly()codeprevnextTop
    Title   : write_assembly
Usage :
Function: not currently available for maq assemblies
Returns : throw
Args :
_basename()codeprevnextTop
 Title   : _basename
Usage : $self->_basename
Function: return the basename of the associate IO file
Returns : scalar string
Args : none
Methods code
next_assemblydescriptionprevnextTop
sub next_assembly {
    my $self = shift;

    my $assembly = Bio::Assembly::Scaffold->new( -progname => $progname );

    # Load contigs and singlets in the scaffold
while ( my $obj = $self->next_contig()) { # Add contig /singlet to assembly
if ($obj->isa('Bio::Assembly::Singlet')) { # a singlet
$assembly->add_singlet($obj); } else { # a contig
$assembly->add_contig($obj); } } return $assembly;
}
next_contigdescriptionprevnextTop
sub next_contig {
    my $self = shift; # object reference
# Read the file of consensus sequences if it has not already been done for
# this Bio:::Assembly::IO stream already
if (not defined $self->_cons) { $self->_parse_cns_file or $self->throw("Associated maq consensus file is not available"); } # Contig and read related
my $contigobj; my %contiginfo; # Loop over all assembly file lines
while ($_ = $self->_readline) { chomp; next if /^$/; # mapview format parsing ; every line is a read...
my %readinfo; @readinfo{ qw(read_name chr posn strand insert_size paired_flag map_qual se_map_qual alt_map_qual num_mm_best_hit sum_qual_mm_best_hit zero_mm_hits one_mm_hits read_len seqstr qualstr) } = split(/\s+/); # sanger conversion
my @qual = map { ord($_)-33 } split('', $readinfo{qualstr}); $readinfo{seq} = Bio::Seq::Quality->new( -id => $readinfo{read_name}, -seq => $readinfo{seqstr}, -qual =>\@ qual ); if ( not defined $contiginfo{start} ) { # First read of new contig or singlet
$contiginfo{'seqnum'} = 1; $contiginfo{'qualobj'} = $self->_next_cons; $contiginfo{'start'} = $contiginfo{'qualobj'}->start; $contiginfo{'end'} = $contiginfo{'qualobj'}->end; $contiginfo{'asmbl_id'} = 'maq_assy['.$self->_basename.']/'.$contiginfo{start}.'-'.$contiginfo{end}; # It may be a singlet, but assume it's a contig for now
$contigobj = $self->_init_contig(\%contiginfo); $self->_store_read(\%readinfo, $contigobj); } else { if ( $readinfo{'posn'} <= $contiginfo{end} ) { # Add read to existing contig
$contiginfo{'seqnum'}++; $self->_store_read(\%readinfo, $contigobj); } else { # Read belongs in a new contig
if ($contiginfo{'seqnum'} > 1) { $self->_store_contig(\%contiginfo, $contigobj); } else { # singlet
# Create a new singlet object from the read info
$contigobj = $self->_store_singlet(\%contiginfo, $contigobj); } # do a pushback
$self->_pushback($_); last; } } } return $contigobj;
}
_init_contigdescriptionprevnextTop
sub _init_contig {
    my ($self, $contiginfo) = @_;
    # Create a contig and attach it to scaffold
my $contigobj = Bio::Assembly::Contig->new( -id => $$contiginfo{'asmbl_id'}, -source => $progname, -strand => 1 ); return $contigobj;
}
_store_contigdescriptionprevnextTop
sub _store_contig {
    my ($self, $contiginfo, $contigobj) = @_;

    $self->throw("Contig object must be defined") unless $contigobj;

    my $consensus_seq = Bio::LocatableSeq->new(
        -id    => $$contiginfo{'asmbl_id'},
        -seq   => $$contiginfo{'qualobj'}->seq,
        -start => 1,
    );
    $contigobj->set_consensus_sequence($consensus_seq);
    my $consensus_qual = Bio::Seq::PrimaryQual->new(
        -id    => $$contiginfo{'asmbl_id'},
        -qual  => $$contiginfo{'qualobj'}->qual,
        -start => 1,
    );
    $contigobj->set_consensus_quality($consensus_qual);

    # Add other misc contig information as features of the contig
# Add other misc read information as subsequence feature
my @other = grep !/asmbl_id|end|qualobj|start/, keys %$contiginfo; my %other; @other{@other} = @$contiginfo{@other}; my $contigtags = Bio::SeqFeature::Generic->new( -primary => '_main_contig_feature', -source => $$contiginfo{'asmbl_id'}, -start => 1, -end => $contigobj->get_consensus_length(), -strand => 1, # dumping ground:
-tag =>\% other ); $contigobj->add_features([ $contigtags ], 1); return $contigobj;
}
_parse_cns_filedescriptionprevnextTop
sub _parse_cns_file {
    my ($self) = @_;
    my @cons;

    $self->{'_cns_parsed'} = 1;
    my $file = $self->file;
    $file =~ s/^[<>+]*//; # byebye parasitic mode chars
my ($fname, $dir, $suf) = fileparse($file, ".maq"); my $cnsf = File::Spec->catdir($dir, "$fname.cns.fastq"); return unless (-e $cnsf ); my $fqio = Bio::SeqIO->new( -file => $cnsf ); my $cns = $fqio->next_seq; # now, infer the contigs on the basis of quality values
# - assuming quality of zero => no coverage
my $qual = $cns->qual; # covered sites
my @sites = grep { $$qual[$_] > 0 } (0..$#$qual); my @ranges = ($sites[0]+1); for my $i (1..$#sites) { if ($sites[$i]-$sites[$i-1]>1) { push @ranges, $sites[$i-1]+1, $sites[$i]+1; } } push @ranges, $sites[-1]; for (my $i = 0; $i<$#ranges; $i+=2) { push @cons, Bio::Seq::Quality->new( -display_id => "${fname}/".$ranges[$i]."-".$ranges[$i+1], -start => $ranges[$i], -end => $ranges[$i+1], -seq => $cns->subseq($ranges[$i], $ranges[$i+1]), -qual => [@{$cns->qual}[$ranges[$i]-1..$ranges[$i+1]-1]] ); } $self->{'_cons'} =\@ cons; return 1;
}
_consdescriptionprevnextTop
sub _cons {
    my $self = shift;
    my $cons = undef;
    if (defined $self->{'_cons'}) {
      $cons = $self->{'_cons'};
    }
    return $cons;
}
_next_consdescriptionprevnextTop
sub _next_cons {
 shift(@{shift->{'_cons'}})
}
_store_readdescriptionprevnextTop
sub _store_read {
   my ($self, $readinfo, $contigobj) = @_;

   # Create an aligned read object
$$readinfo{'strand'} = ($$readinfo{strand} eq '+' ? 1 : -1); my $readobj = Bio::LocatableSeq->new( -display_id => $$readinfo{'read_name'}, -primary_id => $$readinfo{'read_name'}, -seq => $$readinfo{'seqstr'}, -start => 1, -strand => $$readinfo{'strand'}, -alphabet => 'dna' ); # Add read location and sequence to contig (in 'gapped consensus' coordinates)
$$readinfo{'aln_start'} = $$readinfo{'posn'}; $$readinfo{'aln_end'} = $$readinfo{'posn'} + length($$readinfo{'seqstr'})-1; my $alncoord = Bio::SeqFeature::Generic->new( -primary => $readobj->id, -start => $$readinfo{'aln_start'}, -end => $$readinfo{'aln_end'}, -strand => $$readinfo{'strand'}, -qual => join(' ', $$readinfo{seq}->qual), # check here, need to create contigs "by hand"...
-tag => { 'contig' => $contigobj->id() } ); $contigobj->set_seq_coord($alncoord, $readobj); # Add other misc read information as subsequence feature
my @other = grep !/aln_(?:end|start)|seq(?:str)?|strand/, keys %$readinfo; my %other; @other{@other} = @$readinfo{@other}; my $readtags = Bio::SeqFeature::Generic->new( -primary => '_main_read_feature', -source => $readobj->id, -start => $$readinfo{'aln_start'}, -end => $$readinfo{'aln_end'}, -strand => $$readinfo{'strand'}, # dumping ground:
-tag =>\% other ); $contigobj->get_features_collection->add_features([$readtags]); $contigobj->get_features_collection->add_SeqFeature($alncoord, $readtags); return $readobj; } #### revamp for maq
}
_store_singletdescriptionprevnextTop
sub _store_singlet {
    my ($self, $contiginfo, $contigobj) = @_;

    my $contigid = $$contiginfo{'asmbl_id'};
    my $seqref = ($contigobj->each_seq())[0];
    my $singletobj = Bio::Assembly::Singlet->new( -id     => $contigid,
                                                  -seqref => $seqref   );

    # Add other misc contig information as features of the contig
# Add other misc read information as subsequence feature
#my @other = grep !/_sfc|_assembly|_elem/, keys %$contiginfo; # remove the objects; _elem contains a code ref and can't be frozen. Just shooting blind here.
#my %other;
#@other{@other} = @$contiginfo{@other};
#my $contigtags = Bio::SeqFeature::Generic->new(
# -primary => '_main_contig_feature',
# -source => $$contiginfo{asmbl_id},
# -start => 1,
# -end => $singletobj->get_consensus_length(),
# -strand => 1,
# # dumping ground:
# -tag => \%other
#);
#$singletobj->add_features([ $contigtags ], 1);
#$$readinfo{'aln_start'} = $$readinfo{'start'};
#$$readinfo{'aln_end'} = $$readinfo{'end'};
#$$readinfo{'strand'} = ($$readinfo{strand} eq '+' ? 1 : -1);
#my $alncoord = Bio::SeqFeature::Generic->new(
# -primary => '_aligned_coord',
# -source => $$readinfo{read_name},
# -start => $$readinfo{'start'},
# -end => $$readinfo{'end'},
# -strand => $$readinfo{'strand'},
# -tag => { 'contig' => $$contiginfo{asmbl_id} }
# );
#$alncoord->attach_seq($singletobj->seqref);
#$singletobj->add_features([ $alncoord ], 0);
# Add other misc read information as subsequence feature
#my @other = grep !/seqstr|strand/, keys %$readinfo;
#my %other;
#@other{@other} = @$readinfo{@other};
#my $readtags = Bio::SeqFeature::Generic->new(
# -primary => '_main_read_feature',
# -source => $$readinfo{read_name},
# -start => $$readinfo{'aln_start'},
# -end => $$readinfo{'aln_end'},
# -strand => $$readinfo{'strand'},
# # dumping ground:
# -tag => \%other
# );
#$singletobj->get_features_collection->add_features([$readtags]);
#$singletobj->get_features_collection->add_SeqFeature($alncoord, $readtags);
return $singletobj; } ###### writes -- need them??
}
write_assemblydescriptionprevnextTop
sub write_assembly {
    my ($self,@args) = @_;
    $self->throw("Writes not currently available for maq assemblies. Complain to author.")
}
_basenamedescriptionprevnextTop
sub _basename {
    my $self = shift;
    return (fileparse($self->file, ".maq"))[0];
}

1;
__END__
}
General documentation
TODOTop
    *(1)
    Add pod descriptions of maq descriptive data (currently SeqFeatures
added to each contig component)
    *(2)
    Add features describing the aggregate status of reads and contigs based
on the maq "paired flag"
FEEDBACKTop
Mailing ListsTop
User feedback is an integral part of the evolution of this and other
Bioperl modules. Send your comments and suggestions preferably to the
Bioperl mailing lists Your participation is much appreciated.
  bioperl-l@bioperl.org                  - General discussion
http://bioperl.org/wiki/Mailing_lists - About the mailing lists
Support Top
Please direct usage questions or support issues to the mailing list:
bioperl-l@bioperl.org
rather than to the module maintainer directly. Many experienced and
reponsive experts will be able look at the problem and quickly
address it. Please include a thorough description of the problem
with code and data examples if at all possible.
Reporting BugsTop
Report bugs to the BioPerl bug tracking system to help us keep track
the bugs and their resolution. Bug reports can be submitted via email
or the web:
  bioperl-bugs@bio.perl.org
https://redmine.open-bio.org/projects/bioperl/
AUTHOR - Mark A. JensenTop
Email maj -at- fortinbras -dot- us
CONTRIBUTORSTop
Further improvements by Florent Angly
(florent dot angly at gmail dot com)
ACKNOWLEDGEMENTTop
Code and some POD text ripped liberally from Florent Angly's
Bio::Assembly::IO::tigr.
APPENDIXTop
The rest of the documentation details each of the object
methods. Internal methods are usually preceded with a "_".