-
Notifications
You must be signed in to change notification settings - Fork 0
/
transversions.h
103 lines (84 loc) · 2.25 KB
/
transversions.h
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
double transversions_mutations(const char *subject, size_t length,
const char *query)
{
size_t mutations = 0;
const char *sbj = subject, *qry = query;
for (; *sbj; sbj++, qry++) {
if (*sbj != *qry) {
mutations++;
}
}
return static_cast<double>(mutations) / length;
}
double transversions_mutations_length(const char *subject, size_t length,
const char *query)
{
size_t mutations = 0;
for (size_t k = 0; k < length; k++) {
if (subject[k] != query[k]) {
mutations++;
}
}
return static_cast<double>(mutations) / length;
}
double transversions_transversions(const char *subject, size_t length,
const char *query)
{
static const auto is_purine = [](char c) { return c == 'A' || c == 'G'; };
size_t transversions = 0;
const char *sbj = subject, *qry = query;
/*for (; *sbj; sbj++, qry++) {
if ((is_purine(*sbj) && !is_purine(*qry)) ||
(!is_purine(*sbj) && is_purine(*qry))) {
transversions++;
}
}*/
for (size_t k = 0; k < length; k++) {
if ((is_purine(subject[k]) && !is_purine(query[k])) ||
(!is_purine(subject[k]) && is_purine(query[k]))) {
transversions++;
}
}
return static_cast<double>(transversions) / length;
}
double transversions_twiddle_length(const char *subject, size_t length,
const char *query)
{
size_t transversions = 0;
for (size_t k = 0; k < length; k++) {
if (((subject[k] + 1) ^ (query[k] + 1)) & 4) {
transversions++;
}
}
return static_cast<double>(transversions) / length;
}
double transversions_builtin(const char *subject, size_t length,
const char *query)
{
size_t transversions = 0;
for (size_t k = 0; k < length; k++) {
auto c1 = __builtin_popcount(subject[k]);
auto c2 = __builtin_popcount(query[k]);
// if ((c1 == 3 && c2 != 3) || (c1 != 3 && c2 == 3)) {
if ((c1 == 3) ^ (c2 == 3)) {
transversions++;
}
}
return static_cast<double>(transversions) / length;
}
double transversions_table(const char *subject, size_t length,
const char *query)
{
static char table[127];
table['A'] = 1;
table['G'] = 1;
table['C'] = 0;
table['T'] = 0;
size_t transversions = 0;
for (size_t k = 0; k < length; k++) {
if (table[subject[k]] ^ table[query[k]]) {
transversions++;
}
}
return static_cast<double>(transversions) / length;
}