gclib/gclib/GFaSeqGet.cpp

#include "GFaSeqGet.h"
#include "gdna.h"
#include <ctype.h>

void GSubSeq::setup(uint sstart, int slen, int sovl, int qfrom, int qto, uint maxseqlen) {
     if (sovl==0) {
       GFREE(sq);
       sqstart=sstart;
       uint max_len=(maxseqlen>0) ? maxseqlen : MAX_FASUBSEQ;
       sqlen = (slen==0 ? max_len : slen);
       GMALLOC(sq, sqlen);
       return;
       }
  //overlap -- copy the overlapping region
  char* newsq=NULL;
  GMALLOC(newsq, slen);
  memcpy((void*)&newsq[qto], (void*)&sq[qfrom], sovl);
  GFREE(sq);
  sq=newsq;
  sqstart=sstart;
  sqlen=slen;
}

void GFaSeqGet::finit(const char* fn, off_t fofs, bool validate) {
 fh=fopen(fn,"rb");
 if (fh==NULL) {
   GError("Error (GFaSeqGet) opening file '%s'\n",fn);
   }
 fname=Gstrdup(fn);
 initialParse(fofs, validate);
 lastsub=new GSubSeq();
}

GFaSeqGet::GFaSeqGet(const char* faname, uint seqlen, off_t fseqofs, int l_len, int l_blen) {
//for GFastaIndex use mostly -- the important difference is that
//the file offset is to the sequence, not to the defline
  fh=fopen(faname,"rb");
  if (fh==NULL) {
    GError("Error (GFaSeqGet) opening file '%s'\n",faname);
    }
  fname=Gstrdup(faname);
  line_len=l_len;
  line_blen=l_blen;
  seq_len=seqlen;
  if (line_blen<line_len)
       GError("Error (GFaSeqGet): invalid line length info (len=%d, blen=%d)\n",
              line_len, line_blen);
  fseqstart=fseqofs;
  lastsub=new GSubSeq();
}

GFaSeqGet::GFaSeqGet(FILE* f, off_t fofs, bool validate) {
 if (f==NULL) GError("Error (GFaSeqGet) : null file handle!\n");
 seq_len=0;
 fh=f;
 initialParse(fofs, validate);
 lastsub=new GSubSeq();
}

void GFaSeqGet::initialParse(off_t fofs, bool checkall) {
 static const char gfa_ERRPARSE[]="Error (GFaSeqGet): invalid FASTA file format.\n";
 if (fofs!=0) { fseeko(fh,fofs,SEEK_SET); } //e.g. for offsets provided by cdbyank
 //read the first two lines to determine fasta parameters
 fseqstart=fofs;
 int c=getc(fh);
 fseqstart++;
 if (c!='>') GError("Error (GFaSeqGet): not a fasta header?\n");
 while ((c=getc(fh))!=EOF) {
   fseqstart++;
   if (c=='\n' || c=='\r') { break; } //end of defline
   }

 if (c==EOF) GError(gfa_ERRPARSE);
 line_len=0;
 int lendlen=0;
 while ((c=getc(fh))!=EOF) {
  if (c=='\n' || c=='\r') { //end of line encountered
     if (line_len>0) { //end of the first "sequence" line
        lendlen++;
        break;
        }
      else {// another EoL char at the end of defline
        fseqstart++;
        continue;
        }
     }// end-of-line characters
  line_len++;
  }
 //we are at the end of first sequence line
 while ((c=getc(fh))!=EOF) {
   if (c=='\n' || c=='\r') lendlen++;
      else {
       ungetc(c,fh);
       break;
       }
   }
 line_blen=line_len+lendlen;
 if (c==EOF) return;
 // -- you don't need to check it all if you're sure it's safe
 if (checkall) { //validate the rest of the FASTA record
   int llen=0; //last line length
   int elen=0; //length of last line ending
   bool waseol=true;
   while ((c=getc(fh))!=EOF) {
     if (c=='>' && waseol) { ungetc(c,fh); break; }
     if (c=='\n' ||  c=='\r') {
        // eol char
        elen++;
        if (waseol) continue; //2nd eol char
        waseol=true;
        elen=1;
        continue;
        }
     if (c<=32) GError(gfa_ERRPARSE); //invalid character encountered
     //--- on a seq char here:
     if (waseol) {//beginning of a seq line
       if (elen && (llen!=line_len || elen!=lendlen))
           //GError(gfa_ERRPARSE);
         GError("Error: invalid FASTA format for GFaSeqGet; make sure that\n\
  the sequence lines have the same length (except for the last line)");
       waseol=false;
       llen=0;
       elen=0;
       }
     llen++;
     } //while reading chars
   }// FASTA checking was requested
 fseeko(fh,fseqstart,SEEK_SET);
}

