Counting nucleobases in a nucleotide

Question

This question is part of a series solving the Rosalind challenges. For the previous question in this series, see Calculating protein mass ruby. The repository with all my up-to-date solutions so far can be found here.

---

I started the Rosalind challenges roughly a year ago in Ruby. Now I got curious whether I could do the same challenges in C.

Problem: DNA

A string is simply an ordered collection of symbols selected from some alphabet and formed into a word; the length of a string is the number of symbols that it contains. An example of a length 21 DNA string (whose alphabet contains the symbols 'A', 'C', 'G', and 'T') is "ATGCTTCAGAAAGGTCTTACG."

Given:

A DNA string s of length at most 1000 nt.

Return:

Four integers (separated by spaces) counting the respective number of times that the symbols 'A', 'C', 'G', and 'T' occur in s.

Sample Dataset:

AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC

Sample Output:

20 12 17 21

Actual Dataset:

CTCCTCAGATCTCAAACGGCTCTATATTACTAGATAGGAGACACGCCCATACCAGCGACGCGGGGTCACTCATTTTCCCAAGAATCCATGAGTGCGAAGCGCACGTCCATGTGACACAAAATTACTAGAGAGTTTTCAAGTCTGATTACCCGTAGTAAACGACCTTGTGCCGGGTCACTAGTGCAATGAAGAATATGTCAACTATTACTCCCGTGGGATCTATAAAACCAGAAGATCCATTGCACTTGTAGTCGCTGTATAGTCTCTCGTCGTCACCTAGCCGATATGACCGTGCGCGAGTTATCCGGAACCTATAAGTGTTTGCTCTCAACAGTGTCTCAACACATGGAGTCGGTAACCTACTACGAAGCCTGCACCAAGATCGATCAGGGAGAATACCCCCTGACGGTCAACGCCGAAGATCAAAGAGAATGATTCGGCCTAGGGCGATTGGCTATTATCCCGGTCTAACCGCCAGGATACTTCAGTAGATCCCGCTCGACATCTGCCCCCCACAAAGTTATTCAGTTTCGGTGATAATTTCGCTTGAACTCCTATCTATTTAAAAGTTTTCCTATACGATGACTAGTCCCTTGCGAACGATCTTTGCCAGGATGCACGACGGCGAGACAATATTACAATACCGAGTGGAGTGATTGGTATCTACACATACGAAATCTCAATGAGAATGGAAGGTCACACTCGTAACAAACTCCTAAGCGGCGGAGAGCGGAAAGGTATAGTCGAGTCGAAGCCTTTATATCGTGTGGCCAGCAGCTAACACAGAGAAATATGGCGGGAATCATC

Actual output:

231 201 181 194

DNA.c:

#include "DNA.h"

int main()
{
    size_t MAX_LENGTH = 1000;
    char *user_input;

    user_input = (char *) malloc (MAX_LENGTH + 1);
    getline (&user_input, &MAX_LENGTH, stdin);
    output(countACGT(user_input));
}

DNA.h:

#include <ctype.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

typedef struct {
    uint8_t A;
    uint8_t C;
    uint8_t G;
    uint8_t T;
} uint8_t_container;

uint8_t_container countACGT(char *nucleotide)
{
    uint16_t len = strlen(nucleotide);
    char nucleobase;
    uint8_t_container *count = calloc(4, sizeof(uint8_t_container));

    for(uint16_t i = 0; i < len; i++)
    {
        nucleobase = nucleotide[i];
        if isalnum(nucleobase)
        {
            if (nucleobase == 'A')
            {
                count->A += 1;
            }
            else if (nucleobase == 'C')
            {
                count->C += 1;
            }
            else if (nucleobase == 'G')
            {
                count->G += 1;
            }
            else if (nucleobase == 'T')
            {
                count->T += 1;
            }
        }
    }
    return *count;
}

void output(uint8_t_container result)
{
    printf("%i %i %i %i\n", result.A, result.C, result.G, result.T);
}

Explanation:

I prepare the user_input which will contain the nucleotide. The nucleotide (basically an array of nucleobases) will be iterated over. If the input is alphanumeric it will get checked for it's value and the appropriate member value of count (a container of uint8_t) will be increased by one. While the challenge only provides sanitized input data, this should further decrease the likelihood of useless data wreaking havoc. At the end, the container will be retrieved and read. The \n in the printf statement is not required but keeps my terminal sane.

