|
0
|
1 /* reverseComplement.c
|
|
|
2 Program reading a FASTA file as input (one or multiple sequence)
|
|
|
3 and then returning the reverse complementary sequence in the output file.
|
|
|
4 It is important to note that for multiple sequence file the first input
|
|
|
5 sequence is the last one at the output file and the last one in the input
|
|
|
6 is the first one in the output.
|
|
|
7 oscart@uma.es
|
|
|
8 */
|
|
|
9
|
|
|
10 #include <stdio.h>
|
|
|
11 #include <stdlib.h>
|
|
|
12 #include <string.h>
|
|
|
13 #include <ctype.h>
|
|
|
14 #include <inttypes.h>
|
|
|
15 #include "commonFunctions.h"
|
|
|
16
|
|
|
17 #define SEQSIZE 2000000000
|
|
|
18 #define READINPUT 10000
|
|
|
19 #define WSIZE 32
|
|
|
20 #define NREADS 1000000
|
|
|
21
|
|
|
22 int main(int ac, char** av) {
|
|
|
23 FILE *fIn, *fOut;
|
|
|
24 int64_t i, j, nR, seqLen = 0;
|
|
|
25 char *seq, c, toW;
|
|
|
26 long *offset = NULL;
|
|
|
27
|
|
|
28 if (ac != 3)
|
|
|
29 terror("USE: reverseComplement seqFile.IN reverseComplementarySeq.OUT");
|
|
|
30
|
|
|
31 /**
|
|
|
32 * Allocating memory for the sequence
|
|
|
33 * Fixed amount of memory, change it to loadSeqDB?
|
|
|
34 */
|
|
|
35 if ((seq = (char*) malloc(SEQSIZE)) == NULL)
|
|
|
36 terror("memory for Seq");
|
|
|
37
|
|
|
38 if ((fIn = fopen(av[1], "rt")) == NULL)
|
|
|
39 terror("opening IN sequence FASTA file");
|
|
|
40
|
|
|
41 if ((fOut = fopen(av[2], "wt")) == NULL)
|
|
|
42 terror("opening OUT sequence Words file");
|
|
|
43
|
|
|
44 if ((offset = (long *) malloc(sizeof(long)*NREADS)) == NULL)
|
|
|
45 terror("memory for offsets");
|
|
|
46
|
|
|
47 for(i=0;i<NREADS;i++){
|
|
|
48 offset[i]=0;
|
|
|
49 }
|
|
|
50
|
|
|
51 nR = 0;
|
|
|
52 c = fgetc(fIn);
|
|
|
53 while(!feof(fIn)){
|
|
|
54 if(c == '>'){
|
|
|
55 offset[nR++] = ftell(fIn)-1;
|
|
|
56 }
|
|
|
57 c = fgetc(fIn);
|
|
|
58 }
|
|
|
59
|
|
|
60 for(i=nR-1; i>=0; i--){
|
|
|
61 fseek(fIn, offset[i], SEEK_SET);
|
|
|
62 //READ and write header
|
|
|
63 fgets(seq, READINPUT, fIn);
|
|
|
64
|
|
|
65 fprintf(fOut, "%s", seq);
|
|
|
66 //READ and write sequence
|
|
|
67 c = fgetc(fIn);
|
|
|
68 while(c != '>' && !feof(fIn)){
|
|
|
69 if(isupper(c) || islower(c)){
|
|
|
70 seq[seqLen++]=c;
|
|
|
71 }
|
|
|
72 c = fgetc(fIn);
|
|
|
73 }
|
|
|
74 for(j=seqLen-1; j >= 0; j--){
|
|
|
75 switch(seq[j]){
|
|
|
76 case 'A':
|
|
|
77 toW = 'T';
|
|
|
78 break;
|
|
|
79 case 'C':
|
|
|
80 toW = 'G';
|
|
|
81 break;
|
|
|
82 case 'G':
|
|
|
83 toW = 'C';
|
|
|
84 break;
|
|
|
85 case 'T':
|
|
|
86 toW = 'A';
|
|
|
87 break;
|
|
|
88 case 'U':
|
|
|
89 toW = 'A';
|
|
|
90 break;
|
|
|
91 case 'a':
|
|
|
92 toW = 't';
|
|
|
93 break;
|
|
|
94 case 'c':
|
|
|
95 toW = 'g';
|
|
|
96 break;
|
|
|
97 case 'g':
|
|
|
98 toW = 'c';
|
|
|
99 break;
|
|
|
100 case 't':
|
|
|
101 toW = 'a';
|
|
|
102 break;
|
|
|
103 case 'u':
|
|
|
104 toW = 'a';
|
|
|
105 break;
|
|
|
106 default:
|
|
|
107 toW = seq[j];
|
|
|
108 break;
|
|
|
109 }
|
|
|
110 fwrite(&toW, sizeof(char), 1, fOut);
|
|
|
111 }
|
|
|
112 toW='\n';
|
|
|
113 fwrite(&toW, sizeof(char), 1, fOut);
|
|
|
114 seqLen=0;
|
|
|
115 }
|
|
|
116
|
|
|
117 fclose(fIn);
|
|
|
118 fclose(fOut);
|
|
|
119
|
|
|
120 return 0;
|
|
|
121 }
|
|
|
122
|
|
|
123
|
