use strict;

use warnings;

use Carp;

my $opts = parse_params();

do_stats($opts);

exit;

sub error

{

my (@msg) = @_;

if ( scalar @msg )

{

croak join('',@msg);

}

die

"About: Generate some .fa stats. Currently GC content only.\n",

"Usage: ref-stats [OPTIONS]\n",

"Options:\n",

" -r, --refseq <file> The reference sequence.\n",

" -h, -?, --help This help message.\n",

"\n";

}

sub parse_params

{

my $opts = {};

while (my $arg=shift(@ARGV))

{

if ( $arg eq '-r' || $arg eq '--refseq' ) { $$opts{'refseq'} = shift(@ARGV); next }

if ( $arg eq '-?' || $arg eq '-h' || $arg eq '--help' ) { error(); }

error("Unknown parameter \"$arg\". Run -h for help.\n");

}

if ( !exists($$opts{'refseq'}) ) { error("Missing the -r option.\n") }

return $opts;

}

sub do_stats

{

my ($self) = @_;

my $_len = 60;

my %gc_counts = ();

open(my $fh,'<',$$opts{refseq}) or error("$$opts{refseq}: $!");

while (my $line=<$fh>)

{

if ( $line=~/^>/ ) { next; }

chomp($line);

my $len = length($line);

if ( $len ne $_len ) { next; }

my $gc_count = 0;

for (my $i=0; $i<$len; $i++)

{

my $base = substr($line,$i,1);

if ( $base eq 'g' || $base eq 'G' || $base eq 'c' || $base eq 'C' ) { $gc_count++; }

}

$gc_counts{$gc_count}++;

}

close($fh) or error("close $$opts{refseq}: $!");

my %gc_freqs = ();

while (my ($len,$count) = each %gc_counts)

{

my $gc = $len * 100./$_len;

$gc_freqs{$gc} = $count;

}

delete($gc_freqs{0});

use Data::Dumper;

print "# Generated by ref-stats\n";

print Data::Dumper->Dump([\%gc_freqs],[qw(gc_freqs)]);

}