My typing and naming can probably be improved and I'm sure I missed some obvious optimizations. An earlier version used an array inside uint8_t_container, but I decided this probably looked better. Feel free to prove me wrong.

The last time I've written a decent amount of C is a while ago, so I'm especially looking for things that may bite me in the behind when done in production. Alternative solutions are only useful as part of a review or won't be of use. I'm perfectly aware this isn't the most straightforward method of solving the challenge, but I was going for modularity.

Compiled using:

gcc -Wall -Wextra -pedantic DNA.c -o DNA.out

As far as I know that should give me pretty much all warnings available.

Answer 1

There are a few problems with memory management in your original program. Rather than addressing all of them I would suggest you eliminate malloc and calloc altogether in your code. There is no need for them here since you do not need to extend the lifetime of any objects and you are not dealing with unknown sizes or complex data types.

Suggestions

• Naming

  uint8_t_container is not descriptive at all. Change it to nucleobase_frequency or nucleotide_frequency.

• With small problem spaces prefer to use more memory than more computation time

  char can only represent 256 values on most implementations (likely any you will encounter). You can eliminate the if statement that checks whether a character is a valid nucleotide is 'A', 'C', 'G', or 'T' altogether. Instead count the frequency of every character and reduce to the characters you want later.

• Specify pointers to output structs in your parameter list

  The paradigm in C for anything beyond basic data types is to pass pointers to output types instead of returning the type from the function.

  Instead of this

  some_type func (int a, char z);
  

  Prefer this

  void func (some_type* out, int a, char z);
  

  To call this function use

   some_type s;
   func (&s, 23, 'c');
  

• strlen is not required

  As was pointed out in the comments by @Eric Lagergren there is no need to use strlen since we are iterating through the string anyway. Instead you can use the NULL byte as the stopping condition during iteration.

• Use const

  Use const whenever you pass in a pointer to an object or an array that will not change.

• Unnecessary Copy

  The paradigm in C is for functions to accept pointers to structs rather than the structs themselves

  Change

  void output(uint8_t_container result)
  

  to

  void output(const uint8_t_container* result)
  

• Put only function prototypes and not implementation in the .h

Updated Code

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

struct nucleobase_frequency
{
    size_t A;
    size_t C;
    size_t G;
    size_t T;
};

void count_nucleobases (struct nucleobase_frequency* frequency, const char *sequence)
{
    if (frequency == NULL || sequence == NULL)
    {
        return;
    }

    // In practically every environment you use char will
    // be 8-bit which represents 256 possible values.
    // Create and zero-out an array of size 256
    size_t freq[256] = { 0 };

    // Iterate until we reach the NULL byte
    for (size_t i = 0; sequence[i]; i++)
    {
        // Cast the char to unsigned char so we
        // do not access out-of-bounds indices
        unsigned char c = sequence[i];
        freq[c]++;
    }

    // Choose just the frequencies we are interested in
    frequency->A = freq['A'];
    frequency->C = freq['C'];
    frequency->G = freq['G'];
    frequency->T = freq['T'];
}

void print_frequency (const struct nucleobase_frequency* frequency)
{
    if (frequency == NULL)
    {
        return;
    }

    printf("A: %zu\n", frequency->A);
    printf("C: %zu\n", frequency->C);
    printf("G: %zu\n", frequency->G);
    printf("T: %zu\n", frequency->T);
}

Answer 2

Advice 1

typedef struct {
    uint8_t A;
    uint8_t C;
    uint8_t G;
    uint8_t T;
} uint8_t_container;

In real world setting, uint8_t values are not sufficient for representing (absolute) frequencies, consider using at least uint32_t.

Advice 2

if isalnum(nucleobase)

You forgot the parentheses. Would be more funky as

if (isalnum(nucleobase))

Advice 3

uint8_t_container *count = calloc(4, sizeof(uint8_t_container));

You ask 4 times more memory than you need. Consider using

uint8_t_container *count = calloc(1, sizeof(uint8_t_container)); 
                                  ^

Advice 4

return *count;

Above, count is a pointer to a structure, I am not 100% sure, yet it appears to me that since the return value of countACGT is not a pointer to that structure, you convert *count to a value copy, and leave the actual structure being referenced without deallocating: a memory leak.

Advice 5

