static char rcsid[] = "$Id: substring.c 199475 2016-10-23 23:21:59Z twu $";
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "substring.h"
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>		/* For log and exp */

#include "assert.h"
#include "mem.h"
#include "maxent_hr.h"
#include "listdef.h"
#include "list.h"
#include "complement.h"
#include "genome128_hr.h"
#include "mapq.h"
#include "pairdef.h"		/* For State_T */
#include "pair.h"		/* For Pair_print_gsnap */
#include "comp.h"
#include "splice.h"


#define TRIM_MATCH_SCORE 1
#define KNOWN_SPLICESITE_EDGE 1	/* Cannot be 0; otherwise will hit found splice sites */

#define CONSISTENT_TEXT "consistent"
#define TRANSLOCATION_TEXT "translocation"
#define INVERSION_TEXT "inversion"
#define SCRAMBLE_TEXT "scramble"

#define END_SPLICESITE_SEARCH 10
#define MIN_EXON_LENGTH 9
#define END_SPLICESITE_PROB_MATCH 0.90
#define END_SPLICESITE_PROB_MISMATCH 0.95


#ifdef DEBUG
#define debug(x) x
#else
#define debug(x)
#endif

/* mark_mismatches */
#ifdef DEBUG1
#define debug1(x) x
#else
#define debug1(x)
#endif

/* Substring_new */
#ifdef DEBUG2
#define debug2(x) x
#else
#define debug2(x)
#endif

/* Substring_overlap_p and Substring_insert_length */
#ifdef DEBUG3
#define debug3(x) x
#else
#define debug3(x)
#endif


/* splice site probs */
#ifdef DEBUG4
#define debug4(x) x
#else
#define debug4(x)
#endif


/* Substring_convert_to_pairs */
#ifdef DEBUG6
#define debug6(x) x
#else
#define debug6(x)
#endif


/* contains known splicesite */
#ifdef DEBUG7
#define debug7(x) x
#else
#define debug7(x)
#endif


/* trimming.  may also want to turn on DEBUG8 in pair.c */
#ifdef DEBUG8
#define debug8(x) x
#else
#define debug8(x)
#endif

/* bad_stretch_p */
#ifdef DEBUG9
#define debug9(x) x
#else
#define debug9(x)
#endif


/* trim_novel_spliceends */
#ifdef DEBUG13
#define debug13(x) x
#else
#define debug13(x)
#endif


#define LOG_99 -0.01005033585
#define LOG_01 -4.605170186

#if 0
/* Switches on 5 consecutive mismatches */
#define LOG_99_9999 -0.01015034085
#define LOG_99_0001 -9.220390708
#define LOG_25_0001 -10.59663473
#define LOG_25_9999 -1.386394366
#define LOG_01_9999 -4.605270191
#define LOG_01_0001 -13.81551056
#define LOG_75_0001 -9.498022444
#define LOG_75_9999 -0.2877820775
#endif

#if 0
#define LOG_99_999 -0.01105083619
#define LOG_99_001 -6.917805615
#define LOG_25_001 -8.29404964
#define LOG_25_999 -1.387294861
#define LOG_01_999 -4.606170686
#define LOG_01_001 -11.51292546
#define LOG_75_001 -7.195437351
#define LOG_75_999 -0.2886825728
#endif

/* Switches on 4 consecutive mismatches */
#define LOG_99_99 -0.02010067171
#define LOG_99_01 -4.615220522
#define LOG_25_01 -5.991464547
#define LOG_25_99 -1.396344697
#define LOG_01_99 -4.615220522
#define LOG_01_01 -9.210340372
#define LOG_75_01 -4.892852258
#define LOG_75_99 -0.2977324083


static bool print_nsnpdiffs_p;
static bool print_snplabels_p;
static bool show_refdiff_p;

static IIT_T snps_iit;
static int *snps_divint_crosstable;

static IIT_T genes_iit;
static int *genes_divint_crosstable;

static IIT_T splicesites_iit;
static int *splicesites_divint_crosstable;

static int donor_typeint;
static int acceptor_typeint;

static int trim_mismatch_score;
static bool novelsplicingp;
static bool knownsplicingp;
static bool output_sam_p;

static Mode_T mode;

static double genomelength;	/* For BLAST E-value */


char *
Endtype_string (Endtype_T endtype) {
  switch (endtype) {
  case END: return "end";
  case INS: return "ins";
  case DEL: return "del";
  case DON: return "don";
  case ACC: return "acc";
  case AMB_DON: return "amb_don";
  case AMB_ACC: return "amb_acc";
  case TERM: return "term";
  default:
    fprintf(stderr,"Unexpected endtype %d\n",endtype);
    abort();
  }
  return "";
}

char *
Trimaction_string (Trimaction_T trimaction) {
  switch (trimaction) {
  case NO_TRIM: return "NO_TRIM";
  case PRE_TRIMMED: return "PRE_TRIMMED";
  case COMPUTE_TRIM: return "COMPUTE_TRIM";
  default:
    fprintf(stderr,"Unexpected trimaction %d\n",trimaction);
    abort();
  }
  return "";
}



static char complCode[128] = COMPLEMENT_LC;

#if 0
static char *
make_complement_buffered (char *complement, char *sequence, unsigned int length) {
  int i, j;

  /* complement = (char *) CALLOC_OUT(length+1,sizeof(char)); */
  for (i = length-1, j = 0; i >= 0; i--, j++) {
    complement[j] = complCode[(int) sequence[i]];
  }
  complement[length] = '\0';
  return complement;
}
#endif

static char *
make_complement_inplace (char *sequence, unsigned int length) {
  char temp;
  unsigned int i, j;

  for (i = 0, j = length-1; i < length/2; i++, j--) {
    temp = complCode[(int) sequence[i]];
    sequence[i] = complCode[(int) sequence[j]];
    sequence[j] = temp;
  }
  if (i == j) {
    sequence[i] = complCode[(int) sequence[i]];
  }

  return sequence;
}



#define T Substring_T

struct T {
  bool exactp;
  int nmismatches_whole;        /* Over total substring.  Comes from calling procedure */
  int nmismatches_bothdiff;	/* Over region left after trimming */
  int nmismatches_refdiff;	/* Over region left after trimming */
  /* nsnpdiffs = nmismatches_bothdiff - nmismatches_refdiff */

  int nmatches;			/* Not used anymore.  Use nmatches_posttrim instead. */

  int trim_left;
  int trim_right;
  int mandatory_trim_left;	/* Resulting from unalias, plus SOFT_CLIPS_AVOID_CIRCULARIZATION */
  int mandatory_trim_right;	/* Resulting from alias, plus SOFT_CLIPS_AVOID_CIRCULARIZATION */

  bool trim_left_splicep;
  bool trim_right_splicep;

  Chrnum_T chrnum;
  Univcoord_T chroffset;
  Univcoord_T chrhigh;
  Chrpos_T chrlength;

  Univcoord_T left; /* for plus: alignstart - querystart(orig).  for
		       minus: alignend - (querylength -
		       queryend(orig)).  Set when substring is created
		       or made unambiguous, and remains constant */

  Univcoord_T genomicstart;	/* For region corresponding to entire querylength (if extrapolated) */
  Univcoord_T genomicend;

  Endtype_T start_endtype;
  Endtype_T end_endtype;

  int querystart_orig;		/* For part that aligns to genome, pre-trim */
  int queryend_orig;
  int querystart;		/* For part that aligns to genome, post-trim */
  int queryend;
  int amb_splice_pos;		/* Used for ambiguous substrings */
  int querylength;

  Univcoord_T alignstart;	/* For part that aligns to genome, including part that is trimmed (pre-trim) */
  Univcoord_T alignend;

  Univcoord_T alignstart_trim;	/* For part that aligns to genome, excluding part that is trimmed (post-trim) */
  Univcoord_T alignend_trim;

  bool plusp;
  int genestrand;

  char *genomic_bothdiff; /* In same direction as query.  NULL if same
  			     as query.  Has dashes outside of
  			     querystart..(queryend-1) for indels and
  			     splices.  Has lowercase to indicate
  			     trimmed regions, terminal regions,
  			     deletions, splice dinucleotides, and
  			     mismatches from ref and alt.  Use for
  			     MAPQ computations and scoring. */

  char *genomic_refdiff;  /* Same as above, but lowercase for
			     mismatches from ref only.  For
			     non-SNP-tolerant alignment, this is just
			     a pointer to genomic_bothdiff.  Use for
			     NM and MD computations.  */

  float mapq_loglik;

  int sensedir;

  Univcoord_T splicecoord_D;
  int splicesitesD_knowni;	/* Needed for intragenic_splice_p in stage1hr.c */

  bool siteD_knownp;		/* Used for computing Substring_nchimera_known */
  bool siteD_novelp;
  int siteD_pos;
  double siteD_prob;

  /* for shortexon (always use *_1 for acceptor and *_2 for donor) */
  /* for donor/acceptor: the ambiguous position */
  Univcoord_T splicecoord_A;
  int splicesitesA_knowni;

  bool siteA_knownp;
  bool siteA_novelp;
  int siteA_pos;
  double siteA_prob;

  Univcoord_T splicecoord_N;	/* For DNA fusions */
  int siteN_pos;


  bool ambiguous_p;
  int nambcoords;
  Univcoord_T *ambcoords;
  int *amb_knowni;
  int *amb_nmismatches;
  double *amb_probs;
  Endtype_T amb_type;		/* Ambiguous DONs or ACCs */
};


void
Substring_alias_circular (T this) {
  Chrpos_T chrlength;

  if (this != NULL) {
    chrlength = this->chrlength;

    if (this->left + this->querylength > chrlength) {
      if (this->plusp == true) {
	this->mandatory_trim_right = chrlength - (this->left + this->querylength);
	debug2(printf("For alias, setting mandatory_trim_right to be %d\n",this->mandatory_trim_right));
      } else {
	this->mandatory_trim_left = chrlength - (this->left + this->querylength);
	debug2(printf("For alias, setting mandatory_trim_left to be %d\n",this->mandatory_trim_left));
      }
    }

    this->left += chrlength;
    this->genomicstart += chrlength;
    this->genomicend += chrlength;
    this->alignstart += chrlength;
    this->alignend += chrlength;
    this->alignstart_trim += chrlength;
    this->alignend_trim += chrlength;
    this->splicecoord_D += chrlength;
    this->splicecoord_A += chrlength;
    this->splicecoord_N += chrlength;
  }

  return;
}


void
Substring_unalias_circular (T this) {
  Chrpos_T chrlength;

  if (this != NULL) {
    chrlength = this->chrlength;

    if (this->left < chrlength) {
      if (this->plusp == true) {
	this->mandatory_trim_left = chrlength - this->left;
	debug2(printf("For alias, setting mandatory_trim_left to be %d\n",this->mandatory_trim_left));
      } else {
	this->mandatory_trim_right = chrlength - this->left;
	debug2(printf("For alias, setting mandatory_trim_right to be %d\n",this->mandatory_trim_right));
      }
    }

    this->left -= chrlength;
    this->genomicstart -= chrlength;
    this->genomicend -= chrlength;
    this->alignstart -= chrlength;
    this->alignend -= chrlength;
    this->alignstart_trim -= chrlength;
    this->alignend_trim -= chrlength;
    this->splicecoord_D -= chrlength;
    this->splicecoord_A -= chrlength;
    this->splicecoord_N -= chrlength;
  }

  return;
}



static void
fill_w_dashes (char *string, int start, int end) {
  int i;

  for (i = start; i < end; i++) {
    string[i] = '-';
  }
  return;
}


#if 0
static bool
dibase_mismatch_p (char leftnt, char leftcolor, char nt, char rightcolor, char rightnt) {
  char left_predict, right_predict;

  switch (toupper(leftnt)) {
  case 'A':
    switch (leftcolor) {
    case '0': left_predict = 'A'; break;
    case '1': left_predict = 'C'; break;
    case '2': left_predict = 'G'; break;
    case '3': left_predict = 'T'; break;
    }
    break;
  case 'C':
    switch (leftcolor) {
    case '0': left_predict = 'C'; break;
    case '1': left_predict = 'A'; break;
    case '2': left_predict = 'T'; break;
    case '3': left_predict = 'G'; break;
    }
    break;
  case 'G':
    switch (leftcolor) {
    case '0': left_predict = 'G'; break;
    case '1': left_predict = 'T'; break;
    case '2': left_predict = 'A'; break;
    case '3': left_predict = 'C'; break;
    }
    break;
  case 'T':
    switch (leftcolor) {
    case '0': left_predict = 'T'; break;
    case '1': left_predict = 'G'; break;
    case '2': left_predict = 'C'; break;
    case '3': left_predict = 'A'; break;
    }
    break;
  }

  switch (toupper(rightnt)) {
  case 'A':
    switch (rightcolor) {
    case '0': right_predict = 'A'; break;
    case '1': right_predict = 'C'; break;
    case '2': right_predict = 'G'; break;
    case '3': right_predict = 'T'; break;
    }
    break;
  case 'C':
    switch (rightcolor) {
    case '0': right_predict = 'C'; break;
    case '1': right_predict = 'A'; break;
    case '2': right_predict = 'T'; break;
    case '3': right_predict = 'G'; break;
    }
    break;
  case 'G':
    switch (rightcolor) {
    case '0': right_predict = 'G'; break;
    case '1': right_predict = 'T'; break;
    case '2': right_predict = 'A'; break;
    case '3': right_predict = 'C'; break;
    }
    break;
  case 'T':
    switch (rightcolor) {
    case '0': right_predict = 'T'; break;
    case '1': right_predict = 'G'; break;
    case '2': right_predict = 'C'; break;
    case '3': right_predict = 'A'; break;
    }
    break;
  }

  if (left_predict != right_predict) {
    /* printf("left_predict %c != right_predict %c\n",left_predict,right_predict); */
    return false;
  } else if (left_predict == toupper(nt)) {
    return false;
  } else {
    return true;
  }
}
#endif


#if 0
/* query is numbers */
static void
mark_mismatches_dibase (char *genome, char *query, int start, int end, bool plusp) {
  int i;

  debug1(printf("query:  %s\n",query));
  debug1(printf("genome: %s\n",genome));
  debug1(printf("mark:   "));

  if (plusp) {
    /* querypos start */
    debug1(printf("*"));
    if (dibase_mismatch_p(/*leftnt*/genome[start],/*leftcolor*/'0',genome[start],query[start],genome[start+1])) {
      genome[start] = (char) tolower(genome[start]);
    }

    for (i = start+1; i < end-1; i++) {
      debug1(printf("*"));

      if (dibase_mismatch_p(genome[i-1],query[i-1],genome[i],query[i],genome[i+1])) {
	/* printf("Got a dibase mismatch at %d with %c %c %c %c %c\n",
	   i,genome[i-1],query[i-1],genome[i],query[i],genome[i+1]); */
	genome[i] = (char) tolower(genome[i]);
      }
    }

    /* querypos (end-1) */
    if (dibase_mismatch_p(genome[end-2],query[end-2],genome[end-1],
			  /*rightcolor*/'0',/*rightnt*/genome[end-1])) {
      genome[end-1] = (char) tolower(genome[end-1]);
    }

  } else {
    /* querypos start */
    debug1(printf("*"));
    if (dibase_mismatch_p(/*leftnt*/genome[start],/*leftcolor*/'0',genome[start],query[start+1],genome[start+1])) {
      genome[start] = (char) tolower(genome[start]);
    }

    for (i = start+1; i < end-1; i++) {
      debug1(printf("*"));

      if (dibase_mismatch_p(genome[i-1],query[i],genome[i],query[i+1],genome[i+1])) {
	/* printf("Got a dibase mismatch at %d with %c %c %c %c %c\n",
	   i,genome[i-1],query[i],genome[i],query[i+1],genome[i+1]); */
	genome[i] = (char) tolower(genome[i]);
      }
    }

    /* querypos (end-1) */
    if (dibase_mismatch_p(genome[end-2],query[end-1],genome[end-1],
			  /*rightcolor*/'0',/*rightnt*/genome[end-1])) {
      genome[end-1] = (char) tolower(genome[end-1]);
    }
  }

  debug1(printf("\n"));

  return;
}
#endif


#if 0
static int
count_mismatches_dibase (char *genome, char *query, int start, int end, bool plusp) {
  int nmismatches = 0, i;

  debug1(printf("query:  %s\n",query));
  debug1(printf("genome: %s\n",genome));
  debug1(printf("mark:   "));

  if (plusp) {
    /* querypos start */
    debug1(printf("*"));
    if (dibase_mismatch_p(/*leftnt*/genome[start],/*leftcolor*/'0',genome[start],query[start],genome[start+1])) {
      nmismatches++;
    }

    for (i = start+1; i < end-1; i++) {
      debug1(printf("*"));

      if (dibase_mismatch_p(genome[i-1],query[i-1],genome[i],query[i],genome[i+1])) {
	/* printf("Got a dibase mismatch at %d with %c %c %c %c %c\n",
	   i,genome[i-1],query[i-1],genome[i],query[i],genome[i+1]); */
	nmismatches++;
      }
    }

    /* querypos (end-1) */
    if (dibase_mismatch_p(genome[end-2],query[end-2],genome[end-1],
			  /*rightcolor*/'0',/*rightnt*/genome[end-1])) {
      nmismatches++;
    }

  } else {
    /* querypos start */
    debug1(printf("*"));
    if (dibase_mismatch_p(/*leftnt*/genome[start],/*leftcolor*/'0',genome[start],query[start+1],genome[start+1])) {
      nmismatches++;
    }

    for (i = start+1; i < end-1; i++) {
      debug1(printf("*"));

      if (dibase_mismatch_p(genome[i-1],query[i],genome[i],query[i+1],genome[i+1])) {
	/* printf("Got a dibase mismatch at %d with %c %c %c %c %c\n",
	   i,genome[i-1],query[i],genome[i],query[i+1],genome[i+1]); */
	nmismatches++;
      }
    }

    /* querypos (end-1) */
    if (dibase_mismatch_p(genome[end-2],query[end-1],genome[end-1],
			  /*rightcolor*/'0',/*rightnt*/genome[end-1])) {
      nmismatches++;
    }
  }

  debug1(printf("\n"));
  debug1(printf("%d mismatches\n\n",nmismatches));

  return nmismatches;
}
#endif


static int
trim_left_end (int *nmismatches_end, Compress_T query_compress, Univcoord_T left, int querystart, int queryend,
	       int querylength, bool plusp, int genestrand, int trim_mismatch_score) {
  int bestscore, score;
  int trim5, alignlength, pos, prevpos, i;
  int nmismatches;

#if defined(LONG_READLENGTHS)
  int *mismatch_positions = (int *) MALLOC(querylength*sizeof(int));
#elif defined(HAVE_ALLOCA)
  int *mismatch_positions = (int *) ALLOCA(querylength*sizeof(int));
#else
  int mismatch_positions[MAX_READLENGTH];
#endif


  debug8(printf("Entered trim_left_end with querystart %d, queryend %d\n",querystart,queryend));

  alignlength = queryend - querystart;
  bestscore = 0;
  score = 0;
  trim5 = 0;

  if (plusp == true) {
    debug8(printf("Calling Genome_mismatches_right_trim with left %u, pos5 %d, pos3 %d\n",
		  left,querystart,queryend));
    nmismatches = Genome_mismatches_right_trim(mismatch_positions,/*max_mismatches*/alignlength,
					       query_compress,left,/*pos5*/querystart,/*pos3*/queryend,
					       plusp,genestrand);
    debug8(printf("%d mismatches:",nmismatches));
    debug8(
	   for (i = 0; i < nmismatches; i++) {
	     printf(" %d",mismatch_positions[i]);
	   }
	   printf("\n");
	   );

    prevpos = queryend;
    for (i = 0; i < nmismatches; i++) {
      pos = mismatch_positions[i] + 1; /* the position just after the mismatch */
      score += (prevpos - pos)*TRIM_MATCH_SCORE;
      if (score >= bestscore) {	/* want >= and not >, so ties go to end */
	bestscore = score;
	trim5 = pos;
      }
      debug8(printf("Trim left pos %d, score %d, trim5 %d\n",pos,score,trim5));

      score += trim_mismatch_score; /* For the mismatch */
      if (score < 0) {
	score = 0;
      }
      prevpos = pos - 1;		/* On the mismatch */
    }

    score += prevpos*TRIM_MATCH_SCORE;
    if (score >= bestscore) {	/* want >= and not >, so ties go to end */
      /* bestscore = score; */
      trim5 = 0;
    }
    
    *nmismatches_end = 0;
    for (i = 0; i < nmismatches; i++) {
      pos = mismatch_positions[i]; /* the position of the mismatch */
      if (mismatch_positions[i] < trim5) {
	*nmismatches_end += 1;
      }
    }

  } else {
    debug8(printf("Calling Genome_mismatches_left_trim with left %u, pos5 %d - %d, pos3 %d - %d\n",
		  left,querylength,queryend,querylength,querystart));
    nmismatches = Genome_mismatches_left_trim(mismatch_positions,/*max_mismatches*/alignlength,
					      query_compress,left,/*pos5*/querylength - queryend,
					      /*pos3*/querylength - querystart,plusp,genestrand);

    debug8(printf("%d mismatches:",nmismatches));
    debug8(
	   for (i = 0; i < nmismatches; i++) {
	     printf(" %d",querylength - 1 - mismatch_positions[i]);
	   }
	   printf("\n");
	   );

    prevpos = queryend;
    for (i = 0; i < nmismatches; i++) {
      pos = (querylength - 1 - mismatch_positions[i]) + 1; /* the position just after the mismatch */
      score += (prevpos - pos)*TRIM_MATCH_SCORE;
      if (score >= bestscore) {	/* want >= and not >, so ties go to end */
	bestscore = score;
	trim5 = pos;
      }
      debug8(printf("Trim left pos %d, score %d, trim5 %d\n",pos,score,trim5));

      score += trim_mismatch_score; /* For the mismatch */
      if (score < 0) {
	score = 0;
      }
      prevpos = pos - 1;		/* On the mismatch */
    }

    score += prevpos*TRIM_MATCH_SCORE;
    if (score >= bestscore) {	/* want >= and not >, so ties go to end */
      /* bestscore = score; */
      trim5 = 0;
    }
    
    *nmismatches_end = 0;
    for (i = 0; i < nmismatches; i++) {
      pos = querylength - 1 - mismatch_positions[i]; /* the position of the mismatch */
      if (pos < trim5) {
	*nmismatches_end += 1;
      }
    }
  }

#if defined(LONG_READLENGTHS)
  FREE(mismatch_positions);
#elif defined(HAVE_ALLOCA)
  FREEA(mismatch_positions);
#else
  /* Hard-coded use of MAX_READLENGTH */
#endif

  debug8(printf("Trim left pos 0, score %d, trim5 %d, nmismatches_end %d\n",score,trim5,*nmismatches_end));
  debug8(printf("\n"));

  return trim5;
}



/* querystart (typically 0) and queryend (typically querylength) are exclusive */
/* sequences may have had lower case characters marked */
static int
trim_right_end (int *nmismatches_end, Compress_T query_compress, Univcoord_T left, int querystart, int queryend,
		int querylength, bool plusp, int genestrand, int trim_mismatch_score) {
  int bestscore, score;
  int trim3, alignlength, pos, prevpos, i;
  int nmismatches;

#if defined(LONG_READLENGTHS)
  int *mismatch_positions = (int *) MALLOC(querylength*sizeof(int));
#elif defined(HAVE_ALLOCA)
  int *mismatch_positions = (int *) ALLOCA(querylength*sizeof(int));
#else
  int mismatch_positions[MAX_READLENGTH];
#endif

  debug8(printf("Entered trim_right_end with querystart %d, queryend %d\n",querystart,queryend));

  alignlength = queryend - querystart;
  bestscore = 0;
  score = 0;
  trim3 = 0;

  if (plusp == true) {
    nmismatches = Genome_mismatches_left_trim(mismatch_positions,/*max_mismatches*/alignlength,
					      query_compress,left,/*pos5*/querystart,/*pos3*/queryend,
					      plusp,genestrand);

    debug8(printf("%d mismatches:",nmismatches));
    debug8(
	   for (i = 0; i < nmismatches; i++) {
	     printf(" %d",mismatch_positions[i]);
	   }
	   printf("\n");
	   );

    prevpos = querystart - 1;
    for (i = 0; i < nmismatches; i++) {
      pos = mismatch_positions[i] - 1; /* the position just before the mismatch */
      score += (pos - prevpos)*TRIM_MATCH_SCORE;
      if (score >= bestscore) {	/* want >= and not >, so ties go to end */
	bestscore = score;
	trim3 = querylength - pos - 1;
      }
      debug8(printf("Trim right pos %d, score %d, trim3 %d\n",pos,score,trim3));
      
      score += trim_mismatch_score; /* For the mismatch */
      if (score < 0) {
	score = 0;
      }
      prevpos = pos + 1;		/* On the mismatch */
    }

    score += (queryend - 1 - prevpos)*TRIM_MATCH_SCORE;
    if (score >= bestscore) { /* want >= and not >, so ties go to end */
      /* bestscore = score; */
      trim3 = 0;
    }

    *nmismatches_end = 0;
    for (i = 0; i < nmismatches; i++) {
      pos = mismatch_positions[i]; /* the position of the mismatch */
      if (pos > querylength - 1 - trim3) {
	*nmismatches_end += 1;
      }
    }

  } else {
    nmismatches = Genome_mismatches_right_trim(mismatch_positions,/*max_mismatches*/alignlength,
					       query_compress,left,/*pos5*/querylength - queryend,
					       /*pos3*/querylength - querystart,plusp,genestrand);

    debug8(printf("%d mismatches:",nmismatches));
    debug8(
	   for (i = 0; i < nmismatches; i++) {
	     printf(" %d",querylength - 1 - mismatch_positions[i]);
	   }
	   printf("\n");
	   );

    prevpos = querystart - 1;
    for (i = 0; i < nmismatches; i++) {
      pos = (querylength - 1 - mismatch_positions[i]) - 1; /* the position just before the mismatch */
      score += (pos - prevpos)*TRIM_MATCH_SCORE;
      if (score >= bestscore) {	/* want >= and not >, so ties go to end */
	bestscore = score;
	trim3 = querylength - pos - 1;
      }
      debug8(printf("Trim right pos %d, score %d, trim3 %d\n",pos,score,trim3));
      
      score += trim_mismatch_score; /* For the mismatch */
      if (score < 0) {
	score = 0;
      }
      prevpos = pos + 1;		/* On the mismatch */
    }

    score += (queryend - 1 - prevpos)*TRIM_MATCH_SCORE;
    if (score >= bestscore) { /* want >= and not >, so ties go to end */
      /* bestscore = score; */
      trim3 = 0;
    }

    *nmismatches_end = 0;
    for (i = 0; i < nmismatches; i++) {
      pos = querylength - 1 - mismatch_positions[i]; /* the position of the mismatch */
      if (pos > querylength - 1 - trim3) {
	*nmismatches_end += 1;
      }
    }
  }

#if defined(LONG_READLENGTHS)
  FREE(mismatch_positions);
#elif defined(HAVE_ALLOCA)
  FREEA(mismatch_positions);
#else
  /* Hard-coded use of MAX_READLENGTH */
#endif

  debug8(printf("Trim right pos %d, score %d, trim3 %d, nmismatches_end %d\n",queryend-1,score,trim3,*nmismatches_end));
  debug8(printf("\n"));

  return trim3;
}


