[Bioperl-l] Errors using Bio::Assembly

[Robson Francisco de Souza]

	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