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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
|
#include<cmath>
#include "TagAlignments.h"
#include "misc.h"
#include "common.h"
//#define MEM_USAGE
TagAlignments::TagAlignments(bool storeL){//{{{
knowNtotal=false;
knowNreads=false;
Ntotal=0;
Nreads=0;
storeLog = storeL;
}//}}}
void TagAlignments::init(long Nreads,long Ntotal, long M){//{{{
currentRead = 0;
reservedN = 0;
if(Nreads>0){
this->Nreads=Nreads;
knowNreads=true;
readIndex.reserve(Nreads+2);
}
readIndex.push_back(0);
if(Ntotal>0){
this->Ntotal=Ntotal;
knowNtotal=true;
reservedN = Ntotal+1;
trIds.reserve(reservedN);
probs.reserve(reservedN);
}
if(M>0){
this->M=M;
readsInIsoform.assign(M,-1);
}else{
readsInIsoform.clear();
this->M=0;
}
}//}}}
void TagAlignments::pushAlignment(long trId, double prob){//{{{
if(prob<=0)pushAlignmentL(trId, ns_misc::LOG_ZERO);
else pushAlignmentL(trId, log(prob));
}//}}}
void TagAlignments::pushAlignmentL(long trId, double lProb){//{{{
if(trId>=M){
M=trId+1;
readsInIsoform.resize(M,-1);
}
if(readsInIsoform[trId] == currentRead){
// The read has already one alignment to this transcript.
for(long i=readIndex[currentRead];i<(long)trIds.size();i++)
if(trIds[i] == trId){
probs[i] = ns_math::logAddExp(probs[i], lProb);
break;
}
}else{
if(! knowNtotal){
// the size of arrays is unknown try to reserve sensible amount of space if we know Nreads
if(knowNreads && reservedN && ((long)probs.size() == reservedN)){
// we reached the size of reserved space
double dens = (double)probs.size() / currentRead;
dens *= 1.05; //increase it by 5%
reservedN =(long)( reservedN + (dens) * (Nreads - currentRead + 1000.0) );
#ifdef MEM_USAGE
message("TagAlignments:\n size: %ld reserving: %ld capacity before: %ld\n",probs.size(),reservedN,probs.capacity());
#endif
trIds.reserve(reservedN);
probs.reserve(reservedN);
#ifdef MEM_USAGE
message(" capacity after: %ld\n",probs.capacity());
#endif
}else if(knowNreads && (! reservedN) && (currentRead == Nreads / 4 )){
// one quarter in, try to reserve sensible amount of space
double dens = (double)probs.size() / currentRead;
dens *= 1.05; //increase it by 5%
reservedN =(long)((dens) * (Nreads));
#ifdef MEM_USAGE
message("TagAlignments:\n size: %ld reserving: %ld capacity before: %ld\n",probs.size(),reservedN,probs.capacity());
#endif
trIds.reserve(reservedN);
probs.reserve(reservedN);
#ifdef MEM_USAGE
message(" capacity after: %ld\n",probs.capacity());
#endif
}
}
trIds.push_back(trId);
probs.push_back(lProb);
// Mark that transcript trId already has alignment from this read.
readsInIsoform[trId] = currentRead;
}
}//}}}
void TagAlignments::pushRead(){//{{{
// Check whether there were any valid alignments added for this read:
if(readIndex[currentRead] == (int_least32_t) probs.size()){
// If no new alignments, do nothing.
return;
}
// If there are alignments transform from log space if necessary and move to next read.
if(!storeLog){
double logSum = ns_math::logSumExp(probs, readIndex[currentRead], probs.size());
for(long i = readIndex[currentRead]; i<(long)probs.size(); i++)
probs[i] = exp(probs[i]-logSum);
}
// Move to the next read.
currentRead++;
readIndex.push_back(probs.size());
}//}}}
void TagAlignments::finalizeRead(long *M, long *Nreads, long *Ntotal){//{{{
*M = this->M = readsInIsoform.size();
*Nreads = this->Nreads = readIndex.size() - 1;
*Ntotal = this->Ntotal = probs.size();
#ifdef MEM_USAGE
message("TagAlignments: readIndex size: %ld capacity %ld\n",readIndex.size(),readIndex.capacity());
message("TagAlignments: probs size: %ld capacity %ld\n",probs.size(),probs.capacity());
#endif
}//}}}
int_least32_t TagAlignments::getTrId(long i) const {//{{{
if(i<Ntotal)return trIds[i];
return 0;
}//}}}
double TagAlignments::getProb(long i) const {//{{{
if(i<Ntotal)return probs[i];
return 0;
}//}}}
int_least32_t TagAlignments::getReadsI(long i) const {//{{{
if(i<=Nreads)return readIndex[i];
return 0;
}//}}}
|