#if 0
/* Uses hmm */
/* Modified for Stage3_bad_stretch_p */
bool
Substring_bad_stretch_p (T this, Compress_T query_compress_fwd, Compress_T query_compress_rev) {
  int alignlength, startpos, endpos, pos, i;
  float vprob_good, vprob_bad, prev_vprob_good, prev_vprob_bad, good_incr_prob, bad_incr_prob;
#ifdef DEBUG9
  bool result;
  State_T vstate_good[MAX_READLENGTH], vstate_bad[MAX_READLENGTH], state;
#endif

#ifdef HAVE_ALLOCA
  int *mismatch_positions;
#else
  int mismatch_positions[MAX_READLENGTH];
#endif


  debug9(printf("Entered bad_stretch_p with querystart_orig %d, queryend_orig %d, plusp %d\n",
		this->querystart_orig,this->queryend_orig,this->plusp));

  alignlength = this->queryend_orig - this->querystart_orig;
#ifdef HAVE_ALLOCA
  mismatch_positions = (int *) ALLOCA(alignlength*sizeof(int));
#endif

  if (this->plusp == true) {
    startpos = this->querystart_orig;
    endpos = this->queryend_orig;
    debug9(printf("Calling Genome_mismatches_left from pos5 %d to pos3 %d\n",startpos,endpos));
#ifdef DEBUG9
    nmismatches = Genome_mismatches_left(mismatch_positions,/*max_mismatches*/alignlength,
					 query_compress_fwd,this->left,/*pos5*/startpos,/*pos3*/endpos,
					 /*plusp*/true,this->genestrand);
#endif
  } else {
    startpos = this->querylength - this->queryend_orig;
    endpos = this->querylength - this->querystart_orig;
    debug9(printf("Calling Genome_mismatches_left from pos5 %d to pos3 %d\n",startpos,endpos));
#ifdef DEBUG9
    nmismatches = Genome_mismatches_left(mismatch_positions,/*max_mismatches*/alignlength,
					 query_compress_rev,this->left,/*pos5*/startpos,/*pos3*/endpos,
					 /*plusp*/false,this->genestrand);
#endif
  }

  debug9(printf("%d mismatches:",nmismatches));
  debug9(
	 for (i = 0; i < nmismatches; i++) {
	   printf(" %d",mismatch_positions[i]);
	 }
	 printf("\n");
	 );

  /* Initialize priors */
#ifdef COMPUTE_LOG
  prev_vprob_good = log(0.99);
  prev_vprob_bad = log(0.01);
#else
  prev_vprob_good = LOG_99;
  prev_vprob_bad = LOG_01;
#endif

  i = 0;
  for (pos = startpos; pos < endpos; pos++) {
    if (pos != mismatch_positions[i]) {
      /* match */
#ifdef DEBUG9
      printf("querypos %d (match): ",pos);
#endif

      /* state: GOOD */
#ifdef COMPUTE_LOG
      good_incr_prob = log(/*emission_prob*/0.99) + log(/*transition_prob*/0.99);
      bad_incr_prob = log(/*emission_prob*/0.99) + log(/*transition_prob*/0.01);
#else
      good_incr_prob = LOG_99_99;
      bad_incr_prob = LOG_99_01;
#endif

#ifdef DEBUG9
      if (prev_vprob_good + good_incr_prob > prev_vprob_bad + bad_incr_prob) {
	printf("%.2f + %.2f (GOOD->GOOD*) or %.2f + %.2f (BAD->GOOD ).  ",
	       prev_vprob_good,good_incr_prob,prev_vprob_bad,bad_incr_prob);
	vprob_good = prev_vprob_good + good_incr_prob;
	vstate_good[pos] = GOOD;
      } else {
	printf("%.2f + %.2f (GOOD->GOOD ) or %.2f + %.2f (BAD->GOOD!).  ",
	       prev_vprob_good,good_incr_prob,prev_vprob_bad,bad_incr_prob);
	vprob_good = prev_vprob_bad + bad_incr_prob;
	vstate_good[pos] = BAD;
      }
#else
      if (prev_vprob_good + good_incr_prob > prev_vprob_bad + bad_incr_prob) {
	vprob_good = prev_vprob_good + good_incr_prob;
      } else {
	/* vprob_good = prev_vprob_bad + bad_incr_prob; */
	return true;
      }
#endif

      /* state: BAD */
#ifdef COMPUTE_LOG
      good_incr_prob = log(/*emission_prob*/0.25) + log(/*transition_prob*/0.01);
      bad_incr_prob = log(/*emission_prob*/0.25) + log(/*transition_prob*/0.99);
#else
      good_incr_prob = LOG_25_01;
      bad_incr_prob = LOG_25_99;
#endif
      if (prev_vprob_good + good_incr_prob > prev_vprob_bad + bad_incr_prob) {
	vprob_bad = prev_vprob_good + good_incr_prob;
#ifdef DEBUG9
	printf("%.2f + %.2f (GOOD->BAD*) or %.2f + %.2f (BAD->BAD ).\n",
	       prev_vprob_good,good_incr_prob,prev_vprob_bad,bad_incr_prob);
	vstate_bad[pos] = GOOD;
#endif
      } else {
	vprob_bad = prev_vprob_bad + bad_incr_prob;
#ifdef DEBUG9
	printf("%.2f + %.2f (GOOD->BAD ) or %.2f + %.2f (BAD->BAD*).\n",
	       prev_vprob_good,good_incr_prob,prev_vprob_bad,bad_incr_prob);
	vstate_bad[pos] = BAD;
#endif
      }

    } else {
      /* mismatch */
#ifdef DEBUG9
      printf("querypos %d (mismatch): ",pos);
#endif

      /* state: GOOD */
#ifdef COMPUTE_LOG
      good_incr_prob = log(/*emission_prob*/0.01) + log(/*transition_prob*/0.99);
      bad_incr_prob = log(/*emission_prob*/0.01) + log(/*transition_prob*/0.01);
#else
      good_incr_prob = LOG_01_99;
      bad_incr_prob = LOG_01_01;
#endif
#ifdef DEBUG9
      if (prev_vprob_good + good_incr_prob > prev_vprob_bad + bad_incr_prob) {
	printf("%.2f + %.2f (GOOD->GOOD*) or %.2f + %.2f (BAD->GOOD ).  ",
	       prev_vprob_good,good_incr_prob,prev_vprob_bad,bad_incr_prob);
	vprob_good = prev_vprob_good + good_incr_prob;
	vstate_good[pos] = GOOD;
      } else {
	printf("%.2f + %.2f (GOOD->GOOD ) or %.2f + %.2f (BAD->GOOD!).  ",
	       prev_vprob_good,good_incr_prob,prev_vprob_bad,bad_incr_prob);
	vprob_good = prev_vprob_bad + bad_incr_prob;
	vstate_good[pos] = BAD;
      }
#else
      if (prev_vprob_good + good_incr_prob > prev_vprob_bad + bad_incr_prob) {
	vprob_good = prev_vprob_good + good_incr_prob;
      } else {
	/* vprob_good = prev_vprob_bad + bad_incr_prob; */
	return true;
      }
#endif

      /* state: BAD */
#ifdef COMPUTE_LOG
      good_incr_prob = log(/*emission_prob*/0.75) + log(/*transition_prob*/0.01);
      bad_incr_prob = log(/*emission_prob*/0.75) + log(/*transition_prob*/0.99);
#else
      good_incr_prob = LOG_75_01;
      bad_incr_prob = LOG_75_99;
#endif
      if (prev_vprob_good + good_incr_prob > prev_vprob_bad + bad_incr_prob) {
	vprob_bad = prev_vprob_good + good_incr_prob;
#ifdef DEBUG9
	printf("%.2f + %.2f (GOOD->BAD*) or %.2f + %.2f (BAD->BAD ).\n",
	       prev_vprob_good,good_incr_prob,prev_vprob_bad,bad_incr_prob);
	vstate_bad[pos] = GOOD;
#endif
      } else {
	vprob_bad = prev_vprob_bad + bad_incr_prob;
#ifdef DEBUG9
	printf("%.2f + %.2f (GOOD->BAD ) or %.2f + %.2f (BAD->BAD*).\n",
	       prev_vprob_good,good_incr_prob,prev_vprob_bad,bad_incr_prob);
	vstate_bad[pos] = BAD;
#endif
      }

      i++;
    }
    
    prev_vprob_good = vprob_good;
    prev_vprob_bad = vprob_bad;
  }

#ifdef DEBUG9
  result = false;
  if (prev_vprob_good > prev_vprob_bad) {
    printf("Final state: GOOD\n");
    state = GOOD;
  } else {
    printf("Final state: BAD\n");
    state = BAD;
    result = true;
  }

  for (pos = endpos-1; pos >= startpos; pos--) {
    if (state == GOOD) {
      printf("querypos %d: GOOD\n",pos);
      state = vstate_good[pos];
    } else {
      printf("querypos %d: BAD\n",pos);
      state = vstate_bad[pos];
      result = true;
    }
  }
  return result;

#else
  return false;
#endif

}
#endif


void
Substring_free (T *old) {
  if (*old) {
    debug2(printf("Freeing substring %p\n",*old));
    if ((*old)->nambcoords > 0) {
      FREE_OUT((*old)->ambcoords);
      FREE_OUT((*old)->amb_knowni);
      FREE_OUT((*old)->amb_nmismatches);
      FREE_OUT((*old)->amb_probs);
    }
    if ((*old)->genomic_bothdiff != NULL) {
      if ((*old)->genomic_refdiff != (*old)->genomic_bothdiff) {
	FREE_OUT((*old)->genomic_refdiff);
      }
      FREE_OUT((*old)->genomic_bothdiff);
    }

    FREE_OUT(*old);
  }
  return;
}


bool
Substring_contains_p (T this, int querypos) {
  if (querypos >= this->querystart && querypos < this->queryend) {
    return true;
  } else {
    return false;
  }
}


void
Substring_print_ends (T this, int chrnum) {
  if (this == NULL) {
    printf("NA ");
  } else {
    printf("#%d:%llu..%llu ",chrnum,(unsigned long long) this->alignstart,(unsigned long long) this->alignend);
  }
  return;
}


int
Substring_compare (T substring1, T substring2, int alias1, int alias2, Chrpos_T chrlength1, Chrpos_T chrlength2) {
  Univcoord_T alignstart1, alignend1, alignstart2, alignend2;

  if (substring1 == NULL && substring2 == NULL) {
    return 0;
  } else if (substring1 == NULL) {
    return -1;
  } else if (substring2 == NULL) {
    return +1;
  } else {
    alignstart1 = substring1->alignstart;
    alignend1 = substring1->alignend;
    if (alias1 < 0) {
      alignstart1 += chrlength1;
      alignend1 += chrlength1;
    }

    alignstart2 = substring2->alignstart;
    alignend2 = substring2->alignend;
    if (alias2 < 0) {
      alignstart2 += chrlength2;
      alignend2 += chrlength2;
    }

    if (alignstart1 < alignstart2) {
      return -1;
    } else if (alignstart1 > alignstart2) {
      return +1;
    } else if (alignend1 < alignend2) {
      return -1;
    } else if (alignend1 > alignend2) {
      return +1;
    } else {
      return 0;
    }
  }
}


#if 0
bool
Substring_equal_p (T substring1, T substring2) {

  if (substring1 == NULL && substring2 == NULL) {
    return true;
  } else if (substring1 == NULL) {
    return false;
  } else if (substring2 == NULL) {
    return false;
  } else if (substring1->alignstart != substring2->alignstart) {
    return false;
  } else if (substring1->alignend != substring2->alignend) {
    return false;
  } else {
    return true;
  }
}
#endif


bool
Substring_overlap_p (T substring1, T substring2) {
  Univcoord_T low1, high1, low2, high2;

  if (substring1->plusp == true) {
    low1 = substring1->alignstart;
    high1 = substring1->alignend;
  } else {
    low1 = substring1->alignend;
    high1 = substring1->alignstart;
  }

  if (substring2->plusp == true) {
    low2 = substring2->alignstart;
    high2 = substring2->alignend;
  } else {
    low2 = substring2->alignend;
    high2 = substring2->alignstart;
  }

  debug3(printf("Checking overlap between %u..%u and %u..%u",low1,high1,low2,high2));

  if (high2 < low1) {
    debug3(printf(" => no because %u < %u\n",high2,low1));
    return false;
  } else if (low2 > high1) {
    debug3(printf(" => no because %u > %u\n",low2,high1));
    return false;
  } else {
    debug3(printf(" => yes\n"));
    return true;
  }
}


Chrpos_T
Substring_insert_length (T substring5, T substring3) {
  Univcoord_T pos5, pos3;

  pos5 = substring5->genomicstart;
  debug3(printf("pos5 %u\n",substring5->genomicstart));

  pos3 = substring3->genomicend;
  debug3(printf("pos3 %u\n",substring3->genomicend));

  if (pos5 > pos3) {
    return pos5 - pos3;
  } else {
    return pos3 - pos5;
  }
}


bool
Substring_overlap_point_trimmed_p (T substring, Univcoord_T endpos) {
  Univcoord_T low, high;

  if (substring == NULL) {
    return false;
  }

  if (substring->plusp == true) {
    low = substring->alignstart_trim;
    high = substring->alignend_trim;
    debug3(printf("Checking overlap between plus %u..%u and %u",low,high,endpos));
  } else {
    low = substring->alignend_trim;
    high = substring->alignstart_trim;
    debug3(printf("Checking overlap between minus %u..%u and %u",low,high,endpos));
  }


  if (endpos < low) {
    debug3(printf(" => no because %u < %u\n",endpos,low));
    return false;
  } else if (endpos > high) {
    debug3(printf(" => no because %u > %u\n",endpos,high));
    return false;
  } else {
    debug3(printf(" => yes\n"));
    return true;
  }
}


Univcoord_T
Substring_overlap_segment_trimmed (T substring1, T substring2) {
  Univcoord_T maxlow, minhigh;
  Univcoord_T low1, high1, low2, high2;

  if (substring1->plusp == true) {
    low1 = substring1->alignstart_trim;
    high1 = substring1->alignend_trim;
  } else {
    low1 = substring1->alignend_trim;
    high1 = substring1->alignstart_trim;
  }

  if (substring2->plusp == true) {
    low2 = substring2->alignstart_trim;
    high2 = substring2->alignend_trim;
  } else {
    low2 = substring2->alignend_trim;
    high2 = substring2->alignstart_trim;
  }

  debug3(printf("Checking overlap between %u..%u and %u..%u",low1,high1,low2,high2));

  if (high2 < low1) {
    debug3(printf(" => no because %u < %u\n",high2,low1));
    return 0;
  } else if (low2 > high1) {
    debug3(printf(" => no because %u > %u\n",low2,high1));
    return 0;
  } else {
    maxlow = (low1 > low2) ? low1 : low2;
    minhigh = (high1 < high2) ? high1 : high2;
    debug3(printf(" => yes.  maxlow %llu, minhigh %llu.  returning %llu\n",
		  maxlow,minhigh,maxlow + (minhigh - maxlow)/2));
    return maxlow + (minhigh - maxlow)/2;
  }
}



#if 0
Chrpos_T
Substring_insert_length_trimmed (T substring5, T substring3, int hit5_trim_left, int hit5_trim_right, int hit3_trim_left, int hit3_trim_right) {
  Univcoord_T rightpos, leftpos, pos;
  int right_extent, left_extent;
  
  debug3(printf("substring5 %d..%d out of %d.  aligned to %u..%u\n",
		substring5->querystart,substring5->queryend,substring5->querylength,
		substring5->alignstart_trim - substring5->chroffset,substring5->alignend_trim - substring5->chroffset));
  debug3(printf("substring3 %d..%d out of %d.  aligned to %u..%u\n",
		substring3->querystart,substring3->queryend,substring3->querylength,
		substring3->alignstart_trim - substring3->chroffset,substring3->alignend_trim - substring3->chroffset));
  if (substring5->plusp == true) {
    if (substring5->alignend_trim > substring3->alignend_trim) {
      rightpos = substring5->alignend_trim;
      right_extent = substring5->querylength - substring5->queryend - hit5_trim_right;
    } else {
      rightpos = substring3->alignend_trim;
      right_extent = substring3->querylength - substring3->queryend - hit3_trim_right;
    }
    debug3(printf("rightpos = %u, right_extent = %d\n",rightpos,right_extent));

    if (substring5->alignstart_trim < substring3->alignstart_trim) {
      leftpos = substring5->alignstart_trim;
      left_extent = substring5->querystart - hit5_trim_left;
    } else {
      leftpos = substring3->alignstart_trim;
      left_extent = substring3->querystart - hit3_trim_left;
    }
    debug3(printf("leftpos = %u, left_extent = %d\n",leftpos,left_extent));

    return (rightpos - leftpos) + right_extent + left_extent;

  } else {
    if (substring3->alignstart_trim > (pos = substring5->alignend_trim + (substring5->queryend - substring5->querystart /* - hit5_trim_left */) +
				       (substring3->querylength - substring3->querystart - hit3_trim_right))) {
      return 0;
    } else {
      debug3(printf("minus: substring5->alignend_trim %u - substring3->alignstart_trim %u +\n",
		    substring5->alignend_trim - substring5->chroffset,substring3->alignstart_trim - substring3->chroffset));
      debug3(printf("      substring5->queryend %d - substring5->querystart %d +\n",substring5->queryend,substring5->querystart));
      debug3(printf("      substring3->querylength %d - substring3->querystart %d - hit3_trim_right %d\n",
		    substring3->querylength,substring3->querystart,hit3_trim_right));
      return (pos - substring3->alignstart_trim);
    }
  }
}
#endif



static void
mark_mismatches_cmet_gsnap (char *gbuffer, char *query, int start, int end, int genestrand) {
  int i;
  
  debug1(printf("\n"));
  debug1(printf("query:  %s\n",query));
  debug1(printf("genome: %s\n",gbuffer));
  debug1(printf("count:  "));

  if (genestrand == +2) {
    for (i = start; i < end; i++) {
      if (gbuffer[i] == 'G' && query[i] == 'A') {
	debug1(printf("."));
	gbuffer[i] = '.';
      } else if (query[i] != gbuffer[i]) {
	debug1(printf("x"));
	assert(gbuffer[i] != OUTOFBOUNDS);
	gbuffer[i] = (char) tolower(gbuffer[i]);
      } else {
	debug1(printf("*"));
      }
    }

  } else {
    for (i = start; i < end; i++) {
      if (gbuffer[i] == 'C' && query[i] == 'T') {
	debug1(printf("."));
	gbuffer[i] = '.';
      } else if (query[i] != gbuffer[i]) {
	debug1(printf("x"));
	assert(gbuffer[i] != OUTOFBOUNDS);
	gbuffer[i] = (char) tolower(gbuffer[i]);
      } else {
	debug1(printf("*"));
      }
    }
  }

  return;
}


static void
mark_mismatches_cmet_sam (char *gbuffer, char *query, int start, int end, int genestrand) {
  int i;
  
  debug1(printf("query:  %s\n",query));
  debug1(printf("genome: %s\n",gbuffer));
  debug1(printf("count:  "));

  if (genestrand == +2) {
    for (i = start; i < end; i++) {
      if (gbuffer[i] == 'G' && query[i] == 'A') {
	debug1(printf("."));
#if 0
	/* Want to show mismatches */
	gbuffer[i] = 'A';		/* Avoids showing mismatches in MD and NM strings */
#endif
      } else if (query[i] != gbuffer[i]) {
	debug1(printf("x"));
	assert(gbuffer[i] != OUTOFBOUNDS);
	gbuffer[i] = (char) tolower(gbuffer[i]);
      } else {
	debug1(printf("*"));
      }
    }

  } else {
    for (i = start; i < end; i++) {
      if (gbuffer[i] == 'C' && query[i] == 'T') {
	debug1(printf("."));
#if 0
	/* Want to show mismatches */
	gbuffer[i] = 'T';		/* Avoids showing mismatches in MD and NM strings */
#endif
      } else if (query[i] != gbuffer[i]) {
	debug1(printf("x"));
	assert(gbuffer[i] != OUTOFBOUNDS);
	gbuffer[i] = (char) tolower(gbuffer[i]);
      } else {
	debug1(printf("*"));
      }
    }
  }
  
  return;
}



static void
mark_mismatches_atoi_gsnap (char *gbuffer, char *query, int start, int end, int genestrand) {
  int i;
  
  debug1(printf("query:  %s\n",query));
  debug1(printf("genome: %s\n",gbuffer));
  debug1(printf("count:  "));

  if (genestrand == +2) {
    for (i = start; i < end; i++) {
      if (gbuffer[i] == 'T' && query[i] == 'C') {
	debug1(printf("."));
	gbuffer[i] = '.';
      } else if (query[i] != gbuffer[i]) {
	debug1(printf("x"));
	assert(gbuffer[i] != OUTOFBOUNDS);
	gbuffer[i] = (char) tolower(gbuffer[i]);
      } else {
	debug1(printf("*"));
      }
    }

  } else {
    for (i = start; i < end; i++) {
      if (gbuffer[i] == 'A' && query[i] == 'G') {
	debug1(printf("."));
	gbuffer[i] = '.';
      } else if (query[i] != gbuffer[i]) {
	debug1(printf("x"));
	assert(gbuffer[i] != OUTOFBOUNDS);
	gbuffer[i] = (char) tolower(gbuffer[i]);
      } else {
	debug1(printf("*"));
      }
    }
  }
  
  return;
}



static void
mark_mismatches_atoi_sam (char *gbuffer, char *query, int start, int end, int genestrand) {
  int i;
  
  debug1(printf("query:  %s\n",query));
  debug1(printf("genome: %s\n",gbuffer));
  debug1(printf("count:  "));

  if (genestrand == +2) {
    for (i = start; i < end; i++) {
      if (gbuffer[i] == 'T' && query[i] == 'C') {
	debug1(printf("."));
#if 0
	/* Want to show mismatches */
	gbuffer[i] = 'C';		/* Avoids showing mismatches in MD and NM strings */
#endif
      } else if (query[i] != gbuffer[i]) {
	debug1(printf("x"));
	assert(gbuffer[i] != OUTOFBOUNDS);
	gbuffer[i] = (char) tolower(gbuffer[i]);
      } else {
	debug1(printf("*"));
      }
    }

  } else {
    for (i = start; i < end; i++) {
      if (gbuffer[i] == 'A' && query[i] == 'G') {
	debug1(printf("."));
#if 0
	/* Want to show mismatches */
	gbuffer[i] = 'G';		/* Avoids showing mismatches in MD and NM strings */
#endif
      } else if (query[i] != gbuffer[i]) {
	debug1(printf("x"));
	assert(gbuffer[i] != OUTOFBOUNDS);
	gbuffer[i] = (char) tolower(gbuffer[i]);
      } else {
	debug1(printf("*"));
      }
    }
  }

  return;
}


static void
mark_mismatches_ttoc_gsnap (char *gbuffer, char *query, int start, int end, int genestrand) {
  int i;
  
  debug1(printf("query:  %s\n",query));
  debug1(printf("genome: %s\n",gbuffer));
  debug1(printf("count:  "));

  if (genestrand == +2) {
    for (i = start; i < end; i++) {
      if (gbuffer[i] == 'A' && query[i] == 'G') {
	debug1(printf("."));
	gbuffer[i] = '.';
      } else if (query[i] != gbuffer[i]) {
	debug1(printf("x"));
	assert(gbuffer[i] != OUTOFBOUNDS);
	gbuffer[i] = (char) tolower(gbuffer[i]);
      } else {
	debug1(printf("*"));
      }
    }

  } else {
    for (i = start; i < end; i++) {
      if (gbuffer[i] == 'T' && query[i] == 'C') {
	debug1(printf("."));
	gbuffer[i] = '.';
      } else if (query[i] != gbuffer[i]) {
	debug1(printf("x"));
	assert(gbuffer[i] != OUTOFBOUNDS);
	gbuffer[i] = (char) tolower(gbuffer[i]);
      } else {
	debug1(printf("*"));
      }
    }
  }
  
  return;
}



static void
mark_mismatches_ttoc_sam (char *gbuffer, char *query, int start, int end, int genestrand) {
  int i;
  
  debug1(printf("query:  %s\n",query));
  debug1(printf("genome: %s\n",gbuffer));
  debug1(printf("count:  "));

  if (genestrand == +2) {
    for (i = start; i < end; i++) {
      if (gbuffer[i] == 'A' && query[i] == 'G') {
	debug1(printf("."));
#if 0
	/* Want to show mismatches */
	gbuffer[i] = 'G';		/* Avoids showing mismatches in MD and NM strings */
#endif
      } else if (query[i] != gbuffer[i]) {
	debug1(printf("x"));
	assert(gbuffer[i] != OUTOFBOUNDS);
	gbuffer[i] = (char) tolower(gbuffer[i]);
      } else {
	debug1(printf("*"));
      }
    }

  } else {
    for (i = start; i < end; i++) {
      if (gbuffer[i] == 'T' && query[i] == 'C') {
	debug1(printf("."));
#if 0
	/* Want to show mismatches */
	gbuffer[i] = 'C';		/* Avoids showing mismatches in MD and NM strings */
#endif
      } else if (query[i] != gbuffer[i]) {
	debug1(printf("x"));
	assert(gbuffer[i] != OUTOFBOUNDS);
	gbuffer[i] = (char) tolower(gbuffer[i]);
      } else {
	debug1(printf("*"));
      }
    }
  }

  return;
}



void
Substring_setup (bool print_nsnpdiffs_p_in, bool print_snplabels_p_in,
		 bool show_refdiff_p_in, IIT_T snps_iit_in, int *snps_divint_crosstable_in,
		 IIT_T genes_iit_in, int *genes_divint_crosstable_in,
		 IIT_T splicesites_iit_in, int *splicesites_divint_crosstable_in,
		 int donor_typeint_in, int acceptor_typeint_in, int trim_mismatch_score_in,
		 bool novelsplicingp_in, bool knownsplicingp_in,
		 bool output_sam_p_in, Mode_T mode_in, Univcoord_T genomelength_in) {
  print_nsnpdiffs_p = print_nsnpdiffs_p_in;
  print_snplabels_p = print_snplabels_p_in;
  show_refdiff_p = show_refdiff_p_in;

  snps_iit = snps_iit_in;
  snps_divint_crosstable = snps_divint_crosstable_in;

  genes_iit = genes_iit_in;
  genes_divint_crosstable = genes_divint_crosstable_in;

  splicesites_iit = splicesites_iit_in;
  splicesites_divint_crosstable = splicesites_divint_crosstable_in;

  donor_typeint = donor_typeint_in;
  acceptor_typeint = acceptor_typeint_in;

  trim_mismatch_score = trim_mismatch_score_in;
  
  novelsplicingp = novelsplicingp_in;
  knownsplicingp = knownsplicingp_in;

  output_sam_p = output_sam_p_in;
  mode = mode_in;

  genomelength = (double) genomelength_in;

  return;
}



static char *
embellish_genomic (char *genomic_diff, char *query, int querystart, int queryend, int querylength,
		   int extraleft, int extraright, int genestrand) {
  char *result;
  int i, j, k;

  debug1(printf("Entered embellish_genomic with querystart %d, queryend %d, querylength %d, genomic_diff %s\n",
		querystart,queryend,querylength,genomic_diff));

#ifdef DEBUG1
  result = (char *) CALLOC_OUT(querylength+1,sizeof(char));
#else
  result = (char *) MALLOC_OUT((querylength+1) * sizeof(char));
#endif
  result[querylength] = '\0';

  /* Add aligned region with lower-case diffs, surrounded by dashes */
  fill_w_dashes(result,0,querystart);
  debug1(printf("g1: %s (%d..%d) extraleft:%d extraright:%d\n",result,querystart,queryend,extraleft,extraright));

  /* Don't need to know adj anymore, because each substring has its own left */
  debug1(printf("Copying from genomic_diff[%d] to result[%d] for a length of %d - %d\n",querystart,querystart,queryend,querystart));
  strncpy(&(result[querystart]),&(genomic_diff[querystart]),queryend-querystart);
  debug1(printf("g1: %s (%d..%d) extraleft:%d extraright:%d\n",result,querystart,queryend,extraleft,extraright));

  if (mode == STANDARD) {
    /* Skip */
  } else if (mode == CMET_STRANDED || mode == CMET_NONSTRANDED) {
    mark_mismatches_cmet_gsnap(result,query,querystart,queryend,genestrand);
  } else if (mode == ATOI_STRANDED || mode == ATOI_NONSTRANDED) {
    mark_mismatches_atoi_gsnap(result,query,querystart,queryend,genestrand);
  } else if (mode == TTOC_STRANDED || mode == TTOC_NONSTRANDED) {
    mark_mismatches_ttoc_gsnap(result,query,querystart,queryend,genestrand);
  } else {
    abort();
  }

  fill_w_dashes(result,queryend,querylength);
  debug1(printf("g1: %s\n",result));

  /* Add terminal ends as lower-case */
  for (k = 0, i = querystart-1, j = querystart-1; k < extraleft && i >= 0 /*&& j >= 0*/; k++, i--, j--) {
    result[i] = (char) tolower(genomic_diff[j]);
    /* printf("k=%d i=%d result[i]=%c\n",k,i,result[i]); */
    assert(result[i] == 'a' || result[i] == 'c' || result[i] == 'g' || result[i] == 't' || result[i] == 'n');
  }

  for (k = 0, i = queryend, j = queryend; k < extraright && i < querylength /*&& j < genomiclength*/; k++, i++, j++) {
    result[i] = (char) tolower(genomic_diff[j]);
    /* printf("k=%d i=%d result[i]=%c\n",k,i,result[i]); */
#if 0
    assert(result[i] == 'a' || result[i] == 'c' || result[i] == 'g' || result[i] == 't' || result[i] == 'n');
#endif
  }
  debug1(printf("g1: %s\n",result));

  return result;
}


