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

View of /FigMetagenomeTools/coverage.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (download) (as text) (annotate) (vendor branch)
Mon Feb 19 17:15:26 2007 UTC (12 years, 9 months 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;
$0
-i input file
-s only count same bases
-w window size (default = 100)

EOF

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>)
{
    chomp;
    next if (/sequencher/i || /Summary View/i);
    next if (/^\s+\*/ || /^\s+\+/);
    next if (/^\s+\#/);
    next if (/^\s+\.+$/);
    if (/^\s*$/) 
    {
        $rowbeginning+=$length;
        $yposn=0;
	$length=0;
        next;
    } 
    #m/^\s*\S.*\#\S+\s+([\S\:].*)\s*$/; 
    m/(.{43})$/;
    my $seq=$1;
    unless ($seq) 
    {
	    print STDERR "Couldn't parse |$_|\n";
	    next;
    }
    $xposn=$rowbeginning;
    #$seq =~ s/\s+//g;
    if (!$yposn)
    {
	    $length=length($seq)+1;
	    for (my $i=0; $i<=length($seq); $i++)
	    {
		    my $base=substr($seq,$i,1); 
		    $master->[$xposn]=uc($base);
		    $xposn++;
	    }
	    $yposn++;
    }
    else 
    {
	    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;
		    $xposn++;
	    }
	    $yposn++;
    }
}

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";
     $window=0;
     $ma=0;
    }
    $window++;
    $ma+=$coverage->[$i];

    
    #print "$i\t", $coverage->[$i], "\n";
    $line.=$master->[$i];
    ($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