[Bioperl-l] Bio::SeqIO

[Loerch, Patrick]

I am attempting to connect to a database using DBI, bring in a sequence
using a SQL string, and output that sequence using SeqIO, in Fasta format,
to STDOUT. I keep getting an error message which replies "Bad Name after
Fasta' at fetch_exons.pl Line 14". 

I was wondering if this is a problem with the SeqIO (or possibly with the
databse connection). Here is the script below. There are a number of perl
Modules which I am also using. I am finding this to be very frustrating,
cdoes anyone have any suggestions?

Patrick

#!/usr/local/bin/perl -w

use lib '/usr/lib/perl/;
use strict;
use BigFasta;
use JimKent;
use SQLServer;
use Bio::SeqIO;
use Bio::PrimarySeq;

#my $FASTA = "$ENV{HOME}/superlink_masked.fasta";
my $FASTA = "/info/hgapdb/UCSC/Dec0401/chr22.fa";

my $out = Bio::SeqIO->new(-format => 'Fasta', -fh => exons);

my $accession = $ARGV[0] || "NM_005160";

get_database_handle('OLIVE_UCSC');

my $sql = "select a.* from ucsc_refgene_exons a, ucsc_refgene b where
a.refgene_id = b.refgene_id and b.name like '$accession'";


my @exons = select($sql);


my @transcript = ();
my $strand;


for my $e (@exons) {
    
    my ($start, $end) = jk_coords(start => $e->{exonstart},
				  end   => $e->{exonend});
    
    my $exon = rip_substring(fasta => $FASTA, 
			     start => $start,
			     end   => $end);    
    
    my $seq = Bio::PrimarySeq->new(-seq => $exon,
				   -moltype => 'dna'
				   );
    
    push(@transcript, $seq);
    
    $strand = $e->{strand};
}


my $count = 1;

while (@transcript) {

    if ($strand == -1) {

	my $seq = pop(@transcript);
	$seq->display_id("$accession.$count");
	$out->write_seq($seq->revcom);

    } else {

	my $seq = shift(@transcript);
	$seq->display_id("$accession.$count");
	$out->write_seq($seq);    
    }
    
    $count++;

}

release_database_handle();