static char *
embellish_genomic_sam (char *genomic_diff, char *query, int querystart, int queryend, int querylength,
		       int genestrand, bool exactp) {
  char *result;
  int i, j, k;

  result = (char *) MALLOC_OUT((querylength+1) * sizeof(char));
  result[querylength] = '\0';

  strncpy(&(result[querystart]),&(genomic_diff[querystart]),queryend-querystart);

  if (mode == STANDARD) {
    /* Skip */
  } else if (mode == CMET_STRANDED || mode == CMET_NONSTRANDED) {
    mark_mismatches_cmet_sam(result,query,querystart,queryend,genestrand);
  } else if (mode == ATOI_STRANDED || mode == ATOI_NONSTRANDED) {
    mark_mismatches_atoi_sam(result,query,querystart,queryend,genestrand);
  } else if (mode == TTOC_STRANDED || mode == TTOC_NONSTRANDED) {
    mark_mismatches_ttoc_sam(result,query,querystart,queryend,genestrand);
  } else {
    abort();
  }

  /* Add terminal ends as lower-case */
  for (k = 0, i = querystart-1, j = querystart-1; i >= 0 && j >= 0; k++, i--, j--) {
    if (query[i] == genomic_diff[j]) {
      result[i] = genomic_diff[j];
    } else {
      result[i] = (char) tolower(genomic_diff[j]);
    }
    /* printf("k=%d i=%d j=%d result[i]=%c\n",k,i,j,result[i]); */
  }

  if (exactp == true) {
    /* No need to mark mismatches */
  } else {
    for (k = 0, i = queryend, j = queryend; i < querylength /*&& j < genomiclength*/; k++, i++, j++) {
      if (query[i] == genomic_diff[j]) {
	result[i] = genomic_diff[j];
	assert(result[i] == 'A' || result[i] == 'C' || result[i] == 'G' || result[i] == 'T' || result[i] == 'N');
      } else {
	result[i] = (char) tolower(genomic_diff[j]);
	assert(result[i] == 'a' || result[i] == 'c' || result[i] == 'g' || result[i] == 't' || result[i] == 'n');
      }
      /* printf("k=%d i=%d j=%d result[i]=%c\n",k,i,j,result[i]); */
    }
  }

  return result;
}


static int
substring_trim_novel_spliceends (T substring1, T substringN, int *ambig_end_length_5, int *ambig_end_length_3,
				 Splicetype_T *ambig_splicetype_5, Splicetype_T *ambig_splicetype_3,
				 double *ambig_prob_5, double *ambig_prob_3, int orig_sensedir,
				 Univcoord_T chroffset) {
  int new_sensedir;
  Univcoord_T start5, middle5, end5, start3, middle3, end3;
  Univcoord_T genomicstart5, genomicend3;
  bool solve5p, solve3p, plusp;


  debug13(printf("\nEntered substring_trim_novel_spliceends with orig_sensedir %d\n",orig_sensedir));
  *ambig_end_length_5 = 0;
  *ambig_end_length_3 = 0;
  *ambig_prob_5 = 0.0;
  *ambig_prob_3 = 0.0;

  /* start is distal, end is medial */
  if (substringN == NULL) {
    /* Skip 3' end */
    solve3p = false;
    genomicend3 = 0;

  } else if ((plusp = substringN->plusp) == true) {
    solve3p = true;
    genomicend3 = substringN->genomicend;

    middle3 = substringN->alignend_trim + 1;
    if ((start3 = middle3 + END_SPLICESITE_SEARCH) > substringN->genomicend) {
      start3 = substringN->genomicend;
    }
    if ((end3 = middle3 - END_SPLICESITE_SEARCH) < substringN->alignstart_trim + MIN_EXON_LENGTH) {
      end3 = substringN->alignstart_trim + MIN_EXON_LENGTH;
    }
    debug13(printf("\n1 Set end points for 3' trim to be %u..%u..%u\n",
		   start3 - chroffset,middle3 - chroffset,end3 - chroffset));

  } else {
    solve3p = true;
    genomicend3 = substringN->genomicend;

    middle3 = substringN->alignend_trim - 1;
    if ((start3 = middle3 - END_SPLICESITE_SEARCH) < substringN->genomicend) {
      start3 = substringN->genomicend;
    }
    if ((end3 = middle3 + END_SPLICESITE_SEARCH) > substringN->alignstart_trim - MIN_EXON_LENGTH) {
      end3 = substringN->alignstart_trim - MIN_EXON_LENGTH;
    }
    debug13(printf("\n2 Set end points for 3' trim to be %u..%u..%u\n",
		   start3 - chroffset,middle3 - chroffset,end3 - chroffset));
  }

  /* start is distal, end is medial */
  if (substring1 == NULL) {
    /* Skip 5' end */
    solve5p = false;
    genomicstart5 = 0;

  } else if ((plusp = substring1->plusp) == true) {
    solve5p = true;
    genomicstart5 = substring1->genomicstart;

    middle5 = substring1->alignstart_trim - 1;
    if ((start5 = middle5 - END_SPLICESITE_SEARCH) < substring1->genomicstart) {
      start5 = substring1->genomicstart;
    }
    if ((end5 = middle5 + END_SPLICESITE_SEARCH) > substring1->alignend_trim - MIN_EXON_LENGTH) {
      end5 = substring1->alignend_trim - MIN_EXON_LENGTH;
    }
    debug13(printf("\n1 Set end points for 5' trim to be %u..%u..%u\n",
		   start5 - chroffset,middle5 - chroffset,end5 - chroffset));

  } else {
    solve5p = true;
    genomicstart5 = substring1->genomicstart;

    middle5 = substring1->alignstart_trim + 1;
    if ((start5 = middle5 + END_SPLICESITE_SEARCH) > substring1->genomicstart) {
      start5 = substring1->genomicstart;
    }
    if ((end5 = middle5 - END_SPLICESITE_SEARCH) < substring1->alignend_trim + MIN_EXON_LENGTH) {
      end5 = substring1->alignend_trim + MIN_EXON_LENGTH;
    }
    debug13(printf("\n2 Set end points for 5' trim to be %u..%u..%u\n",
		   start5 - chroffset,middle5 - chroffset,end5 - chroffset));
  }

  new_sensedir = Splice_trim_novel_spliceends(&(*ambig_end_length_5),&(*ambig_end_length_3),
					      &(*ambig_splicetype_5),&(*ambig_splicetype_3),
					      &(*ambig_prob_5),&(*ambig_prob_3),orig_sensedir,
					      start5,middle5,end5,solve5p,start3,middle3,end3,solve3p,
					      genomicstart5,genomicend3,chroffset,plusp);

  debug13(printf("Returning ambig_end_length_5 %d and ambig_end_length_3 %d, probs %f and %f, new_sensedir %d\n",
		 *ambig_end_length_5,*ambig_end_length_3,*ambig_prob_5,*ambig_prob_3,new_sensedir));

  return new_sensedir;
}


/* Want querylength and not querylength_adj */
T
Substring_new (int nmismatches_whole, Chrnum_T chrnum, Univcoord_T chroffset,
	       Univcoord_T chrhigh, Chrpos_T chrlength,
	       Compress_T query_compress, Endtype_T start_endtype, Endtype_T end_endtype,
	       int querystart, int queryend, int querylength,
	       Univcoord_T alignstart, Univcoord_T alignend,
	       bool exactp, bool plusp, int genestrand,
	       Trimaction_T trim_left_action, Trimaction_T trim_right_action,
	       int outofbounds_start, int outofbounds_end, int minlength, int orig_sensedir) {
  T new;
  int nmatches;
  /* int nonterminal_trim = 0; */

  int ambig_end_length_5, ambig_end_length_3;
  Splicetype_T ambig_splicetype_5, ambig_splicetype_3;
  double ambig_prob_5, ambig_prob_3;
  int nmismatches_end_left, nmismatches_end_right;
  int new_sensedir = SENSE_NULL;
  int trim;

  
  /* General test for goodness over original region */
  nmatches = queryend - querystart - nmismatches_whole;
  debug2(printf("Initial nmatches %d = queryend %d - querystart %d - nmismatches_whole %d\n",
		nmatches,queryend,querystart,nmismatches_whole));

  if (nmatches - 3*nmismatches_whole < 0) {
    debug2(printf("Substring fails general test 1 for goodness with %d matches and %d mismatches\n",
		  nmatches,nmismatches_whole));
    return (T) NULL;
  } else {
    new = (T) MALLOC_OUT(sizeof(*new));
    debug2(printf("substring %p:\n",new));
  }

  new->exactp = exactp;

  new->chrnum = chrnum;
  new->chroffset = chroffset;
  new->chrhigh = chrhigh;
  new->chrlength = chrlength;

  /* new->genomicstart = genomicstart; */
  /* new->genomicend = genomicend; */
#if 0
  new->genomicstart_adj = genomicstart_adj;
  new->genomicend_adj = genomicend_adj;
#endif

  new->start_endtype = start_endtype;
  new->end_endtype = end_endtype;

  new->querystart_orig = new->querystart = querystart;
  new->queryend_orig = new->queryend = queryend;
  new->amb_splice_pos = 0;
  new->querylength = querylength;

  new->alignstart = new->alignstart_trim = alignstart;
  new->alignend = new->alignend_trim = alignend;

  new->plusp = plusp;
  new->genestrand = genestrand;

  new->splicecoord_D = new->splicecoord_A = new->splicecoord_N = 0;
  new->siteD_pos = new->siteA_pos = new->siteN_pos = 0;

  new->siteD_prob = new->siteA_prob = 0.0;
  new->siteD_knownp = new->siteA_knownp = false;
  new->siteD_novelp = new->siteA_novelp = false;

  debug2(printf("\n***Entered Substring_new with query %d..%d, chrnum %d (chroffset %u, chrhigh %u), plusp %d, outofbounds start %d and end %d\n",
		querystart,queryend,chrnum,chroffset,chrhigh,plusp,outofbounds_start,outofbounds_end));
  /* assert(queryend > querystart); -- Assertion does not hold.  Sometimes queryend == querystart */


  /* Compute coordinates */
  if (plusp == true) {
    new->left = alignstart - querystart;
    new->genomicstart = new->left;
    new->genomicend = new->left + querylength;

    debug2(printf("left is %u\n",new->left));
    debug2(printf("genomicstart is %u, genomicend is %u\n",new->genomicstart,new->genomicend));
    debug2(printf("querylength is %d, alignstart is %u, alignend is %u\n",querylength,alignstart,alignend));
    assert(alignstart + outofbounds_start >= chroffset);
    assert(alignend - outofbounds_end <= chrhigh);
    
  } else {
    new->left = alignend - (querylength - queryend);
    new->genomicend = new->left;
    new->genomicstart = new->left + querylength;

    debug2(printf("left is %u\n",new->left));
    debug2(printf("genomicstart is %u, genomicend is %u\n",new->genomicstart,new->genomicend));
    debug2(printf("querylength is %d, alignstart is %u, alignend is %u\n",querylength,alignstart,alignend));
    assert(alignstart - outofbounds_start <= chrhigh);
    assert(alignend + outofbounds_end >= chroffset);
  }
  
  if (nmismatches_whole < 0) {
    debug2(printf("Counting mismatches from querystart %d to queryend %d\n",querystart,queryend));
    if (plusp == true) {
      nmismatches_whole = 
	Genome_count_mismatches_substring(query_compress,new->left,/*pos5*/querystart,
					  /*pos3*/queryend,/*plusp*/true,genestrand);
    } else {
      nmismatches_whole = 
	Genome_count_mismatches_substring(query_compress,new->left,/*pos5*/querylength - queryend,
					  /*pos3*/querylength - querystart,/*plusp*/false,genestrand);
    }
  }

  /* Assign new->nmismatches_whole */
  new->nmismatches_whole = nmismatches_whole;

  /* Initialize these so an aborted Substring_free knows what to do */
  new->nambcoords = 0;
  new->genomic_bothdiff = (char *) NULL;
  new->genomic_refdiff = (char *) NULL;

  /* Do trimming */
  debug8(printf("trim_left_action %s, trim_right_action %s\n",
		Trimaction_string(trim_left_action),Trimaction_string(trim_right_action)));

  new->mandatory_trim_left = 0;
  new->mandatory_trim_right = 0;

  if (trim_left_action == PRE_TRIMMED) {
    new->trim_left = querystart;

  } else if (trim_left_action == NO_TRIM) {
    if (outofbounds_start > 0) {
      /* Not allowed to trim */
      debug2(printf("outofbounds_start %d > 0, so returning NULL\n",outofbounds_start));
      Substring_free(&new);
      return (T) NULL;
    } else {
      new->trim_left = querystart;
    }

  } else if (new->start_endtype == TERM) {
    /* Accept true terminals generated by GSNAP procedure */
    trim = trim_left_end(&nmismatches_end_left,query_compress,new->left,querystart,queryend,querylength,
			 plusp,genestrand,/*trim_mismatch_score*/-3);
    debug8(printf("trim_left_end: trim_left +%d from querystart %d, nmismatches_end_left = %d\n",
		  trim,querystart,nmismatches_end_left));
    if (outofbounds_start > querystart + trim) {
      trim = outofbounds_start - querystart;
    }
    new->querystart += trim;
    new->trim_left = new->querystart;
    debug8(printf("querystart is now %d\n",new->querystart));

  } else {
    trim = trim_left_end(&nmismatches_end_left,query_compress,new->left,querystart,queryend,querylength,
			 plusp,genestrand,trim_mismatch_score);
    debug8(printf("trim_left_end: trim_left +%d from querystart %d, nmismatches_end_left = %d\n",
		  trim,querystart,nmismatches_end_left));
    if (outofbounds_start > querystart + trim) {
      trim = outofbounds_start - querystart;
    }
    /* nonterminal_trim += new->trim_left; */
    new->querystart += trim;
    new->trim_left = new->querystart;
    debug8(printf("querystart is now %d\n",new->querystart));
  }

  if (trim_right_action == PRE_TRIMMED) {
    new->trim_right = querylength - queryend;

  } else if (trim_right_action == NO_TRIM) {
    if (outofbounds_end > 0) {
      /* Not allowed to trim */
      Substring_free(&new);
      debug2(printf("outofbounds_end %d > 0, so returning NULL\n",outofbounds_end));
      return (T) NULL;
    } else {
      new->trim_right = querylength - queryend;
    }

  } else if (new->end_endtype == TERM) {
    /* Accept true terminals generated by GSNAP procedure */
    trim = trim_right_end(&nmismatches_end_right,query_compress,new->left,querystart,queryend,querylength,
			  plusp,genestrand,/*trim_mismatch_score*/-3);
    debug8(printf("trim_right_end: trim_right +%d from queryend %d, nmismatches_end_right = %d\n",
		  trim,queryend,nmismatches_end_right));
    if (outofbounds_end > queryend - trim) {
      trim = queryend - outofbounds_end;
    }
    new->queryend -= trim;
    new->trim_right = querylength - new->queryend;
    debug8(printf("queryend is now %d\n",new->queryend));

  } else {
    trim = trim_right_end(&nmismatches_end_right,query_compress,new->left,querystart,queryend,querylength,
			  plusp,genestrand,trim_mismatch_score);
    debug8(printf("trim_right_end: trim_right +%d from queryend %d, nmismatches_end_right = %d\n",
		  trim,queryend,nmismatches_end_right));
    if (outofbounds_end > queryend - trim) {
      trim = queryend - outofbounds_end;
    }
    /* nonterminal_trim += new->trim_right; */
    new->queryend -= trim;
    new->trim_right = querylength - new->queryend;
    debug8(printf("queryend is now %d\n",new->queryend));
  }

#if 0
  if (nonterminal_trim >= reject_trimlength) {
    /* Reject non-terminal alignments (including those by sarray search) with excessive trim */
    /* Keep true terminals for now in case they help lead to GMAP alignments */
    debug8(printf("Nonterminal trims %d exceeds reject_trimlength %d, so returning NULL\n",
		  nonterminal_trim,reject_trimlength));
    Substring_free(&new);
    return (T) NULL;
  } else {
#endif
    /* Check for minlength.  Needed to avoid nonsensical terminal alignments */
    if (new->queryend - new->querystart <= minlength) {
      debug8(printf("queryend %d - querystart %d <= minlength %d, so returning NULL\n",
		    new->queryend,new->querystart,minlength));
      Substring_free(&new);
      return (T) NULL;
    }
#if 0
  }
#endif


  new->mandatory_trim_left = 0;
  new->mandatory_trim_right = 0;

  /* ? Should we spend the time to determine trim_left_splicep and
     trim_right_splicep, especially since trimming may not be perfect */
  if (plusp == true) {
    /* nmatches: Counts matches over whole region including trims */
    new->nmatches = (new->alignend - new->alignstart) - new->nmismatches_whole;

    if (trim_left_action == COMPUTE_TRIM) {
      if (querystart == 0 && nmismatches_end_left == 0) {
	debug8(printf("querystart is 0 and nmismatches_end_left is 0, so setting left_end_action to be NO_TRIM\n"));
	trim_left_action = NO_TRIM;
      } else {
	new->alignstart_trim = new->genomicstart + new->trim_left;
      }
    }
    if (trim_right_action == COMPUTE_TRIM) {
      if (queryend == querylength && nmismatches_end_right == 0) {
	debug8(printf("queryend is querylength and nmismatches_end_right is 0, so setting right_end_action to be NO_TRIM\n"));
	trim_right_action = NO_TRIM;
      } else {
	new->alignend_trim = new->genomicend - new->trim_right;
      }
    }
    debug2(printf("Got trims of %d and %d => Revised alignstart_trim and alignend_trim to be %u..%u (%u..%u)\n",
		  new->trim_left,new->trim_right,new->alignstart_trim,new->alignend_trim,
		  new->alignstart_trim - new->chroffset,new->alignend_trim - new->chroffset));
    debug2(printf("genomicstart is %u, genomicend is %u\n",new->genomicstart,new->genomicend));

    new->trim_left_splicep = new->trim_right_splicep = false;
    if (novelsplicingp == true) {
      if (trim_left_action == COMPUTE_TRIM && trim_right_action == COMPUTE_TRIM) {
	new_sensedir = substring_trim_novel_spliceends(/*substring1*/new,/*substringN*/new,&ambig_end_length_5,&ambig_end_length_3,
						       &ambig_splicetype_5,&ambig_splicetype_3,&ambig_prob_5,&ambig_prob_3,
						       orig_sensedir,chroffset);
	if (ambig_end_length_5 > 0) {
	  /* Revise trim_left to be at the splice site */
	  new->trim_left_splicep = true;
	  /* new->querystart += (ambig_end_length_5 - new->trim_left); */
	  /* new->alignstart_trim += (ambig_end_length_5 - new->trim_left); */
	  new->querystart = ambig_end_length_5;
	  new->alignstart_trim = new->genomicstart + ambig_end_length_5;

	  new->trim_left = ambig_end_length_5;
	  if (ambig_splicetype_5 == DONOR || ambig_splicetype_5 == ANTIDONOR) {
	    new->start_endtype = DON;
	    new->siteD_prob = ambig_prob_5;
	  } else {
	    new->start_endtype = ACC;
	    new->siteA_prob = ambig_prob_5;
	  }
	}
	if (ambig_end_length_3 > 0) {
	  /* Revise trim_right to be at the splice site */
	  new->trim_right_splicep = true;
	  /* new->queryend -= (ambig_end_length_3 - new->trim_right); */
	  /* new->alignend_trim -= (ambig_end_length_3 - new->trim_right); */
	  new->queryend = querylength - ambig_end_length_3;
	  new->alignend_trim = new->genomicend - ambig_end_length_3;

	  new->trim_right = ambig_end_length_3;
	  if (ambig_splicetype_3 == DONOR || ambig_splicetype_3 == ANTIDONOR) {
	    new->end_endtype = DON;
	    new->siteD_prob = ambig_prob_3;
	  } else {
	    new->end_endtype = ACC;
	    new->siteA_prob = ambig_prob_3;
	  }
	}

      } else if (trim_left_action == COMPUTE_TRIM) {
	new_sensedir = substring_trim_novel_spliceends(/*substring1*/new,/*substringN*/NULL,&ambig_end_length_5,&ambig_end_length_3,
						       &ambig_splicetype_5,&ambig_splicetype_3,&ambig_prob_5,&ambig_prob_3,
						       orig_sensedir,chroffset);
	if (ambig_end_length_5 > 0) {
	  /* Revise trim_left to be at the splice site */
	  new->trim_left_splicep = true;
	  /* new->querystart += (ambig_end_length_5 - new->trim_left); */
	  /* new->alignstart_trim += (ambig_end_length_5 - new->trim_left); */
	  new->querystart = ambig_end_length_5;
	  new->alignstart_trim = new->genomicstart + ambig_end_length_5;

	  new->trim_left = ambig_end_length_5;
	  if (ambig_splicetype_5 == DONOR || ambig_splicetype_5 == ANTIDONOR) {
	    new->start_endtype = DON;
	    new->siteD_prob = ambig_prob_5;
	  } else {
	    new->start_endtype = ACC;
	    new->siteA_prob = ambig_prob_5;
	  }
	}

      } else if (trim_right_action == COMPUTE_TRIM) {
	new_sensedir = substring_trim_novel_spliceends(/*substring1*/NULL,/*substringN*/new,&ambig_end_length_5,&ambig_end_length_3,
						       &ambig_splicetype_5,&ambig_splicetype_3,&ambig_prob_5,&ambig_prob_3,
						       orig_sensedir,chroffset);
	if (ambig_end_length_3 > 0) {
	  /* Revise trim_right to be at the splice site */
	  new->trim_right_splicep = true;
	  /* new->queryend -= (ambig_end_length_3 - new->trim_right); */
	  /* new->alignend_trim -= (ambig_end_length_3 - new->trim_right); */
	  new->queryend = querylength - ambig_end_length_3;
	  new->alignend_trim = new->genomicend - ambig_end_length_3;

	  new->trim_right = ambig_end_length_3;
	  if (ambig_splicetype_3 == DONOR || ambig_splicetype_3 == ANTIDONOR) {
	    new->end_endtype = DON;
	    new->siteD_prob = ambig_prob_3;
	  } else {
	    new->end_endtype = ACC;
	    new->siteA_prob = ambig_prob_3;
	  }
	}
      }
    }

  } else {
    /* nmatches: Counts matches over whole region including trims */
    new->nmatches = (new->alignstart - new->alignend) - new->nmismatches_whole;

    if (trim_left_action == COMPUTE_TRIM) {
      if (querystart == 0 && nmismatches_end_left == 0) {
	debug8(printf("querystart is 0 and nmismatches_end_left is 0, so setting left_end_action to be NO_TRIM\n"));
	trim_left_action = NO_TRIM;
      } else {
	new->alignstart_trim = new->genomicstart - new->trim_left;
      }
    }
    if (trim_right_action == COMPUTE_TRIM) {
      if (queryend == querylength && nmismatches_end_right == 0) {
	debug8(printf("queryend is querylength and nmismatches_end_right is 0, so setting right_end_action to be NO_TRIM\n"));
	trim_right_action = NO_TRIM;
      } else {
	new->alignend_trim = new->genomicend + new->trim_right;
      }
    }
    debug2(printf("Revised alignstart_trim and alignend_trim to be %u..%u (%u..%u)\n",
		  new->alignstart_trim,new->alignend_trim,
		  new->alignstart_trim - new->chroffset,new->alignend_trim - new->chroffset));

    new->trim_left_splicep = new->trim_right_splicep = false;
    if (novelsplicingp == true) {
      if (trim_left_action == COMPUTE_TRIM && trim_right_action == COMPUTE_TRIM) {
	new_sensedir = substring_trim_novel_spliceends(/*substring1*/new,/*substringN*/new,&ambig_end_length_5,&ambig_end_length_3,
						       &ambig_splicetype_5,&ambig_splicetype_3,&ambig_prob_5,&ambig_prob_3,
						       orig_sensedir,chroffset);
	if (ambig_end_length_5 > 0) {
	  /* Revise trim_left to be at the splice site */
	  new->trim_left_splicep = true;
	  /* new->querystart += (ambig_end_length_5 - new->trim_left); */
	  /* new->alignstart_trim -= (ambig_end_length_5 - new->trim_left); */
	  new->querystart = ambig_end_length_5;
	  new->alignstart_trim = new->genomicstart - ambig_end_length_5;

	  new->trim_left = ambig_end_length_5;
	  if (ambig_splicetype_5 == DONOR || ambig_splicetype_5 == ANTIDONOR) {
	    new->start_endtype = DON;
	    new->siteD_prob = ambig_prob_5;
	  } else {
	    new->start_endtype = ACC;
	    new->siteA_prob = ambig_prob_5;
	  }
	}
	if (ambig_end_length_3 > 0) {
	  /* Revise trim_right to be at the splice site */
	  new->trim_right_splicep = true;
	  /* new->queryend -= (ambig_end_length_3 - new->trim_right); */
	  /* new->alignend_trim += (ambig_end_length_3 - new->trim_right); */
	  new->queryend = querylength - ambig_end_length_3;
	  new->alignend_trim = new->genomicend + ambig_end_length_3;

	  new->trim_right = ambig_end_length_3;
	  if (ambig_splicetype_3 == DONOR || ambig_splicetype_3 == ANTIDONOR) {
	    new->end_endtype = DON;
	    new->siteD_prob = ambig_prob_3;
	  } else {
	    new->end_endtype = ACC;
	    new->siteA_prob = ambig_prob_3;
	  }
	}

      } else if (trim_left_action == COMPUTE_TRIM) {
	new_sensedir = substring_trim_novel_spliceends(/*substring1*/new,/*substringN*/NULL,&ambig_end_length_5,&ambig_end_length_3,
						       &ambig_splicetype_5,&ambig_splicetype_3,&ambig_prob_5,&ambig_prob_3,
						       orig_sensedir,chroffset);
	if (ambig_end_length_5 > 0) {
	  /* Revise trim_left to be at the splice site */
	  new->trim_left_splicep = true;
	  /* new->querystart += (ambig_end_length_5 - new->trim_left); */
	  /* new->alignstart_trim -= (ambig_end_length_5 - new->trim_left); */
	  new->querystart = ambig_end_length_5;
	  new->alignstart_trim = new->genomicstart - ambig_end_length_5;

	  new->trim_left = ambig_end_length_5;
	  if (ambig_splicetype_5 == DONOR || ambig_splicetype_5 == ANTIDONOR) {
	    new->start_endtype = DON;
	    new->siteD_prob = ambig_prob_5;
	  } else {
	    new->start_endtype = ACC;
	    new->siteA_prob = ambig_prob_5;
	  }
	}

      } else if (trim_right_action == COMPUTE_TRIM) {
	new_sensedir = substring_trim_novel_spliceends(/*substring1*/NULL,/*substringN*/new,&ambig_end_length_5,&ambig_end_length_3,
						       &ambig_splicetype_5,&ambig_splicetype_3,&ambig_prob_5,&ambig_prob_3,
						       orig_sensedir,chroffset);
	if (ambig_end_length_3 > 0) {
	  /* Revise trim_right to be at the splice site */
	  new->trim_right_splicep = true;
	  /* new->queryend -= (ambig_end_length_3 - new->trim_right); */
	  /* new->alignend_trim += (ambig_end_length_3 - new->trim_right); */
	  new->queryend = querylength - ambig_end_length_3;
	  new->alignend_trim = new->genomicend + ambig_end_length_3;

	  new->trim_right = ambig_end_length_3;
	  if (ambig_splicetype_3 == DONOR || ambig_splicetype_3 == ANTIDONOR) {
	    new->end_endtype = DON;
	    new->siteD_prob = ambig_prob_3;
	  } else {
	    new->end_endtype = ACC;
	    new->siteA_prob = ambig_prob_3;
	  }
	}
      } 
    }
  }

  new->sensedir = new_sensedir;

  /* nmismatches_bothdiff: Counts matches of trimmed region */
  if (new->trim_left == 0 && new->trim_right == 0) {
    new->nmismatches_bothdiff = nmismatches_whole;

  } else if (plusp == true) {
    new->nmismatches_bothdiff = 
      Genome_count_mismatches_substring(query_compress,new->left,/*pos5*/new->alignstart_trim-new->left,
					/*pos3*/new->alignend_trim - new->left,plusp,genestrand);
  } else {
    new->nmismatches_bothdiff = 
      Genome_count_mismatches_substring(query_compress,new->left,/*pos5*/new->alignend_trim - new->left,
					/*pos3*/new->alignstart_trim - new->left,plusp,genestrand);
  }

  /* General test for goodness over trimmed region */
  nmatches = new->queryend - new->querystart - new->nmismatches_bothdiff;
  debug2(printf("nmatches %d = queryend %d - querystart %d - nmismatches_bothdiff %d\n",
		nmatches,new->queryend,new->querystart,new->nmismatches_bothdiff));
  if (nmatches - 3*new->nmismatches_bothdiff < 0) {
    debug2(printf("Substring fails general test 2 for goodness with %d matches and %d mismatches\n",
		  nmatches,new->nmismatches_bothdiff));
    Substring_free(&new);
    return (T) NULL;
  }

  new->ambiguous_p = false;
  new->nambcoords = 0;
  new->ambcoords = (Univcoord_T *) NULL;
  new->amb_knowni = (int *) NULL;
  new->amb_nmismatches = (int *) NULL;
  new->amb_probs = (double *) NULL;
  new->amb_type = END;

  debug2(printf("Returning substring %p, query %d..%d, trim %d..%d\n",
		new,new->querystart,new->queryend,new->trim_left,new->trim_right));
  return new;
}


