Counting DNA nucleotides in C

Question

I have written code to solve the following Rosalind problem. This is my first time writing in C and I would like a review of my code, particularly in regard to correctness and performance.

#include <stdio.h>
#include <string.h>
int main() {

    /* Input DNA of length at most 1000 nucleotides */
    char dna[1000];
    FILE *input;
    input = fopen("/Users/James/Desktop/rosalind_dna.txt", "r");
    fgets(dna, 1000, input);
    fclose(input);

    /* Count frequency of symbols 'A', 'C', 'G' and 'T' in DNA */
    int i, a, c, g, t;
    for(i=0; i<strlen(dna); i++) {
        if (dna[i] == 'A') {
            a += 1;
        }
        else if (dna[i] == 'C') {
            c += 1;
        }
        else if (dna[i] == 'G') {
            g += 1;
        }
        else if (dna[i] == 'T') {
            t += 1;
        }
    }

    /* Output frequency of DNA nucleotides */
    printf("%i %i %i %i", a, c, g, t);
    return 0;
}

Updated Code

#include <stdio.h>
#include <string.h>
#define READ_BUFFER_SIZE 1000

int main() {

    /* Input DNA of length at most 1000 nucleotides */
    char dna[READ_BUFFER_SIZE];
    FILE *input;    
    input = fopen("/Users/James/Desktop/rosalind_dna.txt", "r");
    if (input == NULL) {
        perror("Error");
    }
    fgets(dna, READ_BUFFER_SIZE, input);
    fclose(input);

    /* Count frequency of symbols 'A', 'C', 'G' and 'T' in DNA */
    int aCount, cCount, gCount, tCount;
    char *current;
    for(current = dna; *current; ++current) {
        switch(*current){
            case 'A':
                ++aCount;
                break;
            case 'C':
                ++cCount;
                break;
            case 'G':
                ++gCount;
                break;
            case 'T':
                ++tCount;
                break;
            default:
                printf("Error: Invalid nucleotide\n");
                break;
        }
    }

    /* Output frequency of DNA nucleotides */
    printf("%i %i %i %i", aCount, cCount, gCount, tCount);
    return 0;
}

Answer 1

I am no C programmer but there are a few things that come to my eye:

Static size

You always read at most 1000 characters, what if the file size is much bigger? Instead you would loop trough the file, reading/processing it chunkwise.

Magic Number

Even if you use the static size you should not enter it at multiple locations. Instead use a constant for it to have to apply changes only in one place:

#define READ_BUFFER_SIZE 1000
//...
char dna[READ_BUFFER_SIZE];
//...
fgets(dna, READ_BUFFER_SIZE, input);

Loop termination

It might be optimized by the compiler but unoptimized your loop end check is n times (for n = strlen(dna)) which itself is takes \$O(n)\$, making your loop \$O(n^2)\$.

Looping

You index through the loop which is more costly than walking a pointer through it like:

char *current;
// loop until the '\0' terminator is encountered
for(current = dna; *current; ++current)

However, note as stated in https://stackoverflow.com/a/2305783 that modern compiler will optimize the indexed version to the same efficient code as the pointer version.

Switching

To find the counter for the current base you are using cascading ifs, which is suboptimal. Instead use a switch statement:

switch(*current) {
default: assert("Invalid base name!");
case 'A':
//...

Incrementing

This is most likely optimized by the compiler but still: You are incrementing so use the (prefix) increment operator ++a instead of a += 1;.

Naming

In this short example the names are nearly okay but in a longer program you would want to expand the variable names e.g.:

• dna -> dnaChunkBuffer (assuming you work in chunks)
• input -> inputFile
• i -> currentBaseIndex
• a, c, g, t -> base*Count

Code-Reuse

As I said I am no C programmer but I could imagine that there is a standard algorithm which helps you count the different characters in a string. If so, I would go for that.

Error Handling

It is just a small example but what happens if there is no file at the location you entered? Your program does not account for that. You should at least test if the fopen function succeeded and given an error otherwise.

Answer 2

The Rosalind problem takes in a string and counts the occurrences of A,C,G,T - reading the string from the file is not part of it so I shall skip over that.

In terms of correctness - does what is intended - it seems fine.

In terms of performance, as @Nobody mentioned, it is most likely that strlen() is being called multiple times, so doing it once outside the loop and using the saved value is a better plan (a compiler my optimize it out but I wouldn't rely on it)

If performance is/becomes an issue you might try removing the if/else/else setup by creating an array of values and indexing into them, but truthfully this is just messing around. For this situation the difference in performance would be negligible.

/* 20 slots, because we only need up as far as T */
int results[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
int strLen = strlen(inputString);
for (int index = 0; index < strLen; index++){
    results[inputString[index] - 'A']++;
}

aCount = results[0];
cCount = results['C'-'A'];
gCount = results['G'-'A'];
tCount = results['T'-'A']; 

Other Points

It is a matter of personal preference but once I move above 3 options I tend to go for a select() rather than the if/else/else

int strLen = strlen(inputString);

for (int index = 0; index < strLen; index++){
    switch (inputString[index]){
        case 'A':
            aCount++;
            break;
        case 'C':
            cCount++;
            break;
        case 'G':
            gCount++;
            break;
        case 'T':
            tCount++;
            break;
    }
}

as it can be a bit easier to follow but this will not impact correctness and AFAIK it does not have better performance than if / else / else.

Answer 3

While the revised version is better, I worry about the 1000 nucleotides limitation. If you ever want to operate on a larger data set, your memory allocation will fail.

It would be better to not buffer it in ram, as the file system already has numerous buffers.

I would change the loop to a:

while (c=fgetc(input)) 

then use an array of 4 for the sums:

unsigned int nucleotides[4];
int n; /* index into nucleotides, per the switch below: 0=a, 1=c, 2=g, 3=t*/

then index into the nucleotides and increment as:

switch (c)
{
    case 'A': n=0; break;
    case 'C': n=1; break;
    case 'G': n=2; break;
    case 'T': n=3; break;
}
nucleotides[n]++;

Finally after the loop:

printf("%i %i %i %i", nucleotides[0], nucleotides[1], nucleotides[2], nucleotides[3]);
return 0;

This will result in a program that is very light on ram, and doesn't spend a lot of time in strlen(); You can successfully run this version as long as there aren't more than 4 billion samples of a nucleotide. If there are more, then change the type of nucleotides to long long.

To make it slightly more readable you can use an enum or #define to map the indexes 0..3 to A..T respectively.