1 /* @source isochore application
2 **
3 ** Finds isochores
4 **
5 ** @author Copyright (C) Peter Rice
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 /* @datastatic PIntarr ******************************************************
29 **
30 ** Integer array
31 **
32 ** @alias SIntarr
33 ** @alias PIntarr
34 **
35 ** @attr Array [ajint*] Integer array
36 ** @attr Size [ajint] Size
37 ** @attr Padding [char[4]] Padding to alignment boundary
38 ******************************************************************************/
39
40 typedef struct SIntarr
41 {
42 ajint* Array;
43 ajint Size;
44 char Padding[4];
45 } OIntarr;
46 #define PIntarr OIntarr*
47
48
49
50
51 /* @datastatic IsochorePFltarr ************************************************
52 **
53 ** Integer array
54 **
55 ** @alias IsochoreSFltarr
56 ** @alias IsochorePFltarr
57 **
58 ** @attr Array [float*] Float array
59 ** @attr Size [ajint] Size
60 ** @attr Padding [char[4]] Padding to alignment boundary
61 ******************************************************************************/
62
63 typedef struct IsochoreSFltarr
64 {
65 float* Array;
66 ajint Size;
67 char Padding[4];
68 } IsochoreOFltarr;
69 #define IsochorePFltarr IsochoreOFltarr*
70
71
72
73 static void isochore_FltarrDel(IsochorePFltarr *parr);
74 static IsochorePFltarr isochore_FltarrNew0(size_t size);
75
76
77
78
79 /* @prog isochore *************************************************************
80 **
81 ** Calculates the G+C content of a DNA sequence
82 ** by sliding a window of size "iwin" in increments of "ishift" bases
83 ** at a time.
84 ** Results are stored in float array "results" with one position for
85 ** each calculated value.
86 **
87 ** Results are also written to an output file with tab delimiters.
88 **
89 ** To plot this, start at the centre of the first window (0 + iwin/2)
90 ** and plot a point from results0>Array[0] every "ishift" bases until
91 ** the end (use isize).
92 **
93 ** Future changes: Users should be able to ask for a sequence range
94 ** and plot just that range. Currently a range can be set on the
95 ** command line but it is ignored. The range is from "ajSeqGetBegin(seq)"
96 ** to "ajSeqGetEnd(seq)".
97 **
98 ******************************************************************************/
99
100
main(int argc,char ** argv)101 int main(int argc, char **argv)
102 {
103 AjPSeq seq;
104 AjPFile out;
105 IsochorePFltarr results;
106
107 AjPGraph plot;
108 AjPGraphdata graphdata;
109
110 ajint iwin;
111 ajint ishift;
112 ajint i;
113 ajint j;
114 ajint k;
115 ajint ipos;
116 ajint isize;
117 const char *sq;
118 ajint igc;
119 ajint imax;
120 ajint ibeg;
121 ajint iend;
122 ajint ilen;
123 float amin = 0.;
124 float amax = 0.;
125
126 embInit("isochore", argc, argv);
127
128 seq = ajAcdGetSeq("sequence");
129 out = ajAcdGetOutfile("outfile");
130 plot = ajAcdGetGraphxy("graph");
131
132 ajGraphInitSeq(plot, seq);
133
134 sq = ajStrGetPtr(ajSeqGetSeqS(seq));
135
136 ibeg = ajSeqGetBegin(seq);
137 iend = ajSeqGetEnd(seq);
138 ilen = ajSeqGetLen(seq);
139
140 iwin = ajAcdGetInt("window");
141 ishift = ajAcdGetInt("shift");
142
143 imax = iend + iwin/2; /* stop at imax */
144 if(imax > ilen)
145 imax = ilen;
146
147 i = ibeg - iwin/2; /* start calculating from i */
148 if(i < 0)
149 i = 0;
150
151 isize = 1 + (imax - iwin - i)/ishift; /* size of results array */
152 results = isochore_FltarrNew0(isize);
153
154 ajDebug("ilen: %d ibeg: %d iend: %d\n", ilen, ibeg, iend);
155 ajDebug("iwin: %d ishift: %d isize: %d imax: %d i: %d\n",
156 iwin, ishift, isize, imax, i);
157
158 ipos = i + iwin/2;
159 ajFmtPrintF(out, "Position\tPercent G+C %d .. %d\n",
160 ibeg, iend);
161
162 for(j=0; j < isize; i+=ishift, j++)
163 { /* sum over window */
164 igc = 0;
165 for(k=0; k < iwin; k++)
166 if(strchr("CcGg", sq[i+k]))
167 igc++;
168
169 results->Array[j] = (float) igc / (float) iwin;
170 ajFmtPrintF(out, "%d\t%.3f\n", ipos, results->Array[j]);
171 ipos += ishift;
172 }
173
174
175 ajFileClose(&out);
176
177
178 #ifndef NO_PLOT
179 i = ibeg - iwin/2; /* start calculating from i */
180 if(i < 0)
181 i = 0;
182 ipos = i + iwin/2;
183
184 /* create the graph */
185
186 graphdata = ajGraphdataNew();
187
188 ajGraphicsCalcRange(results->Array,isize,&amin,&amax);
189
190 ajGraphdataSetTruescale(graphdata,(float)ipos,(float)(ipos+(ishift*isize)),
191 amin,amax);
192 ajGraphdataSetMinmax(graphdata,(float)ipos,(float)(ipos+(ishift*isize)),
193 amin,amax);
194 ajGraphdataSetTypeC(graphdata,"2D Plot");
195 ajGraphdataSetTitleC(graphdata,"");
196
197
198
199 ajGraphDataAdd(plot,graphdata);
200 ajGraphdataCalcXY(graphdata, isize,(float)(ipos),(float)ishift,
201 results->Array);
202
203
204 /* display the region 0 -> 1 for the y axis */
205 ajGraphxySetYstartF(plot,0.0);
206 ajGraphxySetYendF(plot,1.0);
207
208 /* draw the graph */
209 ajGraphxyDisplay(plot,AJTRUE);
210
211 /* Delete the structures and data */
212 ajGraphxyDel(&plot);
213 #endif
214
215
216
217 /*
218 ** plot is an XY graph definition object created by acdSetGraphxy
219 ** something like this to plot the data:
220 ** sequence, start position, increment, array, array size
221 ** Note: seq has Begin and End values which can limit the plot range
222 **
223 ** ioff = ibeg + iwin/2;
224 ** ajPlotInit(plot, seq);
225 ** ajPlotFloat(plot, ioff, ishift, results->Array, isize);
226 **
227 */
228
229 ajSeqDel(&seq);
230 ajFileClose(&out);
231 ajGraphxyDel(&plot);
232 isochore_FltarrDel(&results);
233 embExit();
234
235 return 0;
236 }
237
238
239
240
241
242 /* @funcstatic isochore_FltarrNew0 ********************************************
243 **
244 ** Undocumented.
245 **
246 ** @param [r] size [size_t] Undocumented
247 ** @return [IsochorePFltarr] Undocumented
248 ** @@
249 ******************************************************************************/
250
isochore_FltarrNew0(size_t size)251 static IsochorePFltarr isochore_FltarrNew0(size_t size)
252 {
253 IsochorePFltarr ret;
254
255 AJNEW(ret);
256 ret->Size = (ajint) size;
257 AJCNEW(ret->Array,size);
258
259 return ret;
260 }
261
262
263
264
265 /* @funcstatic isochore_FltarrDel *********************************************
266 **
267 ** Undocumented.
268 **
269 ** @param [d] parr [IsochorePFltarr*] Undocumented
270 ** @return [void]
271 ** @@
272 ******************************************************************************/
273
isochore_FltarrDel(IsochorePFltarr * parr)274 static void isochore_FltarrDel(IsochorePFltarr *parr)
275 {
276 if(!*parr) return;
277 AJFREE((*parr)->Array);
278 AJFREE(*parr);
279
280 return;
281 }
282