Mercurial > repos > aaronquinlan > multi_intersect

comparison BEDTools-Version-2.14.3/src/utils/Fasta/Fasta.cpp @ 0:dfcd8b6c1bda

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author aaronquinlan
date Thu, 03 Nov 2011 10:25:04 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:dfcd8b6c1bda
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