[Bio] / StandaloneTools / find_overlaps_standalone Repository:
ViewVC logotype

View of /StandaloneTools/find_overlaps_standalone

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (download) (annotate) (vendor branch)
Mon Jun 21 21:46:19 2004 UTC (15 years, 4 months ago) by olson
Branch: MAIN, gordon
CVS Tags: gordon1, HEAD
Changes since 1.1: +0 -0 lines
Initial import of gordon's utilities.

#!/usr/bin/perl -w

(($dir = shift) && (-d $dir)) || die "$dir does not exist";

@pegs = ();
if (!-s "$dir/Features/peg/tbl")
{
    die "zero-size $dir/Features/peg/tbl";
}
else
{
    foreach $entry (`cut -f1,2 $dir/Features/peg/tbl`)
    {
	chomp $entry;
	(undef, $loc) = split /\t/, $entry;
	if ($loc =~ /^([^,]+)_(\d+)_(\d+)$/)
	{
	    push(@pegs,[$1,$2,$3]);
	}
    }
    $num_pegs = (scalar @pegs);
    print STDERR "Read $num_pegs pegs from $dir\n";
    
    $overlaps = 0;
    @pegs = sort { ($a->[0] cmp $b->[0]) or (&min($a->[1],$a->[2]) <=> &min($b->[1],$b->[2])) } @pegs;
    for ($i=0; ($i < @pegs); $i++)
    {
	for ($j=$i+1; (($j < @pegs) && ($overlap = &overlap($pegs[$i],$pegs[$j]))); ++$j)
	{
	    if ($overlap > 40) { ++$overlaps; }
	}
    }
    print "$dir\t", int(0.5 + 100*$overlaps/(&max($num_pegs,1))), "%\t$num_pegs\t$overlaps\n";
}

sub overlap {
    my($p1,$p2) = @_;
    
    my ($left1, $right1, $left2, $right2);
    
    if ($p1->[0] ne $p2->[0]) { return 0 }
    
    $left1  = &min($p1->[1],$p1->[2]);
    $left2  = &min($p2->[1],$p2->[2]);
    $right1 = &max($p1->[1],$p1->[2]);
    $right2 = &max($p2->[1],$p2->[2]);
    
    if ($left1 > $left2)
    { 
	($left1, $left2)   = ($left2, $left1); 
	($right1, $right2) = ($right2, $right1);
    }
    
    $ov = 0;
    if ($right1 >= $left2) { $ov = &min($right1,$right2) - $left2 + 1; }
    
    return $ov;
}

sub min {
    my ($x, $y) = @_;
    
    return (($x < $y) ? $x : $y);
}

sub max {
    my ($x, $y) = @_;
    
    return (($x > $y) ? $x : $y);
}

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3