# Data processing task for bioinformatics

I threw together this C program today to handle a bioinformatics data processing task. The program seems to work correctly, but I wanted to know if anyone has suggestions regarding how the input data are parsed and how I've used control structures in the main processing loop.

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define BUFFERSIZE 1024

FILE *fileopen(const char *filename, const char *mode)
{
FILE *fh = fopen(filename, mode);
if(fh == NULL)
{
fprintf(stderr, "error: unable to open file '%s'\n", filename);
exit(1);
}
return fh;
}

void *memalloc(size_t size)
{
void *memory = malloc(size);
if(memory == NULL)
{
fprintf(stderr, "error: unable to allocate memory\n");
exit(1);
}
return memory;
}

int main(int argc, const char **argv)
{
// Parse command line arguments
char seqbuffer1[BUFFERSIZE];
char seqbuffer2[BUFFERSIZE];
char qualbuffer1[BUFFERSIZE];
char qualbuffer2[BUFFERSIZE];
char *seqbufferout;
FILE *seqinfile;
FILE *seqoutfile;
FILE *qualinfile;
FILE *qualoutfile;
int seqlength;
if(argc != 6)
{
fprintf(stderr, "error: 5 arguments required, %d provided\n", argc - 1);
exit(1);
}
seqinfile   = fileopen(argv[1], "r");
qualinfile  = fileopen(argv[2], "r");
seqoutfile  = fileopen(argv[3], "w");
qualoutfile = fileopen(argv[4], "w");
seqlength   = atoi(argv[5]);
assert(seqlength > 0);

seqbufferout = (char *)memalloc( sizeof(char) * (seqlength + 1) );
seqbufferout[seqlength] = '\0';
while(fgets(seqbuffer1, BUFFERSIZE, seqinfile))
{
if(strlen(seqbuffer1) > 0)
{
int ambiguous;
int i;

assert(strncmp(seqbuffer1, ">", 1) == 0);
if(fgets(seqbuffer2, BUFFERSIZE, seqinfile) == NULL)
{
fprintf(stderr, "error: file ends with a fasta header\n");
exit(1);
}

if( fgets(qualbuffer1, BUFFERSIZE, qualinfile) == NULL ||
fgets(qualbuffer2, BUFFERSIZE, qualinfile) == NULL )
{
fprintf(stderr, "error: no quality entry corresponding to sequence '%s'\n", seqbuffer1);
exit(1);
}

ambiguous = 0;
for(i = 0; i < seqlength; i++)
{
if(seqbuffer2[i] == '.')
{
ambiguous = 1;
break;
}
seqbufferout[i] = seqbuffer2[i];
}

if(!ambiguous)
{
fprintf(seqoutfile,  "%s%s\n", seqbuffer1, seqbufferout);
fputs(qualbuffer1, qualoutfile);

i = 0;
char *tok = strtok(qualbuffer2, " ");
while(tok != NULL && i < seqlength - 1)
{
if(i > 0)
fputs(" ", qualoutfile);
fputs(tok, qualoutfile);
tok = strtok(NULL, " ");
i++;
}
fputs("\n", qualoutfile);
}
}
}

fclose(seqinfile);
fclose(qualinfile);
fclose(seqoutfile);
fclose(qualoutfile);
free(seqbufferout);

return 0;
}


Instead of the big indent block you make by doing if(strlen(seqbuffer1) > 0), do: if (strlen(seqbuffer1) <1) continue;, so the entire rest of the loop can be outdented a level, making it more readable.

Similarly, you can turn the if(!ambiguous) into if (ambiguous) continue; and outdent that block as well.

Depending on what version of C you're using, ambiguous could be a boolean instead of an int.

You could be a little more consistent with your NULL comparisons to the result of fgets() - the one in the while() condition is implicitly != NULL but the ones in the if() statements are explicitly == NULL. IMO you should be either explicit everywhere or implicit everywhere - either choice is fine, but mixing styles makes your code a bit more unclear.

assert(strncmp(seqbuffer1, ">", 1) == 0); is the same as the more-easily-understood assert(seqbuffer1[0] == '>'), right?

Also, you've got at least one potential bug: What if BUFFERSIZE is smaller than input line length and then seqlength > BUFFERSIZE?

• While I agree that large statements are hard to read, I think continue is even worse, it leads to spaghetti and is generally considered bad practice (it is banned by MISRA-C for example). Instead, use functions. if(strlen(seqbuffer1) > 0)) parse_string(seqbuffer); This is how you make programs clearer. – Lundin May 29 '12 at 6:24
• continue can really only exacerbate spaghetti code, not cause it: it's only unclear if it's in the middle of a large block... but I'd argue that's a problem with having a large block, not with continue. – pjz May 29 '12 at 18:47
• Here is an example of spaghetti caused solely by the continue keyword: for(i=0; i<condition; i++) { if(something) { continue; } a[count] = b[i]; count++; } The code is hard to read, because it is not obvious to the reader when i and count are increased respectively. – Lundin May 31 '12 at 14:11
• Apart from that, the general definition of spaghetti no matter programming language, is usually "code that jumps both upwards and downwards". Code that only jumps downwards is usually considered acceptable, even if it is using goto. C language continue does not fit in anywhere, except in code jumping upwards. – Lundin May 31 '12 at 14:16
• I think that can be argued both ways; virtualdub.org/blog/pivot/entry.php?id=238 argues opposite to you. – pjz Jun 1 '12 at 19:13

• use perror instead of fprintf on error
• use EXIT_FAILURE as an exit code on failure

main()

• seqbufferout and hence memalloc() are redundant
• but as it is, the call to memalloc should omit sizeof(char), which is 1 by default and should omit the cast.
• fgets calls maybe use sizeof buf instead of BUFFERSIZE
• is fgets(seqbuffer1...) != NULL and strlen(seqbuffer1) == 0 a possible combination? I don't think so - in which case the strlen is redundant
• assert on seqbuffer1[0] == '>', not on a strcmp
• for-loops like this that do a well-defined task should be extracted into a function
• but this for loops seems just to look for a '.'; strchr would be better
• copying seqbuffer2 to seqbufferout seems redundant as you just print the contents of seqbufferout
• the while loop should be extracted into a separate function; it looks as if the loop just strips leading and trailing spaces and multiple spaces from the qualbuffer2 string - is that right?

• so (unless I misunderstood the purpose of the code) the stuff from and including ambiguous = 0 becomes:

if (!strchr(seqbuffer2, '.'))
{
fprintf(seqoutfile,  "%s%s\n", seqbuffer1, seqbuffer2);
fputs(qualbuffer1, qualoutfile);
fputs(strip_spaces(qualbuffer2), qualoutfile);
fputs("\n", qualoutfile);
}

• no checks for error in writing to the output files

• no check for error in seqinfile on exit from the main loop. You need to append the following before closing the files

if (ferror(seqinfile))
{
perror(...);
...
}