« back to all changes in this revision
Viewing changes to src/hhalign/hhfunc-C.h
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Package Import Robot
- Author(s): Olivier Sallou
- Date: 2011-12-14 11:21:56 UTC
- Revision ID: package-import@ubuntu.com-20111214112156-y3chwm0t4yn2nevm
Tags: upstream-1.0.3
Import upstream version 1.0.3
- files added:
- m4
- src
- src/clustal
- src/hhalign
- src/kmpp
- src/squid
added
removed
src/hhalign/hhfunc-C.h
1
/* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
2
3
/*********************************************************************
4
* Clustal Omega - Multiple sequence alignment
5
*
6
* Copyright (C) 2010 University College Dublin
7
*
8
* Clustal-Omega is free software; you can redistribute it and/or
9
* modify it under the terms of the GNU General Public License as
10
* published by the Free Software Foundation; either version 2 of the
11
* License, or (at your option) any later version.
12
*
13
* This file is part of Clustal-Omega.
14
*
15
********************************************************************/
16
17
/*
18
* RCS $Id: hhfunc-C.h 256 2011-06-23 13:57:28Z fabian $
19
*/
20
21
/*
22
* Changelog: Michael Remmert made changes to hhalign stand-alone code
23
* FS implemented some of the changes on 2010-11-11 -> MR1
24
*
25
* did not incorporate SS prediction PSIpred (yet); functions are:
26
* CalculateSS(3),
27
*/
28
29
30
// hhfunc.C
31
32
33
/**
34
* AddBackgroundFrequencies()
35
*
36
* @brief add background frequencies (derived from HMM) to
37
* sequence/profile
38
*
39
* @param[in,out] ppfFreq,
40
* [in] residue frequencies of sequence/profile,
41
* [out] overlayed with HMM background frequencies
42
* @param[in,out] ppfPseudoF,
43
* [in] residue frequencies+pseudocounts of sequence/profile,
44
* [out] overlayed with HMM background frequencies+pseudocounts
45
* @param[in] iSeqLen,
46
* length of sequence/profile (not aligned to HMM)
47
* @pram[in] prHMM,
48
* background HMM
49
* @param[in] ppcSeq,
50
* sequences/profile to be 'softened'
51
* @param[in] pcPrealigned,
52
* sequence aligned to HMM, this is not quite a consensus,
53
* it is identical to 1st sequence but over-writes gap information,
54
* if other sequences in profile have (non-gap) residues
55
* @param[in] iPreCnt
56
* number of sequences pre-aligned (pcPrealigned is 'consensus' of these sequences)
57
* @param[in] pcRepresent
58
* sequence representative of HMM, aligned to pcSeq0
59
*/
60
void
61
AddBackgroundFrequencies(float **ppfFreq, float **ppfPseudoF, float **ppfPseudoTr,
62
int iSeqLen, hmm_light *prHMM,
63
char **ppcSeq, char *pcPrealigned, int iPreCnt, char *pcRepresent)
64
{
65
66
char *pcS = pcPrealigned; /* current residue in pre-aligned sequence */
67
char *pcH = pcRepresent; /* current residue in pre-aligned HMM */
68
int iS = 0; /* position in un-aligned sequence (corresponds to pcS) */
69
int iH = 0; /* position in un-aligned HMM (corresponds to pcH) */
70
int iA; /* residue iterator */
71
int iT; /* transition state iterator */
72
float fFWeight = 0.50 / sqrt((float)(iPreCnt)); /* weight of backgroud frequencies */ /* FIXME: tune value, 0.50 default */
73
//float fFWeight = 0.75;
74
float fFThgiew = UNITY - fFWeight; /* weight of 'true' frequencies */
75
float fGWeight = 0.50 / sqrt((float)(iPreCnt)); /* weight of backgroud transitions */ /* FIXME: tune value, 0.50 default */
76
//float fGWeight = 0.50 /*/ (float)(iPreCnt)*/; /* weight of backgroud transitions */ /* FIXME: tune value, 0.50 default */
77
float fGThgiew = UNITY - fGWeight; /* weight of 'true' transitions */
78
float fAux;
79
/* zf1SeqTrans[] are the (phenomenological) transition probabilities (+ pseudo-counts)
80
for a single sequence. It is almost certain (0.99) to stay in a match state, and very
81
unlikely (0.05) to go from match to insertion/deletion */
82
float zf1SeqTrans[] = {0.98, 0.01, 0.01, 0.25, 0.75, 0.25, 0.75};
83
float zf1SeqInit[] = {0.98, 0.01, 0.01, 0.99, 0.01, 0.99, 0.01};
84
float zf1SeqDel[] = {0.24, 0.01, 0.75, 0.01, 0.01, 0.01, 0.01};
85
float zf1SeqRevrt[] = {0.98, 0.01, 0.01, 0.01, 0.01, 0.99, 0.01};
86
87
if ( (NULL == pcPrealigned) || (NULL == pcRepresent) ){
88
/*printf("%s/%s:%d: WARNING HMM=NULL -- didn't think I would get here (carry on, no danger)\n",
89
__FUNCTION__, __FILE__, __LINE__);*/
90
return;
91
}
92
93
if (NULL == prHMM->p){
94
printf("%s:%s:%d: WARNING -- Background Pseudocounts point to NULL\n"
95
"\tthis is not intended - don't add background but CONTINUE\n",
96
__FUNCTION__, __FILE__, __LINE__);
97
return;
98
99
}
100
101
/* FIXME: should be 0 (FS thinks) but -1 gives better results */
102
iH = iS = 0/*-1*//*+1*/;
103
while ( ('\0' != *pcS) && ('\0' != *pcH) ){
104
105
if ( ('-' != *pcH) && ('-' != *pcS) ){
106
iH++;
107
iS++;
108
109
#if 0
110
/* match state
111
* - HMM had a gap in previous position (now closed)
112
* FIXME: this does not really work */
113
if ((iH > 0) && ('-' == *(pcH-1))){
114
115
for (iT = 0; iT < STATE_TRANSITIONS; iT++){
116
ppfPseudoTr[iS-1][iT] = log2f(zf1SeqRevrt[iT]);
117
}
118
}
119
#endif
120
121
#if 1
122
/* do frequencies -- this is not really useful;
123
frequencies are derived from HMM
124
and contain already pseudocounts (PCs),
125
adding frequencies and then PCs will add PCs _twice_
126
results are better than not to add them,
127
but not as good as PCs */
128
for (iA = 0; iA < AMINOACIDS; iA++){
129
fAux = fFThgiew * ppfFreq[iS][iA] + fFWeight * prHMM->f[iH][iA];
130
ppfFreq[iS][iA] = fAux;
131
}
132
#endif
133
/* do pseudo-counts */
134
for (iA = 0; iA < AMINOACIDS; iA++){
135
fAux =
136
fFThgiew * ppfPseudoF[iS][iA] + fFWeight * prHMM->p[iH][iA];
137
ppfPseudoF[iS][iA] = fAux;
138
}
139
#if 0 /* do state transitions */
140
for (iT = 0; iT < STATE_TRANSITIONS; iT++){
141
#if 1
142
/* this is a very crude method,
143
which averages the logarithms of the transitions,
144
effectively performing a geometric mean -
145
this presumably violates normalisation.
146
however, results are surprisingly good */
147
fAux =
148
fGThgiew * ppfPseudoTr[iS][iT] +
149
fGWeight * prHMM->tr[iH][iT];
150
ppfPseudoTr[iS][iT] = fAux;
151
#else /* crude averaging */
152
/* There are 2 things to consider:
153
(1) one should really blend the probabilities of
154
the transitions, however, by default we have
155
the logarithms thereof.
156
So must exponentiate, blend, and then take log again.
157
This is expensive, and does not seem to lead to better
158
results than blending the logarithms
159
(and violating normalisation)
160
(2) transition probabilities for a single sequence
161
are really easy, there are no insert/delete transitions.
162
However, there is a begin state that is different
163
from the main body.
164
But again, this does not seem to make a blind bit
165
of difference
166
*/
167
if (iS > 0){
168
fAux =
169
fGThgiew * zf1SeqTrans[iT] +
170
fGWeight * prHMM->linTr[iH][iT];
171
ppfPseudoTr[iS][iT] = log2f(fAux);
172
}
173
else {
174
fAux =
175
fGThgiew * zf1SeqInit[iT] +
176
fGWeight * prHMM->linTr[iH][iT];
177
ppfPseudoTr[iS][iT] = log2f(fAux);
178
}
179
#endif /* mixing of linear/log */
180
} /* 0 <= iT < STATE_TRANSITIONS */
181
#endif /* did state transitions */
182
} /* was Match -- neither HMM nor sequence have gap */
183
184
else if ('-' == *pcH){
185
/* sequence opened up gap in HMM */
186
#if 0
187
if ((iH > 0) && ('-' != *(pcH-1)) && (iS > 0)){
188
/* this is the first gap in HMM
189
* FIXME: this does not really work */
190
for (iT = 0; iT < STATE_TRANSITIONS; iT++){
191
ppfPseudoTr[iS-1][iT] = log2f(zf1SeqDel[iT]);
192
}
193
}
194
else {
195
/* do nothing, keep single sequence values exactly as they are*/
196
}
197
#endif
198
iS++;
199
}
200
else if ('-' == *pcS){
201
/* here the single sequence has a gap,
202
and the HMM (as a whole) does not. There may be individual gaps
203
in the HMM at this stage. By ignoring this we say that the HMM
204
dominates the overall behaviour - as in the M2M state as well
205
*/
206
iH++;
207
}
208
209
pcH++;
210
pcS++;
211
212
} /* !EO[seq/hmm] */
213
214
return;
215
216
} /* this is the end of AddBackgroundFrequencies() */
217
218
219
220
/**
221
* ReadAndPrepare()
222
*
223
* @brief Read HMM input file or transfer alignment
224
* and add pseudocounts etc.
225
*
226
* @param[in] iRnPtype
227
* type of read/prepare
228
* enum {INTERN_ALN_2_HMM = 0, READ_ALN_2_HMM, READ_HMM_2_HMM};
229
* @param[in] ppcProf
230
* alignment
231
* @param[in] iCnt
232
* number of seqs in alignment
233
* @param[in,out] prHMM,
234
* [in] if sequences read/prepared, [out] if HMM from file
235
* @param[in] pcPrealigned,
236
* (single) sequence aligned to background HMM
237
* @param[in] pcRepresent,
238
* sequence representing HMM aligned to individual sequence
239
* param[in] pdExWeights
240
* (external) sequence weights, derived from tree
241
* @param[in] infile
242
* name of file with HMM info (formerly also alignment)
243
* @param[out] q
244
* HMM structure with transition probabilities, residue frequencies etc
245
* @param[???] qali
246
* FIXME: what is qali?
247
*
248
* @return FAILURE on error, OK otherwise
249
*/
250
int
251
ReadAndPrepare(int iRnPtype,
252
char **ppcProf, int iCnt, hmm_light *prHMM,
253
char *pcPrealigned, char *pcRepresent, double *pdExWeights,
254
char* infile, HMM& q, Alignment* qali=NULL) {
255
256
//#ifndef CLUSTALO_NOFILE
257
char path[NAMELEN];
258
259
/* NOTE: there are different scenarios:
260
261
(i) ("" != infile) - read HMM from file,
262
transfer frequency/transition/pseudo-count (f/t/p) info into prHMM
263
264
(ii) ('\0' != ppcProf[0]) - transfer sequence/alignment into q/qali,
265
don't save f/t/p into prHMM,
266
on the contrary, if prior f/t/p available then add it to q/qali,
267
this is only done if (1==iCnt)
268
269
(iii) ('\0' == ppcProf[0]) - re-cycle old HMM information
270
recreate a HMM from previous data
271
*/
272
273
/********************************/
274
/*** (o) Recycle internal HMM ***/
275
/********************************/
276
if ( (INTERN_ALN_2_HMM == iRnPtype) && (iCnt <= 0) ){
277
278
/* NOTE: here we are writing into a HMM structure/class;
279
memory has been allocated for this in hhalign.cpp;
280
however, as iCnt<=0, there may not be memory for
281
prHMM->n_display sequences/names.
282
But then, there doesn't have to be.
283
At this stage we are just copying one HMM into another HMM,
284
sequences are irrelevant. The only sequence of (marginal)
285
interest is the consensus sequence */
286
/* FIXME: check that prHMM is valid -- how? */
287
288
const int ciCons = 0;
289
const int ciNoof = ciCons+1;
290
const int ciInvalid = -1;
291
q.n_display = ciNoof; /* only consensus */
292
q.sname = NULL;
293
q.ncons = ciCons;
294
q.nfirst = ciCons;//prHMM->nfirst;
295
q.nss_dssp = ciInvalid;//prHMM->nss_dssp;
296
q.nsa_dssp = ciInvalid;//prHMM->nsa_dssp;
297
q.nss_pred = ciInvalid;//prHMM->nss_pred;
298
q.nss_conf = ciInvalid;//prHMM->nss_conf;
299
q.L = prHMM->L;
300
q.N_in = prHMM->N_in;
301
q.N_filtered = prHMM->N_filtered;
302
/* NOTE: I (FS) think that only ever 1 sequence will be transferred here,
303
that is, the consensus sequence. However, we might want to allow
304
(in the future) to transfer more sequences,
305
hence the awkward for() loop */
306
#if 0
307
for (int i = prHMM->ncons+0; i < prHMM->ncons+q.n_display; i++){
308
/* NOTE: In the original hhalign code the first position
309
is kept open ('\0'). This makes it difficult to use
310
string functions like strlen/strdup etc.
311
Insert a temporary gap (.) to facilitate string operations */
312
char cHead = prHMM->seq[i][0];
313
prHMM->seq[i][0] = '.';
314
q.seq[i] = strdup(prHMM->seq[i]);
315
prHMM->seq[i][0] = q.seq[i][0] = cHead;
316
}
317
#else
318
{
319
char cHead = prHMM->seq[prHMM->ncons][0];
320
prHMM->seq[prHMM->ncons][0] = '.';
321
q.seq[q.ncons] = strdup(prHMM->seq[prHMM->ncons]);
322
prHMM->seq[prHMM->ncons][0] = q.seq[q.ncons][0] = cHead;
323
}
324
#endif
325
const int NEFF_SCORE = 10; /* FIXME: Magic Number */
326
for (int i = 0; i < prHMM->L+1; i++){
327
q.Neff_M[i] = prHMM->Neff_M[i];
328
q.Neff_I[i] = prHMM->Neff_I[i];
329
q.Neff_D[i] = prHMM->Neff_D[i];
330
}
331
q.Neff_HMM = prHMM->Neff_HMM;
332
/* skip longname,name,file,fam,sfam,fold,cl */
333
q.lamda = prHMM->lamda;
334
q.mu = prHMM->mu;
335
336
HMMshadow rShadow = {0}; /* friend of HMM to access private members */
337
rShadow.copyShadowToHMM(q, *prHMM);
338
339
/* skip trans_lin,ss_dssp,sa_dssp,ss_pred,ss_conf,Xcons */
340
/* pav already done in copyShadowToHMM */
341
/* skip pnul */
342
343
return OK;
344
345
346
} /* INTERN_ALN_2_HMM && iCnt<=0 */
347
348
/******************************/
349
/*** (i) Read HMM from file ***/
350
/******************************/
351
char line[LINELEN] = {0}; // input line
352
FILE *inf = NULL;
353
//if ( (0 != strcmp(infile,"")) /*&& (iCnt > 0)*/ )
354
if ( (READ_HMM_2_HMM == iRnPtype) || (READ_ALN_2_HMM == iRnPtype)) {
355
356
if (0 == strcmp(infile,"")){
357
printf("%s:%s:%d:\n"
358
"\texpected to re %s from file but no file specified\n"
359
""
360
, __FUNCTION__, __FILE__, __LINE__
361
, (READ_HMM_2_HMM==iRnPtype?"HMM":"alignment"));
362
return FAILURE;
363
}
364
inf = fopen(infile, "r");
365
if (!inf) OpenFileError(infile);
366
Pathname(path,infile);
367
368
//}
369
//else {
370
//inf = stdin;
371
//if (v>=2) printf("Reading HMM / multiple alignment from standard input ...\n(To get a help list instead, quit and type %s -h.)\n",program_name);
372
//*path='\0';
373
//}
374
375
fgetline(line,LINELEN-1,inf);
376
}
377
378
//if ( (0 != strcmp(infile,"")) && (iCnt > 0) )
379
if ( (READ_HMM_2_HMM == iRnPtype) ){
380
381
if (0 == strcmp(infile,"")){
382
printf("%s:%s:%d: expected to read HMM from file but no file-name\n",
383
__FUNCTION__, __FILE__, __LINE__);
384
}
385
386
// Is infile a HMMER3 file? /* MR1 */
387
if (!strncmp(line,"HMMER3",6))
388
{
389
if (v>=2) {
390
cout<<"Query file is in HMMER3 format\n";
391
cout<<"WARNING! Use of HMMER3 format as input results in dramatically loss of sensitivity!\n";
392
}
393
394
// Read 'query HMMER file
395
rewind(inf);
396
q.ReadHMMer3(inf,path);
397
398
// Don't add transition pseudocounts to query!!
399
400
// Generate an amino acid frequency matrix from f[i][a] with full pseudocount admixture (tau=1) -> g[i][a]
401
q.PreparePseudocounts();
402
403
// DON'T ADD amino acid pseudocounts to query: pcm=0! q.p[i][a] = f[i][a]
404
q.AddAminoAcidPseudocounts(0, par.pca, par.pcb, par.pcc);
405
406
/* further down there is a q.Log2LinTransitionProbs
407
but only if (iCnt>0), however, we still need it it here (i think),
408
there is no danger of doing this twice, as trans_lin is checked
409
FIXME (FS, 2011-01-12) */
410
/* further down there is a q.Log2LinTransitionProbs
411
but only if (iCnt>0), however, we still need it it here (i think),
412
there is no danger of doing this twice, as trans_lin is checked
413
FIXME (FS, 2011-01-12) */
414
//if (par.forward >= 1)
415
{
416
q.Log2LinTransitionProbs(1.0);
417
}
418
419
}
420
// ... or Is infile an old HMMER file?
421
else if (!strncmp(line,"HMMER",5)) {
422
if (v>=2) {
423
cout<<"Query file is in HMMER format\n";
424
cout<<"WARNING! Use of HMMER format as input results in dramatically loss of sensitivity!\n";
425
}
426
427
// Read 'query HMMER file
428
q.ReadHMMer(inf,path);
429
if (v>=2 && q.Neff_HMM>11.0)
430
fprintf(stderr,"WARNING: HMM %s looks too diverse (Neff=%.1f>11). Better check the underlying alignment... \n",q.name,q.Neff_HMM);
431
432
// Don't add transition pseudocounts to query!!
433
434
// Generate an amino acid frequency matrix from f[i][a] with full pseudocount admixture (tau=1) -> g[i][a]
435
q.PreparePseudocounts();
436
437
// DON'T ADD amino acid pseudocounts to query: pcm=0! q.p[i][a] = f[i][a]
438
q.AddAminoAcidPseudocounts(0, par.pca, par.pcb, par.pcc);
439
440
/* further down there is a q.Log2LinTransitionProbs
441
but only if (iCnt>0), however, we still need it it here (i think),
442
there is no danger of doing this twice, as trans_lin is checked
443
FIXME (FS, 2011-01-12) */
444
//if (par.forward >= 1)
445
{
446
q.Log2LinTransitionProbs(1.0);
447
}
448
449
} /* it was a HMMer file */
450
451
// ... or is it an hhm file?
452
else if (!strncmp(line,"NAME",4) || !strncmp(line,"HH",2)) {
453
454
if (v>=2) cout<<"Query file is in HHM format\n";
455
456
// Rewind to beginning of line and read query hhm file
457
rewind(inf);
458
q.Read(inf,path);
459
if (v>=2 && q.Neff_HMM>11.0)
460
fprintf(stderr,"WARNING: HMM %s looks too diverse (Neff=%.1f>11). Better check the underlying alignment... \n",q.name,q.Neff_HMM);
461
462
// Add transition pseudocounts to query -> q.p[i][a]
463
q.AddTransitionPseudocounts();
464
465
// Generate an amino acid frequency matrix from f[i][a] with full pseudocount admixture (tau=1) -> g[i][a]
466
q.PreparePseudocounts();
467
468
// Add amino acid pseudocounts to query: q.p[i][a] = (1-tau)*f[i][a] + tau*g[i][a]
469
q.AddAminoAcidPseudocounts();
470
471
} /* it was a HHM file */
472
else {
473
fprintf(stderr, "%s:%s:%d: Unknown HMM format in infile\n"
474
"infile=%s, #seq=%d\n"
475
, __FUNCTION__, __FILE__, __LINE__
476
, infile, iCnt);
477
return FAILURE;
478
}
479
480
/*fclose(inf);*/
481
482
/*** transfer class info to struct */
483
prHMM->n_display = q.n_display;
484
/* ignore sname*/
485
prHMM->seq = (char **)calloc((q.n_display+1), sizeof(char *));
486
/* FIXME valgrind says bytes get lost in the above calloc during
487
* hmm-iteration
488
*/
489
for (int i = 0; i < q.n_display; i++){
490
/* NOTE: In the original hhalign code the first position
491
is kept open ('\0'). This makes it difficult to use
492
string functions like strlen/strdup etc.
493
Insert a temporary gap (.) to facilitate string operations */
494
char cHead = q.seq[i][0];
495
q.seq[i][0] = '.';
496
prHMM->seq[i] = strdup(q.seq[i]);
497
q.seq[i][0] = prHMM->seq[i][0] = cHead;
498
}
499
prHMM->ncons = q.ncons;
500
prHMM->nfirst = q.nfirst;
501
prHMM->nss_dssp = q.nss_dssp;
502
prHMM->nsa_dssp = q.nsa_dssp;
503
prHMM->nss_pred = q.nss_pred;
504
prHMM->nss_conf = q.nss_conf;
505
prHMM->L = q.L;
506
prHMM->N_in = q.N_in;
507
prHMM->N_filtered = q.N_filtered;
508
const int NEFF_SCORE = 10; /* FIXME: Magic Number */
509
prHMM->Neff_M = (float *)calloc(prHMM->L+1, sizeof(float));
510
prHMM->Neff_I = (float *)calloc(prHMM->L+1, sizeof(float));
511
prHMM->Neff_D = (float *)calloc(prHMM->L+1, sizeof(float));
512
/*for (int i = 0; i < prHMM->L+1; i++){
513
prHMM->Neff_M[i] = prHMM->Neff_I[i] = prHMM->Neff_D[i] = NEFF_SCORE;
514
}*/
515
prHMM->Neff_HMM = q.Neff_HMM;
516
/* skip longname,name,file,fam,sfam,fold,cl */
517
prHMM->lamda = q.lamda;
518
prHMM->mu = q.mu;
519
520
HMMshadow rShadow = {0}; /* friend of HMM to access private members */
521
rShadow.copyHMMtoShadow(q);
522
523
prHMM->f = (float **)calloc(prHMM->L+1, sizeof(float *));
524
prHMM->g = (float **)calloc(prHMM->L+1, sizeof(float *));
525
prHMM->p = (float **)calloc(prHMM->L+1, sizeof(float *));
526
prHMM->tr = (float **)calloc(prHMM->L+1, sizeof(float *));
527
prHMM->linTr = (float **)calloc(prHMM->L+1, sizeof(float *));
528
for (int i = 0; i < prHMM->L+1; i++){
529
prHMM->f[i] = (float *)calloc(AMINOACIDS, sizeof(float));
530
prHMM->g[i] = (float *)calloc(AMINOACIDS, sizeof(float));
531
prHMM->p[i] = (float *)calloc(AMINOACIDS, sizeof(float));
532
for (int j = 0; j < AMINOACIDS; j++){
533
prHMM->f[i][j] = (float)rShadow.f[i][j];
534
prHMM->g[i][j] = (float)rShadow.g[i][j];
535
prHMM->p[i][j] = (float)rShadow.p[i][j];
536
}
537
prHMM->tr[i] = (float *)calloc(STATE_TRANSITIONS, sizeof(float));
538
prHMM->linTr[i] = (float *)calloc(STATE_TRANSITIONS, sizeof(float));
539
for (int j = 0; j< STATE_TRANSITIONS; j++){
540
prHMM->tr[i][j] = (float)rShadow.tr[i][j];
541
prHMM->linTr[i][j] = fpow2(rShadow.tr[i][j]);
542
}
543
} /*0 <= i < prHMM->L+1 */
544
/* skip trans_lin,ss_dssp,sa_dssp,ss_pred,ss_conf,Xcons */
545
for (int j = 0; j < AMINOACIDS; j++){
546
prHMM->pav[j] = (float)rShadow.pav[j];
547
}
548
/* skip pnul */
549
550
} /* have read HHM from file */
551
/*else if ( ((NULL != ppcProf) && (iCnt > 0) && ('\0' != ppcProf[0][0])) ||
552
( (0 != strcmp(infile,"") && (0 == iCnt) )) )*/
553
else if ( (INTERN_ALN_2_HMM == iRnPtype) || (READ_ALN_2_HMM == iRnPtype) ) {
554
555
if ( (INTERN_ALN_2_HMM == iRnPtype) &&
556
( (NULL == ppcProf) || (iCnt <= 0) || ('\0' == ppcProf[0][0]) ) ){
557
printf("%s:%s:%d: was expecting internal alignment but\n"
558
"\tppcProf=%p, #seq=%d, ppcProf[0][0]=%c\n"
559
, __FUNCTION__, __FILE__, __LINE__
560
, ppcProf, iCnt, ppcProf[0][0]);
561
return FAILURE;
562
}
563
else if ( (READ_ALN_2_HMM == iRnPtype) &&
564
(0 == strcmp(infile,"")) ){
565
printf("%s:%s:%d: was expecting to read alignment from file but no filename\n"
566
, __FUNCTION__, __FILE__, __LINE__);
567
return FAILURE;
568
}
569
570
/*******************************/
571
/*** (ii) it is an alignment ***/
572
/*******************************/
573
/* transfer alignment information from clustal character array
574
into pali/q/t classes */
575
/*
576
* NOTE that emissions in HMMer file format contain pseudo-counts.
577
* HHM file format does not contain emission pseudo-counts.
578
* the structure that stores background HMM information does contain pseudo-counts
579
*/
580
Alignment* pali;
581
if (qali==NULL){
582
pali=new Alignment(iCnt); /* FIXME: pass in iCnt to get rid of MAXSEQ */
583
}
584
else{
585
pali=qali;
586
}
587
588
if (par.calibrate) {
589
printf("\nError in %s: only HHM files can be calibrated.\n",program_name);
590
printf("Build an HHM file from your alignment with 'hhmake -i %s' and rerun hhsearch with the hhm file\n\n",infile);
591
exit(1);
592
}
593
594
if (v>=2 && strcmp(infile,"stdin")) cout<<infile<<" is in A2M, A3M or FASTA format\n";
595
596
/* Read alignment from infile into matrix X[k][l] as ASCII
597
(and supply first line as extra argument) */
598
//if (iCnt > 0)
599
if (INTERN_ALN_2_HMM == iRnPtype){
600
pali->Transfer(ppcProf, iCnt);
601
}
602
//else if (0 == iCnt)
603
else if (READ_ALN_2_HMM == iRnPtype){
604
pali->Read(inf, infile, line);
605
}
606
else {
607
printf("%s:%s:%d: FATAL problem\n"
608
"infile = (%s), #seq = %d\n"
609
, __FUNCTION__, __FILE__, __LINE__
610
, infile, iCnt);
611
return FAILURE;
612
}
613
614
/* Convert ASCII to int (0-20), throw out all insert states,
615
record their number in I[k][i]
616
and store marked sequences in name[k] and seq[k] */
617
pali->Compress(infile);
618
619
/* Sort out the nseqdis most dissimilar sequences for display
620
in the output alignments */
621
pali->FilterForDisplay(par.max_seqid, par.coverage, par.qid,
622
par.qsc,par.nseqdis);
623
624
// Filter alignment for min score per column with core query profile, defined by coverage_core and qsc_core
625
//if (par.coresc>-10) pali->HomologyFilter(par.coverage_core, par.qsc_core, par.coresc);
626
627
/* Remove sequences with seq. identity larger than seqid percent
628
(remove the shorter of two) */
629
pali->N_filtered = pali->Filter(par.max_seqid, par.coverage,
630
par.qid, par.qsc, par.Ndiff);
631
632
/* Calculate pos-specific weights,
633
AA frequencies and transitions -> f[i][a], tr[i][a] */
634
pali->FrequenciesAndTransitions(q);
635
if (v>=2 && q.Neff_HMM>11.0)
636
fprintf(stderr,"WARNING: alignment %s looks too diverse (Neff=%.1f>11). Better check it with an alignment viewer... \n",q.name,q.Neff_HMM);
637
638
/*printf("%d %d %f %d (N,Nf,Neff,L) %s:%s:%d\n"
639
, q.N_in, q.N_filtered, q.Neff_HMM, q.L, __FUNCTION__, __FILE__, __LINE__);*/
640
641
// Add transition pseudocounts to query -> p[i][a]
642
q.AddTransitionPseudocounts();
643
644
/* Generate an amino acid frequency matrix from f[i][a]
645
with full pseudocount admixture (tau=1) -> g[i][a] */
646
q.PreparePseudocounts();
647
648
/* Add amino acid pseudocounts to query:
649
p[i][a] = (1-tau)*f[i][a] + tau*g[i][a] */
650
q.AddAminoAcidPseudocounts();
651
652
/* ****** add aligned background pseudocounts ***** */
653
HMMshadow rShadowQ = {0};
654
rShadowQ.copyHMMtoShadow(q);
655
656
AddBackgroundFrequencies(rShadowQ.f, rShadowQ.p, rShadowQ.tr,
657
q.L, prHMM,
658
q.seq, pcPrealigned, iCnt, pcRepresent);
659
660
661
if (qali==NULL){
662
delete(pali); pali = NULL;
663
}
664
665
} /* else if (NULL != ppcProf) // not hmmer/hhalign but alignment */
666
667
//else if ((prHMM->L > 0) && ('\0' == ppcProf[0][0]))
668
else if (INTERN_HMM_2_HMM == iRnPtype){
669
670
/******************************************/
671
/*** (iii) re-cycle old HMM information ***/
672
/******************************************/
673
674
if (prHMM->L <= 0){
675
printf("%s:%s:%d: was expecting to copy HMM structure but L=%d\n"
676
, __FUNCTION__, __FILE__, __LINE__, prHMM->L);
677
}
678
679
printf("%s:%s:%d: RE-CYCLING HMM\n", __FUNCTION__, __FILE__, __LINE__);
680
681
#if 0
682
q.n_display = prHMM->n_display;
683
/* ignore sname*/
684
for (int i = 0; i < q.n_display; i++){
685
/* NOTE: In the original hhalign code the first position
686
is kept open ('\0'). This makes it difficult to use
687
string functions like strlen/strdup etc.
688
Insert a temporary gap (.) to facilitate string operations */
689
char cHead = prHMM->seq[i][0];
690
prHMM->seq[i][0] = '.';
691
q.seq[i] = strdup(prHMM->seq[i]);
692
q.seq[i][0] = prHMM->seq[i][0] = cHead;
693
}
694
q.nfirst = prHMM->nfirst;
695
#else
696
q.n_display = 1;
697
q.nfirst = 0;
698
char cHead = prHMM->seq[prHMM->nfirst][0];
699
prHMM->seq[prHMM->nfirst][0] = '.';
700
q.seq[0] = strdup(prHMM->seq[prHMM->nfirst]);
701
q.seq[q.n_display-1][0] = prHMM->seq[prHMM->nfirst][0] = cHead;
702
#endif
703
q.ncons = prHMM->ncons;
704
q.nss_dssp = prHMM->nss_dssp;
705
q.nsa_dssp = prHMM->nsa_dssp;
706
q.nss_pred = prHMM->nss_pred;
707
q.nss_conf = prHMM->nss_conf;
708
q.L = prHMM->L;
709
q.N_in = prHMM->N_in;
710
q.N_filtered = prHMM->N_filtered;
711
#define NEFF_SCORE 10 /* FIXME: Magic Number */
712
/*for (int i; i < prHMM->L+1; i++){
713
q.Neff_M[i] = q.Neff_I[i] = q.Neff_D[i] = NEFF_SCORE;
714
}*/
715
q.Neff_HMM = prHMM->Neff_HMM;
716
/* skip longname,name,file,fam,sfam,fold,cl */
717
q.lamda = prHMM->lamda;
718
q.mu = prHMM->mu;
719
720
HMMshadow rShadow = {0}; /* friend of HMM to access private members */
721
rShadow.copyShadowToHMM(q, *prHMM);
722
723
}
724
725
if (iCnt > 0){
726
if (par.forward>=1) q.Log2LinTransitionProbs(1.0);
727
}
728
729
if (NULL != inf){
730
fclose(inf);
731
}
732
733
return OK;
734
735
} /*** end: ReadAndPrepare() ***/
736
737
/**
738
* FreeHMMstruct()
739
*
740
* @brief FIXME
741
*
742
* @param[in,out]
743
*/
744
extern "C" void
745
FreeHMMstruct(hmm_light *prHMM)
746
{
747
int i;
748
749
if (NULL != prHMM->f){
750
for (i = 0; i < prHMM->L+1; i++){
751
if (NULL != prHMM->f[i]){
752
free(prHMM->f[i]); prHMM->f[i] = NULL;
753
}
754
} /* 0 <= i < prHMM->L+1 */
755
free(prHMM->f); prHMM->f = NULL;
756
}
757
if (NULL != prHMM->g){
758
for (i = 0; i < prHMM->L+1; i++){
759
if (NULL != prHMM->g[i]){
760
free(prHMM->g[i]); prHMM->g[i] = NULL;
761
}
762
} /* 0 <= i < prHMM->L+1 */
763
free(prHMM->g); prHMM->g = NULL;
764
}
765
if (NULL != prHMM->p){
766
for (i = 0; i < prHMM->L+1; i++){
767
if (NULL != prHMM->p[i]){
768
free(prHMM->p[i]); prHMM->p[i] = NULL;
769
}
770
} /* 0 <= i < prHMM->L+1 */
771
free(prHMM->p); prHMM->p = NULL;
772
}
773
if (NULL != prHMM->tr){
774
for (i = 0; i < prHMM->L+1; i++){
775
if (NULL != prHMM->tr[i]){
776
free(prHMM->tr[i]); prHMM->tr[i] = NULL;
777
}
778
} /* 0 <= i < prHMM->L+1 */
779
free(prHMM->tr); prHMM->tr = NULL;
780
}
781
if (NULL != prHMM->linTr){
782
for (i = 0; i < prHMM->L+1; i++){
783
if (NULL != prHMM->linTr[i]){
784
free(prHMM->linTr[i]); prHMM->linTr[i] = NULL;
785
}
786
} /* 0 <= i < prHMM->L+1 */
787
free(prHMM->linTr); prHMM->linTr = NULL;
788
}
789
790
if (NULL != prHMM->Neff_M){
791
free(prHMM->Neff_M); prHMM->Neff_M = NULL;
792
}
793
if (NULL != prHMM->Neff_I){
794
free(prHMM->Neff_I); prHMM->Neff_I = NULL;
795
}
796
if (NULL != prHMM->Neff_D){
797
free(prHMM->Neff_D); prHMM->Neff_D = NULL;
798
}
799
800
if (NULL != prHMM->seq){
801
for (i = 0; i < prHMM->n_display; i++){
802
if (NULL != prHMM->seq[i]){
803
free(prHMM->seq[i]); prHMM->seq[i] = NULL;
804
}
805
}
806
free(prHMM->seq); prHMM->seq = NULL;
807
}
808
809
memset(prHMM, 0, sizeof(hmm_light));
810
811
} /*** end: FreeHMMstruct() ***/
812
813
814
/**
815
* @brief comparisin function for ascending qsort() of doubles
816
*/
817
int
818
CompDblAsc(const void *pv1, const void *pv2){
819
820
double d1 = *(double *)pv1;
821
double d2 = *(double *)pv2;
822
823
if (d1 > d2) { return +1; }
824
else if (d1 < d2) { return -1; }
825
else { return 0; }
826
827
} /*** end: CompDblAsc() ***/
828
829
830
/**
831
* AlnToHMM2()
832
*
833
* @brief convert alignment into HMM (hhmake)
834
*
835
* @param[out] prHMM
836
* structure with pseudocounts etc
837
* @param[in] pcHMM_input
838
* name of file with HMM
839
*
840
*/
841
extern "C" int
842
AlnToHMM2(hmm_light *rHMM_p,
843
char **ppcSeq, int iN){
844
845
if (0 == par.M){
846
SetDefaults();
847
SetSubstitutionMatrix();
848
par.cons = 1;
849
par.M = 2;
850
par.forward=2;
851
par.Mgaps=100;
852
const int ciGoodMeasureSeq = 10;
853
int iAuxLen = strlen(ppcSeq[0]);
854
par.maxResLen = iAuxLen;
855
par.maxResLen += ciGoodMeasureSeq;
856
par.maxColCnt = iAuxLen + ciGoodMeasureSeq;
857
par.max_seqid=DEFAULT_FILTER;
858
par.loc=0; par.mact=0;
859
/* some minor parameters, affecting alignment (i think) */
860
par.p = 0.0; /* minimum threshold for inclusion in hit list */
861
par.Z = 100; /* minimum threshold for inclusion in hit list and alignment listing */
862
par.z = 1; /* min number of lines in hit list */
863
par.B = 100; /* max number of alignments */
864
par.b = 1; /* min number of alignments */
865
par.E = 1e6; /* maximum threshold for inclusion in hit list and alignment listing */
866
par.altali=1;par.hitrank=0;par.showcons=1; par.showdssp=1;par.showpred=1;par.nseqdis=iN;par.cons=1;
867
}
868
869
const int ciGoodMeasure = 10;
870
int iLen = strlen(ppcSeq[0]) + ciGoodMeasure;
871
if (iLen > par.maxResLen){
872
par.maxResLen = par.maxColCnt = iLen;
873
}
874
HMM rTemp(iN+2, iLen); /* temporary template */
875
Alignment rTempAli(iN+2, iLen); /* temporary alignment */
876
int iParCons = par.cons;
877
878
/*printf(">>>>>>>>>>> %s:%s:%d: there are %d sequences\n", __FUNCTION__, __FILE__, __LINE__, iN);*/
879
880
par.cons = 1;
881
if (OK != ReadAndPrepare(INTERN_ALN_2_HMM,
882
ppcSeq, iN, rHMM_p,
883
NULL/*prealigned*/, NULL/*representative*/, NULL/*weights*/, //YES/*transfer*/,
884
(char*)("")/*in-file*/, rTemp, &rTempAli)) {
885
return FAILURE;
886
}
887
par.cons = iParCons;
888
889
/******/
890
/*** transfer class info to struct */
891
rHMM_p->n_display = rTemp.n_display;
892
rHMM_p->sname = NULL;
893
rHMM_p->seq = (char **)calloc((rTemp.n_display+1), sizeof(char *));
894
895
for (int i = 0; i < rTemp.n_display; i++){
896
/* NOTE: In the original hhalign code the first position
897
is kept open ('\0'). This makes it difficult to use
898
string functions like strlen/strdup etc.
899
Insert a temporary gap (.) to facilitate string operations */
900
char cHead = rTemp.seq[i][0];
901
rTemp.seq[i][0] = '.';
902
rHMM_p->seq[i] = strdup(rTemp.seq[i]);
903
rTemp.seq[i][0] = rHMM_p->seq[i][0] = cHead;
904
}
905
rHMM_p->ncons = rTemp.ncons;
906
rHMM_p->nfirst = rTemp.nfirst;
907
if (-1 == rHMM_p->ncons){
908
/* ncons needed later for alignment of
909
representative sequence and copy of profile.
910
ncons not always set (-1 default),
911
this will cause segmentation fault.
912
nfirst (probably) right index -
913
no problem if not */
914
rHMM_p->ncons = rTemp.nfirst;
915
}
916
rHMM_p->nss_dssp = rTemp.nss_dssp;
917
rHMM_p->nsa_dssp = rTemp.nsa_dssp;
918
rHMM_p->nss_pred = rTemp.nss_pred;
919
rHMM_p->nss_conf = rTemp.nss_conf;
920
rHMM_p->L = rTemp.L;
921
rHMM_p->N_in = rTemp.N_in;
922
rHMM_p->N_filtered = rTemp.N_filtered;
923
#define NEFF_SCORE 10 /* FIXME: Magic Number */ //// get rid of that, FS, 2010-12-22
924
rHMM_p->Neff_M = (float *)calloc(rHMM_p->L+1, sizeof(float));
925
rHMM_p->Neff_I = (float *)calloc(rHMM_p->L+1, sizeof(float));
926
rHMM_p->Neff_D = (float *)calloc(rHMM_p->L+1, sizeof(float));
927
for (int i = 0; i < rHMM_p->L+1; i++){
928
rHMM_p->Neff_M[i] = rHMM_p->Neff_I[i] = rHMM_p->Neff_D[i] = NEFF_SCORE; //// get rid of that, FS, 2010-12-22
929
}
930
rHMM_p->Neff_HMM = rTemp.Neff_HMM;
931
/* skip longname,name,file,fam,sfam,fold,cl */
932
rHMM_p->lamda = rTemp.lamda;
933
rHMM_p->mu = rTemp.mu;
934
935
HMMshadow rShadow = {0}; /* friend of HMM to access private members */
936
rShadow.copyHMMtoShadow(rTemp);
937
938
rHMM_p->Neff_M = (float *)calloc(rHMM_p->L+1, sizeof(float));
939
rHMM_p->Neff_I = (float *)calloc(rHMM_p->L+1, sizeof(float));
940
rHMM_p->Neff_D = (float *)calloc(rHMM_p->L+1, sizeof(float));
941
rHMM_p->f = (float **)calloc(rHMM_p->L+1, sizeof(float *));
942
rHMM_p->g = (float **)calloc(rHMM_p->L+1, sizeof(float *));
943
rHMM_p->p = (float **)calloc(rHMM_p->L+1, sizeof(float *));
944
rHMM_p->tr = (float **)calloc(rHMM_p->L+1, sizeof(float *));
945
rHMM_p->linTr = (float **)calloc(rHMM_p->L+1, sizeof(float *));
946
947
for (int i = 0; i < rHMM_p->L+1; i++){
948
rHMM_p->Neff_M[i] = (float)rShadow.Neff_M[i];
949
rHMM_p->Neff_I[i] = (float)rShadow.Neff_I[i];
950
rHMM_p->Neff_D[i] = (float)rShadow.Neff_D[i];
951
rHMM_p->f[i] = (float *)calloc(AMINOACIDS, sizeof(float));
952
rHMM_p->g[i] = (float *)calloc(AMINOACIDS, sizeof(float));
953
rHMM_p->p[i] = (float *)calloc(AMINOACIDS, sizeof(float));
954
for (int j = 0; j < AMINOACIDS; j++){
955
rHMM_p->f[i][j] = (float)rShadow.f[i][j];
956
rHMM_p->g[i][j] = (float)rShadow.g[i][j];
957
rHMM_p->p[i][j] = (float)rShadow.p[i][j];
958
}
959
rHMM_p->tr[i] = (float *)calloc(STATE_TRANSITIONS, sizeof(float));
960
rHMM_p->linTr[i] = (float *)calloc(STATE_TRANSITIONS, sizeof(float));
961
for (int j = 0; j< STATE_TRANSITIONS; j++){
962
rHMM_p->tr[i][j] = (float)rShadow.tr[i][j];
963
rHMM_p->linTr[i][j] = fpow2(rShadow.tr[i][j]);
964
}
965
} /*0 <= i < prHMM->L+1 */
966
/* skip trans_lin,ss_dssp,sa_dssp,ss_pred,ss_conf,Xcons */
967
for (int j = 0; j < AMINOACIDS; j++){
968
rHMM_p->pav[j] = (float)rShadow.pav[j];
969
}
970
/* skip pnul */
971
972
973
/******/
974
975
976
rTemp.ClobberGlobal();
977
rTempAli.ClobberGlobal();
978
979
return OK;
980
981
} /*** end of AlnToHMM2() ***/
982
983
984
/**
985
* HHMake_Wrapper()
986
*
987
* @brief turn alignment (from file) into HHM (HMM) on file
988
*
989
* @param[out] hmm_out
990
* name of file with HMM info corresponding to tmp_aln
991
* @param[in] tmp_aln
992
* name of file with alignment, to be turned into HMM (HHM)
993
*
994
*/
995
extern "C" int
996
HHMake_Wrapper(char *tmp_aln, char *hmm_out)
997
{
998
999
HMM rTemp; /* temporary template */
1000
Alignment rTempAli; /* temporary alignment */
1001
hmm_light *rHMM_p = NULL;
1002
1003
/** Note:
1004
this is a wrapper for a stand-alone program hhmake.
1005
hhmake uses a different set of parameters than hhalign.
1006
However, all parameters are GLOBAL.
1007
So at this point we register the hhalign settings,
1008
reset them to hhmake settings and set them back
1009
at the end of the function
1010
*/
1011
1012
/* save settings from hhalign */
1013
int iParShowcons=par.showcons;
1014
int iParAppend=par.append;
1015
int iParNseqdis=par.nseqdis;
1016
int iParMark=par.mark;
1017
int iParMaxSeqid=par.max_seqid;
1018
int iParQid=par.qid;
1019
double dParQsc=par.qsc;
1020
int iParCoverage=par.coverage;
1021
int iParNdiff=par.Ndiff;
1022
int iParCoverageCore=par.coverage_core;
1023
double dParQscCore=par.qsc_core;
1024
double dParCoresc=par.coresc;
1025
int iParM=par.M;
1026
int iParMgaps=par.Mgaps;
1027
int iParMatrix=par.matrix;
1028
int iParPcm=par.pcm;
1029
double dParPcw=par.pcw;
1030
double dParGapb=par.gapb;
1031
int iParWg=par.wg;
1032
1033
/* these are settings suitable for hhmake */
1034
par.showcons=1; // write consensus sequence into hhm file
1035
par.append=0; // overwrite output file
1036
par.nseqdis=10; // maximum number of query or template sequences to be recoreded in HMM and diplayed in output alignments
1037
par.mark=0; // 1: only marked sequences (or first) get displayed; 0: most divergent ones get displayed
1038
par.max_seqid=90; // default for maximum sequence identity threshold
1039
par.qid=0; // default for maximum sequence identity threshold
1040
par.qsc=-20.0f; // default for minimum score per column with query
1041
par.coverage=0; // default for minimum coverage threshold
1042
par.Ndiff=100; // pick Ndiff most different sequences from alignment
1043
par.coverage_core=30; // Minimum coverage for sequences in core alignment
1044
par.qsc_core=0.3f; // Minimum score per column with core alignment (HMM)
1045
par.coresc=-20.0f; // Minimum score per column with core alignment (HMM)
1046
par.M=1; // match state assignment is by A2M/A3M
1047
par.Mgaps=50; // above this percentage of gaps, columns are assigned to insert states
1048
par.matrix=0; // Subst.matrix 0: Gonnet, 1: HSDM, 2: BLOSUM50 3: BLOSUM62
1049
par.pcm=0; // no amino acid and transition pseudocounts added
1050
par.pcw=0; // wc>0 weighs columns according to their intra-clomun similarity
1051
par.gapb=0.0; // default values for transition pseudocounts; 0.0: add no transition pseudocounts!
1052
par.wg=0; // 0: use local sequence weights 1: use local ones
1053
1054
1055
if (OK != ReadAndPrepare(READ_ALN_2_HMM,
1056
NULL, 0, rHMM_p, NULL, NULL, NULL,
1057
tmp_aln, rTemp, &rTempAli)) {
1058
return FAILURE;
1059
}
1060
1061
rTemp.WriteToFile(hmm_out);
1062
1063
1064
1065
/* restore settings from hhalign */
1066
par.showcons=iParShowcons;
1067
par.append=iParAppend;
1068
par.nseqdis=iParNseqdis;
1069
par.mark=iParMark;
1070
par.max_seqid=iParMaxSeqid;
1071
par.qid=iParQid;
1072
par.qsc=dParQsc;
1073
par.coverage=iParCoverage;
1074
par.Ndiff=iParNdiff;
1075
par.coverage_core=iParCoverageCore;
1076
par.qsc_core=dParQscCore;
1077
par.coresc=dParCoresc;
1078
par.M=iParM;
1079
par.Mgaps=iParMgaps;
1080
par.matrix=iParMatrix;
1081
par.pcm=iParPcm;
1082
par.pcw=dParPcw;
1083
par.gapb=dParGapb;
1084
par.wg=iParWg;
1085
1086
1087
/* prepare exit */
1088
rTemp.ClobberGlobal();
1089
rTempAli.ClobberGlobal();
1090
1091
return OK;
1092
}
1093
1094
/**
1095
* readHMMWrapper()
1096
*
1097
* @brief read HMM from file, copy (relevant) info into struct
1098
*
1099
* @param[out] prHMM
1100
* structure with pseudocounts etc, probably uninitialised on entry
1101
* @param[in] pcHMM_input
1102
* name of file with HMM
1103
*
1104
*/
1105
extern "C" int
1106
readHMMWrapper(hmm_light *rHMM_p,
1107
char *pcHMM_input)
1108
{
1109
par.maxResLen = 15002;
1110
HMM rTemp(1000,par.maxResLen); /* temporary template */
1111
Alignment rTempAli; /* temporary alignment */
1112
/* AW changed init from {0} to 0 because it failed to compile with
1113
* gcc 4.3.3 with the following error:
1114
* error: braces around initializer for non-aggregate type
1115
*/
1116
/* FS taken out initialiser alltogether */
1117
1118
/* 0th arg of RnP is the type of RnP, ie,
1119
enum {INTERN_ALN_2_HMM = 0, READ_ALN_2_HMM, READ_HMM_2_HMM};*/
1120
/* 1st arg of ReadAndPrepare() is ppcSeqs, 2nd is #seq */
1121
/* FIXME: at this stage the 3rd arg, rHMM_p, only acts as a dummy.
1122
this is rather silly, as it is the actual struct that will
1123
carry the HMM info at the end.
1124
If we write it already in ReadAndPrepare() then we could
1125
dispense with all this friend-class nonsense */
1126
if (OK != ReadAndPrepare(READ_HMM_2_HMM,
1127
NULL, 0, rHMM_p, NULL, NULL, NULL,
1128
pcHMM_input, rTemp, &rTempAli)) {
1129
return FAILURE;
1130
}
1131
1132
/* an important piece of information I want to get out of here
1133
is the consenssus sequence. there are however certain
1134
Pfam HMMs that don't trigger consensus calculation.
1135
at the moment I (FS) don't understand why this is
1136
(or rather why this is not). The proper place to do this
1137
should be inside ReadAndPrepare():ReadHMMer(), but
1138
I have not quite penetrated the logic there.
1139
Therefore I try to catch this condition at this point (here)
1140
and rectify it.
1141
*/
1142
if (-1 == rHMM_p->ncons){
1143
int i, iA;
1144
rHMM_p->ncons = rHMM_p->nfirst;
1145
1146
for (i = 0; i < rHMM_p->L; i++){
1147
double dPmax = 0.00;
1148
int iAmax = -1;
1149
for (iA = 0; iA < AMINOACIDS; iA++){
1150
if (rHMM_p->f[i][iA] > dPmax){
1151
iAmax = iA;
1152
dPmax = rHMM_p->f[i][iA];
1153
}
1154
} /* (0 <= iA < AMINOACIDS) */
1155
rHMM_p->seq[rHMM_p->ncons][i] = i2aa(iAmax);
1156
} /* (0 <= i < rHMM_p->L) */
1157
1158
} /* ncons not set */
1159
1160
1161
rTemp.ClobberGlobal();
1162
rTempAli.ClobberGlobal();
1163
1164
return OK;
1165
1166
} /*** end: readHMMWrapper() ***/
1167
1168
1169
1170
1171
1172
1173
/////////////////////////////////////////////////////////////////////////////
1174
/**
1175
* @brief Do precalculations for q and t to prepare comparison
1176
*/
1177
void
1178
PrepareTemplate(HMM& q, HMM& t, int format)
1179
{
1180
if (format==0) // HHM format
1181
{
1182
// Add transition pseudocounts to template
1183
t.AddTransitionPseudocounts();
1184
1185
// Generate an amino acid frequency matrix from t.f[i][a] with max pseudocounts (tau=1) -> t.g[i][a]
1186
t.PreparePseudocounts();
1187
1188
// Add amino acid pseudocounts to template: t.p[i][a] = (1-tau)*f[i][a] + tau*g[i][a]
1189
t.AddAminoAcidPseudocounts();
1190
}
1191
else // HHMER format
1192
{
1193
// Don't add transition pseudocounts to template
1194
// t.AddTransitionPseudocounts(par.gapd, par.gape, par.gapf, par.gapg, par.gaph, par.gapi, 0.0);
1195
1196
// Generate an amino acid frequency matrix from f[i][a] with full pseudocount admixture (tau=1) -> g[i][a]
1197
t.PreparePseudocounts();
1198
1199
// DON'T ADD amino acid pseudocounts to temlate: pcm=0! t.p[i][a] = t.f[i][a]
1200
t.AddAminoAcidPseudocounts(0, par.pca, par.pcb, par.pcc);
1201
}
1202
1203
// Modify transition probabilities to include SS-dependent penalties
1204
if (par.ssgap) t.UseSecStrucDependentGapPenalties();
1205
1206
if (par.forward>=1) t.Log2LinTransitionProbs(1.0);
1207
1208
// Factor Null model into HMM t
1209
// ATTENTION! t.p[i][a] is divided by pnul[a] (for reasons of efficiency) => do not reuse t.p
1210
t.IncludeNullModelInHMM(q,t); // Can go BEFORE the loop if not dependent on template
1211
1212
return;
1213
}
1214
1215
1216
/*** end of hhfunc-C.h ***/
Loggerhead is a web-based interface for Breezy
Version: 2.0.1