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