1 #include <fstream>
2 #include <iostream>
3 #include <string>
4
5 #include <seqan/arg_parse.h>
6 #include <seqan/store.h>
7 #include "razers.h"
8
9 using namespace seqan;
10
11
main(int argc,const char * argv[])12 int main(int argc, const char * argv[])
13 {
14
15 //////////////////////////////////////////////////////////////////////////////
16 // Define options
17 ArgumentParser parser;
18 RazerSOptions<> options;
19 double delta = 0;
20
21 addUsageLine(parser, "[OPTION]... <reads.fasta>");
22 addArgument(parser, ArgParseArgument(ArgParseArgument::INPUT_FILE, "READS FILE"));
23 addOption(parser, ArgParseOption("mr", "mutation-rate", "set the percent mutation rate", ArgParseOption::DOUBLE));
24 addOption(parser, ArgParseOption("qd", "quality-delta", "add a delta value to qualities", ArgParseOption::DOUBLE));
25
26 ArgumentParser::ParseResult res = parse(parser, argc, argv);
27
28 if (res != ArgumentParser::PARSE_OK)
29 return res == ArgumentParser::PARSE_ERROR;
30
31 FragmentStore<> store;
32 getOptionValue(options.mutationRate, parser, "mutation-rate");
33 getOptionValue(delta, parser, "quality-delta");
34 options.mutationRate /= 100;
35
36 //////////////////////////////////////////////////////////////////////////////
37 // Load reads
38 CharString readsFilename;
39 getArgumentValue(readsFilename, parser, 0);
40
41 SeqFileIn seqFileIn;
42 if (!open(seqFileIn, toCString(readsFilename)) || !loadReads(store, seqFileIn, options))
43 {
44 std::cerr << "Failed to load reads." << std::endl;
45 return 1;
46 }
47
48 //////////////////////////////////////////////////////////////////////////////
49 // Output error distribution
50
51 for (unsigned i = 0; i < length(options.avrgQuality); ++i)
52 {
53 double sequencingError = exp((options.avrgQuality[i] + delta) * log(10.0) / -10.0);
54 double errorProb = 1.0 - (1.0 - sequencingError) * (1.0 - options.mutationRate);
55
56 // std::cout << options.avrgQuality[i] << std::endl;
57 std::cout << i << '\t' << errorProb << std::endl;
58 }
59 return 0;
60 }
61