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