[Bioperl-l] fasta2gff_landmark ...

[Dan Bolser]

Hi,

I am trying to create a script to write one 'GFF landmark' line for
each entry in a fasta file. I'm doing this so that the landmarks show
up in GBrowse (along with my attached similarity 'features').

In a nutshell, I was told to add lines like:

C08HBa0216M19.1    .    contig    1    <length>    .    .    .
ID=C08HBa0216M19.1;Name=C08HBa0216M19.1


to the GFF so that GBrowse would pick up the features which have
C08HBa0216M19.1 as their 'landmark'. This seems to be working great
for the sequences where I added this line manually, but now I want to
write a script to automatically generate these GFF landmarks from a
fasta file.


So far I have a script that looks like this:

<source lang="perl">
#!/usr/bin/perl -w

use strict;

use Bio::SeqIO;
use Bio::SeqFeature::Generic;
use Bio::FeatureIO;

[snip]

my $gffIO =
    new Bio::FeatureIO( -format => 'GFF',
			-version => 3
			);

## Now read in the fasta file to create the landmarks

my $seqIO =
  new Bio::SeqIO( -file => $finished_bacs_repeatmasked,
		  -format => 'fasta'
		);

## Begwin

while(my $seq = $seqIO->next_seq){
  ## Create a new generic feature to hold the properties of the
  ## sequence;

  my $name = $seq->id;
  my $leng = $seq->length;

  my $accn;

  next unless $accn = $finished{$name};

  my $f =
    new Bio::SeqFeature::Generic( -seq_id  => $name,
				  -primary => "SGN",
				  -source  => "genomic_clone",
				  -start   =>     1,
				  -end     => $leng,
				  -tag     => { ID   => $name,
					        Name => $name,
						Dbxref => "GB:$accn",
						Ontology_term => "SO:0000040"
					      }
				);

  print $gffIO->write_feature($f);
  exit;
}
</source>


The script gives me the following error:

Use of uninitialized value in pattern match (m//) at
~/perl5/lib/perl5/Bio/FeatureIO/gff.pm line 107, <FB> line 688.

------------- EXCEPTION: Bio::Root::Exception -------------
MSG: only Bio::SeqFeature::Annotated objects are writeable
STACK: Error::throw
STACK: Bio::Root::Root::throw ~/perl5/lib/perl5/Bio/Root/Root.pm:368
STACK: Bio::FeatureIO::gff::write_feature
~/perl5/lib/perl5/Bio/FeatureIO/gff.pm:272
STACK: ./fasta2gff_landmarks.plx:72
-----------------------------------------------------------


With all due respect to:

* http://www.bioperl.org/wiki/GFF_Refactor and
* http://www.bioperl.org/wiki/GFF_code_audit


I'm a bit stuck.

Any quick fix in the meantime? I realise that using a feature object
is overkill here, but I wanted to learn the 'BioPerl way' of doing
things. Should I drop that approach until the code audit and
refactorization is over, or is this the kind of use-case that can help
with that process?


Thanks for any pointers,
Dan.


P.S. Thanks also to the friendly people on
irc://ircfreenode.net/#bioperl and
irc://irc.freenode.net/#bioinformatics who helped me get this script
together as it is.