C program to find complementary of DNA sequence

Question

[Arindam Ghosh]

I tried to write this C code to find complementary of DNA sequence. Please verify and help to improve.

#include<stdio.h>

main()
{
    int n;

    printf("\n\nEnter length of sequence:");
    scanf("%d",&n);

    char seq[n], com[n];
    int i;

    printf("\nENTER THE SEQUENCE:");
    scanf("%s",&seq);
    for(i;i<n;i++)
    {
        if(seq[i]=='A')
        com[i]='T';
        else if(seq[i]=='T')
        com[i]='A';
        else if(seq[i]=='G')
        com[i]='C';
        else if(seq[i]=='C')
        com[i]='G';
        else if(seq[i]==' ')
        com[i]='_';
        else if(seq[i]!='T' && seq[i]!='A' && seq[i]!='G' && seq[i]!='C')
        com[i]='*';
    }
    printf("\n\n%s\n\n",com);
    printf("\n\n* is non-DNA bases\n\n");
}

[5heikki]

Inputting DNA sequence by hand.. very convenient.

[Arindam Ghosh]

we can copy paste it atleast................

[Anima Mundi]

Just for fun, I made a Python version of the script (outputs reverse complement):

s = raw_input('Enter sequence:\n')

sup = s.upper()

c = ''
for char in sup:
    if char == 'A':
        c += 'T'
    elif char == 'T':
        c += 'A'
    elif char == 'C':
        c += 'G'
    elif char == 'G':
        c += 'C'
    else:
        c += char

print c[::-1]

[John]

A 1-liner might look like:

from string import maketrans
print raw_input('DNA: ').translate(maketrans('AaCcGgTt', 'TtGgCcAa'))[::-1]

Answer 1

[Pierre Lindenbaum]

declare 'main(int argc,char** argv)'

don't print messages like "ENTER SEQUENCE:"

read from stdin or from one or more file

read fasta

use switch/case instead of 'if/else' or even faster, a conversion array char compl[UCHAR_MAX];

a C program returns 0 on success

Answer 2

[Alex Reynolds]

If you want really fast, avoid conditionals and instead declare a static array that uses a character's ASCII decimal value to map to its complementary base:

This takes the reverse complement of a string made up of bases {A, C, T, G, U, a, c, t, g, u, N} by reading the input string backwards and printing mapped bases to standard output.

To compile:

$ gcc -Wall rc.c -o rc

To run, as an example:

$ echo 'AGGTCCA' | ./rc
TGGACCT

[lh3]

FYI: seqtk/bioawk has a more complete complement table. It includes R/Y/D/etc:

https://github.com/lh3/seqtk/blob/master/seqtk.c#L159

[Pierre Lindenbaum]

Thanks Alex, getline is new to me. Is it just a GNU extension or is it part of a recent standard ?

[Pierre Lindenbaum]

memory leak: line should be freed

[Alex Reynolds]

Fixed, thanks. Valgrind suggests the heap is clean:

$ echo -e 'ATTCG\nTTCCA\nGGGAT\nNNaTT' | valgrind -v --track-origins=yes --leak-check=full ./rc
==26625== Memcheck, a memory error detector
==26625== Copyright (C) 2002-2012, and GNU GPL'd, by Julian Seward et al.
==26625== Using Valgrind-3.8.1 and LibVEX; rerun with -h for copyright info
==26625== Command: ./rc
==26625== 
...
CGAAT
TGGAA
ATCCC
AAtNN
==26625== 
==26625== HEAP SUMMARY:
==26625==     in use at exit: 0 bytes in 0 blocks
==26625==   total heap usage: 1 allocs, 1 frees, 120 bytes allocated
==26625== 
==26625== All heap blocks were freed -- no leaks are possible
==26625== 
==26625== For counts of detected and suppressed errors, rerun with: -v
==26625== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 6 from 6)

[Pierre Lindenbaum]

answering my question: http://stackoverflow.com/questions/13112784 "First, getline() is not in the C standard library, but is a POSIX 2008 extension. Normally, it will be available with a POSIX-compatible compiler, as the macros _POSIX_C_SOURCE will be defined with the appropriate values. You possibly have an older compiler from before getline() was standardized, in which case this is a GNU extension, and you must #define _GNU_SOURCE before #include <stdio.h> to enable it, and must be using a GNU-compatible compiler, such as gcc."

[Alex Reynolds]

I didn't have to add an extension to compile under the version of clang that ships in OS X 10.11, so I guess OS X supports that POSIX standard, but on an older Linux box running gcc 4.8.2, I did have to add -std=gnu99 to the build statement. I did not have to add #define _GNU_SOURCE.

[Brian Bushnell]

Don't forget N and U! :)

[Alex Reynolds]

Done and done!

Answer 3

[Alternative]

Unless it is an assignment, I do not see why you won't use existing tools. Here is faRc from Kent utils and much more:

http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/

or

https://github.com/ENCODE-DCC/kentUtils