Mercurial > repos > plus91-technologies-pvt-ltd > softsearch

comparison 2.4/src/cluster.pair.pl @ 13:e3609c8714fb draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author plus91-technologies-pvt-ltd
date Fri, 30 May 2014 03:37:55 -0400
parents
children
comparison
equal deleted inserted replaced
12:980fce54e60a 13:e3609c8714fb
1 #!/usr/bin/perl
2 use strict;
3 use POSIX;
4
5 my $usage = "cluster.pair.pl maxdist\n";
6 my $maxdist = shift or die $usage;
7
8 my %count;
9
10 while (<STDIN>){
11 chomp;
12 my ($sample, $chrstart, $start, $chrend, $end) = split /\t/;
13 my $nstart = floor ($start/$maxdist);
14 my $nend = floor ($end/$maxdist);
15 my $coord = {start=>$start, end=>$end};
16
17 push @{$count{$chrstart}->{$nstart}->{$chrend}->{$nend}->{$sample}}, $coord;
18 }
19
20 print_groups (\%count);
21
22 sub print_groups {
23 my ($rcount) = @_;
24 my %count = %{$rcount};
25
26 foreach my $chrstart (sort {$a<=>$b} keys %count) {
27 foreach my $posstart (sort {$a<=>$b} keys %{$count{$chrstart}}) {
28 my %fcoord = %{$count{$chrstart}->{$posstart}};
29
30 foreach my $chrend (sort {$a<=>$b} keys %fcoord) {
31 foreach my $posend (sort {$a<=>$b} keys %{$fcoord{$chrend}}){
32 my @nsamples = sort {$a cmp $b} (keys %{$fcoord{$chrend}->{$posend}});
33
34 my $cpos = $fcoord{$chrend}->{$posend};
35
36 my @coords;
37 my $totnum=0;
38
39 foreach my $sample (@nsamples) {
40 my ($num, $avgx, $avgy) = calc_moments(@{$cpos->{$sample}});
41 push (@coords, {start=>$avgx, end=>$avgy});
42 $totnum+=$num;
43 }
44
45 my ($num, $avgx, $avgy) = calc_moments(@coords);
46
47 print $chrstart."\t".$avgx."\t".$chrend."\t".$avgy ."\t".$num."\t".$totnum."\t" ;
48
49 print $_."\t" foreach (@nsamples);
50 print "\n";
51 }
52 }
53 }
54 }
55 }
56
57 sub calc_moments {
58 my (@pos) = @_;
59
60 my ($num, $sumx, $sumy) = (0,0,0);
61 foreach my $cpos (@pos) {
62 $num++;
63 $sumx+=$cpos->{start};
64 $sumy+=$cpos->{end};
65 }
66 my $avgx = sprintf ("%d", $sumx/$num);
67 my $avgy = sprintf ("%d", $sumy/$num);
68
69 return ($num, $avgx, $avgy);
70 }