#!/usr/bin/perl

use strict;

my $seq = "tair9.fa";
my $db;
my $chr;
open(SEQ,$seq) || die "Could not open $seq\n";
while(<SEQ>) {
  chomp;
  if ( />(\S+)/ ) {
    $chr =  $1;
    print "chr=$chr\n";
  }	
  else {
    $db->{$chr} .= $_;
  }
}

open(FA, "opa.fa") || die "Could not open opa.fa\n";
my $fa;
my $id;
my $is_snp;
while(<FA>) {
  chomp;
  if ( />(.+)$/ ) {
    $id = $1;
  }
  else {
    $fa->{$id} = $_;

    if ( $id =~ /^(.+)-left/ ) {
      $is_snp->{$1} =1;
    } elsif ( $id =~ /^(.+)-right/ ) {
      $is_snp->{$1}=1;
    }
  }
}


open(BLAT, "tair9.psl") || die "Could not open tair9.psl\n";
for(my$n=0;$n<5;$n++) {
  $_ = <BLAT>;
}
my $best;
my $root;

while(<BLAT>) {
  chomp;
  my ($match, $mismatch, $repmatch, $N, $Qgapcount, $Qgapbases, $Tgapcount, $Tgapbases, $strand, $Qname, $Qsize, $Qstart, $Qend, $Tname, $Tsize, $Tstart, $Tend, $blockCount, $blockSizes, $qstarts,$tStarts ) =  split;
  if ( $best->{$Qname}[0] < $match ) {
    $best->{$Qname} = [$match, $mismatch, $repmatch, $N, $Qgapcount, $Qgapbases, $Tgapcount, $Tgapbases, $strand, $Qname, $Qsize, $Qstart, $Qend, $Tname, $Tsize, $Tstart, $Tend, $blockCount, $blockSizes, $qstarts,$tStarts];
    if ( $Qname =~ /^(.+)-left$/ ) {
      $root->{$Qname} = $1;
    }
    elsif ( $Qname =~ /^(.+)-right$/ ) {
      $root->{$Qname} = $1;
    }
  }
}
my $n=0;

my $snpfile = "/Net/spud/raid/GSCAN/ARABIDOPSIS/DATA/founders.genotypes.txt";
open(SNPS,$snpfile) || die "Could not open $snpfile\n";
$_ = <SNPS>;
chomp;
my @snps = split(/\s+/);
shift @snps;
my $genos;

while(<SNPS> ) {
  chomp;
  my ( $founder, @alleles ) = split;
  for(my $k=0;$k<=$#snps;$k++) {
    my $snp = @snps[$k];
    $genos->{$founder}{$snp} = $alleles[$k];
  }
}
open(MAP, ">map.strand.txt");
my @map;

foreach my $Q ( sort keys %$is_snp ) {
  $n++;
  my $left = "$Q-left";
  my $right = "$Q-right";
  my $pos;
  my $strand;
  if ( exists($best->{$left}) && exists($best->{$right} ) ) {
    my $d;
    if ( $best->{$left}[8] eq "+" ) {
      my $off1 = $best->{$left}[10]-$best->{$left}[12];
      my $off2 = $best->{$right}[11];
      $d = $best->{$right}[15] - $best->{$left}[16]-$off1+$off2;
      $pos = $best->{$left}[16]+2+$off1;
      $strand = +1;
    }
    else {
      my $off1 = $best->{$left}[10]-$best->{$left}[12];
      my $off2 = $best->{$right}[11];
      $d = $best->{$left}[15] - $best->{$right}[16]+$off1-$off2;
      $pos = $best->{$right}[16]+2+$off1;
      $strand = -1;
    }
    
    if ( $best->{$left}[13] eq $best->{$right}[13] && abs($d) < 10 ) {
      my $txt = $d==1 ? "" : "*";
      print "$n $Q $pos $best->{$left}[13] $best->{$left}[8] $best->{$left}[15] $best->{$left}[16]  $best->{$right}[15] $best->{$right}[16] $d $txt\n";
    }
    else {
      print "$n XXX $Q $best->{$left}[13] $best->{$left}[8] $best->{$left}[15] $best->{$left}[16]  $best->{$right}[15] $best->{$right}[16] $d\n";
      if ( exists( ($best->{$left}))) {
	my $off1 = $best->{$left}[10]-$best->{$left}[12];
	$pos = $best->{$left}[16]+2+$off1;
      }
      elsif (exists( ($best->{$right}))) {
	my $off1 = $best->{$right}[10]-$best->{$right}[12];
	$pos = $best->{$right}[15]-1-$off1;
      }
    }
  } 
  else {
    print "badS $Q left $left $fa->{$left} $best->{$left}[13] $best->{$left}[16] right $right $fa->{$right} $best->{$right}[13] $best->{$right}[15]\n";
  }
  push @map, [ $best->{$left}[13], $pos, $Q, $strand] if ( defined($pos));
  
}

@map = sort { my $x = $a->[0] cmp $b->[0]; $x = $a->[1] <=> $b->[1] if ( ! $x ) ; return($x)} @map;

foreach my $ref ( @map ) {
  print join( "\t", @$ref ) . "\n";
  my $pos = $ref->[1];
  my $chr = $ref->[0];
  my $seq = substr($db->{$chr},$pos-1,4);
  my $snp = $ref->[2];
  my $col = $genos->{Col}{$snp};
  print MAP join( "\t", $snp, $chr, $pos, $seq, $col ) . "\n";
}

foreach my $c ( 1,2,3,4,5 ) {
  my $chr = "Chr$c";
  my $out = "chr$c.MAGIC.map";
  open(OUT,">$out" ) || die "Could not write $out\n";
  print OUT "marker\tchromosome\tbuild\tbp\n";
  foreach my $ref ( @map ) {
    if ( $chr eq $ref->[0] ) {
      my $pos = $ref->[1];
      my $snp = $ref->[2];
      print OUT join( "\t", $snp, $c, "tair9", $pos ) . "\n";
    }
  }
  close(OUT);
}
  
