bin/agp_to_fasta.pl

   1 #!/usr/bin/env perl
   2 use strict;
   3 use warnings;
   4
   5 =head1 NAME
   6
   7 agp_to_fasta.pl - assemble an AGP file into fasta
   8
   9 =head1 SYNOPSIS
  10
  11 agp_to_fasta.pl foo.agp source.fasta > result.fasta
  12
  13 =head1 OPTIONS
  14
  15 =head2 -w <num>
  16
  17 Column width of fasta output, with 0 being unlimited (no additional
  18 newlines).  Default 80.
  19
  20 =cut
  21
  22 use File::Temp;
  23 use Getopt::Std;
  24 use Pod::Usage;
  25
  26 use Bio::Index::Fasta;
  27
  28 use CXGN::BioTools::AGP qw/ agp_to_seqs /;
  29
  30 my %opts;
  31 getopts( 'w', \%opts );
  32 $opts{w} = 80 unless defined $opts{w};
  33 pod2usage() unless @ARGV;
  34
  35 @ARGV == 2 or die "must pass exactly 2 file names as arguments\n";
  36
  37 for( @ARGV ) {
  38     -r or die "cannot read file '$_'"
  39 }
  40 my ( $agp_file, $source_seqs ) = @ARGV;
  41
  42 my $index = make_seqs_index( $source_seqs );
  43
  44 my @seqs = agp_to_seqs(
  45     $agp_file,
  46     fetch_default => sub {
  47         my $s = $index->fetch( $_[0] )
  48             or return;
  49         return $s->seq;
  50     }
  51    );
  52
  53 my $o = Bio::SeqIO->new(
  54     -format => 'fasta',
  55     -fh     => \*STDOUT,
  56     ( $opts{w} ? (-width  => 80) : () ),
  57    );
  58
  59 # the while/shift cuts down a bit on disk space usage
  60 while( @seqs ) {
  61     $o->write_seq( shift @seqs );
  62 }
  63
  64 exit;
  65
  66 ################# SUBS ###########3
  67
  68 my $temp;
  69 sub make_seqs_index {
  70     my ( $seqs ) = @_;
  71     $temp = File::Temp->new;
  72     $temp->close;
  73
  74     # index of source sequences, for random access
  75     my $inx = Bio::Index::Fasta->new(
  76         -filename => "$temp",
  77         -write_flag => 1 );
  78     $inx->make_index( $source_seqs );
  79
  80     return $inx;
  81 }
  82