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
|
#include "Transcript.h"
#include "serviceFuns.cpp"
int Transcript::variationAdjust(const Genome &mapGen, char *R)
{
Variation &Var=*mapGen.Var;
if (!Var.yes)
{//no variation
return 0;
};
int dScore=0;//change in the score
uint nMM1=0;
//for each block, check whether it overlaps one or more SNPs
for (uint ie=0; ie<nExons; ie++)
{
//binary search to find nearest SNP
int32 isnp=binarySearch1b <uint> (exons[ie][EX_G], Var.snp.loci, Var.snp.N);
if (isnp>=0)
{
while ((uint)isnp<Var.snp.N && exons[ie][EX_G]+exons[ie][EX_L]>Var.snp.loci[isnp])
{//these SNPs overlap the block
varInd.push_back(isnp); //record snp index
varGenCoord.push_back(Var.snp.loci[isnp]-mapGen.chrStart[Chr]);
varReadCoord.push_back(exons[ie][EX_R]+Var.snp.loci[isnp]-exons[ie][EX_G]);
char ntR=R[varReadCoord.back()];//nt of the read in the SNP position, already trnasformed to + genome strand
uint8 igt;
if (ntR>3) {
igt=4;
} else {
for (igt=1; igt<3; igt++) {//1st or 2nd allele, =3 of none
if (Var.snp.nt[isnp][igt]==ntR) {
break;
};
};
};
//if (ntR == Var.snp.nt[isnp][0])
//{//mark snp that agrees with the reference
// igt*=10;
//};
varAllele.push_back(igt);
if (igt<3 && ntR != Var.snp.nt[isnp][0])
{//non-reference allele, correct nMM and score
++nMM1;
};
++isnp;
};
};
};
#define VAR_noScoreCorrection
#ifndef VAR_noScoreCorrection
if (nMM1>0)
{//one or more mismtaches need to be corrected
uint nMMold=nMM;
alignScore(Read1, G, P);
nMM-=nMM1;
nMatch+=nMM1;
dScore=2*(nMMold-nMM);//score only changes if the number of mismatches is reduced after SNP adjustment
};
#else
//#warning VAR_noScoreCorrection set: no variation score correction
#endif
return dScore;
};
|
