|
0
|
1 /*
|
|
|
2 * Copyright (c) <2008 - 2009>, University of Washington, Simon Fraser University
|
|
|
3 * All rights reserved.
|
|
|
4 *
|
|
|
5 * Redistribution and use in source and binary forms, with or without modification,
|
|
|
6 * are permitted provided that the following conditions are met:
|
|
|
7 *
|
|
|
8 * Redistributions of source code must retain the above copyright notice, this list
|
|
|
9 * of conditions and the following disclaimer.
|
|
|
10 * - Redistributions in binary form must reproduce the above copyright notice, this
|
|
|
11 * list of conditions and the following disclaimer in the documentation and/or other
|
|
|
12 * materials provided with the distribution.
|
|
|
13 * - Neither the name of the <ORGANIZATION> nor the names of its contributors may be
|
|
|
14 * used to endorse or promote products derived from this software without specific
|
|
|
15 * prior written permission.
|
|
|
16 *
|
|
|
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
|
|
18 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
|
|
19 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
|
|
|
20 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
|
|
|
21 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
|
|
|
22 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
|
|
|
23 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
|
|
|
24 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
|
|
|
25 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
|
|
|
26 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
|
|
|
27 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
|
|
28 */
|
|
|
29
|
|
|
30 /*
|
|
|
31 * Author : Faraz Hach
|
|
|
32 * Email : fhach AT cs DOT sfu
|
|
|
33 * Last Update : 2009-01-29
|
|
|
34 */
|
|
|
35
|
|
|
36 #include <stdio.h>
|
|
|
37 #include <stdlib.h>
|
|
|
38 #include <getopt.h>
|
|
|
39 #include <string.h>
|
|
|
40 #include <ctype.h>
|
|
|
41 #include "Common.h"
|
|
|
42 #include "CommandLineParser.h"
|
|
|
43
|
|
|
44 int uniqueMode=1;
|
|
|
45 int indexingMode;
|
|
|
46 int searchingMode;
|
|
|
47 int bisulfiteMode;
|
|
|
48 int pairedEndMode;
|
|
|
49 int pairedEndDiscordantMode;
|
|
|
50 int pairedEndProfilingMode;
|
|
|
51 int seqCompressed;
|
|
|
52 int outCompressed;
|
|
|
53 int cropSize = 0;
|
|
|
54 int progressRep = 0;
|
|
|
55 int minPairEndedDistance=-1;
|
|
|
56 int maxPairEndedDistance=-1;
|
|
|
57 int minPairEndedDiscordantDistance=-1;
|
|
|
58 int maxPairEndedDiscordantDistance=-1;
|
|
|
59 char *seqFile1;
|
|
|
60 char *seqFile2;
|
|
|
61 char *mappingOutput = "output";
|
|
|
62 char *mappingOutputPath = "";
|
|
|
63 char *unmappedOutput = "";
|
|
|
64 char fileName[1000][2][FILE_NAME_LENGTH];
|
|
|
65 int fileCnt;
|
|
|
66 unsigned char errThreshold=2;
|
|
|
67 unsigned char maxHits=0;
|
|
|
68 unsigned char WINDOW_SIZE = 12;
|
|
|
69 unsigned int CONTIG_SIZE;
|
|
|
70 unsigned int CONTIG_MAX_SIZE;
|
|
|
71
|
|
|
72 void printHelp();
|
|
|
73
|
|
|
74 int parseCommandLine (int argc, char *argv[])
|
|
|
75 {
|
|
|
76
|
|
|
77 int o;
|
|
|
78 int index;
|
|
|
79 char *fastaFile = NULL;
|
|
|
80 char *fastaOutputFile = NULL;
|
|
|
81 char *indexFile = NULL;
|
|
|
82 char *batchFile = NULL ;
|
|
|
83 int batchMode = 0;
|
|
|
84 static struct option longOptions[] =
|
|
|
85 {
|
|
|
86
|
|
|
87 // {"bs", no_argument, &bisulfiteMode, 1},
|
|
|
88 {"pe", no_argument, &pairedEndMode, 1},
|
|
|
89 {"discordant-vh", no_argument, &pairedEndDiscordantMode, 1},
|
|
|
90 {"profile", no_argument, &pairedEndProfilingMode, 1},
|
|
|
91 {"seqcomp", no_argument, &seqCompressed, 1},
|
|
|
92 {"outcomp", no_argument, &outCompressed, 1},
|
|
|
93 {"progress", no_argument, &progressRep, 1},
|
|
|
94 {"index", required_argument, 0, 'i'},
|
|
|
95 {"search", required_argument, 0, 's'},
|
|
|
96 {"help", no_argument, 0, 'h'},
|
|
|
97 {"version", no_argument, 0, 'v'},
|
|
|
98 {"seq", required_argument, 0, 'x'},
|
|
|
99 {"seq1", required_argument, 0, 'x'},
|
|
|
100 {"seq2", required_argument, 0, 'y'},
|
|
|
101 {"ws", required_argument, 0, 'w'},
|
|
|
102 {"min", required_argument, 0, 'l'},
|
|
|
103 {"max", required_argument, 0, 'm'},
|
|
|
104 {"crop", required_argument, 0, 'c'}
|
|
|
105
|
|
|
106 };
|
|
|
107
|
|
|
108
|
|
|
109
|
|
|
110 while ( (o = getopt_long ( argc, argv, "f:i:u:o:s:e:n:bhv", longOptions, &index))!= -1 )
|
|
|
111 {
|
|
|
112 switch (o)
|
|
|
113 {
|
|
|
114 case 'i':
|
|
|
115 indexingMode = 1;
|
|
|
116 fastaFile = optarg;
|
|
|
117 break;
|
|
|
118 case 's':
|
|
|
119 searchingMode = 1;
|
|
|
120 fastaFile = optarg;
|
|
|
121 break;
|
|
|
122 case 'b':
|
|
|
123 batchMode = 1;
|
|
|
124 break;
|
|
|
125 case 'c':
|
|
|
126 cropSize = atoi(optarg);
|
|
|
127 break;
|
|
|
128 case 'w':
|
|
|
129 WINDOW_SIZE = atoi(optarg);
|
|
|
130 break;
|
|
|
131 case 'x':
|
|
|
132 seqFile1 = optarg;
|
|
|
133 break;
|
|
|
134 case 'y':
|
|
|
135 seqFile2 = optarg;
|
|
|
136 break;
|
|
|
137 case 'u':
|
|
|
138 unmappedOutput = optarg;
|
|
|
139 break;
|
|
|
140 case 'o':
|
|
|
141 mappingOutput = getMem(FILE_NAME_LENGTH);
|
|
|
142 mappingOutputPath = getMem(FILE_NAME_LENGTH);
|
|
|
143 stripPath (optarg, &mappingOutputPath, &mappingOutput);
|
|
|
144 break;
|
|
|
145 case 'n':
|
|
|
146 maxHits = atoi(optarg);
|
|
|
147 break;
|
|
|
148 case 'e':
|
|
|
149 errThreshold = atoi(optarg);
|
|
|
150 break;
|
|
|
151 case 'l':
|
|
|
152 minPairEndedDistance = atoi(optarg);
|
|
|
153 break;
|
|
|
154 case 'm':
|
|
|
155 maxPairEndedDistance = atoi(optarg);
|
|
|
156 break;
|
|
|
157 case 'h':
|
|
|
158 printHelp();
|
|
|
159 return 0;
|
|
|
160 break;
|
|
|
161 case 'v':
|
|
|
162 fprintf(stdout, "%s.%s\n", versionNumber, versionNumberF);
|
|
|
163 return 0;
|
|
|
164 break;
|
|
|
165 }
|
|
|
166
|
|
|
167 }
|
|
|
168
|
|
|
169 if (indexingMode + searchingMode != 1)
|
|
|
170 {
|
|
|
171 fprintf(stdout, "ERROR: Indexing / Searching mode should be selected\n");
|
|
|
172 return 0;
|
|
|
173 }
|
|
|
174
|
|
|
175 if (WINDOW_SIZE > 14 || WINDOW_SIZE < 8)
|
|
|
176 {
|
|
|
177 fprintf(stdout, "ERROR: Window size should be in [8..14]\n");
|
|
|
178 return 0;
|
|
|
179 }
|
|
|
180
|
|
|
181
|
|
|
182 if ( indexingMode )
|
|
|
183 {
|
|
|
184 CONTIG_SIZE = 15000000;
|
|
|
185 CONTIG_MAX_SIZE = 40000000;
|
|
|
186
|
|
|
187 if (batchMode)
|
|
|
188 {
|
|
|
189 batchFile = fastaFile;
|
|
|
190 fastaFile = NULL;
|
|
|
191 }
|
|
|
192
|
|
|
193 if (batchFile == NULL && fastaFile == NULL)
|
|
|
194 {
|
|
|
195 fprintf(stdout, "ERROR: Reference(s) should be indicated for indexing\n");
|
|
|
196 return 0;
|
|
|
197 }
|
|
|
198
|
|
|
199 if (pairedEndDiscordantMode)
|
|
|
200 {
|
|
|
201 fprintf(stdout, "ERROR: --discordant cannot be used in indexing mode. \n");
|
|
|
202 return 0;
|
|
|
203 }
|
|
|
204
|
|
|
205 }
|
|
|
206
|
|
|
207
|
|
|
208 if ( searchingMode )
|
|
|
209 {
|
|
|
210 CONTIG_SIZE = 300000000;
|
|
|
211 CONTIG_MAX_SIZE = 300000000;
|
|
|
212
|
|
|
213
|
|
|
214 if (batchMode)
|
|
|
215 {
|
|
|
216 batchFile = fastaFile;
|
|
|
217 fastaFile = NULL;
|
|
|
218 }
|
|
|
219
|
|
|
220 if (batchFile == NULL && fastaFile == NULL)
|
|
|
221 {
|
|
|
222 fprintf(stdout, "ERROR: Index File(s) should be indiciated for searching\n");
|
|
|
223 return 0;
|
|
|
224 }
|
|
|
225
|
|
|
226 if (seqFile1 == NULL && seqFile2 == NULL)
|
|
|
227 {
|
|
|
228 fprintf(stdout, "ERROR: Please indicate a sequence file for searching.\n");
|
|
|
229 return 0;
|
|
|
230 }
|
|
|
231
|
|
|
232
|
|
|
233 if (!pairedEndMode && !bisulfiteMode && seqFile2 != NULL)
|
|
|
234 {
|
|
|
235 fprintf(stdout, "ERROR: Second File can be indicated in pairedend/bisulfite mode\n");
|
|
|
236 return 0;
|
|
|
237 }
|
|
|
238
|
|
|
239 if (pairedEndMode && (minPairEndedDistance <0 || maxPairEndedDistance < 0 || minPairEndedDistance > maxPairEndedDistance))
|
|
|
240 {
|
|
|
241 fprintf(stdout, "ERROR: Please enter a valid range for pairedend sequences.\n");
|
|
|
242 return 0;
|
|
|
243 }
|
|
|
244
|
|
|
245 if (pairedEndMode && seqFile1 == NULL)
|
|
|
246 {
|
|
|
247 fprintf(stdout, "ERROR: Please indicate the first file for pairedend search.\n");
|
|
|
248 return 0;
|
|
|
249 }
|
|
|
250
|
|
|
251 if (!pairedEndMode && pairedEndDiscordantMode)
|
|
|
252 {
|
|
|
253 fprintf(stdout, "ERROR: --discordant should be used with --pe");
|
|
|
254 return 0;
|
|
|
255 }
|
|
|
256
|
|
|
257 if (!pairedEndMode && pairedEndProfilingMode)
|
|
|
258 {
|
|
|
259 fprintf(stdout, "ERROR: --profile should be used with --pe");
|
|
|
260 return 0;
|
|
|
261 }
|
|
|
262 }
|
|
|
263
|
|
|
264 int i = 0;
|
|
|
265
|
|
|
266
|
|
|
267 if (batchMode)
|
|
|
268 {
|
|
|
269 FILE *fp = fileOpen(batchFile, "r");
|
|
|
270
|
|
|
271 if (fp == NULL)
|
|
|
272 return 0;
|
|
|
273
|
|
|
274 fileCnt = 0;
|
|
|
275
|
|
|
276 while ( fgets(fileName[fileCnt][0], FILE_NAME_LENGTH, fp))
|
|
|
277 {
|
|
|
278 for (i = strlen(fileName[fileCnt][0])-1; i>=0; i--)
|
|
|
279 if ( !isspace(fileName[fileCnt][0][i]))
|
|
|
280 break;
|
|
|
281 fileName[fileCnt][0][i+1] = '\0';
|
|
|
282
|
|
|
283 if (strcmp(fileName[fileCnt][0], "") != 0)
|
|
|
284 {
|
|
|
285 if (bisulfiteMode)
|
|
|
286 sprintf(fileName[fileCnt][1], "%s.bsindex", fileName[fileCnt][0]);
|
|
|
287 else
|
|
|
288 sprintf(fileName[fileCnt][1], "%s.index", fileName[fileCnt][0]);
|
|
|
289 fileCnt++;
|
|
|
290 }
|
|
|
291 }
|
|
|
292 }
|
|
|
293 else
|
|
|
294 {
|
|
|
295 sprintf(fileName[fileCnt][0], "%s", fastaFile);
|
|
|
296 if (bisulfiteMode)
|
|
|
297 sprintf(fileName[fileCnt][1], "%s.bsindex", fileName[fileCnt][0]);
|
|
|
298 else
|
|
|
299 sprintf(fileName[fileCnt][1], "%s.index", fileName[fileCnt][0]);
|
|
|
300 fileCnt++;
|
|
|
301 }
|
|
|
302
|
|
|
303
|
|
|
304 if (pairedEndProfilingMode)
|
|
|
305 {
|
|
|
306
|
|
|
307 minPairEndedDistance = 0;
|
|
|
308 maxPairEndedDistance = 300000000;
|
|
|
309
|
|
|
310 }
|
|
|
311
|
|
|
312 if (pairedEndDiscordantMode)
|
|
|
313 {
|
|
|
314 minPairEndedDiscordantDistance = minPairEndedDistance;
|
|
|
315 maxPairEndedDiscordantDistance = maxPairEndedDistance;
|
|
|
316
|
|
|
317 minPairEndedDistance = 0;
|
|
|
318 maxPairEndedDistance = 300000000;
|
|
|
319 }
|
|
|
320
|
|
|
321 return 1;
|
|
|
322 }
|
|
|
323
|
|
|
324
|
|
|
325 void printHelp()
|
|
|
326 {
|
|
|
327 char *errorType;
|
|
|
328 if (mrFAST)
|
|
|
329 {
|
|
|
330 fprintf(stdout,"mrFAST : Micro-Read Fast Alignment Search Tool.\n\n");
|
|
|
331 fprintf(stdout,"Usage: mrFAST [options]\n\n");
|
|
|
332 errorType="edit distance";
|
|
|
333 }
|
|
|
334 else
|
|
|
335 {
|
|
|
336 fprintf(stdout,"mrsFAST : Micro-Read Substitutions (only) Fast Alignment Search Tool.\n\n");
|
|
|
337 fprintf(stdout,"mrsFAST is a cache oblivious read mapping tool. mrsFAST capable of mapping\n");
|
|
|
338 fprintf(stdout,"single and paired end reads to the reference genome. Bisulfite treated \n");
|
|
|
339 fprintf(stdout,"sequences are not supported in this version. By default mrsFAST reports \n");
|
|
|
340 fprintf(stdout,"the output in SAM format.\n\n");
|
|
|
341 fprintf(stdout,"Usage: mrsFAST [options]\n\n");
|
|
|
342 errorType="hamming distance";
|
|
|
343 }
|
|
|
344
|
|
|
345 fprintf(stdout,"General Options:\n");
|
|
|
346 fprintf(stdout," -v|--version\t\tCurrent Version.\n");
|
|
|
347 fprintf(stdout," -h\t\t\tShows the help file.\n");
|
|
|
348 fprintf(stdout,"\n\n");
|
|
|
349
|
|
|
350 fprintf(stdout,"Indexing Options:\n");
|
|
|
351 fprintf(stdout," --index [file]\t\tGenerate an index from the specified fasta file. \n");
|
|
|
352 fprintf(stdout," -b\t\t\tIndicates the indexing will be done in batch mode.\n\t\t\tThe file specified in --index should contain the \n\t\t\tlist of fasta files.\n");
|
|
|
353 fprintf(stdout," -ws [int]\t\tSet window size for indexing (default:12-min:8 max:14).\n");
|
|
|
354 // fprintf(stdout," -bs \t\t\tBisulfite mode.");
|
|
|
355 fprintf(stdout,"\n\n");
|
|
|
356
|
|
|
357 fprintf(stdout,"Searching Options:\n");
|
|
|
358 fprintf(stdout," --search [file]\tSearch the specified genome. Index file should be \n\t\t\tin same directory as the fasta file.\n");
|
|
|
359 fprintf(stdout," -b\t\t\tIndicates the mapping will be done in batch mode. \n\t\t\tThe file specified in --search should contain the \n\t\t\tlist of fasta files.\n");
|
|
|
360 fprintf(stdout," --pe \t\t\tSearch will be done in Pairedend mode.\n");
|
|
|
361 // fprintf(stdout," --bs \t\t\tSearch will be done in Bisulfite mode.\n");
|
|
|
362 fprintf(stdout," --seq [file]\t\tInput sequences in fasta/fastq format [file]. If \n\t\t\tpairend reads are interleaved, use this option.\n");
|
|
|
363 fprintf(stdout," --seq1 [file]\t\tInput sequences in fasta/fastq format [file] (First \n\t\t\tfile). Use this option to indicate the first file of \n\t\t\tpair-end reads. You can use this option alone in \n\t\t\tbisulfite mode. \n");
|
|
|
364 fprintf(stdout," --seq2 [file]\t\tInput sequences in fasta/fastq format [file] (Second \n\t\t\tfile). Use this option to indicate the second file of \n\t\t\tpair-end reads. You can use this option alone in \n\t\t\tbisulfite mode. \n");
|
|
|
365 fprintf(stdout," -o [file]\t\tOutput of the mapped sequences. The default is output.\n");
|
|
|
366 fprintf(stdout," --seqcomp \t\tIndicates that the input sequences are compressed(gz).\n");
|
|
|
367 fprintf(stdout," --outcomp \t\tIndicates that output file should be compressed(gz).\n");
|
|
|
368 // fprintf(stdout," -u [file]\t\tSave unmapped sequences to the [file] in fasta format.\n");
|
|
|
369 fprintf(stdout," -n [int]\t\tMaximum number of locations reported for a sequence \n\t\t\t(default 0, all mappings). \n");
|
|
|
370 fprintf(stdout," -e [int]\t\t%s (default 2).\n", errorType);
|
|
|
371 fprintf(stdout," --min [int]\t\tMin inferred distance allowed between two pairend sequences.\n");
|
|
|
372 fprintf(stdout," --max [int]\t\tMax inferred distance allowed between two pairend sequences.\n");
|
|
|
373 }
|
