# Create the Lucy object from an existing Lucy output file 
  @params = ('seqfile' => 'lucy.seq', 'lucy_verbose' => 1);
  $lucyObj = Bio::Tools::Lucy->new(@params);

  # Get names of all sequences
  $names = $lucyObj->get_sequence_names();

  #  Print seq and qual values for sequences >400 bp in order to run CAP3
  foreach $name (@$names) {
      next unless $lucyObj->length_clear($name) > 400;
      print SEQ ">$name\n", $lucyObj->sequence($name), "\n";
      print QUAL ">$name\n", $lucyObj->quality($name), "\n";
  }

  # Get an array of Bio::PrimarySeq objects
  @seqObjs = $lucyObj->get_Seq_Objs();

  @params = ('seqfile' => 'lucy.seq', 'adv_stderr' => 1, 
	     'fwd_desig' => '_F', 'rev_desig' => '_R');
  $lucyObj = Bio::Tools::Lucy->new(@params);

  You should get an array with the sequence names in order to use
  accessor methods.  Note: The Lucy binary program will fail unless
  the sequence names provided as input are unique.

  $names_ref = $lucyObj->get_sequence_names();

  This code snippet will produce a Fasta format file with sequence
  lengths and %GC in the description line.

  foreach $name (@$names) {
      print FILE ">$name\t",
		 $lucyObj->length_clear($name), "\t",
		 $lucyObj->per_GC($name), "\n",
		 $lucyObj->sequence($name), "\n";  
  }

  Print seq and qual values for sequences >400 bp in order to assemble
  them with CAP3 (or other assembler).

  foreach $name (@$names) {
      next unless $lucyObj->length_clear($name) > 400;
      print SEQ ">$name\n", $lucyObj->sequence($name), "\n";
      print QUAL ">$name\n", $lucyObj->quality($name), "\n";
  }

  Get all the sequences as Bio::PrimarySeq objects (eg., for use with
  Bio::Tools::Blast to perform BLAST).

  @seqObjs = $lucyObj->get_Seq_Objs();

  Or use only those sequences that are full length and have a Poly-A
  tail.

  foreach $name (@$names) {
      next unless ($lucyObj->full_length($name) and $lucy->polyA($name));
      push @seqObjs, $lucyObj->get_Seq_Obj($name);
  }

  Get the names of those sequences that were rejected by Lucy.

  $rejects_ref = $lucyObj->get_rejects();

  Print the names of the rejects and 1 letter code for reason they
  were rejected.

  foreach $key (sort keys %$rejects_ref) {
      print "$key:  ", $rejects_ref->{$key};
  }

  There is a lot of other information available about the sequences
  analyzed by Lucy (see APPENDIX).  This module can be used with the
  DBI module to store this sequence information in a database.

    my $self = shift;
    my $attr = $AUTOLOAD;
    $attr =~ s/.*:://;
    $attr = lc $attr;
    $self->throw("Unallowed parameter: $attr !") unless $OK_FIELD{$attr};
    $self->{$attr} = shift if @_;
    return $self->{$attr};

    my ($class,@args) = @_;
    my $self = $class->SUPER::new(@args);
    my ($attr, $value);
    while (@args) {
	$attr = shift @args;
	$attr = lc $attr;
	$value = shift @args;
	$self->{$attr} = $value;
    }
    &_parse($self);	
    return $self;

    my $self = shift;
    $self->{seqfile} =~ /^(\S+)\.\S+$/;
    my $file = $1;
    
    print "Opening $self->{seqfile} for parsing...\n" if $self->{lucy_verbose};
    open SEQ, "$self->{seqfile}" or $self->throw("Could not open sequence file: $self->{seqfile}");
    my ($name, $line);
    my $seq = "";
    my @lines = <SEQ>;
    while ($line = pop @lines) {
	chomp $line;
	if ($line =~ /^>(\S+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)/) {    
            $name = $1;
	    if ($self->{fwd_desig}) {
		$self->{sequences}{$name}{direction} = "F" if $name =~ /^(\S+)($self->{fwd_desig})$/;
	    }
	    if ($self->{rev_desig}) {
                $self->{sequences}{$name}{direction} = "R" if $name =~ /^(\S+)($self->{rev_desig})$/;
            }
	    $self->{sequences}{$name}{min_clone_len} = $2; # this is used for TIGR Assembler, as are $3 and $4
	    $self->{sequences}{$name}{max_clone_len} = $3;
	    $self->{sequences}{$name}{med_clone_len} = $4; 
	    $self->{sequences}{$name}{beg_clear} = $5;
	    $self->{sequences}{$name}{end_clear} = $6;
	    $self->{sequences}{$name}{length_raw} = $seq =~ tr/[AGCTN]//; # from what I've seen, these are the bases Phred calls.  Please let me know if I'm wrong.     
	    my $beg = $5-1; # substr function begins with index 0
	    $seq = $self->{sequences}{$name}{sequence} = substr ($seq, $beg, $6-$beg);
	    my $count = $self->{sequences}{$name}{length_clear} = $seq =~ tr/[AGCTN]//;
	    my $countGC =  $seq =~ tr/[GC]//;
	    $self->{sequences}{$name}{per_GC} = $countGC/$count * 100;
	    $seq = "";
	}
	else {
	    $seq = $line.$seq;
	} 
    }


