text:
/* Merge alignments into a single alignment based on their reference (RF) annotation.
 *
 * Special care is taken to be consistent with missing data '~'
 * annotation in two different contexts.
 *
 * First, from Infernal 1.1rc2 and later cmalign and cmsearch -A
 * output alignments: '~' can occur in GC RF and SS_cons
 * annotation. The '~' after a given nongap/nonmissing RF position x
 * (and before nongap/nonmissing RF position x+1) must all occur 5' of
 * all inserted positions (gap: '.') that occur between x and x+1.
 *
 * Second, from HMMER3, hmmbuild reads and writes '~' in sequences to
 * mark the ends of fragments. Any fragment must begin and end with a
 * contiguous string of 1 or more '~'. '~' characters are not allowed
 * within sequences anywhere except within these leading and trailing
 * contiguous stretches.
 *
 * esl-alimerge will accept as input and produce as output alignments
 * that are consistent with both of these conventions.
 *
 * EPN, Fri Nov 20 16:28:59 2009
 */
#include "esl_config.h"

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include "easel.h"
#include "esl_alphabet.h"
#include "esl_distance.h"
#include "esl_fileparser.h"
#include "esl_getopts.h"
#include "esl_msa.h"
#include "esl_msafile.h"
#include "esl_msafile2.h"
#include "esl_sqio.h"
#include "esl_stopwatch.h"
#include "esl_vectorops.h"

static char banner[] = "merge alignments based on their reference (RF) annotation";
static char usage1[]  = "[options] <alignment file 1> <alignment file 2>";
static char usage2[]  = "[options] --list <file listing n > 1 ali files to merge>\n\
\n\
  Input alignments must be in Stockholm or Pfam format.\n\
  Ouput format choices\n\
  --------------------\n\
  stockholm [default]\n\
  pfam\n\
  a2m\n\
  psiblast\n\
  afa";

static void read_list_file(char *listfile, char ***ret_alifile_list, int *ret_nalifile);
static int  update_maxgap_and_maxmis(ESL_MSA *msa, char *errbuf, int clen, int64_t alen, int *maxgap, int *maxmis);
static int  validate_and_copy_msa_annotation(const ESL_GETOPTS *go, int outfmt, ESL_MSA *mmsa, ESL_MSA **msaA, int clen, int nmsa, int alen_merged, int *maxgap, int *maxmis, char *errbuf);
static int  add_msa(ESL_MSA *mmsa, char *errbuf, ESL_MSA *msa_to_add, int *maxgap, int *maxmis, int clen, int alen_merged);
static int  inflate_string_with_gaps_and_missing(char *src_str, int64_t src_len, int64_t dst_len, int *ngapA, char gapchar, int *nmisA, char mischar, char **ret_dst_str);
static int  inflate_seq_with_gaps(const ESL_ALPHABET *abc, char *src_str, int64_t src_len, int64_t dst_len, int *ngapA, char gapchar, char **ret_dst_str);
static int  validate_no_nongaps_in_rf_gaps(const ESL_ALPHABET *abc, char *rf_str, char *other_str, int64_t len);
static int  determine_gap_columns_to_add(ESL_MSA *msa, int *maxgap, int *maxmis, int clen, int **ret_ngapA, int **ret_nmisA, int **ret_neitherA, char *errbuf);
static void write_pfam_msa_top(FILE *fp, ESL_MSA *msa);
static void write_pfam_msa_gc(FILE *fp, ESL_MSA *msa, int maxwidth);
static int64_t maxwidth(char **s, int n);
static int  rfchar_is_nongap_nonmissing(const ESL_ALPHABET *abc, char rfchar);

static ESL_OPTIONS options[] = {
  /* name         type          default  env   range togs reqs  incomp           help                                                             docgroup */
  { "--list",     eslARG_NONE,    FALSE, NULL, NULL, NULL,NULL, NULL,            "command-line argument is a file that lists ali files to merge",  99 },
  { "-h",         eslARG_NONE,    FALSE, NULL, NULL, NULL,NULL, NULL,            "help; show brief info on version and usage",                     1 },
  { "-o",         eslARG_OUTFILE,  NULL, NULL, NULL, NULL,NULL, NULL,            "output the final alignment to file <f>, not stdout",             1 },
  { "-v",         eslARG_NONE,    FALSE, NULL, NULL, NULL,"-o", NULL,            "print info on merge to stdout; requires -o",                     1 },
  { "--small",    eslARG_NONE,    FALSE, NULL, NULL, NULL,NULL, NULL,            "use minimal RAM (RAM usage will be independent of aln sizes)",   1 },
  { "--rfonly",   eslARG_NONE,    FALSE, NULL, NULL, NULL,NULL, NULL,            "remove all columns that are gaps in GC RF annotation",           1 },
  { "--informat", eslARG_STRING,  FALSE, NULL, NULL, NULL,NULL, NULL,            "NOT YET DISPLAYED",                                              99 },
  { "--outformat",eslARG_STRING,  FALSE, NULL, NULL, NULL,NULL, NULL,            "specify that output aln be format <s> (see choices above)",      1 },
  { "--rna",      eslARG_NONE,    FALSE, NULL, NULL, NULL,NULL,"--amino,--rna",  "alignments to merge are RNA alignments",                         1 },
  { "--dna",      eslARG_NONE,    FALSE, NULL, NULL, NULL,NULL,"--amino,--rna",  "alignments to merge are DNA alignments",                         1 },
  { "--amino",    eslARG_NONE,    FALSE, NULL, NULL, NULL,NULL,"--dna,--rna",    "alignments to merge are protein alignments",                     1 },
  { "--stall",    eslARG_NONE,    FALSE, NULL, NULL, NULL,NULL, NULL,            "arrest after start: for debugging under gdb",                    99 },
  { 0,0,0,0,0,0,0,0,0,0 },
};

