Take the 2-minute tour ×
Code Review Stack Exchange is a question and answer site for peer programmer code reviews. It's 100% free, no registration required.

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;
}
share|improve this question

3 Answers 3

up vote 10 down vote accepted

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 http://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.

share|improve this answer
    
Good points ! However, I have serious doubts about the "You index through the loop which is much more costly than walking a pointer through it like" part. –  Josay Jun 13 at 10:19
    
@Josay: Can you expand on this? In my experience indexing is slower (when unoptimized) because it uses the formula baseAddress + sizeof(element) * index (and increments the index) whereas incrementing a pointer does not need this calculation. So indexing has one more addition and one more multiplication. –  Nobody Jun 13 at 10:25
    
@MatteoItalia: You are right. I discarded this when writing the code but I don't remember why. –  Nobody Jun 13 at 10:38
    
The "looping" section is not completely correct, indexing an array is substantially as cheap as working with pointers. The real problem is that he is potentially calling strlen at each iteration, which is O(n) in the length of the string, thus making the whole thing O(n^2). Also in your example strlen isn't needed (you walk the array twice for no reason), just do for(current = dna; *current; ++current). –  Matteo Italia Jun 13 at 10:39
    
@MatteoItalia: I am not sure you understood me correctly. On the performance of indexing vs pointers: stackoverflow.com/a/2305783 However, I was talking about unoptimized code. –  Nobody Jun 13 at 10:47

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.

share|improve this answer

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.

share|improve this answer
1  
you should indent your case statements inside the switch, other wise a well laid out answer. –  Malachi Jun 13 at 15:44

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.