--- /dev/null
+/* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
+
+/*********************************************************************
+ * Clustal Omega - Multiple sequence alignment
+ *
+ * Copyright (C) 2010 University College Dublin
+ *
+ * Clustal-Omega is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation; either version 2 of the
+ * License, or (at your option) any later version.
+ *
+ * This file is part of Clustal-Omega.
+ *
+ ********************************************************************/
+
+/*
+ * RCS $Id: mymain.c 255 2011-06-22 15:59:07Z fabian $
+ */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <stdio.h>
+#include <string.h>
+#include <assert.h>
+#include <unistd.h>
+#include <argtable2.h>
+#include <ctype.h>
+#include <limits.h>
+#include <libgen.h> /* for basename only */
+
+/* clustal */
+#include "clustal-omega.h"
+
+#include "mymain.h"
+
+/* hhalign (parameters) */
+#include "hhalign/general.h"
+
+typedef struct {
+
+ /* Sequence input
+ */
+ /** sequence type (from cmdline arg) */
+ int iSeqType;
+ /** sequence input file. not directly used by Align() */
+ char *pcSeqInfile;
+ /** Dealign sequences on input. Otherwise we use the alignment
+ * info for background-HMM creation */
+ bool bDealignInputSeqs;
+
+ /* profiles: : pre-aligned sequences, whose alignment will not be changed
+ */
+ /** profile 1: pre-aligned sequence input. not directly used by Align() */
+ char *pcProfile1Infile ;
+ /** profile 2: pre-aligned sequence input. not directly used by Align() */
+ char *pcProfile2Infile;
+
+ /** Input limitations
+ */
+ /** maximum allowed number of input sequences */
+ int iMaxNumSeq;
+ /** maximum allowed input sequence length */
+ int iMaxSeqLen;
+
+ /* Alignment output
+ */
+ /** alignment output file */
+ char *pcAlnOutfile;
+ /** alignment output format */
+ int iAlnOutFormat;
+ /** force overwriting of files */
+ bool bForceFileOverwrite;
+
+ /* multithreading
+ */
+ /** number of threads */
+ int iThreads;
+
+ /* logging
+ */
+ char *pcLogFile;
+
+ opts_t aln_opts;
+
+ /* changes here will have to be reflected in SetDefaultUserOpts(),
+ * FreeUserOpts(), PrintUserOpts() and UserOptsLogicCheck() etc
+ */
+} cmdline_opts_t;
+
+
+
+/* log-file used for non-essential logging in prLog */
+FILE *prLogFile = NULL;
+
+
+
+/**
+ * @brief Sets default user/commandline options
+ *
+ * @param[out] opts
+ * The option structure to initialise
+ *
+ */
+void
+SetDefaultUserOpts(cmdline_opts_t *opts)
+{
+
+ assert(NULL != opts);
+
+ opts->iSeqType = SEQTYPE_UNKNOWN;
+ opts->pcSeqInfile = NULL;
+ opts->bDealignInputSeqs = FALSE;
+ opts->pcProfile1Infile = NULL;
+ opts->pcProfile2Infile = NULL;
+
+ opts->iMaxNumSeq = INT_MAX;
+ opts->iMaxSeqLen = INT_MAX;
+
+ opts->pcAlnOutfile = NULL;
+ opts->iAlnOutFormat = MSAFILE_A2M;
+ opts->bForceFileOverwrite = FALSE;
+
+#ifdef HAVE_OPENMP
+ /* defaults to # of CPUs */
+ opts->iThreads = omp_get_max_threads();
+#else
+ opts->iThreads = 1;
+#endif
+
+ opts->pcLogFile = NULL;
+
+ SetDefaultAlnOpts(& opts->aln_opts);
+}
+/* end of SetDefaultAlnOpts() */
+
+
+
+
+/**
+ * @brief FIXME add doc
+ *
+ */
+void
+PrintUserOpts(FILE *prFile, cmdline_opts_t *opts) {
+
+ /* keep in same order as in struct. FIXME could this be derived from argtable?
+ */
+ /* we only allow protein anyway: fprintf(prFile, "seq-type = %s\n", SeqTypeToStr(opts->iSeqType)); */
+ fprintf(prFile, "option: seq-in = %s\n",
+ NULL != opts->pcSeqInfile? opts->pcSeqInfile: "(null)");
+ fprintf(prFile, "option: dealign = %d\n", opts->bDealignInputSeqs);
+ fprintf(prFile, "option: profile1 = %s\n",
+ NULL != opts->pcProfile1Infile? opts->pcProfile1Infile: "(null)");
+ fprintf(prFile, "option: profile2 = %s\n",
+ NULL != opts->pcProfile2Infile? opts->pcProfile2Infile: "(null)");
+ fprintf(prFile, "option: max-num-seq = %d\n", opts->iMaxNumSeq);
+ fprintf(prFile, "option: max-seq-len = %d\n", opts->iMaxSeqLen);
+ fprintf(prFile, "option: aln-out-file = %s\n",
+ NULL != opts->pcAlnOutfile? opts->pcAlnOutfile: "(null)");
+ fprintf(prFile, "option: aln-out-format = %s\n", SeqfileFormat2String(opts->iAlnOutFormat));
+ fprintf(prFile, "option: force-file-overwrite = %d\n", opts->bForceFileOverwrite);
+ fprintf(prFile, "option: threads = %d\n", opts->iThreads);
+ fprintf(prFile, "option: logFile = %s\n", opts->pcLogFile);
+}
+/* end of PrintUserOpts */
+
+
+
+/**
+ * @brief Frees user opt members allocated during parsing
+ *
+ * @param[out] user_opts
+ * user options whose members are to free
+ *
+ * @see ParseCommandLine()
+ *
+ */
+void
+FreeUserOpts(cmdline_opts_t *user_opts)
+{
+
+ if (NULL != user_opts->pcSeqInfile) {
+ CKFREE(user_opts->pcSeqInfile);
+ }
+ if (NULL != user_opts->pcProfile1Infile) {
+ CKFREE(user_opts->pcProfile1Infile);
+ }
+ if (NULL != user_opts->pcProfile2Infile) {
+ CKFREE(user_opts->pcProfile2Infile);
+ }
+ if (NULL != user_opts->pcAlnOutfile) {
+ CKFREE(user_opts->pcAlnOutfile);
+ }
+ if (NULL != user_opts->pcLogFile) {
+ CKFREE(user_opts->pcLogFile);
+ }
+
+ FreeAlnOpts(& user_opts->aln_opts);
+
+ return;
+}
+/* end of FreeUserOpts() */
+
+
+
+
+/**
+ * @brief Do quick&dirty logic check of used options and call Log(&rLog, LOG_FATAL, ) in case
+ * of any inconsistencies
+ *
+ * @param[in] opts
+ * option structure to check
+ *
+ */
+void
+UserOptsLogicCheck(cmdline_opts_t *opts)
+{
+ /* sequence input
+ *
+ */
+ if (NULL == opts->pcSeqInfile) {
+ if (NULL == opts->pcProfile1Infile && NULL == opts->pcProfile2Infile) {
+ Log(&rLog, LOG_FATAL, "No sequence input was provided. For more information try: --help");
+ }
+ } else {
+ if (NULL != opts->pcProfile1Infile && NULL != opts->pcProfile2Infile) {
+ Log(&rLog, LOG_FATAL, "Can't align two profile alignments AND a 'normal' sequence file");
+ }
+ }
+ /* if a profile was given it should always be no 1, not 2 */
+ if (NULL == opts->pcProfile1Infile && NULL != opts->pcProfile2Infile) {
+ Log(&rLog, LOG_FATAL, "Got a second profile, but no first one.");
+ }
+
+ /* alignment output
+ */
+ if (rLog.iLogLevelEnabled < LOG_WARN && NULL==opts->pcAlnOutfile && NULL==opts->pcLogFile) {
+ Log(&rLog, LOG_FATAL, "%s %s",
+ "You requested alignment output to stdout and verbose logging.",
+ " Alignment and log messages would get mixed up.");
+ }
+ /* distance matrix output impossible in mbed mode, only have clusters, FS, r254 -> */
+ if ( (NULL != opts->aln_opts.pcDistmatOutfile) && (TRUE == opts->aln_opts.bUseMbed) ) {
+ Log(&rLog, LOG_FATAL, "Distance Matrix output not possible in mBed mode.");
+ }
+}
+/* end of UserOptsLogicCheck */
+
+
+
+/**
+ * @brief Parse command line parameters. Will exit if help/usage etc
+ * are called or or call Log(&rLog, LOG_FATAL, ) if an error was detected.
+ *
+ * @param[out] user_opts
+ * User parameter struct, with defaults already set.
+ * @param[in] argc
+ * mains argc
+ * @param[in] argv
+ * mains argv
+ *
+ */
+void
+ParseCommandLine(cmdline_opts_t *user_opts, int argc, char **argv)
+{
+
+ /* argtable command line parsing:
+ * see
+ * http://argtable.sourceforge.net/doc/argtable2-intro.html
+ *
+ * basic structure is: arg_xxxN:
+ * xxx can be int, lit, db1, str, rex, file or date
+ * If N = 0, arguments may appear zero-or-once; N = 1 means
+ * exactly once, N = n means up to maxcount times
+ *
+ *
+ * @note: changes here, might also affect main.cpp:ConvertOldCmdLine()
+ *
+ */
+
+ struct arg_rem *rem_seq_input = arg_rem(NULL, "\nSequence Input:");
+ struct arg_file *opt_seqin = arg_file0("i", "in,infile",
+ "{<file>,-}",
+ "Multiple sequence input file (- for stdin)");
+ struct arg_file *opt_hmm_in = arg_filen(NULL, "hmm-in", "<file>",
+ /*min*/ 0, /*max*/ 128,
+ "HMM input files");
+ struct arg_lit *opt_dealign = arg_lit0(NULL, "dealign",
+ "Dealign input sequences");
+ struct arg_str *opt_seqtype = arg_str0("t", "seqtype",
+ "{Protein,RNA,DNA}",
+ "Force a sequence type (default: auto)");
+ struct arg_file *opt_profile1 = arg_file0(NULL, "profile1,p1",
+ "<file>",
+ "Pre-aligned multiple sequence file (aligned columns will be kept fix)");
+ struct arg_file *opt_profile2 = arg_file0(NULL, "profile2,p2",
+ "<file>",
+ "Pre-aligned multiple sequence file (aligned columns will be kept fix)");
+
+
+ struct arg_rem *rem_guidetree = arg_rem(NULL, "\nClustering:");
+ struct arg_str *opt_pairdist = arg_str0("p", "pairdist",
+ "{ktuple}",
+ "Pairwise alignment distance measure");
+ struct arg_file *opt_distmat_in = arg_file0(NULL, "distmat-in",
+ "<file>",
+ "Pairwise distance matrix input file (skips distance computation)");
+ struct arg_file *opt_distmat_out = arg_file0(NULL, "distmat-out",
+ "<file>",
+ "Pairwise distance matrix output file");
+ struct arg_file *opt_guidetree_in = arg_file0(NULL, "guidetree-in",
+ "<file>",
+ "Guide tree input file (skips distance computation and guide-tree clustering step)");
+ struct arg_file *opt_guidetree_out = arg_file0(NULL, "guidetree-out",
+ "<file>",
+ "Guide tree output file");
+ /* AW: mbed is default since at least R253
+ struct arg_lit *opt_mbed = arg_lit0(NULL, "mbed",
+ "Fast, Mbed-like clustering for guide-tree calculation");
+ struct arg_lit *opt_mbed_iter = arg_lit0(NULL, "mbed-iter",
+ "Use Mbed-like clustering also during iteration");
+ */
+ /* Note: might be better to use arg_str (mbed=YES/NO) but I don't want to introduce an '=' into pipeline, FS, r250 -> */
+ struct arg_lit *opt_full = arg_lit0(NULL, "full",
+ "Use full distance matrix for guide-tree calculation (might be slow; mBed is default)");
+ struct arg_lit *opt_full_iter = arg_lit0(NULL, "full-iter",
+ "Use full distance matrix for guide-tree calculation during iteration (might be slowish; mBed is default)");
+
+ struct arg_str *opt_clustering = arg_str0("c", "clustering",
+ "{UPGMA}",
+ "Clustering method for guide tree");
+
+
+ struct arg_rem *rem_aln_output = arg_rem(NULL, "\nAlignment Output:");
+ struct arg_file *opt_outfile = arg_file0("o", "out,outfile",
+ "{file,-}",
+ "Multiple sequence alignment output file (default: stdout)");
+ struct arg_str *opt_outfmt = arg_str0(NULL, "outfmt",
+ "{a2m=fa[sta],clu[stal],msf,phy[lip],selex,st[ockholm],vie[nna]}",
+ "MSA output file format (default: fasta)");
+
+
+ struct arg_rem *rem_iteration = arg_rem(NULL, "\nIteration:");
+ struct arg_str *opt_num_iterations = arg_str0(NULL, "iterations,iter",
+ /* FIXME "{<n>,auto}", "Number of combined guide-tree/HMM iterations"); */
+ "<n>", "Number of (combined guide-tree/HMM) iterations");
+ struct arg_int *opt_max_guidetree_iterations = arg_int0(NULL, "max-guidetree-iterations",
+ "<n>", "Maximum number guidetree iterations");
+ struct arg_int *opt_max_hmm_iterations = arg_int0(NULL, "max-hmm-iterations",
+ "<n>", "Maximum number of HMM iterations");
+
+
+ struct arg_rem *rem_limits = arg_rem(NULL, "\nLimits (will exit early, if exceeded):");
+ struct arg_int *opt_max_seq = arg_int0(NULL, "maxnumseq", "<n>",
+ "Maximum allowed number of sequences");
+ struct arg_int *opt_max_seqlen = arg_int0(NULL, "maxseqlen", "<l>",
+ "Maximum allowed sequence length");
+
+
+ struct arg_rem *rem_misc = arg_rem(NULL, "\nMiscellaneous:");
+
+ struct arg_lit *opt_autooptions = arg_lit0(NULL, "auto",
+ "Set options automatically (might overwrite some of your options)");
+ struct arg_int *opt_threads = arg_int0(NULL, "threads", "<n>",
+ "Number of processors to use");
+ struct arg_file *opt_logfile = arg_file0("l", "log",
+ "<file>",
+ "Log all non-essential output to this file");
+ struct arg_lit *opt_help = arg_lit0("h", "help",
+ "Print this help and exit");
+ struct arg_lit *opt_version = arg_lit0(NULL, "version",
+ "Print version information and exit");
+ struct arg_lit *opt_long_version = arg_lit0(NULL, "long-version",
+ "Print long version information and exit");
+ struct arg_lit *opt_verbose = arg_litn("v", "verbose",
+ 0, 3,
+ "Verbose output (increases if given multiple times)");
+ struct arg_lit *opt_force = arg_lit0(NULL, "force",
+ "Force file overwriting");
+ struct arg_int *opt_macram = arg_int0(NULL, "MAC-RAM", "<n>", /* keep this quiet for the moment, FS r240 -> */
+ NULL/*"maximum amount of RAM to use for MAC algorithm (in MB)"*/);
+
+
+ struct arg_end *opt_end = arg_end(10); /* maximum number of errors
+ * to store */
+
+ void *argtable[] = {rem_seq_input,
+ opt_seqin,
+ opt_hmm_in,
+ opt_dealign,
+#if 0
+ /* unused since we only support protein for now */
+ opt_seqtype,
+#endif
+ opt_profile1,
+ opt_profile2,
+
+ rem_guidetree,
+#if 0
+ /* no other options then default available or not implemented */
+ opt_pairdist,
+#endif
+ opt_distmat_in,
+ opt_distmat_out,
+ opt_guidetree_in,
+ opt_guidetree_out,
+ opt_full, /* FS, r250 -> */
+ opt_full_iter, /* FS, r250 -> */
+#if 0
+ /* no other options then default available */
+ opt_clustering,
+#endif
+ rem_aln_output,
+ opt_outfile,
+ opt_outfmt,
+
+ rem_iteration,
+ opt_num_iterations,
+ opt_max_guidetree_iterations,
+ opt_max_hmm_iterations,
+
+ rem_limits,
+ opt_max_seq,
+ opt_max_seqlen,
+
+ rem_misc,
+ opt_autooptions,
+ opt_threads,
+ opt_logfile,
+ opt_help,
+ opt_verbose,
+ opt_version,
+ opt_long_version,
+ opt_force,
+ opt_macram, /* FS, r240 -> r241 */
+
+ opt_end};
+ int nerrors;
+
+
+ /* Verify the argtable[] entries were allocated sucessfully
+ */
+ if (arg_nullcheck(argtable)) {
+ Log(&rLog, LOG_FATAL, "insufficient memory (for argtable allocation)");
+ }
+
+ /* Parse the command line as defined by argtable[]
+ */
+ nerrors = arg_parse(argc, argv, argtable);
+
+ /* Special case: '--help' takes precedence over error reporting
+ */
+ if (opt_help->count > 0) {
+ printf("%s - %s (%s)\n", PACKAGE_NAME, PACKAGE_VERSION, PACKAGE_CODENAME);
+
+ printf("\n");
+ printf("Check http://www.clustal.org for more information and updates.\n");
+
+ /*printf("\n");
+ printf("FIXME more info e.g. how it works, pointers to references etc...\n");
+ FIXME which paper to cite etc
+ */
+
+
+ printf("\n");
+ printf("Usage: %s", basename(argv[0]));
+ arg_print_syntax(stdout,argtable, "\n");
+
+ printf("\n");
+ printf("A typical invocation would be: %s -i my-in-seqs.fa -o my-out-seqs.fa -v\n",
+ basename(argv[0]));
+ printf("See below for a list of all options.\n");
+
+ arg_print_glossary(stdout, argtable, " %-25s %s\n");
+ arg_freetable(argtable, sizeof(argtable)/sizeof(argtable[0]));
+ exit(EXIT_SUCCESS);
+ }
+
+ /* Special case: '--version' takes precedence over error reporting
+ */
+ if (opt_version->count > 0) {
+ printf("%s\n", PACKAGE_VERSION);
+ arg_freetable(argtable,sizeof(argtable)/sizeof(argtable[0]));
+ exit(EXIT_SUCCESS);
+ }
+
+ /* Special case: '--long-version' takes precedence over error reporting
+ */
+ if (opt_long_version->count > 0) {
+ char zcLongVersion[1024];
+ PrintLongVersion(zcLongVersion, sizeof(zcLongVersion));
+ printf("%s\n", zcLongVersion);
+ arg_freetable(argtable,sizeof(argtable)/sizeof(argtable[0]));
+ exit(EXIT_SUCCESS);
+ }
+
+ /* If the parser returned any errors then display them and exit
+ */
+ if (nerrors > 0) {
+ /* Display the error details contained in the arg_end struct.*/
+ arg_print_errors(stdout, opt_end, PACKAGE);
+ fprintf(stderr, "For more information try: %s --help\n", argv[0]);
+ arg_freetable(argtable,sizeof(argtable)/sizeof(argtable[0]));
+ exit(EXIT_FAILURE);
+ }
+
+
+ /* ------------------------------------------------------------
+ *
+ * Command line successfully parsed. Now transfer values to
+ * user_opts. While doing so, make sure that given input files
+ * exist and given output files are writable do not exist, or if
+ * they do, should be overwritten.
+ *
+ * No logic checks here! They are done in a different function
+ *
+ * ------------------------------------------------------------*/
+
+
+ /* not part of user_opts because it declared in src/util.h */
+ if (0 == opt_verbose->count) {
+ rLog.iLogLevelEnabled = LOG_WARN;
+ } else if (1 == opt_verbose->count) {
+ rLog.iLogLevelEnabled = LOG_INFO;
+ } else if (2 == opt_verbose->count) {
+ rLog.iLogLevelEnabled = LOG_VERBOSE;
+ } else if (3 == opt_verbose->count) {
+ rLog.iLogLevelEnabled = LOG_DEBUG;
+ }
+
+ user_opts->aln_opts.bAutoOptions = opt_autooptions->count;
+
+ user_opts->bDealignInputSeqs = opt_dealign->count;
+
+ /* NOTE: full distance matrix used to be default - there was
+ --mbed flag but no --full flag after r250 decided that mBed
+ should be default - now need --full flag to turn off mBed.
+ wanted to retain mBed Boolean, so simply added --full flag. if
+ both flags set (erroneously) want --mbed to overwrite --full,
+ therefore do --full 1st, the --mbed, FS, r250 */
+ if (opt_full->count){
+ user_opts->aln_opts.bUseMbed = FALSE;
+ }
+
+ if (opt_full_iter->count){
+ user_opts->aln_opts.bUseMbedForIteration = FALSE;
+ }
+
+ user_opts->bForceFileOverwrite = opt_force->count;
+
+
+
+ /* log-file
+ */
+ if (opt_logfile->count > 0) {
+ user_opts->pcLogFile = CkStrdup(opt_logfile->filename[0]);
+
+ /* warn if already exists or not writable */
+ if (FileExists(user_opts->pcLogFile) && ! user_opts->bForceFileOverwrite) {
+ Log(&rLog, LOG_FATAL, "%s '%s'. %s",
+ "Cowardly refusing to overwrite already existing file",
+ user_opts->pcLogFile,
+ "Use --force to force overwriting.");
+ }
+ if (! FileIsWritable(user_opts->pcLogFile)) {
+ Log(&rLog, LOG_FATAL, "Sorry, I do not have permission to write to file '%s'.",
+ user_opts->pcLogFile);
+ }
+ }
+
+
+ /* normal sequence input (no profile)
+ */
+ if (opt_seqin->count > 0) {
+ user_opts->pcSeqInfile = CkStrdup(opt_seqin->filename[0]);
+ }
+
+ /* Input limitations
+ */
+ /* maximum number of sequences */
+ if (opt_max_seq->count > 0) {
+ user_opts->iMaxNumSeq = opt_max_seq->ival[0];
+ }
+
+ /* maximum sequence length */
+ if (opt_max_seqlen->count > 0) {
+ user_opts->iMaxSeqLen = opt_max_seqlen->ival[0];
+ }
+
+ /* Sequence type
+ */
+ if (opt_seqtype->count > 0) {
+ if (STR_NC_EQ(opt_seqtype->sval[0], "protein")) {
+ user_opts->iSeqType = SEQTYPE_PROTEIN;
+ } else if (STR_NC_EQ(opt_seqtype->sval[0], "rna")) {
+ user_opts->iSeqType = SEQTYPE_RNA;
+ } else if (STR_NC_EQ(opt_seqtype->sval[0], "dna")) {
+ user_opts->iSeqType = SEQTYPE_DNA;
+ } else {
+ Log(&rLog, LOG_FATAL, "Unknown sequence type '%s'", opt_seqtype->sval[0]);
+ }
+ }
+
+ /* Profile input
+ */
+ if (opt_profile1->count > 0) {
+ user_opts->pcProfile1Infile = CkStrdup(opt_profile1->filename[0]);
+ if (! FileExists(user_opts->pcProfile1Infile)) {
+ Log(&rLog, LOG_FATAL, "File '%s' does not exist.", user_opts->pcProfile1Infile);
+ }
+ }
+
+ if (opt_profile2->count > 0) {
+ user_opts->pcProfile2Infile = CkStrdup(opt_profile2->filename[0]);
+ if (! FileExists(user_opts->pcProfile2Infile)) {
+ Log(&rLog, LOG_FATAL, "File '%s' does not exist.", user_opts->pcProfile2Infile);
+ }
+ }
+
+
+ /* HMM input
+ */
+ user_opts->aln_opts.iHMMInputFiles = 0;
+ user_opts->aln_opts.ppcHMMInput = NULL;
+ if (opt_hmm_in->count>0) {
+ int iAux;
+ user_opts->aln_opts.iHMMInputFiles = opt_hmm_in->count;
+ user_opts->aln_opts.ppcHMMInput = (char **) CKMALLOC(
+ user_opts->aln_opts.iHMMInputFiles * sizeof(char*));
+ for (iAux=0; iAux<opt_hmm_in->count; iAux++) {
+ user_opts->aln_opts.ppcHMMInput[iAux] = CkStrdup(opt_hmm_in->filename[iAux]);
+ if (! FileExists(user_opts->aln_opts.ppcHMMInput[iAux])) {
+ Log(&rLog, LOG_FATAL, "File '%s' does not exist.", user_opts->aln_opts.ppcHMMInput[iAux]);
+ }
+ }
+ }
+
+
+ /* Pair distance method
+ */
+ if (opt_pairdist->count > 0) {
+ if (STR_NC_EQ(opt_pairdist->sval[0], "ktuple")) {
+ user_opts->aln_opts.iPairDistType = PAIRDIST_KTUPLE;
+ } else {
+ Log(&rLog, LOG_FATAL, "Unknown pairdist method '%s'", opt_pairdist->sval[0]);
+ }
+ }
+
+
+ /* Distance matrix input
+ */
+ if (opt_distmat_in->count > 0) {
+ user_opts->aln_opts.pcDistmatInfile = CkStrdup(opt_distmat_in->filename[0]);
+ if (! FileExists(user_opts->aln_opts.pcDistmatInfile)) {
+ Log(&rLog, LOG_FATAL, "File '%s' does not exist.", user_opts->aln_opts.pcDistmatInfile);
+ }
+ }
+
+
+ /* Distance matrix output
+ */
+ if (opt_distmat_out->count > 0) {
+ user_opts->aln_opts.pcDistmatOutfile = CkStrdup(opt_distmat_out->filename[0]);
+
+ /* warn if already exists or not writable */
+ if (FileExists(user_opts->aln_opts.pcDistmatOutfile) && ! user_opts->bForceFileOverwrite) {
+ Log(&rLog, LOG_FATAL, "%s '%s'. %s",
+ "Cowardly refusing to overwrite already existing file",
+ user_opts->aln_opts.pcDistmatOutfile,
+ "Use --force to force overwriting.");
+ }
+ if (! FileIsWritable(user_opts->aln_opts.pcDistmatOutfile)) {
+ Log(&rLog, LOG_FATAL, "Sorry, I do not have permission to write to file '%s'.",
+ user_opts->aln_opts.pcDistmatOutfile);
+ }
+ }
+
+ /* Clustering
+ *
+ */
+ if (opt_clustering->count > 0) {
+ if (STR_NC_EQ(opt_clustering->sval[0], "upgma")) {
+ user_opts->aln_opts.iClusteringType = CLUSTERING_UPGMA;
+ } else {
+ Log(&rLog, LOG_FATAL, "Unknown guide-tree clustering method '%s'", opt_clustering->sval[0]);
+ }
+ }
+
+
+ /* Guidetree input
+ */
+ if (opt_guidetree_in->count > 0) {
+ user_opts->aln_opts.pcGuidetreeInfile = CkStrdup(opt_guidetree_in->filename[0]);
+ if (! FileExists(user_opts->aln_opts.pcGuidetreeInfile)) {
+ Log(&rLog, LOG_FATAL, "File '%s' does not exist.", user_opts->aln_opts.pcGuidetreeInfile);
+ }
+ }
+
+
+ /* Guidetree output
+ */
+ if (opt_guidetree_out->count > 0) {
+ user_opts->aln_opts.pcGuidetreeOutfile = CkStrdup(opt_guidetree_out->filename[0]);
+
+ /* warn if already exists or not writable */
+ if (FileExists(user_opts->aln_opts.pcGuidetreeOutfile) && ! user_opts->bForceFileOverwrite) {
+ Log(&rLog, LOG_FATAL, "%s '%s'. %s",
+ "Cowardly refusing to overwrite already existing file",
+ user_opts->aln_opts.pcGuidetreeOutfile,
+ "Use --force to force overwriting.");
+ }
+ if (! FileIsWritable(user_opts->aln_opts.pcGuidetreeOutfile)) {
+ Log(&rLog, LOG_FATAL, "Sorry, I do not have permission to write to file '%s'.",
+ user_opts->aln_opts.pcGuidetreeOutfile);
+ }
+ }
+
+
+ /* max guidetree iterations
+ */
+ if (opt_max_guidetree_iterations->count > 0) {
+ user_opts->aln_opts.iMaxGuidetreeIterations = opt_max_guidetree_iterations->ival[0];
+ }
+
+
+ /* max guidetree iterations
+ */
+ if (opt_max_hmm_iterations->count > 0) {
+ user_opts->aln_opts.iMaxHMMIterations = opt_max_hmm_iterations->ival[0];
+ }
+
+ /* number of iterations
+ */
+ if (opt_num_iterations->count > 0) {
+ if (STR_NC_EQ(opt_num_iterations->sval[0], "auto")) {
+ Log(&rLog, LOG_FATAL, "Automatic iteration not supported at the moment.");
+ user_opts->aln_opts.bIterationsAuto = TRUE;
+
+ } else {
+ int iAux;
+ user_opts->aln_opts.bIterationsAuto = FALSE;
+ for (iAux=0; iAux<(int)strlen(opt_num_iterations->sval[0]); iAux++) {
+ if (! isdigit(opt_num_iterations->sval[0][iAux])) {
+ Log(&rLog, LOG_FATAL, "Couldn't iteration parameter: %s",
+ opt_num_iterations->sval[0]);
+ }
+ }
+ user_opts->aln_opts.iNumIterations = atoi(opt_num_iterations->sval[0]);
+ }
+ }
+
+
+ /* Alignment output
+ */
+ if (opt_outfile->count > 0) {
+ user_opts->pcAlnOutfile = CkStrdup(opt_outfile->filename[0]);
+
+ /* warn if already exists or not writable */
+ if (FileExists(user_opts->pcAlnOutfile) && ! user_opts->bForceFileOverwrite) {
+ Log(&rLog, LOG_FATAL, "%s '%s'. %s",
+ "Cowardly refusing to overwrite already existing file",
+ user_opts->pcAlnOutfile,
+ "Use --force to force overwriting.");
+ }
+ if (! FileIsWritable(user_opts->pcAlnOutfile)) {
+ Log(&rLog, LOG_FATAL, "Sorry, I do not have permission to write to file '%s'.",
+ user_opts->pcAlnOutfile);
+ }
+ }
+
+
+ /* Output format
+ */
+ if (opt_outfmt->count > 0) {
+ /* avoid gcc warning about discarded qualifier */
+ char *tmp = (char *)opt_outfmt->sval[0];
+ user_opts->iAlnOutFormat = String2SeqfileFormat(tmp);
+ if (SQFILE_UNKNOWN == user_opts->iAlnOutFormat) {
+ Log(&rLog, LOG_FATAL, "Unknown output format '%s'", opt_outfmt->sval[0]);
+ }
+ }
+
+ /* Number of threads
+ */
+#ifdef HAVE_OPENMP
+ if (opt_threads->count > 0) {
+ if (opt_threads->ival[0] <= 0) {
+ Log(&rLog, LOG_FATAL, "Changing number of threads to %d doesn't make sense.",
+ opt_threads->ival[0]);
+ }
+ user_opts->iThreads = opt_threads->ival[0];
+ }
+
+#else
+ if (opt_threads->count > 0) {
+ if (opt_threads->ival[0] > 1) {
+ Log(&rLog, LOG_FATAL, "Cannot change number of threads to %d. %s was build without OpenMP support.",
+ opt_threads->ival[0], PACKAGE_NAME);
+ }
+ }
+#endif
+
+
+ /* max MAC RAM (maximum amount of RAM set aside for MAC algorithm)
+ */
+ if (opt_macram->count > 0) { /* FS, r240 -> r241 */
+ user_opts->aln_opts.iMacRam = opt_macram->ival[0];
+ }
+
+
+
+ arg_freetable(argtable,sizeof(argtable)/sizeof(argtable[0]));
+
+ UserOptsLogicCheck(user_opts);
+
+ return;
+}
+/* end of ParseCommandLine() */
+
+
+
+
+/**
+ *
+ * @brief the 'real' main function
+ *
+ */
+int
+MyMain(int argc, char **argv)
+{
+ mseq_t *prMSeq = NULL;
+ mseq_t *prMSeqProfile1 = NULL;
+ mseq_t *prMSeqProfile2 = NULL;
+ cmdline_opts_t cmdline_opts;
+ hhalign_para rHhalignPara = {0};
+
+ /* Must happen first: setup logger */
+ LogDefaultSetup(&rLog);
+
+ /*Log(&rLog, LOG_WARN, "This is a non-public realase of %s. Please do not distribute.", PACKAGE_NAME);*/
+ Log(&rLog, LOG_WARN, "This is a beta version of %s, for protein only.", PACKAGE_NAME); /* FS, r237 -> 238 */
+
+ SetDefaultUserOpts(&(cmdline_opts));
+
+ ParseCommandLine(&cmdline_opts, argc, argv);
+
+ if (NULL != cmdline_opts.pcLogFile) {
+ prLogFile = fopen(cmdline_opts.pcLogFile, "w");
+ LogSetFP(&rLog, LOG_INFO, prLogFile);
+ LogSetFP(&rLog, LOG_VERBOSE, prLogFile);
+ LogSetFP(&rLog, LOG_DEBUG, prLogFile);
+ }
+
+ InitClustalOmega(cmdline_opts.iThreads);
+
+ if (rLog.iLogLevelEnabled < LOG_INFO) {
+ PrintUserOpts(LogGetFP(&rLog, LOG_INFO), & cmdline_opts);
+ PrintAlnOpts(LogGetFP(&rLog, LOG_INFO), & (cmdline_opts.aln_opts));
+ }
+ /* write relevant command-line options for hhalign
+ *
+ */
+ rHhalignPara.iMacRamMB = cmdline_opts.aln_opts.iMacRam;
+
+ /* Read sequence file
+ *
+ */
+ if (NULL != cmdline_opts.pcSeqInfile) {
+ NewMSeq(&prMSeq);
+ if (ReadSequences(prMSeq, cmdline_opts.pcSeqInfile,
+ cmdline_opts.iSeqType,
+ cmdline_opts.iMaxNumSeq, cmdline_opts.iMaxSeqLen)) {
+ Log(&rLog, LOG_FATAL, "Reading sequence file '%s' failed", cmdline_opts.pcSeqInfile);
+ }
+#if TRACE
+ {
+ int iAux;
+ for (iAux=0; iAux<prMSeq->nseqs; iAux++) {
+ Log(&rLog, LOG_FORCED_DEBUG, "seq no %d: seq = %s", iAux, prMSeq->seq[iAux]);
+ LogSqInfo(&prMSeq->sqinfo[iAux]);
+ }
+ }
+#endif
+ }
+ /* k-tuple pairwise distance calculation seg-faults if
+ * only one sequence, simply exit early.
+ * note that for profile/profile alignment prMSeq is NULL
+ * FS, r222->r223 */
+ if (prMSeq && (prMSeq->nseqs <= 1)){
+ Log(&rLog, LOG_FATAL, "File '%s' contains %d sequence%s, nothing to align",
+ cmdline_opts.pcSeqInfile, prMSeq->nseqs, 1==prMSeq->nseqs?"":"s");
+ }
+
+ /* Dealign if requested and neccessary
+ */
+ if (NULL != prMSeq) {
+ if (TRUE == prMSeq->aligned && cmdline_opts.bDealignInputSeqs) {
+ Log(&rLog, LOG_INFO, "Dealigning already aligned input sequences as requested.");
+ DealignMSeq(prMSeq);
+ }
+ }
+
+
+ /* Read profile1
+ *
+ */
+ if (NULL != cmdline_opts.pcProfile1Infile) {
+ NewMSeq(&prMSeqProfile1);
+ if (ReadSequences(prMSeqProfile1, cmdline_opts.pcProfile1Infile,
+ cmdline_opts.iSeqType,
+ cmdline_opts.iMaxNumSeq, cmdline_opts.iMaxSeqLen)) {
+ Log(&rLog, LOG_FATAL, "Reading sequences from profile file '%s' failed",
+ cmdline_opts.pcProfile1Infile);
+ }
+ /* FIXME: commented out. FS, r240 -> r241
+ * for explanation see below */
+ /*if (1==prMSeqProfile1->nseqs) {
+ Log(&rLog, LOG_FATAL, "'%s' contains only one sequence and can therefore not be used as a profile",
+ cmdline_opts.pcProfile1Infile);
+ }*/
+ if (FALSE == prMSeqProfile1->aligned) {
+ Log(&rLog, LOG_FATAL, "Sequences in '%s' are not aligned, i.e. this is not a profile",
+ cmdline_opts.pcProfile1Infile);
+ }
+ }
+
+
+
+ /* Read profile2
+ *
+ */
+ if (NULL != cmdline_opts.pcProfile2Infile) {
+ NewMSeq(&prMSeqProfile2);
+ if (ReadSequences(prMSeqProfile2, cmdline_opts.pcProfile2Infile,
+ cmdline_opts.iSeqType,
+ cmdline_opts.iMaxNumSeq, cmdline_opts.iMaxSeqLen)) {
+ Log(&rLog, LOG_FATAL, "Reading sequences from profile file '%s' failed",
+ cmdline_opts.pcProfile2Infile);
+ }
+ /* FIXME: there is no (clean) way to align a single sequence to a profile.
+ * if we go down the -i route, it causes a seg-fault in the pair-wise
+ * k-tuple distance calculation. However, single sequences can be
+ * understood as 1-profiles. Therefore we have to allow for 1-profiles.
+ * FS, r240 -> r241
+ */
+ /*if (1==prMSeqProfile2->nseqs) {
+ Log(&rLog, LOG_FATAL, "'%s' contains only one sequence and can therefore not be used as a profile",
+ cmdline_opts.pcProfile2Infile);
+ }*/
+ if (FALSE == prMSeqProfile1->aligned) {
+ Log(&rLog, LOG_FATAL, "Sequences in '%s' are not aligned, i.e. this is not a profile",
+ cmdline_opts.pcProfile2Infile);
+ }
+ }
+
+
+ /* Depending on the input we got perform
+ *
+ * (i) normal alignment: seq + optional profile
+ * or
+ * (ii) profile profile alignment
+ *
+ */
+ if (NULL != prMSeq) {
+ if (Align(prMSeq, prMSeqProfile1, & cmdline_opts.aln_opts, rHhalignPara)) {
+ Log(&rLog, LOG_FATAL, "An error occured during the alignment");
+ }
+
+ if (WriteAlignment(prMSeq, cmdline_opts.pcAlnOutfile,
+ cmdline_opts.iAlnOutFormat)) {
+ Log(&rLog, LOG_FATAL, "Could not save alignment to %s", cmdline_opts.pcAlnOutfile);
+ }
+#if 0
+ {
+ bool bSampling = FALSE; /* better set to TRUE for many sequences */
+ bool bReportAll = TRUE;
+ AliStat(prMSeq, bSampling, bReportAll);
+ }
+#endif
+
+
+ } else if (NULL != prMSeqProfile1 && NULL != prMSeqProfile2) {
+ if (AlignProfiles(prMSeqProfile1, prMSeqProfile2, rHhalignPara)) {
+ Log(&rLog, LOG_FATAL, "An error occured during the alignment");
+ }
+ if (WriteAlignment(prMSeqProfile1, cmdline_opts.pcAlnOutfile,
+ cmdline_opts.iAlnOutFormat)) {
+ Log(&rLog, LOG_FATAL, "Could not save alignment to %s", cmdline_opts.pcAlnOutfile);
+ }
+ }
+
+
+ /* cleanup
+ */
+ if (NULL != prMSeq) {
+ FreeMSeq(&prMSeq);
+ }
+ if (NULL != prMSeqProfile1) {
+ FreeMSeq(&prMSeqProfile1);
+ }
+ if (NULL != prMSeqProfile2) {
+ FreeMSeq(&prMSeqProfile2);
+ }
+
+ FreeUserOpts(&cmdline_opts);
+
+ Log(&rLog, LOG_DEBUG, "Successful program exit");
+
+ if (NULL != cmdline_opts.pcLogFile) {
+ fclose(prLogFile);
+ }
+ return EXIT_SUCCESS;
+}
+/* end of MyMain() */
+
+
git://source.jalview.org / jabaws.git / blobdiff