C program to find complementary of DNA sequence
3
1
Entering edit mode
5.9 years ago
Arindam Ghosh ▴ 400

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

}

myposts c dna • 21k views
ADD COMMENT
2
Entering edit mode

Inputting DNA sequence by hand.. very convenient.

ADD REPLY
1
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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]
ADD REPLY
0
Entering edit mode

A 1-liner might look like:

from string import maketrans
print raw_input('DNA: ').translate(maketrans('AaCcGgTt', 'TtGgCcAa'))[::-1]
ADD REPLY
5
Entering edit mode
5.9 years ago
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

ADD COMMENT
2
Entering edit mode
5.9 years ago

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
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

memory leak: line should be freed

ADD REPLY
0
Entering edit mode

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)
ADD REPLY
0
Entering edit mode

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."

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

Don't forget N and U! :)

ADD REPLY
1
Entering edit mode

Done and done!

ADD REPLY
0
Entering edit mode

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

ADD REPLY
1
Entering edit mode
5.9 years ago
Alternative ▴ 260

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

ADD COMMENT

Login before adding your answer.

Traffic: 1561 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6