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
|
#include "ReadAlign.h"
#include "ErrorWarning.h"
void ReadAlign::transformGenome()
{//convert to new genome
alignsGenOut.alBest = trBest; //it will be redefined below if alignment passes mapping filter
if (!mapGen.genomeOut.convYes || mapGen.pGe.transform.type==0 || nTr > P.outFilterMultimapNmax || nTr==0)
return; //no transformation necessary
uint32 nTr1=0;
for (uint32 iTr=0; iTr<nTr; iTr++) {//convert output transcripts into new genome
trMult[iTr]->haploType = ( trMult[iTr]->Chr >= mapGen.nChrReal/2 ? 1 : 2 );
*alignsGenOut.alMult[nTr1]=*trMult[iTr];//copy information before conversion
if (trMult[iTr]->transformGenome(*mapGen.genomeOut.g, *alignsGenOut.alMult[nTr1])) {
if (trBest == trMult[iTr])
alignsGenOut.alBest = alignsGenOut.alMult[nTr1]; //mark the best transcript
++nTr1;
};
};
if (mapGen.pGe.transform.type==2) {//diploid genome, remove identical transcripts
bool keepTr[nTr1];
for (uint32 ia1=0; ia1<nTr1; ia1++)
keepTr[ia1]=true;
for (uint32 ia1=0; ia1<nTr1; ia1++) {
if (!keepTr[ia1])
continue;
for (uint32 ia2=ia1+1; ia2<nTr1; ia2++) {//compare to all previous ones
if (!keepTr[ia1])
continue;
Transcript &a1 = *alignsGenOut.alMult[ia1];
Transcript &a2 = *alignsGenOut.alMult[ia2];
if ( a1.Chr==a2.Chr && a1.Str==a2.Str
&& a1.exons[0][EX_G]-a1.exons[0][EX_R] == a2.exons[0][EX_G]-a2.exons[0][EX_R]
&& a1.exons[a1.nExons][EX_G]+a1.exons[a1.nExons][EX_L]-a1.exons[a1.nExons][EX_R] == a2.exons[a2.nExons][EX_G]+a2.exons[a2.nExons][EX_L]-a2.exons[a2.nExons][EX_R]
) {//matching ends, remove one of the transcripts
alignsGenOut.alMult[ia1]->haploType = 0; //undefined haplotype
alignsGenOut.alMult[ia2]->haploType = 0;
if (a1.maxScore>a2.maxScore) {
keepTr[ia2]=false;
} else {
keepTr[ia1]=false;
};
};
};
};
alignsGenOut.alN=0;
for (uint32 ia1=0; ia1<nTr1; ia1++) {
if (keepTr[ia1]) {
*alignsGenOut.alMult[alignsGenOut.alN] = *alignsGenOut.alMult[ia1];
alignsGenOut.alN++;
};
};
} else {//haploid genome
// for (uint32 iTr=0; iTr<nTr1; iTr++)
// trMult[iTr] = alignsGenOut.alMult[iTr]; //point to new transcsript
alignsGenOut.alN=nTr1;
};
};
|