int
main(int argc, char **argv)
{
  int           status;		               /* easel return code               */
  ESL_GETOPTS  *go      = NULL;	               /* application configuration       */
  ESL_ALPHABET *abc     = NULL;      	       /* biological alphabet             */
  char         *msafile1 = NULL;	       /* msa file 1 (stays NULL if --list) */
  char         *msafile2 = NULL;	       /* msa file 2 (stays NULL if --list) */
  char         *listfile = NULL;	       /* list file name (stays NULL unless --list) */
  int           infmt   = eslMSAFILE_UNKNOWN;  /* format code for input alifiles  */
  int           outfmt  = eslMSAFILE_UNKNOWN;  /* format code for output ali      */
  ESLX_MSAFILE *afp     = NULL;	               /* open alignment file, normal interface */
  ESL_MSAFILE2 *afp2    = NULL;	               /* open alignment file, small-mem interface */
  FILE         *ofp;		               /* output file (default is stdout) */
  char        **alifile_list = NULL;           /* list of alignment files to merge */
  int           nalifile;                      /* size of alifile_list             */
  int           do_stall;                      /* used to stall when debugging     */
  int           fi;                            /* counter over alignment files */
  int           ai, ai2;                       /* counters over alignments */
  int           nali_tot;                      /* number of alignments in all files */
  int          *nali_per_file = NULL;          /* [0..nalifile-1] number of alignments per file */
  int           nseq_tot;                      /* number of sequences in all alignments */
  int           nseq_cur;                      /* number of sequences in current alignment */
  int64_t       alen_cur;                      /* length of current alignment */
  int64_t      *alenA = NULL;                  /* [0..nali_tot-1] alignment length of input msas (even after
						* potentially removingeinserts (--rfonly)) */
  ESL_MSA     **msaA = NULL;                   /* [0..nali_tot-1] all msas read from all files */
  int          *maxgap = NULL;                 /* [0..cpos..cm->clen+1] max number of gap columns
						* before each consensus position in all alignments */
  int          *maxmis = NULL;                 /* [0..cpos..cm->clen+1] max number of missing data columns
						* before each consensus position in all alignments */
  int           nalloc = 0;                    /* current size of msaA */
  int           chunksize = 10;                /* size to increase nalloc by when realloc'ing */
  void         *tmp;                           /* for ESL_RALLOC() */
  int           clen;                          /* consensus length (non-gap #=GC RF length) of all alignments */
  int           cur_clen;                      /* consensus length (non-gap #=GC RF length) of current alignments */
  int           apos;                          /* alignment position */
  ESL_MSA      *mmsa = NULL;                   /* the merged alignment created by merging all alignments in msaA */
  int           alen_mmsa;                     /* number of columns in merged MSA */
  char          errbuf[eslERRBUFSIZE];         /* buffer for error messages */
  char         *tmpstr;                        /* used if -v, for printing file names */
  int         **usemeA = NULL;                 /* [0..nali_tot-1][0..alen]  used only if --rfonly enabled, for removing gap RF columns */
  ESL_STOPWATCH *w  = NULL;                    /* for timing the merge, only used if -o enabled */
  int           do_small;                      /* TRUE if --small, operate in special small memory mode, aln seq data is not stored */
  int           do_rfonly;                     /* TRUE if --rfonly, output only non-gap RF columns (remove all insert columns) */
  int          *ngapA = NULL;              /* [0..alen] number of insert gap columns to add after each alignment column when merging */
  int          *nmisA = NULL;               /* [0..alen] number of missing data ('~') gap columns to add after each alignment column when merging */
  int          *neitherA = NULL;           /* [0..apos..alen] = ngapA[apos] + nmisA[apos] */

  /* output formatting, only relevant if -v */
  char      *namedashes = NULL;                /* string of dashes, an underline */
  int        ni;                               /* counter                        */
  int        namewidth;                        /* max width of file name         */
  int        tmp_len;                          /* current width of merged aln    */

  /* variables only used in small mode (--small) */
  int           ngs_cur;                       /* number of GS lines in current alignment (only used if do_small) */
  int           gs_exists = FALSE;             /* set to TRUE if do_small and any input aln has >= 1 GS line */
  int           maxname, maxgf, maxgc, maxgr;  /* max length of seqname, GF tag, GC tag, GR tag in all input alignments */
  int           maxname_cur, maxgf_cur, maxgc_cur, maxgr_cur; /* max length of seqname, GF tag, GC tag, GR tag in current input alignment */
  int           margin;                        /* total margin length for output msa */

  /***********************************************
   * Parse command line
   ***********************************************/

  go = esl_getopts_Create(options);
  if (esl_opt_ProcessCmdline(go, argc, argv) != eslOK ||
      esl_opt_VerifyConfig(go)               != eslOK)
    {
      printf("Failed to parse command line: %s\n", go->errbuf);
      esl_usage(stdout, argv[0], usage1);
      esl_usage(stdout, argv[0], usage2);
      printf("\nTo see more help on available options, do %s -h\n\n", argv[0]);
      exit(1);
    }

  if (esl_opt_GetBoolean(go, "-h") )
    {
      esl_banner(stdout, argv[0], banner);
      esl_usage (stdout, argv[0], usage1);
      esl_usage (stdout, argv[0], usage2);
      puts("\n where options are:");
      esl_opt_DisplayHelp(stdout, go, 1, 2, 80);
      exit(0);
    }

  if(((! esl_opt_GetBoolean(go, "--list")) && (esl_opt_ArgNumber(go) != 2)) ||
     ((  esl_opt_GetBoolean(go, "--list")) && (esl_opt_ArgNumber(go) != 1)))
    {
      printf("Incorrect number of command line arguments.\n");
      esl_usage(stdout, argv[0], usage1);
      esl_usage(stdout, argv[0], usage2);
      printf("\nTo see more help on available options, do %s -h\n\n", argv[0]);
      exit(1);
    }

  if(esl_opt_GetBoolean(go, "--list")) {
    listfile = esl_opt_GetArg(go, 1);
  }
  else {
    msafile1 = esl_opt_GetArg(go, 1);
    msafile2 = esl_opt_GetArg(go, 2);
  }

  do_small  = (esl_opt_IsOn(go, "--small")) ? TRUE : FALSE;
  do_rfonly = (esl_opt_IsOn(go, "--rfonly"))  ? TRUE : FALSE;

  /* open output file */
  if (esl_opt_GetString(go, "-o") != NULL) {
    if ((ofp = fopen(esl_opt_GetString(go, "-o"), "w")) == NULL)
      esl_fatal("Failed to open -o output file %s\n", esl_opt_GetString(go, "-o"));
  } else ofp = stdout;

  if (esl_opt_IsOn(go, "--informat")) {
    infmt = eslx_msafile_EncodeFormat(esl_opt_GetString(go, "--informat"));
    if (infmt == eslMSAFILE_UNKNOWN) esl_fatal("%s is not a valid input sequence file format for --informat", esl_opt_GetString(go, "--informat"));
    if (do_small && infmt != eslMSAFILE_PFAM) esl_fatal("small memory mode requires Pfam formatted alignments");
  }
  if (esl_opt_IsOn(go, "--outformat")) {
    outfmt = eslx_msafile_EncodeFormat(esl_opt_GetString(go, "--outformat"));
    if (outfmt == eslMSAFILE_UNKNOWN) esl_fatal("%s is not a valid input sequence file format for --outformat", esl_opt_GetString(go, "--outformat"));
    if (do_small && outfmt != eslMSAFILE_PFAM) esl_fatal("we can only output Pfam formatted alignments in small memory mode");
  }
  else outfmt = eslMSAFILE_STOCKHOLM;

  if (do_small) {
    infmt  = eslMSAFILE_PFAM; /* this must be true, else we can't do small memory mode */
    outfmt = eslMSAFILE_PFAM;
  }

  do_stall = esl_opt_GetBoolean(go, "--stall"); /* a stall point for attaching gdb */
  while (do_stall);

  /* determine file names to merge */
  if(listfile != NULL) { /* read list file */
    read_list_file(listfile, &alifile_list, &nalifile);
    if(nalifile == 0) esl_fatal("Failed to read a single alignment file name from %s\n", listfile);
  }
  else { /* we're merging two alignment files from command-line */
    nalifile = 2;
    ESL_ALLOC(alifile_list, sizeof(char *) * nalifile);
    if((status = esl_strdup(msafile1, -1, &(alifile_list[0]))) != eslOK) esl_fatal("Error storing alignment file name %s, error status: %d\n", msafile1, status);
    if((status = esl_strdup(msafile2, -1, &(alifile_list[1]))) != eslOK) esl_fatal("Error storing alignment file name %s, error status: %d\n", msafile2, status);
  }

  /* create and start stopwatch */
  if(ofp != stdout) {
    w = esl_stopwatch_Create();
    esl_stopwatch_Start(w);
  }

  if(esl_opt_GetBoolean(go, "-v")) {
    /* determine the longest file name in alifile_list */
    namewidth = 9; /* length of 'file name' */
    for(fi = 0; fi < nalifile; fi++) {
      if((status = esl_FileTail(alifile_list[fi], FALSE, &tmpstr)) != eslOK) esl_fatal("Memory allocation error.");
      namewidth = ESL_MAX(namewidth, strlen(tmpstr));
      free(tmpstr);
    }

    ESL_ALLOC(namedashes, sizeof(char) * (namewidth+1));
    namedashes[namewidth] = '\0';
    for(ni = 0; ni < namewidth; ni++) namedashes[ni] = '-';
    fprintf(stdout, "# Reading %d alignment files...\n", nalifile);
    fprintf(stdout, "#\n");
    fprintf(stdout, "# %7s  %-*s  %7s  %9s  %9s  %13s  %8s\n", "",        namewidth,"",          "",        "",          "",            "",               "ncols");
    fprintf(stdout, "# %7s  %-*s  %7s  %9s  %9s  %13s  %8s\n", "file #",  namewidth,"file name", "ali #",   "#seq/ali",  "ncols/ali",   "# seq total",    "required");
    fprintf(stdout, "# %7s  %*s  %7s  %9s  %9s  %13s  %8s\n", "-------", namewidth, namedashes,  "-------", "---------", "---------",   "-------------", "--------");
  }

  /* Allocate and initialize */
  nalloc = chunksize;
  ESL_ALLOC(msaA,   sizeof(ESL_MSA *) * nalloc);
  ESL_ALLOC(alenA,  sizeof(int64_t) * nalloc);
  if(do_rfonly) ESL_ALLOC(usemeA, sizeof(int *) * nalloc);

  /* abc handling is weird. We only use alphabet to define gap characters in this miniapp.
   * msa's are actually read in text mode. Thus eslRNA suffices for anything.
   */
  if      (esl_opt_GetBoolean(go, "--amino"))   abc = esl_alphabet_Create(eslAMINO);
  else if (esl_opt_GetBoolean(go, "--dna"))     abc = esl_alphabet_Create(eslDNA);
  else if (esl_opt_GetBoolean(go, "--rna"))     abc = esl_alphabet_Create(eslRNA);
  else if (do_small)                      	abc = esl_alphabet_Create(eslRNA); /* alphabet is only used to define gap characters, so (in this miniapp) we're okay specifying RNA for any alignment (even non-RNA ones) */
  else                                          abc = esl_alphabet_Create(eslRNA); /* ditto */

  /****************************************************************
   *  Read alignments one at a time, storing them all, separately *
   ****************************************************************/

  ai = 0;
  nseq_tot = 0;
  maxname = maxgf = maxgc = maxgr = 0;
  ESL_ALLOC(nali_per_file, sizeof(int) * nalifile);
  esl_vec_ISet(nali_per_file, nalifile, 0);

  for (fi = 0; fi < nalifile; fi++) {

    if (do_small)
      {
	status  = esl_msafile2_Open(alifile_list[fi], NULL, &afp2);
	if      (status == eslENOTFOUND) esl_fatal("Alignment file %s doesn't exist or is not readable\n", alifile_list[fi]);
	else if (status != eslOK)        esl_fatal("Alignment file %s open failed with error %d\n", alifile_list[fi], status);
      }
    else
      {
	status = eslx_msafile_Open(NULL, alifile_list[fi], NULL, infmt, NULL, &afp);
	if (status != eslOK) eslx_msafile_OpenFailure(afp, status);
      }

    /* while loop: while we have an alignment in current alignment file, (statement looks weird b/c we use a different function if --small) */
    while((status = (do_small) ?
	   esl_msafile2_ReadInfoPfam(afp2, NULL, NULL, -1, NULL,NULL, &(msaA[ai]), &nseq_cur, &alen_cur, &ngs_cur, &maxname_cur, &maxgf_cur, &maxgc_cur, &maxgr_cur, NULL, NULL, NULL, NULL, NULL) :
	   eslx_msafile_Read        (afp, &(msaA[ai]))) == eslOK) {

      if(msaA[ai]->rf == NULL) esl_fatal("Error, all alignments must have #=GC RF annotation; alignment %d of file %d does not (%s)\n", nali_per_file[fi], (fi+1), alifile_list[fi]);
      msaA[ai]->abc = abc; /* msa's are read in text mode, so this is (currently) only used to define gap characters, it doesn't even have to be the correct alphabet. if --small, this is set as RNA regardless of input */

      if (do_small) {
	maxname = ESL_MAX(maxname, maxname_cur);
	maxgf   = ESL_MAX(maxgf, maxgf_cur);
	maxgc   = ESL_MAX(maxgc, maxgc_cur);
	maxgr   = ESL_MAX(maxgr, maxgr_cur);
	msaA[ai]->alen = alen_cur;
	if(ngs_cur > 0) gs_exists = TRUE;
      }
      else {
	nseq_cur = msaA[ai]->nseq;
      }
      alenA[ai] = msaA[ai]->alen; /* impt if --small and --rfonly, to remember total width of aln to expect in second pass */
      nali_per_file[fi]++;
      nseq_tot += nseq_cur;

      /* reallocate per-alignment data, if nec */
      if((ai+1) == nalloc) {
	nalloc += chunksize;
	ESL_RALLOC(msaA,   tmp, sizeof(ESL_MSA *) * nalloc);
	ESL_RALLOC(alenA,  tmp, sizeof(int64_t) * nalloc);
	for(ai2 = ai+1; ai2 < nalloc; ai2++) { msaA[ai2] = NULL; }
	if(do_rfonly) {
	  ESL_RALLOC(usemeA, tmp, sizeof(int *) * nalloc);
	  for(ai2 = ai+1; ai2 < nalloc; ai2++) { usemeA[ai2] = NULL; }
	}
      }

      /* either store consensus (non-gap RF) length (if first aln), or verify it is what we expect */
      cur_clen = 0;
      for(apos = 0; apos < (int) msaA[ai]->alen; apos++) {
	if(rfchar_is_nongap_nonmissing(msaA[ai]->abc, msaA[ai]->rf[apos])) cur_clen++;
      }
      if(ai == 0) { /* first alignment, store clen, allocate maxgap, maxmis */
	clen = cur_clen;
	ESL_ALLOC(maxgap, sizeof(int) * (clen+1));
	esl_vec_ISet(maxgap, (clen+1), 0);
	ESL_ALLOC(maxmis, sizeof(int) * (clen+1));
	esl_vec_ISet(maxmis, (clen+1), 0); /* these will all stay 0 unless we see '~' in the alignments */
      }
      else if(cur_clen != clen) {
	esl_fatal("Error, all alignments must have identical non-gap #=GC RF lengths; expected (RF length of first ali read): %d,\nalignment %d of file %d length is %d (%s))\n", clen, nali_per_file[fi], (fi+1), cur_clen, alifile_list[fi]);
      }

      if(do_rfonly) {
	/* Remove all columns that are gaps/missing data in the RF annotation, we keep an array of usemes,
	 * one per aln, in case of --small, so we know useme upon second pass of alignment files */
	ESL_ALLOC(usemeA[ai], sizeof(int) * (msaA[ai]->alen));
	for(apos = 0; apos < msaA[ai]->alen; apos++) { usemeA[ai][apos] = (rfchar_is_nongap_nonmissing(abc, msaA[ai]->rf[apos])) ? TRUE : FALSE; }
	if((status = esl_msa_ColumnSubset(msaA[ai], errbuf, usemeA[ai])) != eslOK) {
	  esl_fatal("status code: %d removing gap RF columns for msa %d from file %s:\n%s", status, (ai+1), alifile_list[fi], errbuf);
	}
      }
      else { /* --rfonly not enabled, determine max number inserts between each position */
	/* update_maxgap_and_maxmis checks to make sure the msa follows our rule for missing data and gaps:
	 * for all positions x..y b/t any two adjacent nongap RF positions (x-1 and y+1 are nongap RF positions)
	 * all missing data columns '~' must come before all gap positions ('.', '-', or '_')
	 */
      	if((status = update_maxgap_and_maxmis(msaA[ai], errbuf, clen, msaA[ai]->alen, maxgap, maxmis)) != eslOK) esl_fatal(errbuf);
      }
      if(esl_opt_GetBoolean(go, "-v")) {
	if((status = esl_FileTail(alifile_list[fi], FALSE, &tmpstr)) != eslOK) esl_fatal("Memory allocation error.");
	tmp_len = clen + esl_vec_ISum(maxgap, clen+1) + esl_vec_ISum(maxmis, clen+1);
	fprintf(stdout, "  %7d  %-*s  %7d  %9d  %9" PRId64 "  %13d  %8d\n", (fi+1),  namewidth, tmpstr, (ai+1), nseq_cur, msaA[ai]->alen, nseq_tot, tmp_len);
	free(tmpstr);
      }
      ai++;
    } /* end of while eslx_msafile_Read() loop */

    if (do_small)
      {
	if      (status == eslEFORMAT) esl_fatal("Alignment file %s, parse error:\n%s\n", alifile_list[fi], afp2->errbuf);
	else if (status == eslEINVAL)  esl_fatal("Alignment file %s, parse error:\n%s\n", alifile_list[fi], afp2->errbuf);
	else if (status != eslEOF)     esl_fatal("Alignment file %s, read failed with error code %d\n", alifile_list[fi], status);
      }
    else if (status != eslEOF) eslx_msafile_ReadFailure(afp, status);

    if(nali_per_file[fi] == 0)     esl_fatal("Failed to read any alignments from file %s\n", alifile_list[fi]);

    if (do_small) esl_msafile2_Close(afp2);
    else          eslx_msafile_Close(afp);
  } /* end of for (fi=0; fi < nalifile; fi++) */
  nali_tot = ai;

  /*********************************************
   *  Merge all alignments into the merged MSA *
   *********************************************/

  /* We allocate space for all sequences, but leave sequences as NULL (nseq = -1).
   * If (do_small) we didn't store the sequences on the first pass through the
   * alignment files, and we'll never allocate space for the sequences in mmsa,
   * we'll just output them as we reread them on another pass through the
   * individual alignments. If we read >= 1 GS line in any of the input alignments,
   * we need to do an additional pass through the files, outputting only GS
   * data. Then, in a final (3rd) pass we'll output aligned data.
   *
   * if (!do_small),  we have the sequences in memory, we'll copy these
   * to the merged alignment, freeing them in the orignal msaA alignments
   * as we go so we never need to allocate the full mmsa while we still have
   * the individual msas (in msaA[]) in memory.
   */
  mmsa = esl_msa_Create(nseq_tot, -1);
  alen_mmsa = clen + esl_vec_ISum(maxgap, clen+1) + esl_vec_ISum(maxmis, clen+1);

  /* Determine what annotation from the input alignments
   * we will include in the merged MSA.
   * See comments in header of validate_and_copy_msa_annotation()
   * for rules on what we include.
   */
  if((status = validate_and_copy_msa_annotation(go, outfmt, mmsa, msaA, nali_tot, clen, alen_mmsa, maxgap, maxmis, errbuf)) != eslOK)
    esl_fatal("Error while checking and copying individual MSA annotation to merged MSA:%s\n", errbuf);

  if (do_small) {
    /* Small memory mode, do up to 2 more passes through the input alignments:
     * pass 2 will output only GS data at top of output alignment file (only performed if >= 1 GS line read in input files
     * pass 3 will output all aligned data at to output alignment file (always performed)
     */

    /* output header, comments, and #=GF data */
    write_pfam_msa_top(ofp, mmsa);

    if(ofp != stdout) {
      if(esl_opt_GetBoolean(go, "-v")) { fprintf(stdout, "#\n"); }
      fprintf(stdout, "# Outputting merged alignment to file %s ... ", esl_opt_GetString(go, "-o"));
      fflush(stdout);
    }

    /* if there was any GS annotation in any of the individual alignments,
     * do second pass through alignment files, outputting GS annotation as we go. */
    if(gs_exists) {
      ai = 0;
      for(fi = 0; fi < nalifile; fi++) {
	/* we're in small memory mode... */
	status = esl_msafile2_Open(alifile_list[fi], NULL, &afp2); /* this should work b/c it did on the first pass */
	if      (status == eslENOTFOUND) esl_fatal("Second pass, alignment file %s doesn't exist or is not readable\n", alifile_list[fi]);
	else if (status == eslEFORMAT)   esl_fatal("Second pass, couldn't determine format of alignment %s\n", alifile_list[fi]);
	else if (status != eslOK)        esl_fatal("Second pass, alignment file %s open failed with error %d\n", alifile_list[fi], status);

	for(ai2 = 0; ai2 < nali_per_file[fi]; ai2++) {
	  status = esl_msafile2_RegurgitatePfam(afp2, ofp,
						maxname, maxgf, maxgc, maxgr, /* max width of a seq name, gf tag, gc tag, gr tag (irrelevant here) */
						FALSE, /* regurgitate stockholm header ? */
						FALSE, /* regurgitate // trailer ? */
						FALSE, /* regurgitate blank lines */
						FALSE, /* regurgitate comments */
						FALSE, /* regurgitate GF ? */
						TRUE,  /* regurgitate GS ? */
						FALSE, /* regurgitate GC ? */
						FALSE, /* regurgitate GR ? */
						FALSE, /* regurgitate aseq ? */
						NULL,  /* regurgitate all seqs, not a subset */
						NULL,  /* regurgitate all seqs, not a subset */
						NULL,  /* we're not keeping a subset of columns */
						NULL,  /* we're not adding all gap columns */
						alenA[ai], /* alignment length, as we read it in first pass (inserts may have been removed since then) */
						'.',   /* gap char, irrelevant */
						NULL,  /* don't return num seqs read */
						NULL); /* don't return num seqs regurgitated */
	  if(status == eslEOF) esl_fatal("Second pass, error out of alignments too soon, when trying to read aln %d of file %s", ai2, alifile_list[fi]);
	  if(status != eslOK)  esl_fatal("Second pass, error reading alignment %d of file %s: %s", ai2, alifile_list[fi], afp2->errbuf);
	  ai++;
	  fflush(ofp);
	}
	esl_msafile2_Close(afp2);
      }
      fprintf(ofp, "\n"); /* a single blank line to separate GS annotation from aligned data */
    }

    /* do another (either second or third) pass through alignment files, outputting aligned sequence data (and GR) as we go */
    ai = 0;
    for(fi = 0; fi < nalifile; fi++) {
      status = esl_msafile2_Open(alifile_list[fi], NULL, &afp2); /* this should work b/c it did on the first pass */
      if      (status == eslENOTFOUND) esl_fatal("Second pass, alignment file %s doesn't exist or is not readable\n", alifile_list[fi]);
      else if (status == eslEFORMAT)   esl_fatal("Second pass, couldn't determine format of alignment %s\n", alifile_list[fi]);
      else if (status != eslOK)        esl_fatal("Second pass, alignment file %s open failed with error %d\n", alifile_list[fi], status);

      for(ai2 = 0; ai2 < nali_per_file[fi]; ai2++) {
	/* determine how many all gap columns to insert after each alignment position
	 * of the child msa when copying it to the merged msa */
	if(! do_rfonly) {
	  if((status = determine_gap_columns_to_add(msaA[ai], maxgap, maxmis, clen, &(ngapA), &(nmisA), &(neitherA), errbuf)) != eslOK)
	    esl_fatal("error determining number of all gap columns to add to alignment %d of file %s", ai2, alifile_list[fi]);
	}
	status = esl_msafile2_RegurgitatePfam(afp2, ofp,
					      maxname, maxgf, maxgc, maxgr, /* max width of a seq name, gf tag, gc tag, gr tag */
					      FALSE, /* regurgitate stockholm header ? */
					      FALSE, /* regurgitate // trailer ? */
					      FALSE, /* regurgitate blank lines */
					      FALSE, /* regurgitate comments */
					      FALSE, /* regurgitate GF ? */
					      FALSE, /* regurgitate GS ? */
					      FALSE, /* regurgitate GC ? */
					      TRUE,  /* regurgitate GR ? */
					      TRUE,  /* regurgitate aseq ? */
					      NULL,  /* regurgitate all seqs, not a subset */
					      NULL,  /* regurgitate all seqs, not a subset */
					      (do_rfonly) ? usemeA[ai] : NULL,
					      (do_rfonly) ? NULL       : neitherA,
					      alenA[ai], /* alignment length, as we read it in first pass (inserts may have been removed since then) */
					      '.',
					      NULL,  /* don't return num seqs read */
					      NULL); /* don't return num seqs regurgitated */

	if(status == eslEOF) esl_fatal("Second pass, error out of alignments too soon, when trying to read aln %d of file %s", ai2, alifile_list[fi]);
	if(status != eslOK)  esl_fatal("Second pass, error reading alignment %d of file %s: %s", ai2, alifile_list[fi], afp2->errbuf);
	free(ngapA);
	free(nmisA);
	free(neitherA);
	esl_msa_Destroy(msaA[ai]);
	msaA[ai] = NULL;
	ai++;
	fflush(ofp);
      }
      esl_msafile2_Close(afp2);
    }
    /* finally, write GC annotation and '//' line */
    margin = maxname + 1;
    if (maxgc > 0 && maxgc+6 > margin) margin = maxgc+6;
    if (maxgr > 0 && maxname+maxgr+7 > margin) margin = maxname+maxgr+7;
    write_pfam_msa_gc(ofp, mmsa, margin);
  } /* end of if(do_small) */

  else { /* ! do_small: for each input alignment in msaA[], add all aligned data to mmsa, then free it  */
    if(esl_opt_GetBoolean(go, "-v")) { fprintf(stdout, "#\n# Merging alignments ... "); fflush(stdout); }
    for(ai = 0; ai < nali_tot; ai++) {
      if((status = add_msa(mmsa, errbuf, msaA[ai], maxgap, maxmis, clen, alen_mmsa)) != eslOK)
	esl_fatal("Error, merging alignment %d of %d:\n%s.", (ai+1), nali_tot, errbuf);
      esl_msa_Destroy(msaA[ai]); /* note: the aligned sequences will have already been freed in add_msa() */
      msaA[ai] = NULL;
    }
    if(esl_opt_GetBoolean(go, "-v")) { fprintf(stdout, "done.\n#\n"); fflush(stdout); }
    mmsa->alen = alen_mmsa; /* it was -1, b/c we filled in each seq as we marched through each msaA[] alignment */
    if(ofp != stdout) { fprintf(stdout, "# Saving alignment to file %s ... ", esl_opt_GetString(go, "-o")); }
    status = eslx_msafile_Write(ofp, mmsa, outfmt);
    if(status != eslOK) esl_fatal("Error, during alignment output; status code: %d\n", status);
  }
  if(ofp != stdout) {
    fflush(stdout);
    fclose(ofp);
    fprintf(stdout, "done\n#\n");
    fflush(stdout);
    esl_stopwatch_Stop(w);
    esl_stopwatch_Display(stdout, w, "# CPU time: ");
  }

  /* clean up and exit */
  if(alifile_list != NULL) {
    for(fi = 0; fi < nalifile; fi++) {
      if(alifile_list[fi] != NULL) free(alifile_list[fi]);
    }
    free(alifile_list);
  }
  if(usemeA != NULL) {
    for(ai = 0; ai < nali_tot; ai++) {
      free(usemeA[ai]);
    }
    free(usemeA);
  }

  if(nali_per_file != NULL) free(nali_per_file);
  if(alenA != NULL)         free(alenA);
  if(namedashes != NULL)    free(namedashes);
  if(msaA != NULL)          free(msaA);
  if(maxgap != NULL)        free(maxgap);
  if(maxmis != NULL)        free(maxmis);
  if(mmsa != NULL)          esl_msa_Destroy(mmsa);
  if(abc  != NULL)          esl_alphabet_Destroy(abc);
  if(w    != NULL)          esl_stopwatch_Destroy(w);
  esl_getopts_Destroy(go);
  return 0;

 ERROR:
  esl_fatal("Out of memory. Reformat to Pfam with esl-reformat and try esl-alimerge --small.");
  return eslEMEM; /*NEVERREACHED*/
}

