[Bio] / FigMetagenomeTools / quality_scores.pl Repository:
ViewVC logotype

View of /FigMetagenomeTools/quality_scores.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (download) (as text) (annotate)
Fri Mar 16 20:51:25 2007 UTC (12 years, 11 months ago) by olson
Branch: MAIN
CVS Tags: mgrast_dev_08112011, mgrast_dev_08022011, mgrast_dev_05262011, mgrast_dev_04082011, mgrast_version_3_2, mgrast_dev_12152011, mgrast_dev_06072011, mgrast_rel_2008_0806, mgrast_dev_10262011, mgrast_dev_02212011, mgrast_rel_2008_0923, mgrast_release_3_0, mgrast_dev_03252011, mgrast_rel_2008_0924, mgrast_rel_2008_1110_v2, mgrast_rel_2008_0625, mgrast_release_3_0_4, mgrast_release_3_0_2, mgrast_release_3_0_3, mgrast_release_3_0_1, mgrast_dev_03312011, mgrast_release_3_1_2, mgrast_release_3_1_1, mgrast_release_3_1_0, mgrast_dev_04132011, mgrast_dev_04012011, mgrast_rel_2008_0919, mgrast_rel_2008_1110, myrast_33, mgrast_rel_2008_0917, mgrast_dev_04052011, mgrast_dev_02222011, HEAD
Changes since 1.1: +0 -1 lines
initial tweaks

#!/usr/bin/perl -w

# count the quality scores for the files from Pedulla

use strict;
use Rob;
my $rob=new Rob;

my $usage=<<EOF;
$0 <options>

-f file
-c cutoff (default=20)

EOF

my ($file, $cutoff)=('', 20);
while (@ARGV) {
 my $test=shift;
 if ($test eq "-f") {$file=shift @ARGV}
 elsif ($test eq "-c") {$cutoff=shift @ARGV}
}

die $usage unless $file;

my $fasta=$rob->read_fasta($file);

print join "\t", "Sequence", "Number of bases", "Mean", "Median", "Standard Deviation", "Longest run of sequences greater than $cutoff", "\n"; 
foreach my $key (keys %$fasta) {
 my $max=0;
 my $current=0;
 my @values=split /\s+/, $fasta->{$key};
 #print STDERR join "|", $key, @values, "\n";
 my @temp;
 # this removes undefs from the array
 foreach (@values) {if (/\d/) {push @temp, $_}}
 @values=@temp;
 foreach my $val (@values) {
  if ($val > 99) {
   print STDERR "Problem with $val\n";
  }
  if ($val < $cutoff) {
   if ($current > $max) {$max=$current}
   $current=0;
  } 
  else {
   $current++;
  }
 }

 my $label=$key;
 $label =~ s/\s*PHD.*//;
 print join "\t", $label, scalar(@values), $rob->mean(\@values), $rob->median(\@values), $rob->stdev(\@values), $max, "\n";
}



MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3