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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
|
#include <seqan/sequence.h>
#include <seqan/bam_io.h>
using namespace seqan2;
int main(int, char const **)
{
CharString bamFileName = getAbsolutePath("demos/tutorial/sam_and_bam_io/example.bam");
CharString baiFileName = getAbsolutePath("demos/tutorial/sam_and_bam_io/example.bam.bai");
CharString rName = "ref";
// Open BamFileIn for reading.
BamFileIn inFile;
if (!open(inFile, toCString(bamFileName)))
{
std::cerr << "ERROR: Could not open " << bamFileName << " for reading.\n";
return 1;
}
// Read BAI index.
BamIndex<Bai> baiIndex;
if (!open(baiIndex, toCString(baiFileName)))
{
std::cerr << "ERROR: Could not read BAI index file " << baiFileName << "\n";
return 1;
}
// Read header.
BamHeader header;
readHeader(header, inFile);
// Translate from reference name to rID.
int rID = 0;
if (!getIdByName(rID, contigNamesCache(context(inFile)), rName))
{
std::cerr << "ERROR: Reference sequence named " << rName << " not known.\n";
return 1;
}
// Translate BEGIN and END arguments to number, 1-based to 0-based.
int beginPos = 9, endPos = 30;
// 1-based to 0-based.
beginPos -= 1;
endPos -= 1;
// Translate number of elements to print to number.
int num = 3;
// Jump the BGZF stream to this position.
bool hasAlignments = false;
if (!jumpToRegion(inFile, hasAlignments, rID, beginPos, endPos, baiIndex))
{
std::cerr << "ERROR: Could not jump to " << beginPos << ":" << endPos << "\n";
return 1;
}
if (!hasAlignments)
return 0; // No alignments here.
// Seek linearly to the selected position.
BamAlignmentRecord record;
int numPrinted = 0;
BamFileOut out(inFile, std::cout, Sam());
while (!atEnd(inFile) && numPrinted < num)
{
readRecord(record, inFile);
// If we are on the next reference or at the end already then we stop.
if (record.rID == -1 || record.rID > rID || record.beginPos >= endPos)
break;
// If we are left of the selected position then we skip this record.
if (record.beginPos < beginPos)
continue;
// Otherwise, we print it to the user.
numPrinted++;
writeRecord(out, record);
}
return 0;
}
|
