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

View of /FigMetagenomeTools/coverage.pl

Parent Directory Parent Directory | Revision Log Revision Log

Revision - (download) (as text) (annotate) (vendor branch)
Mon Feb 19 17:15:26 2007 UTC (13 years, 1 month ago) by olson
Branch: x, 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, y, HEAD
Changes since 1.1: +0 -0 lines
Initial import

#!/usr/bin/perl -w

# count the coverage from a sequencher consensus output file
# basically start with a sequence line, and then in the next lines count if
# the base is the same
# reset if it is a blank line

use strict;

my $usage=<<EOF;
-i input file
-s only count same bases
-w window size (default = 100)


my ($file, $same, $size);
while (@ARGV)
    my $t=shift;
    if ($t eq "-i") {$file=shift}
    if ($t eq "-s") {$same=1}
    if ($t eq "-w") {$size=shift}

unless ($size) {$size=100}
die $usage unless ($file);

my $seq;
my ($yposn, $xposn,$rowbeginning, $length)=(0,0,0,0);
my $master; my $coverage; 
open(IN, $file) || die "can't open $file";
while (<IN>)
    next if (/sequencher/i || /Summary View/i);
    next if (/^\s+\*/ || /^\s+\+/);
    next if (/^\s+\#/);
    next if (/^\s+\.+$/);
    if (/^\s*$/) 
    my $seq=$1;
    unless ($seq) 
	    print STDERR "Couldn't parse |$_|\n";
    #$seq =~ s/\s+//g;
    if (!$yposn)
	    for (my $i=0; $i<=length($seq); $i++)
		    my $base=substr($seq,$i,1); 
	    for (my $i=0; $i<=length($seq); $i++)
		    my $base=substr($seq,$i,1);
		    ($same && ($master->[$xposn] eq uc($base))) ? $coverage->[$xposn]++ : 
			    ($base =~ /[ACGT]/i) ? $coverage->[$xposn]++ : 1;

my $line; my $cov;
my $ma; my $window;
for (my $i=0; $i<=($rowbeginning+$xposn); $i++)
    next unless ($master->[$i]);
    next if ($master->[$i] eq ":");
    next if ($master->[$i] eq " ");
    unless ($coverage->[$i]) {$coverage->[$i]=0}
    if ($window && ($window == $size))
     print "$i\t", $ma/$window, "\n";

    #print "$i\t", $coverage->[$i], "\n";
    ($coverage->[$i]) ? eval{$cov.=$coverage->[$i]} : eval{$cov.=0};
    if ($i && !($i % 60))
    #	print "$line\n$cov\n";
#	($line, $cov)=('','');

#print "$line\n$cov\n";

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3