# Counting nucleobases in a nucleotide

This question is part of a series solving the Rosalind challenges. For the previous question in this series, see Calculating protein mass . 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.

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.

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

• I'm curious why build a 256 byte array and why use strlen? A switch might just get compiled into a jump table and if you increment the array with pointers you get the same effect as strlen just without having to parse the array twice. – Eric Lagergren Dec 10 '16 at 20:14
• @EricLagergren You are absolutely right about strlen being both unnecessary and wasteful. I have edited the code to remove it. – twohundredping Dec 11 '16 at 6:35
• @EricLagergren size_t is 8 bytes on my machine so the array is actually 2KB. Regardless the timing on my machine (Intel Core i7-4810MQ CPU @ 2.80GHz, gcc 5.4.0 -std=c11 -O3) for a sequence of size 1000000000 shows the frequency table approach takes 590ms and the switch statement without a default case takes 5830ms. – twohundredping Dec 11 '16 at 7:11
• Good catch on the unnecessary copy, that's a quite a major flaw. I'm not sure about the size_t freq[256] = { 0 }; part though, I don't think size_t is supposed to be used for such a trick. It also adds a magic number where I was just glad the answers managed to reduce the amount of such numbers. – Mast Dec 11 '16 at 14:03
• @Mast You can use CHAR_BIT to determine the number of bits in a char. From there you define the array to have size 1 << CHAR_BIT. I just didn't want to muddy up the code with such pedantry. Anyway, I am not sure what trick you think I am using. I am just creating an array indexed by all possible values each character can represent. Additionally, nothing is up my sleeves regarding size_t here. It is the appropriate type to use to allow the problem to scale up past the artificial length-1000 limit. – twohundredping Dec 11 '16 at 16:44

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.

if isalnum(nucleobase)


You forgot the parentheses. Would be more funky as

if (isalnum(nucleobase))


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));
^


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.

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.

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.

• You got some interesting points there, I'll definitely take them into account. The calloc asking for 4x as much memory was a stupid mistake, I was actually trying to minimize my memory usage. I simply forgot C even had a switch :-) – Mast Dec 10 '16 at 15:23
• @Mast Well, after all we are all here to learn. :^) – coderodde Dec 10 '16 at 15:24

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.

• I considered using a Makefile overkill for this scenario, but you're right. Regardless of how small a program is, in the real world you'd give it a Makefile anyway. I considered splitting it up into 3 files, except I couldn't get my naming straight since I wanted to keep the name of the final file DNA.c. Is that a stupid idea and inappropriate as a name? Is there an unwritten rule that the final name always is main.c? – Mast Dec 10 '16 at 18:09
• @Mast There is no main.c rule; I was wrong to suggest such name. The name of a .c file containing a main() typically coincides with the name of the app. I don't think the app should be called DNA, for it doesn't convey its goal - on the other hand a file named DNA.c clearly provide handy functions for DNA analysis. – vnp Dec 10 '16 at 19:00
• I went for DNA.c as main file because the code name for the problem is DNA. Which reminds me I did link to the problem set but not to the actual problem. I'll fix that. I think, based on your answer and comment, that I can think of a better naming scheme for the next one. Thanks :-) – Mast Dec 10 '16 at 19:06

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;
}
}
}