const char* GFaSeqGet::subseq(uint cstart, int& clen) {
  //cstart is 1-based genomic coordinate within current fasta sequence
   int maxlen=(seq_len>0)?seq_len : MAX_FASUBSEQ;
   //GMessage("--> call: subseq(%u, %d)\n", cstart, clen);
  if (clen>maxlen) {
    GMessage("Error (GFaSeqGet): subsequence cannot be larger than %d\n", maxlen);
    return NULL;
    }
  if (seq_len>0 && clen+cstart-1>seq_len) {
     GMessage("Error (GFaSeqGet): end coordinate (%d) cannot be larger than sequence length %d\n", clen+cstart-1, seq_len);
     }
  if (lastsub->sq==NULL || lastsub->sqlen==0) {
    lastsub->setup(cstart, clen, 0,0,0,seq_len);
    loadsubseq(cstart, clen);
    lastsub->sqlen=clen;
    return (const char*)lastsub->sq;
    }
  //allow extension up to MAX_FASUBSEQ
  uint bstart=lastsub->sqstart;
  uint bend=lastsub->sqstart+lastsub->sqlen-1;
  uint cend=cstart+clen-1;
  int qlen=0; //only the extra len to be allocated/appended/prepended
  uint qstart=cstart; //start coordinate of the new seq block of length qlen to be read from file
  int newlen=0; //the new total length of the buffered sequence lastsub->sq
  int kovl=0;
  int czfrom=0;//0-based offsets for copying a previously read sequence chunk
  int czto=0;
  uint newstart=cstart;
  if (cstart>=bstart && cend<=bend) { //new reg contained within existing buffer
     return (const char*) &(lastsub->sq[cstart-bstart]) ;
    }
  //extend downward
  uint newend=GMAX(cend, bend);
  if (cstart<bstart) { //requested start < old buffer start
    newstart=cstart;
    newlen=(newend-newstart+1);
    if (newlen>MAX_FASUBSEQ) {
       newlen=MAX_FASUBSEQ;
       newend=cstart+newlen-1; //keep newstart, set newend
       }
    qlen=bstart-cstart;
    if (newend>bstart) { //overlap
       if (newend>bend) {// new region is larger & around the old one - so we have two regions to update
         kovl=bend-bstart+1;
         czfrom=0;
         czto=bstart-cstart;
         lastsub->setup(newstart, newlen, kovl, czfrom, czto, seq_len); //this should realloc and copy the kovl subseq
         qlen=bstart-cstart;
         loadsubseq(newstart, qlen);
         qlen=newend-bend;
         int toread=qlen;
         loadsubseq(bend+1, qlen);
         clen-=(toread-qlen);
         lastsub->sqlen=clen;
         return (const char*)lastsub->sq;
         }
        //newend<=bend
       kovl=newend-bstart+1;
       }
     else { //no overlap with previous buffer
       if (newend>bend) kovl=bend-bstart+1;
                   else kovl=newend-bstart+1;
       }
     qlen=bstart-cstart;
     czfrom=0;
     czto=qlen;
    } //cstart<bstart
   else { //cstart>=bstart, possibly extend upwards
    newstart=bstart;
    newlen=(newend-newstart+1);
    if (newlen>MAX_FASUBSEQ) {
       newstart=bstart+(newlen-MAX_FASUBSEQ);//keep newend, assign newstart
       newlen=MAX_FASUBSEQ;
       if (newstart<=bend) { //overlap with old buffer
          kovl=bend-newstart+1;
          czfrom=newstart-bstart;
          czto=0;
          }
       else { //not overlapping old buffer
         kovl=0;
         }
       } //newstart reassigned
    else { //we can extend the buffer to include the old one
      qlen=newend-bend; //how much to read from file
      qstart=bend+1;
      kovl=bend-bstart+1;
      czfrom=0;
      czto=0;
      }
    }
  lastsub->setup(newstart, newlen, kovl, czfrom, czto, seq_len); //this should realloc but copy any overlapping region
  lastsub->sqlen-=qlen; //appending may result in a premature eof
  int toread=qlen;
  loadsubseq(qstart, qlen); //read the missing chunk, if any
  clen-=(toread-qlen);
  lastsub->sqlen+=qlen;
  return (const char*)(lastsub->sq+(cstart-newstart));
}

