perl/etc/justonce/sperm/doplotkmercmp.pl

   1 #!/bin/env perl
   2 use strict;
   3 use warnings;
   4 use IO::Unread qw(unread);
   5 use Data::Dump qw(ddx);
   6
   7 #die "Usage: $0 <input1> <input2> <output>\n" if @ARGV < 3;
   8 #my ($inf1,$inf2,$outf)=@ARGV;
   9
  10 my (%Files,%Fname2ID);
  11 open L,'<','kmerfreq.lst' or die $!;
  12 while (<L>) {
  13         chomp;
  14         my @t = split /\//,$_;
  15         my @k = split /\./,$t[2];
  16         if ($t[2] =~ /^Donor\./) {
  17                 push @{$Files{$k[1]}{'Donor'}},$_;
  18         } else {
  19                 push @{$Files{$k[1]}{$t[1]}},$_;
  20         }
  21         $Fname2ID{$_} = $k[0];
  22 }
  23 ddx \%Files;
  24 close L;
  25
  26 my $prefix = 'perl plotkmercnt.pl 1e-8 ';
  27 my $outpath = 'out';
  28 my $linending = " &\n";
  29
  30 sub getout($$$$) {
  31         my ($x,$y,$k,$outpath) = @_;
  32         return join(' ',$x,$y,"$outpath/$k-$Fname2ID{$x}-$Fname2ID{$y}");
  33 }
  34
  35 open OUT,'>','run.sh' or die $!;
  36 for my $k ( keys %Files ) {
  37         my @types = keys %{$Files{$k}};
  38         for my $i ( 0 .. ($#types-1) ) {
  39                 for my $j ( ($i+1) .. $#types ) {
  40                         print "$k: $i,$j\n";
  41                         my $arrA = $Files{$k}->{$types[$i]};
  42                         my $arrB = $Files{$k}->{$types[$j]};
  43                         # pairwise
  44                         for my $x ( 0 .. $#$arrA ) {
  45                                 for my $y ( 0 .. $#$arrB ) {
  46                                         print "p $arrA->[$x], $arrB->[$y]\n";
  47                                         print OUT $prefix,getout($arrA->[$x],$arrB->[$y],$k,$outpath),$linending;
  48                                 }
  49                         }
  50                 }
  51         }
  52         for my $t ( @types ) {
  53                 print "$k: $t\n";
  54                 my $arrA = $Files{$k}->{$t};
  55                 # within A
  56                 for my $x ( 0 .. ($#$arrA-1) ) {
  57                         for my $y ( ($x+1) .. $#$arrA ) {
  58                                 print "w $arrA->[$x], $arrA->[$y]\n";
  59                                 print OUT $prefix,getout($arrA->[$x],$arrA->[$y],$k,$outpath),$linending;
  60                         }
  61                 }
  62         }
  63 }
  64 close OUT;
  65
  66 __END__
  67 Well, Combine(7,2) = 21 . I should use this directly.
  68 But we may use different prefix for different tyoes (p or w) later, right ?
  69
  70 For the duplicated ones:
  71 ps -ef |grep perl|perl -lane 'BEGIN {my (%x,%y);} @a=split /\s+/;$b=join(",",@a[10,11,12]);++$x{$b};push @{$y{$b}},$a[1]; END {my @z; for (keys %x) { if ($x{$_}>1) {print join("-",$x{$_},$_,"\t",@{$y{$_}}); push @z,$y{$_}->[0] };  } print join(" ",@z);print $#z }'