-
Notifications
You must be signed in to change notification settings - Fork 0
/
stdaln_interface.cpp
95 lines (82 loc) · 2.31 KB
/
stdaln_interface.cpp
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
/*
* =====================================================================================
*
* Filename: stdaln_interface.cpp
*
* Description: interface to Stdaln lib
*
* Version: 1.0
* Created: 2006年11月21日 19时27分54秒 CST
* Revision: none
* Compiler: gcc
*
* Author: (),
* Company:
*
* =====================================================================================
*/
#include "stdaln_interface.h"
#include "struct.h"
#include <iostream>
using namespace std;
int gap_open = 3; //add by xz
int gap_ext = 1;
int gap_end = 1;
int score_match = 2;
int score_mismatch = -2;
int global_align(const char* str1, int begin1, int end1, const char* str2, int begin2, int end2, int* max_gap)
{
AlnAln *aln_global;
aln_global = aln_stdaln_aux(str1 + begin1, str2 + begin2, &aln_param_nt2nt, 1, end1 - begin1 + 1, end2 - begin2 + 1);
int s = aln_global->score;
int g_count = 0;
int m_count = 0;
*max_gap = 0;
for (int i = 0;i < aln_global->path_len;++i)
{
if (aln_global->outm[i] == ' ') //miss match
{
g_count++;
}
else
{
if (g_count >= 2) //2 or more mismatch is gap, 1 is SNP//???here is unsure,may be we should check '-' in out1&2
{
if (CHECK_PAIRS)
cout << " gap between pairs:" << g_count << endl;
*max_gap += g_count;
}
g_count = 0;
++m_count;
}
}
aln_free_AlnAln(aln_global);
//return s;
return m_count;
}
/* initialize score matrix, only for nt2nt */
void set_score_matrix(int gap_open, int gap_ext, int gap_end, int match_score, int mismatch_score)
{
int i, j, *p;
int row = aln_param_nt2nt.row;
for (i = 0; i < row; ++i)
{
for (j = 0; j < row; ++j)
{
p = aln_param_nt2nt.matrix + i * row + j;
if (*p == 2)
{ /* match */
*p = score_match; //by xz
}
else if (*p == -2)
{
*p = score_mismatch; //by xz, else *p = -2;
}
//printf("%4d",*p);
}
//printf("\n");
}
aln_param_nt2nt.gap_open = gap_open;
aln_param_nt2nt.gap_ext = gap_ext;
aln_param_nt2nt.gap_end = gap_end;
}