[Bio] / FigKernelScripts / FFB2_usable_motifs.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/FFB2_usable_motifs.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 my $usage = "usage: usable_motifs Dir";
2 :     use IO::File;
3 :     use strict;
4 : olson 1.2 use List::Util qw(reduce sum);
5 : overbeek 1.1 use Data::Dumper;
6 :    
7 :     my $dir;
8 :     (
9 :     ($dir = shift @ARGV)
10 :     )
11 :     || die $usage;
12 :    
13 : olson 1.3 my $range = "7-12";
14 :     if (@ARGV)
15 :     {
16 :     $range = shift;
17 :     }
18 :    
19 :     my($min_k, $max_k);
20 :    
21 :     if ($range =~ /^\d+$/)
22 :     {
23 :     $min_k = $max_k = $range;
24 :     }
25 :     elsif ($range =~ /^(\d+)-(\d+)/)
26 :     {
27 :     $min_k = $1;
28 :     $max_k = $2;
29 :     }
30 :     else
31 :     {
32 :     die "Invalid range '$range'\n";
33 :     }
34 :    
35 : overbeek 1.1 if (! -d $dir) { mkdir($dir,0777) || die "could not make $dir" }
36 :     my $i;
37 :    
38 :     my @files;
39 : olson 1.3 for ($i=$min_k; ($i <= $max_k); $i++)
40 : overbeek 1.1 {
41 :     if (! -d "$dir/$i" ) { mkdir("$dir/$i",0777) || die "could not make $dir/$i" }
42 : olson 1.2 $files[$i] = IO::File->new("| gzip -c > $dir/$i/good.oligos.gz");
43 :     # $files[$i] = IO::File->new(">$dir/$i/good.oligos");
44 : overbeek 1.1 }
45 :    
46 :     my $line = <STDIN>;
47 :     while ($line && ($line =~ /^(\S+)\t(\S+)/))
48 :     {
49 : olson 1.3 my $curr = substr($1,0,$min_k);
50 : overbeek 1.1 my @pairs;
51 : olson 1.3 while ($line && ($line =~ /^(\S+)\t(\S+)/) && (substr($1,0,$min_k) eq $curr))
52 : overbeek 1.1 {
53 : olson 1.2 push(@pairs,[$1,$2,length($1)]); # [Oligo,FuncIndex]
54 : overbeek 1.1 $line = <STDIN>;
55 :     }
56 :     &process_set(\@pairs,\@files);
57 :     }
58 :    
59 : olson 1.4 for ($i=$min_k; ($i <= $max_k); $i++)
60 :     {
61 :     close($files[$i]);
62 :     }
63 :    
64 : overbeek 1.1 sub process_set {
65 :     my($pairs,$files) = @_;
66 :    
67 :     my $i;
68 : olson 1.3 for ($i=$min_k; ($i <= $max_k); $i++)
69 : overbeek 1.1 {
70 : olson 1.2 my $fh = $files->[$i];
71 : overbeek 1.1 my %counts;
72 : olson 1.2 # foreach my $pair (@$pairs)
73 :     # {
74 :     # my($oligo,$funcI,$len) = @$pair;
75 :    
76 :     # my $short_oligo;
77 :     # if ($len >= $i)
78 :     # {
79 :     # $short_oligo = substr($oligo,0,$i);
80 :     # $counts{$short_oligo}->{$funcI}++;
81 :     # }
82 :     # }
83 :    
84 :     #
85 :     # The above optimizes to the following.
86 :     #
87 :    
88 :     for (@$pairs)
89 : overbeek 1.1 {
90 : olson 1.2 $counts{substr($_->[0], 0, $i)}->{$_->[1]}++ if $_->[2] >= $i;
91 :     }
92 : overbeek 1.1
93 :    
94 :     foreach my $short_oligo (sort keys(%counts))
95 :     {
96 :     my $funcI_hash = $counts{$short_oligo};
97 : olson 1.2
98 :     # my @tmp = sort { $funcI_hash->{$b} <=> $funcI_hash->{$a}} keys(%$funcI_hash);
99 :     # my $sum = 0;
100 :     # foreach my $funcI (@tmp) { $sum += $funcI_hash->{$funcI} }
101 :     #
102 :     # The above optimizes to the following. No need to sort when we
103 :     # just want the max value.
104 :     #
105 :    
106 :     my($k, $v, $max, $maxv);
107 : overbeek 1.1 my $sum = 0;
108 : olson 1.2 while (($k, $v) = each %$funcI_hash)
109 :     {
110 :     if ($v > $maxv)
111 :     {
112 :     $maxv = $v;
113 :     $max = $k;
114 :     }
115 :     $sum += $v;
116 :     }
117 :    
118 :     # my $sum = sum values %$funcI_hash;
119 :     # my $max = reduce { $funcI_hash->{$a} > $funcI_hash->{$b} ? $a : $b } keys %$funcI_hash;
120 :    
121 :     if ($maxv / $sum > 0.9)
122 :     # if (($funcI_hash->{$tmp[0]} / $sum) > 0.9)
123 : overbeek 1.1 {
124 : olson 1.2 print $fh $short_oligo, "\t", $max, "\n";
125 :     # print $fh "$short_oligo\t$tmp[0]\t" . $funcI_hash->{$tmp[0]}/$sum . "\t" . $funcI_hash->{$tmp[0]} . "\n";
126 :     # print $fh join("\t",($short_oligo,
127 :     # $tmp[0],
128 :     # sprintf("%0.3f",$funcI_hash->{$tmp[0]}/$sum),
129 :     # $funcI_hash->{$tmp[0]})),"\n";
130 : overbeek 1.1 }
131 :     }
132 :     }
133 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3