1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
|
/*****************************************************************
* @LICENSE@
*****************************************************************/
/* revcomp.c
*
* Reverse complement of a IUPAC character string
* CVS $Id: revcomp.c,v 1.6 2003/04/14 16:00:16 eddy Exp $
*/
#include "squidconf.h"
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include "squid.h"
/* Function: revcomp()
*
* Purpose: Reverse complement seq; store in comp.
* Can revcomp "in place" (revcomp(seq, seq)).
*
* Args: comp - destination for reverse complement of seq
* seq - sequence to reverse complement
*
* Returns: NULL on failure; or a (useless) pointer to comp.
*/
char *
revcomp(char *comp, char *seq)
{
char *s;
char c;
if (comp == NULL) return NULL;
if (seq == NULL) return NULL;
StrReverse(comp, seq);
for (s = comp; *s != '\0'; s++)
{
c = *s;
c = sre_toupper(c);
switch (c) {
case 'A': c = 'T'; break;
case 'C': c = 'G'; break;
case 'G': c = 'C'; break;
case 'T': c = 'A'; break;
case 'U': c = 'A'; break;
case 'R': c = 'Y'; break;
case 'Y': c = 'R'; break;
case 'M': c = 'K'; break;
case 'K': c = 'M'; break;
case 'S': c = 'S'; break;
case 'W': c = 'W'; break;
case 'H': c = 'D'; break;
case 'D': c = 'H'; break;
case 'B': c = 'V'; break;
case 'V': c = 'B'; break;
default: break; /* anything else? leave it; it's prob a gap or an X */
}
if (islower((int) *s)) c = (char) sre_tolower((int) c);
*s = c;
}
return comp;
}
#ifdef REVCOMP_TESTDRIVER
/* gcc -g -DREVCOMP_TESTDRIVER revcomp.c sre_string.c shuffle.c sre_math.c sre_ctype.c sqerror.c -lm
*/
int
main(void)
{
float p[4] = {0.25, 0.25, 0.25, 0.25};
char *alphabet = "ACGT";
int len = 10;
char *seq;
seq = RandomSequence(alphabet, p, 4, len);
printf("%s\n", seq);
revcomp(seq, seq);
printf("%s\n", seq);
free(seq);
exit(0);
}
#endif
|
