/*
  bed2bed_display.c

  Convert BED file to BED format suitable for displaying ChIP-seq peaks.
  
  # Arguments:
  #
  #   BED File
  #   Minimal BED score - 5th field (integer)  [-a min]
  #   Maximal BED score - 5th field (integer)  [-b max]
  #
  # Options:
  #
  #   Set Feature name    [-f fname]
  #   BED output formats: [-o 1|2]  def=1
  #
  #     1- BED track with score and strand fields -use score to visualize the track 
  #        the score value will determine the level of gray in which this feature is displayed 
  #        (higher numbers = darker gray).
  #
  #        Track Parameters:
  #
  #        ex: track name=STAT1Peaks description="Robertson 2007, STAT1 peaks" visibility=2 useScore=1
  #
  #     2- BedGraph format
  #        The bedGraph format allows display of continuous-valued data in track format.
  #        chromA  chromStartA  chromEndA  dataValueA
  #        chromB  chromStartB  chromEndB  dataValueB 
  #       
  #        Parameters for bedGraph track definition lines:
  #
  #        track type=bedGraph name=track_label description=center_label visibility=1|2|3 color=r,g,b graphType=bar|points

  Giovanna Ambrosini, EPFL/ISREC, Giovanna.Ambrosini@epfl.ch

  Copyright (c) 2014 EPFL and Swiss Institute of Bioinformatics.

  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program.  If not, see <http://www.gnu.org/licenses/>.

*/
/*
#define DEBUG 
*/
#define _GNU_SOURCE
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <errno.h>
#include <unistd.h>
#include <ctype.h>
#include <getopt.h>
#include <assert.h>
#include <sys/stat.h>
#include "hashtable.h"
#ifdef DEBUG
#include <mcheck.h>
#endif

#include "version.h"

#define BUF_SIZE 8192
#define LINE_SIZE 1024
#define FT_MAX  64
#define CHR_NB 18
#define AC_MAX 18
#define POS_MAX 16
#define SCORE_MAX 12
#define NAME_MAX 128
#define FIELD_MAX 64

typedef struct _options_t {
  int help;
  int debug;
  char *feature;
  int oformat;
  char *autoscale;
  char *always0;
  char *wfunction;
  char *smoothing;
  char *visibility;
  char *trackName;
  char *trackDesc;
  char *trackColor;
} options_t;

static options_t options;

int minScore = -1;
int maxScore = -1;

int rescale_score(int s)  /* Rescale score range to [167,1000] */
{
  int a = 167;
  int b = 1000;
  int val;

  if (maxScore-minScore > 0) {
    val = ((b-a)*(s-minScore)/(maxScore-minScore)) + a;   
  } else {
    val = s;
  }
  return(val);
}

