1 /* @source btwisted application
2 **
3 ** Calculates twist in B DNA
4 ** @author Copyright (C) David Martin (dmartin@hgmp.mrc.ac.uk)
5 ** @author Copyright (C) Alan Bleasby (ableasby@hgmp.mrc.ac.uk)
6 ** @@
7 **
8 ** This program is free software; you can redistribute it and/or
9 ** modify it under the terms of the GNU General Public License
10 ** as published by the Free Software Foundation; either version 2
11 ** of the License, or (at your option) any later version.
12 **
13 ** This program is distributed in the hope that it will be useful,
14 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
15 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 ** GNU General Public License for more details.
17 **
18 ** You should have received a copy of the GNU General Public License
19 ** along with this program; if not, write to the Free Software
20 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
21 ******************************************************************************/
22
23 #include "emboss.h"
24
25
26
27
28 static AjPTable btwisted_getdinucdata(AjPFile anglesfile);
29
30
31
32
33 /* @prog btwisted *************************************************************
34 **
35 ** Calculates the twisting in a B-DNA sequence
36 **
37 ******************************************************************************/
38
main(int argc,char ** argv)39 int main(int argc, char **argv)
40 {
41
42 AjPSeq seq = NULL;
43 AjPFile angles = NULL;
44 AjPFile energies = NULL;
45 AjPFile result = NULL;
46
47 AjPTable angletable = NULL;
48 AjPTable energytable = NULL;
49
50 AjPStr nucs = NULL;
51 const AjPStr valstr = NULL;
52
53 const char * dinuc = NULL;
54 ajint len;
55 ajint begin;
56 ajint end;
57 ajint i;
58 float val;
59 float anglesum = 0.0;
60 float energysum = 0.0;
61 float twists = 0.0;
62
63 float basesperturn = 0.0;
64 float energyperbase = 0.0;
65
66 embInit ("btwisted", argc, argv);
67
68 seq = ajAcdGetSeq ("sequence");
69 angles = ajAcdGetDatafile("angledata");
70 energies = ajAcdGetDatafile("energydata");
71 result = ajAcdGetOutfile ("outfile");
72
73
74 nucs = ajStrNew();
75
76 angletable = btwisted_getdinucdata(angles);
77 energytable = btwisted_getdinucdata(energies);
78
79 ajFileClose(&angles);
80 ajFileClose(&energies);
81
82 begin = ajSeqGetBegin(seq);
83 end = ajSeqGetEnd(seq);
84
85 len = end-begin+1;
86
87 dinuc = ajSeqGetSeqC(seq);
88
89 for(i=begin-1; i<end-1; ++i)
90 {
91 ajStrAssignSubC(&nucs,dinuc,i,i+1);
92 if(!(valstr = ajTableFetchS(angletable, nucs)))
93 ajFatal("Incomplete table");
94
95 ajStrToFloat(valstr,&val);
96 anglesum += val;
97 if(!(valstr = ajTableFetchS(energytable, nucs)))
98 ajFatal("Incomplete table");
99
100 ajStrToFloat(valstr,&val);
101 energysum += val;
102 }
103
104 twists = anglesum / (float)360.0 ;
105 basesperturn = (float) len * (float)360.0 /anglesum;
106 energyperbase = energysum/(float) (len-1);
107
108 ajFmtPrintF(result, "# Output from BTWISTED\n");
109 ajFmtPrintF(result, "# Twisting calculated from %d to %d of %s\n",
110 begin, end, ajSeqGetNameC(seq));
111 ajFmtPrintF(result,"Total twist (degrees): %.1f\n", anglesum);
112 ajFmtPrintF(result,"Total turns : %.2f\n", twists);
113 ajFmtPrintF(result,"Average bases per turn: %.2f\n", basesperturn);
114 ajFmtPrintF(result,"Total stacking energy : %.2f\n", energysum);
115 ajFmtPrintF(result,"Average stacking energy per dinucleotide: %.2f\n",
116 energyperbase);
117
118 ajTablestrFree(&angletable);
119 ajTablestrFree(&energytable);
120
121 ajStrDel(&nucs);
122 ajFileClose(&result);
123 ajSeqDel(&seq);
124
125 embExit ();
126
127 return 0;
128 }
129
130
131
132
133 /* @funcstatic btwisted_getdinucdata ******************************************
134 **
135 ** Undocumented.
136 **
137 ** @param [u] inf [AjPFile] Data file
138 ** @return [AjPTable] Data values table
139 ** @@
140 ******************************************************************************/
141
btwisted_getdinucdata(AjPFile inf)142 static AjPTable btwisted_getdinucdata(AjPFile inf)
143 {
144 AjPStr valstr = NULL;
145 AjPStr key = NULL;
146 AjPStr line = NULL;
147 AjPStrTok token = NULL;
148 AjPTable table = NULL;
149
150 valstr = ajStrNew();
151 line = ajStrNew();
152
153 table = ajTablestrNewCase(20);
154
155 while(ajReadlineTrim(inf,&line))
156 {
157 if(*ajStrGetPtr(line)=='#')
158 continue;
159 token = ajStrTokenNewC(line," \n\t\r");
160 key = ajStrNew();
161 ajStrTokenNextParseC(token," \n\t\r",&key);
162 valstr = ajStrNew();
163 ajStrTokenNextParseC(token," \n\t\r",&valstr);
164 ajTablePut(table,(void *)key,(void *) valstr);
165 ajStrTokenDel(&token);
166 }
167
168
169 ajStrDel(&line);
170
171 return table;
172 }
173