T
Substring_new_ambig_D (int querystart, int queryend, int splice_pos, int querylength,
		       Chrnum_T chrnum, Univcoord_T chroffset,
		       Univcoord_T chrhigh, Chrpos_T chrlength,
		       bool plusp, int genestrand,
#ifdef LARGE_GENOMES
		       Uint8list_T ambcoords,
#else
		       Uintlist_T ambcoords,
#endif
		       Intlist_T amb_knowni, Intlist_T amb_nmismatches, Doublelist_T amb_probs,
		       double amb_common_prob, bool substring1p) {
  int ignore;
  T new = (T) MALLOC_OUT(sizeof(*new));

  debug2(printf("Entered Substring_new_ambig_D with chrnum %d (chroffset %u, chrhigh %u), %d..%d, querylength %d, plusp %d\n",
		chrnum,chroffset,chrhigh,querystart,queryend,querylength,plusp));

  new->exactp = false;
  
  new->chrnum = chrnum;
  new->chroffset = chroffset;
  new->chrhigh = chrhigh;
  new->chrlength = chrlength;

  new->left = 0;
#ifdef LARGE_GENOMES
  if (plusp == true) {
    new->genomicstart = Uint8list_max(ambcoords);
    new->genomicend = Uint8list_min(ambcoords);
  } else {
    new->genomicstart = Uint8list_min(ambcoords);
    new->genomicend = Uint8list_max(ambcoords);
  }
#else
  if (plusp == true) {
    new->genomicstart = Uintlist_max(ambcoords);
    new->genomicend = Uintlist_min(ambcoords);
  } else {
    new->genomicstart = Uintlist_min(ambcoords);
    new->genomicend = Uintlist_max(ambcoords);
  }
#endif
  new->start_endtype = END;
  new->end_endtype = END;

  new->querystart_orig = new->querystart = querystart;
  new->queryend_orig = new->queryend = queryend;
  new->querylength = querylength;

  new->alignstart = new->alignstart_trim = 0;
  new->alignend = new->alignend_trim = 0;

  new->plusp = plusp;
  new->genestrand = genestrand;

  new->siteD_knownp = new->siteA_knownp = false;
  new->siteD_novelp = new->siteA_novelp = false;

  new->siteD_prob = 0.0;
  new->siteA_prob = amb_common_prob;

  new->nmismatches_bothdiff = new->nmismatches_whole = Intlist_min(amb_nmismatches);

#if 0
  if (plusp == true) {
    /* Fails because alignstart and alignend are not known */
    new->nmatches = (new->alignend_trim - new->alignstart_trim) - new->nmismatches_whole;
  } else {
    new->alignoffset = querylength - queryend;
    /* Fails because alignstart and alignend are not known */
    new->nmatches = (new->alignstart_trim - new->alignend_trim) - new->nmismatches_whole;
  }
#endif
  new->nmatches = (queryend - querystart) - new->nmismatches_whole;

  new->genomic_bothdiff = (char *) NULL;
  new->genomic_refdiff = (char *) NULL;
  if (substring1p == true) {
    debug2(printf("substring1p is true, so setting trims to be %d and %d\n",querystart,0));
    new->trim_left = querystart; /* trim_left includes the ambiguous part */
    new->trim_right = 0;
  } else {
    debug2(printf("substring1p is false, so setting trims to be %d and %d\n",0,querylength - queryend));
    new->trim_left = 0;
    new->trim_right = querylength - queryend; /* trim_right includes the ambiguous part */
  }
  new->mandatory_trim_left = 0;
  new->mandatory_trim_right = 0;
  new->trim_left_splicep = new->trim_right_splicep = false;


  new->ambiguous_p = true;
#ifdef LARGE_GENOMES
  new->ambcoords = Uint8list_to_array_out(&new->nambcoords,ambcoords);
#else
  new->ambcoords = Uintlist_to_array_out(&new->nambcoords,ambcoords);
  debug2(printf("ambcoords: %s\n",Uintlist_to_string(ambcoords)));
#endif
  new->amb_knowni = Intlist_to_array_out(&ignore,amb_knowni);
  new->amb_nmismatches = Intlist_to_array_out(&ignore,amb_nmismatches);
  new->amb_probs = Doublelist_to_array_out(&ignore,amb_probs);
  new->amb_splice_pos = splice_pos;
  new->amb_type = DON;

  return new;
}

T
Substring_new_ambig_A (int querystart, int queryend, int splice_pos, int querylength,
		       Chrnum_T chrnum, Univcoord_T chroffset,
		       Univcoord_T chrhigh, Chrpos_T chrlength,
		       bool plusp, int genestrand,
#ifdef LARGE_GENOMES
		       Uint8list_T ambcoords,
#else
		       Uintlist_T ambcoords,
#endif
		       Intlist_T amb_knowni, Intlist_T amb_nmismatches, Doublelist_T amb_probs,
		       double amb_common_prob, bool substring1p) {
  int ignore;
  T new = (T) MALLOC_OUT(sizeof(*new));

  debug2(printf("Entered Substring_new_ambig_A with chrnum %d (chroffset %u, chrhigh %u), %d..%d, querylength %d, plusp %d\n",
		chrnum,chroffset,chrhigh,querystart,queryend,querylength,plusp));

  new->exactp = false;
  
  new->chrnum = chrnum;
  new->chroffset = chroffset;
  new->chrhigh = chrhigh;
  new->chrlength = chrlength;

  new->left = 0;
#ifdef LARGE_GENOMES
  if (plusp == true) {
    new->genomicstart = Uint8list_max(ambcoords);
    new->genomicend = Uint8list_min(ambcoords);
  } else {
    new->genomicstart = Uint8list_min(ambcoords);
    new->genomicend = Uint8list_max(ambcoords);
  }
#else
  if (plusp == true) {
    new->genomicstart = Uintlist_max(ambcoords);
    new->genomicend = Uintlist_min(ambcoords);
  } else {
    new->genomicstart = Uintlist_min(ambcoords);
    new->genomicend = Uintlist_max(ambcoords);
  }
#endif
  new->start_endtype = END;
  new->end_endtype = END;

  new->querystart_orig = new->querystart = querystart;
  new->queryend_orig = new->queryend = queryend;
  new->querylength = querylength;

  new->alignstart = new->alignstart_trim = 0;
  new->alignend = new->alignend_trim = 0;

  new->plusp = plusp;
  new->genestrand = genestrand;

  new->siteD_knownp = new->siteA_knownp = false;
  new->siteD_novelp = new->siteA_novelp = false;

  new->siteA_prob = 0.0;
  new->siteD_prob = amb_common_prob;

  new->nmismatches_bothdiff = new->nmismatches_whole = Intlist_min(amb_nmismatches);

#if 0
  if (plusp == true) {
    /* Fails because alignstart and alignend are not known */
    new->nmatches = (new->alignend_trim - new->alignstart_trim) - new->nmismatches_whole;
  } else {
    new->alignoffset = querylength - queryend;
    /* Fails because alignstart and alignend are not known */
    new->nmatches = (new->alignstart_trim - new->alignend_trim) - new->nmismatches_whole;
  }
#endif
  new->nmatches = (queryend - querystart) - new->nmismatches_whole;

  new->genomic_bothdiff = (char *) NULL;
  new->genomic_refdiff = (char *) NULL;
  if (substring1p == true) {
    debug2(printf("substring1p is true, so setting trims to be %d and %d\n",querystart,0));
    new->trim_left = querystart; /* trim_left includes the ambiguous part */
    new->trim_right = 0;
  } else {
    debug2(printf("substring1p is false, so setting trims to be %d and %d\n",0,querylength - queryend));
    new->trim_left = 0;
    new->trim_right = querylength - queryend; /* trim_right includes the ambiguous part */
  }
  new->mandatory_trim_left = 0;
  new->mandatory_trim_right = 0;
  new->trim_left_splicep = new->trim_right_splicep = false;


  new->ambiguous_p = true;
#ifdef LARGE_GENOMES
  new->ambcoords = Uint8list_to_array_out(&new->nambcoords,ambcoords);
#else
  new->ambcoords = Uintlist_to_array_out(&new->nambcoords,ambcoords);
  debug2(printf("ambcoords: %s\n",Uintlist_to_string(ambcoords)));
#endif
  new->amb_knowni = Intlist_to_array_out(&ignore,amb_knowni);
  new->amb_nmismatches = Intlist_to_array_out(&ignore,amb_nmismatches);
  new->amb_probs = Doublelist_to_array_out(&ignore,amb_probs);
  new->amb_splice_pos = splice_pos;
  new->amb_type = ACC;

  return new;
}


#if 0
void
Substring_unambiguous_bounds (Univcoord_T *genomicstart, Univcoord_T *genomicend, T this, int bingoi) {
  Univcoord_T splicecoord, left;

  splicecoord = this->ambcoords[bingoi];
  
  if (this->plusp == true) {
    left = splicecoord - this->amb_splice_pos;
    *genomicstart = left;
    *genomicend = left + this->querylength;

  } else {
    left = splicecoord - (this->querylength - this->amb_splice_pos);
    *genomicend = left;
    *genomicstart = left + this->querylength;
  }

  return;
}
#endif


Univcoord_T
Substring_set_unambiguous (double *donor_prob, double *acceptor_prob, Univcoord_T *genomicstart, Univcoord_T *genomicend,
			   T this, int bingoi) {

#ifdef DEBUG2
  printf("Entered Substring_set_unambiguous for %d..%d.  plusp %d",this->querystart,this->queryend,this->plusp);
  if (this->amb_type == DON) {
    printf("type DON\n");
  } else {
    printf("type ACC\n");
  }
#endif

  this->nmismatches_whole = this->amb_nmismatches[bingoi];
  
  if (this->plusp == true) {
    if (this->amb_type == DON) {
      *acceptor_prob = this->siteA_prob;
      *donor_prob = this->siteD_prob = this->amb_probs[bingoi];
      this->splicecoord_D = this->ambcoords[bingoi];
      this->splicesitesD_knowni = this->amb_knowni[bingoi];
      this->left = this->splicecoord_D - this->amb_splice_pos;
    } else {
      *donor_prob = this->siteD_prob;
      *acceptor_prob = this->siteA_prob = this->amb_probs[bingoi];
      this->splicecoord_A = this->ambcoords[bingoi];
      this->splicesitesA_knowni = this->amb_knowni[bingoi];
      this->left = this->splicecoord_A - this->amb_splice_pos;
    }

    debug2(printf("left %u\n",this->left));
    *genomicstart = this->genomicstart = this->left;
    *genomicend = this->genomicend = this->left + this->querylength;
    this->alignstart = this->alignstart_trim = this->genomicstart + this->querystart;
    this->alignend = this->alignend_trim =  this->genomicstart + this->queryend;
    this->nmatches = (this->alignend - this->alignstart) - this->nmismatches_whole;

    debug2(printf("querypos %d..%d, alignstart is %u (%u), alignend is %u (%u), genomicstart is %u, genomicend is %u\n",
		  this->querystart,this->queryend,this->alignstart,this->alignstart - this->chroffset,
		  this->alignend,this->alignend - this->chroffset,this->genomicstart,this->genomicend));

  } else {
    if (this->amb_type == DON) {
      *acceptor_prob = this->siteA_prob;
      *donor_prob = this->siteD_prob = this->amb_probs[bingoi];
      this->splicecoord_D = this->ambcoords[bingoi];
      this->splicesitesD_knowni = this->amb_knowni[bingoi];
      this->left = this->splicecoord_D - (this->querylength - this->amb_splice_pos);
    } else {
      *donor_prob = this->siteD_prob;
      *acceptor_prob = this->siteA_prob = this->amb_probs[bingoi];
      this->splicecoord_A = this->ambcoords[bingoi];
      this->splicesitesA_knowni = this->amb_knowni[bingoi];
      this->left = this->splicecoord_A - (this->querylength - this->amb_splice_pos);
    }

    debug2(printf("left %u\n",this->left));
    *genomicend = this->genomicend = this->left;
    *genomicstart = this->genomicstart = this->left + this->querylength;
    this->alignend = this->alignend_trim = this->genomicstart - this->queryend;
    this->alignstart = this->alignstart_trim = this->genomicstart - this->querystart;
    this->nmatches = (this->alignstart - this->alignend) - this->nmismatches_whole;

    debug2(printf("querypos %d..%d, alignstart is %u (%u), alignend is %u (%u), genomicstart is %u, genomicend is %u\n",
		  this->querystart,this->queryend,this->alignstart,this->alignstart - this->chroffset,
		  this->alignend,this->alignend - this->chroffset,this->genomicstart,this->genomicend));
  }

  this->ambiguous_p = false;

  return this->left;
}


/* Look also at Pair_compute_mapq */
float
Substring_compute_mapq (T this, Compress_T query_compress, char *quality_string, bool trim_terminals_p) {
  int mapq_start, mapq_end;
  float best_loglik, loglik;
  Univcoord_T left, splicecoord;
  int i;

  /* mapq */
#if 0
  mapq_start = this->querystart_orig;
  mapq_end = this->queryend_orig;
#else
  /* Need trimmed regions to make this calculation equivalent to Pair_compute_mapq */
  mapq_start = this->querystart;
  mapq_end = this->queryend;
#endif

  /* It appears from simulated reads that it is better not to trim in
     computing MAPQ.  The correct mapping then tends to be selected
     with a higher MAPQ score. */
  /* But if all ends are terminals, then terminal parts should not be
     included in MAPQ scoring */

#if 0
  if (trim_terminals_p == true) {
    if (this->start_endtype == TERM) {
      mapq_start += this->trim_left;
    }
    if (this->end_endtype == TERM) {
      mapq_end -= this->trim_right;
    }
  }
#endif

  if (this->exactp == true) {
    /* this->mapq_loglik = MAPQ_loglik_exact(quality_string,0,querylength); */
    this->mapq_loglik = 0.0;

  } else if (this->ambiguous_p == true) {
    if (this->plusp == true) {
      splicecoord = this->ambcoords[0];
      left = splicecoord - this->amb_splice_pos;
      best_loglik = MAPQ_loglik(query_compress,left,mapq_start,mapq_end,
				this->querylength,quality_string,/*plusp*/true,this->genestrand);
      for (i = 1; i < this->nambcoords; i++) {
	splicecoord = this->ambcoords[i];
	left = splicecoord - this->amb_splice_pos;
	if ((loglik = MAPQ_loglik(query_compress,left,mapq_start,mapq_end,
				  this->querylength,quality_string,/*plusp*/true,this->genestrand)) > best_loglik) {
	  best_loglik = loglik;
	}
      }
    } else {
      splicecoord = this->ambcoords[0];
      left = splicecoord - (this->querylength - this->amb_splice_pos);
      best_loglik = MAPQ_loglik(query_compress,left,mapq_start,mapq_end,
				this->querylength,quality_string,/*plusp*/false,this->genestrand);
      for (i = 1; i < this->nambcoords; i++) {
	splicecoord = this->ambcoords[i];
	left = splicecoord - (this->querylength - this->amb_splice_pos);
	if ((loglik = MAPQ_loglik(query_compress,left,mapq_start,mapq_end,
				  this->querylength,quality_string,/*plusp*/false,this->genestrand)) > best_loglik) {
	  best_loglik = loglik;
	}
      }
    }

    this->mapq_loglik = best_loglik;

  } else {
    debug2(printf("trim_left %d, trim_right %d, mapq_start = %d, mapq_end = %d\n",
		  this->trim_left,this->trim_right,mapq_start,mapq_end));
    this->mapq_loglik = MAPQ_loglik(query_compress,this->left,mapq_start,mapq_end,
				    this->querylength,quality_string,this->plusp,this->genestrand);
    debug2(printf("Substring %u..%u gets loglik %f\n",this->genomicstart - this->chroffset,
		  this->genomicend - this->chroffset,this->mapq_loglik));
  }

  return this->mapq_loglik;
}


/* Note: query needed only for dibase */
int
Substring_display_prep (T this, char *queryuc_ptr, int querylength,
			int extraleft, int extraright,
			Compress_T query_compress_fwd, Compress_T query_compress_rev,
			Genome_T genome) {

#if defined(LONG_READLENGTHS)
  char *genomic_diff, *gbuffer;
#elif defined(HAVE_ALLOCA)
  char *genomic_diff, *gbuffer;
#else
  char *genomic_diff, *gbuffer;
  char gbuffer_alloc[MAX_READLENGTH/*+MAX_END_DELETIONS*/+1];
  bool allocp;
#endif

  /* genomic_bothdiff, genomic_refdiff, and nmismatches_refdiff */
  if (0 && this->exactp == true) {
    this->genomic_bothdiff = (char *) NULL;
    this->genomic_refdiff = (char *) NULL;
    this->nmismatches_refdiff = this->nmismatches_whole;

  } else if (this->plusp == true) {
    if (0 && this->exactp == true && extraleft == 0 && extraright == 0) {
      /* Don't use adj here */
      /* genomic_diff = &(queryuc_ptr[0]); */
      this->genomic_refdiff = (char *) NULL;
      this->nmismatches_refdiff = 0;

    } else {
      /* Used to be this->genomiclength, but doesn't work for large insertions */
#if defined(LONG_READLENGTHS)
      gbuffer = (char *) MALLOC((querylength+1) * sizeof(char));
#elif defined(HAVE_ALLOCA)
      gbuffer = (char *) ALLOCA((querylength+1) * sizeof(char));
#else
      if (querylength < MAX_READLENGTH) {
	gbuffer = gbuffer_alloc;
	allocp = false;
      } else {
	gbuffer = (char *) MALLOC((querylength+1) * sizeof(char));
	allocp = true;
      }
#endif

      debug1(printf("Obtaining genomic_diff from left %u (%u) for querylength %d\n",
		    this->left,this->left - this->chroffset,querylength));
      Genome_fill_buffer_simple(genome,this->left,querylength,gbuffer);
      genomic_diff = gbuffer;

      Genome_mark_mismatches(genomic_diff,querylength,query_compress_fwd,
			     this->left,/*pos5*/this->querystart,/*pos3*/this->queryend,
			     /*plusp*/true,this->genestrand);

    }
    /* Need to perform embellish to put dashes in */
    this->genomic_bothdiff = embellish_genomic(genomic_diff,queryuc_ptr,this->querystart,this->queryend,
					       querylength,extraleft,extraright,this->genestrand);

    if (snps_iit == NULL) {
      this->genomic_refdiff = this->genomic_bothdiff;
      this->nmismatches_refdiff = this->nmismatches_bothdiff;

    } else {
      this->nmismatches_refdiff =
	Genome_count_mismatches_substring_ref(query_compress_fwd,this->left,
					      /*pos5*/this->alignstart_trim - this->left,
					      /*pos3*/this->alignend_trim - this->left,
					      /*plusp*/true,this->genestrand);
	
      Genome_mark_mismatches_ref(genomic_diff,querylength,query_compress_fwd,this->left,
				 /*pos5*/this->querystart,/*pos3*/this->queryend,
				 /*plusp*/true,this->genestrand);
      if (output_sam_p == false) {
	this->genomic_refdiff = embellish_genomic(genomic_diff,queryuc_ptr,this->querystart,this->queryend,
						  querylength,extraleft,extraright,this->genestrand);
      }
    }

    if (output_sam_p == true) {
      this->genomic_refdiff = embellish_genomic_sam(genomic_diff,queryuc_ptr,this->querystart,this->queryend,
						    querylength,this->genestrand,this->exactp);
    }

    if (0 && this->exactp == true && extraleft == 0 && extraright == 0) {
    } else {
#if defined(LONG_READLENGTHS)
      FREE(gbuffer);
#elif defined(HAVE_ALLOCA)
      FREEA(gbuffer);
#else
      if (allocp == true) {
	FREE(gbuffer);
      }
#endif
    }

  } else {
    if (0 && this->exactp == true && extraleft == 0 && extraright == 0) {
      /* Don't use adj here */
      /* genomic_diff = &(queryuc_ptr[0]); -- Not queryrc */
      this->genomic_refdiff = (char *) NULL;
      this->nmismatches_refdiff = 0;

    } else {
      /* Used to be this->genomiclength, but doesn't work for large insertions */
#if defined(LONG_READLENGTHS)
      gbuffer = (char *) MALLOC((querylength+1) * sizeof(char));
#elif defined(HAVE_ALLOCA)
      gbuffer = (char *) ALLOCA((querylength+1) * sizeof(char));
#else
      if (querylength < MAX_READLENGTH) {
	gbuffer = gbuffer_alloc;
	allocp = false;
      } else {
	gbuffer = (char *) MALLOC((querylength+1) * sizeof(char));
	allocp = true;
      }
#endif

      debug1(printf("Obtaining genomic_diff from left %u (%u) for querylength %d, and complemented\n",
		    this->left,this->left - this->chroffset,querylength));
      Genome_fill_buffer_simple(genome,this->left,querylength,gbuffer);
      genomic_diff = make_complement_inplace(gbuffer,querylength);

      Genome_mark_mismatches(genomic_diff,querylength,query_compress_rev,
			     this->left,/*pos5*/querylength - this->queryend,
			     /*pos3*/querylength - this->querystart,
			     /*plusp*/false,this->genestrand);
    }

    /* Need to perform embellish to put dashes in */
    this->genomic_bothdiff = embellish_genomic(genomic_diff,/*not queryrc*/queryuc_ptr,this->querystart,this->queryend,
					       querylength,extraleft,extraright,this->genestrand);

    if (snps_iit == NULL) {
      this->genomic_refdiff = this->genomic_bothdiff;
      this->nmismatches_refdiff = this->nmismatches_bothdiff;
      
    } else {
      this->nmismatches_refdiff = 
	Genome_count_mismatches_substring_ref(query_compress_rev,this->left,
					      /*pos5*/this->alignend_trim - this->left,
					      /*pos3*/this->alignstart_trim - this->left,/*plusp*/false,
					      this->genestrand);
      
      Genome_mark_mismatches_ref(genomic_diff,querylength,query_compress_rev,this->left,
				 /*pos5*/querylength - this->queryend,
				 /*pos3*/querylength - this->querystart,
				 /*plusp*/false,this->genestrand);
      
      if (output_sam_p == false) {
	this->genomic_refdiff = embellish_genomic(genomic_diff,/*not queryrc*/queryuc_ptr,this->querystart,this->queryend,
						  querylength,extraleft,extraright,this->genestrand);
      }
    }

    if (output_sam_p == true) {
      this->genomic_refdiff = embellish_genomic_sam(genomic_diff,/*not queryrc*/queryuc_ptr,this->querystart,this->queryend,
						    querylength,this->genestrand,this->exactp);
    }

    if (0 && this->exactp == true && extraleft == 0 && extraright == 0) {
    } else {
#if defined(LONG_READLENGTHS)
      FREE(gbuffer);
#elif defined(HAVE_ALLOCA)
      FREEA(gbuffer);
#else
      if (allocp == true) {
	FREE(gbuffer);
      }
#endif
    }
  }

  return this->nmismatches_refdiff;
}


Univcoord_T
Substring_left (T this) {
  return this->left;
}


Univcoord_T
Substring_splicecoord_D (T this) {
  return this->splicecoord_D;
}

Univcoord_T
Substring_splicecoord_A (T this) {
  return this->splicecoord_A;
}

Chrpos_T
Substring_chr_splicecoord_D (T this) {
  return (Chrpos_T) (this->splicecoord_D - this->chroffset);
}

Chrpos_T
Substring_chr_splicecoord_A (T this) {
  return (Chrpos_T) (this->splicecoord_A - this->chroffset);
}

int
Substring_splicesitesD_knowni (T this) {
  return this->splicesitesD_knowni;
}

int
Substring_splicesitesA_knowni (T this) {
  return this->splicesitesA_knowni;
}

bool
Substring_plusp (T this) {
  return this->plusp;
}

int
Substring_genestrand (T this) {
  return this->genestrand;
}

char *
Substring_genomic_bothdiff (T this) {
  return this->genomic_bothdiff;
}

char *
Substring_genomic_refdiff (T this) {
  return this->genomic_refdiff;
}


int
Substring_nmismatches_whole (T this) {
  return this->nmismatches_whole;
}


int
Substring_nmismatches_bothdiff (T this) {
  return this->nmismatches_bothdiff;
}


int
Substring_nmismatches_refdiff (T this) {
  return this->nmismatches_refdiff;
}

int
Substring_nmatches_posttrim (T this) {
  if (this->ambiguous_p == true) {
    return 0;
  } else {
    return this->queryend - this->querystart - this->nmismatches_bothdiff;
  }
}

/* nmatches_posttrim plus amb_length */
int
Substring_nmatches (T this) {
  int amb_length;

  if (this->ambiguous_p == false) {
    /* Add the amb_length part */
    amb_length = 0;
    if (this->trim_left_splicep == true) {
      amb_length += this->trim_left;
    }
    if (this->trim_right_splicep == true) {
      amb_length += this->trim_right;
    }
    return amb_length + /*nmatches_posttrim*/ (this->queryend - this->querystart - this->nmismatches_bothdiff);

  } else {
    /* Include the entire ambiguous substring */
    return (this->queryend - this->querystart);
  }
}


void
Substring_set_nmismatches_terminal (T this, int nmismatches_whole, int nmismatches_bothdiff) {
  this->nmismatches_whole = nmismatches_whole;
  this->nmismatches_bothdiff = nmismatches_bothdiff;
  if (this->plusp == true) {
    this->nmatches = (this->alignend_trim - this->alignstart_trim) - nmismatches_whole;
  } else {
    this->nmatches = (this->alignstart_trim - this->alignend_trim) - nmismatches_whole;
  }
  return;
}

Endtype_T
Substring_start_endtype (T this) {
  return this->start_endtype;
}

Endtype_T
Substring_end_endtype (T this) {
  return this->end_endtype;
}

#if 0
void
Substring_set_start_endtype (T this, Endtype_T start_endtype, int start_amb_length) {
  this->start_endtype = start_endtype;
  this->start_amb_length = start_amb_length;
  return;
}

void
Substring_set_end_endtype (T this, Endtype_T end_endtype, int end_amb_length) {
  this->end_endtype = end_endtype;
  this->end_amb_length = end_amb_length;
  return;
}
#endif


float
Substring_mapq_loglik (T this) {
  return this->mapq_loglik;
}

int
Substring_trim_left (T this) {
  return this->trim_left;
}

int
Substring_trim_right (T this) {
  return this->trim_right;
}

bool
Substring_trim_left_splicep (T this) {
  return this->trim_left_splicep;
}

bool
Substring_trim_right_splicep (T this) {
  return this->trim_right_splicep;
}

int
Substring_mandatory_trim_left (T this) {
  return this->mandatory_trim_left;
}

int
Substring_mandatory_trim_right (T this) {
  return this->mandatory_trim_right;
}


int
Substring_querystart (T this) {
  return this->querystart;
}

int
Substring_querystart_orig (T this) {
  return this->querystart_orig;
}

int
Substring_queryend (T this) {
  return this->queryend;
}

int
Substring_queryend_orig (T this) {
  return this->queryend_orig;
}

int
Substring_querylength (T this) {
  return this->querylength;
}

int
Substring_match_length (T this) {
  if (this == NULL) {
    return 0;
  } else {
    return this->queryend - this->querystart;
  }
}

/* Mapped and unmapped */
int
Substring_amb_length (T this) {
  int amb_length;

  if (this->ambiguous_p == false) {
    /* For GSNAP algorithm, we represent trimmed ambiguous ends using these fields */
    amb_length = 0;
    if (this->trim_left_splicep == true) {
      amb_length += this->trim_left;
    }
    if (this->trim_right_splicep == true) {
      amb_length += this->trim_right;
    }
    return amb_length;

  } else {
    /* For substrings algorithm, we represent known ambiguous ends with a substring having ambiguous_p == true */
    return this->queryend - this->querystart;
  }
}


/* Mapped and unmapped */
int
Substring_start_amb_length (T this) {
  if (this->ambiguous_p == false) {
    /* For GSNAP algorithm, we represent trimmed ambiguous ends using these fields */
    if (this->trim_left_splicep == true) {
      return this->trim_left;
    } else {
      return 0;
    }

  } else {
    /* For substrings algorithm, we represent known ambiguous ends with a substring having ambiguous_p == true */
    return this->queryend - this->querystart;
  }
}


/* Mapped and unmapped */
int
Substring_end_amb_length (T this) {
  if (this->ambiguous_p == false) {
    /* For GSNAP algorithm, we represent trimmed ambiguous ends using these fields */
    if (this->trim_right_splicep == true) {
      return this->trim_right;
    } else {
      return 0;
    }

  } else {
    /* For substrings algorithm, we represent known ambiguous ends with a substring having ambiguous_p == true */
    return this->queryend - this->querystart;
  }
}


