Data processing task for bioinformatics

Question

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); // Process reads 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; } 

Answer 1

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?

Answer 2

General comments:

• 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(...); ... }