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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
|
#include "IncludeDefine.h"
#include "Parameters.h"
#include "Transcript.h"
#include "extendAlign.h"
#include "binarySearch2.h"
// #include "stitchGapIndel.cpp"
intScore stitchAlignToTranscript(uint rAend, uint gAend, uint rBstart, uint gBstart, uint L, uint iFragB, uint sjAB, Parameters& P, char* R, Genome &mapGen, Transcript *trA, const uint outFilterMismatchNmaxTotal) {
//stitch together A and B, extend in the gap, returns max score
if (trA->nExons>=MAX_N_EXONS)
return -1000010;
char *G=mapGen.G;
int Score=0;
if (sjAB!=((uint) -1) && trA->exons[trA->nExons-1][EX_sjA]==sjAB \
&& trA->exons[trA->nExons-1][EX_iFrag]==iFragB && rBstart==rAend+1 && gAend+1<gBstart ) {//simple stitching if junction belongs to a database
if (mapGen.sjdbMotif[sjAB]==0 && (L<=mapGen.sjdbShiftRight[sjAB] || trA->exons[trA->nExons-1][EX_L]<=mapGen.sjdbShiftLeft[sjAB]) ) {
return -1000006; //too large repeats around non-canonical junction
};
trA->exons[trA->nExons][EX_L] = L; //new exon length
trA->exons[trA->nExons][EX_R] = rBstart; //new exon r-start
trA->exons[trA->nExons][EX_G] = gBstart; //new exon g-start
trA->canonSJ[trA->nExons-1]=mapGen.sjdbMotif[sjAB]; //mark sj-db
trA->shiftSJ[trA->nExons-1][0]=mapGen.sjdbShiftLeft[sjAB];
trA->shiftSJ[trA->nExons-1][1]=mapGen.sjdbShiftRight[sjAB];
trA->sjAnnot[trA->nExons-1]=1;
trA->sjStr[trA->nExons-1]=mapGen.sjdbStrand[sjAB];;
trA->nExons++;
trA->nMatch+=L;
for (uint ii=rBstart;ii<rBstart+L;ii++) Score+=scoreMatch; //add QS for mapped portions
Score+=P.pGe.sjdbScore;
} else {//general stitching
trA->sjAnnot[trA->nExons-1]=0;
trA->sjStr[trA->nExons-1]=0;
if (trA->exons[trA->nExons-1][EX_iFrag]==iFragB) {//stitch aligns on the same fragment
uint gBend=gBstart+L-1;
uint rBend=rBstart+L-1;
// {//debug
// if (sjAB!=((uint) -1) && trA->exons[trA->nExons-1][EX_sjA]!=((uint) -1) && rBend<=rAend) {//
// Score -= rAend-rBstart+1;
// gAend -= rAend-rBstart+1;
// rAend = rBstart-1;
// trA->exons[trA->nExons-1][EX_L] =rAend-trA->exons[trA->nExons-1][EX_R]+1;
// };
// };
//check if r-overlapping fully and exit
if (rBend<=rAend) return -1000001;
if (gBend<=gAend && trA->exons[trA->nExons-1][EX_iFrag]==iFragB) return -1000002;
//shift the B 5' if overlaps A 3'
if (rBstart<=rAend) {
gBstart+=rAend-rBstart+1;
rBstart=rAend+1;
L=rBend-rBstart+1;
};
for (uint ii=rBstart;ii<=rBend;ii++) Score+=scoreMatch; //add QS for mapped portions
int gGap=gBstart-gAend-1; //could be < 0 for insertions
int rGap=rBstart-rAend-1;//>0 always since we removed overlap
uint nMatch=L;
uint nMM=0;
uint Del=0, Ins=0;
uint nIns=0, nDel=0;
int jR=0; //junction location in R-space
int jCan=999; //canonical junction type
uint gBstart1=gBstart-rGap-1;//the last base of the intron if all read gap belongs to acceptor, i.e. jR=0
// check all the different combinations of gGap and rGap
if ( gGap==0 && rGap==0 ) {//just joined the pieces, w/o stiching or gaps
//do nothing for now
} else if ( gGap>0 && rGap>0 && rGap==gGap ) {//no gaps, just try to fill space
//simple stitching, assuming no insertion in the read
for (int ii=1;ii<=rGap;ii++) {
if (G[gAend+ii]<4 && R[rAend+ii]<4) {//only score genome bases that are not Ns
if ( R[rAend+ii]==G[gAend+ii] ) {
Score+=scoreMatch;
nMatch++;
} else {
Score-=scoreMatch;
nMM++;
};
};
};
} else if ( gGap>rGap ) {//genomic gap (Deletion)
nDel=1;
Del=gGap-rGap; //gGap>0 here
if (Del>P.alignIntronMax && P.alignIntronMax>0) {
return -1000003; //large gaps not allowed
};
int Score1=0;
int jR1=1; //junction location in R-space
do { // 1. move left, until the score for MM is less than canonical advantage
jR1--;
if ( R[rAend+jR1]!=G[gBstart1+jR1] && G[gBstart1+jR1]<4 && R[rAend+jR1]==G[gAend+jR1]) Score1 -= scoreMatch;
} while ( Score1+P.scoreStitchSJshift >= 0 && int(trA->exons[trA->nExons-1][EX_L]) + jR1 > 1);//>=P.alignSJoverhangMin); //also check that we are still within the exon
int maxScore2=-999999;
Score1=0;
int jPen=0;
do { // 2. scan to the right to find the best junction locus
// ?TODO? if genome base is N, how to score?
if ( R[rAend+jR1]==G[gAend+jR1] && R[rAend+jR1]!=G[gBstart1+jR1] ) Score1+=scoreMatch;
if ( R[rAend+jR1]!=G[gAend+jR1] && R[rAend+jR1]==G[gBstart1+jR1] ) Score1-=scoreMatch;
int jCan1=-1; //this marks Deletion
int jPen1=0;
int Score2=Score1;
if (Del>=P.alignIntronMin) {//only check intron motif for large gaps= non-Dels
//check if the intron is canonical, or semi-canonical
if ( G[gAend+jR1+1]==2 && G[gAend+jR1+2]==3 && G[gBstart1+jR1-1]==0 && G[gBstart1+jR1]==2 ) {//GTAG
jCan1=1;
} else if ( G[gAend+jR1+1]==1 && G[gAend+jR1+2]==3 && G[gBstart1+jR1-1]==0 && G[gBstart1+jR1]==1 ) {//CTAC
jCan1=2;
} else if ( G[gAend+jR1+1]==2 && G[gAend+jR1+2]==1 && G[gBstart1+jR1-1]==0 && G[gBstart1+jR1]==2 ) {//GCAG
jCan1=3;
jPen1=P.scoreGapGCAG;
} else if ( G[gAend+jR1+1]==1 && G[gAend+jR1+2]==3 && G[gBstart1+jR1-1]==2 && G[gBstart1+jR1]==1 ) {//CTGC
jCan1=4;
jPen1=P.scoreGapGCAG;
} else if ( G[gAend+jR1+1]==0 && G[gAend+jR1+2]==3 && G[gBstart1+jR1-1]==0 && G[gBstart1+jR1]==1 ) {//ATAC
jCan1=5;
jPen1=P.scoreGapATAC;
} else if ( G[gAend+jR1+1]==2 && G[gAend+jR1+2]==3 && G[gBstart1+jR1-1]==0 && G[gBstart1+jR1]==3 ) {//GTAT
jCan1=6;
jPen1=P.scoreGapATAC;
} else {
jCan1=0;
jPen1=P.scoreGapNoncan;
};
Score2 += jPen1;
};
if (maxScore2 < Score2 ) {//check if the score is the highest. TODO: record the next highest score
maxScore2=Score2;
jR=jR1; //this is the last base of donor
jCan=jCan1;
jPen=jPen1;
};
jR1++;
} while ( jR1 < int(rBend) - int(rAend) );// - int(P.alignSJoverhangMin) );//TODO: do not need to search the full B-transcript, can stop as soon as Score goes down by more than
//repeat length: go back and forth around jR to find repeat length
uint jjL=0,jjR=0;
while ( gAend+jR>=jjL && G[gAend-jjL+jR]==G[gBstart1-jjL+jR] && G[gAend-jjL+jR]<4 && jjL<=MAX_SJ_REPEAT_SEARCH) {//go back
jjL++;
};
while ( gAend+jjR+jR+1<mapGen.nGenome && G[gAend+jjR+jR+1]==G[gBstart1+jjR+jR+1] && G[gAend+jjR+jR+1]<4 && jjR<=MAX_SJ_REPEAT_SEARCH) {//go forward
jjR++;
};
if (jCan<=0) {//flush deletions and non-canonical junction to the left
jR-=jjL;
if (int(trA->exons[trA->nExons-1][EX_L])+jR<1) return -1000005;
jjR+=jjL;
jjL=0;
};
//TODO check here if the internal exon length < minDa, if so exit w/o stitiching
for (int ii=min(1,jR+1);ii<=max(rGap,jR);ii++) {//score donor and acceptor
uint g1=(ii<=jR) ? (gAend+ii):(gBstart1+ii);
if (G[g1]<4 && R[rAend+ii]<4) {//only penalize non-N bases
if ( R[rAend+ii]==G[g1] ) {
if (ii>=1 && ii <=rGap) {//only add +score and matches within the gap
Score+=scoreMatch;
nMatch++;
};
} else {//add -score and MM for all bases
Score-=scoreMatch;
nMM++;
if (ii<1 || ii>rGap) {//subtract previuosly presumed matches
Score-=scoreMatch;
nMatch--;
// if (ii<=jR) nMM--;
};
};
};
};
//score the gap
if (mapGen.sjdbN>0) {//check if the junction is annotated
uint jS=gAend+jR+1, jE=gBstart1+jR;//intron start/end
int sjdbInd=binarySearch2(jS,jE,mapGen.sjdbStart,mapGen.sjdbEnd,mapGen.sjdbN);
if (sjdbInd<0) {
if (Del>=P.alignIntronMin) {
Score += P.scoreGap + jPen; //genome gap penalty + non-canonical penalty
} else {//deletion
Score += Del*P.scoreDelBase + P.scoreDelOpen;
jCan=-1;
trA->sjAnnot[trA->nExons-1]=0;
// jjR-=jjL;
// jR-=jjL;
// jjL=0;
// trA->shiftSJ[trA->nExons-1][0]=0;
// trA->shiftSJ[trA->nExons-1][1]=jjR;
};
} else {//annotated
jCan=mapGen.sjdbMotif[sjdbInd];
if (mapGen.sjdbMotif[sjdbInd]==0) {//shift to match annotations
if (L<=mapGen.sjdbShiftLeft[sjdbInd] || trA->exons[trA->nExons-1][EX_L]<=mapGen.sjdbShiftLeft[sjdbInd]) {
//this is not needed, the check below rAend+jR >= rBend is enough
//it changes the results slightly, so will keep for now. Will remove in STAR3
return -1000006;
};
jR += (int) mapGen.sjdbShiftLeft[sjdbInd];
if ( rAend+jR >= rBend ) {//jR is shifted past rBend
return -1000006;
};
jjL=mapGen.sjdbShiftLeft[sjdbInd];
jjR=mapGen.sjdbShiftRight[sjdbInd];
};
trA->sjAnnot[trA->nExons-1]=1;
trA->sjStr[trA->nExons-1]=mapGen.sjdbStrand[sjdbInd];
Score += P.pGe.sjdbScore;
};
} else {//no annotation
if (Del>=P.alignIntronMin) {//junction, not short deletion
Score += P.scoreGap + jPen;
} else {
Score += Del*P.scoreDelBase + P.scoreDelOpen;
jCan=-1;
trA->sjAnnot[trA->nExons-1]=0;
};
};
trA->shiftSJ[trA->nExons-1][0]=jjL;
trA->shiftSJ[trA->nExons-1][1]=jjR;
trA->canonSJ[trA->nExons-1]=jCan;
if (trA->sjAnnot[trA->nExons-1]==0) {//strand for unannotated junctions
if (jCan>0) {
trA->sjStr[trA->nExons-1]=2-jCan%2; //1=+,2=-
} else {
trA->sjStr[trA->nExons-1]=0;
};
};
} else if ( rGap>gGap ) {//insertion: if also gGap>0, need to stitch
Ins=rGap-gGap;
nIns=1;
if (gGap==0) {//simple insertion, no need to stitch
jR=0;
} else if (gGap<0) {//overlapping seeds: reduce the score
jR=0;
for (int ii=0; ii<-gGap; ii++) {
Score -= scoreMatch;
};
} else {//stitch: define the exon boundary jR
int Score1=0; int maxScore1=0;
for (int jR1=1;jR1<=gGap;jR1++) {//scan to the right to find the best score
if (G[gAend+jR1]<4) {//only penalize goog genome bases
Score1+=( R[rAend+jR1]==G[gAend+jR1] ) ? scoreMatch:-scoreMatch;
Score1+=( R[rAend+Ins+jR1]==G[gAend+jR1] ) ? -scoreMatch:+scoreMatch;
};
if (Score1>maxScore1 || (Score1==maxScore1 && P.alignInsertionFlush.flushRight)) {//equal sign (>=) flushes insertions to the right
maxScore1=Score1;
jR=jR1;
};
};
for (int ii=1;ii<=gGap;ii++) {//score donor and acceptor
uint r1=rAend+ii+(ii<=jR ? 0:Ins);
if (G[gAend+ii]<4 && R[r1]<4) {
if ( R[r1]==G[gAend+ii] ) {
Score+=scoreMatch;
nMatch++;
} else {//add -score and MM for all bases
Score-=scoreMatch;
nMM++;
};
};
};
};
if (P.alignInsertionFlush.flushRight) {
for (; jR<(int)rBend-(int)rAend-(int)Ins; jR++ ){//flush the indel to the right as much as possible
if (R[rAend+jR+1]!=G[gAend+jR+1] || G[gAend+jR+1]==4) {
break;
};
};
if (jR==(int)rBend-(int)rAend-(int)Ins) {//nothing left of the B-piece
return -1000009;
};
};
Score += Ins*P.scoreInsBase + P.scoreInsOpen;
jCan=-2; //marks insertion though it's not used below
}; //different types of gaps selection
#ifdef COMPILE_FOR_LONG_READS
if ( (trA->nMM + nMM)<=outFilterMismatchNmaxTotal )
// if ( Score>0 && nMM<=200 )
#else
if ( (trA->nMM + nMM)<=outFilterMismatchNmaxTotal \
&& ( jCan<0 || (jCan<7 && nMM<= (uint) P.alignSJstitchMismatchNmax[(jCan+1)/2]) ) )
#endif
{//stitching worked only if there no mis-matches for non-GT/AG junctions
trA->nMM += nMM;
trA->nMatch += nMatch;
if (Del>=P.alignIntronMin) {
trA->nGap += nDel;
trA->lGap += Del;
} else {
trA->nDel += nDel;
trA->lDel += Del;
};
//modify exons
if (Del==0 && Ins==0) {//no gap => no new exon, extend the boundary of the previous exon
trA->exons[trA->nExons-1][EX_L] += rBend-rAend;
} else if (Del>0) { //deletion:ca only have Del> or Ins>0
trA->exons[trA->nExons-1][EX_L] += jR; //correct the previous exon boundary
trA->exons[trA->nExons][EX_L] = rBend-rAend-jR; //new exon length
trA->exons[trA->nExons][EX_R] = rAend+jR+1; //new exon r-start
trA->exons[trA->nExons][EX_G] = gBstart1+jR+1; //new exon g-start
trA->nExons++;
} else if (Ins>0) { //Ins>0;
trA->nIns += nIns;
trA->lIns += Ins;
trA->exons[trA->nExons-1][EX_L] += jR; //correct the previous exon boundary
trA->exons[trA->nExons][EX_L] = rBend-rAend-jR-Ins; //new exon length
trA->exons[trA->nExons][EX_R] = rAend+jR+Ins+1; //new exon r-start
trA->exons[trA->nExons][EX_G] = gAend+1+jR; //new exon g-start
trA->canonSJ[trA->nExons-1]=-2; //mark insertion
trA->sjAnnot[trA->nExons-1]=0;
trA->nExons++;
};
} else {//stitching was not accepted
return -1000007;
};
} else if (gBstart+trA->exons[0][EX_R]+P.alignEndsProtrude.nBasesMax >= trA->exons[0][EX_G] || trA->exons[0][EX_G] < trA->exons[0][EX_R]){//if (iFragA==iFragB) stitch aligns from different fragments
//CHECK: this second confdition does not make sense
if (P.alignMatesGapMax>0 && gBstart > trA->exons[trA->nExons-1][EX_G] + trA->exons[trA->nExons-1][EX_L] + P.alignMatesGapMax) {
return -1000004; //gap between mates too large
};
//extend the fragments inside
//note, that this always works, i.e. Score>0
for (uint ii=rBstart;ii<rBstart+L;ii++) Score+=scoreMatch; //add QS for mapped portions
Transcript trExtend;
//TODO: compare extensions to the left and right, pick the best one to be performed first
//otherwise if a large nMM is reached in the 2st extension, it will prevent the 2nd extension
//use the following example:
//>1
//TTCTGTGTCTCCCCCTCCCCCACTGGCTACATGGAGACAGGGGGGGGGGGCCGGGCGGTTCCCGGGCAGAAAAAAA
//>1
//AATATTTGGAACACTTATGTGAAAAATGATTTGTTTTTCTGAAATTTACGTTTCTCTCTGAGTCCTGTAACTGTCC
trExtend.reset();
if ( extendAlign(R, G, rAend+1, gAend+1, 1, 1, DEF_readSeqLengthMax, trA->nMatch, trA->nMM, outFilterMismatchNmaxTotal, P.outFilterMismatchNoverLmax, \
P.alignEndsType.ext[trA->exons[trA->nExons-1][EX_iFrag]][1], &trExtend) ) {
trA->add(&trExtend);
Score += trExtend.maxScore;
trA->exons[trA->nExons-1][EX_L] += trExtend.extendL;
};// if extendAlign for read A
trA->exons[trA->nExons][EX_R] = rBstart;
trA->exons[trA->nExons][EX_G] = gBstart;
trA->exons[trA->nExons][EX_L] = L;
trA->nMatch += L;
trExtend.reset();
//if end extension needs to be forced, use large length. Otherwise, only extend until the beginning of the transcript
uint extlen=P.alignEndsType.ext[iFragB][1] ? DEF_readSeqLengthMax : gBstart-trA->exons[0][EX_G]+trA->exons[0][EX_R];
if ( extendAlign(R, G, rBstart-1, gBstart-1, -1, -1, extlen, trA->nMatch, trA->nMM, outFilterMismatchNmaxTotal, P.outFilterMismatchNoverLmax, \
P.alignEndsType.ext[iFragB][1], &trExtend) ) {
trA->add(&trExtend);
Score += trExtend.maxScore;
trA->exons[trA->nExons][EX_R] -= trExtend.extendL;
trA->exons[trA->nExons][EX_G] -= trExtend.extendL;
trA->exons[trA->nExons][EX_L] += trExtend.extendL;
}; //if extendAlign B
trA->canonSJ[trA->nExons-1]=-3; //mark different fragments junction
trA->sjAnnot[trA->nExons-1]=0;
trA->nExons++;
} else {//no stitching possible
return -1000008;
};
};
trA->exons[trA->nExons-1][EX_iFrag]=iFragB; //the new exon belongs to fragment iFragB
trA->exons[trA->nExons-1][EX_sjA]=sjAB;
return Score;
};
|
