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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
|
#ifndef HITCONTAINER_H_
#define HITCONTAINER_H_
#include<cassert>
#include<iostream>
#include<vector>
#include<algorithm>
#include "utils.h"
#include "GroupInfo.h"
template<class HitType>
class HitContainer {
public:
HitContainer() {
clear();
}
void clear() {
n = nhits = 0;
s.clear();
hits.clear();
s.push_back(0);
}
bool read(std::istream&); // each time a read
void write(std::ostream&); // write all reads' hit out
void push_back(const HitType& hit) {
hits.push_back(hit);
++nhits;
}
//update read information vector etc
void updateRI() {
if (nhits > s.back()) { //Do not change if last read does not have hits
s.push_back(nhits);
++n;
}
}
READ_INT_TYPE getN() { return n; }
HIT_INT_TYPE getNHits() { return nhits; }
READ_INT_TYPE calcNumGeneMultiReads(const GroupInfo&);
READ_INT_TYPE calcNumIsoformMultiReads();
HIT_INT_TYPE getSAt(READ_INT_TYPE pos) { assert(pos >= 0 && pos <= n); return s[pos]; }
HitType& getHitAt(HIT_INT_TYPE pos) { assert(pos >= 0 && pos < nhits); return hits[pos]; }
private:
READ_INT_TYPE n; // n reads in total
HIT_INT_TYPE nhits; // # of hits
std::vector<HIT_INT_TYPE> s;
std::vector<HitType> hits;
};
//Each time only read one read's hits. If you want to start over, must call clear() first!
template<class HitType>
bool HitContainer<HitType>::read(std::istream& in) {
HIT_INT_TYPE tot;
if (!(in>>tot)) return false;
assert(tot > 0);
for (HIT_INT_TYPE i = 0; i < tot; i++) {
HitType hit;
if (!hit.read(in)) return false;
hits.push_back(hit);
}
nhits = nhits + tot;
++n;
s.push_back(nhits);
return true;
}
template<class HitType>
void HitContainer<HitType>::write(std::ostream& out) {
if (n <= 0) return;
for (READ_INT_TYPE i = 0; i < n; i++) {
out<<s[i + 1] - s[i];
for (HIT_INT_TYPE j = s[i]; j < s[i + 1]; j++) {
hits[j].write(out);
}
out<<std::endl;
}
}
template<class HitType>
READ_INT_TYPE HitContainer<HitType>::calcNumGeneMultiReads(const GroupInfo& gi) {
READ_INT_TYPE res = 0;
int *sortgids = NULL;
for (READ_INT_TYPE i = 0; i < n; i++) {
HIT_INT_TYPE num = s[i + 1] - s[i];
sortgids = new int[num];
for (HIT_INT_TYPE j = s[i]; j < s[i + 1]; j++) sortgids[j] = gi.gidAt(hits[j].getSid());
std::sort(sortgids, sortgids + num);
if (std::unique(sortgids, sortgids + num) - sortgids > 1) ++res;
delete[] sortgids;
}
return res;
}
template<class HitType>
READ_INT_TYPE HitContainer<HitType>::calcNumIsoformMultiReads() {
READ_INT_TYPE res = 0;
for (READ_INT_TYPE i = 0; i < n; i++)
if (s[i + 1] - s[i] > 1) ++res;
return res;
}
#endif /* HITCONTAINER_H_ */
|
