7
use Bio::Tools::SeqStats;
14
'f|format:s' => \$format,
17
'a|aggregate!'=> \$aggregate,
21
my $USAGE = "usage: gccalc.pl -f format -i filename\n";
26
$file = shift unless $file;
29
print "Could not open file [$file]\n$USAGE" and exit unless -e $file;
30
$seqin = new Bio::SeqIO(-format => $format,
33
$seqin = new Bio::SeqIO(-format => $format,
36
my ($total_base, $total_gc);
37
while( my $seq = $seqin->next_seq ) {
38
next if( $seq->length == 0 );
39
if( $seq->alphabet eq 'protein' ) {
40
warn("gccalc does not work on amino acid sequences ...skipping this seq");
44
my $seq_stats = Bio::Tools::SeqStats->new('-seq'=>$seq);
45
my $hash_ref = $seq_stats->count_monomers(); # for DNA sequence
46
print "Seq: ", $seq->display_id, " ";
47
print $seq->desc if $seq->desc;
48
print " Len:", $seq->length, "\n";
49
$total_base += $seq->length;
50
$total_gc += $hash_ref->{'G'} + $hash_ref->{'C'};
51
printf "GC content is %.4f\n", ($hash_ref->{'G'} + $hash_ref->{'C'}) /
54
foreach my $base (sort keys %{$hash_ref}) {
55
print "Number of bases of type ", $base, "= ", $hash_ref->{$base},"\n";
60
printf "Total GC content is %.4f out of %d bases\n", $total_gc / $total_base, $total_base;
62
# alternatively one could use code submitted by
67
my @seqarray = split('',$seq);
69
foreach my $base (@seqarray) {
70
$count++ if $base =~ /[G|C]/i;
73
my $len = $#seqarray+1;
82
bp_gccalc - GC content of nucleotide sequences
86
bp_gccalc [-f/--format FORMAT] [-h/--help] filename
88
bp_gccalc [-f/--format FORMAT] < filename
90
bp_gccalc [-f/--format FORMAT] -i filename
94
This scripts prints out the GC content for every nucleotide sequence
99
The default sequence format is fasta.
101
The sequence input can be provided using any of the three methods:
105
=item unnamed argument
111
bp_gccalc -i filename
123
User feedback is an integral part of the evolution of this and other
124
Bioperl modules. Send your comments and suggestions preferably to
125
the Bioperl mailing list. Your participation is much appreciated.
127
bioperl-l@bioperl.org - General discussion
128
http://bioperl.org/wiki/Mailing_lists - About the mailing lists
130
=head2 Reporting Bugs
132
Report bugs to the Bioperl bug tracking system to help us keep track
133
of the bugs and their resolution. Bug reports can be submitted via the
136
https://redmine.open-bio.org/projects/bioperl/
138
=head1 AUTHOR - Jason Stajich
140
Email jason@bioperl.org
144
Based on script code (see bottom) submitted by cckim@stanford.edu
146
Submitted as part of bioperl script project 2001/08/06