# now parse quality values (check for presence of quality file first) 
    if ($self->{qualfile}) {
	open QUAL, "$self->{qualfile}" or $self->throw("Could not open quality file: $self->{qualfile}");
	@lines = <QUAL>;
    }
    elsif (-e "$file.qual") {
	print "You did not set qualfile, but I'm opening $file.qual\n" if $self->{lucy_verbose};
	$self->qualfile("$file.qual");
	open QUAL, "$file.qual" or $self->throw("Could not open quality file: $file.qual");
	@lines = <QUAL>;
    }
    else {
	print "I did not find a quality file.  You will not be able to use all of the accessor methods.\n" if $self->{lucy_verbose};
	@lines = ();
    }

    my (@vals, @slice, $num, $tot, $vals);  
    my $qual = ""; 
    while ($line = pop @lines) {
	chomp $line;
 	if ($line =~ /^>(\S+)/) {
	    $name = $1;
	    @vals = split /\s/ , $qual;
	    @slice = @vals[$self->{sequences}{$name}{beg_clear} .. $self->{sequences}{$name}{end_clear}];
	    $vals = join "\t", @slice;
	    $self->{sequences}{$name}{quality} = $vals;
	    $qual = "";
	    foreach $num (@slice) {
		$tot += $num;
	    }
	    $num = @slice;
	    $self->{sequences}{$name}{avg_quality} = $tot/$num;
	    $tot = 0;
	}
	else {
	    $qual = $line.$qual;
	}
    }

# determine whether reads are full length

    if ($self->{infofile}) {
        open INFO, "$self->{infofile}" or $self->throw("Could not open info file: $self->{infofile}");
	@lines = <INFO>;
    }
    elsif (-e "$file.info") {
	print "You did not set infofile, but I'm opening $file.info\n" if $self->{lucy_verbose};
	$self->infofile("$file.info");
        open INFO, "$file.info" or $self->throw("Could not open info file: $file.info");
	@lines = <INFO>;
    }
    else {
	print "I did not find an info file.  You will not be able to use all of the accessor methods.\n" if $self->{lucy_verbose};
	@lines = ();
    }

    foreach (@lines) {
	/^(\S+).+CLV\s+(\d+)\s+(\d+)$/;
	if ($2>0 && $3>0) {
	    $self->{sequences}{$1}{full_length} = 1 if $self->{sequences}{$1}; # will show cleavage info for rejected sequences too
	}
    }