/* Before trimming */
int
Substring_match_length_orig (T this) {
  if (this == NULL) {
    return 0;
  } else {
    return this->queryend_orig - this->querystart_orig;
  }
}

Chrpos_T
Substring_genomic_alignment_length (T this) {
  if (this == NULL) {
    return 0U;
  } else if (this->alignend_trim > this->alignstart_trim) {
    return this->alignend_trim - this->alignstart_trim;
  } else {
    return this->alignstart_trim - this->alignend_trim;
  }
}


Chrnum_T
Substring_chrnum (T this) {
  return this->chrnum;
}

Univcoord_T
Substring_chroffset (T this) {
  return this->chroffset;
}

Univcoord_T
Substring_chrhigh (T this) {
  return this->chrhigh;
}

Chrpos_T
Substring_chrlength (T this) {
  return this->chrlength;
}

Univcoord_T
Substring_alignstart (T this) {
  return this->alignstart;
}

Univcoord_T
Substring_alignend (T this) {
  return this->alignend;
}

Chrpos_T
Substring_alignstart_chr (T this) {
  if (this->ambiguous_p == true) {
    return 0;
  } else {
    return this->alignstart - this->chroffset;
  }
}

Chrpos_T
Substring_alignend_chr (T this) {
  if (this->ambiguous_p == true) {
    return 0;
  } else {
    return this->alignend - this->chroffset;
  }
}

Chrpos_T
Substring_alignstart_trim_chr (T this) {
  if (this->ambiguous_p == true) {
    return 0;
  } else {
    return this->alignstart_trim - this->chroffset;
  }
}

Chrpos_T
Substring_alignend_trim_chr (T this) {
  if (this->ambiguous_p == true) {
    return 0;
  } else {
    return this->alignend_trim - this->chroffset;
  }
}

Univcoord_T
Substring_alignstart_trim (T this) {
  return this->alignstart_trim;
}

Univcoord_T
Substring_alignend_trim (T this) {
  return this->alignend_trim;
}


#if 0
Univcoord_T
Substring_alignmid_trim (T this) {
  if (this->alignend_trim > this->alignstart_trim) {
    return this->alignstart_trim + (this->alignend_trim - this->alignstart_trim)/2;
  } else {
    return this->alignend_trim + (this->alignstart_trim - this->alignend_trim)/2;
  }
}
#endif


Univcoord_T
Substring_left_genomicseg (T this) {
  return this->left;
}

Univcoord_T
Substring_genomicstart (T this) {
  return this->genomicstart;
}

Chrpos_T
Substring_genomicstart_chr (T this) {
  return this->genomicstart - this->chroffset;
}

Univcoord_T
Substring_genomicend (T this) {
  return this->genomicend;
}

Chrpos_T
Substring_genomicend_chr (T this) {
  return this->genomicend - this->chroffset;
}

double
Substring_amb_donor_prob (T this) {
  double max;
  int i;

  if (this->amb_type == DON) {
    return this->siteD_prob;
  } else {
    max = this->amb_probs[0];
    for (i = 1; i < this->nambcoords; i++) {
      if (this->amb_probs[i] > max) {
	max = this->amb_probs[i];
      }
    }
    return max;
  }
}

double
Substring_amb_acceptor_prob (T this) {
  double max;
  int i;

  if (this->amb_type == ACC) {
    return this->siteA_prob;
  } else {
    max = this->amb_probs[0];
    for (i = 1; i < this->nambcoords; i++) {
      if (this->amb_probs[i] > max) {
	max = this->amb_probs[i];
      }
    }
    return max;
  }
}



double
Substring_siteD_prob (T this) {
  return this->siteD_prob;
}  

double
Substring_siteA_prob (T this) {
  return this->siteA_prob;
}  


int
Substring_siteD_pos (T this) {
  return this->siteD_pos;
}

int
Substring_siteA_pos (T this) {
  return this->siteA_pos;
}

int
Substring_siteN_pos (T this) {
  return this->siteN_pos;
}


int
Substring_nchimera_known (T this) {
  if (this == NULL) {
    return 0;
  } else {
    return (int) this->siteD_knownp + (int) this->siteA_knownp;
  }
}

int
Substring_nchimera_novel (T this) {
  if (this == NULL) {
    return 0;
  } else {
    return (int) this->siteD_novelp + (int) this->siteA_novelp;
  }
}


int
Substring_sensedir (T this) {
  return this->sensedir;
}


bool
Substring_ambiguous_p (T this) {
  assert(this->ambiguous_p == false || this->ambiguous_p == true);
  return this->ambiguous_p;
}

bool
Substring_list_ambiguous_p (List_T list) {
  T this;
  List_T p;

  for (p = list; p != NULL; p = List_next(p)) {
    this = (T) List_head(p);
    if (this->ambiguous_p == true) {
      return true;
    }
  }
  return false;
}

int
Substring_nambcoords (T this) {
  return this->nambcoords;
}

Univcoord_T *
Substring_ambcoords (T this) {
  return this->ambcoords;
}

int *
Substring_amb_knowni (T this) {
  return this->amb_knowni;
}

int *
Substring_amb_nmismatches (T this) {
  return this->amb_nmismatches;
}

double *
Substring_amb_probs (T this) {
  return this->amb_probs;
}




/* circularpos measures query distance from SAM chrlow to origin */
int
Substring_circularpos (T this) {
  if (this == NULL) {
    return -1;

  } else if (this->plusp == true) {
    /* printf("substring plus: looking at %u (trim %u) vs %u+%u\n",
       this->alignend,this->alignend_trim,this->chroffset,this->chrlength); */
    if (this->alignend_trim > this->chroffset + this->chrlength) {
      /* return (this->querystart - this->trim_left) + (this->chroffset + this->chrlength) - this->alignstart; */
      return this->querystart + (this->chroffset + this->chrlength) - this->alignstart_trim;
    } else {
      return -1;
    }

  } else {
    /* printf("substring minus: looking at %u (trim %u) vs %u+%u\n",
       this->alignstart,this->alignstart_trim,this->chroffset,this->chrlength); */
    if (this->alignstart_trim > this->chroffset + this->chrlength) {
      /* return ((this->querylength - this->trim_right) - this->queryend) + (this->chroffset + this->chrlength) - this->alignend; */
      return (this->querylength - this->queryend) + (this->chroffset + this->chrlength) - this->alignend_trim;
    } else {
      return -1;
    }
  }
}


T
Substring_copy (T old) {
  T new;

  if (old == NULL) {
    return (T) NULL;
  } else {
    new = (T) MALLOC_OUT(sizeof(*new));
    debug2(printf("substring %p is a copy of %p\n",new,old));

    new->exactp = old->exactp;
    new->nmismatches_whole = old->nmismatches_whole;
    new->nmismatches_bothdiff = old->nmismatches_bothdiff;
    new->nmismatches_refdiff = old->nmismatches_refdiff;
    new->nmatches = old->nmatches;

    new->trim_left = old->trim_left;
    new->trim_right = old->trim_right;
    new->mandatory_trim_left = old->mandatory_trim_left;
    new->mandatory_trim_right = old->mandatory_trim_right;

    new->trim_left_splicep = old->trim_left_splicep;
    new->trim_right_splicep = old->trim_right_splicep;

    new->chrnum = old->chrnum;
    new->chroffset = old->chroffset;
    new->chrhigh = old->chrhigh;
    new->chrlength = old->chrlength;

    new->left = old->left;
    new->genomicstart = old->genomicstart;
    new->genomicend = old->genomicend;

    new->start_endtype = old->start_endtype;
    new->end_endtype = old->end_endtype;

    new->querystart_orig = old->querystart_orig;
    new->queryend_orig = old->queryend_orig;
    new->querystart = old->querystart;
    new->queryend = old->queryend;
    new->amb_splice_pos = old->amb_splice_pos;
    new->querylength = old->querylength;

    new->alignstart = old->alignstart;
    new->alignend = old->alignend;

    new->alignstart_trim = old->alignstart_trim;
    new->alignend_trim = old->alignend_trim;

    new->plusp = old->plusp;
    new->genestrand = old->genestrand;

    if (old->genomic_bothdiff == NULL) {
      new->genomic_bothdiff = (char *) NULL;
      new->genomic_refdiff = (char *) NULL;
    } else if (old->genomic_refdiff == old->genomic_bothdiff) {
      new->genomic_bothdiff = (char *) CALLOC_OUT(strlen(old->genomic_bothdiff)+1,sizeof(char));
      strcpy(new->genomic_bothdiff,old->genomic_bothdiff);
      new->genomic_refdiff = new->genomic_bothdiff;
    } else {
      new->genomic_bothdiff = (char *) CALLOC_OUT(strlen(old->genomic_bothdiff)+1,sizeof(char));
      strcpy(new->genomic_bothdiff,old->genomic_bothdiff);
      new->genomic_refdiff = (char *) CALLOC_OUT(strlen(old->genomic_refdiff)+1,sizeof(char));
      strcpy(new->genomic_refdiff,old->genomic_refdiff);
    }

    new->mapq_loglik = old->mapq_loglik;

    new->sensedir = old->sensedir;

    new->splicecoord_D = old->splicecoord_D;
    new->splicesitesD_knowni = old->splicesitesD_knowni;
    new->siteD_knownp = old->siteD_knownp;
    new->siteD_novelp = old->siteD_novelp;
    new->siteD_pos = old->siteD_pos;
    new->siteD_prob = old->siteD_prob;

    new->splicecoord_A = old->splicecoord_A;
    new->splicesitesA_knowni = old->splicesitesA_knowni;
    new->siteA_knownp = old->siteA_knownp;
    new->siteA_novelp = old->siteA_novelp;
    new->siteA_pos = old->siteA_pos;
    new->siteA_prob = old->siteA_prob;

    new->splicecoord_N = old->splicecoord_N;
    new->siteN_pos = old->siteN_pos;

    new->ambiguous_p = old->ambiguous_p;
    if (old->nambcoords == 0) {
      new->nambcoords = 0;
      new->ambcoords = (Univcoord_T *) NULL;
      new->amb_knowni = (int *) NULL;
      new->amb_nmismatches = (int *) NULL;
      new->amb_probs = (double *) NULL;
    } else {
      new->nambcoords = old->nambcoords;
      new->ambcoords = (Univcoord_T *) MALLOC_OUT(old->nambcoords * sizeof(Univcoord_T));
      new->amb_knowni = (int *) MALLOC_OUT(old->nambcoords * sizeof(int));
      new->amb_nmismatches = (int *) MALLOC_OUT(old->nambcoords * sizeof(int));
      new->amb_probs = (double *) MALLOC_OUT(old->nambcoords * sizeof(double));

      memcpy(new->ambcoords,old->ambcoords,old->nambcoords * sizeof(Univcoord_T));
      memcpy(new->amb_knowni,old->amb_knowni,old->nambcoords * sizeof(int));
      memcpy(new->amb_nmismatches,old->amb_nmismatches,old->nambcoords * sizeof(int));
      memcpy(new->amb_probs,old->amb_probs,old->nambcoords * sizeof(double));
    }
    new->amb_type = old->amb_type;

    return new;
  }
}


/* Treat as a sense donor, on either plus or minus strand */
T
Substring_new_startfrag (Univcoord_T startfrag_coord, int splice_pos, int nmismatches,
			 Univcoord_T left, Compress_T query_compress,
			 int querylength, bool plusp, int genestrand,
			 Chrnum_T chrnum, Univcoord_T chroffset, Univcoord_T chrhigh, Chrpos_T chrlength) {
  T new;
  int querystart, queryend;
  Univcoord_T genomicstart, alignstart, alignend;
  Endtype_T start_endtype, end_endtype;
  Trimaction_T trim_left_action, trim_right_action;

  /* Previously checked if left >= chroffset + chrlength to exclude
     the duplicate length, but now excluding all translocations to
     circular chromosomes */

  if (chroffset + chrlength < chrhigh) {
    /* Don't splice to circular chromosomes */
    return (T) NULL;

  } else if (plusp == true) {
    genomicstart = left;
    /* genomicend = left + querylength; */

    start_endtype = END;
    end_endtype = FRAG;

    querystart = 0;
    queryend = splice_pos;
    alignstart = genomicstart;
    alignend = genomicstart + splice_pos;
    trim_left_action = COMPUTE_TRIM;	/* querystart == 0 */
    trim_right_action = NO_TRIM;

  } else {
    genomicstart = left + querylength;
    /* genomicend = left; */

    start_endtype = END;
    end_endtype = FRAG;

    querystart = 0;
    queryend = querylength - splice_pos;
    alignstart = genomicstart;
    alignend = genomicstart - (querylength - splice_pos);
    trim_left_action = COMPUTE_TRIM;	/* querystart == 0 */
    trim_right_action = NO_TRIM;
  }

  if ((new = Substring_new(nmismatches,chrnum,chroffset,chrhigh,chrlength,
			   query_compress,start_endtype,end_endtype,querystart,queryend,querylength,
			   alignstart,alignend,/*exactp*/false,plusp,genestrand,
			   trim_left_action,trim_right_action,/*outofbounds_start*/0,/*outofbounds_end*/0,
			   /*minlength*/0,/*orig_sensedir*/SENSE_NULL)) == NULL) {
    return (T) NULL;
  }

  debug2(printf("Making new startfrag with coord %u and left %u, plusp %d, query %d..%d, genome %u..%u\n",
		startfrag_coord,left,plusp,querystart,queryend,alignstart - chroffset,alignend - chroffset));
  new->splicecoord_N = startfrag_coord;
  assert(startfrag_coord == left + splice_pos);

  if (plusp == true) {
    new->siteN_pos = splice_pos;
  } else {
    new->siteN_pos = querylength - splice_pos;
  }

  return new;
}


/* Treat as a sense acceptor, on either plus or minus strand */
T
Substring_new_endfrag (Univcoord_T endfrag_coord, int splice_pos, int nmismatches,
		       Univcoord_T left, Compress_T query_compress,
		       int querylength, bool plusp, int genestrand,
		       Chrnum_T chrnum, Univcoord_T chroffset, Univcoord_T chrhigh, Chrpos_T chrlength) {
  T new;
  int querystart, queryend;
  Univcoord_T genomicstart, alignstart, alignend;
  Endtype_T start_endtype, end_endtype;
  Trimaction_T trim_left_action, trim_right_action;

  /* Previously checked if left >= chroffset + chrlength to exclude
     the duplicate length, but now excluding all translocations to
     circular chromosomes */

  if (chroffset + chrlength < chrhigh) {
    /* Don't splice to circular chromosomes */
    return (T) NULL;

  } else if (plusp == true) {
    genomicstart = left;
    /* genomicend = left + querylength; */

    start_endtype = FRAG;
    end_endtype = END;

    querystart = splice_pos;
    queryend = querylength;
    alignstart = genomicstart + splice_pos;
    alignend = genomicstart + querylength;
    trim_left_action = NO_TRIM;
    trim_right_action = COMPUTE_TRIM; /* queryend == querylength */

  } else {
    genomicstart = left + querylength;
    /* genomicend = left; */

    start_endtype = FRAG;
    end_endtype = END;

    querystart = querylength - splice_pos;
    queryend = querylength;
    alignstart = left + splice_pos;
    alignend = left;
    trim_left_action = NO_TRIM;
    trim_right_action = COMPUTE_TRIM;	/* queryend == querylength */
  }

  if ((new = Substring_new(nmismatches,chrnum,chroffset,chrhigh,chrlength,
			   query_compress,start_endtype,end_endtype,querystart,queryend,querylength,
			   alignstart,alignend,/*exactp*/false,plusp,genestrand,
			   trim_left_action,trim_right_action,/*outofbounds_start*/0,/*outofbounds_end*/0,
			   /*minlength*/0,/*sensedir*/SENSE_NULL)) == NULL) {
    return (T) NULL;
  }

  debug2(printf("Making new endfrag with coord %u and left %u, plusp %d, query %d..%d, genome %u..%u\n",
		endfrag_coord,left,plusp,querystart,queryend,alignstart - chroffset,alignend - chroffset));
  new->splicecoord_N = endfrag_coord;
  assert(endfrag_coord == left + splice_pos);

  if (plusp == true) {
    new->siteN_pos = splice_pos;
  } else {
    new->siteN_pos = querylength - splice_pos;
  }

  return new;
}



T
Substring_new_donor (Univcoord_T donor_coord, int donor_knowni, int donor_pos,
		     int substring_querystart, int substring_queryend, int donor_nmismatches,
		     double donor_prob, Univcoord_T left, Compress_T query_compress,
		     int querylength, bool plusp, int genestrand, int sensedir,
		     Chrnum_T chrnum, Univcoord_T chroffset, Univcoord_T chrhigh, Chrpos_T chrlength) {
  T new;
  int querystart, queryend;
  Univcoord_T genomicstart, alignstart, alignend;
  /* Univcoord_T genomicend; */
  Endtype_T start_endtype, end_endtype;
  Trimaction_T trim_left_action, trim_right_action;

  /* Previously checked if left >= chroffset + chrlength to exclude
     the duplicate length, but now excluding all translocations to
     circular chromosomes */

  if (chroffset + chrlength < chrhigh) {
    /* Don't splice to circular chromosomes */
    return (T) NULL;

  } else if (plusp == true) {
    genomicstart = left;
    /* genomicend = left + querylength; */
    if (sensedir == SENSE_FORWARD) {
      start_endtype = END;
      end_endtype = DON;

      queryend = donor_pos;
#if 0
      querystart = substring_querystart; /* 0, for an end piece */
      if (querystart == 0) {
	trim_left_action = COMPUTE_TRIM;	/* querystart == 0 */
      } else {
	trim_left_action = PRE_TRIMMED;
      }
#else
      querystart = 0;
      trim_left_action = COMPUTE_TRIM;
#endif
      trim_right_action = NO_TRIM;

    } else if (sensedir == SENSE_ANTI) {
      start_endtype = DON;
      end_endtype = END;

      querystart = donor_pos;
      trim_left_action = NO_TRIM;
#if 0
      queryend = substring_queryend; /* querylength, for an end piece */
      if (queryend == querylength) {
	trim_right_action = COMPUTE_TRIM;	/* queryend == querylength */
      } else {
	trim_right_action = PRE_TRIMMED;
      }
#else
      queryend = querylength;
      trim_right_action = COMPUTE_TRIM;
#endif

    } else {
      abort();
    }

    alignstart = genomicstart + querystart;
    alignend = genomicstart + queryend;

  } else {
    genomicstart = left + querylength;
    /* genomicend = left; */
    if (sensedir == SENSE_FORWARD) {
      start_endtype = END;
      end_endtype = DON;

      queryend = querylength - donor_pos;
#if 0
      querystart = substring_querystart; /* 0, for an end piece */
      if (querystart == 0) {
	trim_left_action = COMPUTE_TRIM;	/* querystart == 0 */
      } else {
	trim_left_action = PRE_TRIMMED;
      }
#else
      querystart = 0;
      trim_left_action = COMPUTE_TRIM;
#endif
      trim_right_action = NO_TRIM;

    } else if (sensedir == SENSE_ANTI) {
      start_endtype = DON;
      end_endtype = END;

      querystart = querylength - donor_pos;
      trim_left_action = NO_TRIM;
#if 0
      queryend = substring_queryend; /* querylength, for an end piece */
      if (queryend == querylength) {
	trim_right_action = COMPUTE_TRIM;	/* queryend == querylength */
      } else {
	trim_right_action = PRE_TRIMMED;
      }
#else
      queryend = querylength;
      trim_right_action = COMPUTE_TRIM;
#endif

    } else {
      abort();
    }

    alignstart = genomicstart - querystart;
    alignend = genomicstart - queryend;
  }

  if ((new = Substring_new(donor_nmismatches,chrnum,chroffset,chrhigh,chrlength,
			   query_compress,start_endtype,end_endtype,querystart,queryend,querylength,
			   alignstart,alignend,/*exactp*/false,plusp,genestrand,
			   trim_left_action,trim_right_action,/*outofbounds_start*/0,/*outofbounds_end*/0,
			   /*minlength*/0,sensedir)) == NULL) {
    return (T) NULL;
  }

  debug2(printf("Making new donor with splicesites_i %d, coord %u and left %u, plusp %d, sensedir %d, query %d..%d, trim %d..%d, genome %u..%u\n",
		donor_knowni,donor_coord,left,plusp,sensedir,new->querystart,new->queryend,
		new->trim_left,new->trim_right,alignstart - chroffset,alignend - chroffset));
  debug2(printf("Original bounds were %d..%d\n",substring_querystart,substring_queryend));
  debug2(printf("Setting siteD_prob to be %f\n",donor_prob));

  new->splicecoord_D = donor_coord;
  new->splicesitesD_knowni = donor_knowni;
  assert(donor_coord == left + donor_pos);

  if (new->sensedir == SENSE_NULL) {
    new->sensedir = sensedir;
  }
  if (donor_knowni >= 0) {
    new->siteD_knownp = true;
    /* new->chimera_novelp = false */
  } else {
    /* new->siteD_knownp = false; */
    new->siteD_novelp = true;
  }

  if (plusp == true) {
    new->siteD_pos = donor_pos;
  } else {
    new->siteD_pos = querylength - donor_pos;
  }
  new->siteD_prob = donor_prob;

  return new;
}


T
Substring_new_acceptor (Univcoord_T acceptor_coord, int acceptor_knowni, int acceptor_pos,
			int substring_querystart, int substring_queryend, int acceptor_nmismatches,
			double acceptor_prob, Univcoord_T left, Compress_T query_compress,
			int querylength, bool plusp, int genestrand, int sensedir,
			Chrnum_T chrnum, Univcoord_T chroffset, Univcoord_T chrhigh, Chrpos_T chrlength) {
  T new;
  int querystart, queryend;
  Univcoord_T genomicstart, alignstart, alignend;
  /* Univcoord_T genomicend; */
  Endtype_T start_endtype, end_endtype;
  Trimaction_T trim_left_action, trim_right_action;

  /* Previously checked if left >= chroffset + chrlength to exclude
     the duplicate length, but now excluding all translocations to
     circular chromosomes */

  if (chroffset + chrlength < chrhigh) {
    /* Don't splice to circular chromosomes */
    return (T) NULL;

  } else if (plusp == true) {
    genomicstart = left;
    /* genomicend = left + querylength; */
    if (sensedir == SENSE_FORWARD) {
      start_endtype = ACC;
      end_endtype = END;

      querystart = acceptor_pos;
      trim_left_action = NO_TRIM;
#if 0
      queryend = substring_queryend; /* querylength, for an end piece */
      if (queryend == querylength) {
	trim_right_action = COMPUTE_TRIM;	/* queryend == querylength */
      } else {
	trim_right_action = PRE_TRIMMED;
      }
#else
      queryend = querylength;
      trim_right_action = COMPUTE_TRIM;
#endif

    } else if (sensedir == SENSE_ANTI) {
      start_endtype = END;
      end_endtype = ACC;

      queryend = acceptor_pos;
#if 0
      querystart = substring_querystart; /* 0, for an end piece */
      if (querystart == 0) {
	trim_left_action = COMPUTE_TRIM;	/* querystart == 0 */
      } else {
	trim_left_action = PRE_TRIMMED;
      }
#else
      querystart = 0;
      trim_left_action = COMPUTE_TRIM;
#endif
      trim_right_action = NO_TRIM;

    } else {
      abort();
    }

    alignstart = genomicstart + querystart;
    alignend = genomicstart + queryend;

  } else {
    genomicstart = left + querylength;
    /* genomicend = left; */
    if (sensedir == SENSE_FORWARD) {
      start_endtype = ACC;
      end_endtype = END;

      querystart = querylength - acceptor_pos;
      trim_left_action = NO_TRIM;
#if 0
      queryend = substring_queryend; /* querylength, for an end piece */
      if (queryend == querylength) {
	trim_right_action = COMPUTE_TRIM;	/* queryend == querylength */
      } else {
	trim_right_action = PRE_TRIMMED;
      }
#else
      queryend = querylength;
      trim_right_action = COMPUTE_TRIM;
#endif

    } else if (sensedir == SENSE_ANTI) {
      start_endtype = END;
      end_endtype = ACC;

      queryend = querylength - acceptor_pos;
#if 0
      querystart = substring_querystart; /* 0, for an end piece */
      if (querystart == 0) {
	trim_left_action = COMPUTE_TRIM;	/* querystart == 0 */
      } else {
	trim_left_action = PRE_TRIMMED;
      }
#else
      querystart = 0;
      trim_left_action = COMPUTE_TRIM;
#endif
      trim_right_action = NO_TRIM;

    } else {
      abort();
    }

    alignstart = genomicstart - querystart;
    alignend = genomicstart - queryend;
  }

  if ((new = Substring_new(acceptor_nmismatches,chrnum,chroffset,chrhigh,chrlength,
			   query_compress,start_endtype,end_endtype,querystart,queryend,querylength,
			   alignstart,alignend,/*exactp*/false,plusp,genestrand,
			   trim_left_action,trim_right_action,/*outofbounds_start*/0,/*outofbounds_end*/0,
			   /*minlength*/0,sensedir)) == NULL) {
    return (T) NULL;
  }

  debug2(printf("Making new acceptor with splicesites_i %d, coord %u and left %u, plusp %d, sensedir %d, query %d..%d, trim %d..%d, genome %u..%u\n",
		acceptor_knowni,acceptor_coord,left,plusp,sensedir,new->querystart,new->queryend,
		new->trim_left,new->trim_right,alignstart - chroffset,alignend - chroffset));
  debug2(printf("Original bounds were %d..%d\n",substring_querystart,substring_queryend));
  debug2(printf("Setting siteA_prob to be %f\n",acceptor_prob));

  new->splicecoord_A = acceptor_coord;
  new->splicesitesA_knowni = acceptor_knowni;
  assert(acceptor_coord == left + acceptor_pos);

  if (new->sensedir == SENSE_NULL) {
    new->sensedir = sensedir;
  }
  if (acceptor_knowni >= 0) {
    new->siteA_knownp = true;
    /* new->chimera_novelp = false */
  } else {
    /* new->chimera_knownp = false; */
    new->siteA_novelp = true;
  }

  if (plusp == true) {
    new->siteA_pos = acceptor_pos;
  } else {
    new->siteA_pos = querylength - acceptor_pos;
  }
  new->siteA_prob = acceptor_prob;

  return new;
}



T
Substring_new_shortexon (Univcoord_T acceptor_coord, int acceptor_knowni, Univcoord_T donor_coord, int donor_knowni,
			 int acceptor_pos, int donor_pos, int nmismatches,
			 double acceptor_prob, double donor_prob, Univcoord_T left,
			 Compress_T query_compress, int querylength,
			 bool plusp, int genestrand, int sensedir,
			 bool acceptor_ambp, bool donor_ambp,
			 Chrnum_T chrnum, Univcoord_T chroffset, Univcoord_T chrhigh, Chrpos_T chrlength) {
  T new;
  int querystart, queryend;
  Univcoord_T genomicstart, alignstart, alignend;
  Endtype_T start_endtype, end_endtype;

  if (plusp == true) {
    genomicstart = left;
    /* genomicend = left + querylength; */
    if (sensedir == SENSE_FORWARD) {
      start_endtype = (acceptor_ambp == true) ? AMB_ACC : ACC;
      end_endtype = (donor_ambp == true) ? AMB_DON : DON;
      querystart = acceptor_pos;
      queryend = donor_pos;
      alignstart = genomicstart + acceptor_pos;
      alignend = genomicstart + donor_pos;

    } else if (sensedir == SENSE_ANTI) {
      start_endtype = (donor_ambp == true) ? AMB_DON : DON;
      end_endtype = (acceptor_ambp == true) ? AMB_ACC : ACC;
      querystart = donor_pos;
      queryend = acceptor_pos;
      alignstart = genomicstart + donor_pos;
      alignend = genomicstart + acceptor_pos;

    } else {
      abort();
    }

  } else {
    genomicstart = left + querylength;
    /* genomicend = left; */
    if (sensedir == SENSE_FORWARD) {
      start_endtype = (acceptor_ambp == true) ? AMB_ACC : ACC;
      end_endtype = (donor_ambp == true) ? AMB_DON : DON;
      querystart = querylength - acceptor_pos;
      queryend = querylength - donor_pos;
      alignstart = genomicstart - (querylength - acceptor_pos);
      alignend = genomicstart - (querylength - donor_pos);

    } else if (sensedir == SENSE_ANTI) {
      start_endtype = (donor_ambp == true) ? AMB_DON : DON;
      end_endtype = (acceptor_ambp == true) ? AMB_ACC : ACC;
      querystart = querylength - donor_pos;
      queryend = querylength - acceptor_pos;
      alignstart = genomicstart - (querylength - donor_pos);
      alignend = genomicstart - (querylength - acceptor_pos);

    } else {
      abort();
    }
  }

  if ((new = Substring_new(nmismatches,chrnum,chroffset,chrhigh,chrlength,
			   query_compress,start_endtype,end_endtype,querystart,queryend,querylength,
			   alignstart,alignend,/*exactp*/false,plusp,genestrand,
			   /*trim_left_action*/NO_TRIM,/*trim_right_action*/NO_TRIM,
			   /*outofbounds_start*/0,/*outofbounds_end*/0,/*minlength*/0,sensedir)) == NULL) {
    return (T) NULL;
  }

  debug2(printf("Making new middle with left %u, plusp %d\n",left,plusp));
  new->splicecoord_A = acceptor_coord;
  new->splicesitesA_knowni = acceptor_knowni;
  new->splicecoord_D = donor_coord;
  new->splicesitesD_knowni = donor_knowni;

  if (new->sensedir == SENSE_NULL) {
    new->sensedir = sensedir;
  }

  if (acceptor_knowni >= 0) {
    new->siteA_knownp = true;
    /* new->chimera_novelp = false; */
  } else {
    /* new->chimera_knownp = false; */
    new->siteA_novelp = true;
  }

  if (donor_knowni >= 0) {
    new->siteD_knownp = true;
    /* new->chimera_novelp_2 = false; */
  } else {
    /* new->siteD_knownp_2 = false; */
    new->siteD_novelp = true;
  }

  if (plusp == true) {
    new->siteA_pos = acceptor_pos;
    new->siteD_pos = donor_pos;
  } else {
    new->siteA_pos = querylength - acceptor_pos;
    new->siteD_pos = querylength - donor_pos;
  }

  new->siteA_prob = acceptor_prob;
  new->siteD_prob = donor_prob;

  return new;
}



