|
0
|
1 // ***************************************************************************
|
|
|
2 // FastaIndex.cpp (c) 2010 Erik Garrison <erik.garrison@bc.edu>
|
|
|
3 // Marth Lab, Department of Biology, Boston College
|
|
|
4 // All rights reserved.
|
|
|
5 // ---------------------------------------------------------------------------
|
|
|
6 // Last modified: 9 February 2010 (EG)
|
|
|
7 // ---------------------------------------------------------------------------
|
|
|
8
|
|
|
9 #include "Fasta.h"
|
|
|
10
|
|
|
11 FastaIndexEntry::FastaIndexEntry(string name, int length, long long offset, int line_blen, int line_len)
|
|
|
12 : name(name)
|
|
|
13 , length(length)
|
|
|
14 , offset(offset)
|
|
|
15 , line_blen(line_blen)
|
|
|
16 , line_len(line_len)
|
|
|
17 {}
|
|
|
18
|
|
|
19 FastaIndexEntry::FastaIndexEntry(void) // empty constructor
|
|
|
20 { clear(); }
|
|
|
21
|
|
|
22 FastaIndexEntry::~FastaIndexEntry(void)
|
|
|
23 {}
|
|
|
24
|
|
|
25 void FastaIndexEntry::clear(void)
|
|
|
26 {
|
|
|
27 name = "";
|
|
|
28 length = NULL;
|
|
|
29 offset = -1; // no real offset will ever be below 0, so this allows us to
|
|
|
30 // check if we have already recorded a real offset
|
|
|
31 line_blen = NULL;
|
|
|
32 line_len = NULL;
|
|
|
33 }
|
|
|
34
|
|
|
35 ostream& operator<<(ostream& output, const FastaIndexEntry& e) {
|
|
|
36 // just write the first component of the name, for compliance with other tools
|
|
|
37 output << split(e.name, ' ').at(0) << "\t" << e.length << "\t" << e.offset << "\t" <<
|
|
|
38 e.line_blen << "\t" << e.line_len;
|
|
|
39 return output; // for multiple << operators.
|
|
|
40 }
|
|
|
41
|
|
|
42 FastaIndex::FastaIndex(void)
|
|
|
43 {}
|
|
|
44
|
|
|
45 void FastaIndex::readIndexFile(string fname) {
|
|
|
46 string line;
|
|
|
47 long long linenum = 0;
|
|
|
48 indexFile.open(fname.c_str(), ifstream::in);
|
|
|
49 if (indexFile.is_open()) {
|
|
|
50 while (getline (indexFile, line)) {
|
|
|
51 ++linenum;
|
|
|
52 // the fai format defined in samtools is tab-delimited, every line being:
|
|
|
53 // fai->name[i], (int)x.len, (long long)x.offset, (int)x.line_blen, (int)x.line_len
|
|
|
54 vector<string> fields = split(line, '\t');
|
|
|
55 if (fields.size() == 5) { // if we don't get enough fields then there is a problem with the file
|
|
|
56 // note that fields[0] is the sequence name
|
|
|
57 char* end;
|
|
|
58 string name = split(fields[0], " \t").at(0); // key by first token of name
|
|
|
59 sequenceNames.push_back(name);
|
|
|
60 this->insert(make_pair(name, FastaIndexEntry(fields[0], atoi(fields[1].c_str()),
|
|
|
61 strtoll(fields[2].c_str(), &end, 10),
|
|
|
62 atoi(fields[3].c_str()),
|
|
|
63 atoi(fields[4].c_str()))));
|
|
|
64 } else {
|
|
|
65 cerr << "Warning: malformed fasta index file " << fname <<
|
|
|
66 "does not have enough fields @ line " << linenum << endl;
|
|
|
67 cerr << line << endl;
|
|
|
68 exit(1);
|
|
|
69 }
|
|
|
70 }
|
|
|
71 } else {
|
|
|
72 cerr << "could not open index file " << fname << endl;
|
|
|
73 exit(1);
|
|
|
74 }
|
|
|
75 }
|
|
|
76
|
|
|
77 // for consistency this should be a class method
|
|
|
78 bool fastaIndexEntryCompare ( FastaIndexEntry a, FastaIndexEntry b) { return (a.offset<b.offset); }
|
|
|
79
|
|
|
80 ostream& operator<<(ostream& output, FastaIndex& fastaIndex) {
|
|
|
81 vector<FastaIndexEntry> sortedIndex;
|
|
|
82 for(vector<string>::const_iterator it = fastaIndex.sequenceNames.begin(); it != fastaIndex.sequenceNames.end(); ++it)
|
|
|
83 {
|
|
|
84 sortedIndex.push_back(fastaIndex[*it]);
|
|
|
85 }
|
|
|
86 sort(sortedIndex.begin(), sortedIndex.end(), fastaIndexEntryCompare);
|
|
|
87 for( vector<FastaIndexEntry>::iterator fit = sortedIndex.begin(); fit != sortedIndex.end(); ++fit) {
|
|
|
88 output << *fit << endl;
|
|
|
89 }
|
|
|
90 return output;
|
|
|
91 }
|
|
|
92
|
|
|
93 void FastaIndex::indexReference(string refname) {
|
|
|
94 // overview:
|
|
|
95 // for line in the reference fasta file
|
|
|
96 // track byte offset from the start of the file
|
|
|
97 // if line is a fasta header, take the name and dump the last sequnece to the index
|
|
|
98 // if line is a sequence, add it to the current sequence
|
|
|
99 //cerr << "indexing fasta reference " << refname << endl;
|
|
|
100 string line;
|
|
|
101 FastaIndexEntry entry; // an entry buffer used in processing
|
|
|
102 entry.clear();
|
|
|
103 int line_length = 0;
|
|
|
104 long long offset = 0; // byte offset from start of file
|
|
|
105 long long line_number = 0; // current line number
|
|
|
106 bool mismatchedLineLengths = false; // flag to indicate if our line length changes mid-file
|
|
|
107 // this will be used to raise an error
|
|
|
108 // if we have a line length change at
|
|
|
109 // any line other than the last line in
|
|
|
110 // the sequence
|
|
|
111 bool emptyLine = false; // flag to catch empty lines, which we allow for
|
|
|
112 // index generation only on the last line of the sequence
|
|
|
113 ifstream refFile;
|
|
|
114 refFile.open(refname.c_str());
|
|
|
115 if (refFile.is_open()) {
|
|
|
116 while (getline(refFile, line)) {
|
|
|
117 ++line_number;
|
|
|
118 line_length = line.length();
|
|
|
119 if (line[0] == ';') {
|
|
|
120 // fasta comment, skip
|
|
|
121 } else if (line[0] == '+') {
|
|
|
122 // fastq quality header
|
|
|
123 getline(refFile, line);
|
|
|
124 line_length = line.length();
|
|
|
125 offset += line_length + 1;
|
|
|
126 // get and don't handle the quality line
|
|
|
127 getline(refFile, line);
|
|
|
128 line_length = line.length();
|
|
|
129 } else if (line[0] == '>' || line[0] == '@') { // fasta /fastq header
|
|
|
130 // if we aren't on the first entry, push the last sequence into the index
|
|
|
131 if (entry.name != "") {
|
|
|
132 mismatchedLineLengths = false; // reset line length error tracker for every new sequence
|
|
|
133 emptyLine = false;
|
|
|
134 flushEntryToIndex(entry);
|
|
|
135 entry.clear();
|
|
|
136 }
|
|
|
137 entry.name = line.substr(1, line_length - 1);
|
|
|
138 } else { // we assume we have found a sequence line
|
|
|
139 if (entry.offset == -1) // NB initially the offset is -1
|
|
|
140 entry.offset = offset;
|
|
|
141 entry.length += line_length;
|
|
|
142 if (entry.line_len) {
|
|
|
143 //entry.line_len = entry.line_len ? entry.line_len : line_length + 1;
|
|
|
144 if (mismatchedLineLengths || emptyLine) {
|
|
|
145 if (line_length == 0) {
|
|
|
146 emptyLine = true; // flag empty lines, raise error only if this is embedded in the sequence
|
|
|
147 } else {
|
|
|
148 if (emptyLine) {
|
|
|
149 cerr << "ERROR: embedded newline";
|
|
|
150 } else {
|
|
|
151 cerr << "ERROR: mismatched line lengths";
|
|
|
152 }
|
|
|
153 cerr << " at line " << line_number << " within sequence " << entry.name <<
|
|
|
154 endl << "File not suitable for fasta index generation." << endl;
|
|
|
155 exit(1);
|
|
|
156 }
|
|
|
157 }
|
|
|
158 // this flag is set here and checked on the next line
|
|
|
159 // because we may have reached the end of the sequence, in
|
|
|
160 // which case a mismatched line length is OK
|
|
|
161 if (entry.line_len != line_length + 1) {
|
|
|
162 mismatchedLineLengths = true;
|
|
|
163 if (line_length == 0) {
|
|
|
164 emptyLine = true; // flag empty lines, raise error only if this is embedded in the sequence
|
|
|
165 }
|
|
|
166 }
|
|
|
167 } else {
|
|
|
168 entry.line_len = line_length + 1; // first line
|
|
|
169 }
|
|
|
170 entry.line_blen = entry.line_len - 1;
|
|
|
171 }
|
|
|
172 offset += line_length + 1;
|
|
|
173 }
|
|
|
174 // we've hit the end of the fasta file!
|
|
|
175 // flush the last entry
|
|
|
176 flushEntryToIndex(entry);
|
|
|
177 } else {
|
|
|
178 cerr << "could not open reference file " << refname << " for indexing!" << endl;
|
|
|
179 exit(1);
|
|
|
180 }
|
|
|
181 }
|
|
|
182
|
|
|
183 void FastaIndex::flushEntryToIndex(FastaIndexEntry& entry) {
|
|
|
184 string name = split(entry.name, " \t").at(0); // key by first token of name
|
|
|
185 sequenceNames.push_back(name);
|
|
|
186 this->insert(make_pair(name, FastaIndexEntry(entry.name, entry.length,
|
|
|
187 entry.offset, entry.line_blen,
|
|
|
188 entry.line_len)));
|
|
|
189
|
|
|
190 }
|
|
|
191
|
|
|
192 void FastaIndex::writeIndexFile(string fname) {
|
|
|
193 //cerr << "writing fasta index file " << fname << endl;
|
|
|
194 ofstream file;
|
|
|
195 file.open(fname.c_str());
|
|
|
196 if (file.is_open()) {
|
|
|
197 file << *this;
|
|
|
198 } else {
|
|
|
199 cerr << "could not open index file " << fname << " for writing!" << endl;
|
|
|
200 exit(1);
|
|
|
201 }
|
|
|
202 }
|
|
|
203
|
|
|
204 FastaIndex::~FastaIndex(void) {
|
|
|
205 indexFile.close();
|
|
|
206 }
|
|
|
207
|
|
|
208 FastaIndexEntry FastaIndex::entry(string name) {
|
|
|
209 FastaIndex::iterator e = this->find(name);
|
|
|
210 if (e == this->end()) {
|
|
|
211 cerr << "unable to find FASTA index entry for '" << name << "'" << endl;
|
|
|
212 exit(1);
|
|
|
213 } else {
|
|
|
214 return e->second;
|
|
|
215 }
|
|
|
216 }
|
|
|
217
|
|
|
218 string FastaIndex::indexFileExtension() { return ".fai"; }
|
|
|
219
|
|
|
220 void FastaReference::open(string reffilename, bool usemmap) {
|
|
|
221 filename = reffilename;
|
|
|
222 if (!(file = fopen(filename.c_str(), "r"))) {
|
|
|
223 cerr << "could not open " << filename << endl;
|
|
|
224 exit(1);
|
|
|
225 }
|
|
|
226 index = new FastaIndex();
|
|
|
227 struct stat stFileInfo;
|
|
|
228 string indexFileName = filename + index->indexFileExtension();
|
|
|
229 // if we can find an index file, use it
|
|
|
230 if(stat(indexFileName.c_str(), &stFileInfo) == 0) {
|
|
|
231 index->readIndexFile(indexFileName);
|
|
|
232 } else { // otherwise, read the reference and generate the index file in the cwd
|
|
|
233 cerr << "index file " << indexFileName << " not found, generating..." << endl;
|
|
|
234 index->indexReference(filename);
|
|
|
235 index->writeIndexFile(indexFileName);
|
|
|
236 }
|
|
|
237 if (usemmap) {
|
|
|
238 usingmmap = true;
|
|
|
239 int fd = fileno(file);
|
|
|
240 struct stat sb;
|
|
|
241 if (fstat(fd, &sb) == -1)
|
|
|
242 cerr << "could not stat file" << filename << endl;
|
|
|
243 filesize = sb.st_size;
|
|
|
244 // map the whole file
|
|
|
245 filemm = mmap(NULL, filesize, PROT_READ, MAP_SHARED, fd, 0);
|
|
|
246 }
|
|
|
247 }
|
|
|
248
|
|
|
249 FastaReference::~FastaReference(void) {
|
|
|
250 fclose(file);
|
|
|
251 if (usingmmap) {
|
|
|
252 munmap(filemm, filesize);
|
|
|
253 }
|
|
|
254 delete index;
|
|
|
255 }
|
|
|
256
|
|
|
257 string FastaReference::getSequence(string seqname) {
|
|
|
258 FastaIndexEntry entry = index->entry(seqname);
|
|
|
259 int newlines_in_sequence = entry.length / entry.line_blen;
|
|
|
260 int seqlen = newlines_in_sequence + entry.length;
|
|
|
261 char* seq = (char*) calloc (seqlen + 1, sizeof(char));
|
|
|
262 if (usingmmap) {
|
|
|
263 memcpy(seq, (char*) filemm + entry.offset, seqlen);
|
|
|
264 } else {
|
|
|
265 fseek64(file, entry.offset, SEEK_SET);
|
|
|
266 fread(seq, sizeof(char), seqlen, file);
|
|
|
267 }
|
|
|
268 seq[seqlen] = '\0';
|
|
|
269 char* pbegin = seq;
|
|
|
270 char* pend = seq + (seqlen/sizeof(char));
|
|
|
271 pend = remove(pbegin, pend, '\n');
|
|
|
272 pend = remove(pbegin, pend, '\0');
|
|
|
273 string s = seq;
|
|
|
274 free(seq);
|
|
|
275 s.resize((pend - pbegin)/sizeof(char));
|
|
|
276 return s;
|
|
|
277 }
|
|
|
278
|
|
|
279 // TODO cleanup; odd function. use a map
|
|
|
280 string FastaReference::sequenceNameStartingWith(string seqnameStart) {
|
|
|
281 try {
|
|
|
282 return (*index)[seqnameStart].name;
|
|
|
283 } catch (exception& e) {
|
|
|
284 cerr << e.what() << ": unable to find index entry for " << seqnameStart << endl;
|
|
|
285 exit(1);
|
|
|
286 }
|
|
|
287 }
|
|
|
288
|
|
|
289 string FastaReference::getSubSequence(string seqname, int start, int length) {
|
|
|
290 FastaIndexEntry entry = index->entry(seqname);
|
|
|
291 if (start < 0 || length < 1) {
|
|
|
292 cerr << "Error: cannot construct subsequence with negative offset or length < 1" << endl;
|
|
|
293 exit(1);
|
|
|
294 }
|
|
|
295 // we have to handle newlines
|
|
|
296 // approach: count newlines before start
|
|
|
297 // count newlines by end of read
|
|
|
298 // subtracting newlines before start find count of embedded newlines
|
|
|
299 int newlines_before = start > 0 ? (start - 1) / entry.line_blen : 0;
|
|
|
300 int newlines_by_end = (start + length - 1) / entry.line_blen;
|
|
|
301 int newlines_inside = newlines_by_end - newlines_before;
|
|
|
302 int seqlen = length + newlines_inside;
|
|
|
303 char* seq = (char*) calloc (seqlen + 1, sizeof(char));
|
|
|
304 if (usingmmap) {
|
|
|
305 memcpy(seq, (char*) filemm + entry.offset + newlines_before + start, seqlen);
|
|
|
306 } else {
|
|
|
307 fseek64(file, (off_t) (entry.offset + newlines_before + start), SEEK_SET);
|
|
|
308 fread(seq, sizeof(char), (off_t) seqlen, file);
|
|
|
309 }
|
|
|
310 seq[seqlen] = '\0';
|
|
|
311 char* pbegin = seq;
|
|
|
312 char* pend = seq + (seqlen/sizeof(char));
|
|
|
313 pend = remove(pbegin, pend, '\n');
|
|
|
314 pend = remove(pbegin, pend, '\0');
|
|
|
315 string s = seq;
|
|
|
316 free(seq);
|
|
|
317 s.resize((pend - pbegin)/sizeof(char));
|
|
|
318 return s;
|
|
|
319 }
|
|
|
320
|
|
|
321 long unsigned int FastaReference::sequenceLength(string seqname) {
|
|
|
322 FastaIndexEntry entry = index->entry(seqname);
|
|
|
323 return entry.length;
|
|
|
324 }
|
|
|
325
|
