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 consambig *************************************************************
40 **
41 ** Creates an ambiguous 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 const char *p;
51 AjPSeqset seqset;
52 AjPSeqout seqout;
53 AjPSeq seqo;
54 AjPStr name = NULL;
55 AjPStr cons;
56
57 embInit ("consambig", argc, argv);
58
59 seqset = ajAcdGetSeqset ("sequence");
60 seqout = ajAcdGetSeqout("outseq");
61 name = ajAcdGetString ("name");
62
63 nseqs = ajSeqsetGetSize(seqset);
64 if(nseqs<2)
65 ajFatal("Insufficient sequences (%d) to create a matrix",nseqs);
66
67 mlen = ajSeqsetGetLen(seqset);
68 for(i=0;i<nseqs;++i) /* check sequences are same length */
69 {
70 p = ajSeqsetGetseqSeqC(seqset,i);
71 if(strlen(p)!=mlen)
72 {
73 ajWarn("Sequence lengths are not equal: expect %d '%S' is %d",
74 mlen, ajSeqsetGetseqNameS(seqset, i), p);
75 break;
76 }
77 }
78
79 ajSeqsetFmtUpper(seqset);
80
81 cons = ajStrNewRes(mlen+1);
82
83 ajAlignConsAmbig (seqset, &cons);
84
85 /* write out consensus sequence */
86 seqo = ajSeqNew();
87 ajSeqAssignSeqS(seqo,cons);
88
89 if (ajSeqsetIsNuc(seqset))
90 ajSeqSetNuc(seqo);
91 else
92 ajSeqSetProt(seqo);
93
94 if(name == NULL)
95 ajSeqAssignNameS(seqo,ajSeqsetGetNameS(seqset));
96 else
97 ajSeqAssignNameS(seqo,name);
98
99 ajSeqoutWriteSeq(seqout,seqo);
100 ajSeqoutClose(seqout);
101
102 ajStrDel(&cons);
103 ajSeqDel(&seqo);
104 ajSeqsetDel(&seqset);
105 ajSeqoutDel(&seqout);
106 ajStrDel(&name);
107
108 embExit();
109
110 return 0;
111 }
112