[Bio] / StandaloneTools / find_overlaps_standalone.pl Repository:
ViewVC logotype

Annotation of /StandaloneTools/find_overlaps_standalone.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (view) (download) (as text)

1 : parrello 1.1 #!/usr/bin/perl -w
2 :    
3 :     (($dir = shift) && (-d $dir)) || die "$dir does not exist";
4 :    
5 :     @pegs = ();
6 :     if (!-s "$dir/Features/peg/tbl")
7 :     {
8 :     die "zero-size $dir/Features/peg/tbl";
9 :     }
10 :     else
11 :     {
12 :     foreach $entry (`cut -f1,2 $dir/Features/peg/tbl`)
13 :     {
14 :     chomp $entry;
15 :     (undef, $loc) = split /\t/, $entry;
16 :     if ($loc =~ /^([^,]+)_(\d+)_(\d+)$/)
17 :     {
18 :     push(@pegs,[$1,$2,$3]);
19 :     }
20 :     }
21 :     $num_pegs = (scalar @pegs);
22 :     print STDERR "Read $num_pegs pegs from $dir\n";
23 :    
24 :     $overlaps = 0;
25 :     @pegs = sort { ($a->[0] cmp $b->[0]) or (&min($a->[1],$a->[2]) <=> &min($b->[1],$b->[2])) } @pegs;
26 :     for ($i=0; ($i < @pegs); $i++)
27 :     {
28 :     for ($j=$i+1; (($j < @pegs) && ($overlap = &overlap($pegs[$i],$pegs[$j]))); ++$j)
29 :     {
30 :     if ($overlap > 40) { ++$overlaps; }
31 :     }
32 :     }
33 :     print "$dir\t", int(0.5 + 100*$overlaps/(&max($num_pegs,1))), "%\t$num_pegs\t$overlaps\n";
34 :     }
35 :    
36 :     sub overlap {
37 :     my($p1,$p2) = @_;
38 :    
39 :     my ($left1, $right1, $left2, $right2);
40 :    
41 :     if ($p1->[0] ne $p2->[0]) { return 0 }
42 :    
43 :     $left1 = &min($p1->[1],$p1->[2]);
44 :     $left2 = &min($p2->[1],$p2->[2]);
45 :     $right1 = &max($p1->[1],$p1->[2]);
46 :     $right2 = &max($p2->[1],$p2->[2]);
47 :    
48 :     if ($left1 > $left2)
49 :     {
50 :     ($left1, $left2) = ($left2, $left1);
51 :     ($right1, $right2) = ($right2, $right1);
52 :     }
53 :    
54 :     $ov = 0;
55 :     if ($right1 >= $left2) { $ov = &min($right1,$right2) - $left2 + 1; }
56 :    
57 :     return $ov;
58 :     }
59 :    
60 :     sub min {
61 :     my ($x, $y) = @_;
62 :    
63 :     return (($x < $y) ? $x : $y);
64 :     }
65 :    
66 :     sub max {
67 :     my ($x, $y) = @_;
68 :    
69 :     return (($x > $y) ? $x : $y);
70 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3