/* Function: read_list_file
 * Date:     EPN, Fri Nov 20 16:41:32 2009
 *
 * Read a file listing alignment files to merge.
 * Store file names in *ret_alifile_list and return it,
 * return number of files in ret_nalifile and return it.
 * Each white-space delimited token is considered a
 * different alignment name.
 *
 * Dies if we encounter an error.
 *
 * Returns: void.
 */
void
read_list_file(char *listfile, char ***ret_alifile_list, int *ret_nalifile)
{
  int             status;
  ESL_FILEPARSER *efp;
  char           *tok;
  int nalloc     = 10;
  int chunksize  = 10;
  char **alifile_list = NULL;
  int n = 0;
  void *tmp;

  ESL_ALLOC(alifile_list, sizeof(char *) * nalloc);
  status = esl_fileparser_Open(listfile, NULL,  &efp);
  if     (status == eslENOTFOUND) esl_fatal("List file %s does not exist or is not readable\n", listfile);
  else if(status == eslEMEM)      esl_fatal("Ran out of memory when opening list file %s\n", listfile);
  else if(status != eslOK)        esl_fatal("Error opening list file %s\n", listfile);

  esl_fileparser_SetCommentChar(efp, '#');
  while((status = esl_fileparser_GetToken(efp, &tok, NULL)) != eslEOF) {
    if(n == nalloc) { nalloc += chunksize; ESL_RALLOC(alifile_list, tmp, sizeof(char *) * nalloc); }
    if((status = esl_strdup(tok, -1, &(alifile_list[n++]))) != eslOK) {
      esl_fatal("Error storing alignment file name while reading list file %s, error status: %d\n", listfile, status);
    }
  }
  esl_fileparser_Close(efp);
  *ret_alifile_list = alifile_list;
  *ret_nalifile = n;

  return;

 ERROR:
  esl_fatal("Out of memory.");
  return; /*NOTREACHED*/
}

/* Function: update_maxgap_and_maxmis
 * Date:     EPN, Sun Nov 22 09:40:48 2009
 *           (stolen from Infernal's cmalign.c)
 *
 * Update maxgap[] and maxmis[], arrays that keeps track of the max
 * number of gap columns and missing data columns ('~' gap #=GC RF)
 * columns before each cpos (consensus (non-gap #=GC RF) column)
 *
 * We require that all missing columns between cpos x and x+1 occur before
 * all gap columns between cpos x and x+1.
 *
 * Consensus columns are index [0..cpos..clen].
 *
 * max{gap,mis}[0]      is number of gap/missing columns before 1st cpos.
 * max{gap,mis}[clen-1] is number of gap/missing columns before final cpos.
 * max{gap,mis}[clen]   is number of gap/missing columns after  final cpos.
 *
 * Returns: eslOK on success.
 *          eslEINVAL if msa->rf is NULL, nongap RF length is not clen of
 *                    for any two cpos x and x+1, a gap column precedes a
 *                    missing data column.
 *
 */
int
update_maxgap_and_maxmis(ESL_MSA *msa, char *errbuf, int clen, int64_t alen, int *maxgap, int *maxmis)
{
  int apos;
  int cpos = 0;
  int ngap = 0;
  int nmis = 0;

  for(apos = 0; apos < alen; apos++) {
    if(esl_abc_CIsGap(msa->abc, msa->rf[apos])) {
      ngap++;
    }
    else if (esl_abc_CIsMissing(msa->abc, msa->rf[apos])) {
      nmis++;
      if(ngap > 0) ESL_FAIL(eslEINVAL, errbuf, "after nongap RF pos %d, %d gap columns precede a missing data column (none should)", cpos, ngap);
    }
    else {
      maxgap[cpos] = ESL_MAX(maxgap[cpos], ngap);
      maxmis[cpos] = ESL_MAX(maxmis[cpos], nmis);
      cpos++;
      ngap = 0;
      nmis = 0;
    }
  }

  /* update final value, max{ins,el}[clen+1], the number of inserts
   * after the final consensus position */
  maxgap[cpos] = ESL_MAX(maxgap[cpos], ngap);
  maxmis[cpos] = ESL_MAX(maxmis[cpos], nmis);
  if(cpos != clen) ESL_FAIL(eslEINVAL, errbuf, "second pass expected clen (%d) not equal to actual clen (%d).\n", clen, cpos);

  return eslOK;
}

/* Function: validate_and_copy_msa_annotation
 * Date:     EPN, Tue Nov 24 05:35:50 2009
 *
 * Decide what individual MSA annotation from
 * the input alignments in msaA[], if any, will be
 * included in the merged alignment (mmsa) and
 * add that info to it.
 *
 * Rules for what to include in mmsa:
 *
 * We include name,desc,acc,author annotation in merged alignment
 * if it is identical in all msaA[] input alignments.
 *
 * We include comments and per-file (GF) annotation if they
 * are present and identical in all input msaA[] alignments.
 *
 * We include per-column (GC) annotation if it is present and
 * identical *with-respect-to* #=GC RF annotation AND all
 * the annotation in gap  #=GC RF columns in all msaA[] are gaps
 * ('.'). This also pertains to the following parsed per-column
 * annotation: ss_cons, sa_cons, pp_cons, and rf. With rf,
 * de-gapped rf annotation must be identical in all input
 * alignments period, if it is not we'll die with an error message.
 *
 * Per-sequence information and per-residue information is always
 * included in merged alignment. This is done by add_msa() function.
 *
 * Returns: eslOK on success.
 *          eslEINCONCEIVABLE if input/msa is corrupt in some way (example: ngf>0 but gf_tag[0] == NULL)
 *          eslEMEM on memory error
 *          if !eslOK, errbuf is filled before return
 */
