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
|
#include <fstream>
#include <iostream>
#include <string>
#include <seqan/arg_parse.h>
#include <seqan/store.h>
#include "razers.h"
using namespace seqan2;
int main(int argc, const char * argv[])
{
//////////////////////////////////////////////////////////////////////////////
// Define options
ArgumentParser parser;
RazerSOptions<> options;
double delta = 0;
addUsageLine(parser, "[OPTION]... <reads.fasta>");
addArgument(parser, ArgParseArgument(ArgParseArgument::INPUT_FILE, "READS FILE"));
addOption(parser, ArgParseOption("mr", "mutation-rate", "set the percent mutation rate", ArgParseOption::DOUBLE));
addOption(parser, ArgParseOption("qd", "quality-delta", "add a delta value to qualities", ArgParseOption::DOUBLE));
ArgumentParser::ParseResult res = parse(parser, argc, argv);
if (res != ArgumentParser::PARSE_OK)
return res == ArgumentParser::PARSE_ERROR;
FragmentStore<> store;
getOptionValue(options.mutationRate, parser, "mutation-rate");
getOptionValue(delta, parser, "quality-delta");
options.mutationRate /= 100;
//////////////////////////////////////////////////////////////////////////////
// Load reads
CharString readsFilename;
getArgumentValue(readsFilename, parser, 0);
SeqFileIn seqFileIn;
if (!open(seqFileIn, toCString(readsFilename)) || !loadReads(store, seqFileIn, options))
{
std::cerr << "Failed to load reads." << std::endl;
return 1;
}
//////////////////////////////////////////////////////////////////////////////
// Output error distribution
for (unsigned i = 0; i < length(options.avrgQuality); ++i)
{
double sequencingError = exp((options.avrgQuality[i] + delta) * log(10.0) / -10.0);
double errorProb = 1.0 - (1.0 - sequencingError) * (1.0 - options.mutationRate);
// std::cout << options.avrgQuality[i] << std::endl;
std::cout << i << '\t' << errorProb << std::endl;
}
return 0;
}
|