int
process_bed(FILE *input, char *iFile) 
{
  unsigned long start, end;
  int score;
  int new_score;
  int count = 1;
  char *s, *res, *buf;
  size_t bLen = LINE_SIZE;

  if (options.debug)
    fprintf(stderr, " Processing BED file %s\n", iFile);
  if ((s = malloc(bLen * sizeof(char))) == NULL) {
    perror("process_bed: malloc");
    exit(1);
  }
  /* Print Track Header */
  if (options.oformat == 1) {
    printf("track name=\"%s\" description=\"%s\" visibility=%s useScore=1\n", options.trackName, options.trackDesc, options.visibility);
  } else if (options.oformat == 2) {
    printf("track type=bedGraph name=\"%s\" description=\"%s\" visibility=%s color=%s autoScale=%s alwaysZero=%s maxHeightPixels=100:50:20 graphType=bar priority=30 windowingFunction=%s smoothingWindow=%s\n", options.trackName, options.trackDesc, options.visibility, options.trackColor, options.autoscale, options.always0, options.wfunction, options.smoothing);
  }
#ifdef DEBUG
  int c = 1; 
#endif
  while ((res = fgets(s, (int) bLen, input)) != NULL) {
    size_t cLen = strlen(s);
    char bed_fld[12][FIELD_MAX]; 
    char strand = '\0';
    char field[32];
    int i = 0;

    while (cLen + 1 == bLen && s[cLen - 1] != '\n') {
      bLen *= 2;
      if ((s = realloc(s, bLen)) == NULL) {
        perror("process_file: realloc");
        exit(1);
      }
      res = fgets(s + cLen, (int) (bLen - cLen), input);
      cLen = strlen(s);
    }
    if (s[cLen - 1] == '\n')
      s[cLen - 1] = 0;
    else if (s[cLen - 1] == '\n' && s[cLen - 2] == '\r')
      s[cLen - 2] = 0;

    buf = s;
    /* printf ("BED line : %s\n", s); */
    /* Check line */
    /* Get first word/field */
    if ( sscanf(s, "%31[^ ]", field) == 1 ) {
      if ( (strcmp(field, "track") == 0) || (strcmp(field, "browser") == 0) || (strcmp(field, "itemRgb") == 0) )
        continue;
    }
    /* Get BED fields */
    /* Chrom NB */
    while (*buf != 0 && !isspace(*buf)) {
      if (i >= CHR_NB) {
        bed_fld[0][i] = 0;
        fprintf(stderr, "Chrom NB too long \"%s\" Skipping line...\n", bed_fld[0]);
        break;
      }
      bed_fld[0][i++] = *buf++;
    }
    if (i < AC_MAX)
      bed_fld[0][i] = 0;
    else
      continue;
    while (isspace(*buf))
      buf++;
    /* Start Position */
    i = 0;
    while (isdigit(*buf)) {
      if (i >= POS_MAX) {
        fprintf(stderr, "Start position too large \"%s\" \n", buf);
        exit(1);
      }
      bed_fld[1][i++] = *buf++;
    }
    bed_fld[1][i] = 0;
    start = (unsigned long)atoi(bed_fld[1]);
    while (isspace(*buf))
      buf++;
    /* End Position */
    i = 0;
    while (isdigit(*buf)) { 
      if (i >= POS_MAX) {
        fprintf(stderr, "End position too large \"%s\" \n", buf);
        exit(1);
      }
      bed_fld[2][i++] = *buf++;
    }
    bed_fld[2][i] = 0;
    end = (unsigned long)atoi(bed_fld[2]);
    while (isspace(*buf))
      buf++;
    /* Name/Feature */
    i = 0;
    while (*buf != 0 && !isspace(*buf)) {
      if (i >= FIELD_MAX) {
        fprintf(stderr, "Field Name too long \"%s\" \n", buf);
        fclose(input);
        exit(1);
      }
      bed_fld[3][i++] = *buf++;
    }
    bed_fld[3][i] = 0;
    while (isspace(*buf))
      buf++;
    /* Score */
    i = 0;
    while (isdigit(*buf) || *buf == '+' || *buf == '-' || *buf == '.') { 
      if (i >= SCORE_MAX) {
        fprintf(stderr, "Score too large \"%s\" \n", buf);
        exit(1);
      }
      bed_fld[4][i++] = *buf++;
    }
    bed_fld[4][i] = 0;
    score = atoi(bed_fld[4]);
    while (isspace(*buf))
      buf++;
    /* Field 6/Strand */
    i = 0;
    while (*buf != 0 && !isspace(*buf)) {
      if (i >= FIELD_MAX) {
        fprintf(stderr, "Field 6 too long \"%s\" \n", buf);
        fclose(input);
        exit(1);
      }
      bed_fld[5][i++] = *buf++;
    }
    bed_fld[5][i] = 0;
    while (isspace(*buf))
      buf++;
   /* Field 7 */
    i = 0;
    while (*buf != 0 && !isspace(*buf)) {
      if (i >= FIELD_MAX) {
        fprintf(stderr, "Field 7 too long \"%s\" \n", buf);
        fclose(input);
        exit(1);
      }
      bed_fld[6][i++] = *buf++;
    }
    bed_fld[6][i] = 0;
    while (isspace(*buf))
      buf++;
   /* Field 8 */
    i = 0;
    while (*buf != 0 && !isspace(*buf)) {
      if (i >= FIELD_MAX) {
        fprintf(stderr, "Field 8 too long \"%s\" \n", buf);
        fclose(input);
        exit(1);
      }
      bed_fld[7][i++] = *buf++;
    }
    bed_fld[7][i] = 0;
    while (isspace(*buf))
      buf++;
   /* Field 9 */
    i = 0;
    while (*buf != 0 && !isspace(*buf)) {
      if (i >= FIELD_MAX) {
        fprintf(stderr, "Field 9 too long \"%s\" \n", buf);
        fclose(input);
        exit(1);
      }
      bed_fld[8][i++] = *buf++;
    }
    bed_fld[8][i] = 0;
    while (isspace(*buf))
      buf++;
   /* Field 10 */
    i = 0;
    while (*buf != 0 && !isspace(*buf)) {
      if (i >= FIELD_MAX) {
        fprintf(stderr, "Field 10 too long \"%s\" \n", buf);
        fclose(input);
        exit(1);
      }
      bed_fld[9][i++] = *buf++;
    }
    bed_fld[9][i] = 0;
    while (isspace(*buf))
      buf++;
   /* Field 11 */
    i = 0;
    while (*buf != 0 && !isspace(*buf)) {
      if (i >= FIELD_MAX) {
        fprintf(stderr, "Field 11 too long \"%s\" \n", buf);
        fclose(input);
        exit(1);
      }
      bed_fld[10][i++] = *buf++;
    }
    bed_fld[10][i] = 0;
    while (isspace(*buf))
      buf++;
   /* Field 12 */
    i = 0;
    while (*buf != 0 && !isspace(*buf)) {
      if (i >= FIELD_MAX) {
        fprintf(stderr, "Field 12 too long \"%s\" \n", buf);
        fclose(input);
        exit(1);
      }
      bed_fld[11][i++] = *buf++;
    }
    bed_fld[11][i] = 0;
    while (isspace(*buf))
      buf++;
#ifdef DEBUG
    printf(" [%d]   Chr nb: %s  Start: %lu  End: %lu  Name: %s  Score: %d  Strand: %c  Field #7 %s  Field #8 %s  Field #9 %s  Field #10 %s  Field #11 %s  Field #12 %s\n", c++, bed_fld[0], start, end, bed_fld[3], score, bed_fld[5], bed_fld[6], bed_fld[7], bed_fld[8], bed_fld[9], bed_fld[10], bed_fld[11]);
#endif
    /* Set Strand field  */
    if (bed_fld[5][0] == '+' || bed_fld[5][0] == '-')
      strand = bed_fld[5][0];
    else
      strand = '.';
    /* Set Feature/Name field */
    if (options.feature == NULL) { 
      if (bed_fld[3][0] != '\0') {
        options.feature = malloc((strlen(bed_fld[3])+1) * sizeof(char)); 
        strcpy(options.feature, bed_fld[3]);
      }
    }
    /* Print out BED line */
    if (options.oformat == 1) {
      /* Rescale Score */
      new_score = rescale_score(score);
      printf("%s\t%lu\t%lu\t%s%d\t%d\t%c\n", bed_fld[0], start, end, options.feature, count, new_score, strand);
      count++;
    } else if (options.oformat == 2) {
      printf("%s\t%lu\t%lu\t%d\n", bed_fld[0], start, end, score);
    }
  } /* End of While */
  if (input != stdin) {
    fclose(input);
  }
  return 0;
}

