gene_matrix3

#!/usr/bin/perl

use strict;

my $in=shift;
my $locus=shift;
my $vcds=shift;
my $vdist=shift;
my $vdhs=shift;
my $thresh=shift;

open(O, "> $locus.out");

my $e = `zgrep $locus $in > tmp.qtl.3`;

my %ma; my %p; my %loc; my $x=2; my %rloc;
open(IN, "< tmp.qtl.3");
while(my $i=<IN>) {
  chomp($i);  
  my @l=split(/\t/,$i);
  if ($l[9]<$thresh) { $p{"$l[2]-$l[5]-$l[6]"}=1; }
}

open(IN, "< tmp.qtl.3");
while(my $i=<IN>) {
  chomp($i);
  my @l=split(/\t/,$i);
  if ($l[2] eq $locus) {
  if (exists $p{"$l[2]-$l[5]-$l[6]"}) {
  if ($l[9] ne "NA") {
    if (!(exists $loc{"$l[5]-$l[6]"})) { $loc{"$l[5]-$l[6]"}=$x; $rloc{$x}="$l[5]-$l[6]"; $x++; }
    $ma{$l[0]}[$loc{"$l[5]-$l[6]"}]=$l[9];
    $ma{$l[0]}[0]=$l[3];
    $ma{$l[0]}[1]=$l[1];
  }}}
}

open(C, "< $vcds");
my %cds; my %xc;
while(my $i=<C>) {
  chomp($i);
  my @l=split(/\t/,$i);
  $l[7]=~/LOCUS=([^;]*)/;
  if ($1 eq $locus) {
  if (exists $ma{$l[2]}) {
    $cds{$l[2]}=$l[12];
    if (!(exists $loc{"$l[12]-CDS"})) { $loc{"$l[12]-CDS"}=$x; $rloc{$x}="$l[12]-CDS"; $xc{$x}=1; $x++; }
    $ma{$l[2]}[$loc{"$l[12]-CDS"}]=1;
  }
  }
}

open(P, "< $vdist");
my %prox; my %xp;
while(my $i=<P>) {
  chomp($i);
  my @l=split(/\t/,$i);
  $l[7]=~/LOCUS=([^;]*)/;
  if ($1 eq $locus) {
  if (exists $ma{$l[2]}) {
    $prox{$l[2]}=$l[11];
    if (!(exists $loc{"$l[11]-DIST"})) { $loc{"$l[11]-DIST"}=$x; $rloc{$x}="$l[11]-DIST"; $xp{$x}=1; $x++; }
    $ma{$l[2]}[$loc{"$l[11]-DIST"}]=1; # $l[12]
  }
  }
}

open(D, "< $vdhs");
my %dhs; my %xd;
while(my $i=<D>) {
  chomp($i);
  my @l=split(/\t/,$i);
  $l[7]=~/LOCUS=([^;]*)/;
  if ($1 eq $locus) {
  if (exists $ma{$l[2]}) {
    $dhs{$l[2]}=$l[12];
    if (!(exists $loc{"$l[12]-DHS"})) { $loc{"$l[12]-DHS"}=$x; $rloc{$x}="$l[12]-DHS"; $xd{$x}=1; $x++; }
    $ma{$l[2]}[$loc{"$l[12]-DHS"}]=1;
  }
  }
}

my $c=0; my $rst;
print O "P\tMAF"; $rst="m\$P ~ m\$MAF"; while($c<$x) { if (exists $rloc{$c}) { print O "\t$rloc{$c}"; $rloc{$c}=~s/-/\./g; $rst.=" + m\$".$rloc{$c}; } $c++; } print O "\n"; 

foreach my $gene (keys %ma) {
  my $y=0; while($y<$x) { 
    if (exists $xc{$y}) { if (!($ma{$gene}[$y])) { $ma{$gene}[$y]=0; } }
    elsif (exists $xd{$y}) { if (!($ma{$gene}[$y])) { $ma{$gene}[$y]=0; } }
    elsif (exists $xp{$y}) { if (!($ma{$gene}[$y])) { $ma{$gene}[$y]=0; } }
    elsif ($y==1) { if (!($ma{$gene}[$y])) { $ma{$gene}[$y]=0; } }
    else { 
      if (!($ma{$gene}[$y])) { $ma{$gene}[$y]=0; } 
      else { $ma{$gene}[$y]=-log($ma{$gene}[$y])/log(10); }
    }
    $y++; 
  }

  print O join("\t",@{$ma{$gene}})."\n";
}

#open(R, "> tmp.R");
#print R "m <- read.table(\"$locus.out\",header=T)\n";
#print R "ml <- lm($rst)\n";
#print R "summary(ml)\n";
#my $go = `Rscript tmp.R`;

my $go = `Rscript lasso.R $locus.out`;

print $go;