[Bioperl-l] Errors using Bio::Assembly

[anunberg at oriongenomics.com]

Thanks for the quick response, I will be at the OReilly conference next 
week so no rush.  I would add the bit about the qual files in the docs.  I 
was wondering how you were calling high qual discrenpancies, and i thought 
it was just the default upper/lower case you were using in the trace 
sequences in the ace file.  Also, the high_quality_discrepancies returns 
an array of Bio::SeqFeatures::Generic objects yes? but in the web 
documentation it says a SeqFeature Collection object.

My ace files are not in the context of phredphrap. I only save the ace 
files themselves.  I think I will attempt to right a method that will do 
high qual discrepancies based on the ace file itself.

Thanks again, this does help alot
Andy
 On Fri, 31 Jan 2003, 
Robson Francisco de Souza wrote:

> 
> 	Hi Andy,
> 
> 	You are trying to use a Bio::Assembly::Contig method that depends
> on quality values for the read sequences, but you haven't the appropriate
> quality objects first. An ACE file has no data about read base qualities,
> those are loaded by consed from the phd file in ../phd_dir. To use
> high_quality_discrepancies from Bio::Assembly::ContigAnalysis, you must
> first add quality objects using set_seq_qual. See my changes to you code
> below. See my comments on the other errors after the code.
> 
> use lib '/loginhome/anunberg/bioperl-1.2';
> use strict;
> use Bio::Assembly::IO;
> use Bio::Assembly::ContigAnalysis;
> use Carp qw(cluck);
> my $in=Bio::Assembly::IO->new(-file=>$ARGV[0], -format=>'ace');
> my $assembly=$in->next_assembly;
> my @contigs=$assembly->all_contigs;
> 
> my $PhdDir = '../phd_dir'; # change #1: phd directory
> 
> foreach my $contig(@contigs){
> 	my $analysis=Bio::Assembly::ContigAnalysis->new(-contig=>$contig);
> 
> 	# change #2: this block sets quality values for reads
> 	#
> 	foreach my $seqID ($contig->get_seq_ids) {
> 	  my $swq = Bio::SeqIO->new(-file=>"$PhdDir/${seqID}.phd.1",
>                       		    -format=>'phd');
> 	  $swq = $swq->next_seq();
>           my $seq = $contig->get_seq_by_name($seqID);
> 	  $contig->set_seq_qual($seq,$swq);
>         }
> 	# change #2: end
> 
> 	my @hqd=$analysis->high_quality_discrepancies;
> #	my @ss=$analysis->single_strand; # does not work. See below.
> 	my @lcd=$analysis->low_consensus_quality;
> 	my $id=$contig->id;
> 	print "$id\t";
> 	print scalar @hqd."\t";
> 	print scalar @lcd."\n";
> }
> 
> > >My first error is when calling the next_assembly method.  Some parsing
> > >error occurs, I think.  I got mutlitiple warnings all the same
> > >Argument "" isn't numeric in numeric comparison (<=>) at
> > >/loginhome/anunberg/bioperl-1.2/Bio/SeqFeature/Collection.pm line 436,
> > ><GEN0> line 54.
> > >the lines correspond to the QA tag lines in the ace file.  Sometimes the
> > >same complaint occurs one line above it if its blank, but not all blank
> > >lines above the QA tag produce this warning
> > >when using the method high_quality_discrepancies, i get this warning
> > >Argument "" isn't numeric in numeric comparison (<=>) at
> > >/loginhome/anunberg/bioperl-1.2/Bio/SeqFeature/Collection.pm line 436,
> > ><GEN0> line 358.
> > >which is a blank line at the end of the ace file.
> 
> 	This one I could not repreduce. The example in t/data is working
> on your system? Maybe you can send me an sample file latter, so that I can
> reproduce your error, but the errors refer to a sort method in
> Bio::SeqFeature::Collection. It looks like the module did not load the
> alignment coordinate for your reads...
> 
> > >I may be using the method wrong, but using an ace file where i know to have
> > >a high quality descrepancie, i still get 0 .
> > >The documentation is a little off for this method.  The description says
> > >that a Bio::SeqFeature::Collection is returned, but the code and the
> > >example has an array of Bio::SeqFeature::Generic objects being returned(if
> > >i am reading the code correctly).
> 
> 	Ok. My patch to your code fix this one. Maybe I should add this
> example to the docs?
> 
> > >the single_strand method gives me the following exception
> > >------------- EXCEPTION Bio::Root::NotImplemented -------------
> > >MSG: Abstract method
> > >"Bio::Assembly::ContigAnalysis::_merge_overlapping_features" is not
> > >implemented by package Bio::Assembly::ContigAnalysis.
> > >This is not your fault - author of Bio::Assembly::ContigAnalysis should be
> > >blamed!
> > >
> > >STACK Bio::Root::RootI::throw_not_implemented
> > >/loginhome/anunberg/bioperl-1.2/Bio/Root/RootI.pm:529
> > >STACK Bio::Assembly::ContigAnalysis::_merge_overlapping_features
> > >/loginhome/anunberg/bioperl-1.2/Bio/Assembly/ContigAnalysis.pm:476
> > >STACK Bio::Assembly::ContigAnalysis::single_strand
> > >/loginhome/anunberg/bioperl-1.2/Bio/Assembly/ContigAnalysis.pm:442
> > >STACK toplevel ../testing_BioAssembly.pl:31
> 
> 	As the error says: it is my fault and my only. I did not had time
> to complete the single_strand method yet. This method needs a complete
> rewrite. As I'm quite a lot busy wright now, I'll need a few weeks to
> finish it.
> 
> > >Thanks to anyone who can help.
> > >Andy
> 
> 	Hope that helps.
> 	Cheers,
> 			Robson
> 
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at bioperl.org
> http://bioperl.org/mailman/listinfo/bioperl-l
>