← Back to branch summary

~ubuntu-branches/ubuntu/trusty/tigr-glimmer/trusty

« back to all changes in this revision

Viewing changes to SimpleMake/anomaly.cc

  • Committer: Bazaar Package Importer
  • Author(s): Andreas Tille, Charles Plessy, Andreas Tille
  • Date: 2008-04-22 11:59:07 UTC
  • mfrom: (3.1.1 lenny)
  • Revision ID: james.westby@ubuntu.com-20080422115907-rpsy3kgts0vu6lxq
Tags: 3.02-1
http://bugs.debian.org/385258
[ Charles Plessy ]
* debian/watch:
  - Replaced by the new one written by Nelson (Closes: #385258)

[ Andreas Tille ]
*  New upstream version
* Group maintenance by Debian-Med team
  - DM-Upload-Allowed: Yes
  - Vcs tags
  - Use correct address as Uploader: Steffen Moeller <moeller@debian.org>
* Standards-Version: 3.7.3 (no changes needed)
* debhelper >= 5
* Moved Homepage from long description to control fields
* Removed [Biology] from short description

Show diffs side-by-side

added added

removed removed

Lines of Context:
SimpleMake/anomaly.cc
 
1
//  A. L. Delcher
 
2
//
 
3
//  File:  anomaly.cc
 
4
//
 
5
//  Last Modified:  Fri May 19 17:10:05 EDT 2006
 
6
//
 
7
//  This program reads the sequence in the first command-line
 
8
//  file and then takes the start and end positions specified in
 
9
//  the second command-line file and checks for anomalous
 
10
//  start/stop codons and frame shifts.
 
11
 
 
12
 
 
13
#include  "anomaly.hh"
 
14
 
 
15
 
 
16
// Global variables
 
17
 
 
18
static bool  Check_Previous_Stop = false;
 
19
  // Determines whether to check if codon before start coordinate
 
20
  // is a valid stop codon
 
21
static bool  Check_Start_Codon = true;
 
22
  // Determines whether to check if first codon is a valid start
 
23
static char  * Coord_File_Name = NULL;
 
24
  // From the second command-line argument
 
25
static int  Num_Start_Codons;
 
26
  // Number of different start codon patterns
 
27
static int  Num_Stop_Codons;
 
28
  // Number of different stop codon patterns
 
29
static char  * Sequence_File_Name = NULL;
 
30
  // From the first command-line argument
 
31
static vector <const char *>  Start_Codon;
 
32
  // Sequences assumed to be start codons
 
33
static vector <const char *>  Stop_Codon;
 
34
  // Sequences assumed to be stop codons
 
35
 
 
36
 
 
37
int  main
 
38
    (int argc, char * argv [])
 
39
 
 
40
  {
 
41
   FILE  * fp;
 
42
   string  Data, hdr;
 
43
   char  * Buffer, Line [MAX_LINE], Name [MAX_LINE];
 
44
   char  Codon [4] = "tag";
 
45
   int  Direction, Frame_Shift;
 
46
   long int  Buffer_Len, Gene_Len;
 
47
   long int  i, j, Begin, End, Len, Start, Stop;
 
48
   int  problem_ct = 0, ok_ct = 0;
 
49
   
 
50
   try
 
51
     {
 
52
      Parse_Command_Line (argc, argv);
 
53
 
 
54
      Set_Start_And_Stop_Codons ();
 
55
 
 
56
      fp = File_Open (Sequence_File_Name, "r");
 
57
 
 
58
      Buffer = (char *) Safe_malloc (INIT_SIZE);
 
59
      Buffer_Len = INIT_SIZE;
 
60
 
 
61
      Fasta_Read (fp, Data, hdr);
 
62
 
 
63
      fclose (fp);
 
64
 
 
65
      Len = Data . length ();
 
66
      Data . insert (Data . begin (), 'x');
 
67
        // Put a dummy character at the front of Data so subscripts
 
68
        // will start at 1
 
69
 
 
70
      fp = File_Open (Coord_File_Name, "r");
 
71
 
 
72
      while  (fgets (Line, MAX_LINE, fp) != NULL)
 
73
        {
 
74
         bool  problem = false;
 
75
 
 
76
         if  (sscanf (Line, "%s %ld %ld", Name, & Start, & End) != 3)
 
77
             {
 
78
              printf ("Bad line:  %s\n...Skipping\n", Line);
 
79
              continue;
 
80
             }
 
81
 
 
82
         if  (Start < End && End - Start <= Len / 2 || Start - End > Len / 2)
 
83
             {
 
84
              Direction = +1;
 
85
              Gene_Len = 1 + End - Start;
 
86
              if  (Gene_Len < 0)
 
87
                  Gene_Len += Len;
 
88
 
 
89
              if  (Buffer_Len < Gene_Len + 1)
 
90
                  Buffer = (char *) Safe_realloc (Buffer, 1 + Gene_Len);
 
91
              Buffer_Len = 1 + Gene_Len;
 
92
              for  (i = 0;  i < Gene_Len;  i ++)
 
93
                {
 
94
                 if  (Start + i <= Len)
 
95
                     j = Start + i;
 
96
                   else
 
97
                     j = Start + i - Len;
 
98
                 Buffer [i] = tolower (Data [j]);
 
99
                }
 
100
              Buffer [i] = '\0';
 
101
             }
 
102
           else
 
103
             {
 
104
              Direction = -1;
 
105
              Gene_Len = 1 + Start - End;
 
106
              if  (Gene_Len < 0)
 
107
                  Gene_Len += Len;
 
108
 
 
109
              if  (Buffer_Len < Gene_Len + 1)
 
110
                  Buffer = (char *) Safe_realloc (Buffer, 1 + Gene_Len);
 
111
              Buffer_Len = 1 + Gene_Len;
 
112
              for  (i = 0;  i < Gene_Len;  i ++)
 
113
                {
 
114
                 if  (Start - i >= 1)
 
115
                     j = Start - i;
 
116
                   else
 
117
                     j = Start - i + Len;
 
118
                 Buffer [i] = Complement (tolower (Data [j]));
 
119
                }
 
120
              Buffer [i] = '\0';
 
121
             }
 
122
 
 
123
         if  (Check_Previous_Stop)
 
124
             {
 
125
              if  (Direction == +1)
 
126
                  {
 
127
                   for  (i = 3;  i > 0;  i --)
 
128
                     if  (Start - i < 1)
 
129
                         Codon [i] = tolower (Data [Start - i + Len]);
 
130
                       else
 
131
                         Codon [i] = tolower (Data [Start - i]);
 
132
                  }
 
133
                else
 
134
                  {
 
135
                   for  (i = 3;  i > 0;  i --)
 
136
                     if  (Start + i > Len)
 
137
                         Codon [i] = Complement (tolower (Data [Start + i - Len]));
 
138
                       else
 
139
                         Codon [i] = Complement (tolower (Data [Start + i]));
 
140
                  }
 
141
              if  (! Is_Stop_Codon (Codon))
 
142
                  {
 
143
                   printf ("%-10s %8ld %8ld no stop before start\n",
 
144
                                  Name, Start, End);
 
145
                   problem = true;
 
146
                  }
 
147
             }
 
148
         if  (Check_Start_Codon && ! Is_Start_Codon (Buffer))
 
149
             {
 
150
              printf ("%-10s has bad start codon = %.3s\n", Name, Buffer);
 
151
              problem = true;
 
152
             }
 
153
         if  (! Is_Stop_Codon (Buffer + Gene_Len - 3))
 
154
             {
 
155
              printf ("%-10s has bad stop codon = %s\n", Name, Buffer + Gene_Len - 3);
 
156
              problem = true;
 
157
              for  (j = Gene_Len;  j < Len;  j += 3)
 
158
                {
 
159
                 for  (i = 0;  i < 3;  i ++)
 
160
                   if  (Direction == +1)
 
161
                       {
 
162
                        if  (Start + i + j > Len)
 
163
                            Codon [i] = tolower (Data [Start + i + j - Len]);
 
164
                          else
 
165
                            Codon [i] = tolower (Data [Start + i + j]);
 
166
                       }
 
167
                     else
 
168
                       {
 
169
                        if  (Start - i - j < 1)
 
170
                            Codon [i] = Complement (tolower (Data [Start - i - j + Len]));
 
171
                          else
 
172
                            Codon [i] = Complement (tolower (Data [Start - i - j]));
 
173
                       }
 
174
                 if  (Is_Stop_Codon (Codon))
 
175
                     break;
 
176
                }
 
177
              assert (j < Len);
 
178
              printf ("           next stop occurs at offset %ld"
 
179
                   "  Gene_Len = %ld  diff = %+ld\n",
 
180
                   j, Gene_Len, j - Gene_Len + 3);
 
181
             }
 
182
 
 
183
         Frame_Shift = (Gene_Len % 3);
 
184
         if  (Frame_Shift)
 
185
             {
 
186
              printf ("%-10s %8ld %8ld has %+d frame shift\n",
 
187
                              Name, Start, End, Frame_Shift);
 
188
              problem = true;
 
189
 
 
190
              for  (i = 0;  i < Gene_Len - 3;  i += 3)
 
191
                if  (Is_Stop_Codon (Buffer + i))
 
192
                    break;
 
193
              if  (i < Gene_Len - 3)
 
194
                  {
 
195
                   Stop = Start + Direction * (i - 1);
 
196
                   if  (Stop < 1)
 
197
                       Stop += Len;
 
198
                   else if  (Stop > Len)
 
199
                       Stop -= Len;
 
200
                   printf ("   Best prefix is %8ld %8ld   Len = %ld\n",
 
201
                                Start, Stop, i);
 
202
                  }
 
203
                else
 
204
                  {
 
205
                   printf ("   No stop found in start frame\n");
 
206
                   continue;
 
207
                  }
 
208
 
 
209
              for  (i = Gene_Len - 6;  i >= 0;  i -= 3)
 
210
                if  (Is_Stop_Codon (Buffer + i))
 
211
                    break;
 
212
              i += 3;
 
213
              Begin = Start + Direction * i;
 
214
              if  (Begin < 1)
 
215
                  Begin += Len;
 
216
              else if  (Stop > Len)
 
217
                  Begin -= Len;
 
218
              printf ("   Best suffix is %8ld %8ld   Len = %ld\n",
 
219
                           Begin, End, Gene_Len - i - 3);
 
220
 
 
221
             }
 
222
           else
 
223
             {
 
224
              for  (i = 0;  i < Gene_Len - 3;  i += 3)
 
225
                if  (Is_Stop_Codon (Buffer + i))
 
226
                    {
 
227
                     printf ("%-10s has stop codon %.3s at offset %ld"
 
228
                          "  Gene_Len = %ld  diff = %+ld\n",
 
229
                          Name, Buffer + i, i, Gene_Len, Gene_Len - 3 - i);
 
230
                     problem = true;
 
231
                    }
 
232
             }
 
233
         if  (problem)
 
234
             problem_ct ++;
 
235
           else
 
236
             ok_ct ++;
 
237
        }
 
238
 
 
239
      fprintf (stderr, "     OK orfs = %7d\n", ok_ct);
 
240
      fprintf (stderr, "Problem orfs = %7d\n", problem_ct);
 
241
     }
 
242
   catch (std :: exception & e)
 
243
     {
 
244
      cerr << "** Standard Exception **" << endl;
 
245
      cerr << e << endl;
 
246
      exit (EXIT_FAILURE);
 
247
     }
 
248
 
 
249
   return  0;
 
250
  }
 
251
 
 
252
 
 
253
 
 
254
static bool  Is_Start_Codon
 
255
    (const char * p)
 
256
 
 
257
//  Return true iff the first 3 characters of p match a
 
258
//  string in global  Start_Codon .
 
259
 
 
260
  {
 
261
   int  i;
 
262
 
 
263
   for  (i = 0;  i < Num_Start_Codons;  i ++)
 
264
     if  (strncmp (p, Start_Codon [i], 3) == 0)
 
265
         return  true;
 
266
 
 
267
   return  false;
 
268
  }
 
269
 
 
270
 
 
271
 
 
272
static bool  Is_Stop_Codon
 
273
    (const char * p)
 
274
 
 
275
//  Return true iff the first 3 characters of p match a
 
276
//  string in global  Stop_Codon .
 
277
 
 
278
  {
 
279
   int  i;
 
280
 
 
281
   for  (i = 0;  i < Num_Stop_Codons;  i ++)
 
282
     if  (strncmp (p, Stop_Codon [i], 3) == 0)
 
283
         return  true;
 
284
 
 
285
   return  false;
 
286
  }
 
287
 
 
288
 
 
289
 
 
290
static void  Parse_Command_Line
 
291
    (int argc, char * argv [])
 
292
 
 
293
//  Get options and parameters from command line with  argc
 
294
//  arguments in  argv [0 .. (argc - 1)] .
 
295
 
 
296
  {
 
297
   char  * p, * q;
 
298
   bool  errflg = false;
 
299
   int  ch;
 
300
 
 
301
   optarg = NULL;
 
302
 
 
303
   while  (! errflg && ((ch = getopt (argc, argv, "A:stZ:")) != EOF))
 
304
     switch  (ch)
 
305
       {
 
306
        case  'A' :
 
307
          Start_Codon . clear ();
 
308
          for  (p = strtok (optarg, ",");  p != NULL;  p = strtok (NULL, ","))
 
309
            {
 
310
             q = strdup (p);
 
311
             Make_Lower_Case (q);
 
312
             Start_Codon . push_back (q);
 
313
            }
 
314
          break;
 
315
 
 
316
        case  's' :
 
317
          Check_Start_Codon = FALSE;
 
318
          break;
 
319
 
 
320
        case  't' :
 
321
          Check_Previous_Stop = TRUE;
 
322
          break;
 
323
 
 
324
        case  'Z' :
 
325
          Stop_Codon . clear ();
 
326
          for  (p = strtok (optarg, ",");  p != NULL;  p = strtok (NULL, ","))
 
327
            {
 
328
             q = strdup (p);
 
329
             Make_Lower_Case (q);
 
330
             Stop_Codon . push_back (q);
 
331
            }
 
332
          break;
 
333
 
 
334
        case  '?' :
 
335
          fprintf (stderr, "Unrecognized option -%c\n", optopt);
 
336
 
 
337
        default :
 
338
          errflg = true;
 
339
       }
 
340
 
 
341
   if  (errflg)
 
342
       {
 
343
        Usage ();
 
344
        exit (EXIT_FAILURE);
 
345
       }
 
346
 
 
347
   if  (optind > argc - 2)
 
348
       {
 
349
        Usage ();
 
350
        exit (EXIT_FAILURE);
 
351
       }
 
352
 
 
353
   Sequence_File_Name = argv [optind ++];
 
354
   Coord_File_Name = argv [optind ++];
 
355
 
 
356
   return;
 
357
  }
 
358
 
 
359
 
 
360
 
 
361
static void  Set_Start_And_Stop_Codons
 
362
    (void)
 
363
 
 
364
//  Set globals  Start_Codon  and  Stop_Codon  to the sequences
 
365
//  that are allowed to be start and stop codons for genes.
 
366
 
 
367
  {
 
368
   int  i, n;
 
369
 
 
370
   if  (Start_Codon . size () == 0)
 
371
       {
 
372
        n = sizeof (DEFAULT_START_CODON) / sizeof (char *);
 
373
        for  (i = 0;  i < n;  i ++)
 
374
          Start_Codon . push_back (DEFAULT_START_CODON [i]);
 
375
       }
 
376
 
 
377
   if  (Stop_Codon . size () == 0)
 
378
       {
 
379
        n = sizeof (DEFAULT_STOP_CODON) / sizeof (char *);
 
380
        for  (i = 0;  i < n;  i ++)
 
381
          Stop_Codon . push_back (DEFAULT_STOP_CODON [i]);
 
382
       }
 
383
 
 
384
   Num_Start_Codons = Start_Codon . size ();
 
385
   Num_Stop_Codons = Stop_Codon . size ();
 
386
 
 
387
   return;
 
388
  }
 
389
 
 
390
 
 
391
 
 
392
static void  Usage
 
393
    (void)
 
394
 
 
395
//  Print to stderr description of options and command line for
 
396
//  this program.
 
397
 
 
398
  {
 
399
   fprintf (stderr,
 
400
       "USAGE:  anomaly [options] <sequence-file> <coord-file>\n"
 
401
       "\n"
 
402
       "Read DNA sequence in <sequence-file> and for each region specified\n"
 
403
       "by the coordinates in <coord-file>, check whether the region\n"
 
404
       "represents a normal gene, i.e., it begins with a start codon, ends\n"
 
405
       "with a stop codon, and has no frame shifts.\n"
 
406
       "Output goes to standard output.\n"
 
407
       "\n"
 
408
       "Options:\n"
 
409
       " -A <codon-list>\n"
 
410
       "    Use comma-separated list of codons as start codons\n"
 
411
       "    Sample format:  -A atg,gtg\n"
 
412
       "    Default start codons are atg,gtg,ttg\n"
 
413
       " -s\n"
 
414
       "    Omit the check that the first codon is a start codon.\n"
 
415
       " -t\n"
 
416
       "    Check whether the codon preceding the start coordinate position\n"
 
417
       "    is a stop codon.  This is useful if the coordinates represent\n"
 
418
       "    the entire region between stop codons.\n"
 
419
       " -Z <codon-list>\n"
 
420
       "    Use comma-separated list of codons as stop codons\n"
 
421
       "    Sample format:  -Z tag,tga,taa\n"
 
422
       "\n");
 
423
 
 
424
   return;
 
425
  }
 
426
 
 
427
 
 
428