int
validate_and_copy_msa_annotation(const ESL_GETOPTS *go, int outfmt, ESL_MSA *mmsa, ESL_MSA **msaA, int nmsa, int clen, int alen_merged, int *maxgap, int *maxmis, char *errbuf)
{
  int status;
  int j;                    /* counter over alignment annotations */
  int j2;                   /* counter over alignment annotations */
  int ai;                   /* counter over alignments */
  char *dealigned   = NULL; /* a temporary, dealigned string */
  char *dealigned2  = NULL; /* another temporary, dealigned string */
  char *gapped_out  = NULL; /* a temporary string with gaps added to fit into merged aln */
  int *ngapA    = NULL; /* [0..alen] number of insert gap columns to add after each alignment column when merging */
  int *nmisA     = NULL; /* [0..alen] number of missing data ('~') gap columns to add after each alignment column when merging */
  int *neitherA = NULL; /* [0..apos..alen] = ngapA[apos] + nmisA[apos] */
  int do_add;
  int found_tag;
  int be_verbose = FALSE;

  /* we only print info about annotation if -v AND we'll actually
   * output it (as stockholm or pfam)
   * (we actually don't even need to be in this function if we're
   * not output in stockholm or pfam...)
   */
  if((esl_opt_GetBoolean(go, "-v")) &&
     (outfmt == eslMSAFILE_STOCKHOLM || outfmt == eslMSAFILE_PFAM))
    { be_verbose = TRUE; }

  if(be_verbose) fprintf(stdout, "#\n");

  if(nmsa == 0) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "in validate_and_copy_msa_annotation(): zero child alignments.");

  /* First, determine how many all gap columns to insert after each alignment position
   * of the first child msa, so we can (possibly) gap out GC,SS_cons,SA_cons,PP_cons annotation
   * to appropriate length when adding it to the merged MSA. */
  if((status = determine_gap_columns_to_add(msaA[0], maxgap, maxmis, clen, &(ngapA), &(nmisA), &(neitherA), errbuf)) != eslOK)
    return status;

  /* Note: esl_strcmp() can handle NULL strings (they are not identical to non-NULL strings) */

  /*********************************************************************/
  /* Check if name annotation is identical in all alignments */
  do_add = TRUE; /* until proven otherwise */
  if(msaA[0]->name != NULL) {
    for(ai = 1; ai < nmsa; ai++) {
      if(esl_strcmp(msaA[0]->name, msaA[ai]->name) != 0) { do_add = FALSE; break; }
    }
    if(do_add) {
      if(be_verbose) fprintf(stdout, "# Identical name annotation from all alignments transferred to merged alignment.\n");
      if((status = esl_strdup(msaA[0]->name, -1, &(mmsa->name))) != eslOK) goto ERROR;
    }
    else if(be_verbose) fprintf(stdout, "# Name annotation is not identical in all alignments; not included in merged alignment.\n");
  }
  else if(be_verbose) fprintf(stdout, "# Name annotation absent from (at least) first alignment; not included in merged alignment.\n");

  /*********************************************************************/
  /* Check if description annotation is identical in all alignments */
  do_add = TRUE; /* until proven otherwise */
  if(msaA[0]->desc != NULL) {
    for(ai = 1; ai < nmsa; ai++) {
      if(esl_strcmp(msaA[0]->desc, msaA[ai]->desc) != 0) { do_add = FALSE; break; }
    }
    if(do_add) {
      if(be_verbose) fprintf(stdout, "# Identical description annotation from all alignments transferred to merged alignment.\n");
      if((status = esl_strdup(msaA[0]->desc, -1, &(mmsa->desc))) != eslOK) goto ERROR;
    }
    else if(be_verbose) fprintf(stdout, "# Description annotation is not identical in all alignments; not included in merged alignment.\n");
  }
  else if(be_verbose) fprintf(stdout, "# Description annotation absent from (at least) first alignment; not included in merged alignment.\n");

  /*********************************************************************/
  /* Check if accession annotation is identical in all alignments */
  do_add = TRUE; /* until proven otherwise */
  if(msaA[0]->acc != NULL) {
    for(ai = 1; ai < nmsa; ai++) {
      if(esl_strcmp(msaA[0]->acc, msaA[ai]->acc) != 0) { do_add = FALSE; break; }
    }
    if(do_add) {
      if(be_verbose) fprintf(stdout, "# Identical accession annotation from all alignments transferred to merged alignment.\n");
      if((status = esl_strdup(msaA[0]->acc, -1, &(mmsa->acc))) != eslOK) goto ERROR;
    }
    else if(be_verbose) fprintf(stdout, "# Accession annotation is not identical in all alignments; not included in merged alignment.\n");
  }
  else if(be_verbose) fprintf(stdout, "# Accession annotation absent from (at least) first alignment; not included in merged alignment.\n");

  /*********************************************************************/
  /* Check if author annotation is identical in all alignments */
  do_add = TRUE; /* until proven otherwise */
  if(msaA[0]->au != NULL) {
    for(ai = 1; ai < nmsa; ai++) {
      if(esl_strcmp(msaA[0]->au, msaA[ai]->au) != 0) { do_add = FALSE; break; }
    }
    if(do_add) {
      if(be_verbose) fprintf(stdout, "# Identical author annotation from all alignments transferred to merged alignment.\n");
      if((status = esl_strdup(msaA[0]->au, -1, &(mmsa->au))) != eslOK) goto ERROR;
    }
    else if(be_verbose) fprintf(stdout, "# Author annotation is not identical in all alignments; not included in merged alignment.\n");
  }
  else if(be_verbose) fprintf(stdout, "# Author annotation absent from (at least) first alignment; not included in merged alignment.\n");

  /*********************************************************************/
  /* Check per-file (GF) annotation, must be present and identical in all msaA[] alignments to be included */
  if(msaA[0]->ngf > 0) {
    for(j = 0; j < msaA[0]->ngf; j++) {
      do_add = TRUE; /* until proven otherwise */
      /* verify that what we think is true is true */
      if(msaA[0]->gf_tag[j] == NULL) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "unexpectedly, GF tag %d of msaA[0] is NULL, but msaA[0]->ngf is %d.\n", j, msaA[0]->ngf);
      if(msaA[0]->gf[j]    == NULL)  ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "unexpectedly, GF annotation %d of msaA[0] is NULL, but msaA[0]->ngf is %d.\n", j, msaA[0]->ngf);
      for(ai = 1; ai < nmsa; ai++) {
	found_tag = FALSE;
	for(j2 = 0; j2 < msaA[ai]->ngf; j2++) {
	  if(esl_strcmp(msaA[0]->gf_tag[j], msaA[ai]->gf_tag[j2]) == 0) {
	    found_tag = TRUE;
	    if(esl_strcmp(msaA[0]->gf[j], msaA[ai]->gf_tag[j2]) != 0) {
	      do_add = FALSE;
	    }
	    break; /* if we found a match, do_add remains TRUE */
	  }
	}
	if(found_tag && do_add) {
	  if(be_verbose) fprintf(stdout, "# Identical GF tag %s annotation from all alignments transferred to merged alignment.\n", msaA[0]->gf_tag[j]);
	  if((status = esl_msa_AddGF(mmsa, msaA[0]->gf_tag[j], -1,  msaA[0]->gf[j], -1)) != eslOK) goto ERROR;
	}
	else {
	  if(be_verbose) fprintf(stdout, "# GF tag %s annotation from first alignment absent from >= 1 other alignments; not included in merged alignment.\n", msaA[0]->gf_tag[j]);
	}
      }
    }
  }
  else if(be_verbose) fprintf(stdout, "# Unparsed GF annotation absent from (at least) first alignment; not included in merged alignment.\n");

  /*********************************************************************/
  /* Check comments, all must be identically ordered and identical in all msaA[] aligments to include them */
  if(msaA[0]->ncomment > 0) {
    do_add = TRUE; /* until proven otherwise */
    /* make sure all alignments have same number of comments */
    for(ai = 1; ai < nmsa; ai++) {
      if(msaA[ai]->ncomment != msaA[0]->ncomment) {
	do_add = FALSE;
	break;
      }
    }
    if(do_add) {
      /* make sure all alignments have identical comments */
      for(j = 0; j < msaA[0]->ncomment; j++) {
	/* verify that what we think is true is true */
	if(msaA[0]->comment[j] == NULL) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "unexpectedly, comment %d of msaA[0] is NULL, but msaA[0]->ncomment is %d.\n", j, msaA[0]->ncomment);
	for(ai = 1; ai < nmsa; ai++) {
	  if(esl_strcmp(msaA[0]->comment[j], msaA[ai]->comment[j]) != 0) { /* comment doesn't match */
	    do_add = FALSE;
	    break;
	  }
	}
      }
    }
    if(do_add) {
      for(j = 0; j < msaA[0]->ncomment; j++) {
	if((status = esl_msa_AddComment(mmsa, msaA[0]->comment[j], -1))!= eslOK) goto ERROR;
      }
      if(be_verbose) fprintf(stdout, "# All alignments have identical comments in the same order. These were transferred to merged alignment.\n");
    }
    else {
      if(be_verbose) fprintf(stdout, "# Comments are not identical in all alignments; not included in merged alignment.\n");
    }
  }
  else if(be_verbose) fprintf(stdout, "# No comments in (at least) first alignment; not included in merged alignment.\n");

  /*********************************************************************/
  /* Check unparsed per-column (GC) annotation, it must include all gaps in gap RF columns and
   * be identical once gap RF columns are removed in all msaA[] alignments to be included. */

  if(msaA[0]->ngc > 0) {
    for(j = 0; j < msaA[0]->ngc; j++) {
      do_add = TRUE; /* until proven otherwise */
      /* verify that what we think is true is true */
      if(msaA[0]->gc_tag[j] == NULL) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "unexpectedly, GC tag %d of msaA[0] is NULL, but msaA[0]->ngf is %d.\n", j, msaA[0]->ngc);
      if(msaA[0]->gc[j]     == NULL) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "unexpectedly, GC annotation %d of msaA[0] is NULL, but msaA[0]->ngf is %d.\n", j, msaA[0]->ngc);

      /* ensure it does not have non-gaps in gap RF columns */
      if(validate_no_nongaps_in_rf_gaps(msaA[0]->abc, msaA[0]->rf, msaA[0]->gc[j], msaA[0]->alen)) { /* returns TRUE if gc[j] has 0 non-gap characters in gap columns of RF annotation */
	/* dealign gc line */
	if((status = esl_strdup(msaA[0]->gc[j], msaA[0]->alen, &(dealigned))) != eslOK) goto ERROR;
	if((status = esl_strdealign(dealigned, msaA[0]->rf, "-_.~", NULL)) != eslOK) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "unexpected error dealigning GC tag %s of msaA[0]", msaA[0]->gc_tag[j]);

	for(ai = 1; ai < nmsa; ai++) {
	  found_tag = FALSE;
	  for(j2 = 0; j2 < msaA[ai]->ngc; j2++) {
	    if(esl_strcmp(msaA[0]->gc_tag[j], msaA[ai]->gc_tag[j2]) == 0) {
	      found_tag = TRUE;
	      /* ensure it does not have non-gaps in gap RF columns */
	      if(validate_no_nongaps_in_rf_gaps(msaA[ai]->abc, msaA[ai]->rf, msaA[ai]->gc[j2], msaA[ai]->alen)) { /* returns TRUE if gc[j2] has 0 non-gap characters in gap columns of RF annotation */
		/* dealign */
		if((status = esl_strdup(msaA[ai]->gc[j2], msaA[ai]->alen, &(dealigned2))) != eslOK) goto ERROR;
		if((status = esl_strdealign(dealigned2, msaA[ai]->rf, "-_.~", NULL)) != eslOK) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "unexpected error dealigning GC tag %s of msaA[%d]", msaA[ai]->gc_tag[j2], ai);
		/* check identity */
		if(esl_strcmp(dealigned, dealigned2) != 0) { do_add = FALSE; }
		free(dealigned2);
		dealigned2 = NULL;
		break; /* if we matched, do_add remains TRUE */
	      }
	    }
	  }
	} /* end of (for(ai = 1...)) */
	if(dealigned != NULL) { free(dealigned); dealigned = NULL; }
	if(found_tag && do_add) {
	  /* gap out the the GC annotation to fit in merged alignment */
	  if((status = inflate_string_with_gaps_and_missing(msaA[0]->gc[j], msaA[0]->alen, alen_merged, neitherA, '.', NULL, '~', &(gapped_out))) != eslOK) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "Error adding gaps to create GC tag %s annotation for merged alignment.", msaA[0]->gc_tag[j]);
	  if((status = esl_msa_AppendGC(mmsa, msaA[0]->gc_tag[j], gapped_out)) != eslOK) goto ERROR;
	  free(gapped_out);
	  gapped_out = NULL;
	  if(be_verbose) fprintf(stdout, "# Identical GC tag %s annotation from all alignments transferred to merged alignment.\n", msaA[0]->gc_tag[j]);
	}
	else {
	  if(be_verbose) fprintf(stdout, "# GC tag %s annotation from first alignment absent from or different in >= 1 other alignments; not included in merged alignment.\n", msaA[0]->gc_tag[j]);
	}
      }
    } /* end of for(j = 0 j < msaA[0]->ngc... */
  } /* end of if(msaA[0]->ngc > 0) */
  else if(be_verbose) fprintf(stdout, "# Unparsed GC annotation absent from (at least) first alignment; not included in merged alignment.\n");

  /*********************************************************************/
  /* Check ss_cons: it must include all gaps in gap RF columns and be identical once gap RF columns are removed in all
   * msaA[] alignments to be included. (Same requirements as unparsed GC annotation, so code block below is analogous to one above). */
  if(msaA[0]->ss_cons != NULL) {
    do_add = TRUE; /* until proven otherwise */
    /* ensure it does not have non-gaps in gap RF columns */
    if(validate_no_nongaps_in_rf_gaps(msaA[0]->abc, msaA[0]->rf, msaA[0]->ss_cons, msaA[0]->alen)) { /* returns TRUE if ss_cons has 0 non-gap characters in gap columns of RF annotation */
      /* dealign ss_cons */
      if((status = esl_strdup(msaA[0]->ss_cons, msaA[0]->alen, &(dealigned))) != eslOK) goto ERROR;
      if((status = esl_strdealign(dealigned, msaA[0]->rf, "-_.~", NULL)) != eslOK) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "unexpected error dealigning ss_cons of msaA[0]");
      for(ai = 1; ai < nmsa; ai++) {
	if(msaA[ai]->ss_cons == NULL) {
	  do_add = FALSE;
	  break;
	}
	/* ss_cons != NULL, ensure it does not have non-gaps in gap RF columns */
	if(validate_no_nongaps_in_rf_gaps(msaA[ai]->abc, msaA[ai]->rf, msaA[ai]->ss_cons, msaA[ai]->alen)) { /* returns TRUE if ss_cons has 0 non-gap characters in gap columns of RF annotation */
	  /* dealign */
	  if((status = esl_strdup(msaA[ai]->ss_cons, msaA[ai]->alen, &(dealigned2))) != eslOK) goto ERROR;
	  if((status = esl_strdealign(dealigned2, msaA[ai]->rf, "-_.~", NULL)) != eslOK) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "unexpected error dealigning ss_cons of msaA[%d]", ai);
	  /* check identity */
	  if(esl_strcmp(dealigned, dealigned2) != 0) { do_add = FALSE; }
	  free(dealigned2);
	  dealigned2 = NULL;
	  break; /* if we matched, do_add remains TRUE */
	}
      } /* end of (for(ai = 1...)) */
      if(dealigned != NULL) { free(dealigned); dealigned = NULL; }
      if(do_add) {
	/* gap out the the ss_cons to fit in merged alignment,
	 * as a special case, we use '~' in SS_cons as we add missing data and gaps
	 * like we do in RF (for all other annotation we only add gaps) */
	if((status = inflate_string_with_gaps_and_missing(msaA[0]->ss_cons, msaA[0]->alen, alen_merged, ngapA, '.', nmisA, '~', &(gapped_out))) != eslOK) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "Error adding gaps to create SS_cons annotation for merged alignment.");
	if(mmsa->ss_cons != NULL) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "Error adding SS_cons to merged alignment, it is already non-NULL.");
	if((status = esl_strdup(gapped_out, alen_merged, &(mmsa->ss_cons))) != eslOK) goto ERROR;
	free(gapped_out);
	gapped_out = NULL;
	if(be_verbose) fprintf(stdout, "# Identical SS_cons annotation from all alignments transferred to merged alignment.\n");
      }
      else {
	if(be_verbose) fprintf(stdout, "# SS_cons annotation from first alignment absent from or different in >= 1 other alignments; not included in merged alignment.\n");
      }
    }
  } /* end of if(msaA[0]->ss_cons != NULL) */
  else if(be_verbose) fprintf(stdout, "# SS_cons annotation absent from (at least) first alignment; not included in merged alignment.\n");

  /*********************************************************************/
  /* Check sa_cons: it must include all gaps in gap RF columns and be identical once gap RF columns are removed in all
   * msaA[] alignments to be included. (Same requirements as unparsed GC annotation, so code block below is analogous to one above). */
  if(msaA[0]->sa_cons != NULL) {
    do_add = TRUE; /* until proven otherwise */
    /* ensure it does not have non-gaps in gap RF columns */
    if(validate_no_nongaps_in_rf_gaps(msaA[0]->abc, msaA[0]->rf, msaA[0]->sa_cons, msaA[0]->alen)) { /* returns TRUE if sa_cons has 0 non-gap characters in gap columns of RF annotation */
      /* dealign sa_cons */
      if((status = esl_strdup(msaA[0]->sa_cons, msaA[0]->alen, &(dealigned))) != eslOK) goto ERROR;
      if((status = esl_strdealign(dealigned, msaA[0]->rf, "-_.~", NULL)) != eslOK) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "unexpected error dealigning sa_cons of msaA[0]");
      for(ai = 1; ai < nmsa; ai++) {
	if(msaA[ai]->sa_cons == NULL) {
	  do_add = FALSE;
	  break;
	}
	/* sa_cons != NULL, ensure it does not have non-gaps in gap RF columns */
	if(validate_no_nongaps_in_rf_gaps(msaA[ai]->abc, msaA[ai]->rf, msaA[ai]->sa_cons, msaA[ai]->alen)) { /* returns TRUE if sa_cons has 0 non-gap characters in gap columns of RF annotation */
	  /* dealign */
	  if((status = esl_strdup(msaA[ai]->sa_cons, msaA[ai]->alen, &(dealigned2))) != eslOK) goto ERROR;
	  if((status = esl_strdealign(dealigned2, msaA[ai]->rf, "-_.~", NULL)) != eslOK) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "unexpected error dealigning sa_cons of msaA[%d]", ai);
	  /* check identity */
	  if(esl_strcmp(dealigned, dealigned2) != 0) { do_add = FALSE; }
	  free(dealigned2);
	  dealigned2 = NULL;
	  break; /* if we matched, do_add remains TRUE */
	}
      } /* end of (for(ai = 1...)) */
      if(dealigned != NULL) { free(dealigned); dealigned = NULL; }
      if(do_add) {
	/* gap out the the sa_cons to fit in merged alignment */
	if((status = inflate_string_with_gaps_and_missing(msaA[0]->sa_cons, msaA[0]->alen, alen_merged, neitherA, '.', NULL, '~', &(gapped_out))) != eslOK) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "Error adding gaps to create SA_cons annotation for merged alignment.");
	if(mmsa->sa_cons != NULL) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "Error adding SA_cons to merged alignment, it is already non-NULL.");
	if((status = esl_strdup(gapped_out, alen_merged, &(mmsa->sa_cons))) != eslOK) goto ERROR;
	free(gapped_out);
	gapped_out = NULL;
	if(be_verbose) fprintf(stdout, "# Identical SA_cons annotation from all alignments transferred to merged alignment.\n");
      }
      else {
	if(be_verbose) fprintf(stdout, "# SA_cons annotation from first alignment absent from or different in >= 1 other alignments; not included in merged alignment.\n");
      }
    }
  } /* end of if(msaA[0]->sa_cons != NULL) */
  else if(be_verbose) fprintf(stdout, "# SA_cons annotation absent from (at least) first alignment; not included in merged alignment.\n");

  /*********************************************************************/
  /* Check pp_cons: it must include all gaps in gap RF columns and be identical once gap RF columns are removed in all
   * msaA[] alignments to be included. (Same requirements as unparsed GC annotation, so code block below is analogous to one above). */
  if(msaA[0]->pp_cons != NULL) {
    do_add = TRUE; /* until proven otherwise */
    /* ensure it does not have non-gaps in gap RF columns */
    if(validate_no_nongaps_in_rf_gaps(msaA[0]->abc, msaA[0]->rf, msaA[0]->pp_cons, msaA[0]->alen)) { /* returns TRUE if pp_cons has 0 non-gap characters in gap columns of RF annotation */
      /* dealign pp_cons */
      if((status = esl_strdup(msaA[0]->pp_cons, msaA[0]->alen, &(dealigned))) != eslOK) goto ERROR;
      if((status = esl_strdealign(dealigned, msaA[0]->rf, "-_.~", NULL)) != eslOK) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "unexpected error dealigning pp_cons of msaA[0]");
      for(ai = 1; ai < nmsa; ai++) {
	if(msaA[ai]->pp_cons == NULL) {
	  do_add = FALSE;
	  break;
	}
	/* pp_cons != NULL, ensure it does not have non-gaps in gap RF columns */
	if(validate_no_nongaps_in_rf_gaps(msaA[ai]->abc, msaA[ai]->rf, msaA[ai]->pp_cons, msaA[ai]->alen)) { /* returns TRUE if pp_cons has 0 non-gap characters in gap columns of RF annotation */
	  /* dealign */
	  if((status = esl_strdup(msaA[ai]->pp_cons, msaA[ai]->alen, &(dealigned2))) != eslOK) goto ERROR;
	  if((status = esl_strdealign(dealigned2, msaA[ai]->rf, "-_.~", NULL)) != eslOK) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "unexpected error dealigning pp_cons of msaA[%d]", ai);
	  /* check identity */
	  if(esl_strcmp(dealigned, dealigned2) != 0) { do_add = FALSE; }
	  free(dealigned2);
	  dealigned2 = NULL;
	  break; /* if we matched, do_add remains TRUE */
	}
      } /* end of (for(ai = 1...)) */
      if(dealigned != NULL) { free(dealigned); dealigned = NULL; }
      if(do_add) {
	/* gap out the the pp_cons to fit in merged alignment */
	if((status = inflate_string_with_gaps_and_missing(msaA[0]->pp_cons, msaA[0]->alen, alen_merged, neitherA, '.', NULL, '~', &(gapped_out))) != eslOK) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "Error adding gaps to create PP_cons annotation for merged alignment.");
	if(mmsa->pp_cons != NULL) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "Error adding PP_cons to merged alignment, it is already non-NULL.");
	if((status = esl_strdup(gapped_out, alen_merged, &(mmsa->pp_cons))) != eslOK) goto ERROR;
	free(gapped_out);
	gapped_out = NULL;
	if(be_verbose) fprintf(stdout, "# Identical PP_cons annotation from all alignments transferred to merged alignment.\n");
      }
      else {
	if(be_verbose) fprintf(stdout, "# PP_cons annotation from first alignment absent from or different in >= 1 other alignments; not included in merged alignment.\n");
      }
    }
  } /* end of if(msaA[0]->pp_cons != NULL) */
  else if(be_verbose) fprintf(stdout, "# PP_cons annotation absent from (at least) first alignment; not included in merged alignment.\n");

  /*********************************************************************/
  /* Finally, validate that RF annotation is identical in all alignments after removing gaps. */

  if(msaA[0]->rf == NULL) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "All alignments must have #= GC RF annotation.");
  /* dealign rf */
  if((status = esl_strdup(msaA[0]->rf, msaA[0]->alen, &(dealigned))) != eslOK) goto ERROR;
  if((status = esl_strdealign(dealigned, dealigned, "-_.~", NULL)) != eslOK) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "unexpected error dealigning RF of msaA[0]");
  for(ai = 1; ai < nmsa; ai++) {
    if(msaA[ai]->rf == NULL) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "All alignments must have #= GC RF annotation.");
    /* dealign */
    if((status = esl_strdup(msaA[ai]->rf, msaA[ai]->alen, &(dealigned2))) != eslOK) goto ERROR;
    if((status = esl_strdealign(dealigned2, dealigned2, "-_.~", NULL)) != eslOK) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "unexpected error dealigning RF of msaA[%d]", ai);
    /* check identity */
    if(esl_strcmp(dealigned, dealigned2) != 0) {
      printf("%s\n%s\n", dealigned, dealigned2);
      ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "All alignments must have identical #=GC RF annotation, once gaps (\".\",\"-\",\"_\") are removed.\nAlignment %d de-gapped RF annotation differs from that of alignment 1.\n%s\n%s", ai+1, dealigned, dealigned2);
    }
    if(dealigned2 != NULL) { free(dealigned2); dealigned2 = NULL; }
  }
  if(dealigned  != NULL) { free(dealigned);  dealigned = NULL; }
  /* gap out the the RF to fit in merged alignment */
  if((status = inflate_string_with_gaps_and_missing(msaA[0]->rf, msaA[0]->alen, alen_merged, ngapA, '.', nmisA, '~', &(gapped_out))) != eslOK) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "Error adding gaps to create RF annotation for merged alignment.");
  if(mmsa->rf != NULL) ESL_XFAIL(eslEINCONCEIVABLE, errbuf, "Error adding RF to merged alignment, it is already non-NULL.");
  if((status = esl_strdup(gapped_out, alen_merged, &(mmsa->rf))) != eslOK) goto ERROR;
  free(gapped_out);
  gapped_out = NULL;
  if(be_verbose) fprintf(stdout, "# Identical RF annotation from all alignments transferred to merged alignment.\n");

  if(dealigned  != NULL) free(dealigned);
  if(dealigned2 != NULL) free(dealigned2);
  if(gapped_out != NULL) free(gapped_out);
  if(ngapA      != NULL) free(ngapA);
  if(nmisA      != NULL) free(nmisA);
  if(neitherA   != NULL) free(neitherA);
  return eslOK;

 ERROR:
  if(dealigned  != NULL) free(dealigned);
  if(dealigned2 != NULL) free(dealigned2);
  if(gapped_out != NULL) free(gapped_out);
  if(ngapA      != NULL) free(ngapA);
  if(nmisA      != NULL) free(nmisA);
  if(neitherA   != NULL) free(neitherA);
  return status;
}

