|
14 | 14 | #include <map> // for operator!=, __map_iterator, operator== |
15 | 15 | #include <vector> // for vector |
16 | 16 |
|
17 | | -static int decode_nt16_acgt(uint8_t nt16) { |
18 | | - if (nt16 == 1) return 0; // A |
19 | | - if (nt16 == 2) return 1; // C |
20 | | - if (nt16 == 4) return 2; // G |
21 | | - if (nt16 == 8) return 3; // T |
22 | | - return -1; |
23 | | -} |
24 | | - |
25 | | -double dust_score_nt16(const uint8_t *seq, int32_t l, int32_t window) { |
26 | | - const int32_t WLEN = 3; |
27 | | - const int32_t WTOT = 64; |
28 | | - const int32_t WMASK = WTOT - 1; |
29 | | - if (window < WLEN || window > 64) |
30 | | - return -1.0; |
31 | | - |
32 | | - int32_t wCount[WTOT]; |
33 | | - int32_t wSeq[64]; |
34 | | - memset(wCount, 0, sizeof(wCount)); |
35 | | - memset(wSeq, 0, sizeof(wSeq)); |
36 | | - |
37 | | - int64_t score = 0; |
38 | | - int64_t maxScore = 0; |
39 | | - int32_t t = 0; |
40 | | - int32_t n = -WLEN; |
41 | | - |
42 | | - for (int32_t i = 0; i < l; ++i) { |
43 | | - const int b = decode_nt16_acgt(bam_seqi(seq, i)); |
44 | | - if (b < 0) |
45 | | - continue; // ignore N/ambiguous |
46 | | - |
47 | | - t = ((t << 2) | b) & WMASK; |
48 | | - if (++n >= 0) { |
49 | | - const int32_t k = n % window; |
50 | | - if (n >= window) { |
51 | | - const int32_t x = wSeq[k]; |
52 | | - if (wCount[x] > 0) |
53 | | - score -= --wCount[x]; |
54 | | - score += wCount[t]++; |
55 | | - if (score > maxScore) |
56 | | - maxScore = score; |
57 | | - } else { |
58 | | - score += wCount[t]++; |
59 | | - } |
60 | | - wSeq[k] = t; |
61 | | - } |
62 | | - } |
63 | | - |
64 | | - if (n <= 0) |
65 | | - return 0.0; |
66 | | - if (n >= window) |
67 | | - return (200.0 * (double)maxScore) / (window * (window - 1)); |
68 | | - return (200.0 * (double)score) / (n * (n + 1)); |
69 | | -} |
70 | | - |
71 | 17 |
|
72 | 18 | int file_is_older(const char *file1, const char *file2) { |
73 | 19 | struct stat st1, st2; |
|
0 commit comments