if (nucleobase == 'A')
{
    count->A += 1;
}
else if (nucleobase == 'C')
{
    count->C += 1;
}
else if (nucleobase == 'G')
{
    count->G += 1;
}
else if (nucleobase == 'T')
{
    count->T += 1;
}

This is asking for a switch.

Advice 6

Finally, more idiomatic C would be having a char pointer sliding through the nucleotides and halting at zero terminator. That way, no need for strlen.

Summa summarum

All in all, I had this in mind:

typedef struct {
    uint32_t a_count;
    uint32_t c_count;
    uint32_t g_count;
    uint32_t t_count;
} uint32_t_container;

uint32_t_container count_nucleotides(char *nucleotides)
{
    uint32_t_container result;

    result.a_count = 0;
    result.c_count = 0;
    result.g_count = 0;
    result.t_count = 0;

    for (char *c = nucleotides;; c++)
    {
        switch (*c)
        {
            case 'A':
                result.a_count++;
                break;

            case 'C':
                result.c_count++;
                break;

            case 'G':
                result.g_count++;
                break;

            case 'T':
                result.t_count++;
                break;

            case 0:
                return result;
        }
    }
}

Hope that helps.

Answer 3

A .h file containing function definitions achieves nothing but trouble. For starters you cannot #include such file in more than one source (the only excuse to place a function definition in an include file is static inline qualifier). A correct modularization is a header containing declarations, and another .c file providing definitions:

// DNA.h
#ifndef DNA_H
#define DNA_H
typedef struct {
    uint8_t A;
    uint8_t C;
    uint8_t G;
    uint8_t T;
} uint8_t_container;

uint8_t_container countACGT(char *nucleotide);
void output(uint8_t_container result);
#endif

// DNA.c
#include "DNA.h"
uint8_t_container countACGT(char *nucleotide)
{
    ....
}
void output(uint8_t_container result)
{
    ....
}

// main.c
#include "DNA.h"
int main()
{
    ....
}

and build the application with

gcc -o nucleobase_counter main.c dna.c

In real life you would also have a Makefile to automate the build.

Answer 4

You've already had some good suggestions covering most of the main ideas, however there's always something else to think about.

Consider 'const correctness'

If you're not going to be modifying a parameter (such as your nucleotide), then consider declaring it as const. This makes it clear to the caller that they can rely on the information not changing. It also helps to stop you accidentally overwriting values, which brings me on to:

Consider following a literals left approach

In C, having your literals on the right can be error prone (a small typo and == becomes =). Unlike languages like C#, both this:

if (nucleobase == 'A')

And this:

if (nucleobase = 'A')

Will compile quite happily (although you may get a warning depending on your compiler / options). The difference being that in the second statement nucleobase is actually updated every time to 'A', so the statement always triggers. Following a literals left approach, you can't make this mistake. Whilst this:

if ('A' == nucleobase)

Compiles, this throws an error:

if ('A' = nucleobase)

memory leak in main

In addition to the memory leak in your count ACGT method, you're also not cleaning up memory in main (you malloc for user_input and never free it). Whilst not strictly necessary, I tend to free all malloced buffers before exiting.

Check return values

You're making calls to library functions that return status values as part of the function call. Whilst it's unlikely that you're going to get failures from either malloc or getline for this application getting into the habit of checking the return values to validate them will mean that you don't forget when it gets to your production code.

Naming

In my experience, all caps variable names are usually reserved for constants / macros. Your MAX_LENGTH variable defined in main isn't declared as const however because of it's name may be interpreted as being a constant on first inspection (I know that's how I interpreted it). As it stands, the value of MAX_LENGTH may be updated by your getline call. Again, in your current application this is unlikely to be an issue however in production code it might do (particularly if it was declared at global scope).

struct initializer

If you go with the approach suggested by @coderodde, rather than calloc, I'd suggest using the struct initializer {0} rather than initialising each member separately.

Modifying their code for the above would turn it into (note I've also used the nucleotides parameter, rather than declaring a new local variable because to me it looks cleaner):

uint32_t_container count_nucleotides(const char *nucleotides)
{
    uint32_t_container result={0};

    for (;; nucleotides++)
    {
        switch (*nucleotides)
        {
            case 'A':
                result.a_count++;
                break;

            case 'C':
                result.c_count++;
                break;

            case 'G':
                result.g_count++;
                break;

            case 'T':
                result.t_count++;
                break;

            case 0:
                return result;
        }
    }
}