binaries/src/clustalo/src/mymain.c

   1 /* -*- mode: c; tab-width: 4; c-basic-offset: 4;  indent-tabs-mode: nil -*- */
   2
   3 /*********************************************************************
   4  * Clustal Omega - Multiple sequence alignment
   5  *
   6  * Copyright (C) 2010 University College Dublin
   7  *
   8  * Clustal-Omega is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU General Public License as
  10  * published by the Free Software Foundation; either version 2 of the
  11  * License, or (at your option) any later version.
  12  *
  13  * This file is part of Clustal-Omega.
  14  *
  15  ********************************************************************/
  16
  17 /*
  18  *  RCS $Id: mymain.c 255 2011-06-22 15:59:07Z fabian $
  19  */
  20
  21 #ifdef HAVE_CONFIG_H
  22 #include "config.h"
  23 #endif
  24
  25 #include <stdio.h>
  26 #include <string.h>
  27 #include <assert.h>
  28 #include <unistd.h>
  29 #include <argtable2.h>
  30 #include <ctype.h>
  31 #include <limits.h>
  32 #include <libgen.h> /* for basename only */
  33
  34 /* clustal */
  35 #include "clustal-omega.h"
  36
  37 #include "mymain.h"
  38
  39 /* hhalign (parameters) */
  40 #include "hhalign/general.h"
  41
  42 typedef struct {
  43
  44     /* Sequence input
  45      */
  46     /** sequence type (from cmdline arg) */
  47     int iSeqType;
  48     /** sequence input file. not directly used by Align() */
  49     char *pcSeqInfile;
  50     /** Dealign sequences on input. Otherwise we use the alignment
  51      * info for background-HMM creation */
  52     bool bDealignInputSeqs;
  53
  54     /* profiles: : pre-aligned sequences, whose alignment will not be changed
  55      */
  56     /** profile 1: pre-aligned sequence input. not directly used by Align() */
  57     char *pcProfile1Infile ;
  58     /** profile 2: pre-aligned sequence input. not directly used by Align() */
  59     char *pcProfile2Infile;
  60
  61     /** Input limitations
  62      */
  63     /** maximum allowed number of input sequences */
  64     int iMaxNumSeq;
  65     /** maximum allowed input sequence length */
  66     int iMaxSeqLen;
  67
  68     /* Alignment output
  69      */
  70     /** alignment output file */
  71     char *pcAlnOutfile;
  72     /** alignment output format */
  73     int iAlnOutFormat;
  74     /** force overwriting of files */
  75     bool bForceFileOverwrite;
  76
  77     /* multithreading
  78      */
  79     /** number of threads */
  80     int iThreads;
  81
  82     /* logging
  83      */
  84     char *pcLogFile;
  85
  86     opts_t aln_opts;
  87
  88     /* changes here will have to be reflected in SetDefaultUserOpts(),
  89      * FreeUserOpts(), PrintUserOpts() and UserOptsLogicCheck() etc
  90      */
  91 } cmdline_opts_t;
  92
  93
  94
  95 /* log-file used for non-essential logging in prLog */
  96 FILE *prLogFile = NULL;
  97
  98
  99
 100 /**
 101  * @brief Sets default user/commandline options
 102  *
 103  * @param[out] opts
 104  * The option structure to initialise
 105  *
 106  */
 107 void
 108 SetDefaultUserOpts(cmdline_opts_t *opts)
 109 {
 110
 111     assert(NULL != opts);
 112
 113     opts->iSeqType = SEQTYPE_UNKNOWN;
 114     opts->pcSeqInfile = NULL;
 115     opts->bDealignInputSeqs = FALSE;
 116     opts->pcProfile1Infile = NULL;
 117     opts->pcProfile2Infile = NULL;
 118
 119     opts->iMaxNumSeq = INT_MAX;
 120     opts->iMaxSeqLen = INT_MAX;
 121
 122     opts->pcAlnOutfile = NULL;
 123     opts->iAlnOutFormat =  MSAFILE_A2M;
 124     opts->bForceFileOverwrite = FALSE;
 125
 126 #ifdef HAVE_OPENMP
 127     /* defaults to # of CPUs */
 128     opts->iThreads = omp_get_max_threads();
 129 #else
 130     opts->iThreads = 1;
 131 #endif
 132
 133     opts->pcLogFile = NULL;
 134
 135     SetDefaultAlnOpts(& opts->aln_opts);
 136 }
 137 /* end of SetDefaultAlnOpts() */
 138
 139
 140
 141
 142 /**
 143  * @brief FIXME add doc
 144  *
 145  */
 146 void
 147 PrintUserOpts(FILE *prFile, cmdline_opts_t *opts) {
 148
 149     /* keep in same order as in struct. FIXME could this be derived from argtable?
 150      */
 151     /* we only allow protein anyway: fprintf(prFile, "seq-type = %s\n", SeqTypeToStr(opts->iSeqType)); */
 152     fprintf(prFile, "option: seq-in = %s\n",
 153             NULL != opts->pcSeqInfile? opts->pcSeqInfile: "(null)");
 154     fprintf(prFile, "option: dealign = %d\n", opts->bDealignInputSeqs);
 155     fprintf(prFile, "option: profile1 = %s\n",
 156             NULL != opts->pcProfile1Infile? opts->pcProfile1Infile: "(null)");
 157     fprintf(prFile, "option: profile2 = %s\n",
 158             NULL != opts->pcProfile2Infile? opts->pcProfile2Infile: "(null)");
 159     fprintf(prFile, "option: max-num-seq = %d\n", opts->iMaxNumSeq);
 160     fprintf(prFile, "option: max-seq-len = %d\n", opts->iMaxSeqLen);
 161     fprintf(prFile, "option: aln-out-file = %s\n",
 162             NULL != opts->pcAlnOutfile? opts->pcAlnOutfile: "(null)");
 163     fprintf(prFile, "option: aln-out-format = %s\n", SeqfileFormat2String(opts->iAlnOutFormat));
 164     fprintf(prFile, "option: force-file-overwrite = %d\n", opts->bForceFileOverwrite);
 165     fprintf(prFile, "option: threads = %d\n", opts->iThreads);
 166     fprintf(prFile, "option: logFile = %s\n", opts->pcLogFile);
 167 }
 168 /* end of PrintUserOpts */
 169
 170
 171
 172 /**
 173  * @brief Frees user opt members allocated during parsing
 174  *
 175  * @param[out] user_opts
 176  * user options whose members are to free
 177  *
 178  * @see ParseCommandLine()
 179  *
 180  */
 181 void
 182 FreeUserOpts(cmdline_opts_t *user_opts)
 183 {
 184
 185     if (NULL != user_opts->pcSeqInfile) {
 186         CKFREE(user_opts->pcSeqInfile);
 187     }
 188     if (NULL != user_opts->pcProfile1Infile) {
 189         CKFREE(user_opts->pcProfile1Infile);
 190     }
 191     if (NULL != user_opts->pcProfile2Infile) {
 192         CKFREE(user_opts->pcProfile2Infile);
 193     }
 194     if (NULL != user_opts->pcAlnOutfile) {
 195         CKFREE(user_opts->pcAlnOutfile);
 196     }
 197     if (NULL != user_opts->pcLogFile) {
 198         CKFREE(user_opts->pcLogFile);
 199     }
 200
 201     FreeAlnOpts(& user_opts->aln_opts);
 202
 203     return;
 204 }
 205 /* end of FreeUserOpts() */
 206
 207
 208
 209
 210 /**
 211  * @brief Do quick&dirty logic check of used options and call Log(&rLog, LOG_FATAL, ) in case
 212  * of any inconsistencies
 213  *
 214  * @param[in] opts
 215  * option structure to check
 216  *
 217  */
 218 void
 219 UserOptsLogicCheck(cmdline_opts_t *opts)
 220 {
 221     /* sequence input
 222      *
 223      */
 224     if (NULL == opts->pcSeqInfile) {
 225         if (NULL == opts->pcProfile1Infile && NULL == opts->pcProfile2Infile) {
 226             Log(&rLog, LOG_FATAL, "No sequence input was provided. For more information try: --help");
 227         }
 228     } else {
 229         if (NULL != opts->pcProfile1Infile && NULL != opts->pcProfile2Infile) {
 230             Log(&rLog, LOG_FATAL, "Can't align two profile alignments AND a 'normal' sequence file");
 231         }
 232     }
 233     /* if a profile was given it should always be no 1, not 2 */
 234     if (NULL == opts->pcProfile1Infile && NULL != opts->pcProfile2Infile) {
 235         Log(&rLog, LOG_FATAL, "Got a second profile, but no first one.");
 236     }
 237
 238     /* alignment output
 239      */
 240     if (rLog.iLogLevelEnabled < LOG_WARN && NULL==opts->pcAlnOutfile && NULL==opts->pcLogFile) {
 241         Log(&rLog, LOG_FATAL, "%s %s",
 242               "You requested alignment output to stdout and verbose logging.",
 243               " Alignment and log messages would get mixed up.");
 244     }
 245     /* distance matrix output impossible in mbed mode, only have clusters, FS, r254 ->  */
 246     if ( (NULL != opts->aln_opts.pcDistmatOutfile) && (TRUE == opts->aln_opts.bUseMbed) ) {
 247         Log(&rLog, LOG_FATAL, "Distance Matrix output not possible in mBed mode.");
 248     }
 249 }
 250 /* end of UserOptsLogicCheck */
 251
 252
 253
 254 /**
 255  * @brief Parse command line parameters. Will exit if help/usage etc
 256  * are called or or call Log(&rLog, LOG_FATAL, ) if an error was detected.
 257  *
 258  * @param[out] user_opts
 259  * User parameter struct, with defaults already set.
 260  * @param[in] argc
 261  * mains argc
 262  * @param[in] argv
 263  * mains argv
 264  *
 265  */
 266 void
 267 ParseCommandLine(cmdline_opts_t *user_opts, int argc, char **argv)
 268 {
 269
 270      /* argtable command line parsing:
 271      * see
 272      * http://argtable.sourceforge.net/doc/argtable2-intro.html
 273      *
 274      * basic structure is: arg_xxxN:
 275      * xxx can be int, lit, db1, str, rex, file or date
 276      * If N = 0, arguments may appear zero-or-once; N = 1 means
 277      * exactly once, N = n means up to maxcount times
 278      *
 279      *
 280      * @note: changes here, might also affect main.cpp:ConvertOldCmdLine()
 281      *
 282      */
 283
 284     struct arg_rem  *rem_seq_input  = arg_rem(NULL, "\nSequence Input:");
 285     struct arg_file *opt_seqin = arg_file0("i", "in,infile",
 286                                             "{<file>,-}",
 287                                             "Multiple sequence input file (- for stdin)");
 288     struct arg_file *opt_hmm_in = arg_filen(NULL, "hmm-in", "<file>",
 289                                             /*min*/ 0, /*max*/ 128,
 290                                             "HMM input files");
 291     struct arg_lit *opt_dealign = arg_lit0(NULL, "dealign",
 292                                            "Dealign input sequences");
 293     struct arg_str *opt_seqtype = arg_str0("t", "seqtype",
 294                                            "{Protein,RNA,DNA}",
 295                                            "Force a sequence type (default: auto)");
 296     struct arg_file *opt_profile1 = arg_file0(NULL, "profile1,p1",
 297                                               "<file>",
 298                                               "Pre-aligned multiple sequence file (aligned columns will be kept fix)");
 299     struct arg_file *opt_profile2 = arg_file0(NULL, "profile2,p2",
 300                                               "<file>",
 301                                               "Pre-aligned multiple sequence file (aligned columns will be kept fix)");
 302
 303
 304     struct arg_rem  *rem_guidetree  = arg_rem(NULL, "\nClustering:");
 305     struct arg_str *opt_pairdist = arg_str0("p", "pairdist",
 306                                              "{ktuple}",
 307                                              "Pairwise alignment distance measure");
 308     struct arg_file *opt_distmat_in = arg_file0(NULL, "distmat-in",
 309                                                 "<file>",
 310                                                 "Pairwise distance matrix input file (skips distance computation)");
 311     struct arg_file *opt_distmat_out = arg_file0(NULL, "distmat-out",
 312                                                  "<file>",
 313                                                  "Pairwise distance matrix output file");
 314     struct arg_file *opt_guidetree_in = arg_file0(NULL, "guidetree-in",
 315                                                   "<file>",
 316                                                   "Guide tree input file (skips distance computation and guide-tree clustering step)");
 317     struct arg_file *opt_guidetree_out = arg_file0(NULL, "guidetree-out",
 318                                                    "<file>",
 319                                                    "Guide tree output file");
 320     /* AW: mbed is default since at least R253
 321        struct arg_lit *opt_mbed = arg_lit0(NULL, "mbed",
 322        "Fast, Mbed-like clustering for guide-tree calculation");
 323        struct arg_lit *opt_mbed_iter = arg_lit0(NULL, "mbed-iter",
 324        "Use Mbed-like clustering also during iteration");
 325     */
 326     /* Note: might be better to use arg_str (mbed=YES/NO) but I don't want to introduce an '=' into pipeline, FS, r250 -> */
 327     struct arg_lit *opt_full = arg_lit0(NULL, "full",
 328                                         "Use full distance matrix for guide-tree calculation (might be slow; mBed is default)");
 329     struct arg_lit *opt_full_iter = arg_lit0(NULL, "full-iter",
 330                                         "Use full distance matrix for guide-tree calculation during iteration (might be slowish; mBed is default)");
 331
 332     struct arg_str *opt_clustering = arg_str0("c", "clustering",
 333                                               "{UPGMA}",
 334                                               "Clustering method for guide tree");
 335
 336
 337     struct arg_rem *rem_aln_output  = arg_rem(NULL, "\nAlignment Output:");
 338     struct arg_file *opt_outfile = arg_file0("o", "out,outfile",
 339                                              "{file,-}",
 340                                              "Multiple sequence alignment output file (default: stdout)");
 341     struct arg_str *opt_outfmt = arg_str0(NULL, "outfmt",
 342                                             "{a2m=fa[sta],clu[stal],msf,phy[lip],selex,st[ockholm],vie[nna]}",
 343                                             "MSA output file format (default: fasta)");
 344
 345
 346     struct arg_rem *rem_iteration  = arg_rem(NULL, "\nIteration:");
 347     struct arg_str *opt_num_iterations = arg_str0(NULL, "iterations,iter",
 348                                                   /* FIXME "{<n>,auto}", "Number of combined guide-tree/HMM iterations"); */
 349                                                   "<n>", "Number of (combined guide-tree/HMM) iterations");
 350     struct arg_int *opt_max_guidetree_iterations = arg_int0(NULL, "max-guidetree-iterations",
 351                                                             "<n>", "Maximum number guidetree iterations");
 352     struct arg_int *opt_max_hmm_iterations = arg_int0(NULL, "max-hmm-iterations",
 353                                                       "<n>", "Maximum number of HMM iterations");
 354
 355
 356     struct arg_rem *rem_limits  = arg_rem(NULL, "\nLimits (will exit early, if exceeded):");
 357     struct arg_int *opt_max_seq = arg_int0(NULL, "maxnumseq", "<n>",
 358                                            "Maximum allowed number of sequences");
 359     struct arg_int *opt_max_seqlen = arg_int0(NULL, "maxseqlen", "<l>",
 360                                               "Maximum allowed sequence length");
 361
 362
 363     struct arg_rem *rem_misc  = arg_rem(NULL, "\nMiscellaneous:");
 364
 365     struct arg_lit *opt_autooptions = arg_lit0(NULL, "auto",
 366                                          "Set options automatically (might overwrite some of your options)");
 367     struct arg_int *opt_threads = arg_int0(NULL, "threads", "<n>",
 368                                               "Number of processors to use");
 369     struct arg_file *opt_logfile = arg_file0("l", "log",
 370                                              "<file>",
 371                                              "Log all non-essential output to this file");
 372     struct arg_lit *opt_help = arg_lit0("h", "help",
 373                                          "Print this help and exit");
 374     struct arg_lit *opt_version = arg_lit0(NULL, "version",
 375                                            "Print version information and exit");
 376     struct arg_lit *opt_long_version = arg_lit0(NULL, "long-version",
 377                                            "Print long version information and exit");
 378     struct arg_lit *opt_verbose = arg_litn("v", "verbose",
 379                                            0, 3,
 380                                            "Verbose output (increases if given multiple times)");
 381     struct arg_lit *opt_force = arg_lit0(NULL, "force",
 382                                          "Force file overwriting");
 383     struct arg_int *opt_macram = arg_int0(NULL, "MAC-RAM", "<n>", /* keep this quiet for the moment, FS r240 -> */
 384                                           NULL/*"maximum amount of RAM to use for MAC algorithm (in MB)"*/);
 385
 386
 387     struct arg_end *opt_end = arg_end(10); /* maximum number of errors
 388                                             * to store */
 389
 390     void *argtable[] = {rem_seq_input,
 391                         opt_seqin,
 392                         opt_hmm_in,
 393                         opt_dealign,
 394 #if 0
 395                         /* unused since we only support protein for now */
 396                         opt_seqtype,
 397 #endif
 398                         opt_profile1,
 399                         opt_profile2,
 400
 401                         rem_guidetree,
 402 #if 0
 403                         /* no other options then default available or not implemented */
 404                         opt_pairdist,
 405 #endif
 406                         opt_distmat_in,
 407                         opt_distmat_out,
 408                         opt_guidetree_in,
 409                         opt_guidetree_out,
 410                         opt_full, /* FS, r250 -> */
 411                         opt_full_iter, /* FS, r250 -> */
 412 #if 0
 413                         /* no other options then default available */
 414                         opt_clustering,
 415 #endif
 416                         rem_aln_output,
 417                         opt_outfile,
 418                         opt_outfmt,
 419
 420                         rem_iteration,
 421                         opt_num_iterations,
 422                         opt_max_guidetree_iterations,
 423                         opt_max_hmm_iterations,
 424
 425                         rem_limits,
 426                         opt_max_seq,
 427                         opt_max_seqlen,
 428
 429                         rem_misc,
 430                         opt_autooptions,
 431                         opt_threads,
 432                         opt_logfile,
 433                         opt_help,
 434                         opt_verbose,
 435                         opt_version,
 436                         opt_long_version,
 437                         opt_force,
 438                         opt_macram, /* FS, r240 -> r241 */
 439
 440                         opt_end};
 441     int nerrors;
 442
 443
 444     /* Verify the argtable[] entries were allocated sucessfully
 445      */
 446     if (arg_nullcheck(argtable)) {
 447         Log(&rLog, LOG_FATAL, "insufficient memory (for argtable allocation)");
 448     }
 449
 450     /* Parse the command line as defined by argtable[]
 451      */
 452     nerrors = arg_parse(argc, argv, argtable);
 453
 454     /* Special case: '--help' takes precedence over error reporting
 455      */
 456     if (opt_help->count > 0) {
 457         printf("%s - %s (%s)\n", PACKAGE_NAME, PACKAGE_VERSION, PACKAGE_CODENAME);
 458
 459         printf("\n");
 460         printf("Check http://www.clustal.org for more information and updates.\n");
 461
 462         /*printf("\n");
 463           printf("FIXME more info e.g. how it works, pointers to references etc...\n");
 464           FIXME which paper to cite etc
 465         */
 466
 467
 468         printf("\n");
 469         printf("Usage: %s", basename(argv[0]));
 470         arg_print_syntax(stdout,argtable, "\n");
 471
 472         printf("\n");
 473         printf("A typical invocation would be: %s -i my-in-seqs.fa -o my-out-seqs.fa -v\n",
 474                basename(argv[0]));
 475         printf("See below for a list of all options.\n");
 476
 477         arg_print_glossary(stdout, argtable, "  %-25s %s\n");
 478         arg_freetable(argtable, sizeof(argtable)/sizeof(argtable[0]));
 479         exit(EXIT_SUCCESS);
 480     }
 481
 482     /* Special case: '--version' takes precedence over error reporting
 483      */
 484     if (opt_version->count > 0) {
 485         printf("%s\n", PACKAGE_VERSION);
 486         arg_freetable(argtable,sizeof(argtable)/sizeof(argtable[0]));
 487         exit(EXIT_SUCCESS);
 488     }
 489
 490     /* Special case: '--long-version' takes precedence over error reporting
 491      */
 492     if (opt_long_version->count > 0) {
 493         char zcLongVersion[1024];
 494         PrintLongVersion(zcLongVersion, sizeof(zcLongVersion));
 495         printf("%s\n", zcLongVersion);
 496         arg_freetable(argtable,sizeof(argtable)/sizeof(argtable[0]));
 497         exit(EXIT_SUCCESS);
 498     }
 499
 500     /* If the parser returned any errors then display them and exit
 501      */
 502     if (nerrors > 0) {
 503         /* Display the error details contained in the arg_end struct.*/
 504         arg_print_errors(stdout, opt_end, PACKAGE);
 505         fprintf(stderr, "For more information try: %s --help\n", argv[0]);
 506         arg_freetable(argtable,sizeof(argtable)/sizeof(argtable[0]));
 507         exit(EXIT_FAILURE);
 508     }
 509
 510
 511     /* ------------------------------------------------------------
 512      *
 513      * Command line successfully parsed. Now transfer values to
 514      * user_opts. While doing so, make sure that given input files
 515      * exist and given output files are writable do not exist, or if
 516      * they do, should be overwritten.
 517      *
 518      * No logic checks here! They are done in a different function
 519      *
 520      * ------------------------------------------------------------*/
 521
 522
 523     /* not part of user_opts because it declared in src/util.h */
 524     if (0 == opt_verbose->count) {
 525         rLog.iLogLevelEnabled = LOG_WARN;
 526     } else if (1 == opt_verbose->count) {
 527         rLog.iLogLevelEnabled = LOG_INFO;
 528     } else if (2 == opt_verbose->count) {
 529         rLog.iLogLevelEnabled = LOG_VERBOSE;
 530     } else if (3 == opt_verbose->count) {
 531         rLog.iLogLevelEnabled = LOG_DEBUG;
 532     }
 533
 534     user_opts->aln_opts.bAutoOptions = opt_autooptions->count;
 535
 536     user_opts->bDealignInputSeqs = opt_dealign->count;
 537
 538     /* NOTE: full distance matrix used to be default - there was
 539        --mbed flag but no --full flag after r250 decided that mBed
 540        should be default - now need --full flag to turn off mBed.
 541        wanted to retain mBed Boolean, so simply added --full flag. if
 542        both flags set (erroneously) want --mbed to overwrite --full,
 543        therefore do --full 1st, the --mbed, FS, r250 */
 544     if (opt_full->count){
 545         user_opts->aln_opts.bUseMbed = FALSE;
 546     }
 547
 548     if (opt_full_iter->count){
 549         user_opts->aln_opts.bUseMbedForIteration = FALSE;
 550     }
 551
 552     user_opts->bForceFileOverwrite = opt_force->count;
 553
 554
 555
 556     /* log-file
 557      */
 558     if (opt_logfile->count > 0) {
 559         user_opts->pcLogFile = CkStrdup(opt_logfile->filename[0]);
 560
 561         /* warn if already exists or not writable */
 562         if (FileExists(user_opts->pcLogFile) && ! user_opts->bForceFileOverwrite) {
 563             Log(&rLog, LOG_FATAL, "%s '%s'. %s",
 564                   "Cowardly refusing to overwrite already existing file",
 565                   user_opts->pcLogFile,
 566                   "Use --force to force overwriting.");
 567         }
 568         if (! FileIsWritable(user_opts->pcLogFile)) {
 569             Log(&rLog, LOG_FATAL, "Sorry, I do not have permission to write to file '%s'.",
 570                   user_opts->pcLogFile);
 571         }
 572     }
 573
 574
 575     /* normal sequence input (no profile)
 576      */
 577     if (opt_seqin->count > 0) {
 578         user_opts->pcSeqInfile = CkStrdup(opt_seqin->filename[0]);
 579     }
 580
 581     /* Input limitations
 582      */
 583     /* maximum number of sequences */
 584     if (opt_max_seq->count > 0) {
 585         user_opts->iMaxNumSeq = opt_max_seq->ival[0];
 586     }
 587
 588     /* maximum sequence length */
 589     if (opt_max_seqlen->count > 0) {
 590         user_opts->iMaxSeqLen = opt_max_seqlen->ival[0];
 591     }
 592
 593     /* Sequence type
 594      */
 595     if (opt_seqtype->count > 0) {
 596         if (STR_NC_EQ(opt_seqtype->sval[0], "protein")) {
 597             user_opts->iSeqType = SEQTYPE_PROTEIN;
 598         } else if (STR_NC_EQ(opt_seqtype->sval[0], "rna")) {
 599             user_opts->iSeqType = SEQTYPE_RNA;
 600         } else if  (STR_NC_EQ(opt_seqtype->sval[0], "dna")) {
 601             user_opts->iSeqType = SEQTYPE_DNA;
 602         } else {
 603             Log(&rLog, LOG_FATAL, "Unknown sequence type '%s'", opt_seqtype->sval[0]);
 604         }
 605     }
 606
 607     /* Profile input
 608      */
 609     if (opt_profile1->count > 0) {
 610         user_opts->pcProfile1Infile = CkStrdup(opt_profile1->filename[0]);
 611         if (! FileExists(user_opts->pcProfile1Infile)) {
 612             Log(&rLog, LOG_FATAL, "File '%s' does not exist.", user_opts->pcProfile1Infile);
 613         }
 614     }
 615
 616     if (opt_profile2->count > 0) {
 617         user_opts->pcProfile2Infile = CkStrdup(opt_profile2->filename[0]);
 618         if (! FileExists(user_opts->pcProfile2Infile)) {
 619             Log(&rLog, LOG_FATAL, "File '%s' does not exist.", user_opts->pcProfile2Infile);
 620         }
 621     }
 622
 623
 624     /* HMM input
 625      */
 626     user_opts->aln_opts.iHMMInputFiles = 0;
 627     user_opts->aln_opts.ppcHMMInput = NULL;
 628     if (opt_hmm_in->count>0) {
 629         int iAux;
 630         user_opts->aln_opts.iHMMInputFiles = opt_hmm_in->count;
 631         user_opts->aln_opts.ppcHMMInput = (char **) CKMALLOC(
 632             user_opts->aln_opts.iHMMInputFiles * sizeof(char*));
 633         for (iAux=0; iAux<opt_hmm_in->count; iAux++) {
 634             user_opts->aln_opts.ppcHMMInput[iAux] = CkStrdup(opt_hmm_in->filename[iAux]);
 635             if (! FileExists(user_opts->aln_opts.ppcHMMInput[iAux])) {
 636                 Log(&rLog, LOG_FATAL, "File '%s' does not exist.", user_opts->aln_opts.ppcHMMInput[iAux]);
 637             }
 638         }
 639     }
 640
 641
 642     /* Pair distance method
 643      */
 644     if (opt_pairdist->count > 0) {
 645         if (STR_NC_EQ(opt_pairdist->sval[0], "ktuple")) {
 646             user_opts->aln_opts.iPairDistType = PAIRDIST_KTUPLE;
 647         } else {
 648             Log(&rLog, LOG_FATAL, "Unknown pairdist method '%s'", opt_pairdist->sval[0]);
 649         }
 650     }
 651
 652
 653     /* Distance matrix input
 654      */
 655     if (opt_distmat_in->count > 0) {
 656         user_opts->aln_opts.pcDistmatInfile = CkStrdup(opt_distmat_in->filename[0]);
 657         if (! FileExists(user_opts->aln_opts.pcDistmatInfile)) {
 658             Log(&rLog, LOG_FATAL, "File '%s' does not exist.", user_opts->aln_opts.pcDistmatInfile);
 659         }
 660     }
 661
 662
 663     /* Distance matrix output
 664      */
 665     if (opt_distmat_out->count > 0) {
 666         user_opts->aln_opts.pcDistmatOutfile = CkStrdup(opt_distmat_out->filename[0]);
 667
 668         /* warn if already exists or not writable */
 669         if (FileExists(user_opts->aln_opts.pcDistmatOutfile) && ! user_opts->bForceFileOverwrite) {
 670             Log(&rLog, LOG_FATAL, "%s '%s'. %s",
 671                   "Cowardly refusing to overwrite already existing file",
 672                   user_opts->aln_opts.pcDistmatOutfile,
 673                   "Use --force to force overwriting.");
 674         }
 675         if (! FileIsWritable(user_opts->aln_opts.pcDistmatOutfile)) {
 676             Log(&rLog, LOG_FATAL, "Sorry, I do not have permission to write to file '%s'.",
 677                 user_opts->aln_opts.pcDistmatOutfile);
 678         }
 679     }
 680
 681     /* Clustering
 682      *
 683      */
 684     if (opt_clustering->count > 0) {
 685         if (STR_NC_EQ(opt_clustering->sval[0], "upgma")) {
 686             user_opts->aln_opts.iClusteringType = CLUSTERING_UPGMA;
 687         } else {
 688             Log(&rLog, LOG_FATAL, "Unknown guide-tree clustering method '%s'", opt_clustering->sval[0]);
 689         }
 690     }
 691
 692
 693     /* Guidetree input
 694      */
 695     if (opt_guidetree_in->count > 0) {
 696         user_opts->aln_opts.pcGuidetreeInfile = CkStrdup(opt_guidetree_in->filename[0]);
 697         if (! FileExists(user_opts->aln_opts.pcGuidetreeInfile)) {
 698             Log(&rLog, LOG_FATAL, "File '%s' does not exist.", user_opts->aln_opts.pcGuidetreeInfile);
 699         }
 700     }
 701
 702
 703     /* Guidetree output
 704      */
 705     if (opt_guidetree_out->count > 0) {
 706         user_opts->aln_opts.pcGuidetreeOutfile = CkStrdup(opt_guidetree_out->filename[0]);
 707
 708         /* warn if already exists or not writable */
 709         if (FileExists(user_opts->aln_opts.pcGuidetreeOutfile) && ! user_opts->bForceFileOverwrite) {
 710             Log(&rLog, LOG_FATAL, "%s '%s'. %s",
 711                   "Cowardly refusing to overwrite already existing file",
 712                   user_opts->aln_opts.pcGuidetreeOutfile,
 713                   "Use --force to force overwriting.");
 714         }
 715         if (! FileIsWritable(user_opts->aln_opts.pcGuidetreeOutfile)) {
 716             Log(&rLog, LOG_FATAL, "Sorry, I do not have permission to write to file '%s'.",
 717                   user_opts->aln_opts.pcGuidetreeOutfile);
 718         }
 719     }
 720
 721
 722     /* max guidetree iterations
 723      */
 724     if (opt_max_guidetree_iterations->count > 0) {
 725         user_opts->aln_opts.iMaxGuidetreeIterations = opt_max_guidetree_iterations->ival[0];
 726     }
 727
 728
 729     /* max guidetree iterations
 730      */
 731     if (opt_max_hmm_iterations->count > 0) {
 732         user_opts->aln_opts.iMaxHMMIterations = opt_max_hmm_iterations->ival[0];
 733     }
 734
 735      /* number of iterations
 736       */
 737      if (opt_num_iterations->count > 0) {
 738         if (STR_NC_EQ(opt_num_iterations->sval[0], "auto")) {
 739             Log(&rLog, LOG_FATAL, "Automatic iteration not supported at the moment.");
 740             user_opts->aln_opts.bIterationsAuto = TRUE;
 741
 742         } else {
 743             int iAux;
 744             user_opts->aln_opts.bIterationsAuto = FALSE;
 745             for (iAux=0; iAux<(int)strlen(opt_num_iterations->sval[0]); iAux++) {
 746                 if (! isdigit(opt_num_iterations->sval[0][iAux])) {
 747                     Log(&rLog, LOG_FATAL, "Couldn't iteration parameter: %s",
 748                           opt_num_iterations->sval[0]);
 749                 }
 750             }
 751             user_opts->aln_opts.iNumIterations = atoi(opt_num_iterations->sval[0]);
 752         }
 753     }
 754
 755
 756     /* Alignment output
 757      */
 758     if (opt_outfile->count > 0) {
 759         user_opts->pcAlnOutfile = CkStrdup(opt_outfile->filename[0]);
 760
 761         /* warn if already exists or not writable */
 762         if (FileExists(user_opts->pcAlnOutfile) && ! user_opts->bForceFileOverwrite) {
 763             Log(&rLog, LOG_FATAL, "%s '%s'. %s",
 764                   "Cowardly refusing to overwrite already existing file",
 765                   user_opts->pcAlnOutfile,
 766                   "Use --force to force overwriting.");
 767         }
 768         if (! FileIsWritable(user_opts->pcAlnOutfile)) {
 769             Log(&rLog, LOG_FATAL, "Sorry, I do not have permission to write to file '%s'.",
 770                   user_opts->pcAlnOutfile);
 771         }
 772     }
 773
 774
 775     /* Output format
 776      */
 777     if (opt_outfmt->count > 0) {
 778         /* avoid gcc warning about discarded qualifier */
 779         char *tmp = (char *)opt_outfmt->sval[0];
 780         user_opts->iAlnOutFormat = String2SeqfileFormat(tmp);
 781         if (SQFILE_UNKNOWN == user_opts->iAlnOutFormat) {
 782             Log(&rLog, LOG_FATAL, "Unknown output format '%s'", opt_outfmt->sval[0]);
 783         }
 784     }
 785
 786     /* Number of threads
 787      */
 788 #ifdef HAVE_OPENMP
 789     if (opt_threads->count > 0) {
 790         if (opt_threads->ival[0] <= 0) {
 791             Log(&rLog, LOG_FATAL, "Changing number of threads to %d doesn't make sense.",
 792                   opt_threads->ival[0]);
 793         }
 794         user_opts->iThreads = opt_threads->ival[0];
 795     }
 796
 797 #else
 798     if (opt_threads->count > 0) {
 799         if (opt_threads->ival[0] > 1) {
 800             Log(&rLog, LOG_FATAL, "Cannot change number of threads to %d. %s was build without OpenMP support.",
 801                   opt_threads->ival[0], PACKAGE_NAME);
 802         }
 803     }
 804 #endif
 805
 806
 807     /* max MAC RAM (maximum amount of RAM set aside for MAC algorithm)
 808      */
 809     if (opt_macram->count > 0) { /* FS, r240 -> r241 */
 810         user_opts->aln_opts.iMacRam = opt_macram->ival[0];
 811     }
 812
 813
 814
 815     arg_freetable(argtable,sizeof(argtable)/sizeof(argtable[0]));
 816
 817     UserOptsLogicCheck(user_opts);
 818
 819     return;
 820 }
 821 /* end of ParseCommandLine() */
 822
 823
 824
 825
 826 /**
 827  *
 828  * @brief the 'real' main function
 829  *
 830  */
 831 int
 832 MyMain(int argc, char **argv)
 833 {
 834     mseq_t *prMSeq = NULL;
 835     mseq_t *prMSeqProfile1 = NULL;
 836     mseq_t *prMSeqProfile2 = NULL;
 837     cmdline_opts_t cmdline_opts;
 838     hhalign_para rHhalignPara = {0};
 839
 840     /* Must happen first: setup logger */
 841     LogDefaultSetup(&rLog);
 842
 843     /*Log(&rLog, LOG_WARN, "This is a non-public realase of %s. Please do not distribute.", PACKAGE_NAME);*/
 844     Log(&rLog, LOG_WARN, "This is a beta version of %s, for protein only.", PACKAGE_NAME); /* FS, r237 -> 238 */
 845
 846     SetDefaultUserOpts(&(cmdline_opts));
 847
 848     ParseCommandLine(&cmdline_opts, argc, argv);
 849
 850     if (NULL != cmdline_opts.pcLogFile) {
 851         prLogFile = fopen(cmdline_opts.pcLogFile, "w");
 852         LogSetFP(&rLog, LOG_INFO, prLogFile);
 853         LogSetFP(&rLog, LOG_VERBOSE, prLogFile);
 854         LogSetFP(&rLog, LOG_DEBUG, prLogFile);
 855     }
 856
 857     InitClustalOmega(cmdline_opts.iThreads);
 858
 859     if (rLog.iLogLevelEnabled < LOG_INFO) {
 860         PrintUserOpts(LogGetFP(&rLog, LOG_INFO), & cmdline_opts);
 861         PrintAlnOpts(LogGetFP(&rLog, LOG_INFO), & (cmdline_opts.aln_opts));
 862     }
 863     /* write relevant command-line options for hhalign
 864      *
 865      */
 866     rHhalignPara.iMacRamMB = cmdline_opts.aln_opts.iMacRam;
 867
 868     /* Read sequence file
 869      *
 870      */
 871     if (NULL != cmdline_opts.pcSeqInfile) {
 872         NewMSeq(&prMSeq);
 873         if (ReadSequences(prMSeq, cmdline_opts.pcSeqInfile,
 874                           cmdline_opts.iSeqType,
 875                           cmdline_opts.iMaxNumSeq, cmdline_opts.iMaxSeqLen)) {
 876             Log(&rLog, LOG_FATAL, "Reading sequence file '%s' failed", cmdline_opts.pcSeqInfile);
 877         }
 878 #if TRACE
 879         {
 880             int iAux;
 881             for (iAux=0; iAux<prMSeq->nseqs; iAux++) {
 882                 Log(&rLog, LOG_FORCED_DEBUG, "seq no %d: seq = %s", iAux, prMSeq->seq[iAux]);
 883                 LogSqInfo(&prMSeq->sqinfo[iAux]);
 884             }
 885         }
 886 #endif
 887     }
 888     /* k-tuple pairwise distance calculation seg-faults if
 889      * only one sequence, simply exit early.
 890      * note that for profile/profile alignment prMSeq is NULL
 891      * FS, r222->r223 */
 892     if (prMSeq && (prMSeq->nseqs <= 1)){
 893         Log(&rLog, LOG_FATAL, "File '%s' contains %d sequence%s, nothing to align",
 894               cmdline_opts.pcSeqInfile, prMSeq->nseqs, 1==prMSeq->nseqs?"":"s");
 895     }
 896
 897     /* Dealign if requested and neccessary
 898      */
 899     if (NULL != prMSeq) {
 900         if (TRUE == prMSeq->aligned && cmdline_opts.bDealignInputSeqs) {
 901             Log(&rLog, LOG_INFO, "Dealigning already aligned input sequences as requested.");
 902             DealignMSeq(prMSeq);
 903         }
 904     }
 905
 906
 907     /* Read profile1
 908      *
 909      */
 910     if (NULL != cmdline_opts.pcProfile1Infile) {
 911         NewMSeq(&prMSeqProfile1);
 912         if (ReadSequences(prMSeqProfile1, cmdline_opts.pcProfile1Infile,
 913                           cmdline_opts.iSeqType,
 914                           cmdline_opts.iMaxNumSeq, cmdline_opts.iMaxSeqLen)) {
 915             Log(&rLog, LOG_FATAL, "Reading sequences from profile file '%s' failed",
 916                   cmdline_opts.pcProfile1Infile);
 917         }
 918         /* FIXME: commented out. FS, r240 -> r241
 919          * for explanation see below */
 920         /*if (1==prMSeqProfile1->nseqs) {
 921           Log(&rLog, LOG_FATAL, "'%s' contains only one sequence and can therefore not be used as a profile",
 922           cmdline_opts.pcProfile1Infile);
 923           }*/
 924         if (FALSE == prMSeqProfile1->aligned) {
 925             Log(&rLog, LOG_FATAL, "Sequences in '%s' are not aligned, i.e. this is not a profile",
 926                   cmdline_opts.pcProfile1Infile);
 927         }
 928     }
 929
 930
 931
 932     /* Read profile2
 933      *
 934      */
 935     if (NULL != cmdline_opts.pcProfile2Infile) {
 936         NewMSeq(&prMSeqProfile2);
 937         if (ReadSequences(prMSeqProfile2, cmdline_opts.pcProfile2Infile,
 938                           cmdline_opts.iSeqType,
 939                           cmdline_opts.iMaxNumSeq, cmdline_opts.iMaxSeqLen)) {
 940             Log(&rLog, LOG_FATAL, "Reading sequences from profile file '%s' failed",
 941                   cmdline_opts.pcProfile2Infile);
 942         }
 943         /* FIXME: there is no (clean) way to align a single sequence to a profile.
 944          * if we go down the -i route, it causes a seg-fault in the pair-wise
 945          * k-tuple distance calculation. However, single sequences can be
 946          * understood as 1-profiles. Therefore we have to allow for 1-profiles.
 947          * FS, r240 -> r241
 948          */
 949         /*if (1==prMSeqProfile2->nseqs) {
 950           Log(&rLog, LOG_FATAL, "'%s' contains only one sequence and can therefore not be used as a profile",
 951           cmdline_opts.pcProfile2Infile);
 952           }*/
 953         if (FALSE == prMSeqProfile1->aligned) {
 954             Log(&rLog, LOG_FATAL, "Sequences in '%s' are not aligned, i.e. this is not a profile",
 955                   cmdline_opts.pcProfile2Infile);
 956         }
 957     }
 958
 959
 960     /* Depending on the input we got perform
 961      *
 962      * (i) normal alignment: seq + optional profile
 963      * or
 964      * (ii) profile profile alignment
 965      *
 966      */
 967     if (NULL != prMSeq) {
 968         if (Align(prMSeq, prMSeqProfile1, & cmdline_opts.aln_opts, rHhalignPara)) {
 969             Log(&rLog, LOG_FATAL, "An error occured during the alignment");
 970         }
 971
 972         if (WriteAlignment(prMSeq, cmdline_opts.pcAlnOutfile,
 973                            cmdline_opts.iAlnOutFormat)) {
 974             Log(&rLog, LOG_FATAL, "Could not save alignment to %s", cmdline_opts.pcAlnOutfile);
 975         }
 976 #if 0
 977         {
 978             bool bSampling = FALSE; /* better set to TRUE for many sequences */
 979             bool bReportAll = TRUE;
 980             AliStat(prMSeq, bSampling, bReportAll);
 981         }
 982 #endif
 983
 984
 985     } else if (NULL != prMSeqProfile1 && NULL != prMSeqProfile2) {
 986         if (AlignProfiles(prMSeqProfile1, prMSeqProfile2, rHhalignPara)) {
 987             Log(&rLog, LOG_FATAL, "An error occured during the alignment");
 988         }
 989         if (WriteAlignment(prMSeqProfile1, cmdline_opts.pcAlnOutfile,
 990                            cmdline_opts.iAlnOutFormat)) {
 991             Log(&rLog, LOG_FATAL, "Could not save alignment to %s", cmdline_opts.pcAlnOutfile);
 992         }
 993     }
 994
 995
 996     /* cleanup
 997      */
 998     if (NULL != prMSeq) {
 999         FreeMSeq(&prMSeq);
1000     }
1001     if (NULL != prMSeqProfile1) {
1002         FreeMSeq(&prMSeqProfile1);
1003     }
1004     if (NULL != prMSeqProfile2) {
1005         FreeMSeq(&prMSeqProfile2);
1006     }
1007
1008     FreeUserOpts(&cmdline_opts);
1009
1010     Log(&rLog, LOG_DEBUG, "Successful program exit");
1011
1012     if (NULL != cmdline_opts.pcLogFile) {
1013         fclose(prLogFile);
1014     }
1015     return EXIT_SUCCESS;
1016 }
1017 /*  end of MyMain() */
1018
1019