/* Function: add_msa
 * Date:     EPN, Mon Nov 23 05:54:37 2009
 *
 * Add a "child" MSA we read from a file to the merged
 * MSA - the merged alignment that we'll eventually
 * output. We free each string in the child as soon
 * as we've added it to the merged, to save memory.
 *
 * We add all sequence data (aseq), and per sequence
 * annotation, including sqname, sqdesc, sqacc, pp, ss,
 * sa, as well as non-parsed GS and GR annotation.
 *
 * <maxgap>[0..clen] is an array specifying the
 * number of inserted columns necessary between
 * each consensus position.
 *
 * max{gap,mis}[0]      is number of gaps/missing columns before 1st cpos.
 * max{gap,mis}[clen-1] is number of gaps/missing columns before final cpos.
 * max{gap,mis}[clen]   is number of gaps/missing columns after  final cpos.
 *
 * <alen_merged> is the number of columns in the merged
 * alignment. This is the non-gap RF length plus the
 * sum of the maxgap vector.
 *
 * Returns: eslOK on success.
 *          eslEMEM on memory allocation failure.
 */
int
add_msa(ESL_MSA *mmsa, char *errbuf, ESL_MSA *msa_to_add, int *maxgap, int *maxmis, int clen, int alen_merged)
{
  int status;
  int i;                 /* counter over sequences in msa_to_add */
  int j;                 /* counter over alignment annotations */
  int mi;                /* counter over sequences in mmsa */
  void *tmp;             /* for reallocations */
  char *tmpstr;          /* used for copying GR annotation */
  int  nseq_existing;    /* number of sequences already added to mmsa, by previous calls of this function */
  int *ngapA = NULL;     /* [0..alen] number of insert gap columns to add after each alignment column when merging */
  int *nmisA = NULL;     /* [0..alen] number of missing data ('~') gap columns to add after each alignment column when merging */
  int *neitherA = NULL;  /* [0..apos..alen] = ngapA[apos] + nmisA[apos] */

  nseq_existing = mmsa->nseq;

  /* determine how many all gap columns to insert after each alignment position
   * of the child msa when copying it to the merged msa */
  if((status = determine_gap_columns_to_add(msa_to_add, maxgap, maxmis, clen, &(ngapA), &(nmisA), &(neitherA), errbuf)) != eslOK)
    return status;

  nseq_existing = mmsa->nseq;

  /* Append msa_to_add's sequence data and per-sequence annotation to mmsa after adding necessary gaps */
  /* sequence names and aligned sequence data (digitized) */
  for(i = 0, mi = nseq_existing; i < msa_to_add->nseq; i++, mi++) {
      esl_strdup(msa_to_add->sqname[i], -1, &(mmsa->sqname[mi]));

      status = inflate_seq_with_gaps(msa_to_add->abc, msa_to_add->aseq[i], msa_to_add->alen, alen_merged, neitherA, '.', &(mmsa->aseq[mi]));
      if     (status == eslEMEM)   ESL_XFAIL(status, errbuf, "Out of memory adding sequence %s", mmsa->sqname[mi]);
      else if(status != eslOK)     ESL_XFAIL(status, errbuf, "Found internal missing data symbols in seq: %s", mmsa->sqname[mi]);
      free(msa_to_add->aseq[i]); /* free immediately */
      msa_to_add->aseq[i] = NULL;
  }

  /* parsed annotation that is optional */
  /* sqacc */
  if(msa_to_add->sqacc != NULL) {
    if(mmsa->sqacc == NULL) { /* allocate for all sequences, even ones added in previous calls to add_msa() */
      ESL_ALLOC(mmsa->sqacc, sizeof(char *) * (nseq_existing + msa_to_add->nseq));
      for(mi = 0; mi < nseq_existing; mi++) { mmsa->sqacc[mi] = NULL; }
    }
    else { /* reallocate; to add space for new seqs */
      ESL_RALLOC(mmsa->sqacc, tmp, sizeof(char *) * (nseq_existing + msa_to_add->nseq));
    }
    for(i = 0, mi = nseq_existing; i < msa_to_add->nseq; i++, mi++) {
      if(msa_to_add->sqacc[i] != NULL) {
	if((status = esl_strdup(msa_to_add->sqacc[i], -1, &(mmsa->sqacc[mi]))) != eslOK)
	  ESL_XFAIL(status, errbuf, "Memory allocation error when copying sequence accession number %d.\n", i+1);
	free(msa_to_add->sqacc[i]); /* free immediately */
	msa_to_add->sqacc[i] = NULL;
      }
      else {
	mmsa->sqacc[mi] = NULL;
      }
    }
  } /* end of if(msa_to_add->sqacc != NULL */
  else if(mmsa->sqacc != NULL) {
    /* msa_to_add->sqacc == NULL, but mmsa->sqacc != NULL, reallocate mmsa->sqacc, and set new ones to NULL */
    ESL_RALLOC(mmsa->sqacc, tmp, sizeof(char *) * (nseq_existing + msa_to_add->nseq));
    for(mi = nseq_existing; mi < nseq_existing + msa_to_add->nseq; mi++) {
      mmsa->sqacc[mi] = NULL;
    }
  }

  /* sqdesc */
  if(msa_to_add->sqdesc != NULL) {
    if(mmsa->sqdesc == NULL) { /* allocate for all sequences, even ones added in previous calls to add_msa() */
      ESL_ALLOC(mmsa->sqdesc, sizeof(char *) * (nseq_existing + msa_to_add->nseq));
      for(mi = 0; mi < nseq_existing; mi++) { mmsa->sqdesc[mi] = NULL; }
    }
    else { /* reallocate; to add space for new seqs */
      ESL_RALLOC(mmsa->sqdesc, tmp, sizeof(char *) * (nseq_existing + msa_to_add->nseq));
    }
    for(i = 0, mi = nseq_existing; i < msa_to_add->nseq; i++, mi++) {
      if(msa_to_add->sqdesc[i] != NULL) {
	if((status = esl_strdup(msa_to_add->sqdesc[i], -1, &(mmsa->sqdesc[mi]))) != eslOK)
	  ESL_XFAIL(status, errbuf, "Memory allocation error when copying sequence description number %d.\n", i+1);
	free(msa_to_add->sqdesc[i]); /* free immediately */
	msa_to_add->sqdesc[i] = NULL;
      }
      else {
	mmsa->sqdesc[mi] = NULL;
      }
    }
  } /* end of if(msa_to_add->sqdesc != NULL */
  else if(mmsa->sqdesc != NULL) {
    /* msa_to_add->sqdesc == NULL, but mmsa->sqdesc != NULL, reallocate mmsa->sqdesc, and set new ones to NULL */
    ESL_RALLOC(mmsa->sqdesc, tmp, sizeof(char *) * (nseq_existing + msa_to_add->nseq));
    for(mi = nseq_existing; mi < nseq_existing + msa_to_add->nseq; mi++) {
      mmsa->sqdesc[mi] = NULL;
    }
  }

  /* per-seq posterior probabilities */
  if(msa_to_add->pp != NULL) {
    if(mmsa->pp == NULL) { /* allocate for all sequences, even ones added in previous calls to add_msa() */
      ESL_ALLOC(mmsa->pp, sizeof(char *) * (nseq_existing + msa_to_add->nseq));
      for(mi = 0; mi < nseq_existing; mi++) { mmsa->pp[mi] = NULL; }
    }
    else { /* reallocate; to add space for new seqs */
      ESL_RALLOC(mmsa->pp, tmp, sizeof(char *) * (nseq_existing + msa_to_add->nseq));
    }
    for(i = 0, mi = nseq_existing; i < msa_to_add->nseq; i++, mi++) {
      if(msa_to_add->pp[i] != NULL) {
	if((status = inflate_string_with_gaps_and_missing(msa_to_add->pp[i], msa_to_add->alen, alen_merged, neitherA, '.', NULL, '~', &(mmsa->pp[mi]))) != eslOK)
	  ESL_XFAIL(status, errbuf, "Memory allocation error when copying sequence number %d posterior probabilities.\n", i+1);
	free(msa_to_add->pp[i]); /* free immediately */
	msa_to_add->pp[i] = NULL;
      }
      else {
	mmsa->pp[mi] = NULL;
      }
    }
  } /* end of if(msa_to_add->pp != NULL */
  else if(mmsa->pp != NULL) {
    /* msa_to_add->pp == NULL, but mmsa->pp != NULL, reallocate mmsa->pp, and set new ones to NULL */
    ESL_RALLOC(mmsa->pp, tmp, sizeof(char *) * (nseq_existing + msa_to_add->nseq));
    for(mi = nseq_existing; mi < nseq_existing + msa_to_add->nseq; mi++) {
      mmsa->pp[mi] = NULL;
    }
  }

  /* per-seq secondary structure */
  if(msa_to_add->ss != NULL) {
    if(mmsa->ss == NULL) { /* allocate for all sequences, even ones added in previous calls to add_msa() */
      ESL_ALLOC(mmsa->ss, sizeof(char *) * (nseq_existing + msa_to_add->nseq));
      for(mi = 0; mi < nseq_existing; mi++) { mmsa->ss[mi] = NULL; }
    }
    else { /* reallocate; to add space for new seqs */
      ESL_RALLOC(mmsa->ss, tmp, sizeof(char *) * (nseq_existing + msa_to_add->nseq));
    }
    for(i = 0, mi = nseq_existing; i < msa_to_add->nseq; i++, mi++) {
      if(msa_to_add->ss[i] != NULL) {
	if((status = inflate_string_with_gaps_and_missing(msa_to_add->ss[i], msa_to_add->alen, alen_merged, neitherA, '.', NULL, '~', &(mmsa->ss[mi]))) != eslOK)
	  ESL_XFAIL(status, errbuf, "Memory allocation error when copying sequence number %d secondary structure.\n", i+1);
	free(msa_to_add->ss[i]); /* free immediately */
	msa_to_add->ss[i] = NULL;
      }
      else {
	mmsa->ss[mi] = NULL;
      }
    }
  } /* end of if(msa_to_add->ss != NULL */
  else if(mmsa->ss != NULL) {
    /* msa_to_add->ss == NULL, but mmsa->ss != NULL, reallocate mmsa->ss, and set new ones to NULL */
    ESL_RALLOC(mmsa->ss, tmp, sizeof(char *) * (nseq_existing + msa_to_add->nseq));
    for(mi = nseq_existing; mi < nseq_existing + msa_to_add->nseq; mi++) {
      mmsa->ss[mi] = NULL;
    }
  }

  /* per-seq surface accessibility */
  if(msa_to_add->sa != NULL) {
    if(mmsa->sa == NULL) { /* allocate for all sequences, even ones added in previous calls to add_msa() */
      ESL_ALLOC(mmsa->sa, sizeof(char *) * (nseq_existing + msa_to_add->nseq));
      for(mi = 0; mi < nseq_existing; mi++) { mmsa->sa[mi] = NULL; }
    }
    else { /* reallocate; to add space for new seqs */
      ESL_RALLOC(mmsa->sa, tmp, sizeof(char *) * (nseq_existing + msa_to_add->nseq));
    }
    for(i = 0, mi = nseq_existing; i < msa_to_add->nseq; i++, mi++) {
      if(msa_to_add->sa[i] != NULL) {
	if((status = inflate_string_with_gaps_and_missing(msa_to_add->sa[i], msa_to_add->alen, alen_merged, neitherA, '.', NULL, '~', &(mmsa->sa[mi]))) != eslOK)
	  ESL_XFAIL(status, errbuf, "Memory allocation error when copying sequence number %d surface accessibility.\n", i+1);
	free(msa_to_add->sa[i]); /* free immediately */
	msa_to_add->sa[i] = NULL;
      }
      else {
	mmsa->sa[mi] = NULL;
      }
    }
  } /* end of if(msa_to_add->sa != NULL */
  else if(mmsa->sa != NULL) {
    /* msa_to_add->sa == NULL, but mmsa->sa != NULL, reallocate mmsa->sa, and set new ones to NULL */
    ESL_RALLOC(mmsa->sa, tmp, sizeof(char *) * (nseq_existing + msa_to_add->nseq));
    for(mi = nseq_existing; mi < nseq_existing + msa_to_add->nseq; mi++) {
      mmsa->sa[mi] = NULL;
    }
  }

  /* Unparsed per-sequence (GS) annotation */
  if(msa_to_add->ngs > 0) {
    for(j = 0; j < msa_to_add->ngs; j++) {
      for(i = 0, mi = nseq_existing; i < msa_to_add->nseq; i++, mi++) {
	if(msa_to_add->gs[j][i] != NULL)
	  if((status =esl_msa_AddGS(mmsa, msa_to_add->gs_tag[j], -1, mi, msa_to_add->gs[j][i], -1)) != eslOK)
	    ESL_XFAIL(status, errbuf, "Memory allocation error when copying sequence number %d GS annotation.\n", i+1);
      }
      free(msa_to_add->gs[j][i]); /* free immediately */
      msa_to_add->gs[j][i] = NULL;
    }
  }
  /* caller will free the rest of gs via esl_msa_Destroy() */

  /* unparsed per-residue (GR) annotation */
  if(msa_to_add->gr != NULL) {
    for(j = 0; j < msa_to_add->ngr; j++) {
      for(i = 0, mi = nseq_existing; i < msa_to_add->nseq; i++, mi++) {
	if(msa_to_add->gr[j][i] != NULL) {
	  if((status = inflate_string_with_gaps_and_missing(msa_to_add->gr[j][i], msa_to_add->alen, alen_merged, neitherA, '.', NULL, '~', &(tmpstr))) != eslOK)
	    ESL_XFAIL(status, errbuf, "Memory allocation error when copying sequence number %d GR annotation.\n", i+1);
	  if((status = esl_msa_AppendGR(mmsa, msa_to_add->gr_tag[j], mi, tmpstr)) != eslOK)
	    ESL_XFAIL(status, errbuf, "Memory allocation error when copying sequence number %d GR annotation.\n", i+1);
	  free(tmpstr);
	  free(msa_to_add->gr[j][i]); /* free immediately */
	  msa_to_add->gr[j][i] = NULL;
	}
      }
    }
  }
  /* caller will free the rest of gr in esl_msa_Destroy() */

  /* msa_to_add should be destroyed by caller */

  /* update nseq in mmsa */
  mmsa->nseq += msa_to_add->nseq;

  if(ngapA    != NULL) free(ngapA);
  if(nmisA     != NULL) free(nmisA);
  if(neitherA != NULL) free(neitherA);
  return eslOK;

 ERROR:
  if(ngapA != NULL) free(ngapA);
  return status;
}