char* GFaSeqGet::copyRange(uint cstart, uint cend, bool revCmpl, bool upCase) {
  if (cstart>cend) { Gswap(cstart, cend); }
  int clen=cend-cstart+1;
  const char* gs=subseq(cstart, clen);
  if (gs==NULL) return NULL;
  char* r=NULL;
  GMALLOC(r,clen+1);
  r[clen]=0;
  memcpy((void*)r,(void*)gs, clen);
  if (revCmpl) reverseComplement(r,clen);
  if (upCase) {
       for (int i=0;i<clen;i++)
            r[i]=toupper(r[i]);
       }
  return r;
 }

const char* GFaSeqGet::loadsubseq(uint cstart, int& clen) {
  //assumes enough lastsub->sq space allocated previously
  //only loads the requested clen chars from file, at offset &lastsub->sq[cstart-lastsub->sqstart]
  int sofs=cstart-lastsub->sqstart;
  int lendlen=line_blen-line_len;
  char* seqp=lastsub->sq+sofs;
  //find the proper file offset and read the appropriate lines
  uint seqofs=cstart-1;
  uint startlno = seqofs/line_len;
  int lineofs = seqofs % line_len;
  off_t fstart=fseqstart + (startlno*line_blen);
  fstart+=lineofs;
  
  fseeko(fh, fstart, SEEK_SET);
  int toread=clen;
  int maxlen=(seq_len>0)? seq_len-cstart+1 : MAX_FASUBSEQ ;
  if (toread==0) toread=maxlen; //read max allowed, or to the end of file
  int actualrlen=0;
  int sublen=0;
  if (lineofs>0) { //read the partial first line
    int reqrlen=line_len-lineofs;
    if (reqrlen>toread) reqrlen=toread; //in case we need to read just a few chars
    actualrlen=fread((void*)seqp, 1, reqrlen, fh);
    if (actualrlen<reqrlen) { //eof reached prematurely
      while (seqp[actualrlen-1]=='\n' || seqp[actualrlen-1]=='\r') actualrlen--;
      //check for new sequences in between
      clen=actualrlen;
      sublen+=actualrlen;
      return (const char*)seqp;
      }
    toread-=reqrlen;
    sublen+=reqrlen;
    fseeko(fh, lendlen, SEEK_CUR);
    }
  //read the rest of the lines
  while (toread>=line_len) {
    char* rseqp=&(seqp[sublen]);
    actualrlen=fread((void*)rseqp, 1, line_len, fh);
    /*
    char dbuf[256];dbuf[255]=0;
    strncpy(dbuf,rseqp, actualrlen);
    dbuf[actualrlen]=0;
    GMessage("<<<read line: %s\n",dbuf);
    */
    if (actualrlen<line_len) {
      while (rseqp[actualrlen-1]=='\n' || rseqp[actualrlen-1]=='\r') actualrlen--;
      sublen+=actualrlen;
      clen=sublen;
      return (const char*)seqp;
      }
    toread-=actualrlen;
    sublen+=actualrlen;
    fseeko(fh, lendlen, SEEK_CUR);
    }
  // read the last partial line, if any
  if (toread>0) {
    char* rseqp=&(seqp[sublen]);
    actualrlen=fread((void*)rseqp, 1, toread, fh);
    if (actualrlen<toread) {
      while (rseqp[actualrlen-1]=='\n' || rseqp[actualrlen-1]=='\r')
          actualrlen--;
      }
    sublen+=actualrlen;
    }
  //lastsub->sqlen+=sublen;
  clen=sublen;
  
  return (const char*)seqp;
  }
Revision:	144
Committed:	Thu Dec 22 19:15:46 2011 UTC (12 years, 8 months ago) by gpertea
File size:	10196 byte(s)
Log Message:	added GBitVec.h; removed operator> requirement for sorted vectors; swap() is now template Gswap()
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) { Gswap(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	}