1 /* @source cons application
2 **
3 ** Calculates a consensus
4 ** @author Copyright (C) Tim Carver (tcarver@hgmp.mrc.ac.uk)
5 ** @@
6 **
7 **
8 ** -plurality - defines no. of +ve scoring matches below
9 ** which there is no consensus.
10 **
11 ** -identity - defines the number of identical symbols
12 ** requires in an alignment column for it to
13 ** included in the consensus.
14 **
15 ** -setcase - upper/lower case given if score above/below
16 ** user defined +ve matching threshold.
17 **
18 **
19 ** This program is free software; you can redistribute it and/or
20 ** modify it under the terms of the GNU General Public License
21 ** as published by the Free Software Foundation; either version 2
22 ** of the License, or (at your option) any later version.
23 **
24 ** This program is distributed in the hope that it will be useful,
25 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
26 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27 ** GNU General Public License for more details.
28 **
29 ** You should have received a copy of the GNU General Public License
30 ** along with this program; if not, write to the Free Software
31 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
32 ******************************************************************************/
33
34 #include "emboss.h"
35
36
37
38
39 /* @prog cons *****************************************************************
40 **
41 ** Creates a consensus from multiple alignments
42 **
43 ******************************************************************************/
44
main(int argc,char ** argv)45 int main(int argc, char **argv)
46 {
47 ajint nseqs;
48 ajuint mlen;
49 ajint i;
50 ajint identity;
51 float fplural;
52 float setcase;
53 const char *p;
54 AjPSeqset seqset;
55 AjPSeqout seqout;
56 AjPSeq seqo;
57 AjPStr name = NULL;
58 AjPStr cons;
59 AjPMatrix cmpmatrix = 0;
60
61
62 embInit ("cons", argc, argv);
63
64 seqset = ajAcdGetSeqset ("sequence");
65 cmpmatrix = ajAcdGetMatrix("datafile");
66 fplural = ajAcdGetFloat("plurality");
67 setcase = ajAcdGetFloat("setcase");
68 identity = ajAcdGetInt("identity");
69 seqout = ajAcdGetSeqout("outseq");
70 name = ajAcdGetString ("name");
71
72 nseqs = ajSeqsetGetSize(seqset);
73 if(nseqs<2)
74 ajFatal("Insufficient sequences (%d) to create a matrix",nseqs);
75
76 mlen = ajSeqsetGetLen(seqset);
77 for(i=0;i<nseqs;++i) /* check sequences are same length */
78 {
79 p = ajSeqsetGetseqSeqC(seqset,i);
80 if(strlen(p)!=mlen)
81 {
82 ajWarn("Sequence lengths are not equal: expect %d '%S' is %d",
83 mlen, ajSeqsetGetseqNameS(seqset, i), p);
84 break;
85 }
86 }
87
88 ajSeqsetFmtUpper(seqset);
89
90 cons = ajStrNewRes(mlen+1);
91 embConsCalc (seqset, cmpmatrix, nseqs, mlen,
92 fplural, setcase, identity, ajFalse, &cons);
93
94 /* write out consensus sequence */
95 seqo = ajSeqNew();
96 ajSeqAssignSeqS(seqo,cons);
97
98 if (ajSeqsetIsNuc(seqset))
99 ajSeqSetNuc(seqo);
100 else
101 ajSeqSetProt(seqo);
102
103 if(name == NULL)
104 ajSeqAssignNameS(seqo,ajSeqsetGetNameS(seqset));
105 else
106 ajSeqAssignNameS(seqo,name);
107
108 ajSeqoutWriteSeq(seqout,seqo);
109 ajSeqoutClose(seqout);
110
111 ajStrDel(&cons);
112 ajSeqDel(&seqo);
113 ajSeqsetDel(&seqset);
114 ajMatrixDel(&cmpmatrix);
115 ajSeqoutDel(&seqout);
116 ajStrDel(&name);
117
118 embExit();
119
120 return 0;
121 }
122