/* inflate_string_with_gaps_and_missing
 *
 * Given a string, create a new one that is a copy of it,
 * but with missing data and gaps added before each position (apos)
 * as specified by n{gap,mis}A[0..apos..len]. <gapchar> and <mischar>
 * specify the gap character and missing data character. Either
 * ngapA or nmisA can be either be NULL, if so they're treated as if
 * they're all zeroes (no gaps/missing data will be added anywhere).
 *
 * n{gap,mis}A[0]       - number of gaps/missing to add before first posn
 * n{gap,mis}A[apos]    - number of gaps/missing to add before posn apos
 * n{gap,mis}A[src_len] - number of gaps/missing to add after  final posn
 *
 * By convention, missing data symbols are always put before gap symbols
 * when more than 1 of each are to be placed after the some nongap RF position.
 *
 * ret_str is allocated here.
 *
 * Returns eslOK on success.
 *         eslEMEM on memory error.
 */
int
inflate_string_with_gaps_and_missing(char *src_str, int64_t src_len, int64_t dst_len, int *ngapA, char gapchar, int *nmisA, char mischar, char **ret_dst_str)
{
  int status;
  int src_apos = 0;
  int dst_apos  = 0;
  int i;
  char *dst_str;

  ESL_ALLOC(dst_str, sizeof(char) * (dst_len+1));
  dst_str[dst_len] = '\0';

  /* add gaps before first position */
  if(nmisA) for(i = 0; i < nmisA[0]; i++) dst_str[dst_apos++] = mischar;
  if(ngapA) for(i = 0; i < ngapA[0]; i++) dst_str[dst_apos++] = gapchar;

  /* add gaps after every position */
  for(src_apos = 0; src_apos < src_len; src_apos++) {
    dst_str[dst_apos++] = src_str[src_apos];
    if(nmisA) for(i = 0; i < nmisA[(src_apos+1)]; i++) dst_str[dst_apos++] = mischar;
    if(ngapA) for(i = 0; i < ngapA[(src_apos+1)]; i++) dst_str[dst_apos++] = gapchar;
  }

  *ret_dst_str = dst_str;
  return eslOK;

 ERROR:
  return eslEMEM;
}

