-
Notifications
You must be signed in to change notification settings - Fork 0
/
read_quality.h
104 lines (73 loc) · 3.66 KB
/
read_quality.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
104
#ifndef _READ_QUALITY_H_
#define _READ_QUALITY_H_
#include "load_data.h"
#include "seed.h"
/* Read quality */
#define MAX_READ_QUALITY 93
#define READ_QUALITY_LEVELS 4
#define MAX_READ_QUALITY_LEVEL 3
#define READ_QUALITY_LEVEL0_UB 2
#define READ_QUALITY_LEVEL1_UB 7
#define READ_QUALITY_LEVEL2_UB 10
#define READ_QUALITY_LEVEL3_UB MAX_READ_QUALITY
extern const short READ_QUALITY_LEVEL_UPPER_BOUNDS[READ_QUALITY_LEVELS];
extern int read_quality_min_symbol_code;
/*
* Qualities
*/
typedef unsigned char QUAL_TYPE;
#define QUAL_SIZE_BITS (sizeof(QUAL_TYPE))
#define QUAL_SIZE_USED_BITS 2
#define QMASK (0x3)
#define SHIFTED_QMASK (0xC0)
#define QUAL_SEQUENCE_BLOCK_BITS (QUAL_SIZE_BITS << 3)
/* sizeof(QUAL_TYPE) * 8) */
#define QUALS_PER_SEQUENCE_BLOCK (QUAL_SIZE_BITS << 2)
/*(sizeof(QUAL_TYPE) * 8 / QUAL_SIZE_USED_BITS)*/
#define QMASK_SHIFT ((QUAL_SIZE_BITS) - (QUAL_SIZE_USED_BITS))
/* Actual allocated length for keeping len codes in a compressed sequence*/
#define COMPRESSED_Q_LEN(len) (((len) >> 2) + (((len) & 3) ? 1 : 0))
/* For QUALS_PER_SEQUENCE_BLOCK = 4:
* ceiling(len / (float)QUALS_PER_SEQUENCE_BLOCK) = len / QUALS_PER_SEQUENCE_BLOCK + ((len % QUALS_PER_SEQUENCE_BLOCK) > 0 ? 1 : 0)
*/
#define COMPRESSED_Q_IDX(i) ((i) >> 2)
/* For QUALS_PER_SEQUENCE_BLOCK = 4: len /QUALS_PER_SEQUENCE_BLOCK) */
#define COMPRESSED_Q_OFFSET(i) ((i) & 3)
/* For QUALS_PER_SEQUENCE_BLOCK = 4: len % QUALS_PER_SEQUENCE_BLOCK) */
#define TO_NTH_QUAL(sequence, n, val) (sequence[COMPRESSED_Q_IDX(n)] |= ((val) & QMASK) << (COMPRESSED_Q_OFFSET(n) << 1));
/* sequence[n/QUALS_PER_SEQUENCE_BLOCK] |= (((val & QMASK) << (QMASK_SHIFT)) >> ((n % QUALS_PER_SEQUENCE_BLOCK)*QUAL_SIZE_USED_BITS)); */
#define NTH_QUAL(sequence, n) ((sequence[COMPRESSED_Q_IDX(n)] & (QMASK << (COMPRESSED_Q_OFFSET(n) << 1))) >> (COMPRESSED_Q_OFFSET(n) << 1))
/* For QUALS_PER_SEQUENCE_BLOCK = 4, QMASK = 3, QMASK_SHIFT = 6, QUAL_SIZE_USED_BITS = 2:
* ((((QMASK << (QMASK_SHIFT)) >> ((n % QUALS_PER_SEQUENCE_BLOCK)*QUAL_SIZE_USED_BITS)) & sequence[n/QUALS_PER_SEQUENCE_BLOCK]) >> (QMASK_SHIFT - ((n%QUALS_PER_SEQUENCE_BLOCK)*QUAL_SIZE_USED_BITS)))
*/
/*
* Seed quality mask
*/
#define QSEED_MASK (0x1f) /* 0x1f = 31 to detect value that are multiple of 32 */
#define QSEED_MASK_SHIFT (0x5) /* 2^5 = 32 to get true coordinates */
#define COMPRESSED_QSEED_LEN(len) (((len) >> QSEED_MASK_SHIFT) + (((len) & QSEED_MASK) ? 1 : 0))
#define COMPRESSED_QSEED_IDX(i) ((i) >> QSEED_MASK_SHIFT)
#define COMPRESSED_QSEED_OFFSET(i) ((i) & QSEED_MASK)
#define SET_NTH_QSEED(sequence, n) (sequence[COMPRESSED_QSEED_IDX(n)] |= (1) << (COMPRESSED_QSEED_OFFSET(n)))
#define GET_NTH_QSEED(sequence, n) ((sequence[COMPRESSED_QSEED_IDX(n)] & (QSEED_MASK << COMPRESSED_QSEED_OFFSET(n))) >> (COMPRESSED_QSEED_OFFSET(n)))
/* Read quality handling */
/*(READ_QUALITY_LEVELS - 1) * MIN (Q, 40) / 40 => index of the score */
/* @Deprecated */
/* #define QUALITY_LEVEL(q) ((int) ( ( (READ_QUALITY_LEVELS - 1) * MIN((q), MAX_READ_QUALITY)) / MAX_READ_QUALITY + .5))*/
int get_quality_level(short quality);
#define QUALITY_LEVEL(q) (get_quality_level(q))
int fastq_parse_quality_from_char(char c);
#define FASTQ_PARSE_QUALITY_FROM_CHAR(c) (fastq_parse_quality_from_char(c))
/**
* Allocate and compress the list of "short" qualities
*/
QUAL_TYPE* compress_quality_sequence(const short* qual, int len, int shift);
/**
* Obtain the reverse quality sequence (used for RC alignments)
*/
#ifdef NUCLEOTIDES
void quality__reverse_compressed(const QUAL_TYPE* sequence, QUAL_TYPE* dest, int len);
#else
void quality__reverse_compressed(const QUAL_TYPE first, const QUAL_TYPE* sequence, QUAL_TYPE* dest, int len);
#endif
#endif /* _READ_QUALITY_H_ */