Get the length of each line in file with C and write in output file

Question

I am a biology student and I am trying to learn perl, python and C and also use the scripts in my work. So, I have a file as follows:

>sequence1
ATCGATCGATCG
>sequence2
AAAATTTT
>sequence3
CCCCGGGG  

The output should look like this, that is the name of each sequence and the count of characters in each line and printing the total number of sequences in the end of the file.

sequence1 12
sequence2 8
sequence3 8
Total number of sequences = 3

I could make the perl and python scripts work, this is the python script as an example:

#!/usr/bin/python

import sys

my_file = open(sys.argv[1]) #open the file
my_output = open(sys.argv[2], "w") #open output file

total_sequence_counts = 0

for line in my_file: 
    if line.startswith(">"):
        sequence_name = line.rstrip('\n').replace(">","") 
        total_sequence_counts += 1 
        continue    
    dna_length = len(line.rstrip('\n')) 
    my_output.write(sequence_name + " " + str(dna_length) + '\n')
my_output.write("Total number of sequences = " + str(total_sequence_counts) + '\n')

Now, I want to write the same script in C, this is what I have achieved so far:

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

int main(int argc, char *argv[])
{
input = FILE *fopen(const char *filename, "r");
output = FILE *fopen(const char *filename, "w"); 

double total_sequence_counts = 0;
char sequence_name[];

char line [4095]; // set a temporary line length
char buffer = (char *) malloc (sizeof(line) +1); // allocate some memory

while (fgets(line, sizeof(line), filename) != NULL) { // read until new line character is not found in line

    buffer = realloc(*buffer, strlen(line) + strlen(buffer) + 1); // realloc buffer to adjust buffer size
    if (buffer == NULL) { // print error message if memory allocation fails
        printf("\n Memory error");
        return 0;
    }
    if (line[0] == ">") {
        sequence_name = strcpy(sequence_name, &line[1]); 
        total_sequence_counts += 1
        }
        else {
            double length = strlen(line);
            fprintf(output, "%s \t %ld", sequence_name, length);
        }
    fprintf(output, "%s \t %ld", "Total number of sequences = ", total_sequence_counts);
}
    int fclose(FILE *input); // when you are done working with a file, you should close it using this function. 
    return 0;
    int fclose(FILE *output);
    return 0;
}

But this code, of course is full of mistakes, my problem is that despite studying a lot, I still can't properly understand and use the memory allocation and pointers so I know I especially have mistakes in that part. It would be great if you could comment on my code and see how it can turn into a script that actually work. By the way, in my actual data, the length of each line is not defined so I need to use malloc and realloc for that purpose.

Answer 1

For a simple program like this, where you look at short lines one at a time, you shouldn't worry about dynamic memory allocation. It is probably good enough to use local buffers of a reasonable size.

Another thing is that C isn't particularly suited for quick-and-dirty string processing. For example, there isn't a strstrip function in the standard library. You usually end up implementing such behaviour yourself.

An example implementation looks like this:

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



#define MAXLEN 80       /* Maximum line length, including null terminator */

int main(int argc, char *argv[])
{
    FILE *in;
    FILE *out;

    char line[MAXLEN];          /* Current line buffer */
    char ref[MAXLEN] = "";      /* Sequence reference buffer */
    int nseq = 0;               /* Sequence counter */

    if (argc != 3) {
        fprintf(stderr, "Usage: %s infile outfile\n", argv[0]);
        exit(1);
    }

    in = fopen(argv[1], "r");
    if (in == NULL) {
        fprintf(stderr, "Couldn't open %s.\n", argv[1]);
        exit(1);        
    }

    out = fopen(argv[2], "w");
    if (in == NULL) {
        fprintf(stderr, "Couldn't open %s for writing.\n", argv[2]);
        exit(1);        
    }

    while (fgets(line, sizeof(line), in)) {
        int len = strlen(line);

        /* Strip whitespace from end */
        while (len > 0 && isspace(line[len - 1])) len--;
        line[len] = '\0';

        if (line[0] == '>') {
            /* First char is '>': copy from second char in line */
            strcpy(ref, line + 1);
        } else {
            /* Other lines are sequences */
            fprintf(out, "%s: %d\n", ref, len);
            nseq++;
        }
    }

    fprintf(out, "Total number of sequences. %d\n", nseq);

    fclose(in);
    fclose(out);

    return 0;
}