void
Substring_assign_donor_prob (T donor) {

  if (donor == NULL) {
    return;

  } else if (donor->siteD_knownp == false) {
    /* Prob already assigned */

  } else if (donor->sensedir == SENSE_FORWARD) {
    if (donor->plusp == true) {
      donor->siteD_prob = Maxent_hr_donor_prob(donor->splicecoord_D,donor->chroffset);
    } else {
      donor->siteD_prob = Maxent_hr_antidonor_prob(donor->splicecoord_D,donor->chroffset);
    }

  } else if (donor->sensedir == SENSE_ANTI) {
    if (donor->plusp == true) {
      donor->siteD_prob = Maxent_hr_antidonor_prob(donor->splicecoord_D,donor->chroffset);
    } else {
      donor->siteD_prob = Maxent_hr_donor_prob(donor->splicecoord_D,donor->chroffset);
    }

  } else {
    /* SENSE_NULL */
    donor->siteD_prob = 0.0;
  }

  return;
}

void
Substring_assign_acceptor_prob (T acceptor) {

  if (acceptor == NULL) {
    return;

  } else if (acceptor->siteA_knownp == false) {
    /* Prob already assigned */

  } else if (acceptor->sensedir == SENSE_FORWARD) {
    if (acceptor->plusp == true) {
      acceptor->siteA_prob = Maxent_hr_acceptor_prob(acceptor->splicecoord_A,acceptor->chroffset);
    } else {
      acceptor->siteA_prob = Maxent_hr_antiacceptor_prob(acceptor->splicecoord_A,acceptor->chroffset);
    }
    
  } else if (acceptor->sensedir == SENSE_ANTI) {
    if (acceptor->plusp == true) {
      acceptor->siteA_prob = Maxent_hr_antiacceptor_prob(acceptor->splicecoord_A,acceptor->chroffset);
    } else {
      acceptor->siteA_prob = Maxent_hr_acceptor_prob(acceptor->splicecoord_A,acceptor->chroffset);
    }

  } else {
    /* SENSE_NULL */
    acceptor->siteA_prob = 0.0;
  }

  return;
}


void
Substring_assign_shortexon_prob (T shortexon) {

  if (shortexon->siteA_knownp == false) {
    /* Prob1 already assigned */

  } else if (shortexon->sensedir == SENSE_FORWARD) {
    if (shortexon->plusp == true) {
      shortexon->siteA_prob = Maxent_hr_acceptor_prob(shortexon->splicecoord_A,shortexon->chroffset);
    } else {
      shortexon->siteA_prob = Maxent_hr_antiacceptor_prob(shortexon->splicecoord_A,shortexon->chroffset);
    }

  } else if (shortexon->sensedir == SENSE_ANTI) {
    if (shortexon->plusp == true) {
      shortexon->siteA_prob = Maxent_hr_antiacceptor_prob(shortexon->splicecoord_A,shortexon->chroffset);
    } else {
      shortexon->siteA_prob = Maxent_hr_acceptor_prob(shortexon->splicecoord_A,shortexon->chroffset);
    }

  } else {
    abort();
  }

  if (shortexon->siteD_knownp == false) {
    /* Prob2 already assigned */

  } else if (shortexon->sensedir == SENSE_FORWARD) {
    if (shortexon->plusp == true) {
      shortexon->siteD_prob = Maxent_hr_donor_prob(shortexon->splicecoord_D,shortexon->chroffset);
    } else {
      shortexon->siteD_prob = Maxent_hr_antidonor_prob(shortexon->splicecoord_D,shortexon->chroffset);
    }

  } else if (shortexon->sensedir == SENSE_ANTI) {
    if (shortexon->plusp == true) {
      shortexon->siteD_prob = Maxent_hr_antidonor_prob(shortexon->splicecoord_D,shortexon->chroffset);
    } else {
      shortexon->siteD_prob = Maxent_hr_donor_prob(shortexon->splicecoord_D,shortexon->chroffset);
    }

  } else {
    abort();
  }

  return;
}



static int
ascending_siteD_pos_cmp (const void *a, const void *b) {
  T x = * (T *) a;
  T y = * (T *) b;

  if (x->siteD_pos < y->siteD_pos) {
    return -1;
  } else if (x->siteD_pos > y->siteD_pos) {
    return +1;
  } else if (x->genomicstart < y->genomicstart) {
    return -1;
  } else if (x->genomicstart > y->genomicstart) {
    return +1;
  } else if (x->siteD_knownp == true && y->siteD_knownp == false) {
    return -1;
  } else if (y->siteD_knownp == true && x->siteD_knownp == false) {
    return +1;
  } else {
    return 0;
  }
}

static int
ascending_siteA_pos_cmp (const void *a, const void *b) {
  T x = * (T *) a;
  T y = * (T *) b;

  if (x->siteA_pos < y->siteA_pos) {
    return -1;
  } else if (x->siteA_pos > y->siteA_pos) {
    return +1;
  } else if (x->genomicstart < y->genomicstart) {
    return -1;
  } else if (x->genomicstart > y->genomicstart) {
    return +1;
  } else if (x->siteA_knownp == true && y->siteA_knownp == false) {
    return -1;
  } else if (y->siteA_knownp == true && x->siteA_knownp == false) {
    return +1;
  } else {
    return 0;
  }
}

static int
ascending_siteN_pos_cmp (const void *a, const void *b) {
  T x = * (T *) a;
  T y = * (T *) b;

  if (x->siteN_pos < y->siteN_pos) {
    return -1;
  } else if (x->siteN_pos > y->siteN_pos) {
    return +1;
  } else if (x->genomicstart < y->genomicstart) {
    return -1;
  } else if (x->genomicstart > y->genomicstart) {
    return +1;
  } else {
    return 0;
  }
}

static int
descending_siteD_pos_cmp (const void *a, const void *b) {
  T x = * (T *) a;
  T y = * (T *) b;

  if (x->siteD_pos < y->siteD_pos) {
    return -1;
  } else if (x->siteD_pos > y->siteD_pos) {
    return +1;
  } else if (x->genomicstart > y->genomicstart) {
    return -1;
  } else if (x->genomicstart < y->genomicstart) {
    return +1;
  } else if (x->siteD_knownp == true && y->siteD_knownp == false) {
    return -1;
  } else if (y->siteD_knownp == true && x->siteD_knownp == false) {
    return +1;
  } else {
    return 0;
  }
}

static int
descending_siteA_pos_cmp (const void *a, const void *b) {
  T x = * (T *) a;
  T y = * (T *) b;

  if (x->siteA_pos < y->siteA_pos) {
    return -1;
  } else if (x->siteA_pos > y->siteA_pos) {
    return +1;
  } else if (x->genomicstart > y->genomicstart) {
    return -1;
  } else if (x->genomicstart < y->genomicstart) {
    return +1;
  } else if (x->siteA_knownp == true && y->siteA_knownp == false) {
    return -1;
  } else if (y->siteA_knownp == true && x->siteA_knownp == false) {
    return +1;
  } else {
    return 0;
  }
}

static int
descending_siteN_pos_cmp (const void *a, const void *b) {
  T x = * (T *) a;
  T y = * (T *) b;

  if (x->siteN_pos < y->siteN_pos) {
    return -1;
  } else if (x->siteN_pos > y->siteN_pos) {
    return +1;
  } else if (x->genomicstart > y->genomicstart) {
    return -1;
  } else if (x->genomicstart < y->genomicstart) {
    return +1;
  } else {
    return 0;
  }
}

List_T
Substring_sort_siteD_halves (List_T hitlist, bool ascendingp) {
  List_T sorted = NULL;
  T x, *hits;
  int n, i, j;
  bool *eliminate;

  n = List_length(hitlist);
  debug(printf("Checking %d spliceends for duplicates...",n));
  if (n == 0) {
    debug(printf("\n"));
    return NULL;
  }

  hits = (T *) MALLOCA(n * sizeof(T));
  List_fill_array_and_free((void **) hits,&hitlist);

  if (ascendingp == true) {
    qsort(hits,n,sizeof(T),ascending_siteD_pos_cmp);
  } else {
    qsort(hits,n,sizeof(T),descending_siteD_pos_cmp);
  }

  /* Check for duplicates */
  eliminate = (bool *) CALLOCA(n,sizeof(bool));
  for (i = 0; i < n; i++) {
    x = hits[i];
    j = i+1;
    while (j < n && hits[j]->siteD_pos == x->siteD_pos && hits[j]->genomicstart == x->genomicstart) {
      eliminate[j] = true;
      j++;
    }
  }

  debug(j = 0);
  for (i = n-1; i >= 0; i--) {
    x = hits[i];
    if (eliminate[i] == false) {
      sorted = List_push(sorted,x);
    } else {
      Substring_free(&x);
      debug(j++);
    }
  }
  debug(printf("%d eliminated\n",j));

  FREEA(hits);
  FREEA(eliminate);

  return sorted;
}

List_T
Substring_sort_siteA_halves (List_T hitlist, bool ascendingp) {
  List_T sorted = NULL;
  T x, *hits;
  int n, i, j;
  bool *eliminate;

  n = List_length(hitlist);
  debug(printf("Checking %d spliceends for duplicates...",n));
  if (n == 0) {
    debug(printf("\n"));
    return NULL;
  }

  hits = (T *) MALLOCA(n * sizeof(T));
  List_fill_array_and_free((void **) hits,&hitlist);

  if (ascendingp == true) {
    qsort(hits,n,sizeof(T),ascending_siteA_pos_cmp);
  } else {
    qsort(hits,n,sizeof(T),descending_siteA_pos_cmp);
  }

  /* Check for duplicates */
  eliminate = (bool *) CALLOCA(n,sizeof(bool));
  for (i = 0; i < n; i++) {
    x = hits[i];
    j = i+1;
    while (j < n && hits[j]->siteA_pos == x->siteA_pos && hits[j]->genomicstart == x->genomicstart) {
      eliminate[j] = true;
      j++;
    }
  }

  debug(j = 0);
  for (i = n-1; i >= 0; i--) {
    x = hits[i];
    if (eliminate[i] == false) {
      sorted = List_push(sorted,x);
    } else {
      Substring_free(&x);
      debug(j++);
    }
  }
  debug(printf("%d eliminated\n",j));

  FREEA(hits);
  FREEA(eliminate);

  return sorted;
}

List_T
Substring_sort_siteN_halves (List_T hitlist, bool ascendingp) {
  List_T sorted = NULL;
  T x, *hits;
  int n, i, j;
  bool *eliminate;

  n = List_length(hitlist);
  debug(printf("Checking %d spliceends for duplicates...",n));
  if (n == 0) {
    debug(printf("\n"));
    return NULL;
  }

  hits = (T *) MALLOCA(n * sizeof(T));
  List_fill_array_and_free((void **) hits,&hitlist);

  if (ascendingp == true) {
    qsort(hits,n,sizeof(T),ascending_siteN_pos_cmp);
  } else {
    qsort(hits,n,sizeof(T),descending_siteN_pos_cmp);
  }

  /* Check for duplicates */
  eliminate = (bool *) CALLOCA(n,sizeof(bool));
  for (i = 0; i < n; i++) {
    x = hits[i];
    j = i+1;
    while (j < n && hits[j]->siteN_pos == x->siteN_pos && hits[j]->genomicstart == x->genomicstart) {
      eliminate[j] = true;
      j++;
    }
  }

  debug(j = 0);
  for (i = n-1; i >= 0; i--) {
    x = hits[i];
    if (eliminate[i] == false) {
      sorted = List_push(sorted,x);
    } else {
      Substring_free(&x);
      debug(j++);
    }
  }
  debug(printf("%d eliminated\n",j));

  FREEA(hits);
  FREEA(eliminate);

  return sorted;
}


static void
print_snp_labels (Filestring_T fp, T this, Shortread_T queryseq) {
  int *snps, nsnps, querypos, i;
  char *label, *seq1, *seq2;
  bool allocp, printp = false;
  Interval_T interval;
  Chrpos_T position;

  if (this->plusp == true) {
    snps = IIT_get_with_divno(&nsnps,snps_iit,
			      snps_divint_crosstable[this->chrnum],
			      this->alignstart - this->chroffset + 1,this->alignend - this->chroffset,
			      /*sortp*/false);
  } else {
    snps = IIT_get_with_divno(&nsnps,snps_iit,
			      snps_divint_crosstable[this->chrnum],
			      this->alignend - this->chroffset + 1,this->alignstart - this->chroffset,
			      /*sortp*/false);
  }

  FPRINTF(fp,",snps:");

  seq1 = Shortread_fullpointer_uc(queryseq);
  if (this->genomic_bothdiff == NULL) {
    seq2 = Shortread_fullpointer(queryseq);
  } else {
    seq2 = this->genomic_bothdiff;
  }

  if (this->plusp) {

#if 0
    for (i = 0; i < nsnps; i++) {
      interval = IIT_interval(snps_iit,snps[i]);
      position = Interval_low(interval);
      querypos = position - (this->genomicstart-this->chroffset) - 1;
      printf("%d ",querypos);
    }
    printf("\n");
#endif

    for (i = 0; i < nsnps; i++) {
      interval = IIT_interval(snps_iit,snps[i]);
      position = Interval_low(interval);
      querypos = position - (this->genomicstart-this->chroffset) - 1;
      /* assert(querypos >= 0 && querypos < this->genomiclength); */

#if 0
      alleles = IIT_typestring(snps_iit,Interval_type(interval));
      if (c == alleles[0] || c == alleles[1]) ;
#endif

      if (isupper(seq2[querypos]) && seq1[querypos] != seq2[querypos]) {
	label = IIT_label(snps_iit,snps[i],&allocp);
	if (printp) {
	  FPRINTF(fp,"|");
	}
	FPRINTF(fp,"%d@",querypos+1);
	FPRINTF(fp,"%s",label);
	printp = true;
	if (allocp) FREE(label);
      }
    }

  } else {

    for (i = nsnps-1; i >= 0; i--) {
      interval = IIT_interval(snps_iit,snps[i]);
      position = Interval_low(interval);
      querypos = (this->genomicstart-this->chroffset) - position;
      /* assert(querypos >= 0 && querypos < this->genomiclength); */

#if 0
      /* printf("\n%d%c\n",querypos,c); */
      alleles = IIT_typestring(snps_iit,Interval_type(interval));
      if (c == alleles[0] || c == alleles[1]) ;
#endif

      if (isupper(seq2[querypos]) && seq1[querypos] != seq2[querypos]) {
	label = IIT_label(snps_iit,snps[i],&allocp);
	if (printp) {
	  FPRINTF(fp,"|");
	}
	FPRINTF(fp,"%d@",querypos+1);
	FPRINTF(fp,"%s",label);
	printp = true;
	if (allocp) FREE(label);
      }
    }

  }

  FREE(snps);

  return;
}



bool
Substring_contains_known_splicesite (T this) {
  Univcoord_T left;
  bool method2;
  int length;
#ifdef DEBUG7
  Chrpos_T start, end;
  bool method1;
#endif

  if (splicesites_iit == NULL) {
    return false;
  } else if (this->plusp == true) {
#ifdef DEBUG7
    start = this->alignstart+KNOWN_SPLICESITE_EDGE - this->chroffset + 1;
    end = this->alignend-KNOWN_SPLICESITE_EDGE - this->chroffset;
    method1 = IIT_exists_with_divno(splicesites_iit,splicesites_divint_crosstable[this->chrnum],start,end);
    debug7(printf("Checking if #%d:%u..%u has a known splice site: %d\n",this->chrnum,start,end,method1));
#endif

    left = this->alignstart;
    length = this->alignend - left;
    method2 = Splicetrie_splicesite_p(left,1,length-1);
    debug7(printf("Checking if #%d:%u (%d..%d) has a known splice site: %d\n",
		  this->chrnum,left - this->chroffset,1,length-1,method2));

    /* assert(method2 == method1); */
    /* if (method2 == true) exit(0); */
    return method2;

  } else {
#ifdef DEBUG7
    start = this->alignend+KNOWN_SPLICESITE_EDGE - this->chroffset + 1;
    end = this->alignstart-KNOWN_SPLICESITE_EDGE - this->chroffset;
    method1 = IIT_exists_with_divno(splicesites_iit,splicesites_divint_crosstable[this->chrnum],start,end);
    debug7(printf("Checking if #%d:%u..%u has a known splice site: %d\n",this->chrnum,start,end,method1));
#endif

    left = this->alignend;
    length = this->alignstart - left;
    method2 = Splicetrie_splicesite_p(left,1,length-1);
    debug7(printf("Checking if #%d:%u (%d..%d) has a known splice site: %d\n",
		  this->chrnum,left - this->chroffset,1,length-1,method2));

    /* assert(method2 == method1); */
    /* if (method2 == true) exit(0); */
    return method2;
  }
}



/* Has a copy in pair.c */
static void
print_splicesite_labels (Filestring_T fp, T this, int typeint, int chimera_pos, char *tag) {
  Chrpos_T splicesitepos;
  int *splicesites, nsplicesites, i;
  char *label;
  bool allocp;

  if (this->plusp == true) {
    splicesitepos = this->genomicstart - this->chroffset + chimera_pos;
  } else {
    splicesitepos = this->genomicstart - this->chroffset - chimera_pos;
  }

  /* Note: this->chimera_knownp might not be set for GMAP alignments */
  if (this->siteD_knownp == true) {
    /* Note: IIT_get_typed_signed_with_divno does not work here */
    splicesites = IIT_get_exact_multiple_with_divno(&nsplicesites,splicesites_iit,
						    splicesites_divint_crosstable[this->chrnum],
						    splicesitepos,splicesitepos+1U,typeint);
    if (nsplicesites == 0) {
#if 0
      /* Can happen if splicesites_iit are for introns, and not splicesites */
      fprintf(stderr,"Supposed to have a splicesite at chrnum %d, %u..%u, type %d\n",
	      this->chrnum,splicesitepos,splicesitepos+1U,typeint);
#endif
    } else {
      FPRINTF(fp,",%s:",tag);
      label = IIT_label(splicesites_iit,splicesites[0],&allocp);
      FPRINTF(fp,"%s",label);
      if (allocp) FREE(label);

      for (i = 1; i < nsplicesites; i++) {
	label = IIT_label(splicesites_iit,splicesites[i],&allocp);
	FPRINTF(fp,"|%s",label);
	if (allocp) FREE(label);
      }
      FREE(splicesites);
    }
  }

  return;
}


#if 0
/* replaced by intragenic_splice_p in stage1hr.c */
/* donor_genomicpos and acceptor_genomicpos are Univcoord_T */
bool
Substring_intragenic_splice_p (Genomicpos_T splicedistance, Chrnum_T chrnum,
			       Genomicpos_T donor_genomicpos, Genomicpos_T acceptor_genomicpos,
			       Genomicpos_T chroffset) {
  Genomicpos_T obsdistance;
  Genomicpos_T chrpos;
  int *splicesites, nsplicesites, i;
  char *annot, *restofheader;
  bool alloc_header_p;

  debug5(printf("Entered Substring_intragenic_splice_p with donor %u and acceptor %u, chroffset %u => chrpos %u and %u\n",
		donor_genomicpos,acceptor_genomicpos,chroffset,donor_genomicpos-chroffset,acceptor_genomicpos-chroffset));

  chrpos = donor_genomicpos - chroffset;
  splicesites = IIT_get_typed_with_divno(&nsplicesites,splicesites_iit,
					 splicesites_divint_crosstable[chrnum],
					 chrpos,chrpos,donor_typeint,/*sortp*/false);
  for (i = 0; i < nsplicesites; i++) {
    annot = IIT_annotation(&restofheader,splicesites_iit,splicesites[i],&alloc_header_p);
    debug5(printf("Comparing obsdistance %s at donor %u with splicedistance %u\n",restofheader,chrpos,splicedistance));
    if (sscanf(restofheader,"%u",&obsdistance) == 1 &&
	obsdistance >= splicedistance) {
      if (alloc_header_p == true) {
	FREE(restofheader);
      }
      FREE(splicesites);
      return true;
    }
    if (alloc_header_p == true) {
      FREE(restofheader);
    }
  }
  FREE(splicesites);

  chrpos = acceptor_genomicpos - chroffset;
  splicesites = IIT_get_typed_with_divno(&nsplicesites,splicesites_iit,
					 splicesites_divint_crosstable[chrnum],
					 acceptor_genomicpos,acceptor_genomicpos,donor_typeint,
					 /*sortp*/false);
  for (i = 0; i < nsplicesites; i++) {
    annot = IIT_annotation(&restofheader,splicesites_iit,splicesites[i],&alloc_header_p);
    debug5(printf("Comparing obsdistance %s at acceptor %u with splicedistance %u\n",restofheader,chrpos,splicedistance));
    if (sscanf(restofheader,"%u",&obsdistance) == 1 &&
	obsdistance >= splicedistance) {
      if (alloc_header_p == true) {
	FREE(restofheader);
      }
      FREE(splicesites);
      return true;
    }
    if (alloc_header_p == true) {
      FREE(restofheader);
    }
  }
  FREE(splicesites);

  return false;
}
#endif


Chrpos_T
Substring_compute_chrpos (T this, int hardclip_low, bool hide_soft_clips_p) {
  Chrpos_T chrpos;

  if (hide_soft_clips_p == true) {
    if (this->plusp == true) {
      /* Add 1 to report in 1-based coordinates */
      chrpos = this->genomicstart - this->chroffset + 1U;
      chrpos += hardclip_low;

    } else {
      /* Add 1 to report in 1-based coordinates */
      chrpos = this->genomicend - this->chroffset + 1U;
      chrpos += hardclip_low;
    }

  } else {
    if (this->plusp == true) {
      chrpos = this->genomicstart - this->chroffset + 1U;
      if (this->querystart > hardclip_low) {
	chrpos += this->querystart; /* not querystart_orig */
      } else {
	chrpos += hardclip_low;
      }

    } else {
      chrpos = this->genomicend - this->chroffset + 1U;
      if (this->querylength - this->queryend > hardclip_low) {
	chrpos += this->querylength - this->queryend; /* not queryend_orig */
      } else {
	chrpos += hardclip_low;
      }
    }
  }
    
  return chrpos;
}



/* Taken from NCBI Blast 2.2.29, algo/blast/core/blast_stat.c */
/* Karlin-Altschul formula: m n exp(-lambda * S + log k) = k m n exp(-lambda * S) */
/* Also in pair.c */

static double
blast_evalue (int alignlength, int nmismatches) {
  double k = 0.1;
  double lambda = 1.58;		/* For a +1, -1 scoring scheme */
  double score;
  
  score = (double) ((alignlength - nmismatches) /* scored as +1 */ - nmismatches /* scored as -1 */);

  return k * (double) alignlength * genomelength * exp(-lambda * score);
}

static double
blast_bitscore (int alignlength, int nmismatches) {
  double k = 0.1;
  double lambda = 1.58;		/* For a +1, -1 scoring scheme */
  double score;
  
  score = (double) ((alignlength - nmismatches) /* scored as +1 */ - nmismatches /* scored as -1 */);
  return (score * lambda - log(k)) / log(2.0);
}


double
Substring_evalue (T substring) {
  int alignlength_trim;

#if 0
  /* Doesn't work for ambiguous substrings */
  if (substring->plusp == true) {
    alignlength_trim = (int) (substring->alignend_trim - substring->alignstart_trim);
  } else {
    alignlength_trim = (int) (substring->alignstart_trim - substring->alignend_trim);
  }
  assert(alignlength_trim == substring->queryend - substring->querystart);
#else
  alignlength_trim = substring->queryend - substring->querystart;
#endif

  return blast_evalue(alignlength_trim,substring->nmismatches_bothdiff);
}


void
Substring_print_m8 (Filestring_T fp, T substring, Shortread_T headerseq, char *acc_suffix,
		    char *chr, bool invertp) {
  double identity;
  int alignlength_trim;

  FPRINTF(fp,"%s%s",Shortread_accession(headerseq),acc_suffix); /* field 0: accession */

  FPRINTF(fp,"\t%s",chr);	/* field 1: chr */

  /* field 2: identity */
  if (substring->plusp == true) {
    alignlength_trim = (int) (substring->alignend_trim - substring->alignstart_trim);
  } else {
    alignlength_trim = (int) (substring->alignstart_trim - substring->alignend_trim);
  }

  identity = (double) (alignlength_trim - substring->nmismatches_bothdiff)/(double) alignlength_trim;
  FPRINTF(fp,"\t%.1f",100.0*identity);


  FPRINTF(fp,"\t%d",alignlength_trim); /* field 3: query length */

  FPRINTF(fp,"\t%d",substring->nmismatches_bothdiff); /* field 4: nmismatches */

  FPRINTF(fp,"\t0");		/* field 5: gap openings */

  FPRINTF(fp,"\t%d",substring->querystart + 1); /* field 6: query start */

  FPRINTF(fp,"\t%d",substring->queryend); /* field 7: query end */

  /* fields 8 and 9: chr start and end */
  if (substring->plusp == true) {
    if (invertp == false) {
      FPRINTF(fp,"\t%u\t%u",substring->alignstart_trim - substring->chroffset + 1U,
	      substring->alignend_trim - substring->chroffset);
    } else {
      FPRINTF(fp,"\t%u\t%u",substring->alignend_trim - substring->chroffset,
	      substring->alignstart_trim - substring->chroffset + 1U);
    }
  } else {
    if (invertp == false) {
      FPRINTF(fp,"\t%u\t%u",substring->alignstart_trim - substring->chroffset,
	      substring->alignend_trim - substring->chroffset + 1U);
    } else {
      FPRINTF(fp,"\t%u\t%u",substring->alignend_trim - substring->chroffset + 1U,
	      substring->alignstart_trim - substring->chroffset);
    }
  }

  /* field 10: E value */
  FPRINTF(fp,"\t%.2g",blast_evalue(alignlength_trim,substring->nmismatches_bothdiff));

 /* field 11: bit score */
  FPRINTF(fp,"\t%.1f",blast_bitscore(alignlength_trim,substring->nmismatches_bothdiff));
  
  FPRINTF(fp,"\n");

  return;
}



static void
print_forward (Filestring_T fp, char *string, int n) {
  int i;
  
  FPRINTF(fp,"%.*s",n,string);
  return;
}




static void
print_lc (Filestring_T fp, char *string, int n) {
  int i;
  
  for (i = 0; i < n; i++) {
    FPRINTF(fp,"%c",(char) tolower(string[i]));
  }
  return;
}



static void
print_revcomp (Filestring_T fp, char *nt, int len) {
  int i;

  FPRINTF(fp,"%.*R",len,nt);
  return;
}

static void
print_revcomp_lc (Filestring_T fp, char *nt, int len) {
  int i;

  for (i = len-1; i >= 0; --i) {
    FPRINTF(fp,"%c",(char) tolower(complCode[(int) nt[i]]));
  }
  return;
}



