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
|
#include <pbdata/FASTQReader.hpp>
#include <climits>
#include <cmath>
#include <cstdio>
FASTQReader::FASTQReader() : FASTAReader() { endOfReadDelim = '\n'; }
int FASTQReader::GetNext(FASTASequence &seq) { return ((FASTAReader *)this)->GetNext(seq); }
unsigned char FASTQReader::phredQVtoPacbioQV(unsigned char phredQV)
{
int qual = floor(100.0 * log10(pow(10.0, phredQV / 10.0) - 1.0) + 0.5);
qual = qual > 250 ? 250 : qual;
qual = qual < 1 ? 1 : qual;
return (unsigned char)qual;
}
int FASTQReader::GetNext(FASTQSequence &seq)
{
seq.Free(); // Free seq before being reused.
char c;
while (curPos < fileSize and
((c = filePtr[curPos]) == ' ' or c == '\t' or c == '\n' or c == '\r')) {
curPos++;
}
if (curPos >= fileSize) {
return false;
}
GenomeLength p = curPos;
AdvanceToTitleStart(p, '@');
CheckValidTitleStart(p, '@');
ReadTitle(p, seq);
// Title ends on '\n', consume that;
p++;
GenomeLength p2;
p2 = p;
while (p2 < fileSize and filePtr[p2] != '\n') {
p2++;
}
if (p2 - p > UINT_MAX) {
std::cout
<< "ERROR! Reading sequences stored in more than 4Gbytes of space is not supported."
<< std::endl;
std::exit(EXIT_FAILURE);
}
seq.length = p2 - p;
GenomeLength seqPos;
if (seq.length > 0) {
seq.seq = ProtectedNew<Nucleotide>(seq.length);
p2 = p;
seqPos = 0;
while (p2 < fileSize and filePtr[p2] != '\n') {
seq.seq[seqPos] = filePtr[p2];
p2++;
seqPos++;
}
} else {
seq.seq = 0;
}
p = p2;
AdvanceToTitleStart(p, '+');
CheckValidTitleStart(p, '+');
while (p < fileSize and filePtr[p] != '\n') {
p++;
}
p++; // skip '\n'
p2 = p;
while (p2 < fileSize and filePtr[p2] != '\n') {
p2++;
}
seq.length = p2 - p;
if (seq.length > 0) {
seq.qual.Allocate(seq.length);
p2 = p;
GenomeLength seqPos = 0;
while (p2 < fileSize and filePtr[p2] != '\n') {
seq.qual[seqPos] = filePtr[p2] - FASTQSequence::charToQuality;
p2++;
seqPos++;
}
} else {
seq.qual.data = NULL;
}
curPos = p2;
seq.deleteOnExit = true;
return 1;
}
int FASTQReader::Advance(int nSteps)
{
// An advance of a FASTQ file is simply twice the number of
// advances of FASTA, since each nucleotide sequence has a quality
// sequence.
return ((FASTAReader *)this)->Advance(nSteps * 2);
}
|