A lot of code is about enforcing arguments and opening and closing files. (You could cut out a lot of code if you used stdin and stdout with file redirections.)

The core is the big while loop. Things to note:

• fgets returns NULL on error or when the end of file is reached.
• The first lines determine the length of the line and then remove white-space from the end.
• It is not enough to decrement length, at the end the stripped string must be terminated with the null character '\0'
• When you check the first character in the line, you should check against a char, not a string. In C, single and double quotes are not interchangeable. ">" is a string literal of two characters, '>' and the terminating '\0'.
• When dealing with countable entities like chars in a string, use integer types, not floating-point numbers. (I've used (signed) int here, but because there can't be a negative number of chars in a line, it might have been better to have used an unsigned type.)
• The notation line + 1 is equivalent to &line[1].
• The code I've shown doesn't check that there is always one reference per sequence. I'll leave this as exercide to the reader.

For a beginner, this can be quite a lot to keep track of. For small text-processing tasks like yours, Python and Perl are definitely better suited.

Edit: The solution above won't work for long sequences; it is restricted to MAXLEN characters. But you don't need dynamic allocation if you only need the length, not the contents of the sequences.

Here's an updated version that doesn't read lines, but read characters instead. In '>' context, it stored the reference. Otherwise it just keeps a count:

#include <stdlib.h>
#include <stdio.h>
#include <ctype.h>      /* for isspace() */



#define MAXLEN 80       /* Maximum line length, including null terminator */

int main(int argc, char *argv[])
{
    FILE *in;
    FILE *out;

    int nseq = 0;               /* Sequence counter */
    char ref[MAXLEN];           /* Reference name */

    in = fopen(argv[1], "r");
    out = fopen(argv[2], "w");

    /* Snip: Argument and file checking as above */

    while (1) {
        int c = getc(in);

        if (c == EOF) break;

        if (c == '>') {
            int n = 0;

            c = fgetc(in);
            while (c != EOF && c != '\n') {
                if (n < sizeof(ref) - 1) ref[n++] = c;
                c = fgetc(in);
            }
            ref[n] = '\0';
        } else {
            int len = 0;
            int n = 0;

            while (c != EOF && c != '\n') {
                n++;
                if (!isspace(c)) len = n;
                c = fgetc(in);
            }

            fprintf(out, "%s: %d\n", ref, len);
            nseq++;
        }
    }

    fprintf(out, "Total number of sequences. %d\n", nseq);

    fclose(in);
    fclose(out);

    return 0;
}

Notes:

• fgetc reads a single byte from a file and returns this byte or EOF when the file has ended. In this implementation, that's the only reading function used.
• Storing a reference string is implemented via fgetc here too. You could probably use fgets after skipping the initial angle bracket, too.
• The counting just reads bytes without storing them. n is the total count, len is the count up to the last non-space. (Your lines probably consist only of ACGT without any trailing space, so you could skip the test for space and use n instead of len.)

Answer 2

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

int main(int argc, char *argv[]){
    FILE *my_file = fopen(argv[1], "r");
    FILE *my_output = fopen(argv[2], "w");
    int total_sequence_coutns = 0;
    char *sequence_name;
    int dna_length;
    char *line = NULL;
    size_t size = 0;

    while(-1 != getline(&line, &size, my_file)){
        if(line[0] == '>'){
            sequence_name = strdup(strtok(line, ">\n"));
            total_sequence_coutns +=1;
            continue;
        }
        dna_length = strlen(strtok(line, "\n"));
        fprintf(my_output, "%s %d\n", sequence_name, dna_length);
        free(sequence_name);
    }
    fprintf(my_output, "Total number of sequences = %d\n", total_sequence_coutns);
    fclose(my_file);
    fclose(my_output);
    free(line);
    return (0);
}