/* inflate_seq_with_gaps()
 *
 * Given an aligned sequence, create a new one that is a copy of it,
 * but with gaps added before each position (apos) as specified by
 * ngapA[0..apos..len]. <gapchar> specifies the gap character.
 *
 * This function follows the HMMER3 convention of using '~' to mark
 * fragment sequences. If a sequence is a fragment, the first and final
 * string of gaps (contiguous gaps) will be marked as '~'. That is,
 * the 1st column will necessarily be a '~' and the last will be a '~'.
 * This function takes care to add '~' as gap characters in fragments,
 * where appropriate, to follow HMMER3 convention.
 *
 * This function is importantly different from inflate_string_with_gaps_and_missing()
 * in that it does not allow internal '~' missing data characters to
 * be added (which are dictated by nmisA[] in inflate_string_with_gaps_and_missing())
 * Those '~' are possible in Infernal output alignments only in
 * #=GC RF and #=GC SS_cons alignment.
 *
 * It only adds '~' as necessary to the beginning and ends of fragments.
 *
 * ngapA[0]       - number of gaps/missing to add before first posn
 * ngapA[apos]    - number of gaps/missing to add before posn apos
 * ngapA[src_len] - number of gaps/missing to add after  final posn
 *
 * ret_str is allocated here.
 *
 * Returns eslOK on success.
 *         eslEMEM on memory error.
 *         eslEINVAL if sequence (src_str) has an internal missing
 *           data symbol that violates HMMER3 convention.
 */
int
inflate_seq_with_gaps(const ESL_ALPHABET *abc, char *src_str, int64_t src_len, int64_t dst_len, int *ngapA, char gapchar, char **ret_dst_str)
{
  int status;
  int src_apos = 0;
  int dst_apos  = 0;
  int src_flpos = 0; /* position of rightmost '~' in contiguous string that begins at position 0 */
  int src_frpos = 0; /* position of leftmost  '~' in contiguous string that ends   at position src_len-1 */
  int i;
  char *dst_str;
  char char2add;
  int  i_am_fragment = FALSE;

  ESL_ALLOC(dst_str, sizeof(char) * (dst_len+1));
  dst_str[dst_len] = '\0';

  /* determine if sequence is a fragment, and determine src_flpos and src_frpos */
  src_flpos = 0;
  while(src_flpos < src_len && esl_abc_CIsMissing(abc, src_str[src_flpos])) { src_flpos++; }
  src_flpos--; /* we overshot by 1 */
  /* src_flpos is now rightmost '~' in stretch of '~' that begin at position 0, so
   * it's -1 if sequence is not a fragment (there are no leading '~' in src_str)
   */

  src_frpos = src_len-1;
  while(src_frpos > -1 && esl_abc_CIsMissing(abc, src_str[src_frpos])) { src_frpos--; }
  src_frpos++; /* we overshot by 1 */
  /* src_frpos is now leftmost '~' in stretch of '~' that end at position src_len-1, so
   * it's src_len if sequence is not a fragment (there are no trailing '~' in src_str)
   */
  i_am_fragment = (src_flpos != -1 || src_frpos != src_len) ? TRUE : FALSE;

  /* Now verify that we don't have any internal '~' characters */
  for(src_apos = src_flpos+1; src_apos <= src_frpos-1; src_apos++) {
    if(esl_abc_CIsMissing(abc, src_str[src_apos])) return eslEINVAL;
  }

  /* add gaps before first position */
  src_apos = 0;
  if(i_am_fragment) {
    char2add = ((src_flpos >= src_apos) || (src_frpos <= src_apos)) ? esl_abc_CGetMissing(abc) : gapchar;
  }
  else {
    char2add = gapchar;
  }
  if(ngapA) for(i = 0; i < ngapA[0]; i++) dst_str[dst_apos++] = char2add;

  /* add gaps or missing character after every position */
  for(src_apos = 0; src_apos < src_len; src_apos++) {
    dst_str[dst_apos++] = src_str[src_apos];
    if(i_am_fragment) {
      char2add = ((src_flpos >= src_apos) || (src_frpos <= src_apos)) ? esl_abc_CGetMissing(abc) : gapchar;
    }
    else {
      char2add = gapchar;
    }
    for(i = 0; i < ngapA[(src_apos+1)]; i++) dst_str[dst_apos++] = char2add;
  }

  *ret_dst_str = dst_str;
  return eslOK;

 ERROR:
  return eslEMEM;
}

/* validate_no_nongaps_in_rf_gaps
 *
 * Given an RF string with gaps defined as by alphabet <abc>
 * and another string of same length. Make sure none of the
 * positions that are gaps in the RF string are non-gaps in the
 * other string. Return TRUE if none are. Return FALSE if at
 * least one is.
 *
 * Returns TRUE if 0 characters in <other_str> in same position
 *         as a gap in <rf_str> are non-gaps. FALSE otherwise.
 */
int
validate_no_nongaps_in_rf_gaps(const ESL_ALPHABET *abc, char *rf_str, char *other_str, int64_t len)
{
  int64_t i;
  for(i = 0; i < len; i++) {
    if((! rfchar_is_nongap_nonmissing(abc, rf_str[i])) && (rfchar_is_nongap_nonmissing(abc, other_str[i]))) return FALSE;
  }
  return TRUE;
}

/* determine_gap_columns_to_add
 * (stolen and slightly modified from Infernal 1.1rc1's cmalign.c)
 *
 * Given <maxgap> and <maxmis>, two arrays of the number of gap RF
 * positions and '~' RF (missing data) after each non-gap RF
 * (consensus) position in the eventual final merged alignment,
 * calculate how many inserts and missing data inserts
 * we need to add at each position of <msa> to expand it out to the
 * appropriate size of the eventual merged alignment.
 *
 * max{gap,mis}[0]      is number of gaps/missing before 1st cpos in merged aln
 * max{gap,mis}[cpos]   is number of gaps/missing after cpos in merged aln
 *                             for cpos = 1..clen
 * clen is the number of non-gap RF positions in msa (and in eventual merged msa).
 *
 * We allocate fill and return ret_ngapA[0..msa->alen], ret_nmisA[0..msa->alen],
 * and ret_neitherA[0..msa->alen] here.
 *
 * ret_n{gap,mis}gapA[0]      is number of gaps/missing to add before 1st position of msa
 * ret_n{gap,mis}gapA[apos]   is number of gaps/missing to add after alignment position apos
 *                             for apos = 1..msa->alen
 *
 * ret_neitherA[apos] = ngapA[apos] + nmisA[apos]
 *
 * This is similar to the esl_msa.c helper function of the same name,
 * but that function does not bother with missing data '~'.
 *
 * Returns eslOK on success.
 *         eslEMEM on memory alloaction error
 *         eslEINVAL if gaps occur before missing data b/t any two nongap RF posns
 *         eslERANGE if a value exceeds what we expected (based on earlier
 *                   checks before this function was entered).
 *         if !eslOK, errbuf if filled.
 */
int
determine_gap_columns_to_add(ESL_MSA *msa, int *maxgap, int *maxmis, int clen, int **ret_ngapA, int **ret_nmisA, int **ret_neitherA, char *errbuf)
{
  int status;
  int apos;
  int prv_cpos = 0;  /* alignment position corresponding to consensus position cpos-1 */
  int cpos = 0;
  int ngap = 0;
  int nmis = 0;
  int *ngapA    = NULL;
  int *nmisA    = NULL;
  int *neitherA = NULL;

  /* contract check */
  if(maxmis[0] != 0) ESL_FAIL(eslEINVAL, errbuf, "missing characters exist prior to first cpos, this shouldn't happen.\n");

  ESL_ALLOC(ngapA,    sizeof(int) * (msa->alen+1));
  ESL_ALLOC(nmisA,    sizeof(int) * (msa->alen+1));
  ESL_ALLOC(neitherA, sizeof(int) * (msa->alen+1));
  esl_vec_ISet(ngapA,    (msa->alen+1), 0);
  esl_vec_ISet(nmisA,    (msa->alen+1), 0);
  esl_vec_ISet(neitherA, (msa->alen+1), 0);

  for(apos = 0; apos < msa->alen; apos++) {
    if(esl_abc_CIsMissing(msa->abc, msa->rf[apos])) {
      nmis++;
      if(ngap > 0) ESL_FAIL(eslEINVAL, errbuf, "after nongap RF pos %d, %d gap columns precede a missing data column (none should)", cpos, ngap);
    }
    else if(esl_abc_CIsGap(msa->abc, msa->rf[apos])) {
      ngap++;
    }
    else { /* a consensus position */
      /* a few sanity checks */
      if(ngap  > maxgap[cpos]) ESL_FAIL(eslEINCONCEIVABLE, errbuf, "%d inserts before cpos %d greater than max expected (%d).\n", ngap, cpos, maxgap[cpos]);
      if(nmis  > maxmis[cpos]) ESL_FAIL(eslEINCONCEIVABLE, errbuf, "%d EL inserts before cpos %d greater than max expected (%d).\n", nmis, cpos, maxmis[cpos]);

      if (cpos == 0) {
	/* gaps and/or missing data before first consensus position: flush right */
	if(maxgap[cpos] > 0 && maxmis[cpos] == 0) { /* most common case */
	  ngapA[prv_cpos]  = maxgap[cpos] - ngap; /* flush right gaps (no missing) */
	}
	else if(maxgap[cpos] == 0 && maxmis[cpos] > 0) {
	  nmisA[prv_cpos]  = maxmis[cpos] - nmis; /* flush right missing (no gaps) */
	}
	else if(maxgap[cpos] > 0 && maxmis[cpos] > 0) {
	  /* missing data (ELs) is always 5' of gaps */
	  nmisA[prv_cpos]      = maxmis[cpos] - nmis; /* flush right */
	  ngapA[prv_cpos+nmis] = maxgap[cpos] - ngap; /* flush right, after missing */
	}
      }
      else {
	/* gaps and/missing data after interior consensus position (i.e. not before 1st cpos or after last)
	 * Determine where to place inserts and/or missing data.
	 * Handle each of 4 possibilities separately, note that if
	 * maxgap[cpos] == 0 then ngap == 0, and if maxmis[cpos] == 0 then nmis == 0 (see sanity check above).
	 */
	if(maxgap[cpos] >  0 && maxmis[cpos] == 0) { /* most common case */
	  ngapA[prv_cpos + 1 + (ngap/2)] = maxgap[cpos] - ngap; /* internal cpos: split */
	}
	else if(maxgap[cpos] == 0 && maxmis[cpos] > 0) {
	  nmisA[prv_cpos + 1 + (nmis/2)] = maxmis[cpos] - nmis; /* internal cpos: split */
	}
	else if(maxgap[cpos] >  0 && maxmis[cpos] > 0) {
	  /* missing data is always 5' of gaps */
	  nmisA[prv_cpos + 1 +        (nmis/2)] = maxmis[cpos] - nmis; /* internal cpos: split */
	  ngapA[prv_cpos + 1 + nmis + (ngap/2)] = maxgap[cpos] - ngap; /* internal cpos: split */
	}
	/* final case is if (maxgap[cpos] == 0 && maxmis[cpos] == 0)
	 * in this case we do nothing.
	 */
      }
      cpos++;
      prv_cpos = apos;
      ngap = 0;
      nmis = 0;
    }
  }
  /* deal with gaps and missing data after final consensus position */

  /* first, validate that clen is what it should be */
  if(cpos != clen) {
    if(ngapA != NULL) free(ngapA);
    if(nmisA != NULL) free(nmisA);
    if(neitherA != NULL) free(neitherA);
    ESL_FAIL(eslEINCONCEIVABLE, errbuf, "consensus length (%d) is not the expected length (%d).", cpos, clen);
  }

  if(maxgap[cpos] > 0 && maxmis[cpos] == 0) { /* most common case */
    ngapA[prv_cpos + 1 + ngap] = maxgap[cpos] - ngap; /* flush left gaps (no missing) */
  }
  else if(maxgap[cpos] == 0 && maxmis[cpos] > 0) {
    nmisA[prv_cpos + 1 + nmis] = maxmis[cpos] - nmis; /* flush left missing (no gaps) */
  }
  else if(maxgap[cpos] > 0 && maxmis[cpos] > 0) {
    /* missing data (ELs) is always 5' of gaps */
    nmisA[prv_cpos + 1 + nmis]        = maxmis[cpos] - nmis; /* flush left */
    ngapA[prv_cpos + 1 + nmis + ngap] = maxgap[cpos] - ngap; /* flush left, after missing */
  }

  /* determine neitherA[], the number of gaps due to either inserts or missing data after each apos */
  for(apos = 0; apos <= msa->alen; apos++) {
    neitherA[apos] = ngapA[apos] + nmisA[apos];
  }

  *ret_ngapA  = ngapA;
  *ret_nmisA = nmisA;
  *ret_neitherA = neitherA;

  return eslOK;

 ERROR:
  if(ngapA    != NULL) free(ngapA);
  if(nmisA    != NULL) free(nmisA);
  if(neitherA != NULL) free(neitherA);
  ESL_FAIL(status, errbuf, "Memory allocation error.");
  return status; /*NEVERREACHED*/
}

