#! /usr/bin/gawk -f
# Last edited on 2011-12-09 01:10:24 by stolfi

BEGIN{
  abort = -1;
  split("", ps_scaff);   # Scaffold of probeset, indexed by probeset name.
  split("", ps_start);   # Start of probeset, indexed by probeset name.
  split("", ps_end);     # End of probeset, indexed by probeset name.
  np = read_probeset_table("probeset.txt",ps_scaff,ps_start,ps_end);
  split("", ps_khnames); # KH genes overlapping clones, indexed by probeset name.
}

function read_probeset_table(fname,ps_scaff,ps_start,ps_end,    ntbl,nlin,lin,fld,nfld,ps_name,tmp)
{
  ntbl=0; # Number of entries in table.
  nlin=0; # Number of lines read.
  while((getline lin < fname) > 0) { 
    nlin++;
    if (! match(lin, /^[ \011]*([\#]|$)/))
      { gsub(/[\011]/, " ", lin);
        gsub(/^[ ]+/, "", lin);
        gsub(/[ ]+$/, "", lin);
        gsub(/[ ][ ]+/, " ", lin);
        nfld = split(lin, fld, /[ ]+/);
        if (nfld != 13)
          { tbl_error(fname, nlin, lin, "bad number of fields in table"); }
        if (fld[2] !~ /^Scaffold_[0-9]+$/)
          { tbl_error(fname, nlin, lin, "bad scaffold name"); }
        if ((fld[3] !~ /^[0-9]+$/) || (fld[4] !~ /^[0-9]+$/))
          { tbl_error(fname, nlin, lin, ("bad position range = \"" fld[3] " " fld[4] "\"")); }
        ps_name = fld[5];
        if (ps_name in ps_start) 
          { tbl_error(fname, nlin, lin, ("repeated probeset name = \"" ps_name "\"")); }
        ps_scaff[ps_name] = fld[2];
        ps_start[ps_name] = fld[3];
        ps_end[ps_name] = fld[4];
        # printf "%s %s %d %d\n", ps_name, ps_scaff[ps_name], ps_start[ps_name], ps_end[ps_name] > "/dev/stderr";
        ntbl++;
      }
  }
  if (ERRNO != "0") { tbl_error(fname, nlin, ERRNO); }
  close (fname);
  if (nlin == 0) { arg_error(("file \"" fname "\" empty or missing")); }
  printf "loaded %6d map pairs from %s\n", ntbl, fname > "/dev/stderr"
  return ntbl;
}

(abort >= 0) { exit abort; }

# Discard comment lines:
/^[ \011]*([\#]|$)/ { next; }

# Cleanup:
// { 
  gsub(/[\015]/, "");
  gsub(/[\011]/, " ");
  gsub(/^[ ]+/, "");
  gsub(/[ ]+$/, "");
  gsub(/[ ][ ]+/, " ");
}

# Parse a KH gene model line:
($3 ~ /^Scaffold_[0-9]+$/) {
  if (NF != 11) { data_error(("wrong number of fields")); }
  kh_name = $2;
  kh_scaff = $3;
  nex = $9;
  kh_starts = $10;
  gsub(/[,]$/, "", kh_starts);
  nstarts = split(kh_starts, kh_start, ",");
  kh_ends = $11;
  gsub(/[,]$/, "", kh_ends);
  nends = split(kh_ends, kh_end, ",");
  if (nstarts != nex) { data_error(("wrong number of exon starts")); }
  if (nends != nex)   { data_error(("wrong number of exon ends")); }
  for (i = 1; i <= nex; i++) {
    if (kh_start[i] !~ /^[0-9]+$/) { data_error(("bad exon start = \"" kh_start[i] "\"")); }
    if (kh_end[i]   !~ /^[0-9]+$/) { data_error(("bad exon end = \"" kh_start[i] "\"")); }
  }
        
  kh_psnames = "";
  # printf "[kh] %-30s %3d %s\n", kh_name, nex, kh_scaff > "/dev/stderr";
  for (ps_name in ps_scaff) {
    # printf "  [ps] %-20s %d %d\n", ps_scaff[ps_name], ps_start[ps_name], ps_end[ps_name] > "/dev/stderr";
    if (kh_scaff == ps_scaff[ps_name]) {
      for (i = 1; i <= nex; i++) {
        # printf "    [ex] %d %d\n", kh_start[i], kh_end[i] > "/dev/stderr";
        over_a = (kh_start[i]+0 <= ps_end[ps_name]+0);
        over_b = (kh_end[i]+0 >= ps_start[ps_name]+0);
        if (over_a && over_b) {
          # Overlap detected:
          ps_khnames[ps_name] = (ps_khnames[ps_name] kh_name ",");
          kh_psnames = (kh_psnames ps_name ",");
        }
      }
    }
  }
  kh_psnames = clean_blanks(kh_psnames);
  printf "%-30s %-20s %s\n", kh_name, kh_scaff, (kh_psnames != "" ? kh_psnames : "(NONE)") > "ps-KH-inv.txt";
  next;
}

// { data_error(("bad line format")); }

END{
  close("ps-KH-inv.txt");  printf "ps-KH-inv.txt closed\n" > "/dev/stderr";
  if (abort >= 0) { exit(abort); }
  printf "writing ps-KH-dir.txt:\n" > "/dev/stderr";
  for (ps_name in ps_scaff) {
    khns = clean_blanks(ps_khnames[ps_name]);
    printf "%-30s %-20s %12d %12d %s\n", \
      ps_name, ps_scaff[ps_name], ps_start[ps_name], ps_end[ps_name], \
      (khns != "" ? khns : "(NONE)") > "ps-KH-dir.txt";
  }
  close("ps-KH-dir.txt"); printf "ps-KH-dir.txt closed\n" > "/dev/stderr";
}

function clean_blanks(x)
{
  gsub(/^[ ]+/, "", x);
  gsub(/[ ]+$/, "", x);
  gsub(/[ ][ ]+/, " ", x);
  gsub(/^[(]NONE[)]$/, "", x);
  return x;
}

function data_warning(msg) {
  printf "%s:%d: !! %s\n", FILENAME, FNR, msg > "/dev/stderr";
  printf "  %s\n", $0  > "/dev/stderr";
}

function data_error(msg) {
  printf "%s:%d: ** %s\n", FILENAME, FNR, msg > "/dev/stderr";
  printf "  %s\n", $0  > "/dev/stderr";
  abort = 1;
  exit abort;
}

function tbl_error(f,n,lin,msg)
{ 
  printf "%s:%d: %s\n", f, n, msg > "/dev/stderr";
  printf "  %s\n", lin  > "/dev/stderr";
  abort = 1;
  exit 1
}