perl/etc/crep/gen_cluster_input.pl

   1 #!/usr/bin/perl -w
   2 use strict;
   3 use warnings;
   4 use Data::Dump;
   5
   6 open IA,'<','crep_all_tsv_new.txt.up2.txt' or die $!;
   7 open IB,'<','crep_all_tsv_new.txt.up2.rev.txt' or die $!;
   8 open IR,'<','resLOC.rpratio' or die $!;
   9
  10 open OA,'>','ricexpro.up2.clu.txt' or die $!;
  11 open OB,'>','ricexpro.up2.rev.clu.txt' or die $!;
  12 open OS,'>','ricexpro.up2.stat.txt' or die $!;
  13
  14 my (%dA,%dB,%dFlag,%cntdflag);
  15 while (<IA>) {
  16         my ($id,$cnt,$value) = split /\t/;
  17         $id =~ s/^LOC_Os/t/;
  18         $dA{$id} = log($value) / log(2);
  19         ++$dFlag{$id};
  20 }
  21 close IA;
  22 while (<IB>) {
  23         my ($id,$cnt,$value) = split /\t/;
  24         $id =~ s/^LOC_Os/t/;
  25         $dB{$id} = log($value) / log(2);
  26         $dFlag{$id} += 100;
  27 }
  28 close IB;
  29
  30 for my $k (keys %dFlag) {
  31         push @{ $cntdflag{$dFlag{$k}} },$k;
  32 }
  33 for my $k (keys %cntdflag) {
  34         print "$k -> ",scalar @{ $cntdflag{$k} },': ';
  35         print "\n";     # join(',',@{ $cntdflag{$k} }),
  36 }
  37 # 1 -> 5332:
  38 # 101 -> 531:
  39 # 100 -> 8173:
  40 my ($max,$min,$sum,$n,$ss)=(-999,999,0,0,0);
  41 for my $id (@{ $cntdflag{'101'} }) {
  42         #print OS "# $id: $dA{$id}, $dB{$id}\n";
  43         $sum += $dA{$id}+$dB{$id};
  44         $ss += $dA{$id}*$dA{$id} + $dB{$id}*$dB{$id};
  45         $n += 2;
  46         $max = $dA{$id} if $max < $dA{$id};
  47         $min = $dA{$id} if $min > $dA{$id};
  48         $max = $dB{$id} if $max < $dB{$id};
  49         $min = $dB{$id} if $min > $dB{$id};
  50 }
  51 my $mean = $sum/$n;
  52 print "Range:[$min,$max] Mean: $mean SD:",sqrt($ss/$n - $mean*$mean)," Cnt:$n\n";
  53 =pod
  54 1 -> 5332:
  55 101 -> 531:
  56 100 -> 8173:
  57 Range:[-5.28337053983209,5.82964609302319] Mean: -0.575308843011672 SD:1.62278420169048 Cnt:1062
  58 =cut
  59 my $outHead = <IR>;
  60 print OS "$outHead";
  61 print OA "$outHead";
  62 print OB "$outHead";
  63 while (<IR>) {
  64         my ($id,$acc,@dat) = split /\t/;
  65         my $LOCid = $id;
  66         $id =~ s/^LOC_Os/t/;
  67         my $tmp = join('', map {my $a; if (defined $_) {$a = int($_); $a="+$a" if $_>=0; $a="-$a" if $_<0 and $a==0; $a; } } ($dA{$id},$dB{$id}));
  68         my $newLOCid = $LOCid . $tmp;
  69         my $newid = "${id}_${acc}$tmp";
  70         if ($dFlag{$id} == '101') {
  71                 print OS join("\t",$newid,   $dA{$id} .';'. $dB{$id}   ,@dat);
  72         } elsif ($dFlag{$id} == '1') {
  73                 print OA join("\t",$newid,$newLOCid,@dat);
  74         } elsif ($dFlag{$id} == '100') {
  75                 print OB join("\t",$newid,$newLOCid,@dat);
  76         } else { die; }
  77 }
  78 close IR;
  79 close OS;
  80 close OA;
  81 close OB;
  82
  83 __END__
  84 $ wc -l crep_all_tsv_new.txt.up2.*txt ricexpro.up2.* resLOC.rpratio
  85    8704 crep_all_tsv_new.txt.up2.rev.txt
  86    5863 crep_all_tsv_new.txt.up2.txt
  87    4756 ricexpro.up2.clu.txt
  88    7344 ricexpro.up2.rev.clu.txt
  89     401 ricexpro.up2.stat.txt
  90   12499 resLOC.rpratio
  91   39567 total
  92
  93 ./cluster -f ricexpro.up2.clu.txt -l -cg a -ng -g 2 -e 2
  94 ./cluster -f ricexpro.up2.rev.clu.txt -l -cg a -ng -g 2 -e 2