/* write_pfam_msa_top
 *
 * Highly specialized function for printing out the 'top' of a
 * Stockholm alignment file, i.e. the Stockholm header line, comments
 * and GF annotation to an msa file. This function is necessary when
 * printing the merged alignment file in small memory mode (when
 * --small enabled). <msa> is an alignment with no sequence
 * information (no aseq, ax, GS, or GR data).
 */
void
write_pfam_msa_top(FILE *fp, ESL_MSA *msa)
{
  int i, maxgf;

  /* rest of code in this function was stolen verbatim from actually_write_stockholm in esl_msa.c */

  /* Magic Stockholm header
   */
  fprintf(fp, "# STOCKHOLM 1.0\n");

  /* Free text comments
   */
  for (i = 0;  i < msa->ncomment; i++)
    fprintf(fp, "# %s\n", msa->comment[i]);
  if (msa->ncomment > 0) fprintf(fp, "\n");

  maxgf = maxwidth(msa->gf_tag, msa->ngf);
  if (maxgf < 2) maxgf = 2;

  /* GF section: per-file annotation
   */
  if (msa->name != NULL) fprintf(fp, "#=GF %-*s %s\n", maxgf, "ID", msa->name);
  if (msa->acc  != NULL) fprintf(fp, "#=GF %-*s %s\n", maxgf, "AC", msa->acc);
  if (msa->desc != NULL) fprintf(fp, "#=GF %-*s %s\n", maxgf, "DE", msa->desc);
  if (msa->au   != NULL) fprintf(fp, "#=GF %-*s %s\n", maxgf, "AU", msa->au);

  /* Thresholds are hacky. Pfam has two. Rfam has one.
   */
  if      (msa->cutset[eslMSA_GA1] && msa->cutset[eslMSA_GA2])
    fprintf(fp, "#=GF %-*s %.1f %.1f\n",
	    maxgf, "GA", msa->cutoff[eslMSA_GA1], msa->cutoff[eslMSA_GA2]);
  else if (msa->cutset[eslMSA_GA1])
    fprintf(fp, "#=GF %-*s %.1f\n",
	    maxgf, "GA", msa->cutoff[eslMSA_GA1]);

  if      (msa->cutset[eslMSA_NC1] && msa->cutset[eslMSA_NC2])
    fprintf(fp, "#=GF %-*s %.1f %.1f\n",
	    maxgf, "NC", msa->cutoff[eslMSA_NC1], msa->cutoff[eslMSA_NC2]);
  else if (msa->cutset[eslMSA_NC1])
    fprintf(fp, "#=GF %-*s %.1f\n",
	    maxgf, "NC", msa->cutoff[eslMSA_NC1]);

  if      (msa->cutset[eslMSA_TC1] && msa->cutset[eslMSA_TC2])
    fprintf(fp, "#=GF %-*s %.1f %.1f\n",
	    maxgf, "TC", msa->cutoff[eslMSA_TC1], msa->cutoff[eslMSA_TC2]);
  else if (msa->cutset[eslMSA_TC1])
    fprintf(fp, "#=GF %-*s %.1f\n",
	    maxgf, "TC", msa->cutoff[eslMSA_TC1]);

  for (i = 0; i < msa->ngf; i++)
    fprintf(fp, "#=GF %-*s %s\n", maxgf, msa->gf_tag[i], msa->gf[i]);
  fprintf(fp, "\n");

  return;
}

/* write_pfam_msa_gc
 *
 * Highly specialized function for printing out the GC annotation to a
 * Pfam Stockholm alignment file (1 line/seq). This function is
 * necessary when printing the merged alignment file in small memory
 * mode (when --small enabled). <msa> is an alignment with no sequence
 * information (no aseq, ax, GS, nor GR data).
 */
void
write_pfam_msa_gc(FILE *fp, ESL_MSA *msa, int margin)
{
  int i;
  if (msa->ss_cons != NULL) { fprintf(fp, "#=GC %-*s %s\n", margin-6, "SS_cons", msa->ss_cons); }
  if (msa->sa_cons != NULL) { fprintf(fp, "#=GC %-*s %s\n", margin-6, "SA_cons", msa->sa_cons); }
  if (msa->pp_cons != NULL) { fprintf(fp, "#=GC %-*s %s\n", margin-6, "PP_cons", msa->pp_cons); }
  if (msa->rf != NULL)      { fprintf(fp, "#=GC %-*s %s\n", margin-6, "RF", msa->rf); }
  for (i = 0; i < msa->ngc; i++) { fprintf(fp, "#=GC %-*s %s\n", margin-6, msa->gc_tag[i], msa->gc[i]); }
  fprintf(fp, "//\n");
}

/* maxwidth()
 * Return the length of the longest string in
 * an array of strings.
 */
static int64_t
maxwidth(char **s, int n)
{
  int64_t max,len;
  int     i;

  max = 0;
  for (i = 0; i < n; i++)
    if (s[i] != NULL)
      {
	len = strlen(s[i]);
	if (len > max) max = len;
      }
  return max;
}

/* rfchar_is_nongap_nonmissing()
 * Return FALSE if a character from RF annotation is
 * a gap or a missing character, else return TRUE.
 */
static int
rfchar_is_nongap_nonmissing(const ESL_ALPHABET *abc, char rfchar) {
  if(esl_abc_CIsGap    (abc, rfchar)) return FALSE;
  if(esl_abc_CIsMissing(abc, rfchar)) return FALSE;
  return TRUE;
}

redpajama_set_name: RedPajamaGithub


We debated whether to participate in the craziness of Black Friday deals, but ultimately decided that we were shopping them anyway, so we might as well save our favorite finds and share them with you! Hopefully this makes your gift giving that much easier!
Today we're sharing holiday party wear. Those special pieces that adds festivity to even the most simple everyday looks. We've got 50% off pieces from Gap, Jcrew Factory, Express, Nordstrom and awesome deals from all over the rest of the web.
and for the kiddos too!
The only code you should need for these is Jcrew, enter "HOLIDAY" for 30% off at check-out!


You might be bubbling with literary excitement about a new story idea. It is very difficult to sleep in this situation.
You might feel like you have more 'tweets' left in them. Sometimes it's hard for a writer to say goodnight to Twitter.
Your excitable creative muse is refusing to go to sleep.
You can't stop thinking about how much you dislike your latest draft. Why is it that our writer demons always come out to play in the early hours?
You are questioning your entire literary career.
Your blog stats are spiking due to your latest piece of flash fiction.
Your blog stats have flatlined ever since you posted your latest piece of flash fiction.
You are waiting 'patiently' for an international writer friend, who is living in a different time zone, to answer your 'writer in literary distress' email.
Before you went to bed you submitted a piece of fiction into a competition and are busy visualising how you will react when you are informed your short story has won.
You are riddled with writer envy, after a writer acquaintance tweeted earlier to say that they have a new literary agent and this particular agent is treating them to a large KFC bucket, fries and a Sprite. All you can think about is how your acquaintance went a step too far with the KFC treat revelation.
The book you read before bedtime was so good it has given you a bad dose of writer envy.
You can't get an old unfinished story out of your head.
You can't stop thinking about a piece of negative feedback you got earlier in the day.
You are lying in bed deliberating over whether or not to kill off one of your characters.
You are trying to work out how you can kill a character.
You pinned too many motivational and inspiring quotes on Pinterest and are now charged!
You are feeling wired after a late night writing session.
It is the night before your book is launched. How does one sleep in this situation?
You can't stop thinking about an attractive and dreamy fictional character.
You are waiting for the literary elves to appear and magically finish their latest draft.
Your characters won't stay quiet.
You have been suffering from Writer's Block and are wondering whether you will ever write properly again.
After a day of editing every time you close your eyes you see your own typos and badly worded sentences.
You can't stop thinking about a catastrophic plot hole which you uncovered earlier.
They are suffering from Writer's Insomnia.
Soon many reasons us writers get insomnia!!!
Its a wonder writers sleep. I sometimes wish for a little insomnia so i can get some lone writing time!
Ha ha. I hope you're compiling a book of these lists for writers. No. 8 is number one for me.
Usually it's #3, #12, or #21 that get me.
Pingback: 28 Reasons Why A Writer Might Be Having A Sleepless Night #Writer #AmWriting — BlondeWriteMore – Keri L.
You forgot the one where they had too much caffeine and it was keeping them awake and jittery!
Oh my goodness a fav writer moment for me!


text:
package vod

//Licensed under the Apache License, Version 2.0 (the "License");
//you may not use this file except in compliance with the License.
//You may obtain a copy of the License at
//
//http://www.apache.org/licenses/LICENSE-2.0
//
//Unless required by applicable law or agreed to in writing, software
//distributed under the License is distributed on an "AS IS" BASIS,
//WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
//See the License for the specific language governing permissions and
//limitations under the License.
//
// Code generated by Alibaba Cloud SDK Code Generator.
// Changes may cause incorrect behavior and will be lost if the code is regenerated.

// ThumbnailListInGetVideoInfo is a nested struct in vod response
type ThumbnailListInGetVideoInfo struct {
	Thumbnail []Thumbnail `json:"Thumbnail" xml:"Thumbnail"`
}


VO5 Shea Cashmere condition is enriched with apple extract so it naturally leaves your hair shiny, soft and smooth. Each case contains 6 bottles. Case buy only. Save Money! Buy for a group or your extended family. Allow 14 days for shipment from warehouse. Delivery time to your destination is defined by the shipper according to the delivery option you selected. All shipments are by individual case. If you order a multiple number of cases they may not all arrive on the same day.


text:
Actress Bhumika looking stunning in the photoshoot
A post shared by Bhumika Chawla (@bhumika_chawla_t)
Bhumika Chawla once played a shining role as a top heroine in Telugu. She stopped doing acting in movies after her wedding. After her re-entry, She also made a name for herself in the film industry with her performance in the characters. She is very active on social media for her fans. She shared a photo from a recent photoshoot through her Instagram account. Bhumika Chawla in a white and pink dress looked very beautiful in this photo. Her fans commented that She is looking stunning in this photo. Bhumika is currently working on the movie Edhe Maa Katha. Already The film unit has released the first look of the film. Sumanth, Bhumika, Srikanth, Tanya Hope are all wearing helmets along the mountain roads and the look of riding the bullet bike is captivating to everyone. Gurupavan is directing the film which is being produced under the banner of Gurappa Parameshwara Productions. Music is provided by Sunil Kashyap. Currently, the shooting of this film is in full swing. The movie team is trying to complete the shooting of this movie soon and release it soon.
Ileana says no plans for the New Year
Ram's new look goes viral

redpajama_set_name: RedPajamaCommonCrawl


Students will beg to play these engaging games on any interactive whiteboard or computer. They can race the clock or each other, or they can "play" on their own. Teachers won't have to tell them that what they're really doing is practicing number bonds (aka fact families) and developing mathematical reasoning. Use these games with a Singapore Math program or with any other math curriculum, and make math come alive their students!