6
use Bio::Tools::SeqStats;
13
'f|format:s' => \$format,
16
'a|aggregate!'=> \$aggregate,
20
my $USAGE = "usage: gccalc.pl -f format -i filename\n";
25
$file = shift unless $file;
28
print "Could not open file [$file]\n$USAGE" and exit unless -e $file;
29
$seqin = new Bio::SeqIO(-format => $format,
32
$seqin = new Bio::SeqIO(-format => $format,
35
my ($total_base, $total_gc);
36
while( my $seq = $seqin->next_seq ) {
37
next if( $seq->length == 0 );
38
if( $seq->alphabet eq 'protein' ) {
39
warn("gccalc does not work on amino acid sequences ...skipping this seq");
43
my $seq_stats = Bio::Tools::SeqStats->new('-seq'=>$seq);
44
my $hash_ref = $seq_stats->count_monomers(); # for DNA sequence
45
print "Seq: ", $seq->display_id, " ";
46
print $seq->desc if $seq->desc;
47
print " Len:", $seq->length, "\n";
48
$total_base += $seq->length;
49
$total_gc += $hash_ref->{'G'} + $hash_ref->{'C'};
50
printf "GC content is %.4f\n", ($hash_ref->{'G'} + $hash_ref->{'C'}) /
53
foreach my $base (sort keys %{$hash_ref}) {
54
print "Number of bases of type ", $base, "= ", $hash_ref->{$base},"\n";
59
printf "Total GC content is %.4f out of %d bases\n", $total_gc / $total_base, $total_base;
61
# alternatively one could use code submitted by
66
my @seqarray = split('',$seq);
68
foreach my $base (@seqarray) {
69
$count++ if $base =~ /[G|C]/i;
72
my $len = $#seqarray+1;
81
gccalc - GC content of nucleotide sequences
85
gccalc [-f/--format FORMAT] [-h/--help] filename
87
gccalc [-f/--format FORMAT] < filename
89
gccalc [-f/--format FORMAT] -i filename
93
This scripts prints out the GC content for every nucleotide sequence
98
The default sequence format is fasta.
100
The sequence input can be provided using any of the three methods:
104
=item unnamed argument
122
User feedback is an integral part of the evolution of this and other
123
Bioperl modules. Send your comments and suggestions preferably to
124
the Bioperl mailing list. Your participation is much appreciated.
126
bioperl-l@bioperl.org - General discussion
127
http://bioperl.org/wiki/Mailing_lists - About the mailing lists
129
=head2 Reporting Bugs
131
Report bugs to the Bioperl bug tracking system to help us keep track
132
of the bugs and their resolution. Bug reports can be submitted via the
135
https://redmine.open-bio.org/projects/bioperl/
137
=head1 AUTHOR - Jason Stajich
139
Email jason@bioperl.org
143
Based on script code (see bottom) submitted by cckim@stanford.edu
145
Submitted as part of bioperl script project 2001/08/06