static void
print_genomic (Filestring_T fp, T substring, char *deletion, int deletionlength, bool invertp,
	       Shortread_T queryseq) {
  int i;

  if (invertp == false) {
    if (substring->genomic_bothdiff == NULL) {
      /* Exact match */
      Shortread_print_oneline_uc(fp,queryseq);

    } else if (show_refdiff_p == true) {
      print_forward(fp,substring->genomic_refdiff,substring->queryend);
      if (deletion != NULL) {
	print_lc(fp,deletion,deletionlength);
      }
      print_forward(fp,&(substring->genomic_refdiff[substring->queryend]),substring->querylength - deletionlength - substring->queryend);
    } else {
      print_forward(fp,substring->genomic_bothdiff,substring->queryend);
      if (deletion != NULL) {
	print_lc(fp,deletion,deletionlength);
      }
      print_forward(fp,&(substring->genomic_bothdiff[substring->queryend]),substring->querylength - deletionlength - substring->queryend);
    }

    for (i = 0; i < Shortread_choplength(queryseq); i++) {
      FPRINTF(fp,"*");
    }
    FPRINTF(fp,"\t");
    FPRINTF(fp,"%d..%d",1 + substring->querystart,substring->queryend);

  } else {
    if (substring->genomic_bothdiff == NULL) {
      /* Exact match */
      Shortread_print_oneline_revcomp_uc(fp,queryseq);

    } else if (show_refdiff_p == true) {
      print_revcomp(fp,&(substring->genomic_refdiff[substring->querystart]),substring->querylength - substring->querystart);
      if (deletion != NULL) {
	print_revcomp_lc(fp,deletion,deletionlength);
      }
      print_revcomp(fp,&(substring->genomic_refdiff[deletionlength]),substring->querystart - deletionlength);

    } else {
      print_revcomp(fp,&(substring->genomic_bothdiff[substring->querystart]),substring->querylength - substring->querystart);
      if (deletion != NULL) {
	print_revcomp_lc(fp,deletion,deletionlength);
      }
      print_revcomp(fp,&(substring->genomic_bothdiff[deletionlength]),substring->querystart - deletionlength);
    }
    for (i = 0; i < Shortread_choplength(queryseq); i++) {
      FPRINTF(fp,"*");
    }
    FPRINTF(fp,"\t");
    FPRINTF(fp,"%d..%d",1 + substring->querylength - substring->queryend,
	   substring->querylength - substring->querystart);
  }
  return;
}


static void
print_coordinates (Filestring_T fp, T substring, char *chr, bool invertp) {

  if (substring->plusp == true) {
    if (invertp == false) {
      FPRINTF(fp,"+%s:%u..%u",chr,substring->alignstart_trim - substring->chroffset + 1U,
	     substring->alignend_trim - substring->chroffset);
    } else {
      FPRINTF(fp,"-%s:%u..%u",chr,substring->alignend_trim - substring->chroffset,
	     substring->alignstart_trim - substring->chroffset + 1U);
    }
  } else {
    if (invertp == false) {
      FPRINTF(fp,"-%s:%u..%u",chr,substring->alignstart_trim - substring->chroffset,
	     substring->alignend_trim - substring->chroffset + 1U);
    } else {
      FPRINTF(fp,"+%s:%u..%u",chr,substring->alignend_trim - substring->chroffset + 1U,
	     substring->alignstart_trim - substring->chroffset);
    }
  }

  return;
}



void
Substring_print_alignment (Filestring_T fp, Junction_T pre_junction, T substring, Junction_T post_junction,
			   Shortread_T queryseq, Genome_T genome, char *chr, bool invertp) {
  char *deletion_string;
  int deletion_length;
  Junctiontype_T type1, type2;
  Chrpos_T splice_distance_1, splice_distance_2;

  if (post_junction == NULL) {
    deletion_string = (char *) NULL;
    deletion_length = 0;
  } else if (Junction_type(post_junction) != DEL_JUNCTION) {
    deletion_string = (char *) NULL;
    deletion_length = 0;
  } else {
    deletion_string = Junction_deletion_string(post_junction,genome,substring->plusp);
    deletion_length = Junction_nindels(post_junction);
  }

  print_genomic(fp,substring,deletion_string,deletion_length,invertp,queryseq);
  FREE(deletion_string);
  FPRINTF(fp,"\t");
  print_coordinates(fp,substring,chr,invertp);

  FPRINTF(fp,"\t");
  if (pre_junction == NULL) {
    type1 = NO_JUNCTION;
    /* Handle result of substring_trim_novel_spliceends */
    if (invertp == false) {
      if (substring->start_endtype == DON) {
	FPRINTF(fp,"donor:%.2f",substring->siteD_prob);
      } else if (substring->start_endtype == ACC) {
	FPRINTF(fp,"acceptor:%.2f",substring->siteA_prob);
      } else {
	FPRINTF(fp,"start:%d",substring->trim_left);
      }
    } else {
      if (substring->end_endtype == DON) {
	FPRINTF(fp,"donor:%.2f",substring->siteD_prob);
      } else if (substring->end_endtype == ACC) {
	FPRINTF(fp,"acceptor:%.2f",substring->siteA_prob);
      } else {
	FPRINTF(fp,"start:%d",substring->trim_right);
      }
    }
  } else if ((type1 = Junction_type(pre_junction)) == INS_JUNCTION) {
    FPRINTF(fp,"ins:%d",Junction_nindels(pre_junction));
  } else if (type1 == DEL_JUNCTION) {
    FPRINTF(fp,"del:%d",Junction_nindels(pre_junction));
  } else if (type1 == SPLICE_JUNCTION) {
    if (invertp == false) {
      if (Junction_sensedir(pre_junction) == SENSE_ANTI) {
	FPRINTF(fp,"donor:%.2f",Junction_donor_prob(pre_junction));
      } else {
	FPRINTF(fp,"acceptor:%.2f",Junction_acceptor_prob(pre_junction));
      }
    } else {
      if (Junction_sensedir(pre_junction) == SENSE_ANTI) {
	FPRINTF(fp,"acceptor:%.2f",Junction_acceptor_prob(pre_junction));
      } else {
	FPRINTF(fp,"donor:%.2f",Junction_donor_prob(pre_junction));
      }
    }
  } else if (type1 == CHIMERA_JUNCTION) {
    FPRINTF(fp,"distant:%u",Junction_splice_distance(pre_junction));
  } else {
    abort();
  }

  FPRINTF(fp,"..");

  if (post_junction == NULL) {
    type2 = NO_JUNCTION;
    /* Handle result of substring_trim_novel_spliceends */
    if (invertp == false) {
      if (substring->end_endtype == DON) {
	FPRINTF(fp,"donor:%.2f",substring->siteD_prob);
      } else if (substring->end_endtype == ACC) {
	FPRINTF(fp,"acceptor:%.2f",substring->siteA_prob);
      } else {
	FPRINTF(fp,"end:%d",substring->trim_right);
      }
    } else {
      if (substring->start_endtype == DON) {
	FPRINTF(fp,"donor:%.2f",substring->siteD_prob);
      } else if (substring->start_endtype == ACC) {
	FPRINTF(fp,"acceptor:%.2f",substring->siteA_prob);
      } else {
	FPRINTF(fp,"end:%d",substring->trim_left);
      }
    }
  } else if ((type2 = Junction_type(post_junction)) == INS_JUNCTION) {
    FPRINTF(fp,"ins:%d",Junction_nindels(post_junction));
  } else if (type2 == DEL_JUNCTION) {
    FPRINTF(fp,"del:%d",Junction_nindels(post_junction));
  } else if (type2 == SPLICE_JUNCTION) {
    if (invertp == false) {
      if (Junction_sensedir(post_junction) == SENSE_ANTI) {
	FPRINTF(fp,"acceptor:%.2f",Junction_acceptor_prob(post_junction));
      } else {
	FPRINTF(fp,"donor:%.2f",Junction_donor_prob(post_junction));
      }
    } else {
      if (Junction_sensedir(post_junction) == SENSE_ANTI) {
	FPRINTF(fp,"donor:%.2f",Junction_donor_prob(post_junction));
      } else {
	FPRINTF(fp,"acceptor:%.2f",Junction_acceptor_prob(post_junction));
      }
    }
  } else if (type2 == CHIMERA_JUNCTION) {
    FPRINTF(fp,"distant:%u",Junction_splice_distance(post_junction));
  } else {
    abort();
  }

  FPRINTF(fp,",matches:%d,sub:%d",Substring_nmatches_posttrim(substring),substring->nmismatches_bothdiff);
  if (print_nsnpdiffs_p) {
    FPRINTF(fp,"+%d=%d",substring->nmismatches_refdiff - substring->nmismatches_bothdiff,substring->nmismatches_refdiff);
    if (print_snplabels_p && substring->nmismatches_refdiff > substring->nmismatches_bothdiff) {
      print_snp_labels(fp,substring,queryseq);
    }
  }

  if (type1 == SPLICE_JUNCTION && type2 == SPLICE_JUNCTION) {
    if (invertp == false) {
      if (Junction_sensedir(pre_junction) == SENSE_FORWARD) {
	FPRINTF(fp,",dir:sense");
      } else if (Junction_sensedir(pre_junction) == SENSE_ANTI) {
	FPRINTF(fp,",dir:antisense");
      } else {
	FPRINTF(fp,",dir:unknown");
      }
    } else {
      if (Junction_sensedir(pre_junction) == SENSE_FORWARD) {
	FPRINTF(fp,",dir:antisense");
      } else if (Junction_sensedir(pre_junction) == SENSE_ANTI) {
	FPRINTF(fp,",dir:sense");
      } else {
	FPRINTF(fp,",dir:unknown");
      }
    }
    splice_distance_1 = Junction_splice_distance(pre_junction);
    splice_distance_2 = Junction_splice_distance(post_junction);
    if (splice_distance_1 == 0 && splice_distance_2 == 0) {
      /* Skip */
    } else if (splice_distance_1 == 0) {
      FPRINTF(fp,",splice_type:consistent");
      FPRINTF(fp,",splice_dist_2:%u",splice_distance_2);
    } else if (splice_distance_2 == 0) {
      FPRINTF(fp,",splice_type:consistent");
      FPRINTF(fp,",splice_dist_1:%u",splice_distance_1);
    } else {
      FPRINTF(fp,",splice_type:consistent");
      FPRINTF(fp,",splice_dist_1:%u",splice_distance_1);
      FPRINTF(fp,",splice_dist_2:%u",splice_distance_2);
    }

  } else if (type1 == SPLICE_JUNCTION) {
    if (invertp == false) {
      if (Junction_sensedir(pre_junction) == SENSE_FORWARD) {
	FPRINTF(fp,",dir:sense");
      } else if (Junction_sensedir(pre_junction) == SENSE_ANTI) {
	FPRINTF(fp,",dir:antisense");
      } else {
	FPRINTF(fp,",dir:unknown");
      }
    } else {
      if (Junction_sensedir(pre_junction) == SENSE_FORWARD) {
	FPRINTF(fp,",dir:antisense");
      } else if (Junction_sensedir(pre_junction) == SENSE_ANTI) {
	FPRINTF(fp,",dir:sense");
      } else {
	FPRINTF(fp,",dir:unknown");
      }
    }
    if ((splice_distance_1 = Junction_splice_distance(pre_junction)) > 0) {
      FPRINTF(fp,",splice_type:consistent");
      FPRINTF(fp,",splice_dist_1:%u",splice_distance_1);
    }

  } else if (type2 == SPLICE_JUNCTION) {
    if (invertp == false) {
      if (Junction_sensedir(post_junction) == SENSE_FORWARD) {
	FPRINTF(fp,",dir:sense");
      } else if (Junction_sensedir(post_junction) == SENSE_ANTI) {
	FPRINTF(fp,",dir:antisense");
      } else {
	FPRINTF(fp,",dir:unknown");
      }
    } else {
      if (Junction_sensedir(post_junction) == SENSE_FORWARD) {
	FPRINTF(fp,",dir:antisense");
      } else if (Junction_sensedir(post_junction) == SENSE_ANTI) {
	FPRINTF(fp,",dir:sense");
      } else {
	FPRINTF(fp,",dir:unknown");
      }
    }
    if ((splice_distance_2 = Junction_splice_distance(post_junction)) > 0) {
      FPRINTF(fp,",splice_type:consistent");
      FPRINTF(fp,",splice_dist_2:%u",splice_distance_2);
    }
  }

  return;
}


#if 0
void
Substring_print_single (Filestring_T fp, T substring, Shortread_T queryseq,
			char *chr, bool invertp) {

  print_genomic(fp,substring,/*deletion*/(char *) NULL,/*deletionlength*/0,invertp,queryseq);
  FPRINTF(fp,"\t");
  print_coordinates(fp,substring,chr,invertp);

  FPRINTF(fp,"\t");
  if (invertp == false) {
    switch (substring->start_endtype) {
    case END: FPRINTF(fp,"start:%d",substring->trim_left); break;
    case TERM: FPRINTF(fp,"term:%d",substring->trim_left); break;
    default: fprintf(stderr,"start_endtype is %d\n",substring->start_endtype); abort(); break;
    }
  } else {
    switch (substring->end_endtype) {
    case END: FPRINTF(fp,"start:%d",substring->trim_right); break;
    case TERM: FPRINTF(fp,"term:%d",substring->trim_right); break;
    default: fprintf(stderr,"end_endtype is %d\n",substring->end_endtype); abort(); break;
    }
  }

  FPRINTF(fp,"..");

  if (invertp == false) {
    switch (substring->end_endtype) {
    case END: FPRINTF(fp,"end:%d",substring->trim_right); break;
    case TERM: FPRINTF(fp,"term:%d",substring->trim_right); break;
    default: fprintf(stderr,"end_endtype is %d\n",substring->end_endtype); abort(); break;
    }
  } else {
    switch (substring->start_endtype) {
    case END: FPRINTF(fp,"end:%d",substring->trim_left); break;
    case TERM: FPRINTF(fp,"term:%d",substring->trim_left); break;
    default: fprintf(stderr,"start_endtype is %d\n",substring->start_endtype); abort(); break;
    }
  }

  FPRINTF(fp,",matches:%d,sub:%d",substring->nmatches,substring->nmismatches_bothdiff);
  if (print_nsnpdiffs_p) {
    FPRINTF(fp,"+%d=%d",substring->nmismatches_refdiff - substring->nmismatches_bothdiff,substring->nmismatches_refdiff);
    if (print_snplabels_p && substring->nmismatches_refdiff > substring->nmismatches_bothdiff) {
      print_snp_labels(fp,substring,queryseq);
    }
  }

  return;
}
#endif


#if 0
void
Substring_print_insertion_1 (Filestring_T fp, T substring1, T substring2, int nindels, 
			     Shortread_T queryseq, char *chr, bool invertp) {
  T substring;

  if (invertp == false) {
    substring = substring1;
    print_genomic(fp,substring1,/*deletion*/NULL,/*deletionlength*/0,/*invertp*/false,
		  queryseq);

  } else {
    substring = substring2;
    print_genomic(fp,substring2,/*deletion*/NULL,/*deletionlength*/0,/*invertp*/true,
		  queryseq);
  }

  FPRINTF(fp,"\t");

  print_coordinates(fp,substring,chr,invertp);


  FPRINTF(fp,"\t");
  if (invertp == false) {
    FPRINTF(fp,"start:%d..ins:%d,matches:%d,sub:%d",
	    substring->trim_left,nindels,substring->nmatches,substring->nmismatches_bothdiff);
  } else {
    FPRINTF(fp,"start:%d..ins:%d,matches:%d,sub:%d",
	    substring->trim_right,nindels,substring->nmatches,substring->nmismatches_bothdiff);
  }
  if (print_nsnpdiffs_p) {
    FPRINTF(fp,"+%d=%d",substring->nmismatches_refdiff - substring->nmismatches_bothdiff,substring->nmismatches_refdiff);
    if (print_snplabels_p && substring->nmismatches_refdiff > substring->nmismatches_bothdiff) {
      print_snp_labels(fp,substring,queryseq);
    }
  }


  return;
}
#endif

#if 0
void
Substring_print_insertion_2 (Filestring_T fp, T substring1, T substring2, int nindels,
			     Shortread_T queryseq, char *chr, bool invertp) {
  T substring;

  if (invertp == false) {
    substring = substring2;
    print_genomic(fp,substring2,/*deletion*/NULL,/*deletionlength*/0,/*invertp*/false,
		  queryseq);

  } else {
    substring = substring1;
    print_genomic(fp,substring1,/*deletion*/NULL,/*deletionlength*/0,/*invertp*/true,
		  queryseq);
  }

  FPRINTF(fp,"\t");

  print_coordinates(fp,substring,chr,invertp);


  FPRINTF(fp,"\t");
  if (invertp == false) {
    FPRINTF(fp,"ins:%d..end:%d,matches:%d,sub:%d",
	    nindels,substring->trim_right,substring->nmatches,substring->nmismatches_bothdiff);
  } else {
    FPRINTF(fp,"ins:%d..end:%d,matches:%d,sub:%d",
	    nindels,substring->trim_left,substring->nmatches,substring->nmismatches_bothdiff);
  }
  if (print_nsnpdiffs_p) {
    FPRINTF(fp,"+%d=%d",substring->nmismatches_refdiff - substring->nmismatches_bothdiff,substring->nmismatches_refdiff);
    if (print_snplabels_p && substring->nmismatches_refdiff > substring->nmismatches_bothdiff) {
      print_snp_labels(fp,substring,queryseq);
    }
  }

  return;
}
#endif


#if 0
void
Substring_print_deletion_1 (Filestring_T fp, T substring1, T substring2, int nindels, 
			    char *deletion, Shortread_T queryseq, char *chr, 
			    bool invertp) {
  T substring;

  if (invertp == false) {
    substring = substring1;
    print_genomic(fp,substring1,deletion,nindels,/*invertp*/false,queryseq);
  } else {
    substring = substring2;
    print_genomic(fp,substring2,deletion,nindels,/*invertp*/true,queryseq);
  }

  FPRINTF(fp,"\t");

  print_coordinates(fp,substring,chr,invertp);


  FPRINTF(fp,"\t");
  if (invertp == false) {
    FPRINTF(fp,"start:%d..del:%d,matches:%d,sub:%d",
	    substring->trim_left,nindels,substring->nmatches,substring->nmismatches_bothdiff);
  } else {
    FPRINTF(fp,"start:%d..del:%d,matches:%d,sub:%d",
	    substring->trim_right,nindels,substring->nmatches,substring->nmismatches_bothdiff);
  }
  if (print_nsnpdiffs_p) {
    FPRINTF(fp,"+%d=%d",substring->nmismatches_refdiff - substring->nmismatches_bothdiff,substring->nmismatches_refdiff);
    if (print_snplabels_p && substring->nmismatches_refdiff > substring->nmismatches_bothdiff) {
      print_snp_labels(fp,substring,queryseq);
    }
  }

  return;
}
#endif


#if 0
void
Substring_print_deletion_2 (Filestring_T fp, T substring1, T substring2, int nindels, 
			    Shortread_T queryseq, char *chr, bool invertp) {
  T substring;

  if (invertp == false) {
    substring = substring2;
    print_genomic(fp,substring2,/*deletion*/NULL,/*deletionlength*/0,/*invertp*/false,
		  queryseq);

  } else {
    substring = substring1;
    print_genomic(fp,substring1,/*deletion*/NULL,/*deletionlength*/0,/*invertp*/true,
		  queryseq);
  }

  FPRINTF(fp,"\t");

  print_coordinates(fp,substring,chr,invertp);

  FPRINTF(fp,"\t");
  if (invertp == false) {
    FPRINTF(fp,"del:%d..end:%d,matches:%d,sub:%d",
	    nindels,substring->trim_right,substring->nmatches,substring->nmismatches_bothdiff);
  } else {
    FPRINTF(fp,"del:%d..end:%d,matches:%d,sub:%d",
	    nindels,substring->trim_left,substring->nmatches,substring->nmismatches_bothdiff);
  }
  if (print_nsnpdiffs_p) {
    FPRINTF(fp,"+%d=%d",substring->nmismatches_refdiff - substring->nmismatches_bothdiff,substring->nmismatches_refdiff);
    if (print_snplabels_p && substring->nmismatches_refdiff > substring->nmismatches_bothdiff) {
      print_snp_labels(fp,substring,queryseq);
    }
  }
  
  return;
}
#endif


/* This logic used in splice part of SAM_print */
static void
print_splice_distance (Filestring_T fp, T donor, T acceptor, Chrpos_T distance, int sensedir, char *tag) {
  bool normalp = true;

  if (donor == NULL || acceptor == NULL) {
    /* Don't print anything */
  } else if (distance == 0U) {
    FPRINTF(fp,",splice_type:%s",TRANSLOCATION_TEXT);
  } else {
    if (donor->plusp != acceptor->plusp) {
      FPRINTF(fp,",splice_type:%s",INVERSION_TEXT);
      normalp = false;
    } else if (donor->plusp == true) {
      if (sensedir != SENSE_ANTI) {
	/* SENSE_FORWARD or SENSE_NULL */
	if (acceptor->genomicstart < donor->genomicstart) {
	  FPRINTF(fp,",splice_type:%s",SCRAMBLE_TEXT);
	  normalp = false;
	}
      } else {
	if (donor->genomicstart < acceptor->genomicstart) {
	  FPRINTF(fp,",splice_type:%s",SCRAMBLE_TEXT);
	  normalp = false;
	}
      }
    } else {
      if (sensedir != SENSE_ANTI) {
	/* SENSE_FORWARD or SENSE_NULL */
	if (donor->genomicstart < acceptor->genomicstart) {
	  FPRINTF(fp,",splice_type:%s",SCRAMBLE_TEXT);
	  normalp = false;
	}
      } else {
	if (acceptor->genomicstart < donor->genomicstart) {
	  FPRINTF(fp,",splice_type:%s",SCRAMBLE_TEXT);
	  normalp = false;
	}
      }
    }
    if (normalp == true) {
      FPRINTF(fp,",splice_type:%s",CONSISTENT_TEXT);
    }
    FPRINTF(fp,",%s:%u",tag,distance);
  }

  return;
}

static void
print_shortexon_splice_distances (Filestring_T fp, Chrpos_T distance1, Chrpos_T distance2) {
  if (distance1 == 0U || distance2 == 0U) {
    /* Skip */
  } else {
    FPRINTF(fp,",splice_type:%s",CONSISTENT_TEXT);
    FPRINTF(fp,",splice_dist_1:%u",distance1);
    FPRINTF(fp,",splice_dist_2:%u",distance2);
  }

  return;
}



void
Substring_print_donor (Filestring_T fp, T donor, int sensedir, bool invertp, Shortread_T queryseq,
		       Univ_IIT_T chromosome_iit, T acceptor, Chrpos_T chimera_distance) {
  char *chr, *label_tag, *splice_dist_tag;
  bool allocp;

#ifdef CHECK_KNOWNI
  Chrpos_T splicesitepos;
  int *splicesites;
  int nsplicesites;
#endif

  print_genomic(fp,donor,/*deletion*/NULL,/*deletionlength*/0,invertp,queryseq);

  FPRINTF(fp,"\t");
  chr = Univ_IIT_label(chromosome_iit,donor->chrnum,&allocp);
  print_coordinates(fp,donor,chr,invertp);

  /* printf("donor chimera_pos is %d\n",donor->chimera_pos); */
  FPRINTF(fp,"\t");
  if (sensedir == SENSE_FORWARD) {
    if (invertp == false) {
      FPRINTF(fp,"start:%d..donor:%.2f",donor->trim_left,donor->siteD_prob);
      label_tag = "label_2";
      splice_dist_tag = "splice_dist_2";
    } else {
      FPRINTF(fp,"donor:%.2f..end:%d",donor->siteD_prob,donor->trim_left);
      label_tag = "label_1";
      splice_dist_tag = "splice_dist_1";
    }
  } else if (sensedir == SENSE_ANTI) {
    if (invertp == false) {
      FPRINTF(fp,"donor:%.2f..end:%d",donor->siteD_prob,donor->trim_right);
      label_tag = "label_1";
      splice_dist_tag = "splice_dist_1";
    } else {
      FPRINTF(fp,"start:%d..donor:%.2f",donor->trim_right,donor->siteD_prob);
      label_tag = "label_2";
      splice_dist_tag = "splice_dist_2";
    }
  } else {
    /* SENSE_NULL */
    if (invertp == false) {
      FPRINTF(fp,"start:%d..splice:%.2f",donor->trim_left,donor->siteD_prob);
      label_tag = "label_2";
      splice_dist_tag = "splice_dist_2";
    } else {
      FPRINTF(fp,"splice:%.2f..end:%d",donor->siteD_prob,donor->trim_left);
      label_tag = "label_1";
      splice_dist_tag = "splice_dist_1";
    }
  }

  FPRINTF(fp,",matches:%d,sub:%d",donor->nmatches,donor->nmismatches_bothdiff);
  if (print_nsnpdiffs_p) {
    FPRINTF(fp,"+%d=%d",donor->nmismatches_refdiff - donor->nmismatches_bothdiff,donor->nmismatches_refdiff);
    if (print_snplabels_p && donor->nmismatches_refdiff > donor->nmismatches_bothdiff) {
      print_snp_labels(fp,donor,queryseq);
    }
  }

  if (sensedir == SENSE_FORWARD) {
    if (invertp == false) {
      FPRINTF(fp,",dir:sense");
    } else {
      FPRINTF(fp,",dir:antisense");
    }
  } else if (sensedir == SENSE_ANTI) {
    if (invertp == false) {
      FPRINTF(fp,",dir:antisense");
    } else {
      FPRINTF(fp,",dir:sense");
    }
  } else {
    /* SENSE_NULL */
    FPRINTF(fp,",dir:unknown");
  }

  if (acceptor != NULL) {
    print_splice_distance(fp,donor,acceptor,chimera_distance,sensedir,splice_dist_tag);
  }

#ifdef CHECK_KNOWNI
  if (donor->siteD_knownp == false && splicesites_iit) {
    if (donor->plusp == true) {
      splicesitepos = donor->genomicstart - donor->chroffset + donor->chimera_pos;
    } else {
      splicesitepos = donor->genomicstart - donor->chroffset - donor->chimera_pos;
    }
    splicesites = IIT_get_exact_multiple_with_divno(&nsplicesites,splicesites_iit,
						    splicesites_divint_crosstable[donor->chrnum],
						    splicesitepos,splicesitepos+1U,donor_typeint);
    assert(nsplicesites == 0);
  }
#endif

  if (donor->siteD_knownp && splicesites_iit) {
    print_splicesite_labels(fp,donor,donor_typeint,donor->siteD_pos,label_tag);
  }

  if (allocp == true) {
    FREE(chr);
  }

  return;
}

void 
Substring_print_acceptor (Filestring_T fp, T acceptor, int sensedir, bool invertp, Shortread_T queryseq, 
			  Univ_IIT_T chromosome_iit, T donor, Chrpos_T chimera_distance) {
  char *chr, *label_tag, *splice_dist_tag;
  bool allocp;

#ifdef CHECK_KNOWNI
  Chrpos_T splicesitepos;
  int *splicesites;
  int nsplicesites;
#endif

  print_genomic(fp,acceptor,/*deletion*/NULL,/*deletionlength*/0,invertp,queryseq);

  FPRINTF(fp,"\t");
  chr = Univ_IIT_label(chromosome_iit,acceptor->chrnum,&allocp);
  print_coordinates(fp,acceptor,chr,invertp);

  /* printf("acceptor chimera_pos is %d\n",acceptor->chimera_pos); */
  FPRINTF(fp,"\t");
  if (sensedir == SENSE_FORWARD) {
    if (invertp == false) {
      FPRINTF(fp,"acceptor:%.2f..end:%d",acceptor->siteA_prob,acceptor->trim_right);
      label_tag = "label_1";
      splice_dist_tag = "splice_dist_1";
    } else {
      FPRINTF(fp,"start:%d..acceptor:%.2f",acceptor->trim_right,acceptor->siteA_prob);
      label_tag = "label_2";
      splice_dist_tag = "splice_dist_2";
    }
  } else if (sensedir == SENSE_ANTI) {
    if (invertp == false) {
      FPRINTF(fp,"start:%d..acceptor:%.2f",acceptor->trim_left,acceptor->siteA_prob);
      label_tag = "label_2";
      splice_dist_tag = "splice_dist_2";
    } else {
      FPRINTF(fp,"acceptor:%.2f..end:%d",acceptor->siteA_prob,acceptor->trim_left);
      label_tag = "label_1";
      splice_dist_tag = "splice_dist_1";
    }
  } else {
    /* SENSE_NULL */
    if (invertp == false) {
      FPRINTF(fp,"splice:%.2f..end:%d",acceptor->siteA_prob,acceptor->trim_right);
      label_tag = "label_1";
      splice_dist_tag = "splice_dist_1";
    } else {
      FPRINTF(fp,"start:%d..splice:%.2f",acceptor->trim_right,acceptor->siteA_prob);
      label_tag = "label_2";
      splice_dist_tag = "splice_dist_2";
    }
  }

  FPRINTF(fp,",matches:%d,sub:%d",acceptor->nmatches,acceptor->nmismatches_bothdiff);
  if (print_nsnpdiffs_p) {
    FPRINTF(fp,"+%d=%d",acceptor->nmismatches_refdiff - acceptor->nmismatches_bothdiff,acceptor->nmismatches_refdiff);
    if (print_snplabels_p && acceptor->nmismatches_refdiff > acceptor->nmismatches_bothdiff) {
      print_snp_labels(fp,acceptor,queryseq);
    }
  }

  if (sensedir == SENSE_FORWARD) {
    if (invertp == false) {
      FPRINTF(fp,",dir:sense");
    } else {
      FPRINTF(fp,",dir:antisense");
    }
  } else if (sensedir == SENSE_ANTI) {
    if (invertp == false) {
      FPRINTF(fp,",dir:antisense");
    } else {
      FPRINTF(fp,",dir:sense");
    }
  } else {
    /* SENSE_NULL */
    FPRINTF(fp,",dir:unknown");
  }

  if (donor != NULL) {
    print_splice_distance(fp,donor,acceptor,chimera_distance,sensedir,splice_dist_tag);
  }

#ifdef CHECK_KNOWNI
  if (acceptor->chimera_knownp == false && splicesites_iit) {
    if (acceptor->plusp == true) {
      splicesitepos = acceptor->genomicstart - acceptor->chroffset + acceptor->siteA_pos;
    } else {
      splicesitepos = acceptor->genomicstart - acceptor->chroffset - acceptor->siteA_pos;
    }
    splicesites = IIT_get_exact_multiple_with_divno(&nsplicesites,splicesites_iit,
						    splicesites_divint_crosstable[acceptor->chrnum],
						    splicesitepos,splicesitepos+1U,acceptor_typeint);
    assert(nsplicesites == 0);
  }
#endif


  if (acceptor->siteA_knownp && splicesites_iit) {
    print_splicesite_labels(fp,acceptor,acceptor_typeint,acceptor->siteA_pos,label_tag);
  }


  if (allocp == true) {
    FREE(chr);
  }

  return;
}


