[Bio] / FigWebServices / ma_to_tf_nr.cgi Repository:
ViewVC logotype

Diff of /FigWebServices/ma_to_tf_nr.cgi

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1, Fri Jan 20 16:14:37 2006 UTC revision 1.3, Thu Feb 9 04:30:58 2006 UTC
# Line 78  Line 78 
78    
79  sub find_tfs  sub find_tfs
80  {  {
81      print STDERR "find_tfs called\n";
82    my ($fig,$cgi,$html,$spotid_to_ratio)=@_;    my ($fig,$cgi,$html,$spotid_to_ratio)=@_;
83    $new_html = [];    $new_html = [];
84    $dir = "/home/mkubal/public_html";    $dir = "/home/mkubal/public_html";
# Line 86  Line 87 
87    
88    open(OUT2,">$dir/tfs_to_ratio.txt");    open(OUT2,">$dir/tfs_to_ratio.txt");
89    
90      open(OUT3,">$dir/refseq_in_exp.txt");
91    
92    open(IN,"$dir/spotid_to_refseq.txt");    open(IN,"$dir/spotid_to_refseq.txt");
93    %spotid_to_refseq;    %spotid_to_refseq;
94    while ($_ = <IN>){    while ($_ = <IN>){
95        chomp($_);        chomp($_);
       #print STDERR "spot_to_refseq:$_\n";  
96        @temp = split("\t",$_);        @temp = split("\t",$_);
97        $spotid_to_refseq{$temp[0]} = $temp[1];        if($spotid_to_refseq{$temp[0]}){
98              $ref = $spotid_to_refseq{$temp[0]};
99              push(@$ref,$temp[1]);
100              print STDERR "old spotid:$temp[1]\n";
101          }
102          else{
103              $spotid_to_refseq{$temp[0]} = [$temp[1]];
104              print STDERR "new spotid:$temp[1]\n";
105          }
106    }    }
107    
108    open(IN2,"$dir/refseq_to_transfactor.txt.nonredundant");    open(IN2,"$dir/refseq_to_transfactor.txt.nonredundant");
109    %refseq_to_tf;    %refseq_to_tf;
110    while ($_ = <IN2>){    while ($_ = <IN2>){
111        chomp($_);        chomp($_);
112        #print STDERR "refseq_to_transfactor:$_\n";        print STDERR "refseq_to_transfactor:$_\n";
113        @temp = split("\t",$_);        @temp = split("\t",$_);
114        if($refseq_to_tf{$temp[0]}){        if($refseq_to_tf{$temp[0]}){
115            $ref = $refseq_to_tf{$temp[0]};            $ref = $refseq_to_tf{$temp[0]};
# Line 109  Line 119 
119           $refseq_to_tf{$temp[0]} = [$temp[1]];           $refseq_to_tf{$temp[0]} = [$temp[1]];
120        }        }
121    }    }
122      close(IN2);
123    
124    push(@$new_html,"<HTML><HEAD>    push(@$new_html,"<HTML><HEAD>
125    <TITLE>strep</TITLE>    <TITLE>strep</TITLE>
# Line 122  Line 133 
133    my $row_string = "";    my $row_string = "";
134    foreach my $id (@ids)    foreach my $id (@ids)
135    {    {
136          print STDERR "ids here:$id\n";
137        my $ratio = $spotid_to_ratio->{$id};        my $ratio = $spotid_to_ratio->{$id};
138        my $refseq = $spotid_to_refseq{$id};        my $refseqs = $spotid_to_refseq{$id};
139          foreach $refseq (@$refseqs){
140              print OUT3 "$refseq\n";
141    
142        my $tfs_ref = $refseq_to_tf{$refseq};        my $tfs_ref = $refseq_to_tf{$refseq};
143    
144        foreach $tfs (@$tfs_ref){        foreach $tfs (@$tfs_ref){
# Line 136  Line 151 
151                 print OUT2 "$tfs\t$ratio\n";                 print OUT2 "$tfs\t$ratio\n";
152            }            }
153        }        }
154          }
155    
156    }    }
157    
158    close(OUT2);    close(OUT2);
159      close(OUT3);
160    
161    push(@$new_html,"</TABLE>");    push(@$new_html,"</TABLE>");
162    #push(@$new_html,    #push(@$new_html,
# Line 156  Line 173 
173  {  {
174    my ($fig,$cgi,$html)=@_;    my ($fig,$cgi,$html)=@_;
175    $new_html = [];    $new_html = [];
176      %significant;
177    $dir = "/home/mkubal/public_html";    $dir = "/home/mkubal/public_html";
178    
179    print STDERR "made it here\n";    open(IN2,"$dir/refseq_to_transfactor.txt.nonredundant");
180      %tfs_to_refseq;
181      while ($_ = <IN2>){
182          chomp($_);
183          @temp = split("\t",$_);
184          if($tfs_to_refseq{$temp[1]}){
185              $ref = $tfs_to_refseq{$temp[1]};
186              push(@$ref,$temp[0]);
187          }
188          else{
189             $tfs_to_refseq{$temp[1]} = [$temp[0]];
190          }
191      }
192      close(IN2);
193    
194      open(IN3,"$dir/refseq_in_exp.txt");
195      %refseq_in_exp;
196      while ($_ = <IN3>){
197          chomp($_);
198          $refseq_in_exp{$_} = 1;
199      }
200      close(IN3);
201    
202      #print STDERR "made it here\n";
203    open(IN,"$dir/tfs_to_ratio.txt");    open(IN,"$dir/tfs_to_ratio.txt");
204    %tfs_combinations;    %tfs_combinations;
205    %tfs_counts;    %tfs_counts;
   
206    %tfs_expected;    %tfs_expected;
207      %tfs_total;
208    open(IN3,"$dir/stats.txt.nonredundant");    open(IN3,"$dir/stats.txt.nonredundant");
209    open(SUMMARY,">$dir/nonredundant_stats.summary");    open(SUMMARY,">$dir/nonredundant_stats.summary");
210    while ($_ = <IN3>){    while ($_ = <IN3>){
211       @temp = split("\t",$_);       @temp = split("\t",$_);
212       $expected = $temp[2];       $expected = $temp[2];
213          $total = $temp[1];
214       chomp($expected);       chomp($expected);
215       $tfs_expected{$temp[0]} = $expected;       $tfs_expected{$temp[0]} = $expected;
216          $tfs_total{$temp[0]} = $total;
217    }    }
218    close(IN3);    close(IN3);
219    
# Line 186  Line 229 
229    
230    @negative_exclusives;    @negative_exclusives;
231    @positive_exclusives;    @positive_exclusives;
232      @therest;
233    
234    foreach my $k (keys(%tfs_combinations)){    foreach my $k (keys(%tfs_combinations)){
235       print STDERR "k:$k\n";       #print STDERR "k:$k\n";
236       $sign = "not_set";       $sign = "not_set";
237       $exclusive = 1;       $exclusive = 1;
238       foreach $l (@lines){       foreach $l (@lines){
# Line 218  Line 262 
262          }          }
263          else{push(@negative_exclusives,$k)}          else{push(@negative_exclusives,$k)}
264        }        }
265         else{push(@therest,$k)}
266   }   }
267       push(@$new_html,"<HTML><HEAD>       push(@$new_html,"<HTML><HEAD>
268       <TITLE>strep</TITLE>       <TITLE>strep</TITLE>
# Line 230  Line 275 
275    
276       foreach my $tfs (@positive_exclusives){       foreach my $tfs (@positive_exclusives){
277          my $count = $tfs_counts{$tfs};          my $count = $tfs_counts{$tfs};
278          my $observed = $count * (1/2439);          my $observed = $count * (1/311);
279          my $expected = $tfs_expected{$tfs};          my $expected = $tfs_expected{$tfs};
280          my $ratio = $observed/$expected;          my $ratio = $observed/$expected;
281            if($ratio > 10 ){ $significant{$tfs} = $ratio};
282          my $row_string = "<TR><TD>$tfs</TD><TD>$count</TD><TD>$ratio</TD></TR>";          my $row_string = "<TR><TD>$tfs</TD><TD>$count</TD><TD>$ratio</TD></TR>";
283          push(@$new_html,$row_string);          push(@$new_html,$row_string);
284          print SUMMARY "UP\t$tfs\t$observed\t$expected\t$ratio\n";          print SUMMARY "UP\t$tfs\t$observed\t$expected\t$ratio\n";
# Line 248  Line 294 
294       foreach my $tfs (@negative_exclusives)       foreach my $tfs (@negative_exclusives)
295       {       {
296          my $count = $tfs_counts{$tfs};          my $count = $tfs_counts{$tfs};
297          my $observed = $count * (1/2439);          my $observed = $count * (1/311);
298          my $expected = $tfs_expected{$tfs};          my $expected = $tfs_expected{$tfs};
299          my $ratio = $observed/$expected;          my $ratio = $observed/$expected;
300            if($ratio > 10 ){ $significant{$tfs} = $ratio};
301          my $row_string = "<TR><TD>$tfs</TD><TD>$count</TD><TD>$ratio</TD></TR>";          my $row_string = "<TR><TD>$tfs</TD><TD>$count</TD><TD>$ratio</TD></TR>";
302          push(@$new_html,$row_string);          push(@$new_html,$row_string);
303          print SUMMARY "DOWN\t$tfs\t$observed\t$expected\t$ratio\n";          print SUMMARY "DOWN\t$tfs\t$observed\t$expected\t$ratio\n";
304        }
305    
306      push(@$new_html,"</TABLE>");
307    
308      push(@$new_html,"<br><br>");
309    
310      push(@$new_html,"<TABLE border><TR><TH>MIXED Transcription Factor Combinations</TH><TH>Number of Genes</TH><TH>Observed/Expected</TH></TR>");
311    
312      foreach my $tfs (@therest)
313      {
314          my $count = $tfs_counts{$tfs};
315          my $observed = $count * (1/311);
316          my $expected = $tfs_expected{$tfs};
317          my $ratio = $observed/$expected;
318          if($ratio > 10 ){ $significant{$tfs} = $ratio};
319          my $row_string = "<TR><TD>$tfs</TD><TD>$count</TD><TD>$ratio</TD></TR>";
320          push(@$new_html,$row_string);
321          print SUMMARY "MIXED\t$tfs\t$observed\t$expected\t$ratio\n";
322       }       }
323    
324       push(@$new_html,"</TABLE>");       push(@$new_html,"</TABLE>");
325    
326       &HTML::show_page($cgi,$new_html);    push(@$new_html,"<br><br>");
327       exit;  
328      push(@$new_html,"<TABLE border><TR><TH>SIGNIFICANT Transcription Factor Combinations</TH><TH>RefSeq in Exp</TH><TH>Other RefSeq</TH><TH>Observed/Expected</TH><TH>ORPV</TH></TR>");
329    
330    
331      @sorted = sort {$b <=> $a} values(%significant);
332      %sorted_hash;
333      foreach $s (values(%significant)){
334          $sorted_hash{$s} = 1;
335  }  }
336    
337      @final;
338      %already_in;
339      foreach $s (@sorted){
340          foreach $tfs (keys(%significant)){
341              if($already_in{$tfs}){ $do_nothing = 1}
342              else{
343                  $ratio = $significant{$tfs};
344                  if( $ratio >= $s){
345                      push(@final,$tfs);
346                      $already_in{$tfs} = 1;
347                  }
348              }
349          }
350      }
351    
352      foreach my $tfs (@final)
353      {
354          my $ids = $tfs_to_refseq{$tfs};
355          @in_id_strings = ();
356          @other_id_strings = ();
357          foreach $id (@$ids){
358              if($refseq_in_exp{$id}){
359                  $id_string ="<a href='http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?db=NucCore&cmd=search&term=$id'>$id</a>";
360                  push(@in_id_strings,$id_string);
361             }
362             else {
363                 $id_string ="<a href='http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?db=NucCore&cmd=search&term=$id'>$id</a>";
364                 push(@other_id_strings,$id_string);
365             }
366          }
367    
368          $in_ids_line = join(",",@in_id_strings);
369          $other_ids_line = join(",",@other_id_strings);
370          my $ratio = $significant{$tfs};
371          my $count = $tfs_counts{$tfs};
372    
373          my $d = 311;
374          my $m = $tfs_total{$tfs};
375          my $n = 20659;
376          my $orpv = &orpv($count,$d,$m,$n);
377    
378          my $row_string = "<TR><TD>$tfs</TD><TD>$in_ids_line</TD><TD>$other_ids_line</TD><TD>$ratio</TD><TD>$orpv</TD></TR>";
379          push(@$new_html,$row_string);
380          print SUMMARY "SIGNIFICANT\t$tfs\t$id\t$ratio\n";
381      }
382    
383      push(@$new_html,"</TABLE>");
384    
385      &HTML::show_page($cgi,$new_html);
386      exit;
387    }
388    
389  sub parse_inputs  sub parse_inputs
390  {  {
# Line 288  Line 409 
409    
410  }  }
411    
412    sub gammaln{
413        my($x) = @_;
414        $x = $x -1;
415        if($x){
416            my $result = ( $x *log($x) ) - $x;
417            print STDERR "gammaln result: $result\n";
418            return $result;
419    
420        }
421        else{ return 0}
422    }
423    
424    sub lchoose{
425        my($n,$x) = @_;
426        return &gammaln($n+1) - &gammaln($x+1) - &gammaln($n-$x+1);
427    }
428    
429    sub ldhypg {
430        my($x,$d,$m,$n)= @_;
431        if ($x>$m || $d-$x > $n-$m || $d>$n || $x>$d) {
432            print STDERR "x,d,m,n:$x,$d,$m,$n\n";
433            return(-inf)
434        }
435        else
436        {
437            my ($p, $q, $r) = (&lchoose($m,$x), &lchoose($n-$m,$d-$x), &lchoose($n,$d));
438            print STDERR "p, q, r = $p\t$q\t$r\n";
439            return($p + $q - $r);
440            #return(&lchoose($m,$x) + &lchoose($n-$m,$d-$x) - &lchoose($n,$d));
441        }
442    }
443    
444    sub orpv{
445        print STDERR "orpv called\n";
446        my($x,$d,$m,$n) = @_;
447        $pval = 0;
448        for ( $y=$x; $y <= $n; $y++) {
449            $result = &ldhypg($y,$d,$m,$n);
450            $exp_of_result = exp($result);
451            if($exp_of_result == 1) {$do_nothing = 1}
452            else{$pval= $pval+ $exp_of_result}
453        }
454        print STDERR "pval: $pval\n";
455        return($pval);
456    }

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3