ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/gclib/GFaSeqGet.cpp
Revision: 16
Committed: Mon Jul 18 20:56:02 2011 UTC (8 years, 4 months ago) by gpertea
File size: 10195 byte(s)
Log Message:
sync with local source

Line File contents
1 #include "GFaSeqGet.h"
2 #include "gdna.h"
3 #include <ctype.h>
4
5 void GSubSeq::setup(uint sstart, int slen, int sovl, int qfrom, int qto, uint maxseqlen) {
6 if (sovl==0) {
7 GFREE(sq);
8 sqstart=sstart;
9 uint max_len=(maxseqlen>0) ? maxseqlen : MAX_FASUBSEQ;
10 sqlen = (slen==0 ? max_len : slen);
11 GMALLOC(sq, sqlen);
12 return;
13 }
14 //overlap -- copy the overlapping region
15 char* newsq=NULL;
16 GMALLOC(newsq, slen);
17 memcpy((void*)&newsq[qto], (void*)&sq[qfrom], sovl);
18 GFREE(sq);
19 sq=newsq;
20 sqstart=sstart;
21 sqlen=slen;
22 }
23
24 void GFaSeqGet::finit(const char* fn, off_t fofs, bool validate) {
25 fh=fopen(fn,"rb");
26 if (fh==NULL) {
27 GError("Error (GFaSeqGet) opening file '%s'\n",fn);
28 }
29 fname=Gstrdup(fn);
30 initialParse(fofs, validate);
31 lastsub=new GSubSeq();
32 }
33
34 GFaSeqGet::GFaSeqGet(const char* faname, uint seqlen, off_t fseqofs, int l_len, int l_blen) {
35 //for GFastaIndex use mostly -- the important difference is that
36 //the file offset is to the sequence, not to the defline
37 fh=fopen(faname,"rb");
38 if (fh==NULL) {
39 GError("Error (GFaSeqGet) opening file '%s'\n",faname);
40 }
41 fname=Gstrdup(faname);
42 line_len=l_len;
43 line_blen=l_blen;
44 seq_len=seqlen;
45 if (line_blen<line_len)
46 GError("Error (GFaSeqGet): invalid line length info (len=%d, blen=%d)\n",
47 line_len, line_blen);
48 fseqstart=fseqofs;
49 lastsub=new GSubSeq();
50 }
51
52 GFaSeqGet::GFaSeqGet(FILE* f, off_t fofs, bool validate) {
53 if (f==NULL) GError("Error (GFaSeqGet) : null file handle!\n");
54 seq_len=0;
55 fh=f;
56 initialParse(fofs, validate);
57 lastsub=new GSubSeq();
58 }
59
60 void GFaSeqGet::initialParse(off_t fofs, bool checkall) {
61 static const char gfa_ERRPARSE[]="Error (GFaSeqGet): invalid FASTA file format.\n";
62 if (fofs!=0) { fseeko(fh,fofs,SEEK_SET); } //e.g. for offsets provided by cdbyank
63 //read the first two lines to determine fasta parameters
64 fseqstart=fofs;
65 int c=getc(fh);
66 fseqstart++;
67 if (c!='>') GError("Error (GFaSeqGet): not a fasta header?\n");
68 while ((c=getc(fh))!=EOF) {
69 fseqstart++;
70 if (c=='\n' || c=='\r') { break; } //end of defline
71 }
72
73 if (c==EOF) GError(gfa_ERRPARSE);
74 line_len=0;
75 int lendlen=0;
76 while ((c=getc(fh))!=EOF) {
77 if (c=='\n' || c=='\r') { //end of line encountered
78 if (line_len>0) { //end of the first "sequence" line
79 lendlen++;
80 break;
81 }
82 else {// another EoL char at the end of defline
83 fseqstart++;
84 continue;
85 }
86 }// end-of-line characters
87 line_len++;
88 }
89 //we are at the end of first sequence line
90 while ((c=getc(fh))!=EOF) {
91 if (c=='\n' || c=='\r') lendlen++;
92 else {
93 ungetc(c,fh);
94 break;
95 }
96 }
97 line_blen=line_len+lendlen;
98 if (c==EOF) return;
99 // -- you don't need to check it all if you're sure it's safe
100 if (checkall) { //validate the rest of the FASTA record
101 int llen=0; //last line length
102 int elen=0; //length of last line ending
103 bool waseol=true;
104 while ((c=getc(fh))!=EOF) {
105 if (c=='>' && waseol) { ungetc(c,fh); break; }
106 if (c=='\n' || c=='\r') {
107 // eol char
108 elen++;
109 if (waseol) continue; //2nd eol char
110 waseol=true;
111 elen=1;
112 continue;
113 }
114 if (c<=32) GError(gfa_ERRPARSE); //invalid character encountered
115 //--- on a seq char here:
116 if (waseol) {//beginning of a seq line
117 if (elen && (llen!=line_len || elen!=lendlen))
118 //GError(gfa_ERRPARSE);
119 GError("Error: invalid FASTA format for GFaSeqGet; make sure that\n\
120 the sequence lines have the same length (except for the last line)");
121 waseol=false;
122 llen=0;
123 elen=0;
124 }
125 llen++;
126 } //while reading chars
127 }// FASTA checking was requested
128 fseeko(fh,fseqstart,SEEK_SET);
129 }
130
131 const char* GFaSeqGet::subseq(uint cstart, int& clen) {
132 //cstart is 1-based genomic coordinate within current fasta sequence
133 int maxlen=(seq_len>0)?seq_len : MAX_FASUBSEQ;
134 //GMessage("--> call: subseq(%u, %d)\n", cstart, clen);
135 if (clen>maxlen) {
136 GMessage("Error (GFaSeqGet): subsequence cannot be larger than %d\n", maxlen);
137 return NULL;
138 }
139 if (seq_len>0 && clen+cstart-1>seq_len) {
140 GMessage("Error (GFaSeqGet): end coordinate (%d) cannot be larger than sequence length %d\n", clen+cstart-1, seq_len);
141 }
142 if (lastsub->sq==NULL || lastsub->sqlen==0) {
143 lastsub->setup(cstart, clen, 0,0,0,seq_len);
144 loadsubseq(cstart, clen);
145 lastsub->sqlen=clen;
146 return (const char*)lastsub->sq;
147 }
148 //allow extension up to MAX_FASUBSEQ
149 uint bstart=lastsub->sqstart;
150 uint bend=lastsub->sqstart+lastsub->sqlen-1;
151 uint cend=cstart+clen-1;
152 int qlen=0; //only the extra len to be allocated/appended/prepended
153 uint qstart=cstart; //start coordinate of the new seq block of length qlen to be read from file
154 int newlen=0; //the new total length of the buffered sequence lastsub->sq
155 int kovl=0;
156 int czfrom=0;//0-based offsets for copying a previously read sequence chunk
157 int czto=0;
158 uint newstart=cstart;
159 if (cstart>=bstart && cend<=bend) { //new reg contained within existing buffer
160 return (const char*) &(lastsub->sq[cstart-bstart]) ;
161 }
162 //extend downward
163 uint newend=GMAX(cend, bend);
164 if (cstart<bstart) { //requested start < old buffer start
165 newstart=cstart;
166 newlen=(newend-newstart+1);
167 if (newlen>MAX_FASUBSEQ) {
168 newlen=MAX_FASUBSEQ;
169 newend=cstart+newlen-1; //keep newstart, set newend
170 }
171 qlen=bstart-cstart;
172 if (newend>bstart) { //overlap
173 if (newend>bend) {// new region is larger & around the old one - so we have two regions to update
174 kovl=bend-bstart+1;
175 czfrom=0;
176 czto=bstart-cstart;
177 lastsub->setup(newstart, newlen, kovl, czfrom, czto, seq_len); //this should realloc and copy the kovl subseq
178 qlen=bstart-cstart;
179 loadsubseq(newstart, qlen);
180 qlen=newend-bend;
181 int toread=qlen;
182 loadsubseq(bend+1, qlen);
183 clen-=(toread-qlen);
184 lastsub->sqlen=clen;
185 return (const char*)lastsub->sq;
186 }
187 //newend<=bend
188 kovl=newend-bstart+1;
189 }
190 else { //no overlap with previous buffer
191 if (newend>bend) kovl=bend-bstart+1;
192 else kovl=newend-bstart+1;
193 }
194 qlen=bstart-cstart;
195 czfrom=0;
196 czto=qlen;
197 } //cstart<bstart
198 else { //cstart>=bstart, possibly extend upwards
199 newstart=bstart;
200 newlen=(newend-newstart+1);
201 if (newlen>MAX_FASUBSEQ) {
202 newstart=bstart+(newlen-MAX_FASUBSEQ);//keep newend, assign newstart
203 newlen=MAX_FASUBSEQ;
204 if (newstart<=bend) { //overlap with old buffer
205 kovl=bend-newstart+1;
206 czfrom=newstart-bstart;
207 czto=0;
208 }
209 else { //not overlapping old buffer
210 kovl=0;
211 }
212 } //newstart reassigned
213 else { //we can extend the buffer to include the old one
214 qlen=newend-bend; //how much to read from file
215 qstart=bend+1;
216 kovl=bend-bstart+1;
217 czfrom=0;
218 czto=0;
219 }
220 }
221 lastsub->setup(newstart, newlen, kovl, czfrom, czto, seq_len); //this should realloc but copy any overlapping region
222 lastsub->sqlen-=qlen; //appending may result in a premature eof
223 int toread=qlen;
224 loadsubseq(qstart, qlen); //read the missing chunk, if any
225 clen-=(toread-qlen);
226 lastsub->sqlen+=qlen;
227 return (const char*)(lastsub->sq+(cstart-newstart));
228 }
229
230 char* GFaSeqGet::copyRange(uint cstart, uint cend, bool revCmpl, bool upCase) {
231 if (cstart>cend) { swap(cstart, cend); }
232 int clen=cend-cstart+1;
233 const char* gs=subseq(cstart, clen);
234 if (gs==NULL) return NULL;
235 char* r=NULL;
236 GMALLOC(r,clen+1);
237 r[clen]=0;
238 memcpy((void*)r,(void*)gs, clen);
239 if (revCmpl) reverseComplement(r,clen);
240 if (upCase) {
241 for (int i=0;i<clen;i++)
242 r[i]=toupper(r[i]);
243 }
244 return r;
245 }
246
247 const char* GFaSeqGet::loadsubseq(uint cstart, int& clen) {
248 //assumes enough lastsub->sq space allocated previously
249 //only loads the requested clen chars from file, at offset &lastsub->sq[cstart-lastsub->sqstart]
250 int sofs=cstart-lastsub->sqstart;
251 int lendlen=line_blen-line_len;
252 char* seqp=lastsub->sq+sofs;
253 //find the proper file offset and read the appropriate lines
254 uint seqofs=cstart-1;
255 uint startlno = seqofs/line_len;
256 int lineofs = seqofs % line_len;
257 off_t fstart=fseqstart + (startlno*line_blen);
258 fstart+=lineofs;
259
260 fseeko(fh, fstart, SEEK_SET);
261 int toread=clen;
262 int maxlen=(seq_len>0)? seq_len-cstart+1 : MAX_FASUBSEQ ;
263 if (toread==0) toread=maxlen; //read max allowed, or to the end of file
264 int actualrlen=0;
265 int sublen=0;
266 if (lineofs>0) { //read the partial first line
267 int reqrlen=line_len-lineofs;
268 if (reqrlen>toread) reqrlen=toread; //in case we need to read just a few chars
269 actualrlen=fread((void*)seqp, 1, reqrlen, fh);
270 if (actualrlen<reqrlen) { //eof reached prematurely
271 while (seqp[actualrlen-1]=='\n' || seqp[actualrlen-1]=='\r') actualrlen--;
272 //check for new sequences in between
273 clen=actualrlen;
274 sublen+=actualrlen;
275 return (const char*)seqp;
276 }
277 toread-=reqrlen;
278 sublen+=reqrlen;
279 fseeko(fh, lendlen, SEEK_CUR);
280 }
281 //read the rest of the lines
282 while (toread>=line_len) {
283 char* rseqp=&(seqp[sublen]);
284 actualrlen=fread((void*)rseqp, 1, line_len, fh);
285 /*
286 char dbuf[256];dbuf[255]=0;
287 strncpy(dbuf,rseqp, actualrlen);
288 dbuf[actualrlen]=0;
289 GMessage("<<<read line: %s\n",dbuf);
290 */
291 if (actualrlen<line_len) {
292 while (rseqp[actualrlen-1]=='\n' || rseqp[actualrlen-1]=='\r') actualrlen--;
293 sublen+=actualrlen;
294 clen=sublen;
295 return (const char*)seqp;
296 }
297 toread-=actualrlen;
298 sublen+=actualrlen;
299 fseeko(fh, lendlen, SEEK_CUR);
300 }
301 // read the last partial line, if any
302 if (toread>0) {
303 char* rseqp=&(seqp[sublen]);
304 actualrlen=fread((void*)rseqp, 1, toread, fh);
305 if (actualrlen<toread) {
306 while (rseqp[actualrlen-1]=='\n' || rseqp[actualrlen-1]=='\r')
307 actualrlen--;
308 }
309 sublen+=actualrlen;
310 }
311 //lastsub->sqlen+=sublen;
312 clen=sublen;
313
314 return (const char*)seqp;
315 }