void
Substring_print_shortexon (Filestring_T fp, T shortexon, int sensedir, bool invertp, Shortread_T queryseq,
			   Univ_IIT_T chromosome_iit, Chrpos_T distance1, Chrpos_T distance2) {
  char *chr;
  bool allocp;

#ifdef CHECK_KNOWNI
  Chrpos_T splicesitepos;
  int *splicesites;
  int nsplicesites;
#endif

#ifdef CHECK_KNOWNI
    if (shortexon->chimera_knownp == false && splicesites_iit) {
      if (shortexon->plusp == true) {
	splicesitepos = shortexon->genomicstart - shortexon->chroffset + shortexon->chimera_pos;
      } else {
	splicesitepos = shortexon->genomicstart - shortexon->chroffset - shortexon->chimera_pos;
      }
      splicesites = IIT_get_exact_multiple_with_divno(&nsplicesites,splicesites_iit,
						      splicesites_divint_crosstable[shortexon->chrnum],
						      splicesitepos,splicesitepos+1U,acceptor_typeint);
      assert(nsplicesites == 0);
    }

    if (shortexon->chimera_knownp_2 == false && splicesites_iit) {
      if (shortexon->plusp == true) {
	splicesitepos = shortexon->genomicstart - shortexon->chroffset + shortexon->chimera_pos_2;
      } else {
	splicesitepos = shortexon->genomicstart - shortexon->chroffset - shortexon->chimera_pos_2;
      }
      splicesites = IIT_get_exact_multiple_with_divno(&nsplicesites,splicesites_iit,
						      splicesites_divint_crosstable[shortexon->chrnum],
						      splicesitepos,splicesitepos+1U,donor_typeint);
      assert(nsplicesites == 0);
    }
#endif

  print_genomic(fp,shortexon,/*deletion*/NULL,/*deletionlength*/0,invertp,queryseq);

  FPRINTF(fp,"\t");

  chr = Univ_IIT_label(chromosome_iit,shortexon->chrnum,&allocp);
  print_coordinates(fp,shortexon,chr,invertp);

  FPRINTF(fp,"\t");
  if (sensedir == SENSE_FORWARD && invertp == false) {
    FPRINTF(fp,"acceptor:%.2f..donor:%.2f",shortexon->siteA_prob,shortexon->siteD_prob);
  } else if (sensedir == SENSE_FORWARD && invertp == true) {
    FPRINTF(fp,"donor:%.2f..acceptor:%.2f",shortexon->siteD_prob,shortexon->siteA_prob);
  } else if (sensedir == SENSE_ANTI && invertp == false) {
    FPRINTF(fp,"donor:%.2f..acceptor:%.2f",shortexon->siteD_prob,shortexon->siteA_prob);
  } else if (sensedir == SENSE_ANTI && invertp == true) {
    FPRINTF(fp,"acceptor:%.2f..donor:%.2f",shortexon->siteA_prob,shortexon->siteD_prob);
  }

  FPRINTF(fp,",matches:%d,sub:%d",shortexon->nmatches,shortexon->nmismatches_bothdiff);
  if (print_nsnpdiffs_p) {
    FPRINTF(fp,"+%d=%d",shortexon->nmismatches_refdiff - shortexon->nmismatches_bothdiff,shortexon->nmismatches_refdiff);
    if (print_snplabels_p && shortexon->nmismatches_refdiff > shortexon->nmismatches_bothdiff) {
      print_snp_labels(fp,shortexon,queryseq);
    }
  }

  if (sensedir == SENSE_FORWARD && invertp == false) {
    FPRINTF(fp,",dir:sense");
    print_shortexon_splice_distances(fp,distance1,distance2);

    if (shortexon->siteA_knownp && splicesites_iit) {
      print_splicesite_labels(fp,shortexon,acceptor_typeint,
			      shortexon->siteA_pos,/*tag*/"label_1");
    }
    if (shortexon->siteD_knownp && splicesites_iit) {
      print_splicesite_labels(fp,shortexon,donor_typeint,
			      shortexon->siteD_pos,/*tag*/"label_2");
    }

  } else if (sensedir == SENSE_FORWARD && invertp == true) {
    FPRINTF(fp,",dir:antisense");
    print_shortexon_splice_distances(fp,distance1,distance2);

    if (shortexon->siteD_knownp && splicesites_iit) {
      print_splicesite_labels(fp,shortexon,donor_typeint,
			      shortexon->siteD_pos,/*tag*/"label_1");
    }
    if (shortexon->siteA_knownp && splicesites_iit) {
      print_splicesite_labels(fp,shortexon,acceptor_typeint,
			      shortexon->siteA_pos,/*tag*/"label_2");
    }

  } else if (sensedir == SENSE_ANTI && invertp == false) {
    FPRINTF(fp,",dir:antisense");
    print_shortexon_splice_distances(fp,distance1,distance2);



    if (shortexon->siteD_knownp && splicesites_iit) {
      print_splicesite_labels(fp,shortexon,donor_typeint,
			      shortexon->siteD_pos,/*tag*/"label_1");
    }

    if (shortexon->siteA_knownp && splicesites_iit) {
      print_splicesite_labels(fp,shortexon,acceptor_typeint,
			      shortexon->siteA_pos,/*tag*/"label_2");
    }

  } else if (sensedir == SENSE_ANTI && invertp == true) {
    FPRINTF(fp,",dir:sense");
    print_shortexon_splice_distances(fp,distance1,distance2);
    if (shortexon->siteA_knownp && splicesites_iit) {
      print_splicesite_labels(fp,shortexon,acceptor_typeint,
			      shortexon->siteA_pos,/*tag*/"label_1");
    }
    if (shortexon->siteD_knownp && splicesites_iit) {
      print_splicesite_labels(fp,shortexon,donor_typeint,
			      shortexon->siteD_pos,/*tag*/"label_2");
    }
  }

  if (allocp == true) {
    FREE(chr);
  }

  return;
}


/* Needs to be here to access splicesites_iit */
void
Substring_print_gmap (Filestring_T fp, struct Pair_T *pairs, int npairs, int nsegments, bool invertedp,
		      Endtype_T start_endtype, Endtype_T end_endtype,
		      Chrnum_T chrnum, Univcoord_T chroffset, Univcoord_T chrhigh,
		      int querylength, bool watsonp, int cdna_direction, int score,
		      int insertlength, int pairscore, int mapq_score, Univ_IIT_T chromosome_iit,
		      bool pairedp, GMAP_source_T gmap_source) {
  Pair_print_gsnap(fp,pairs,npairs,nsegments,invertedp,
		   start_endtype,end_endtype,chrnum,chroffset,chrhigh,querylength,
		   watsonp,cdna_direction,score,insertlength,pairscore,mapq_score,
		   chromosome_iit,splicesites_iit,splicesites_divint_crosstable,
		   donor_typeint,acceptor_typeint,pairedp,gmap_source);
  return;
}


/************************************************************************
 *   Resolve against known genes
 ************************************************************************/

Overlap_T
Substring_gene_overlap (T this, bool favor_multiexon_p) {
  Chrpos_T start, end;

  if (genes_iit == NULL) {
    return false;
  } else if (this->plusp == true) {
    start = this->alignstart - this->chroffset + 1;
    end = this->alignend - this->chroffset;
    debug7(printf("Checking if #%d:%u..%u overlaps an exon in a known gene\n",this->chrnum,start,end));
    return IIT_gene_overlap(genes_iit,genes_divint_crosstable[this->chrnum],start,end,favor_multiexon_p);

  } else {
    start = this->alignend - this->chroffset + 1;
    end = this->alignstart - this->chroffset;
    debug7(printf("Checking if #%d:%u..%u overlaps an exon in a known gene\n",this->chrnum,start,end));
    return IIT_gene_overlap(genes_iit,genes_divint_crosstable[this->chrnum],start,end,favor_multiexon_p);
  }
}



/************************************************************************
 *   Multiclean
 ************************************************************************/

static char *
get_total_tally (long int *tally, char *ptr) {
  int n;
  char *end;

  if ((end = index(ptr,'\n')) == NULL) {
    fprintf(stderr,"Premature end of line %s\n",ptr);
    return 0;
  }
  /* fprintf(stderr,"Getting tally for %.*s\n",end-ptr,ptr); */

  while (ptr < end) {
    while (ptr < end && !isdigit((int) *ptr)) {
      ptr++;
    }
    if (ptr < end) {
      sscanf(ptr,"%d",&n);
#if 0
      debug(if (n > 0) printf(" %d",n));
#endif
      (*tally) += n;
      while (ptr < end && !isspace(*ptr)) {
	ptr++;
      }
      while (ptr < end && isspace(*ptr)) {
	ptr++;
      }
    }
  }

  return ptr;
}

long int
Substring_tally (T this, IIT_T tally_iit, int *tally_divint_crosstable) {
  long int total = 0U;
  Interval_T interval;
  char *annotation, *restofheader, *ptr;
#if 0
  bool alloc_chr_p;
#endif
  bool allocp;
  Chrpos_T chrpos, intervalend;

  Chrpos_T coordstart, coordend, pos5, pos3;
  int *matches;
  int nmatches, i;
  
  pos5 = this->alignstart_trim - this->chroffset;
  pos3 = this->alignend_trim - this->chroffset;

  if (pos5 < pos3) {
    coordstart = pos5;
    coordend = pos3;
  } else {
    coordstart = pos3;
    coordend = pos5;
  }
  coordstart += 1U;		/* Because tally IIT is 1-based */
  debug(printf("coordstart = %u, coordend = %u\n",coordstart,coordend));

#if 0
  chr = Univ_IIT_label(chromosome_iit,this->chrnum,&alloc_chr_p);
#endif
  matches = IIT_get_with_divno(&nmatches,tally_iit,tally_divint_crosstable[this->chrnum],
			       coordstart,coordend,/*sortp*/false);

  for (i = 0; i < nmatches; i++) {
    annotation = IIT_annotation(&restofheader,tally_iit,matches[i],&allocp);

    interval = IIT_interval(tally_iit,matches[i]);
    chrpos = Interval_low(interval);
    intervalend = Interval_high(interval);

    ptr = annotation;

    while (chrpos < coordstart) {
      if ((ptr = index(ptr,'\n')) == NULL) {
	fprintf(stderr,"Premature end of tally from %u to %u\n",
		Interval_low(interval),Interval_high(interval));
	return total;
      } else {
	ptr++;
      }
      chrpos++;
    }

    while (chrpos <= intervalend && chrpos <= coordend) {
      ptr = get_total_tally(&total,ptr);
      ptr++;
      chrpos++;
    }

    if (allocp == true) {
      FREE(restofheader);
    }
  }

  FREE(matches);

#if 0
  if (alloc_chr_p) {
    FREE(chr);
  }
#endif

  debug(printf("Subtotal = %ld\n",total));
  return total;
}


bool
Substring_runlength_p (T this, IIT_T runlength_iit, int *runlength_divint_crosstable) {
  Chrpos_T coordstart, coordend, pos5, pos3;
  
  pos5 = this->alignstart_trim - this->chroffset;
  pos3 = this->alignend_trim - this->chroffset;

  if (pos5 < pos3) {
    coordstart = pos5;
    coordend = pos3;
  } else {
    coordstart = pos3;
    coordend = pos5;
  }
  coordstart += 1U;		/* Because runlength IIT is 1-based */
  debug(printf("coordstart = %u, coordend = %u\n",coordstart,coordend));

  /* chr = Univ_IIT_label(chromosome_iit,this->chrnum,&alloc_chr_p); */
  return IIT_exists_with_divno(runlength_iit,runlength_divint_crosstable[this->chrnum],
				coordstart,coordend);
}


int
Substring_count_mismatches_region (T this, int trim_left, int trim_right,
				   Compress_T query_compress_fwd, Compress_T query_compress_rev) {
  int left_bound, right_bound;

  if (this == NULL) {
    return 0;
  } else if (this->ambiguous_p == true) {
    return this->nmismatches_whole; /* Could refine to test each left in ambcoords */
  }

  left_bound = trim_left;
  right_bound = this->querylength - trim_right;

  if (this->queryend_orig < left_bound) {
    return 0;
  } else if (this->querystart_orig > right_bound) {
    return 0;
  } else {
    if (this->querystart_orig > left_bound) {
      left_bound = this->querystart_orig;
    }
    if (this->queryend_orig < right_bound) {
      right_bound = this->queryend_orig;
    }

    if (this->plusp) {
      return Genome_count_mismatches_substring(query_compress_fwd,this->left,/*pos5*/left_bound,
					       /*pos3*/right_bound,/*plusp*/true,this->genestrand);
    } else {
      return Genome_count_mismatches_substring(query_compress_rev,this->left,
					       /*pos5*/this->querylength - right_bound,
					       /*pos3*/this->querylength - left_bound,
					       /*plusp*/false,this->genestrand);
    }
  }
}


/************************************************************************
 *   Conversion to Pair_T format
 ************************************************************************/

List_T
Substring_convert_to_pairs (List_T pairs, T substring, int querylength, Shortread_T queryseq,
			    int hardclip_low, int hardclip_high, int queryseq_offset) {
  int querystart, queryend, querypos, i;
  Chrpos_T chrpos;
  char *seq1;
  char genome;

  if (substring == NULL) {
    return pairs;
  }

  debug6(printf("*** Entered Substring_convert_to_pairs with querylength %d, hardclip_low %d, hardclip_high %d, queryseq_offset %d\n",
		querylength,hardclip_low,hardclip_high,queryseq_offset));

  seq1 = Shortread_fullpointer_uc(queryseq);
  if (substring->plusp == true) {
    if (hardclip_low > substring->querystart) {
      querystart = hardclip_low;
    } else {
      querystart = substring->querystart;
    }

    if (querylength - hardclip_high < substring->queryend) {
      queryend = querylength - hardclip_high;
    } else {
      queryend = substring->queryend;
    }
    /* Pairs are all zero-based, so do not add 1 */
#if 0
    chrpos = substring->genomicstart_adj + querystart - substring->chroffset /*+ 1U*/;
#else
    chrpos = substring->genomicstart + querystart - substring->chroffset /*+ 1U*/;
#endif

    debug6(printf("plus conversion\n"));
    debug6(printf("querystart %d, queryend %d, plusp %d\n",querystart,queryend,substring->plusp));
    debug6(printf("alignstart %u, alignend %u\n",substring->alignstart_trim - substring->chroffset,
		  substring->alignend_trim - substring->chroffset));
    debug6(printf("chrpos %u\n",chrpos));

    if (substring->genomic_bothdiff == NULL) {
      /* Exact match */
      for (i = querystart, querypos = queryseq_offset + querystart; i < queryend; i++, querypos++) {
	pairs = List_push_out(pairs,(void *) Pair_new_out(querypos,/*genomepos*/chrpos++,
							  seq1[i],/*comp*/MATCH_COMP,seq1[i]));
      }
    } else if (show_refdiff_p == true) {
      for (i = querystart, querypos = queryseq_offset + querystart; i < queryend; i++, querypos++) {
	if (isupper(genome = substring->genomic_refdiff[i])) {
	  assert(seq1[i] == genome || seq1[i] == 'N');
	  pairs = List_push_out(pairs,(void *) Pair_new_out(querypos,/*genomepos*/chrpos++,
							    seq1[i],/*comp*/MATCH_COMP,genome));
	} else {
	  assert(seq1[i] != toupper(genome));
	  pairs = List_push_out(pairs,(void *) Pair_new_out(querypos,/*genomepos*/chrpos++,
							    seq1[i],/*comp*/MISMATCH_COMP,toupper(genome)));
	}
      }
    } else {
      /* printf("querystart %d, queryend %d\n",querystart,queryend); */
      /* printf("seq1   %s\n",seq1); */
      /* printf("genome %s\n",substring->genomic_bothdiff); */
      for (i = querystart, querypos = queryseq_offset + querystart; i < queryend; i++, querypos++) {
	if (isupper(genome = substring->genomic_bothdiff[i])) {
	  assert(seq1[i] == genome || seq1[i] == 'N');
	  pairs = List_push_out(pairs,(void *) Pair_new_out(querypos,/*genomepos*/chrpos++,
							    seq1[i],/*comp*/MATCH_COMP,genome));
	} else {
	  assert(seq1[i] != toupper(genome));
	  pairs = List_push_out(pairs,(void *) Pair_new_out(querypos,/*genomepos*/chrpos++,
							    seq1[i],/*comp*/MISMATCH_COMP,toupper(genome)));
	}
      }
    }

  } else {
    if (hardclip_high > substring->querystart) {
      querystart = hardclip_high;
    } else {
      querystart = substring->querystart;
    }

    if (querylength - hardclip_low < substring->queryend) {
      queryend = querylength - hardclip_low;
    } else {
      queryend = substring->queryend;
    }
    /* For minus, to get 0-based coordinates, subtract 1 */
#if 0
    chrpos = substring->genomicstart_adj - querystart - substring->chroffset - 1U;
#else
    chrpos = substring->genomicstart - querystart - substring->chroffset - 1U;
#endif

    debug6(printf("minus conversion\n"));
    debug6(printf("querystart %d, queryend %d, plusp %d\n",querystart,queryend,substring->plusp));
    debug6(printf("alignstart %u, alignend %u\n",substring->alignstart_trim - substring->chroffset,
		  substring->alignend_trim - substring->chroffset));
    debug6(printf("chrpos %u\n",chrpos));

    if (substring->genomic_bothdiff == NULL) {
      /* Exact match */
      for (i = querystart, querypos = queryseq_offset + querystart; i < queryend; i++, querypos++) {
	pairs = List_push_out(pairs,(void *) Pair_new_out(querypos,/*genomepos*/chrpos--,
							  seq1[i],/*comp*/MATCH_COMP,seq1[i]));
      }
    } else if (show_refdiff_p == true) {
      for (i = querystart, querypos = queryseq_offset + querystart; i < queryend; i++, querypos++) {
	if (isupper(genome = substring->genomic_refdiff[i])) {
	  assert(seq1[i] == genome || seq1[i] == 'N');
	  pairs = List_push_out(pairs,(void *) Pair_new_out(querypos,/*genomepos*/chrpos--,
							    seq1[i],/*comp*/MATCH_COMP,genome));
	} else {
	  assert(seq1[i] != toupper(genome));
	  pairs = List_push_out(pairs,(void *) Pair_new_out(querypos,/*genomepos*/chrpos--,
							    seq1[i],/*comp*/MISMATCH_COMP,toupper(genome)));
	}
      }
    } else {
      for (i = querystart, querypos = queryseq_offset + querystart; i < queryend; i++, querypos++) {
	if (isupper(genome = substring->genomic_bothdiff[i])) {
	  /* assert(seq1[i] == genome || seq1[i] == 'N'); */
	  pairs = List_push_out(pairs,(void *) Pair_new_out(querypos,/*genomepos*/chrpos--,
							    seq1[i],/*comp*/MATCH_COMP,genome));
	} else {
	  /* assert(seq1[i] != toupper(genome)); */
	  pairs = List_push_out(pairs,(void *) Pair_new_out(querypos,/*genomepos*/chrpos--,
							    seq1[i],/*comp*/MISMATCH_COMP,toupper(genome)));
	}
      }
    }
  }

  debug6(Pair_dump_list(pairs,true));
  return pairs;
}


List_T
Substring_add_insertion (List_T pairs, T substringA, T substringB, int querylength,
			 int insertionlength, Shortread_T queryseq,
			 int hardclip_low, int hardclip_high, int queryseq_offset) {
  int querystartA, queryendA, querystartB, queryendB, querypos, i;
  Chrpos_T chrendA;
  char *seq1;


  if (substringA->plusp == true) {
    if (hardclip_low > substringA->querystart) {
      querystartA = hardclip_low;
    } else {
      querystartA = substringA->querystart;
    }

    if (querylength - hardclip_high < substringA->queryend) {
      queryendA = querylength - hardclip_high;
    } else {
      queryendA = substringA->queryend;
    }

    if (hardclip_low > substringB->querystart) {
      querystartB = hardclip_low;
    } else {
      querystartB = substringB->querystart;
    }

    if (querylength - hardclip_high < substringB->queryend) {
      queryendB = querylength - hardclip_high;
    } else {
      queryendB = substringB->queryend;
    }

    /* Pairs are all zero-based, so do not add 1 */
#if 0
    chrendA = substringA->genomicstart_adj + queryendA - substringA->chroffset /*+ 1U*/;
#else
    chrendA = substringA->genomicstart + queryendA - substringA->chroffset /*+ 1U*/;
#endif

  } else {
    if (hardclip_high > substringA->querystart) {
      querystartA = hardclip_high;
    } else {
      querystartA = substringA->querystart;
    }

    if (querylength - hardclip_low < substringA->queryend) {
      queryendA = querylength - hardclip_low;
    } else {
      queryendA = substringA->queryend;
    }

    if (hardclip_high > substringB->querystart) {
      querystartB = hardclip_high;
    } else {
      querystartB = substringB->querystart;
    }

    if (querylength - hardclip_low < substringB->queryend) {
      queryendB = querylength - hardclip_low;
    } else {
      queryendB = substringB->queryend;
    }

    /* Pairs are all zero-based, so subtract 1 */
#if 0
    chrendA = substringA->genomicstart_adj - queryendA - substringA->chroffset - 1U;
#else
    chrendA = substringA->genomicstart - queryendA - substringA->chroffset - 1U;
#endif
  }

  if (querystartA <= queryendA && querystartB <= queryendB) {
    seq1 = Shortread_fullpointer_uc(queryseq);
    querypos = queryendA + queryseq_offset;
    i = queryendA;
    while (--insertionlength >= 0) {
      pairs = List_push_out(pairs,(void *) Pair_new_out(querypos++,/*genomepos*/chrendA,
							seq1[i++],/*comp*/INDEL_COMP,' '));
    }
  }

  return pairs;
}


List_T
Substring_add_deletion (List_T pairs, T substringA, T substringB, int querylength,
			char *deletion, int deletionlength,
			int hardclip_low, int hardclip_high, int queryseq_offset) {
  int querystartA, queryendA, querystartB, queryendB, querypos, k;
  Chrpos_T chrendA;

  if (substringA->plusp == true) {
    if (hardclip_low > substringA->querystart) {
      querystartA = hardclip_low;
    } else {
      querystartA = substringA->querystart;
    }

    if (querylength - hardclip_high < substringA->queryend) {
      queryendA = querylength - hardclip_high;
    } else {
      queryendA = substringA->queryend;
    }

    if (hardclip_low > substringB->querystart) {
      querystartB = hardclip_low;
    } else {
      querystartB = substringB->querystart;
    }

    if (querylength - hardclip_high < substringB->queryend) {
      queryendB = querylength - hardclip_high;
    } else {
      queryendB = substringB->queryend;
    }

    /* Pairs are all zero-based, so do not add 1 */
#if 0
    chrendA = substringA->genomicstart_adj + queryendA - substringA->chroffset /*+ 1U*/;
#else
    chrendA = substringA->genomicstart + queryendA - substringA->chroffset /*+ 1U*/;
#endif

    if (querystartA < queryendA && querystartB < queryendB) {
      querypos = queryendA + queryseq_offset;
      for (k = 0; k < deletionlength; k++) {
	pairs = List_push_out(pairs,(void *) Pair_new_out(querypos,/*genomepos*/chrendA++,
							  ' ',/*comp*/INDEL_COMP,deletion[k]));
      }
    }

  } else {
    if (hardclip_high > substringA->querystart) {
      querystartA = hardclip_high;
    } else {
      querystartA = substringA->querystart;
    }

    if (querylength - hardclip_low < substringA->queryend) {
      queryendA = querylength - hardclip_low;
    } else {
      queryendA = substringA->queryend;
    }

    if (hardclip_high > substringB->querystart) {
      querystartB = hardclip_high;
    } else {
      querystartB = substringB->querystart;
    }

    if (querylength - hardclip_low < substringB->queryend) {
      queryendB = querylength - hardclip_low;
    } else {
      queryendB = substringB->queryend;
    }

    /* Pairs are all zero-based, so subtract 1 */
#if 0
    chrendA = substringA->genomicstart_adj - queryendA - substringA->chroffset - 1U;
#else
    chrendA = substringA->genomicstart - queryendA - substringA->chroffset - 1U;
#endif

    if (querystartA <= queryendA && querystartB <= queryendB) {
      querypos = queryendA + queryseq_offset;
      for (k = 0; k < deletionlength; k++) {
	pairs = List_push_out(pairs,(void *) Pair_new_out(querypos,/*genomepos*/chrendA--,
							  ' ',/*comp*/INDEL_COMP,deletion[k]));
      }
    }
  }

  return pairs;
}



List_T
Substring_add_intron (List_T pairs, T substringA, T substringB, int querylength,
		      int hardclip_low, int hardclip_high, int queryseq_offset) {
  int querystartA, queryendA, querystartB, queryendB, querypos;
  Chrpos_T chrendA;

  if (substringA->plusp == true) {
    if (hardclip_low > substringA->querystart) {
      querystartA = hardclip_low;
    } else {
      querystartA = substringA->querystart;
    }

    if (querylength - hardclip_high < substringA->queryend) {
      queryendA = querylength - hardclip_high;
    } else {
      queryendA = substringA->queryend;
    }

    if (hardclip_low > substringB->querystart) {
      querystartB = hardclip_low;
    } else {
      querystartB = substringB->querystart;
    }

    if (querylength - hardclip_high < substringB->queryend) {
      queryendB = querylength - hardclip_high;
    } else {
      queryendB = substringB->queryend;
    }

    /* Pairs are all zero-based, so do not add 1 */
#if 0
    chrendA = substringA->genomicstart_adj + queryendA - substringA->chroffset /*+ 1U*/;
#else
    chrendA = substringA->genomicstart + queryendA - substringA->chroffset /*+ 1U*/;
#endif

  } else {
    if (hardclip_high > substringA->querystart) {
      querystartA = hardclip_high;
    } else {
      querystartA = substringA->querystart;
    }

    if (querylength - hardclip_low < substringA->queryend) {
      queryendA = querylength - hardclip_low;
    } else {
      queryendA = substringA->queryend;
    }

    if (hardclip_high > substringB->querystart) {
      querystartB = hardclip_high;
    } else {
      querystartB = substringB->querystart;
    }

    if (querylength - hardclip_low < substringB->queryend) {
      queryendB = querylength - hardclip_low;
    } else {
      queryendB = substringB->queryend;
    }

    /* Pairs are all zero-based, so subtract 1 */
#if 0
    chrendA = substringA->genomicstart_adj - queryendA - substringA->chroffset - 1U;
#else
    chrendA = substringA->genomicstart - queryendA - substringA->chroffset - 1U;
#endif
  }

  if (querystartA <= queryendA && querystartB <= queryendB) {
    /* Add gapholder */
    /* All we really need for Pair_print_sam is to set gapp to be true */
    querypos = queryendA + queryseq_offset;
    pairs = List_push_out(pairs,(void *) Pair_new_out(querypos,/*genomepos*/chrendA,
						      ' ',/*comp*/FWD_CANONICAL_INTRON_COMP,' '));
  }

  return pairs;
}



