diff --git a/src/Makefile.in b/src/Makefile.in index 6beceb19..152ecef5 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -49,12 +49,14 @@ MYINCDIRS = -I../${ESLDIR} \ PROGS = bathsearch\ + bathalign\ bathbuild\ bathconvert\ bathfetch\ bathstat PROGOBJS = bathsearch.o\ + bathalign.o\ bathbuild.o\ bathconvert.o\ bathfetch.o\ diff --git a/src/bathalign.c b/src/bathalign.c new file mode 100644 index 00000000..54fb8278 --- /dev/null +++ b/src/bathalign.c @@ -0,0 +1,288 @@ +/* bathalign: align protein sequences to a BATH profile HMM + */ +#include "p7_config.h" + +#include +#include +#include + +#include "easel.h" +#include "esl_alphabet.h" +#include "esl_getopts.h" +#include "esl_msa.h" +#include "esl_msafile.h" +#include "esl_sq.h" +#include "esl_sqio.h" +#include "esl_vectorops.h" + +#include "hmmer.h" + +static int map_alignment(const char *msafile, const P7_HMM *hmm, ESL_SQ ***ret_sq, P7_TRACE ***ret_tr, int *ret_ntot); + + +#define ALPHOPTS "--amino,--dna,--rna" /* Exclusive options for alphabet choice */ + +static ESL_OPTIONS options[] = { + /* name type default env range toggles reqs incomp help docgroup*/ + { "-h", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "show brief help on version and usage", 1 }, + { "-o", eslARG_OUTFILE, NULL, NULL, NULL, NULL, NULL, NULL, "output alignment to file , not stdout", 1 }, + + { "--mapali", eslARG_INFILE, NULL, NULL, NULL, NULL, NULL, NULL, "include alignment in file (same ali that HMM came from)", 2 }, + { "--trim", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "trim terminal tails of nonaligned residues from alignment", 2 }, + { "--amino", eslARG_NONE, FALSE, NULL, NULL, ALPHOPTS, NULL, NULL, "assert , both protein: no autodetection", 2 }, + { "--dna", eslARG_NONE, FALSE, NULL, NULL, ALPHOPTS, NULL, NULL, "assert , both DNA: no autodetection", 2 }, + { "--rna", eslARG_NONE, FALSE, NULL, NULL, ALPHOPTS, NULL, NULL, "assert , both RNA: no autodetection", 2 }, + { "--informat", eslARG_STRING, NULL, NULL, NULL, NULL, NULL, NULL, "assert is in format : no autodetection", 2 }, + { "--outformat", eslARG_STRING, "Stockholm", NULL, NULL, NULL, NULL, NULL, "output alignment in format ", 2 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, +}; + +static char usage[] = "[-options] "; +static char banner[] = "align protein sequences to a profile HMM"; + + + + +static void +cmdline_failure(char *argv0, char *format, ...) +{ + va_list argp; + printf("\nERROR: "); + va_start(argp, format); + vfprintf(stdout, format, argp); + va_end(argp); + esl_usage(stdout, argv0, usage); + printf("\nTo see more help on available options, do %s -h\n\n", argv0); + exit(1); +} + +static void +cmdline_help(char *argv0, ESL_GETOPTS *go) +{ + p7_banner (stdout, argv0, banner); + esl_usage (stdout, argv0, usage); + puts("\nBasic options:"); + esl_opt_DisplayHelp(stdout, go, 1, 2, 80); + puts("\nLess common options:"); + esl_opt_DisplayHelp(stdout, go, 2, 2, 80); + puts("\nSequence input formats include: FASTA, EMBL, GenBank, UniProt"); + puts("Alignment output formats include: Stockholm, Pfam, A2M, PSIBLAST\n"); + exit(0); +} + + + +int +main(int argc, char **argv) +{ + ESL_GETOPTS *go = NULL; /* application configuration */ + char *hmmfile = NULL; /* HMM file name */ + char *seqfile = NULL; /* sequence file name */ + char *mapfile = NULL; /* optional mapped MSA file name */ + int infmt = eslSQFILE_UNKNOWN; + int outfmt = eslMSAFILE_STOCKHOLM; + P7_HMMFILE *hfp = NULL; /* open HMM file */ + ESL_SQFILE *sqfp = NULL; /* open sequence file */ + char *outfile = NULL; /* output filename */ + FILE *ofp = stdout; /* output stream */ + ESL_SQ **sq = NULL; /* array of sequences */ + void *p = NULL; /* tmp ptr for reallocation */ + int nseq = 0; /* # of sequences in */ + int mapseq = 0; /* # of sequences in mapped MSA */ + int totseq = 0; /* # of seqs in all sources */ + ESL_ALPHABET *abc = NULL; /* alphabet (set from the HMM file)*/ + P7_HMM *hmm = NULL; + P7_TRACE **tr = NULL; /* array of tracebacks */ + ESL_MSA *msa = NULL; /* resulting multiple alignment */ + int msaopts = 0; /* flags to p7_tracealign_Seqs() */ + int idx; /* counter over seqs, traces */ + int status; /* easel/hmmer return code */ + char errbuf[eslERRBUFSIZE]; + + /* Parse the command line + */ + go = esl_getopts_Create(options); + if (esl_opt_ProcessCmdline(go, argc, argv) != eslOK) cmdline_failure(argv[0], "Failed to parse command line: %s\n", go->errbuf); + if (esl_opt_VerifyConfig(go) != eslOK) cmdline_failure(argv[0], "Error in configuration: %s\n", go->errbuf); + if (esl_opt_GetBoolean(go, "-h") ) cmdline_help (argv[0], go); + if (esl_opt_ArgNumber(go) != 2) cmdline_failure(argv[0], "Incorrect number of command line arguments.\n"); + + hmmfile = esl_opt_GetArg(go, 1); + seqfile = esl_opt_GetArg(go, 2); + + if (strcmp(hmmfile, "-") == 0 && strcmp(seqfile, "-") == 0) + cmdline_failure(argv[0], "Either or may be '-' (to read from stdin), but not both.\n"); + + msaopts |= p7_ALL_CONSENSUS_COLS; /* default as of 3.1 */ + if (esl_opt_GetBoolean(go, "--trim")) msaopts |= p7_TRIM; + + /* If caller declared an input format, decode it + */ + if (esl_opt_IsOn(go, "--informat")) { + infmt = esl_sqio_EncodeFormat(esl_opt_GetString(go, "--informat")); + if (infmt == eslSQFILE_UNKNOWN) cmdline_failure(argv[0], "%s is not a recognized input sequence file format\n", esl_opt_GetString(go, "--informat")); + } + + /* Determine output alignment file format */ + outfmt = esl_msafile_EncodeFormat(esl_opt_GetString(go, "--outformat")); + if (outfmt == eslMSAFILE_UNKNOWN) cmdline_failure(argv[0], "%s is not a recognized output MSA file format\n", esl_opt_GetString(go, "--outformat")); + + /* Open output stream */ + if ( (outfile = esl_opt_GetString(go, "-o")) != NULL) + { + if ((ofp = fopen(outfile, "w")) == NULL) + cmdline_failure(argv[0], "failed to open -o output file %s for writing\n", outfile); + } + + + /* If caller forced an alphabet on us, create the one the caller wants + */ + if (esl_opt_GetBoolean(go, "--amino")) abc = esl_alphabet_Create(eslAMINO); + else if (esl_opt_GetBoolean(go, "--dna")) abc = esl_alphabet_Create(eslDNA); + else if (esl_opt_GetBoolean(go, "--rna")) abc = esl_alphabet_Create(eslRNA); + + /* Read one HMM, and make sure there's only one. + */ + status = p7_hmmfile_OpenE(hmmfile, NULL, &hfp, errbuf); + if (status == eslENOTFOUND) p7_Fail("File existence/permissions problem in trying to open HMM file %s.\n%s\n", hmmfile, errbuf); + else if (status == eslEFORMAT) p7_Fail("File format problem in trying to open HMM file %s.\n%s\n", hmmfile, errbuf); + else if (status != eslOK) p7_Fail("Unexpected error %d in opening HMM file %s.\n%s\n", status, hmmfile, errbuf); + + status = p7_hmmfile_Read(hfp, &abc, &hmm); + if (status == eslEFORMAT) p7_Fail("Bad file format in HMM file %s:\n%s\n", hfp->fname, hfp->errbuf); + else if (status == eslEINCOMPAT) p7_Fail("HMM in %s is not in the expected %s alphabet\n", hfp->fname, esl_abc_DecodeType(abc->type)); + else if (status == eslEOF) p7_Fail("Empty HMM file %s? No HMM data found.\n", hfp->fname); + else if (status != eslOK) p7_Fail("Unexpected error in reading HMMs from %s\n", hfp->fname); + + status = p7_hmmfile_Read(hfp, &abc, NULL); + if (status != eslEOF) p7_Fail("HMM file %s does not contain just one HMM\n", hfp->fname); + p7_hmmfile_Close(hfp); + + + /* We're going to build up two arrays: sequences and traces. + * If --mapali option is chosen, the first set of sequences/traces is from the provided alignment + */ + if ( (mapfile = esl_opt_GetString(go, "--mapali")) != NULL) + { + map_alignment(mapfile, hmm, &sq, &tr, &mapseq); + } + totseq = mapseq; + + /* Read digital sequences into an array (possibly concat'ed onto mapped seqs) + */ + status = esl_sqfile_OpenDigital(abc, seqfile, infmt, NULL, &sqfp); + if (status == eslENOTFOUND) p7_Fail("Failed to open sequence file %s for reading\n", seqfile); + else if (status == eslEFORMAT) p7_Fail("Sequence file %s is empty or misformatted\n", seqfile); + else if (status != eslOK) p7_Fail("Unexpected error %d opening sequence file %s\n", status, seqfile); + + ESL_RALLOC(sq, p, sizeof(ESL_SQ *) * (totseq + 1)); + sq[totseq] = esl_sq_CreateDigital(abc); + nseq = 0; + while ((status = esl_sqio_Read(sqfp, sq[totseq+nseq])) == eslOK) + { + nseq++; + ESL_RALLOC(sq, p, sizeof(ESL_SQ *) * (totseq+nseq+1)); + sq[totseq+nseq] = esl_sq_CreateDigital(abc); + } + if (status == eslEFORMAT) esl_fatal("Parse failed (sequence file %s):\n%s\n", + sqfp->filename, esl_sqfile_GetErrorBuf(sqfp)); + else if (status != eslEOF) esl_fatal("Unexpected error %d reading sequence file %s", status, sqfp->filename); + esl_sqfile_Close(sqfp); + totseq += nseq; + + + /* Remaining initializations, including trace array allocation + */ + ESL_RALLOC(tr, p, sizeof(P7_TRACE *) * totseq); + for (idx = mapseq; idx < totseq; idx++) + tr[idx] = p7_trace_CreateWithPP(); + + p7_tracealign_computeTraces(hmm, sq, mapseq, totseq - mapseq, tr); + + p7_tracealign_Seqs(sq, tr, totseq, hmm->M, msaopts, hmm, &msa); + + esl_msafile_Write(ofp, msa, outfmt); + + for (idx = 0; idx <= totseq; idx++) esl_sq_Destroy(sq[idx]); /* including sq[nseq] because we overallocated */ + for (idx = 0; idx < totseq; idx++) p7_trace_Destroy(tr[idx]); + free(sq); + free(tr); + esl_msa_Destroy(msa); + p7_hmm_Destroy(hmm); + if (ofp != stdout) fclose(ofp); + esl_alphabet_Destroy(abc); + esl_getopts_Destroy(go); + return eslOK; + + ERROR: + return status; +} + + + +/***************************************************************** + * Internal functions used by main and API + *****************************************************************/ + +static int +map_alignment(const char *msafile, const P7_HMM *hmm, ESL_SQ ***ret_sq, P7_TRACE ***ret_tr, int *ret_ntot) +{ + ESL_SQ **sq = NULL; + P7_TRACE **tr = NULL; + ESL_MSAFILE *afp = NULL; + ESL_MSA *msa = NULL; + ESL_ALPHABET *abc = (ESL_ALPHABET *) hmm->abc; /* removing const'ness to make compiler happy. Safe. */ + int *matassign = NULL; + uint32_t chksum = 0; + int i,k; + int status; + + status = esl_msafile_Open(&abc, msafile, NULL, eslMSAFILE_UNKNOWN, NULL, &afp); + if (status != eslOK) esl_msafile_OpenFailure(afp, status); + + status = esl_msafile_Read(afp, &msa); + if (status != eslOK) esl_msafile_ReadFailure(afp, status); + + if (! (hmm->flags & p7H_CHKSUM) ) esl_fatal("HMM has no checksum. --mapali unreliable without it."); + if (! (hmm->flags & p7H_MAP) ) esl_fatal("HMM has no map. --mapali can't work without it."); + esl_msa_Checksum(msa, &chksum); + if (hmm->checksum != chksum) esl_fatal("--mapali MSA %s isn't same as the one HMM came from (checksum mismatch)", msafile); + + ESL_ALLOC(sq, sizeof(ESL_SQ *) * msa->nseq); + ESL_ALLOC(tr, sizeof(P7_TRACE *) * msa->nseq); + ESL_ALLOC(matassign, sizeof(int) * (msa->alen + 1)); + + esl_vec_ISet(matassign, msa->alen+1, 0); + for (k = 1; k <= hmm->M; k++) matassign[hmm->map[k]] = 1; + + p7_trace_FauxFromMSA(msa, matassign, p7_DEFAULT, tr); + + /* The 'faux' core traces constructed by FauxFromMSA() may contain + * D->I and I->D transitions. They may *only* now be passed to + * p7_tracealign_Seqs(), which can deal with these 'illegal' + * transitions, in order to exactly reproduce the input --mapali + * alignment. + */ + + for (i = 0; i < msa->nseq; i++) + esl_sq_FetchFromMSA(msa, i, &(sq[i])); + + *ret_ntot = msa->nseq; + *ret_tr = tr; + *ret_sq = sq; + + esl_msafile_Close(afp); + esl_msa_Destroy(msa); + free(matassign); + return eslOK; + + ERROR: + *ret_ntot = 0; + *ret_tr = NULL; + *ret_sq = NULL; + if (afp != NULL) esl_msafile_Close(afp); + if (msa != NULL) esl_msa_Destroy(msa); + if (matassign != NULL) free(matassign); + return status; +} + diff --git a/src/bathbuild.c b/src/bathbuild.c index 7aa5ae5f..2228eef0 100644 --- a/src/bathbuild.c +++ b/src/bathbuild.c @@ -67,6 +67,7 @@ static ESL_OPTIONS options[] = { { "-o", eslARG_OUTFILE,FALSE, NULL, NULL, NULL, NULL, NULL, "direct summary output to file , not stdout", 1 }, { "-O", eslARG_OUTFILE,FALSE, NULL, NULL, NULL, NULL, NULL, "resave annotated, possibly modified MSA to file ", 1 }, { "--ct", eslARG_INT, "1", NULL, NULL, NULL, NULL, NULL, "use alt genetic code of NCBI transl table ", 1 }, + { "--fs", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "calculate FS3/FS5 stats, for bathsearch --fs/--fsonly", 1 }, /* Alternate model construction strategies */ { "--fast", eslARG_NONE, "default",NULL, NULL, CONOPTS, NULL, NULL, "assign cols w/ >= symfrac residues as consensus", 3 }, @@ -264,6 +265,7 @@ output_header(const ESL_GETOPTS *go, const struct cfg_s *cfg) if (fprintf(cfg->ofp, "# input file: %s\n", cfg->infile) < 0) ESL_EXCEPTION_SYS(eslEWRITE, "write failed"); if (fprintf(cfg->ofp, "# output HMM file: %s\n", cfg->hmmfile) < 0) ESL_EXCEPTION_SYS(eslEWRITE, "write failed"); + if (fprintf(cfg->ofp, "# frameshift stats calculated: %s\n", (esl_opt_GetBoolean(go, "--fs") ? "YES" : "NO")) < 0) ESL_EXCEPTION_SYS(eslEWRITE, "write failed"); if (esl_opt_IsUsed(go, "-n") && fprintf(cfg->ofp, "# name (the single) HMM: %s\n", esl_opt_GetString(go, "-n")) < 0) ESL_EXCEPTION_SYS(eslEWRITE, "write failed"); if (esl_opt_IsUsed(go, "-o") && fprintf(cfg->ofp, "# output directed to file: %s\n", esl_opt_GetString(go, "-o")) < 0) ESL_EXCEPTION_SYS(eslEWRITE, "write failed"); @@ -608,6 +610,8 @@ usual_master(const ESL_GETOPTS *go, struct cfg_s *cfg) /* special arguments for hmmbuild */ info[i].bld->fsprob = p7P_FSPROB; + //frameshift stats are opt-in, as in bathsearch; bathconvert adds them later if needed + info[i].bld->fs = esl_opt_GetBoolean(go, "--fs"); info[i].bld->w_len = (go != NULL && esl_opt_IsOn (go, "--w_length")) ? esl_opt_GetInteger(go, "--w_length"): -1; info[i].bld->w_beta = (go != NULL && esl_opt_IsOn (go, "--w_beta")) ? esl_opt_GetReal (go, "--w_beta") : p7_DEFAULT_WINDOW_BETA; if ( info[i].bld->w_beta < 0 || info[i].bld->w_beta > 1 ) esl_fatal("Invalid window-length beta value\n"); diff --git a/src/bathsearch.c b/src/bathsearch.c index 12224f04..fd0cf53d 100644 --- a/src/bathsearch.c +++ b/src/bathsearch.c @@ -745,9 +745,10 @@ serial_master(ESL_GETOPTS *go, struct cfg_s *cfg) om = NULL; /* optimized query profile */ if(esl_opt_IsUsed(go, "--fs") || esl_opt_IsUsed(go, "--fsonly")) { //check that HMM is properly formated for bathsearch - if(!(hmm->fsprob && hmm->ct)) p7_Fail("HMM file %s not formated for this version bathsearch. Please run 'bathconvert new_file.bhmm old_file.bhmm'.\n", cfg->queryfile); - if( hmm->evparam[p7_FTAUFS3] == p7_EVPARAM_UNSET ) p7_Fail("HMM file %s not formated for this version bathsearch. Please run 'bathconvert new_file.bhmm old_file.bhmm'.\n", cfg->queryfile); - if( hmm->evparam[p7_FTAUFS5] == p7_EVPARAM_UNSET ) p7_Fail("HMM file %s not formated for this version bathsearch. Please run 'bathconvert new_file.bhmm old_file.bhmm'.\n", cfg->queryfile); + if( !(hmm->fsprob && hmm->ct) || + hmm->evparam[p7_FTAUFS3] == p7_EVPARAM_UNSET || + hmm->evparam[p7_FTAUFS5] == p7_EVPARAM_UNSET ) + p7_Fail("HMM file %s has no frameshift statistics, which --fs/--fsonly require.\nRebuild with 'bathbuild --fs', or add them with 'bathconvert new_file.bhmm %s'.\n", cfg->queryfile, cfg->queryfile); } else { diff --git a/testsuite/i9-optional-annotation.pl b/testsuite/i9-optional-annotation.pl index 46940aac..bb0f0366 100755 --- a/testsuite/i9-optional-annotation.pl +++ b/testsuite/i9-optional-annotation.pl @@ -58,7 +58,7 @@ BEGIN close ALI1; close SEQ1; -@output = `$builddir/src/bathbuild $tmppfx.bhmm $tmppfx.sto 2>&1`; +@output = `$builddir/src/bathbuild --fs $tmppfx.bhmm $tmppfx.sto 2>&1`; if ($? != 0) { die "FAIL: bathbuild failed\n"; } @output = `$builddir/src/bathsearch --fs --tblout $tmppfx.tbl $tmppfx.bhmm $tmppfx.seq 2>&1`;