Wrapper for Clustal Omega.
[jabaws.git] / 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