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
|
/*
* hit_set.cpp
*
* Created on: Jul 31, 2009
* Author: Ben Langmead
*/
#include <iostream>
#include <vector>
#include <seqan/sequence.h>
#include "alphabet.h"
#include "hit_set.h"
using namespace std;
using namespace seqan;
/**
* Report up to 'khits' hits from this HitSet.
*/
void HitSet::reportUpTo(ostream& os, int khits) {
khits = min(khits, (int)size());
String<Dna5> seqrc;
String<char> qualr;
for(int i = 0; i < khits; i++) {
const HitSetEnt& h = ents[i];
if(!h.fw && seqan::empty(seqrc)) {
// Lazily initialize seqrc and qualr
seqrc = seq;
reverseComplementInPlace(seqrc, color);
assert_eq(seqan::length(seqrc), seqan::length(seq));
qualr = qual;
reverseInPlace(qualr);
assert_eq(seqan::length(qualr), seqan::length(qual));
}
os << name << '\t'
<< (h.fw ? '+' : '-') << '\t'
<< h.h.first << '\t'
<< h.h.second << '\t'
<< (h.fw ? seq : seqrc) << '\t'
<< (h.fw ? qual : qualr) << '\t'
<< h.oms << '\t';
for(size_t i = 0; i < h.edits.size(); i++) {
const Edit& e = h.edits[i];
os << e.pos;
if(e.type == EDIT_TYPE_SNP) os << "S";
os << ":" << (char)e.chr << ">" << (e.qchr != 0 ? (char)e.qchr : (char)seq[e.pos]);
if(i < h.edits.size()-1 || !h.cedits.empty()) os << ",";
}
for(size_t i = 0; i < h.cedits.size(); i++) {
const Edit& e = h.cedits[i];
os << e.pos;
if(e.type == EDIT_TYPE_SNP) os << "S";
os << ":" << (char)e.chr << ">" << (e.qchr != 0 ? (char)e.qchr : (char)seq[e.pos]);
if(i < h.cedits.size()-1) os << ",";
}
os << endl;
}
}
ostream& operator << (ostream& os, const HitSetEnt& hs) {
os << "\t" << hs.h.first << ":" << hs.h.second;
return os;
}
ostream& operator << (ostream& os, const HitSet& hs) {
os << hs.name << ":" << hs.seq << ":" << hs.qual << endl;
vector<HitSetEnt>::const_iterator it;
for(it = hs.ents.begin(); it != hs.ents.end(); it++) {
os << (*it);
}
return os;
}
|