# parse rejects (and presence of poly-A if Lucy has been modified)

    if ($self->{stderrfile}) {
        open STDERR_LUCY, "$self->{stderrfile}" or $self->throw("Could not open quality file: $self->{stderrfile}");
	@lines = <STDERR_LUCY>;

    }
    elsif (-e "$file.stderr") {
	print "You did not set stderrfile, but I'm opening $file.stderr\n" if $self->{lucy_verbose};
	$self->stderrfile("$file.stderr");
        open STDERR_LUCY, "$file.stderr" or $self->throw("Could not open quality file: $file.stderr");
	@lines = <STDERR_LUCY>;
    }
    else {
	print "I did not find a standard error file.  You will not be able to use all of the accessor methods.\n" if $self->{lucy_verbose};
        @lines = ();
    }

    if ($self->{adv_stderr}) {
	foreach (@lines) {
	    $self->{reject}{$1} = "Q" if /dropping\s+(\S+)/;
	    $self->{reject}{$1} = "V" if /Vector: (\S+)/;
	    $self->{reject}{$1} = "E" if /Empty: (\S+)/;
	    $self->{reject}{$1} = "S" if /Short: (\S+)/;
	    $self->{sequences}{$1}{polyA} = 1 if /(\S+) has PolyA/;
	    if (/Dropped PolyA: (\S+)/) {
		$self->{reject}{$1} = "P";		
		delete $self->{sequences}{$1};
	    }
	}
    }
    else {
	foreach (@lines) {
	    $self->{reject}{$1} = "R" if /dropping\s+(\S+)/;
	}
    }

    my $self = shift;
    my($seqobj, @seqobjs);
    foreach my $key (sort keys %{$self->{sequences}}) {
	$seqobj = Bio::PrimarySeq->new( -seq => "$self->{sequences}{$key}{sequence}",
					-id => "$key");
	push @seqobjs, $seqobj;
    }
    return\@ seqobjs;

    my ($self, $key) = @_;
    my $seqobj = Bio::PrimarySeq->new( -seq => "$self->{sequences}{$key}{sequence}",
                                    -id => "$key");
    return $seqobj;

    my $self = shift;
    my @keys = sort keys %{$self->{sequences}};
    return\@ keys;

    my ($self, $key) = @_;
    return $self->{sequences}{$key}{sequence};

    my($self, $key) = @_;
    return $self->{sequences}{$key}{quality};

    my($self, $key) = @_;
    return $self->{sequences}{$key}{avg_quality};

    my($self, $key) = @_;
    return $self->{sequences}{$key}{direction} if $self->{sequences}{$key}{direction}; 
    return "";

    my($self, $key) = @_;
    return $self->{sequences}{$key}{length_raw};

    my($self, $key) = @_;
    return $self->{sequences}{$key}{length_clear};

    my($self, $key) = @_;
    return $self->{sequences}{$key}{beg_clear};

    my($self, $key) = @_;
    return $self->{sequences}{$key}{end_clear};

    my($self, $key) = @_;
    return $self->{sequences}{$key}{per_GC};

    my($self, $key) = @_;
    return 1 if $self->{sequences}{$key}{full_length};
    return 0;

    my($self, $key) = @_;
    return 1 if $self->{sequences}{$key}{polyA};
    return 0;

    my $self = shift;
    return $self->{reject};

    bioperl-l@bioperl.org             - General discussion
    http://bio.perl.org/MailList.html - About the mailing lists

    bioperl-bugs@bio.perl.org
    http://bugzilla.bioperl.org/

  352a353,354
  >       /* AGW added next line */
  >       fprintf(stderr, "Empty: %s\n", seqs[i].name);
  639a642,643
  > 	    /* AGW added next line */
  > 	    fprintf(stderr, "Short/ no insert: %s\n", seqs[i].name);
  678c682,686
  < 	if (left) seqs[i].left+=left;
  ---
  > 	if (left) {
  > 	  seqs[i].left+=left;
  > 	  /*  AGW added next line */
  > 	  fprintf(stderr, "%s has PolyA (left).\n", seqs[i].name);
  > 	}
  681c689,693
  < 	if (right) seqs[i].right-=right;
  ---
  > 	if (right) {
  > 	  seqs[i].right-=right;
  > 	  /* AGW added next line */
  > 	  fprintf(stderr, "%s has PolyA (right).\n", seqs[i].name);
  > 	}
  682a695,696
  > 	  /* AGW added next line */
  > 	  fprintf(stderr, "Dropped PolyA: %s\n", seqs[i].name);	
  734a749,750
  > 	  /* AGW added next line */
  > 	  fprintf(stderr, "Vector: %s\n", seqs[i].name);