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

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

dna c • 30k views
2
Entering edit mode

Inputting DNA sequence by hand.. very convenient.

1
Entering edit mode

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

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]

0
Entering edit mode

A 1-liner might look like:

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

5
Entering edit mode
8.0 years ago
declare 'main(int argc,char** argv)'


don't print messages like "ENTER SEQUENCE:"

read from stdin or from one or more file

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

a C program returns 0 on success

3
Entering edit mode
8.0 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
1
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

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 ?

0
Entering edit mode

memory leak: line should be freed

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

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.

0
Entering edit mode

Don't forget N and U! :)

1
Entering edit mode

Done and done!

2
Entering edit mode
8.0 years ago
Alternative ▴ 270

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:

or

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