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
|
#include <iostream>
#include <seqan/sequence.h>
#include <seqan/seq_io.h>
#include <seqan/stream.h>
using namespace seqan2;
int main(int argc, char const ** argv)
{
if (argc != 5)
{
std::cerr << "USAGE: query_fai FILE.fa SEQ BEGIN END\n";
return 0;
}
// Try to load index and create on the fly if necessary.
FaiIndex faiIndex;
if (!open(faiIndex, argv[1]))
{
if (!build(faiIndex, argv[1]))
{
std::cerr << "ERROR: Index could not be loaded or built.\n";
return 0;
}
if (!save(faiIndex)) // Name is stored from when reading.
{
std::cerr << "ERROR: Index could not be written do disk.\n";
return 0;
}
}
// Translate sequence name to index.
unsigned idx = 0;
if (!getIdByName(idx, faiIndex, argv[2]))
{
std::cerr << "ERROR: Index does not know about sequence " << argv[2] << "\n";
return 0;
}
// Convert positions into integers.
unsigned beginPos = 0, endPos = 0;
if (!lexicalCast(beginPos, argv[3]))
{
std::cerr << "ERROR: Cannot cast " << argv[3] << " into an unsigned.\n";
return 0;
}
if (!lexicalCast(endPos, argv[4]))
{
std::cerr << "ERROR: Cannot cast " << argv[4] << " into an unsigned.\n";
return 0;
}
// Make sure begin and end pos are on the sequence and begin <= end.
if (beginPos > sequenceLength(faiIndex, idx))
beginPos = sequenceLength(faiIndex, idx);
if (endPos > sequenceLength(faiIndex, idx))
endPos = sequenceLength(faiIndex, idx);
if (beginPos > endPos)
endPos = beginPos;
// Finally, get infix of sequence.
Dna5String sequenceInfix;
readRegion(sequenceInfix, faiIndex, idx, beginPos, endPos);
std::cout << sequenceInfix << "\n";
return 0;
}
|