int
main(int argc, char *argv[])
{
#ifdef DEBUG
  mcheck(NULL);
  mtrace();
#endif
  FILE *input;
  int i = 0;
  int usage = 0;
  options.oformat = 1;

  static struct option long_options[] =
    {
      /* These options may or may not set a flag.
         We distinguish them by their indices. */
      {"debug",      no_argument,       0, 'd'},
      {"help",       no_argument,       0, 'h'},
      {"feature",    required_argument, 0, 'f'},
      {"oformat",    required_argument, 0, 'o'},
      {"minscore",   required_argument, 0, 'a'},
      {"maxscore",   required_argument, 0, 'b'},
      {"name",       required_argument, 0,  0 },
      {"desc",       required_argument, 0,  0 },
      {"color",      required_argument, 0,  0 },
      {"autoscale",  required_argument, 0,  0 },
      {"always0",    required_argument, 0,  0 },
      {"wfunction",  required_argument, 0,  0 },
      {"smoothing",  required_argument, 0,  0 },
      {"visibility", required_argument, 0,  0 },
      {0, 0, 0, 0}
    };

  int option_index = 0;

  while (1) {
    int c = getopt_long(argc, argv, "dhf:a:b:o:", long_options, &option_index);
    if (c == -1)
      break;
    switch (c) {
      case 'd':
        options.debug = 1;
        break;
      case 'h':
        options.help = 1;
        break;
      case 'f':
        options.feature = optarg;
	break;
      case 'a':
        minScore = atoi(optarg);
        break;
      case 'b':
        maxScore = atoi(optarg);
        break;
      case 'o':
        options.oformat = atoi(optarg);
        break;
      case 0:
        /* This option is to set the annotation track line */
        if (strcmp(long_options[option_index].name, "name") == 0) {
          options.trackName = optarg;
        }
        if (strcmp(long_options[option_index].name, "desc") == 0) {
          options.trackDesc = optarg;
        }
        if (strcmp(long_options[option_index].name, "color") == 0) {
          options.trackColor = optarg;
        }
        if (strcmp(long_options[option_index].name, "autoscale") == 0) {
          options.autoscale = optarg;
        }
        if (strcmp(long_options[option_index].name, "always0") == 0) {
          options.always0 = optarg;
        }
        if (strcmp(long_options[option_index].name, "wfunction") == 0) {
          options.wfunction = optarg;
        }
        if (strcmp(long_options[option_index].name, "smoothing") == 0) {
          options.smoothing = optarg;
        }
        if (strcmp(long_options[option_index].name, "visibility") == 0) {
          options.visibility = optarg;
        }
        break;
      default:
        printf ("?? getopt returned character code 0%o ??\n", c);
        usage = 1;
    }
  }
  if (optind > argc || options.help == 1 || usage || minScore == -1 || maxScore == -1) {
    fprintf(stderr, "Usage: %s [options] [-a <min BED score>] [-b <max BED score>] [<] <BED file|stdin>\n"
             "      - version %s\n"
             "      where options are:\n"
	     "  \t\t -d|--debug                Produce Debug information\n"
	     "  \t\t -h|--help                 Show this Help text\n"
             "  \t\t -f|--feature <ft>         Set Feature name <ft>\n"
             "  \t\t -o|--oformat <1|2>        Set output Format (1:BED[Def] 2:BedGraph)\n"
             "  \t\t --name <name>             Set name for track name field [def. name=SGA-feature]\n"
             "  \t\t --desc <desc>             Set track description field [def. desc=\"ChIP-Seq Custom data\"]\n"
             "  \t\t --color <col>             Define the track color in comma-separated RGB values [def. 100,100,100]\n"
             "  \t\t --autoscale  <on|off>     Data viewing paramenter: set auto-scale to UCSC data view [def=OFF]\n"
             "  \t\t --always0    <on|off>     Data viewing paramenter: always include zero [def=OFF]\n"
             "  \t\t --wfunction  <func>       Data viewing paramenter: windowing function [def=mean+whiskers|maximum|mean|minimum]\n"
             "  \t\t --smoothing  <grade>      Data viewing paramenter: smoothing window [def=OFF[0], <grade>=0,2..16]\n"
             "  \t\t --visibility <mode>       Display mode: [def=dense|full|hide] or [def=1|2|3]\n"
	     "\n\tConvert BED file to BED format suitable for displaying ChIP-seq peaks.\n"
	     " \tTwo formats are supported: 1) BED Track with Score field; 2) BedGraph.\n"
	     " \tFor BED track format, only the graphical parameter 'visibility' can be set.\n\n",
	     argv[0], VERSION);
      return 1;
  }
  if (argc > optind) {
      if(!strcmp(argv[optind],"-")) {
          input = stdin;
      } else {
          input = fopen(argv[optind], "r");
          if (NULL == input) {
              fprintf(stderr, "Unable to open '%s': %s(%d)\n",
                  argv[optind], strerror(errno), errno);
             exit(EXIT_FAILURE);
          }
          if (options.debug)
             fprintf(stderr, "Processing file %s\n", argv[optind]);
      }
  } else {
      input = stdin;
  }
  if (options.feature != NULL) {
    char *s = options.feature;
    i = 0;
    while (*s != 0  && !isspace(*s)) {
      if (i >= FT_MAX) {
        fprintf(stderr, "Feature Name too long \"%s\" \n", options.feature);
        return 1;
      }
      s++;
    }
  }
  if (options.trackName == NULL) {
    if (options.feature != NULL) {
      if ((options.trackName =  malloc((strlen(options.feature) + 1) * sizeof(char))) == NULL) {
        perror("process_sga: malloc Feature");
        exit(1);
      }
      strcpy(options.trackName, options.feature);
    } else {
      if ((options.trackDesc = malloc(24 * sizeof(char))) == NULL) {
        perror("main: malloc trackDesc");
        exit(1);
      }
      strcpy(options.trackName, "ChIP-SeqPeaks");
    }
  }
  if (options.trackDesc == NULL) {
    if ((options.trackDesc = malloc(32 * sizeof(char))) == NULL) {
      perror("main: malloc trackDesc");
      exit(1);
    }
    strcpy(options.trackDesc, "ChIP-Seq Custom data");
  }
  if (options.trackColor == NULL) {
    if ((options.trackColor = malloc(12 * sizeof(char))) == NULL) {
      perror("main: malloc trackColor");
      exit(1);
    }
    strcpy(options.trackColor, "0,200,100");
  }
  if (options.autoscale == NULL) {
    if ((options.autoscale = malloc(4 * sizeof(char))) == NULL) {
      perror("main: malloc autoscale");
      exit(1);
    }
    strcpy(options.autoscale, "off");
  }
  if (options.always0 == NULL) {
    if ((options.always0 = malloc(4 * sizeof(char))) == NULL) {
      perror("main: malloc always0");
      exit(1);
    }
    strcpy(options.always0, "off");
  }
  if (options.wfunction == NULL) {
    if ((options.wfunction = malloc(5 * sizeof(char))) == NULL) {
      perror("main: malloc wfunction");
      exit(1);
    }
    strcpy(options.wfunction, "mean");
  }
  if (options.smoothing == NULL) {
    if ((options.smoothing = malloc(4 * sizeof(char))) == NULL) {
      perror("main: malloc smoothing");
      exit(1);
    }
    strcpy(options.smoothing, "off");
  }
  if (options.visibility == NULL) {
    if ((options.visibility = malloc(5 * sizeof(char))) == NULL) {
      perror("main: malloc visibility");
      exit(1);
    }
    if (options.oformat == 1) {
      strcpy(options.visibility, "1");
    } else if (options.oformat == 2) {
      strcpy(options.visibility, "dense");
    }
  }
  if (options.debug) {
    fprintf(stderr, " Arguments:\n");
    if (options.feature != NULL)
      fprintf(stderr, " Feature : %s\n", options.feature);
    fprintf(stderr, " Minimal BED Score: %d\n", minScore);
    fprintf(stderr, " maximal BED score: %d\n", maxScore);
    if (options.oformat == 1) {
      fprintf(stderr, " BED Track Display options:\n");
      fprintf(stderr, " Visibility: %s\n", options.visibility);
      fprintf(stderr, " useScore=1\n");
    } else if (options.oformat == 2) {
      fprintf(stderr, " BedGraph Track Display options:\n");
      fprintf(stderr, " -----------------------------\n");
      fprintf(stderr, " Track name: %s\n", options.trackName);
      fprintf(stderr, " Track description: %s\n", options.trackDesc);
      fprintf(stderr, " Track color: %s\n", options.trackColor);
      fprintf(stderr, " Autoscale: %s\n", options.autoscale);
      fprintf(stderr, " Always0: %s\n", options.always0);
      fprintf(stderr, " Wfunction: %s\n", options.wfunction);
      fprintf(stderr, " Smoothing: %s\n", options.smoothing);
      fprintf(stderr, " Visibility: %s\n", options.visibility);
    }
  }
  if (process_bed(input, argv[optind++]) != 0) {
    return 1;
  }
  return 0;
}