+++ /dev/null
-/* initfa.c */
-
-/* $Name: fa_34_26_5 $ - $Id: initfa.c,v 1.148 2007/04/26 18:40:58 wrp Exp $ */
-
-/* copyright (c) 1996, 1997, 1998 William R. Pearson and the U. of Virginia */
-
-/* init??.c files provide function specific initializations */
-
-/* h_init() - called from comp_lib.c, comp_thr.c to initialize pstruct ppst
- which includes the alphabet, and pam matrix
-
- alloc_pam() - allocate pam matrix space
- init_pam2() - convert from 1D to 2D pam
-
- init_pamx() - convert from 1D to 2D pam
-
- f_initenv() - set up mngmsg and pstruct defaults
- f_getopt() - read fasta specific command line options
- f_getarg() - read ktup
-
- resetp() - reset the parameters, scoring matrix for DNA-DNA/DNA-prot
-
- query_parm() - ask for ktup
- last_init() - some things must be done last
-
- f_initpam() - set some parameters based on the pam matrix
-*/
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <ctype.h>
-#include <string.h>
-#include <math.h>
-
-#ifdef UNIX
-#include <sys/types.h>
-#include <sys/stat.h>
-#endif
-
-#include "defs.h"
-#include "structs.h"
-#include "param.h"
-
-#ifndef PCOMPLIB
-#include "mw.h"
-#else
-#include "p_mw.h"
-#endif
-
-#define XTERNAL
-#include "upam.h"
-#include "uascii.h"
-#undef XTERNAL
-
-#define MAXWINDOW 32
-
-int initpam(char *, struct pstruct *);
-void init_pam2 (struct pstruct *ppst);
-void extend_pssm(unsigned char *aa0, int n0, struct pstruct *ppst);
-void build_xascii(int *qascii, char *save_str);
-void ann_ascii(int *qascii, char *ann_arr);
-void re_ascii(int *qascii, int *pascii);
-extern int nrand(int);
-
-/* at some point, all the defaults should be driven from this table */
-/*
-#pgm q_seq l_seq p_seq matrix g_open g_ext fr_shft e_cut ktup
-# -n/-p -s -e -f -h/-j -E argv[3]
-fasta prot(0) prot(0) prot(0) bl50 -10 -2 - 10.0 2
-fasta dna(1) dna(1) dna(1) +5/-4 -14 -4 - 2.0 6
-ssearch prot(0) prot(0) prot(0) bl50 -10 -2 - 10.0 -
-ssearch dna(1) dna(1) dna(1) +5/-4 -14 -4 - 2.0 -
-fastx dna(1) prot(0) prot(0) BL50 -12 -2 -20 5.0 2
-fasty dna(1) prot(0) prot(0) BL50 -12 -2 -20/-24 5.0 2
-tfastx dna(1) prot(0) prot(0) BL50 -14 -2 -20 5.0 2
-tfasty dna(1) prot(0) prot(0) BL50 -14 -2 -20/-24 5.0 2
-fasts prot(0) prot(0) prot(0) MD20-MS - - - 5.0 -
-fasts dna(1) dna(1) dna(1) +2/-4 - - - 5.0 1
-tfasts prot(0) dna(1) prot(0) MD10-MS - - - 2.0 1
-fastf prot(0) prot(0) prot(0) MD20 - - - 2.0 1
-tfastf prot(0) dna(1) prot(0) MD10 - - - 1.0 1
-fastm prot(0) prot(0) prot(0) MD20 - - - 5.0 1
-fastm dna(1) dna(1) dna(1) +2/-4 - - - 2.0 1
-tfastm prot(0) dna(1) prot(0) MD10 - - - 2.0 1
-*/
-
-struct pgm_def_str {
- int pgm_id;
- char *prog_func;
- char *pgm_abbr;
- char *iprompt0;
- char *ref_str;
- int PgmDID;
- char *smstr;
- int g_open_mod;
- int gshift;
- int hshift;
- int e_cut;
- int ktup;
-};
-
-char *ref_str_a[]={
- "\nPlease cite:\n W.R. Pearson & D.J. Lipman PNAS (1988) 85:2444-2448\n",
- "\nPlease cite:\n T. F. Smith and M. S. Waterman, (1981) J. Mol. Biol. 147:195-197; \n W.R. Pearson (1991) Genomics 11:635-650\n",
- "\nPlease cite:\n Pearson et al, Genomics (1997) 46:24-36\n",
- "\nPlease cite:\n Mackey et al. Mol. Cell. Proteomics (2002) 1:139-147\n",
- "\nPlease cite:\n W.R. Pearson (1996) Meth. Enzymol. 266:227-258\n"
-};
-
-#define FA_PID 1
-#define SS_PID 2
-#define FX_PID 3
-#define FY_PID 4
-#define FS_PID 5
-#define FF_PID 6
-#define FM_PID 7
-#define RSS_PID 8
-#define RFX_PID 9
-#define SSS_PID 10 /* old (slow) non-PG Smith-Waterman */
-#define TFA_PID FA_PID+10
-#define TFX_PID FX_PID+10
-#define TFY_PID FY_PID+10
-#define TFS_PID FS_PID+10
-#define TFF_PID FF_PID+10
-#define TFM_PID FM_PID+10
-
-struct pgm_def_str
-pgm_def_arr[20] = {
- {0, "", "", "", NULL, 400, "", 0, 0, 0, 1.0, 0 }, /* 0 */
- {FA_PID, "FASTA", "fa",
- "FASTA searches a protein or DNA sequence data bank",
- NULL, 401, "BL50", 0, 0, 0, 10.0, 2}, /* 1 - FASTA */
- {SS_PID, "SSEARCH","gsw","SSEARCH searches a sequence data bank",
- NULL, 404, "BL50", 0, 0, 0, 10.0, 0}, /* 2 - SSEARCH */
- {FX_PID, "FASTX","fx",
- "FASTX compares a DNA sequence to a protein sequence data bank",
- NULL, 405, "BL50", -2, -20, 0, 5.0, 2}, /* 3 - FASTX */
- {FY_PID, "FASTY", "fy",
- "FASTY compares a DNA sequence to a protein sequence data bank",
- NULL, 405, "BL50", -2, -20, -24, 5.0, 2}, /* 4 - FASTY */
- {FS_PID, "FASTS", "fs",
- "FASTS compares linked peptides to a protein data bank",
- NULL, 400, "MD20-MS", 0, 0, 0, 5.0, 1}, /* 5 - FASTS */
- {FF_PID, "FASTF", "ff",
- "FASTF compares mixed peptides to a protein databank",
- NULL, 400, "MD20", 0, 0, 0, 2.0, 1 }, /* 6 - FASTF */
- {FM_PID, "FASTM", "fm",
- "FASTM compares ordered peptides to a protein data bank",
- NULL, 400, "MD20", 0, 0, 0, 5.0, 1 }, /* 7 - FASTM */
- {RSS_PID, "PRSS", "rss",
- "PRSS evaluates statistical signficance using Smith-Waterman",
- NULL, 401, "BL50", 0, 0, 0, 1000.0, 0 }, /* 8 - PRSS */
- {RFX_PID,"PRFX", "rfx",
- "PRFX evaluates statistical signficance using FASTX",
- NULL, 401, "BL50", -2, -20, -24, 1000.0, 2 }, /* 9 - PRFX */
- {SSS_PID, "OSEARCH","ssw","OSEARCH searches a sequence data bank",
- NULL, 404, "BL50", 0, 0, 0, 10.0, 0}, /* 2 - OSEARCH */
- {TFA_PID, "TFASTA", "tfa",
- "TFASTA compares a protein to a translated DNA data bank",
- NULL, 402, "BL50", -2, 0, 0, 5.0, 2 },
- {0, "", "", "", NULL, 400, "", 0, 0, 0, 1.0, 0 }, /* 0 */
- {TFX_PID, "TFASTX", "tfx",
- "TFASTX compares a protein to a translated DNA data bank",
- NULL, 406, "BL50", -2, -20, 0, 2.0, 2},
- {TFY_PID, "TFASTY", "tfy",
- "TFASTY compares a protein to a translated DNA data bank",
- NULL, 406, "BL50", -2, -20, -24, 2.0, 2},
- {TFS_PID, "TFASTS", "tfs",
- "TFASTS compares linked peptides to a translated DNA data bank",
- NULL, 400, "MD10-MS", 0, 0, 0, 2.0, 2 },
- {TFF_PID, "TFASTF", "tff",
- "TFASTF compares mixed peptides to a protein databank",
- NULL, 400, "MD10", 0, 0, 0, 1.0, 1 },
- {TFM_PID, "TFASTM", "tfm",
- "TFASTM compares ordered peptides to a translated DNA databank",
- NULL, 400, "MD10", 0, 0, 0, 1.0, 1 }
-};
-
-struct msg_def_str {
- int pgm_id;
- int q_seqt;
- int l_seqt;
- int p_seqt;
- int sw_flag;
- int stages;
- int qframe;
- int nframe;
- int nrelv, srelv, arelv;
- char *f_id0, *f_id1, *label;
-};
-
-/* pgm_id q_seqt l_seqt p_seqt sw_f st qf nf nrv srv arv s_ix */
-struct msg_def_str msg_def_arr[20] = {
- {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, "", "", ""}, /* ID=0 */
- {FA_PID, SEQT_UNK, SEQT_PROT, SEQT_PROT, 1, 1, 1, -1, 3, 1, 3,
- "fa","sw", "opt"},
- {SS_PID, SEQT_UNK, SEQT_PROT, SEQT_PROT, 1, 1, 1, -1, 1, 1, 1,
- "sw","sw", "s-w"},
- {FX_PID, SEQT_DNA, SEQT_PROT, SEQT_PROT, 1, 1, 2, -1, 3, 1, 3,
- "fx","sx", "opt"},
- {FY_PID, SEQT_DNA, SEQT_PROT, SEQT_PROT, 1, 1, 2, -1, 3, 1, 3,
- "fy","sy", "opt"},
- {FS_PID, SEQT_UNK, SEQT_PROT, SEQT_PROT, 1, 1, 1, -1, 3, 2, 3,
- "fs","fs", "initn init1"},
- {FF_PID, SEQT_PROT,SEQT_PROT, SEQT_PROT, 1, 1, 1, -1, 3, 2, 3,
- "ff","ff", "initn init1"},
- {FM_PID, SEQT_PROT,SEQT_PROT, SEQT_PROT, 1, 1, 1, -1, 3, 2, 3,
- "fm","fm","initn init1"},
- {RSS_PID, SEQT_UNK,SEQT_PROT, SEQT_PROT, 0, 1, 1, -1, 1, 1, 1,
- "rss","sw","s-w"},
- {RFX_PID, SEQT_DNA,SEQT_PROT, SEQT_PROT, 0, 1, 2, -1, 3, 1, 3,
- "rfx","sx","opt"},
- {SSS_PID, SEQT_UNK,SEQT_PROT, SEQT_PROT, 1, 1, 1, -1, 1, 1, 1,
- "sw","sw", "s-w"},
- {TFA_PID, SEQT_PROT,SEQT_DNA, SEQT_PROT, 0, 1, 1, 6, 3, 1, 3,
- "tfa","fa","initn init1"},
- {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, "", "", ""}, /* ID=12 */
- {TFX_PID, SEQT_PROT,SEQT_DNA, SEQT_PROT, 1, 1, 1, 2, 3, 2, 3,
- "tfx","sx","initn opt"},
- {TFY_PID, SEQT_PROT,SEQT_DNA, SEQT_PROT, 1, 1, 1, 2, 3, 2, 3,
- "tfy","sy","initn opt"},
- {TFS_PID, SEQT_PROT,SEQT_DNA, SEQT_PROT, 1, 1, 1, 6, 3, 2, 3,
- "tfs","fs","initn init1"},
- {TFF_PID, SEQT_PROT,SEQT_DNA, SEQT_PROT, 1, 1, 1, 6, 3, 2, 3,
- "tff","ff","initn init1"},
- {TFM_PID, SEQT_PROT,SEQT_DNA, SEQT_PROT, 1, 1, 1, 6, 3, 2, 3,
- "tfm","fm","initn init1"}
-};
-
-int
-get_pgm_id() {
-
- int rval=0;
-
-#ifdef FASTA
-#ifndef TFAST
- pgm_def_arr[FA_PID].ref_str = ref_str_a[0];
- rval=FA_PID;
-#else
- pgm_def_arr[TFA_PID].ref_str = ref_str_a[0];
- rval=TFA_PID;
-#endif
-#endif
-
-#ifdef FASTX
-#ifndef TFAST
-#ifndef PRSS
- pgm_def_arr[FX_PID].ref_str = ref_str_a[2];
- rval=FX_PID;
-#else
- pgm_def_arr[RFX_PID].ref_str = ref_str_a[2];
- rval=RFX_PID;
-#endif
-#else
- pgm_def_arr[TFX_PID].ref_str = ref_str_a[2];
- rval=TFX_PID;
-#endif
-#endif
-
-#ifdef FASTY
-#ifndef TFAST
- pgm_def_arr[FY_PID].ref_str = ref_str_a[2];
- rval=FY_PID;
-#else
- pgm_def_arr[TFY_PID].ref_str = ref_str_a[2];
- rval=TFY_PID;
-#endif
-#endif
-
-#ifdef FASTS
-#ifndef TFAST
- pgm_def_arr[FS_PID].ref_str = ref_str_a[3];
- rval=FS_PID;
-#else
- pgm_def_arr[TFS_PID].ref_str = ref_str_a[3];
- rval=TFS_PID;
-#endif
-#endif
-
-#ifdef FASTF
-#ifndef TFAST
- pgm_def_arr[FF_PID].ref_str = ref_str_a[3];
- rval=FF_PID;
-#else
- pgm_def_arr[TFF_PID].ref_str = ref_str_a[3];
- rval=TFF_PID;
-#endif
-#endif
-
-#ifdef FASTM
-#ifndef TFAST
- pgm_def_arr[FM_PID].ref_str = ref_str_a[3];
- rval=FM_PID;
-#else
- pgm_def_arr[TFM_PID].ref_str = ref_str_a[3];
- rval=TFM_PID;
-#endif
-#endif
-
-#ifdef SSEARCH
- pgm_def_arr[SS_PID].ref_str = ref_str_a[1];
- rval=SS_PID;
-#endif
-
-#ifdef OSEARCH
- pgm_def_arr[SSS_PID].ref_str = ref_str_a[1];
- rval=SSS_PID;
-#endif
-
-#ifdef PRSS
-#ifndef FASTX
- pgm_def_arr[RSS_PID].ref_str = ref_str_a[4];
- rval=RSS_PID;
-#endif
-#endif
-
- return rval;
-}
-
-char *iprompt1=" test sequence file name: ";
-char *iprompt2=" database file name: ";
-
-char *verstr="version 34.26.5 April 26, 2007";
-
-char *s_optstr = "13Ac:f:g:h:j:k:nopP:r:s:St:Ux:y:";
-
-static int mktup=2;
-static int ktup_set = 0;
-static int gap_set=0;
-static int del_set=0;
-static int mshuff_set = 0;
-static int prot2dna = 0;
-
-extern int max_workers;
-
-extern void s_abort(char *, char *);
-extern void init_ascii(int ext_sq, int *sascii, int dnaseq);
-extern int standard_pam(char *smstr, struct pstruct *ppst,
- int del_set, int gap_set);
-extern void mk_n_pam(int *arr,int siz, int mat, int mis);
-extern int karlin(int , int, double *, double *, double *);
-extern void init_karlin_a(struct pstruct *, double *, double **);
-extern int do_karlin_a(int **, struct pstruct *, double *,
- double *, double *, double *, double *);
-
-#if defined(TFAST) || defined(FASTX) || defined(FASTY)
-extern void aainit(int tr_type, int debug);
-#endif
-
-char *iprompt0, *prog_func, *refstr;
-
-
-/* Sets defaults assuming a protein sequence */
-void h_init (struct pstruct *ppst, struct mngmsg *m_msp, char *pgm_abbr)
-{
- struct pgm_def_str pgm_def;
- int i, pgm_id;
-
- ppst->pgm_id = pgm_id = get_pgm_id();
- pgm_def = pgm_def_arr[pgm_id];
-
- /* check that pgm_def_arr[] is valid */
- if (pgm_def.pgm_id != pgm_id) {
- fprintf(stderr,
- "**pgm_def integrity failure: def.pgm_id %d != pgm_id %d**\n",
- pgm_def.pgm_id, pgm_id);
- exit(1);
- }
-
- /* check that msg_def_arr[] is valid */
- if (msg_def_arr[pgm_id].pgm_id != pgm_id) {
- fprintf(stderr,
- "**msg_def integrity failure: def.pgm_id %d != pgm_id %d**\n",
- msg_def_arr[pgm_id].pgm_id, pgm_id);
- exit(1);
- }
-
- strncpy(pgm_abbr,pgm_def.pgm_abbr,MAX_SSTR);
- iprompt0 = pgm_def.iprompt0;
- refstr = pgm_def.ref_str;
- prog_func = pgm_def.prog_func;
-
- /* MAXTOT = MAXTST + MAXLIB for everything except TFAST,
- where it is MAXTST + MAXTRN */
- m_msp->max_tot = MAXTOT;
-
- /* set up DNA query sequence if required*/
- if (msg_def_arr[pgm_id].q_seqt == SEQT_DNA) {
- memcpy(qascii,nascii,sizeof(qascii));
- m_msp->qdnaseq = SEQT_DNA;
- }
- else { /* when SEQT_UNK, start with protein */
- memcpy(qascii,aascii,sizeof(qascii));
- m_msp->qdnaseq = msg_def_arr[pgm_id].q_seqt;
- }
-
-#if defined(FASTF) || defined(FASTS) || defined(FASTM)
- qascii[','] = ESS;
- /* also initialize aascii, nascii for databases */
- qascii['*'] = NA;
-#endif
-
- /* initialize a pam matrix */
- strncpy(ppst->pamfile,pgm_def.smstr,MAX_FN);
- standard_pam(ppst->pamfile,ppst,del_set,gap_set);
- ppst->have_pam2 = 0;
-
- /* this is always protein by default */
- ppst->nsq = naa;
- ppst->nsqx = naax;
- for (i=0; i<=ppst->nsqx; i++) {
- ppst->sq[i] = aa[i];
- ppst->hsq[i] = haa[i];
- ppst->sqx[i]=aax[i]; /* sq = aa */
- ppst->hsqx[i]=haax[i]; /* hsq = haa */
- }
- ppst->sq[ppst->nsqx+1] = ppst->sqx[ppst->nsqx+1] = '\0';
-
- /* set up the c_nt[] mapping */
-
-#if defined(FASTS) || defined(FASTF) || defined(FASTM)
- ppst->c_nt[ESS] = ESS;
-#endif
- ppst->c_nt[0]=0;
- for (i=1; i<=nnt; i++) {
- ppst->c_nt[i]=gc_nt[i];
- ppst->c_nt[i+nnt]=gc_nt[i]+nnt;
- }
-}
-
-/*
- * alloc_pam(): allocates memory for the 2D pam matrix as well
- * as for the integer array used to transmit the pam matrix
- */
-void
-alloc_pam (int d1, int d2, struct pstruct *ppst)
-{
- int i, *d2p;
- char err_str[128];
-
- if ((ppst->pam2[0] = (int **) malloc (d1 * sizeof (int *))) == NULL) {
- sprintf(err_str,"Cannot allocate 2D pam matrix: %d",d1);
- s_abort (err_str,"");
- }
-
- if ((ppst->pam2[1] = (int **) malloc (d1 * sizeof (int *))) == NULL) {
- sprintf(err_str,"Cannot allocate 2D pam matrix: %d",d1);
- s_abort (err_str,"");
- }
-
- if ((d2p = pam12 = (int *) calloc (d1 * d2, sizeof (int))) == NULL) {
- sprintf(err_str,"Cannot allocate 2D pam matrix: %d",d1);
- s_abort (err_str,"");
- }
-
- for (i = 0; i < d1; i++, d2p += d2)
- ppst->pam2[0][i] = d2p;
-
- if ((d2p=pam12x= (int *) malloc (d1 * d2 * sizeof (int))) == NULL) {
- sprintf(err_str,"Cannot allocate 2d pam matrix: %d",d2);
- s_abort (err_str,"");
- }
-
- for (i = 0; i < d1; i++, d2p += d2)
- ppst->pam2[1][i] = d2p;
-
- ppst->have_pam2 = 1;
-}
-
-/*
- * init_pam2(struct pstruct pst): Converts 1-D pam matrix to 2-D
- */
-void
-init_pam2 (struct pstruct *ppst) {
- int i, j, k, nsq;
-
- nsq = ppst->nsq;
-
- ppst->pam2[0][0][0] = -BIGNUM;
- ppst->pam_h = -1; ppst->pam_l = 1;
-
- k = 0;
- for (i = 1; i <= nsq; i++) {
- ppst->pam2[0][0][i] = ppst->pam2[0][i][0] = -BIGNUM;
- for (j = 1; j <= i; j++) {
- ppst->pam2[0][j][i] = ppst->pam2[0][i][j] = pam[k++] - ppst->pamoff;
- if (ppst->pam_l > ppst->pam2[0][i][j]) ppst->pam_l =ppst->pam2[0][i][j];
- if (ppst->pam_h < ppst->pam2[0][i][j]) ppst->pam_h =ppst->pam2[0][i][j];
- }
- }
-}
-
-void
-init_pamx (struct pstruct *ppst) {
- int i, j, k, nsq, pam_xx, pam_xm;
- int sa_x, sa_t, tmp;
-
- nsq = ppst->nsq;
-
- ppst->nt_align = (ppst->dnaseq== SEQT_DNA || ppst->dnaseq == SEQT_RNA);
-
- if (ppst->nt_align) {
- sa_x = pascii['N'];
- sa_t = sa_x;
- }
- else {
- sa_x = pascii['X'];
- sa_t = pascii['*'];
- }
-
- if (ppst->dnaseq == SEQT_RNA) {
- tmp = ppst->pam2[0][nascii['G']][nascii['G']] - 1;
- ppst->pam2[0][nascii['A']][nascii['G']] =
- ppst->pam2[0][nascii['C']][nascii['T']] =
- ppst->pam2[0][nascii['C']][nascii['U']] = tmp;
- }
-
- if (ppst->pam_x_set) {
- for (i=1; i<=nsq; i++) {
- ppst->pam2[0][sa_x][i] = ppst->pam2[0][i][sa_x]=ppst->pam_xm;
- ppst->pam2[0][sa_t][i] = ppst->pam2[0][i][sa_t]=ppst->pam_xm;
- }
- ppst->pam2[0][sa_x][sa_x]=ppst->pam_xx;
- ppst->pam2[0][sa_t][sa_t]=ppst->pam_xx;
- }
- else {
- ppst->pam_xx = ppst->pam2[0][sa_x][sa_x];
- ppst->pam_xm = ppst->pam2[0][1][sa_x];
- }
-
- pam_xx = ppst->pam_xx;
- pam_xm = ppst->pam_xm;
-
- if (ppst->ext_sq_set) { /* using extended alphabet */
- /* fill in pam2[1] matrix */
- ppst->pam2[1][0][0] = -BIGNUM;
- /* fill in additional parts of the matrix */
- for (i = 1; i <= nsq; i++) {
-
- /* -BIGNUM to all matches vs 0 */
- ppst->pam2[0][0][i+nsq] = ppst->pam2[0][i+nsq][0] =
- ppst->pam2[1][0][i+nsq] = ppst->pam2[1][i+nsq][0] =
- ppst->pam2[1][0][i] = ppst->pam2[1][i][0] = -BIGNUM;
-
- for (j = 1; j <= nsq; j++) {
-
- /* replicate pam2[0] to i+nsq, j+nsq */
- ppst->pam2[0][i+nsq][j] = ppst->pam2[0][i][j+nsq] =
- ppst->pam2[0][i+nsq][j+nsq] = ppst->pam2[1][i][j] =
- ppst->pam2[0][i][j];
-
- /* set the high portion of pam2[1] to the corresponding value
- of pam2[1][sa_x][j] */
-
- ppst->pam2[1][i+nsq][j] = ppst->pam2[1][i][j+nsq]=
- ppst->pam2[1][i+nsq][j+nsq]=ppst->pam2[0][sa_x][j];
- }
- }
- }
-}
-
-/* function specific initializations */
-void
-f_initenv (struct mngmsg *m_msp, struct pstruct *ppst, unsigned char **aa0) {
- struct msg_def_str m_msg_def;
- int pgm_id;
-
- pgm_id = ppst->pgm_id;
- m_msg_def = msg_def_arr[pgm_id];
-
- m_msp->last_calc_flg=0;
-
- strncpy(m_msp->f_id0,m_msg_def.f_id0,sizeof(m_msp->f_id0));
- strncpy(m_msp->f_id1,m_msg_def.f_id1,sizeof(m_msp->f_id1));
- strncpy (m_msp->label, m_msg_def.label, sizeof(m_msp->label));
-
-#ifndef SSEARCH
- strncpy (m_msp->alab[0],"initn",20);
- strncpy (m_msp->alab[1],"init1",20);
- strncpy (m_msp->alab[2],"opt",20);
-#else
- strncpy (m_msp->alab[0],"s-w opt",20);
-#endif
-
- ppst->gdelval += pgm_def_arr[pgm_id].g_open_mod;
- ppst->sw_flag = m_msg_def.sw_flag;
- m_msp->e_cut=pgm_def_arr[pgm_id].e_cut;
-
- ppst->score_ix = 0;
- ppst->histint = 2;
- m_msp->qframe = m_msg_def.qframe;
- ppst->sw_flag = m_msg_def.sw_flag;
- m_msp->nframe = m_msg_def.nframe;
- m_msp->nrelv = m_msg_def.nrelv;
- m_msp->srelv = m_msg_def.srelv;
- m_msp->arelv = m_msg_def.arelv;
- m_msp->stages = m_msg_def.stages;
-#if defined(PRSS)
- m_msp->shuff_wid = 0;
- m_msp->shuff_max = 200;
-#endif
-
- /* see param.h for the definition of all these */
-
- m_msp->qshuffle = 0;
- m_msp->nm0 = 1;
- m_msp->escore_flg = 0;
-
- /* pam information */
- ppst->pam_pssm = 0;
-#if defined(FASTS) || defined(FASTF) || defined(FASTM)
- ppst->pam_xx = ppst->pam_xm = 0;
-#else
- ppst->pam_xx = 1; /* set >0 to use pam['X']['X'] value */
- ppst->pam_xm = -1; /* set >0 to use pam['X']['A-Z'] value */
-#endif
- ppst->pam_x_set = 0;
- ppst->pam_set = 0;
- ppst->pam_pssm = 0;
- ppst->p_d_set = 0;
- ppst->pamoff = 0;
- ppst->ext_sq_set = 0;
-
- if (pgm_def_arr[ppst->pgm_id].ktup > 0) {
- mktup = 2;
- ppst->param_u.fa.bestscale = 300;
- ppst->param_u.fa.bestoff = 36;
- ppst->param_u.fa.bkfact = 6;
- ppst->param_u.fa.scfact = 3;
- ppst->param_u.fa.bktup = 2;
- ppst->param_u.fa.ktup = 0;
- ppst->param_u.fa.bestmax = 50;
- ppst->param_u.fa.pamfact = 1;
- ppst->param_u.fa.altflag = 0;
- ppst->param_u.fa.optflag = 1;
- ppst->param_u.fa.iniflag = 0;
- ppst->param_u.fa.optcut = 0;
- ppst->param_u.fa.optcut_set = 0;
- ppst->param_u.fa.cgap = 0;
- ppst->param_u.fa.optwid = MAXWINDOW;
- }
-
-}
-
-/* switches for fasta only */
-
-static int shift_set=0;
-static int subs_set=0;
-static int sw_flag_set=0;
-static int nframe_set=0;
-static int wid_set=0;
-
-void
-f_getopt (char copt, char *optarg,
- struct mngmsg *m_msg, struct pstruct *ppst)
-{
- int pgm_id;
- char *bp;
-
- pgm_id = ppst->pgm_id;
-
- switch (copt) {
- case '1':
- if (pgm_def_arr[pgm_id].ktup > 0) {
- ppst->param_u.fa.iniflag=1;
- }
- break;
- case '3':
- nframe_set = 1;
- if (pgm_id == TFA_PID) {
- m_msg->nframe = 3; break;
- }
- else {
- m_msg->nframe = 1; /* for TFASTXY */
- m_msg->qframe = 1; /* for FASTA, FASTX */
- }
- break;
- case 'A':
- ppst->sw_flag= 1;
- sw_flag_set = 1;
- break;
- case 'c':
- if (pgm_def_arr[pgm_id].ktup > 0) {
- sscanf (optarg, "%d", &ppst->param_u.fa.optcut);
- ppst->param_u.fa.optcut_set = 1;
- }
- break;
- case 'f':
- sscanf (optarg, "%d", &ppst->gdelval);
- if (ppst->gdelval > 0) ppst->gdelval = -ppst->gdelval;
- del_set = 1;
- break;
- case 'g':
- sscanf (optarg, "%d", &ppst->ggapval);
- if (ppst->ggapval > 0) ppst->ggapval = -ppst->ggapval;
- gap_set = 1;
- break;
- case 'h':
- sscanf (optarg, "%d", &ppst->gshift);
- if (ppst->gshift > 0) ppst->gshift = -ppst->gshift;
- shift_set = 1;
- break;
- case 'j':
- sscanf (optarg, "%d", &ppst->gsubs);
- subs_set = 1;
- break;
- case 'k':
- sscanf (optarg, "%d", &m_msg->shuff_max);
- mshuff_set = 1;
- break;
- case 'n':
- m_msg->qdnaseq = SEQT_DNA;
- re_ascii(qascii,nascii);
- strncpy(m_msg->sqnam,"nt",4);
- prot2dna = 1;
- break;
- case 'o':
- if (pgm_def_arr[pgm_id].ktup > 0) {
- ppst->param_u.fa.optflag = 0;
- msg_def_arr[pgm_id].nrelv = m_msg->nrelv = 2;
- }
- break;
- case 'p':
- m_msg->qdnaseq = SEQT_PROT;
- ppst->dnaseq = SEQT_PROT;
- strncpy(m_msg->sqnam,"aa",4);
- break;
- case 'P':
- strncpy(ppst->pgpfile,optarg,MAX_FN);
- if ((bp=strchr(ppst->pgpfile,' '))!=NULL) {
- *bp='\0';
- ppst->pgpfile_type = atoi(bp+1);
- }
- else ppst->pgpfile_type = 0;
- ppst->pgpfile[MAX_FN-1]='\0';
- ppst->pam_pssm = 1;
- break;
- case 'r':
- sscanf(optarg,"%d/%d",&ppst->p_d_mat,&ppst->p_d_mis);
- if (ppst->p_d_mat > 0 && ppst->p_d_mis < 0) {
- ppst->p_d_set = 1;
- strncpy(ppst->pamfile,optarg,40);
- }
- break;
- case 's':
- strncpy (ppst->pamfile, optarg, 120);
- ppst->pamfile[120-1]='\0';
- if (!standard_pam(ppst->pamfile,ppst,del_set, gap_set)) {
- initpam (ppst->pamfile, ppst);
- }
- ppst->pam_set=1;
- break;
- case 'S': /* turn on extended alphabet for seg */
- ppst->ext_sq_set = 1;
- break;
- case 't':
- if (tolower(optarg[0])=='t') {
- m_msg->term_code = aascii['*']; optarg++;
- }
- if (*optarg) {sscanf (optarg, "%d", &ppst->tr_type);}
- break;
- case 'U':
- m_msg->qdnaseq = SEQT_RNA;
- memcpy(qascii,nascii,sizeof(qascii));
- strncpy(m_msg->sqnam,"nt",4);
- nt[nascii['T']]='U';
- prot2dna=1;
- break;
- case 'x':
- if (strchr(optarg,',')!=NULL) {
- sscanf (optarg,"%d,%d",&ppst->pam_xx, &ppst->pam_xm);
- }
- else {
- sscanf (optarg,"%d",&ppst->pam_xx);
- ppst->pam_xm = ppst->pam_xx;
- }
- ppst->pam_x_set=1;
- break;
- case 'y':
- if (pgm_def_arr[pgm_id].ktup > 0) {
- sscanf (optarg, "%d", &ppst->param_u.fa.optwid);
- wid_set = 1;
- }
- break;
- }
-}
-
-void
-f_lastenv (struct mngmsg *m_msg, struct pstruct *ppst)
-{
- char save_str[MAX_SSTR];
-
-#if !defined(FASTM) && !defined(FASTS) && !defined(FASTF)
- strncpy(save_str,"*",sizeof(save_str));
-#else
- strncpy(save_str,",",sizeof(save_str));
-#endif
-
- if (m_msg->qdnaseq == SEQT_UNK) {
- build_xascii(qascii,save_str);
- if (m_msg->ann_flg) ann_ascii(qascii,m_msg->ann_arr);
- }
-
-/* this check allows lc DNA sequence queries with FASTX */
-#if defined(FASTA) && !defined(FASTS) && !defined(FASTM) && !defined(FASTF)
- else
- init_ascii(ppst->ext_sq_set,qascii,m_msg->qdnaseq);
-#endif
-}
-
-void
-f_getarg (int argc, char **argv, int optind,
- struct mngmsg *m_msg, struct pstruct *ppst)
-{
-
- if (pgm_def_arr[ppst->pgm_id].ktup > 0) {
- if (argc - optind >= 4) {
- sscanf (argv[optind + 3], "%d", &ppst->param_u.fa.ktup);
- ktup_set = 1;
- }
- else
- ppst->param_u.fa.ktup = -ppst->param_u.fa.bktup;
- }
-
- if (ppst->pgm_id == RSS_PID && argc - optind > 3) {
- sscanf (argv[optind + 3], "%d", &m_msg->shuff_max);
- }
-
- if (ppst->pgm_id == RFX_PID && argc - optind > 4) {
- sscanf (argv[optind + 4], "%d", &m_msg->shuff_max);
- }
-}
-
-/* fills in the query ascii mapping from the parameter
- ascii mapping.
-*/
-
-void
-re_ascii(int *qascii, int *pascii) {
- int i;
-
- for (i=0; i < 128; i++) {
- if (qascii[i] > '@' || qascii[i] < ESS) {
- qascii[i] = pascii[i];
- }
- }
-}
-
-
-/* recode has become function specific to accommodate FASTS/M */
-/* modified 28-Dec-2004 to ensure that all mapped characters
- are valid */
-int
-recode(unsigned char *seq, int n, int *qascii, int nsqx) {
- int i,j;
- char save_c;
-
-#if defined(FASTS) || defined(FASTM)
- qascii[',']=ESS;
-#endif
-
- for (i=0; i < n; i++) {
- save_c = seq[i];
- if (seq[i] > '@') seq[i] = qascii[seq[i]];
- if (seq[i] > nsqx && seq[i]!=ESS) {
- fprintf(stderr, "*** Warning - unrecognized residue at %d:%c - %2d\n",
- i,save_c,save_c);
- seq[i] = qascii['X'];
- }
- }
- seq[i]=EOSEQ;
- return i;
-}
-
-/* here we have the query sequence, all the command line options,
- but we need to set various parameter options based on the type
- of the query sequence (m_msg->qdnaseq = 0:protein/1:DNA) and
- the function (FASTA/FASTX/TFASTA)
-*/
-
-/* this resetp is for conventional a FASTA/TFASTXYZ search */
-void
-resetp (struct mngmsg *m_msg, struct pstruct *ppst) {
- int i, pgm_id;
-
- pgm_id = ppst->pgm_id;
-
-#if defined(TFAST)
- if (m_msg->qdnaseq == SEQT_DNA || m_msg->qdnaseq == SEQT_RNA) {
- fprintf(stderr," %s compares a protein to a translated\n\
-DNA sequence library. Do not use a DNA query/scoring matrix.\n",prog_func);
- exit(1);
- }
-#else
-#if (defined(FASTX) || defined(FASTY))
- if (!(m_msg->qdnaseq == SEQT_DNA || m_msg->qdnaseq == SEQT_RNA)) {
- fprintf(stderr," FASTX/Y compares a DNA sequence to a protein database\n");
- fprintf(stderr," Use a DNA query\n");
- exit(1);
- }
-#endif
-#endif
-
-/* this code changes parameters for programs (FA_PID, SS_PID, FS_PID,
- RSS_PID) that can examine either protein (initial state) or DNA
- Modified May, 2006 to reset e_cut for DNA comparisons.
-*/
-
- if (msg_def_arr[pgm_id].q_seqt == SEQT_UNK) {
- if (m_msg->qdnaseq == SEQT_DNA || m_msg->qdnaseq == SEQT_RNA) {
- msg_def_arr[pgm_id].q_seqt = m_msg->qdnaseq;
- msg_def_arr[pgm_id].p_seqt = SEQT_DNA;
- msg_def_arr[pgm_id].l_seqt = SEQT_DNA;
- if (m_msg->qdnaseq == SEQT_DNA) msg_def_arr[pgm_id].qframe = 2;
- pgm_def_arr[pgm_id].e_cut /= 5.0;
- }
- else {
- msg_def_arr[pgm_id].q_seqt = SEQT_PROT;
- }
- }
-
- ppst->dnaseq = msg_def_arr[pgm_id].p_seqt;
- if (!sw_flag_set) ppst->sw_flag = msg_def_arr[pgm_id].sw_flag;
- if (!m_msg->e_cut_set) m_msg->e_cut=pgm_def_arr[pgm_id].e_cut;
-
- if (ppst->dnaseq == SEQT_DNA && m_msg->qdnaseq==SEQT_RNA) {
- ppst->dnaseq = SEQT_RNA;
- ppst->nt_align = 1;
- }
- if (ppst->dnaseq==SEQT_DNA) pascii = &nascii[0];
- else if (ppst->dnaseq==SEQT_RNA) {
- pascii = &nascii[0];
- ppst->sq[nascii['T']] = 'U';
- }
- else pascii = &aascii[0];
- m_msg->ldnaseq = msg_def_arr[pgm_id].l_seqt;
- if (m_msg->ldnaseq & SEQT_DNA) {
- memcpy(lascii,nascii,sizeof(lascii));
-#ifndef TFAST
-#ifdef DNALIB_LC
- init_ascii(ppst->ext_sq_set,lascii,m_msg->ldnaseq);
-#endif
-#else
- /* no init_ascii() because we translate lower case library sequences */
-#endif
- }
- else {
- memcpy(lascii,aascii,sizeof(lascii)); /* initialize lib mapping */
-
-#if defined(FASTF) || defined(FASTS) || defined(FASTM)
- lascii['*'] = NA;
-#endif
- init_ascii(ppst->ext_sq_set,lascii,m_msg->ldnaseq);
- }
-
- if (!nframe_set) {
- m_msg->qframe = msg_def_arr[pgm_id].qframe;
- m_msg->nframe = msg_def_arr[pgm_id].nframe;
- }
-
- /* the possibilities:
- -i -3 qframe revcomp
- FA_D/FX - - 2 0
- FA_D/FX + - 2 1
- FA_D/FX - + 1 0
- FA_D/FX + + 2 1
- */
-
- if (m_msg->qdnaseq == SEQT_DNA) {
- m_msg->nframe = 1;
- if (m_msg->qframe == 1 && m_msg->revcomp==1) {
- m_msg->qframe = m_msg->revcomp+1;
- }
- }
- else if (m_msg->qdnaseq == SEQT_RNA) {
- m_msg->qframe = m_msg->revcomp+1;
- m_msg->nframe = 1;
- }
-
- /* change settings for DNA search */
- if (ppst->dnaseq == SEQT_DNA || ppst->dnaseq == SEQT_RNA) {
- ppst->histint = 4;
-
- if (!del_set) {
-#ifdef OLD_FASTA_GAP
- ppst->gdelval = -16; /* def. del penalty */
-#else
- ppst->gdelval = -12; /* def. open penalty */
-#endif
- }
- if (!gap_set) ppst->ggapval = -4; /* def. gap penalty */
-
- if (pgm_def_arr[pgm_id].ktup > 0) {
- /* these parameters are used to scale optcut, they should be replaced
- by statistically based parameters */
- if (!wid_set) ppst->param_u.fa.optwid = 16;
- ppst->param_u.fa.bestscale = 80;
- ppst->param_u.fa.bkfact = 5;
- ppst->param_u.fa.scfact = 1;
- ppst->param_u.fa.bktup = 6;
- ppst->param_u.fa.bestmax = 80;
- ppst->param_u.fa.bestoff = 45;
-
- if (!sw_flag_set) {
- ppst->sw_flag = 0;
- strncpy(m_msg->f_id1,"bs",sizeof(m_msg->f_id1));
- }
-
- /* largest ktup */
- mktup = 6;
-
- if (ppst->param_u.fa.pamfact >= 0) ppst->param_u.fa.pamfact = 0;
- if (ppst->param_u.fa.ktup < 0)
- ppst->param_u.fa.ktup = -ppst->param_u.fa.bktup;
- }
-
- ppst->nsq = nnt;
- ppst->nsqx = nntx;
- for (i=0; i<=ppst->nsqx; i++) {
- ppst->hsq[i] = hnt[i];
- ppst->sq[i] = nt[i];
- ppst->hsqx[i] = hntx[i];
- ppst->sqx[i] = ntx[i];
- }
- ppst->sq[ppst->nsqx+1] = ppst->sqx[ppst->nsqx+1] = '\0';
-
- if (!ppst->pam_set) {
- if (ppst->p_d_set)
- mk_n_pam(npam,nnt,ppst->p_d_mat,ppst->p_d_mis);
-#if !defined(FASTS) && !defined(FASTM)
- else if (ppst->pamfile[0]=='\0' || strncmp(ppst->pamfile,"BL50",4)==0) {
- strncpy (ppst->pamfile, "+5/-4", sizeof(ppst->pamfile));
- }
-#else
- else if (strncmp(ppst->pamfile,"MD20",4)==0) {
- strncpy (ppst->pamfile, "+2/-2", sizeof(ppst->pamfile));
- ppst->p_d_mat = +2;
- ppst->p_d_mis = -2;
- mk_n_pam(npam,nnt,ppst->p_d_mat,ppst->p_d_mis);
- }
-#endif
- pam = npam;
- }
-
- strncpy (m_msg->sqnam, "nt",sizeof(m_msg->sqnam));
- strncpy (m_msg->sqtype, "DNA",sizeof(m_msg->sqtype));
- } /* end DNA reset */
-
- else { /* other parameters for protein comparison */
- if (pgm_def_arr[pgm_id].ktup > 0) {
- if (!wid_set) {
- if (ppst->param_u.fa.ktup==1) ppst->param_u.fa.optwid = 32;
- else ppst->param_u.fa.optwid = 16;
- }
- }
- if (!del_set) {ppst->gdelval += pgm_def_arr[pgm_id].g_open_mod;}
- if (!shift_set) {ppst->gshift = pgm_def_arr[pgm_id].gshift;}
- if (!subs_set) {ppst->gsubs = pgm_def_arr[pgm_id].hshift;}
- }
-
-}
-
-/* query_parm() this function asks for any additional parameters
- that have not been provided. Could be null. */
-void
-query_parm (struct mngmsg *m_msp, struct pstruct *ppst)
-{
- char qline[40];
-
- if (pgm_def_arr[ppst->pgm_id].ktup > 0) {
- if (ppst->param_u.fa.ktup < 0)
- ppst->param_u.fa.ktup = -ppst->param_u.fa.ktup;
-
- if (ppst->param_u.fa.ktup == 0) {
- printf (" ktup? (1 to %d) [%d] ", mktup, ppst->param_u.fa.bktup);
- if (fgets (qline, sizeof(qline), stdin) == NULL) exit (0);
- else sscanf(qline,"%d",&ppst->param_u.fa.ktup);
- }
- if (ppst->param_u.fa.ktup == 0)
- ppst->param_u.fa.ktup = ppst->param_u.fa.bktup;
- else ktup_set = 1;
- }
-
-#if defined(PRSS)
- if (m_msp->shuff_max < 10) m_msp->shuff_max = 200;
-
- if (!mshuff_set) {
- printf(" number of shuffles [%d]? ",m_msp->shuff_max);
- fflush(stdout);
- if (fgets (qline, sizeof(qline), stdin) == NULL) exit (0);
- else sscanf(qline,"%d",&m_msp->shuff_max);
- }
-
- if (ppst->zs_win == 0) {
- printf (" local (window) (w) or uniform (u) shuffle [u]? ");
- if (fgets (qline, sizeof(qline), stdin) == NULL) exit (0);
- else if (qline[0]=='w' || qline[0]=='W') {
- m_msp->shuff_wid = 20;
- printf(" local shuffle window size [%d]? ",m_msp->shuff_wid);
- if (fgets (qline, sizeof(qline), stdin) == NULL) exit (0);
- else sscanf(qline,"%d",&m_msp->shuff_wid);
- }
- }
-#endif
-}
-
-/* last_init() cannot look at aa0, n0, because it is only run once,
- it is not run before each new aa0 search */
-void
-last_init (struct mngmsg *m_msg, struct pstruct *ppst
-#ifdef PCOMPLIB
- ,int nnodes
-#endif
- )
-{
- int ix_l, ix_i, i, pgm_id;
- double *kar_p;
- double aa0_f[MAXSQ];
-
- pgm_id = ppst->pgm_id;
-
-#if defined(FASTF) || defined(FASTS) || defined(FASTM)
- m_msg->nohist = 1;
- m_msg->shuff_max = 2000;
-#ifndef PCOMPLIB
- ppst->shuff_node = m_msg->shuff_max/max_workers;
-#else
- ppst->shuff_node = m_msg->shuff_max/nnodes;
-#endif
-#endif
-
- if (m_msg->aln.llen < 1) {
- m_msg->aln.llen = 60;
- }
-
-#ifndef PCOMPLIB
-#if defined(FASTX) || defined(FASTY) || defined(TFAST)
- /* set up translation tables: faatran.c */
- aainit(ppst->tr_type,ppst->debug_lib);
-#endif
-#endif
-
-/* a sanity check */
-#if !defined(TFAST)
- if (m_msg->revcomp && m_msg->qdnaseq!=SEQT_DNA && m_msg->qdnaseq!=SEQT_RNA) {
- fprintf(stderr," cannot reverse complement protein\n");
- m_msg->revcomp = 0;
- }
-#endif
-
- if (pgm_def_arr[pgm_id].ktup > 0) {
-
- if (ppst->param_u.fa.ktup < 0)
- ppst->param_u.fa.ktup = -ppst->param_u.fa.ktup;
-
- if (ppst->param_u.fa.ktup < 1 || ppst->param_u.fa.ktup > mktup) {
- fprintf(stderr," warning ktup = %d out of range [1..%d], reset to %d\n",
- ppst->param_u.fa.ktup, mktup, ppst->param_u.fa.bktup);
- ppst->param_u.fa.ktup = ppst->param_u.fa.bktup;
- }
- }
-
- if (pgm_id == TFA_PID) {
- m_msg->revcomp *= 3;
- if (m_msg->nframe == 3) m_msg->nframe += m_msg->revcomp;
- }
- else if (pgm_id == TFX_PID || pgm_id == TFY_PID) {
- if (m_msg->nframe == 1) m_msg->nframe += m_msg->revcomp;
- }
-
-#if !defined(TFAST)
- /* for fasta/fastx searches, itt iterates the the query strand */
- m_msg->nitt1 = m_msg->qframe-1;
-#else
- /* for tfasta/tfastxy searches, itt iterates the library frames */
- m_msg->nitt1 = m_msg->nframe-1;
-#endif
-
- if (pgm_def_arr[pgm_id].ktup > 0) {
- if (ppst->param_u.fa.ktup>=2 && !wid_set) {
- ppst->param_u.fa.optwid=16;
- switch (pgm_id) {
- case FA_PID:
- m_msg->thr_fact = 32;
- break;
- case FX_PID:
- case FY_PID:
- m_msg->thr_fact = 16;
- break;
- case TFA_PID:
- case TFX_PID:
- case TFY_PID:
- m_msg->thr_fact = 8;
- break;
- default:
- m_msg->thr_fact = 4;
- }
- }
- else { m_msg->thr_fact = 4;}
- }
- else m_msg->thr_fact = 4;
-
-#if defined(PRSS)
- if (m_msg->shuff_max < 10) m_msg->shuff_max = 200;
- if (ppst->zsflag < 10) ppst->zsflag += 10;
- if (ppst->zs_win > 0) {
- m_msg->shuff_wid = ppst->zs_win;
- }
-#endif
-
- if (pgm_def_arr[ppst->pgm_id].ktup > 0) {
- if (ppst->param_u.fa.iniflag) {
- ppst->score_ix = 1;
- strncpy (m_msg->label, "initn init1", sizeof(m_msg->label));
- }
- else if (ppst->param_u.fa.optflag) {
- ppst->score_ix = 2;
- m_msg->stages = 1;
- }
- }
-
- if (!ppst->have_pam2) {
- alloc_pam (MAXSQ, MAXSQ, ppst);
- init_pam2(ppst);
- }
- init_pamx(ppst);
-
- if (ppst->pam_ms) {
- if (m_msg->qdnaseq == SEQT_PROT) {
- /* code to make 'L'/'I' identical scores */
- ix_l = pascii['L'];
- ix_i = pascii['I'];
- ppst->pam2[0][ix_l][ix_i] = ppst->pam2[0][ix_i][ix_l] =
- ppst->pam2[0][ix_l][ix_l] = ppst->pam2[0][ix_i][ix_i] =
- (ppst->pam2[0][ix_l][ix_l]+ppst->pam2[0][ix_i][ix_i]+1)/2;
- for (i=1; i<=ppst->nsq; i++) {
- ppst->pam2[0][i][ix_i] = ppst->pam2[0][i][ix_l] =
- (ppst->pam2[0][i][ix_l]+ppst->pam2[0][i][ix_i]+1)/2;
- ppst->pam2[0][ix_i][i] = ppst->pam2[0][ix_l][i] =
- (ppst->pam2[0][ix_i][i]+ppst->pam2[0][ix_l][i]+1)/2;
- }
-
- /* code to make 'Q'/'K' identical scores */
- if (!shift_set) {
- ix_l = pascii['Q'];
- ix_i = pascii['K'];
- ppst->pam2[0][ix_l][ix_i] = ppst->pam2[0][ix_i][ix_l] =
- ppst->pam2[0][ix_l][ix_l] = ppst->pam2[0][ix_i][ix_i] =
- (ppst->pam2[0][ix_l][ix_l]+ppst->pam2[0][ix_i][ix_i]+1)/2;
- for (i=1; i<=ppst->nsq; i++) {
- ppst->pam2[0][i][ix_i] = ppst->pam2[0][i][ix_l] =
- (ppst->pam2[0][i][ix_l]+ppst->pam2[0][i][ix_i]+1)/2;
- ppst->pam2[0][ix_i][i] = ppst->pam2[0][ix_l][i] =
- (ppst->pam2[0][ix_i][i]+ppst->pam2[0][ix_l][i]+1)/2;
- }
- }
- }
- }
-
- /*
- print_pam(ppst);
- */
-
- /* once we have a complete pam matrix, we can calculate Lambda and K
- for "average" sequences */
- kar_p = NULL;
- init_karlin_a(ppst, aa0_f, &kar_p);
- do_karlin_a(ppst->pam2[0], ppst, aa0_f,
- kar_p, &m_msg->Lambda, &m_msg->K, &m_msg->H);
- free(kar_p);
-
-#if defined(FASTF) || defined(FASTS) || defined(FASTM)
- if (ppst->ext_sq_set) {
- fprintf(stderr," -S not available on [t]fast[fs]\n");
- ppst->ext_sq_set = 0;
-
- /* reset sascii to ignore -S, map lc */
- init_ascii(0,lascii,0);
- }
-#endif
-}
-
-/* this function is left over from the older FASTA format scoring
- matrices that allowed additional parameters (bktup, bkfact) to be
- set in the scoring matrix. It is no longer used. A modern version
- would set parameters based on lambda and K.
-*/
-/*
-void
-f_initpam (line, ppst)
-char *line;
-struct pstruct *ppst;
-{
- if (sscanf (line, " %d %d %d %d %d %d %d", &ppst->param_u.fa.scfact,
- &ppst->param_u.fa.bestoff, &ppst->param_u.fa.bestscale,
- &ppst->param_u.fa.bkfact, &ppst->param_u.fa.bktup,
- &ppst->param_u.fa.bestmax, &ppst->histint) != 7)
- {
- printf (" bestcut parameters - bad format\n");
- exit (1);
- }
-}
-*/
-
-/* alloc_pam2 creates a profile structure */
-int **
-alloc_pam2p(int len, int nsq) {
- int i;
- int **pam2p;
-
- if ((pam2p = (int **)calloc(len,sizeof(int *)))==NULL) {
- fprintf(stderr," Cannot allocate pam2p: %d\n",len);
- return NULL;
- }
-
- if((pam2p[0] = (int *)calloc((nsq+1)*len,sizeof(int)))==NULL) {
- fprintf(stderr, "Cannot allocate pam2p[0]: %d\n", (nsq+1)*len);
- free(pam2p);
- return NULL;
- }
-
- for (i=1; i<len; i++) {
- pam2p[i] = pam2p[0] + (i*(nsq+1));
- }
-
- return pam2p;
-}
-
-void free_pam2p(int **pam2p) {
- if (pam2p) {
- free(pam2p[0]);
- free(pam2p);
- }
-}
-
-/* sortbest has now become comparison function specific so that we can use
- a different comparison for fasts/f
-*/
-#if !defined(FASTS) && !defined (FASTF) && !defined(FASTM)
-#ifndef PCOMPLIB
-void
-qshuffle() {}
-#endif
-
-int
-last_calc(unsigned char *aa0, unsigned char *aa1, int maxn,
- struct beststr **bestp_arr, int nbest,
- struct mngmsg *m_msg, struct pstruct *pst,
- void **f_str, void *rs_str)
-{
- return nbest;
-}
-
-void sortbest (bptr, nbest, irelv)
-struct beststr **bptr;
-int nbest, irelv;
-{
- int gap, i, j;
- struct beststr *tmp;
-
- for (gap = nbest/2; gap > 0; gap /= 2)
- for (i = gap; i < nbest; i++)
- for (j = i - gap; j >= 0; j-= gap) {
- if (bptr[j]->score[irelv] >= bptr[j + gap]->score[irelv]) break;
- tmp = bptr[j];
- bptr[j] = bptr[j + gap];
- bptr[j + gap] = tmp;
- }
-}
-
-void show_aux(FILE *fp, struct beststr *bptr) {}
-void header_aux(FILE *fp) {}
-
-#else
-void sortbest (bptr, nbest, irelv)
-struct beststr **bptr;
-int nbest, irelv;
-{
- int gap, i, j;
- struct beststr *tmp;
-
- for (gap = nbest/2; gap > 0; gap /= 2)
- for (i = gap; i < nbest; i++)
- for (j = i - gap; j >= 0; j-= gap) {
- if (bptr[j]->escore < bptr[j + gap]->escore) break;
- tmp = bptr[j];
- bptr[j] = bptr[j + gap];
- bptr[j + gap] = tmp;
- }
-}
-
-#if defined(FASTS) || defined(FASTM)
-
-#ifndef PCOMPLIB
-/* this shuffle is for FASTS */
-/* convert ',' -> '\0', shuffle each of the substrings */
-void
-qshuffle(unsigned char *aa0, int n0, int nm0) {
-
- unsigned char **aa0start, *aap, tmp;
- int i,j,k, ns;
-
- if ((aa0start=(unsigned char **)calloc(nm0+1,
- sizeof(unsigned char *)))==NULL) {
- fprintf(stderr,"cannot calloc for qshuffle %d\n",nm0);
- exit(1);
- }
-
- aa0start[0]=aa0;
- for (k=1,i=0; i<n0; i++) {
- if (aa0[i]==EOSEQ || aa0[i]==ESS) {
- aa0[i]='\0';
- aa0start[k++] = &aa0[i+1];
- }
- }
-
- /* aa0start has the beginning of each substring */
- for (k=0; k<nm0; k++) {
- aap=aa0start[k];
- ns = strlen((char *)aap);
- for (i=ns; i>1; i--) {
- j = nrand(i);
- tmp = aap[j];
- aap[j] = aap[i-1];
- aap[i-1] = tmp;
- }
- aap[ns] = 0;
- }
-
- for (k=1; k<nm0; k++) {
-/* aap = aa0start[k];
- while (*aap) fputc(pst.sq[*aap++],stderr);
- fputc('\n',stderr);
-*/
- aa0start[k][-1]=ESS;
- }
-
- free(aa0start);
-}
-#endif
-#endif
-
-#ifdef FASTF
-#ifndef PCOMPLIB
-void qshuffle(unsigned char *aa0, int n0, int nm0) {
-
- int i, j, k, nmpos;
- unsigned char tmp;
- int nmoff;
-
- nmoff = (n0 - nm0 - 1)/nm0 + 1;
-
- for (i = nmoff-1 ; i > 0 ; i--) {
-
- /* j = nrand(i); if (i == j) continue;*/ /* shuffle columns */
- j = (nmoff -1 ) - i;
- if (i <= j) break; /* reverse columns */
-
- /* swap all i'th column residues for all j'th column residues */
- for(nmpos = 0, k = 0 ; k < nm0 ; k++, nmpos += nmoff+1 ) {
- tmp = aa0[nmpos + i];
- aa0[nmpos + i] = aa0[nmpos + j];
- aa0[nmpos + j] = tmp;
- }
- }
-}
-#endif
-#endif
-
-
-/* show additional best_str values */
-void show_aux(FILE *fp, struct beststr *bptr) {
- fprintf(fp," %2d %3d",bptr->segnum,bptr->seglen);
-}
-
-void header_aux(FILE *fp) {
- fprintf(fp, " sn sl");
-}
-#endif
-
-void
-fill_pam(int **pam2p, int n0, int nsq, double **freq2d, double scale) {
- int i, j;
- double freq;
-
- /* fprintf(stderr, "scale: %g\n", scale); */
-
- /* now fill in the pam matrix: */
- for (i = 0 ; i < n0 ; i++) {
- for (j = 1 ; j <=20 ; j++) {
- freq = scale * freq2d[i][j-1];
- if ( freq < 0.0) freq -= 0.5;
- else freq += 0.5;
- pam2p[i][j] = (int)(freq);
- }
- }
-}
-
-double
-get_lambda(int **pam2p, int n0, int nsq, unsigned char *query) {
- double lambda, H;
- double *pr, tot, sum;
- int i, ioff, j, min, max;
-
- /* get min and max scores */
- min = BIGNUM;
- max = -BIGNUM;
- if(pam2p[0][1] == -BIGNUM) {
- ioff = 1;
- n0++;
- } else {
- ioff = 0;
- }
-
- for (i = ioff ; i < n0 ; i++) {
- for (j = 1; j <= nsq ; j++) {
- if (min > pam2p[i][j])
- min = pam2p[i][j];
- if (max < pam2p[i][j])
- max = pam2p[i][j];
- }
- }
-
- /* fprintf(stderr, "min: %d\tmax:%d\n", min, max); */
-
- if ((pr = (double *) calloc(max - min + 1, sizeof(double))) == NULL) {
- fprintf(stderr, "Couldn't allocate memory for score probabilities: %d\n", max - min + 1);
- exit(1);
- }
-
- tot = (double) rrtotal * (double) rrtotal * (double) n0;
- for (i = ioff ; i < n0 ; i++) {
- for (j = 1; j <= nsq ; j++) {
- pr[pam2p[i][j] - min] +=
- (double) ((double) rrcounts[aascii[query[i]]] * (double) rrcounts[j]) / tot;
- }
- }
-
- sum = 0.0;
- for(i = 0 ; i <= max-min ; i++) {
- sum += pr[i];
- /* fprintf(stderr, "%3d: %g %g\n", i+min, pr[i], sum); */
- }
- /* fprintf(stderr, "sum: %g\n", sum); */
-
- for(i = 0 ; i <= max-min ; i++) { pr[i] /= sum; }
-
- if (!karlin(min, max, pr, &lambda, &H)) {
- fprintf(stderr, "Karlin lambda estimation failed\n");
- }
-
- /* fprintf(stderr, "lambda: %g\n", lambda); */
- free(pr);
-
- return lambda;
-}
-
-/*
- *aa0 - query sequence
- n0 - length
- pamscale - scaling for pam matrix - provided by apam.c, either
- 0.346574 = ln(2)/2 (P120, BL62) or
- 0.231049 = ln(2)/3 (P250, BL50)
-*/
-
-void
-scale_pssm(int **pssm2p, double **freq2d,
- unsigned char *query, int n0,
- int **pam2, double pamscale);
-
-static unsigned char ustandard_aa[] ="\0ARNDCQEGHILKMFPSTWYV";
-
-void
-read_pssm(unsigned char *aa0, int n0, int nsq,
- double pamscale,
- FILE *fp, int pgpf_type, struct pstruct *ppst) {
- int i, j, len, k;
- int qi, rj; /* qi - index query; rj - index residues (1-20) */
- int **pam2p;
- int first, too_high;
- unsigned char *query, ctmp;
- char dline[512];
- double freq, **freq2d, lambda, new_lambda;
- double scale, scale_high, scale_low;
-
- pam2p = ppst->pam2p[0];
-
- if (pgpf_type == 0) {
-
- if(1 != fread(&len, sizeof(int), 1, fp)) {
- fprintf(stderr, "error reading from checkpoint file: %d\n", len);
- exit(1);
- }
-
- if(len != n0) {
- fprintf(stderr, "profile length (%d) and query length (%d) don't match!\n",
- len,n0);
- exit(1);
- }
-
- /* read over query sequence stored in BLAST profile */
- if(NULL == (query = (unsigned char *) calloc(len+2, sizeof(char)))) {
- fprintf(stderr, "Couldn't allocate memory for query!\n");
- exit(1);
- }
-
- if(len != fread(query, sizeof(char), len, fp)) {
- fprintf(stderr, "Couldn't read query sequence from profile: %s\n", query);
- exit(1);
- }
- }
- else if (pgpf_type == 1) {
-
- if ((fgets(dline,sizeof(dline),fp) == NULL) ||
- (1 != sscanf(dline, "%d",&len))) {
- fprintf(stderr, "error reading from checkpoint file: %d\n", len);
- exit(1);
- }
-
- if(len != n0) {
- fprintf(stderr, "profile length (%d) and query length (%d) don't match!\n",
- len,n0);
- exit(1);
- }
-
- /* read over query sequence stored in BLAST profile */
- if(NULL == (query = (unsigned char *) calloc(len+2, sizeof(char)))) {
- fprintf(stderr, "Couldn't allocate memory for query!\n");
- exit(1);
- }
-
- if (fgets((char *)query,len+2,fp)==NULL) {
- fprintf(stderr, "Couldn't read query sequence from profile: %s\n", query);
- exit(1);
- }
- }
- else {
- fprintf(stderr," Unrecognized PSSM file type: %d\n",pgpf_type);
- exit(1);
- }
-
- /* currently we don't do anything with query; ideally, we should
- check to see that it actually matches aa0 ... */
-
- /* quick 2d array alloc: */
- if((freq2d = (double **) calloc(n0, sizeof(double *))) == NULL) {
- fprintf(stderr, "Couldn't allocate memory for frequencies!\n");
- exit(1);
- }
-
- if((freq2d[0] = (double *) calloc(n0 * 20, sizeof(double))) == NULL) {
- fprintf(stderr, "Couldn't allocate memory for frequencies!\n");
- exit(1);
- }
-
- /* a little pointer arithmetic to fill out 2d array: */
- for (i = 1 ; i < n0 ; i++) {
- freq2d[i] = freq2d[i-1] + 20;
- }
-
- if (pgpf_type == 0) {
- for (qi = 0 ; qi < n0 ; qi++) {
- for (rj = 0 ; rj < 20 ; rj++) {
- if(1 != fread(&freq, sizeof(double), 1, fp)) {
- fprintf(stderr, "Error while reading frequencies!\n");
- exit(1);
- }
- freq2d[qi][rj] = freq;
- }
- }
- }
- else {
- for (qi = 0 ; qi < n0 ; qi++) {
- if ((fgets(dline,sizeof(dline),fp) ==NULL) ||
- (k = sscanf(dline,"%c %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg\n",
- &ctmp, &freq2d[qi][0], &freq2d[qi][1], &freq2d[qi][2], &freq2d[qi][3], &freq2d[qi][4],
- &freq2d[qi][5], &freq2d[qi][6], &freq2d[qi][7], &freq2d[qi][8], &freq2d[qi][9],
- &freq2d[qi][10], &freq2d[qi][11], &freq2d[qi][12], &freq2d[qi][13], &freq2d[qi][14],
- &freq2d[qi][15], &freq2d[qi][16], &freq2d[qi][17], &freq2d[qi][18], &freq2d[qi][19]))<1) {
- fprintf(stderr, "Error while reading frequencies: %d read!\n",k);
- exit(1);
- }
- for (rj=0; rj<20; rj++) { freq2d[qi][rj] /= 10.0; } /* reverse scaling */
- }
- }
-
- scale_pssm(ppst->pam2p[0], freq2d, query, n0, ppst->pam2[0],pamscale);
-
- free(freq2d[0]);
- free(freq2d);
-
- free(query);
-}
-
-void
-scale_pssm(int **pssm2p, double **freq2d, unsigned char *query, int n0, int **pam2, double pamscale) {
- int i, qi, rj;
- double freq, new_lambda, lambda;
- int first, too_high;
- double scale, scale_high, scale_low;
-
- for (qi = 0 ; qi < n0 ; qi++) {
- for (rj = 0 ; rj < 20 ; rj++) {
- if (freq2d[qi][rj] > 1e-20) {
- freq = log(freq2d[qi][rj] /((double) (rrcounts[rj+1])/(double) rrtotal));
- freq /= pamscale; /* this gets us close to originial pam scores */
- freq2d[qi][rj] = freq;
- }
- else {
- /* when blastpgp decides to leave something out, it puts 0's in all the frequencies
- in the binary checkpoint file. In the ascii version, however, it uses BLOSUM62
- values. I will put in scoring matrix values as well */
-
- freq2d[qi][rj] = pam2[aascii[query[qi]]][rj+1];
- }
- }
- }
-
- /* now figure out the right scale */
- scale = 1.0;
- lambda = get_lambda(pam2, 20, 20, ustandard_aa);
-
- /* should be near 1.0 because of our initial scaling by ppst->pamscale */
- /* fprintf(stderr, "real_lambda: %g\n", lambda); */
-
- /* get initial high/low scale values: */
- first = 1;
- while (1) {
- fill_pam(pssm2p, n0, 20, freq2d, scale);
- new_lambda = get_lambda(pssm2p, n0, 20, query);
-
- if (new_lambda > lambda) {
- if (first) {
- first = 0;
- scale = scale_high = 1.0 + 0.05;
- scale_low = 1.0;
- too_high = 1;
- } else {
- if (!too_high) break;
- scale = (scale_high += scale_high - 1.0);
- }
- } else if (new_lambda > 0) {
- if (first) {
- first = 0;
- scale_high = 1.0;
- scale = scale_low = 1.0 - 0.05;
- too_high = 0;
- } else {
- if (too_high) break;
- scale = (scale_low += scale_low - 1.0);
- }
- } else {
- fprintf(stderr, "new_lambda (%g) <= 0; matrix has positive average score", new_lambda);
- exit(1);
- }
- }
-
- /* now do binary search between low and high */
- for (i = 0 ; i < 10 ; i++) {
- scale = 0.5 * (scale_high + scale_low);
- fill_pam(pssm2p, n0, 20, freq2d, scale);
- new_lambda = get_lambda(pssm2p, n0, 20, query);
-
- if (new_lambda > lambda) scale_low = scale;
- else scale_high = scale;
- }
-
- scale = 0.5 * (scale_high + scale_low);
- fill_pam(pssm2p, n0, 20, freq2d, scale);
-
- /*
- fprintf(stderr, "final scale: %g\n", scale);
-
- for (qi = 0 ; qi < n0 ; qi++) {
- fprintf(stderr, "%4d %c: ", qi+1, query[qi]);
- for (rj = 1 ; rj <= 20 ; rj++) {
- fprintf(stderr, "%4d", pssm2p[qi][rj]);
- }
- fprintf(stderr, "\n");
- }
- */
-}
-
-#if defined(SSEARCH) || (defined(PRSS) && !defined(FASTX))
-int
-parse_pssm_asn_fa(FILE *afd, int *n_rows, int *n_cols,
- unsigned char **query, double ***freqs,
- char *matrix, int *gap_open, int *gap_extend,
- double *lambda);
-
-/* the ASN.1 pssm includes information about the scoring matrix used
- (though not the gap penalty in the current version PSSM:2) The PSSM
- scoring matrix and gap penalties should become the default if they
- have not been set explicitly.
-*/
-
-int
-read_asn_pssm(unsigned char *aa0, int n0, int nsq,
- double pamscale, FILE *fp, struct pstruct *ppst) {
-
- int i, j, len, k;
- int qi, rj; /* qi - index query; rj - index residues (1-20) */
- int **pam2p;
- int first, too_high;
- unsigned char *query, ctmp;
- char dline[512];
- char matrix[MAX_SSTR];
- double psi2_lambda;
- double freq, **freq2d, lambda, new_lambda;
- double scale, scale_high, scale_low;
- int gap_open, gap_extend;
- int n_rows, n_cols;
-
- pam2p = ppst->pam2p[0];
-
- if (parse_pssm_asn_fa(fp, &n_rows, &n_cols, &query, &freq2d,
- matrix, &gap_open, &gap_extend, &psi2_lambda)<=0) {
- return -1;
- }
-
- if (!gap_set) {
- if (gap_open) {
- if (gap_open > 0) {gap_open = -gap_open;}
- ppst->gdelval = gap_open;
- }
- else if (strncmp(matrix,"BLOSUM62",8)==0) {
- ppst->gdelval = -11;
- }
- gap_set = 1;
- }
- if (!del_set) {
- if (gap_extend) {
- if (gap_extend > 0) {gap_extend = -gap_extend;}
- ppst->ggapval = gap_extend;
- }
- else if (strncmp(matrix,"BLOSUM62",8)==0) {
- ppst->ggapval = -1;
- }
- del_set = 1;
- }
-
- if (strncmp(matrix, "BLOSUM62", 8)== 0 && !ppst->pam_set) {
- strncpy(ppst->pamfile, "BL62", 120);
- standard_pam(ppst->pamfile,ppst,del_set, gap_set);
- if (!ppst->have_pam2) {
- alloc_pam (MAXSQ, MAXSQ, ppst);
- }
- init_pam2(ppst);
- ppst->pam_set = 1;
- }
-
- if (n_cols < n0) {
- fprintf(stderr, " query length: %d != n_cols: %d\n",n0, n_cols);
- exit(1);
- }
-
- scale_pssm(ppst->pam2p[0], freq2d, query, n0, ppst->pam2[0],pamscale);
-
- free(freq2d[0]);
- free(freq2d);
-
- free(query);
- return 1;
-}
-#endif
-
-void
-last_params(unsigned char *aa0, int n0,
- struct mngmsg *m_msg,
- struct pstruct *ppst
-#ifdef PCOMPLIB
- , struct qmng_str *qm_msg
-#endif
- ) {
- int i, nsq;
- FILE *fp;
-
- if (n0 < 0) { return;}
-
- ppst->n0 = m_msg->n0;
-
- if (ppst->ext_sq_set) { nsq = ppst->nsqx; }
- else {nsq = ppst->nsq;}
-
-/* currently, profiles are only available for SSEARCH, PRSS */
-#if defined(SSEARCH) || defined(PRSS)
-
- ppst->pam2p[0] = alloc_pam2p(n0,nsq);
- ppst->pam2p[1] = alloc_pam2p(n0,nsq);
-
- if (ppst->pam_pssm) {
- if ((ppst->pgpfile_type == 0) && (fp=fopen(ppst->pgpfile,"rb"))) {
- read_pssm(aa0, n0, ppst->nsq, ppst->pamscale, fp, 0, ppst);
- extend_pssm(aa0, n0, ppst);
- }
- else if ((ppst->pgpfile_type == 1) && (fp=fopen(ppst->pgpfile,"r"))) {
- read_pssm(aa0, n0, ppst->nsq, ppst->pamscale, fp, 1, ppst);
- extend_pssm(aa0, n0, ppst);
- }
-#if defined(SSEARCH) || (defined(PRSS) && !defined(FASTX))
- else if ((ppst->pgpfile_type == 2) && (fp=fopen(ppst->pgpfile,"rb"))) {
- if (read_asn_pssm(aa0, n0, ppst->nsq, ppst->pamscale, fp, ppst)>0) {
- extend_pssm(aa0, n0, ppst);
- }
- else {
- fprintf(stderr," Could not parse PSSM file: %s\n",ppst->pgpfile);
- ppst->pam_pssm = 0;
- return;
- }
- }
-#endif
- else {
- fprintf(stderr," Could not open PSSM file: %s\n",ppst->pgpfile);
- ppst->pam_pssm = 0;
- return;
- }
- }
-#endif
-
-#if defined(FASTF) || defined(FASTS) || defined(FASTM)
- m_msg->nm0 = 1;
- for (i=0; i<n0; i++)
- if (aa0[i]==EOSEQ || aa0[i]==ESS) m_msg->nm0++;
-
-/*
- for FASTS, we can do statistics in one of two different ways
- if there are <= 10 query fragments, then we calculate probabilistic
- scores for every library sequence. If there are > 10 fragments, this
- takes much too long and too much memory, so we use the old fashioned
- raw score only z-score normalized method initially, and then calculate
- the probabilistic scores for the best hits. To scale those scores, we
- also need a set of random probabilistic scores. So we do the qshuffle
- to get them.
-
- For FASTF, precalculating probabilities is prohibitively expensive,
- so we never do it; FASTF always acts like FASTS with nfrags>10.
-
-*/
-
-#if defined(FASTS) || defined(FASTM)
- if (m_msg->nm0 > 10) m_msg->escore_flg = 0;
- else m_msg->escore_flg = 1;
-#endif
-
- if (m_msg->escore_flg && (ppst->zsflag&1)) {
- m_msg->last_calc_flg = 0;
- m_msg->qshuffle = 0;
- }
- else { /* need random query, second set of 2000 scores */
- m_msg->last_calc_flg = 1;
- m_msg->qshuffle = 1;
- }
-#else
- m_msg->last_calc_flg = 0;
- m_msg->qshuffle = 0;
- m_msg->escore_flg = 0;
- m_msg->nm0 = 1;
-#endif
-
-/* adjust the ktup if appropriate */
-
- if (!ktup_set && pgm_def_arr[ppst->pgm_id].ktup > 0) {
- if (m_msg->qdnaseq == SEQT_PROT) {
- ppst->param_u.fa.ktup = pgm_def_arr[ppst->pgm_id].ktup;
-#if defined(FASTS) || defined(FASTM)
- if (n0 > 100) ppst->param_u.fa.ktup = 2;
-#endif
- if (n0 < 40) ppst->param_u.fa.ktup = 1;
- }
- else if (m_msg->qdnaseq == SEQT_DNA || m_msg->qdnaseq == SEQT_RNA) {
- if (n0 < 20) ppst->param_u.fa.ktup = 1;
-#if defined(FASTS) || defined(FASTM)
- /* with the current (April 12 2005) dropfs2.c - ktup cannot be > 2 */
- else ppst->param_u.fa.ktup = 2;
-#else
- else if (n0 < 50) ppst->param_u.fa.ktup = 2;
- else if (n0 < 100) ppst->param_u.fa.ktup = 3;
-#endif
- }
- }
-
-#ifdef PCOMPLIB
- qm_msg->nm0 = m_msg->nm0;
- qm_msg->escore_flg = m_msg->escore_flg;
- qm_msg->qshuffle = m_msg->qshuffle;
- qm_msg->pam_pssm = 0;
-#endif
-}
-
-/* given a good profile in ppst->pam2p[0], make an extended profile
- in ppst->pam2p[1]
-*/
-void
-extend_pssm(unsigned char *aa0, int n0, struct pstruct *ppst) {
-
- int i, j, nsq;
- int sa_x, sa_t, sa_b, sa_z;
- int **pam2p0, **pam2p1;
-
- nsq = ppst->nsq;
-
- pam2p0 = ppst->pam2p[0];
- pam2p1 = ppst->pam2p[1];
-
- sa_x = pascii['X'];
- sa_t = pascii['*'];
- sa_b = pascii['B'];
- sa_z = pascii['Z'];
-
- /* fill in boundaries, B, Z, *, X */
- for (i=0; i<n0; i++) {
- pam2p0[i][0] = -BIGNUM;
- pam2p0[i][sa_b] = (int)
- (((float)pam2p0[i][pascii['N']]+(float)pam2p0[i][pascii['D']]+0.5)/2.0);
- pam2p0[i][sa_z] = (int)
- (((float)pam2p0[i][pascii['Q']]+(float)pam2p0[i][pascii['E']]+0.5)/2.0);
- pam2p0[i][sa_x] = ppst->pam_xm;
- pam2p0[i][sa_t] = ppst->pam_xm;
- }
-
- /* copy pam2p0 into pam2p1 */
- for (i=0; i<n0; i++) {
- pam2p1[i][0] = -BIGNUM;
- for (j=1; j<=ppst->nsq; j++) {
- pam2p1[i][j] = pam2p0[i][j];
- }
- }
-
- /* then fill in extended characters, if necessary */
- if (ppst->ext_sq_set) {
- for (i=0; i<n0; i++) {
- for (j=1; j<=ppst->nsq; j++) {
- pam2p0[i][nsq+j] = pam2p0[i][j];
- pam2p1[i][nsq+j] = ppst->pam_xm;
- }
- }
- }
-}
git://source.jalview.org / jabaws.git / blobdiff