---

gnGBKSource.cpp

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
// File:            gnGBKSource.h
// Purpose:         Implements gnBaseSource for GenBank sequences
// Description:     
// Changes:        
// Version:         libGenome 0.5.1 
// Author:          Aaron Darling 
// Modified by:     
// Copyright:       (c) Aaron Darling 
// Licenses:        See COPYING file for details
/////////////////////////////////////////////////////////////////////////////
#include "gn/gnFilter.h"
#include "gn/gnFeature.h"
#include "gn/gnGBKSource.h"
#include "gn/gnSourceSpec.h"
#include "gn/gnSourceHeader.h"
#include "gn/gnSourceQualifier.h"
#include "gn/gnLocation.h"
#include "gn/gnStringTools.h"
#include "gn/gnDebug.h"
#include <string>

gnGBKSource::gnGBKSource()
{
        m_openString = "";
        m_pFilter = gnFilter::fullDNASeqFilter();
        if(m_pFilter == NULL){
                DebugMsg("Error using static sequence filters.");
        }
}
gnGBKSource::gnGBKSource( const gnGBKSource& s ) : gnFileSource(s)
{
        vector< gnFileContig* >::const_iterator iter = s.m_contigList.begin();
        for( ; iter != s.m_contigList.end(); ++iter )
        {
                m_contigList.push_back( (*iter)->Clone() );
        }
}
gnGBKSource::~gnGBKSource()
{
        m_ifstream.close();
        vector< gnFileContig* >::iterator iter = m_contigList.begin();
        for( ; iter != m_contigList.end(); ++iter )
        {
                gnFileContig* fg = *iter;
                *iter = 0;
                delete fg;
        }
}
boolean gnGBKSource::HasContig( const string& name ) const
{
        for(uint32 i = 0 ; i <= m_contigList.size(); i++ )
        {
                if( name == m_contigList[i]->GetName() )
                        return true;
        }
        return false;
}
uint32 gnGBKSource::GetContigID( const string& name ) const
{
        for(uint32 i = 0 ; i <= m_contigList.size(); i++ )
        {
                if( name == m_contigList[i]->GetName() )
                        return i;
        }
        return ALL_CONTIGS;
}
string gnGBKSource::GetContigName( const uint32 i ) const
{
        if( i < m_contigList.size() )
        {
                return m_contigList[i]->GetName();
        }
        return "";
}
gnSeqI gnGBKSource::GetContigSeqLength( const uint32 i ) const
{
        if( i == ALL_CONTIGS)
                return m_spec->GetLength();
        if( i < m_contigList.size() )
        {
                return m_contigList[i]->GetSeqLength();
        }
        return GNSEQI_ERROR;
}

boolean gnGBKSource::SeqRead( const gnSeqI start, char* buf, uint32& bufLen, const uint32 contigI ){
        uint64 startPos = 0;
        uint64 readableBytes = 0;
        if( !SeqSeek( start, contigI, startPos, readableBytes ) )
        {
                bufLen = 0;
                return false;
        }
        
        if( contigI == ALL_CONTIGS )
        {
                uint32 curLen = 0;
                uint64 bytesRead = 0;
                while (curLen < bufLen)
                {
//SeqSeek to start, Figure out how much can be read before SeqSeeking again.
                        if(readableBytes <= 0)  //Look out for zero length contigs!  IMPLEMENT ME
                                if( !SeqSeek( start + curLen, contigI, startPos, readableBytes ) ){
                                        bufLen = curLen;
                                        return true;
                                }
                        //readLen is the amount to read on this pass
                        uint64 readLen = (bufLen - curLen) < readableBytes ? (bufLen - curLen) : readableBytes; 
                        Array<gnSeqC> array_buf( readLen );
                        gnSeqC* tmpBuf = array_buf.data;

                        // read chars and filter
                        m_ifstream.read(tmpBuf, readLen);
                        uint64 gc = m_ifstream.gcount();
                        bytesRead += gc;
                        readableBytes -= gc;
                        for(uint32 i=0; i < gc; i++){
                                if( m_pFilter->IsValid(tmpBuf[i]) ){
                                        buf[curLen] = tmpBuf[i];
                                        curLen++;
                                }
                        }
                        if(m_ifstream.eof()){   //we hit the end of the file.  bail out.
                                m_ifstream.clear();
                                bufLen = curLen;
                                return true;
                        }
                }
                bufLen = curLen;
        }
        else if( contigI < m_contigList.size() )
        {
                uint32 curLen = 0;
                //check to see if the buffer is bigger than the contig.  if so truncate it.
                gnSeqI contigSize = m_contigList[contigI]->GetSeqLength();
                bufLen = bufLen < contigSize ? bufLen : contigSize;
                while (curLen < bufLen)
                {
                        uint64 readLen = bufLen - curLen;       //the amount to read on this pass
                        Array<gnSeqC> array_buf( readLen );
                        gnSeqC* tmpBuf = array_buf.data;

                        // read chars and filter
                        m_ifstream.read(tmpBuf, readLen);
                        uint64 gc = m_ifstream.gcount();
//                      cout << "Read " << gc << " chars from " << m_openString << "\n";
//                      cout << "Checking character validity on: " << tmpBuf << "\n";
                        for(uint32 i=0; i < gc; i++){
                                if( m_pFilter->IsValid(tmpBuf[i]) ){
                                        buf[curLen] = tmpBuf[i];
                                        curLen++;
                                }
                        }
                        if(m_ifstream.eof()){   //we hit the end of the file.  bail out.
                                m_ifstream.clear();
                                bufLen = curLen;
                                return true;
                        }
                }
                bufLen = curLen;
        }
        return true;

}
// private:
// figures out which contig the sequence starts at then calls SeqStartPos to get the offset within that contig
// returns startPos, the file offset where the sequence starts
// returns true if successful, false otherwise
boolean gnGBKSource::SeqSeek( const gnSeqI start, const uint32& contigI, uint64& startPos, uint64& readableBytes )
{
        if( contigI == ALL_CONTIGS )
        {
                // find first contig
                gnSeqI curIndex = 0;
                vector< gnFileContig* >::iterator iter = m_contigList.begin();
                for( ; iter != m_contigList.end(); ++iter )
                {
                        uint64 len = (*iter)->GetSeqLength();
                        if( (curIndex + len) > start )
                                break;
                        curIndex += len;
                }
                if( iter == m_contigList.end() )
                        return false;
                // seek to start
                gnSeqI startIndex = start - curIndex;  //startIndex is starting pos. within the contig
                return SeqStartPos( startIndex, *(*iter), startPos, readableBytes );
        }
        else if( contigI < m_contigList.size() )
        {
                return SeqStartPos( start, *(m_contigList[contigI]), startPos, readableBytes );
        }
        return false;
}
//Returns startPos, the file offset where the sequence starts.
boolean gnGBKSource::SeqStartPos( const gnSeqI start, gnFileContig& contig, uint64& startPos, uint64& readableBytes )
{
        readableBytes = 0;
        uint32 curLen = 0;
        //seek to the file offset where the contig starts
        startPos = contig.GetSectStartEnd(gnContigSequence).first;      //set startPos to start where the contig starts
        m_ifstream.seekg( startPos, ios::beg );
        if( m_ifstream.eof() ){
                ErrorMsg("ERROR in gnGBKSource::Incorrect contig start position, End of file reached!\n");
                return false;
        }
        while( true )
        {
                  // READ the rest of the contig skipping over invalid characters until we get to the starting base pair.
                  // startPos will contain the file offset with the starting base pair
                uint32 tmpbufsize = contig.GetSectStartEnd(gnContigSequence).second - startPos;
                if(tmpbufsize == 0){
                        ErrorMsg("ERROR in gnGBKSource: stored contig size is incorrect.");
                        return false;
                }
                uint64 startOffset = start;
                if(contig.HasRepeatSeqGap()){   //check for sequence integrity
                        startOffset += (9 + m_newlineSize) * (start / 60 + 1) + start / 10 + 1;
                        startPos+=startOffset;
                        m_ifstream.seekg(startPos , ios::beg);
                        readableBytes = contig.GetSectStartEnd(gnContigSequence).second - startPos;
                        return true;
                }

                //sequence is corrupt.  read in base by base
                tmpbufsize = tmpbufsize < BUFFER_SIZE ? tmpbufsize : BUFFER_SIZE;  //read in the smaller of the two.
                Array<char> array_buf( tmpbufsize );
                char* tmpbuf = array_buf.data;

                m_ifstream.read( tmpbuf, tmpbufsize );
                if( m_ifstream.eof() ){
                        ErrorMsg("ERROR in gnGBKSource::Read End of file reached!\n");
                        return false;
                }
                for( uint32 i=0; i < tmpbufsize; ++i ){
                        if( m_pFilter->IsValid(tmpbuf[i]) ){
                                if( curLen >= start ){ //stop when we reach the starting offset within this contig
                                        startPos += i;
                                        m_ifstream.seekg( startPos, ios::beg );  //seek to startPos
                                        readableBytes = contig.GetSectStartEnd(gnContigSequence).second - startPos;
                                        return true;
                                }
                                ++curLen;  //each time we read a valid b.p., increment the sequence length
                        }
                }
                startPos += tmpbufsize;
        }
        return true;
}

void gnGBKSource::FormatString(string& data, uint32 offset, uint32 width){
        //first remove newlines and corresponding whitespace
        string::size_type newline_loc = data.find_first_of('\n', 0);
        while(newline_loc != string::npos){
                if(data[newline_loc-1] == '\r')
                        newline_loc--;
                string::size_type text_loc = newline_loc;
                while((data[text_loc] == ' ') ||(data[text_loc] == '    ')||(data[text_loc] == '\n')||(data[text_loc] == '\r')){
                        text_loc++;
                        if(text_loc+1 == data.length())
                                break;
                }
                data = (data.substr(0, newline_loc) + " " + data.substr(text_loc));
                newline_loc = data.find_first_of('\n', 0);
        }
        //now reformat with newlines and whitespace, observing word boundaries...
        string output_string = "";
        for(uint32 charI = 0; charI < data.length();){
                //get the substring to append and increment charI
                string::size_type base_loc = charI;
                string append_string;
                while(base_loc - charI <= width){
                        string::size_type space_loc = data.find_first_of(' ', base_loc+1);
                        if(space_loc - charI < width)
                                base_loc = space_loc;
                        else if(base_loc == charI){
                                //word is too big for one line.  split it.
                                append_string = data.substr(charI, width);
                                charI+=width;
                        }else{
                                append_string = data.substr(charI, base_loc - charI);
                                charI = base_loc;
                        }
                }
                output_string += string(offset, ' ') + append_string;
                if(charI + width < data.length())
                        output_string += "\r\n";
        }
        data = output_string;
}

void gnGBKSource::WriteHeader(gnMultiSpec* spec, const string& hdr, ofstream& m_ofstream){
        gnBaseHeader* gpbh = NULL;
        uint32 header_index = 0;
        do{
                gpbh = spec->GetHeader(hdr, header_index);
                if(gpbh != NULL)
                        m_ofstream << gpbh->GetHeader();
                header_index++;
        }while(gpbh != NULL);
}

boolean gnGBKSource::Write(gnSequence& seq, const string& filename){
        ofstream m_ofstream(filename.c_str(), ios::out | ios::binary);
        if(!m_ofstream.is_open())
                return false;

        string newline = "\r\n";
        gnGenomeSpec* spec = seq.GetSpec();

        // output general file header first if one exists.
        if(spec->GetHeaderListLength() == 1){
                gnBaseHeader *gpbh = spec->GetHeader(0);
                string name = gpbh->GetHeaderName();
                //IMPLEMENT ME!  Is platform specific newline substitution necessary?
                if(string::npos != name.find(".SEQ")){
                        string header = gpbh->GetHeader();
                        m_ofstream << header;
                }
        }

        Array<gnSeqC> array_buf( BUFFER_SIZE );
        gnSeqC *bases = array_buf.data;

        for(uint32 specI = 0; specI < spec->GetSpecListLength(); specI++){
                gnFragmentSpec* subSpec = spec->GetSpec(specI);
                
                //write out contig headers.  start with LOCUS...
                m_ofstream << "LOCUS       ";
                //write Locus Name
                string contigName = subSpec->GetName();
                if(contigName.length() > SEQ_LOCUS_NAME_LENGTH)
                        contigName = contigName.substr(0, SEQ_LOCUS_NAME_LENGTH);
                uint32 filler_size = SEQ_LOCUS_NAME_LENGTH - contigName.length();
                m_ofstream << contigName << string(filler_size, ' ');
                //write Locus Length
                string length_string = uintToString(subSpec->GetLength());
                filler_size = SEQ_LOCUS_SIZE_LENGTH - length_string.size();
                m_ofstream << string(filler_size, ' ') << length_string << " bp ";
                //write dnatype
                string dnatype = string(SEQ_LOCUS_DNATYPE_LENGTH, ' ');
                uint32 head_look_i = 0;
                gnBaseHeader* gpbh = subSpec->GetHeader("LOCUS", head_look_i);
                if(gpbh != NULL)
                        dnatype = gpbh->GetHeader().substr(SEQ_LOCUS_DNATYPE_OFFSET, SEQ_LOCUS_DNATYPE_LENGTH);
                m_ofstream << dnatype << string(2, ' ');
                //write circularity
                string circular = subSpec->IsCircular() ? string("circular  ") : string(10, ' ');
                m_ofstream << circular;
                //write division code
                string division = string(SEQ_LOCUS_DIVCODE_LENGTH, ' ');
                if(gpbh != NULL)
                        division = gpbh->GetHeader().substr(SEQ_LOCUS_DIVCODE_OFFSET, SEQ_LOCUS_DIVCODE_LENGTH);
                m_ofstream << division;
                //write date -- IMPLEMENT ME!  format the real date to spec! dd-mmm-yyyy
                string date = string(SEQ_LOCUS_DATE_LENGTH, ' ');
                if(gpbh != NULL)
                        date = gpbh->GetHeader().substr(SEQ_LOCUS_DATE_OFFSET, SEQ_LOCUS_DATE_LENGTH);
                m_ofstream << string(7, ' ') << date << "\r\n";
                
                //write out the rest of the headers if they were supplied!
                WriteHeader(subSpec, "DEFINITION", m_ofstream);
                WriteHeader(subSpec, "ACCESSION", m_ofstream);
                WriteHeader(subSpec, "VERSION", m_ofstream);
                WriteHeader(subSpec, "KEYWORDS", m_ofstream);
                WriteHeader(subSpec, "SEGMENT", m_ofstream);
                WriteHeader(subSpec, "SOURCE", m_ofstream);
                WriteHeader(subSpec, "REFERENCE", m_ofstream);
                WriteHeader(subSpec, "COMMENT", m_ofstream);

                //write out feature table!
                m_ofstream << "FEATURES             Location/Qualifiers" << "\r\n";
                for(uint32 featureI = 0; featureI < subSpec->GetFeatureListLength(); featureI++){
                        //write a feature tag
                        gnBaseFeature *gpmf = subSpec->GetFeature(featureI);
                        string featureName = gpmf->GetName();
                        m_ofstream << string(SEQ_SUBTAG_COLUMN, ' ') << featureName;
                        m_ofstream << string(SEQ_FEATURE_LOC_OFFSET - featureName.length() - SEQ_SUBTAG_COLUMN, ' ');
                        //ready to output location b.s.
                        uint32 location_count = gpmf->GetLocationListLength();
                        uint32 line_pos = SEQ_FEATURE_LOC_OFFSET;
                        uint32 parenthesis_count = 0;
                        if(location_count > 1){
                                m_ofstream << "join(";
                                line_pos += 5;
                                parenthesis_count++;
                        }
                        gnLocation::gnLocationType loc_type = gpmf->GetLocationType();
                        switch(loc_type){
                                case gnLocation::LT_Standard:
                                        break;
                                case gnLocation::LT_Complement:
                                        m_ofstream << "complement(";
                                        line_pos += 11;
                                        parenthesis_count++;
                                        break;
                                case gnLocation::LT_Order:
                                        m_ofstream << "order(";
                                        line_pos += 6;
                                        parenthesis_count++;
                                        break;
                                case gnLocation::LT_Group:
                                        m_ofstream << "group(";
                                        parenthesis_count++;
                                        line_pos += 6;
                                        break;
                                case gnLocation::LT_OneOf:
                                        m_ofstream << "one-of(";
                                        parenthesis_count++;
                                        line_pos += 7;
                                        break;                          
                                default:
                                        break;
                        }
                        //create the location string, then see if it will fit on the line
                        string location;
                        for(uint32 locationI = 0; locationI < location_count; locationI++){
                                gnLocation gpl = gpmf->GetLocation(locationI);
                                if(gpl.IsStartBoundLonger())
                                        location += ">";
                                if(gpl.IsStartBoundShorter())
                                        location += "<";
                                location += uintToString(gpl.GetStart());
                                gnSeqI end_loc = gpl.GetEnd();
                                if(end_loc != 0){
                                        switch(gpl.GetType()){
                                                case gnLocation::LT_BetweenBases:
                                                        location += "^";
                                                        break;
                                                case gnLocation::LT_OneOf:
                                                        location += ".";
                                                        break;
                                                default:
                                                        location += "..";
                                                        break;
                                        }
                                        if(gpl.IsEndBoundShorter())
                                                location += "<";
                                        if(gpl.IsEndBoundLonger())
                                                location += ">";
                                        location+= uintToString(end_loc);
                                }
                                if(locationI +1 < location_count)
                                        location += ",";
                                else{   //append necessary parenthesis
                                        for(;parenthesis_count > 0; parenthesis_count--)
                                                location += ")";
                                }
                                //put it on this line if it fits.  otherwise make a new line.
                                if(line_pos + location.length() < SEQ_COLUMN_WIDTH){
                                        m_ofstream << location;
                                        line_pos += location.length();
                                }else{
                                        m_ofstream << "\r\n" << string(SEQ_FEATURE_LOC_OFFSET, ' ') << location;
                                        line_pos = SEQ_FEATURE_LOC_OFFSET + location.length();
                                }
                                location = "";
                        }
                        m_ofstream << "\r\n";
                        //now output qualifiers!  yaay!
                        //god damn this is a big ugly piece of code.
                        
                        uint32 qualifier_count = gpmf->GetQualifierListLength();
                        for(uint32 qualifierI = 0; qualifierI < qualifier_count; qualifierI++){
                                m_ofstream << string(SEQ_FEATURE_LOC_OFFSET, ' ');
                                gnBaseQualifier* qualifier = gpmf->GetQualifier(qualifierI);
                                m_ofstream << "/" << qualifier->GetName() << "=";
                                //IMPLEMENT ME! do a better word wrap on this bitch.
                                string qually = string(qualifier->GetValue());
//                              FormatString(qually, SEQ_FEATURE_LOC_OFFSET, 80 - SEQ_FEATURE_LOC_OFFSET);
//                              qually = qually.substr(SEQ_FEATURE_LOC_OFFSET);
                                m_ofstream << qually << "\r\n";
                        }
                        if(gpmf != NULL)
                                delete gpmf;
                }
                
                //get information about the sequence we're writing out.
                gnSeqI readOffset = seq.contigStart(specI);
                gnSeqI readLength = seq.contigLength(specI);

                //finally - output base count and origin
                m_ofstream << "BASE COUNT ";
                gnSeqI a_count=0, c_count=0, g_count=0, t_count=0, other_count=0;
                gnSeqI countLen = readLength + readOffset;
                for(gnSeqI countI = readOffset; countI < countLen;){
                        gnSeqI writeLen = countLen - countI < BUFFER_SIZE ? countLen - countI : BUFFER_SIZE;
                        if(!seq.ToArray(bases, writeLen, countI))
                                return false;
                        gnSeqI a, c, g, t, other;
                        BaseCount(string(bases, writeLen), a, c, g, t, other);
                        a_count += a;
                        c_count += c;
                        g_count += g;
                        t_count += t;
                        other_count += other;
                        countI += writeLen;
                }
                m_ofstream << uintToString(a_count) << " a ";
                m_ofstream << uintToString(c_count) << " c ";
                m_ofstream << uintToString(g_count) << " g ";
                m_ofstream << uintToString(t_count) << " t ";
                m_ofstream << uintToString(other_count) << " others" << "\r\n";

                string origin = "ORIGIN\r\n";
                head_look_i = 0;
                gpbh = subSpec->GetHeader("ORIGIN", head_look_i);
                if(gpbh != NULL)
                        origin = gpbh->GetHeader();
                m_ofstream << origin;
                
                //write out the sequence
                gnSeqI contig_bases = 0;
                while(readLength > 0){  //buffer the read/writes
                        gnSeqI writeLen = readLength < BUFFER_SIZE + 20 ? readLength : BUFFER_SIZE + 20;
                        boolean success = seq.ToArray(bases, writeLen, readOffset);
                        if(!success)
                                return false;
                        //print each 60 on their own lines...
                        for(gnSeqI curbaseI = 0; curbaseI < writeLen; curbaseI += 60){
                                string baseIndexStr = uintToString(contig_bases + curbaseI +1);
                                m_ofstream << string(SEQ_BASES_INDEX_END - baseIndexStr.length(), ' ');
                                m_ofstream << baseIndexStr;
                                for(gnSeqI base_offset = 0; base_offset <= 50; base_offset+=10){
                                        if(writeLen <= curbaseI + base_offset)
                                                break;
                                        int64 print_length = writeLen - (curbaseI + base_offset);
                                        print_length = print_length > 10 ? 10 : print_length;
                                        m_ofstream << ' ' << string(bases + curbaseI + base_offset, print_length);
                                }
                                m_ofstream << "\r\n";
                        }
                        readLength -= writeLen;
                        readOffset += writeLen;
                        contig_bases += writeLen;
                }
                m_ofstream << "//\r\n";
        }
        
        m_ofstream.close();
        return true;
}

gnFileContig* gnGBKSource::GetFileContig( const uint32 contigI ) const{
        if(m_contigList.size() > contigI)
                return m_contigList[contigI];
        return NULL;
}

//File parsing access routine
boolean gnGBKSource::ParseStream( istream& fin )
{
        // INIT temp varables
        uint32 readState = 0;
        uint32 lineStart = 0;
        // INIT buffer
        uint32 sectionStart = 0;
        uint64 streamPos = 0;
        uint64 bufReadLen = 0;
        uint64 remainingBuffer = 0;
        Array<char> array_buf( BUFFER_SIZE );
        char* buf = array_buf.data;
        gnFragmentSpec* curFrag = 0;
        gnSourceSpec* curSpec = 0;
        gnSourceHeader *curHeader;
        gnFeature* curFeature;
        gnFileContig* curContig = 0;
        gnLocation::gnLocationType curBaseLocationType;
        gnSeqI curLocationStart;
        int32 curStartLength = 0;
        int32 curEndLength = 0;
        string curLocContig = "";
        string curQualifierName;
        uint64 curQualifierStart;
        string curContigName = "";
        gnSeqI seqLength = 0;
        gnSeqI seqChunk, seqChunkCount, gapChunk;
        boolean corruptWarning = false;
        
        //decide what type of newlines we have
        DetermineNewlineType();

        m_spec = new gnGenomeSpec();
        while( !fin.eof() )
        {
                if(sectionStart > 0){
                        if(readState == 14)
                                sectionStart = lineStart;
                        remainingBuffer = bufReadLen - sectionStart;
                        memmove(buf, buf+sectionStart, remainingBuffer);
                }
                  // read chars
                fin.read( buf + remainingBuffer, BUFFER_SIZE - remainingBuffer);
                streamPos -= remainingBuffer;
                lineStart -= sectionStart;
                sectionStart = 0;
                bufReadLen = fin.gcount();
                bufReadLen += remainingBuffer;
                
                for( uint32 i=remainingBuffer ; i < bufReadLen ; i++ )
                {
                        char ch = buf[i];
                        switch( readState )
                        {
                                case 0:         //Assume we are in header at the start of a new line.  
                                                        //Look for keywords starting in column 1
                                        if((ch == '\n')&&(buf[lineStart] != ' ')&&(buf[lineStart] != '  ')){  //not equal to space or tab
                                                if(curSpec == NULL){
                                                        curSpec = new gnSourceSpec(this, m_spec->GetSpecListLength());
                                                        curFrag = new gnFragmentSpec();
                                                        curFrag->AddSpec(curSpec);
                                                        curSpec->SetSourceName(m_openString);
                                                        m_spec->AddSpec(curFrag);
                                                }
                                                if(lineStart != sectionStart){  //Add the previous header to our list
                                                        uint32 j = SEQ_HEADER_NAME_LENGTH-1;
                                                        for(; j > 0; j--)       
                                                                if((buf[sectionStart+j] != ' ')&&(buf[sectionStart+j] != '      '))
                                                                        break;
                                                        string header_name = string(buf+sectionStart, j+1);
                                                        curHeader = new gnSourceHeader(this, header_name, sectionStart + streamPos, lineStart - sectionStart);
                                                        //if this is header info _before_ a locus statement then its a general file header.
                                                        if(strncmp(&buf[lineStart], "LOCUS", 5) == 0)
                                                                m_spec->AddHeader(curHeader);
                                                        else    //otherwise its a fragment header.
                                                                curFrag->AddHeader(curHeader);
                                                        sectionStart = lineStart;
                                                }
                                                
                                                if(strncmp(&buf[lineStart], "FEATURES", 8) == 0){
                                                        sectionStart = i + 1;
                                                        readState = 1;  //read in features
                                                }else if(strncmp(&buf[lineStart], "ORIGIN", 6) == 0){
                                                        curHeader = new gnSourceHeader(this, string("ORIGIN"), sectionStart + streamPos, i - sectionStart + 1);
                                                        curFrag->AddHeader(curHeader);
                                                        curContig = new gnFileContig();
                                                        curContig->SetName(curContigName);
                                                        curContigName = "";
                                                        readState = 13;  //read in base pairs
                                                }else if(strncmp(&buf[lineStart], "LOCUS", 5) == 0){
                                                        if(strncmp(&buf[lineStart+SEQ_LOCUS_CIRCULAR_COLUMN-1], "circular", 8) == 0)
                                                                curFrag->SetCircular(true);
                                                        uint32 j = SEQ_LOCUS_NAME_LENGTH;
                                                        for(; j >= 0; j--)      
                                                                if((buf[lineStart+SEQ_LOCUS_NAME_COLUMN+j-1] != ' ')&&(buf[sectionStart+SEQ_LOCUS_NAME_COLUMN+j-1] != ' '))
                                                                        break;
                                                        curContigName = string(buf+lineStart+SEQ_LOCUS_NAME_COLUMN-1, j+1);
                                                        curFrag->SetName(curContigName);
                                                }
                                        }
                                        if(ch == '\n'){
                                                lineStart = i + 1;
                                        }
                                        break;
                                case 1: //look for feature tag in column six.  ignore whitespace before feature.
                                        if((ch == ' ')||(ch == '        ')){
                                                break;
                                        }else if(ch == '\n'){
                                                lineStart = i + 1;
                                                sectionStart = i + 1;
                                                break;
                                        }else if(sectionStart == i){ //there was no whitespace, we hit a TAG instead
                                                i--;
                                                readState = 0; //Deal with a Header TAG
                                                sectionStart = i + 1;
                                                break;
                                        }else if((i - lineStart == SEQ_SUBTAG_COLUMN)||((buf[lineStart]=='      ')&&(i==lineStart+1))){
                                                sectionStart = i;
                                                readState = 2;
                                        } //
                                case 2:  //Get the feature name.  stop on whitespace
                                        if((ch == ' ')||(ch == '        ')){
                                                string featureName(buf+sectionStart, i - sectionStart);
                                                curFeature = new gnFeature(featureName);
                                                curFrag->AddFeature(curFeature);
                                                sectionStart = i + 1;
                                                readState = 3;
                                        }
                                        break;
                                case 3:   //Ignore whitespace before feature location
                                        if((ch == ' ')||(ch == '        ')){
                                                break;
                                        }else if((ch == '\r')||(ch == '\n')){
                                                lineStart = i+1;
                                                break;
                                        }
                                        sectionStart = i;
                                        readState = 4;
                                // KNOWN BUG HERE!
                                // if JOIN is outside a group of complemented coordinates the feature type
                                // is incorrectly set to LT_COMPLEMENT.  Instead each location's type should be set to LT_COMPLEMENT 
                                case 4:         //Read a location start.  stop on (<.:^ and whitespace
                                        if((ch == ' ')||(ch == '        ')||(ch == '(')||(ch == '.')||(ch=='^')||(ch==':')){
                                                string starter(buf+sectionStart, i - sectionStart);
                                                if(ch == '('){
                                                        if(starter == "complement")
                                                                curFeature->SetLocationType(gnLocation::LT_Complement);
                                                        else if(starter == "order")
                                                                curFeature->SetLocationType(gnLocation::LT_Order);
                                                        else if(starter == "group")
                                                                curFeature->SetLocationType(gnLocation::LT_Group);
                                                        else if(starter == "one-of")
                                                                curFeature->SetLocationType(gnLocation::LT_OneOf);
                                                        sectionStart = i + 1;   //ignore join since it is default.
                                                        break;
                                                }else if(ch == ':'){
                                                        curLocContig = starter;
                                                        sectionStart = i + 1;
                                                        break;
                                                }
                                                curLocationStart = atoi(starter.c_str());
                                                readState = 6;  //read in end base by default.
                                                if(ch == '.'){
                                                        //go to special state to look for another one.
                                                        readState = 5;
                                                        sectionStart = i + 1;
                                                        break;
                                                }else if(ch == '^'){
                                                        curBaseLocationType = gnLocation::LT_BetweenBases;
                                                }else if((ch == ' ')||(ch == '  ')){
                                                        //no end location go to qualifier
                                                        gnLocation curLocation(curLocationStart, curLocationStart);
                                                        curFeature->AddLocation(curLocation, curFeature->GetLocationListLength());
                                                        readState = 7;
                                                }
                                                sectionStart = i + 1;

                                        }else if(ch == '<'){
                                                curStartLength = -1;
                                                sectionStart = i + 1;
                                        }else if(ch == '>'){
                                                curStartLength = 1;
                                                sectionStart = i + 1;
                                        }
                                        break;
                                case 5: //look for another period or location start.
                                        if(ch == '.'){
                                                curBaseLocationType = gnLocation::LT_Standard;
                                                readState = 6;
                                                sectionStart = i + 1;
                                                break;
                                        }
                                        curBaseLocationType = gnLocation::LT_OneOf;
                                case 6: //see if there's a second location value.  stop on >, and whitespace
                                        if(ch == '>'){
                                                curEndLength = 1;
                                                sectionStart = i + 1;
                                        }else if(ch == '<'){
                                                curEndLength = -1;
                                                sectionStart = i + 1;
                                        }else if((ch == ' ')||(ch == '  ')||(ch == ',')){
                                                //read end location
                                                string ender(buf+sectionStart, i - sectionStart);
                                                gnSeqI curLocationEnd = atoi(ender.c_str());
                                                gnLocation curLocation(curLocationStart, curStartLength, curLocationEnd, curEndLength, curBaseLocationType);
                                                curEndLength = 0;
                                                curStartLength = 0;
                                                curFeature->AddLocation(curLocation, curFeature->GetLocationListLength());
                                                readState = ch == ',' ? 3 : 7;  //read another loc if we need to.
                                                sectionStart = i+1;
                                        }
                                        break;
                                case 7:  //skip to start of qualifier
                                        if((ch != ' ')&&(ch != '        ')&&(lineStart == i)){
                                                sectionStart = i;       // Hit a tag.  go deal with it.
                                                readState = 0;
                                                i--;
                                        }else if((ch != ' ')&&(ch != '  ')&&((lineStart == i - SEQ_SUBTAG_COLUMN)||((buf[lineStart]=='  ')&&(i==lineStart+1)))){
                                                sectionStart = i;       // Hit a feature.  go deal with it.
                                                readState = 2;
                                                i--;
                                        }else if(ch == ','){  //oops!  another location to read!
                                                sectionStart = i+1;
                                                readState = 3;
                                        }else if(ch == '/'){  //finally, a qualifier.
                                                sectionStart = i+1;
                                                readState = 8;
                                        }else if(ch == '\n')
                                                lineStart = i + 1;
                                        break;
                                case 8:         //get a qualifier, stop on =
                                        if(ch == '='){
                                                curQualifierName = string(buf+sectionStart, i - sectionStart);
                                                readState = 9;
                                                sectionStart = i+1;
                                        }
                                        break;
                                case 9:         //are we getting a string? look for " or [
                                        if(ch == '"'){
                                                readState = 10;
                                                sectionStart = i;
                                                curQualifierStart = i + streamPos;
                                        }else if(ch == '['){
                                                readState = 11;
                                                sectionStart = i;
                                        }else if((ch == '\r')||(ch == '\n')){
                                                curFeature->AddQualifier(new gnSourceQualifier(this, curQualifierName, sectionStart + streamPos, i - sectionStart));
                                                sectionStart = i+1;
                                                readState = 7; //look for another qualifier
                                        }
                                        break;
                                case 10:                //read until the end of the quotation. look out for escaped quotes
                                        if(ch == '"')
                                                readState = 11;
                                        if(ch == '\n'){
                                                lineStart = i + 1;
                                        }
                                        break;
                                case 11:
                                        if(ch != '"'){
                                                curFeature->AddQualifier(new gnSourceQualifier(this, curQualifierName, curQualifierStart, i - sectionStart));
                                                sectionStart = i+1;
                                                readState = 7;  //look for another qualifier.
                                                if(ch == '\n')
                                                        lineStart = i + 1;
                                        }else
                                                readState = 10;  //quote was escaped.  look for another.
                                        break;
                                case 12:
                                        if(ch == ']'){
                                                curFeature->AddQualifier(new gnSourceQualifier(this, curQualifierName, sectionStart + streamPos, i - sectionStart));
                                                sectionStart = i+1;
                                                readState = 7;  //look for another qualifier.
                                        }
                                        break;
                                case 13:        //start the sequence read.
                                        curContig->SetSectStart(gnContigSequence, i - 1 + streamPos);
                                        curContig->SetRepeatSeqGap(true);
                                        seqChunk = 0;
                                        seqChunkCount = 0;
                                        gapChunk = m_newlineSize + 1;
                                        readState = 14;
                                        break;
                                case 14:
                                        while(i < bufReadLen){
                                                ch = buf[i];
                                                if((ch == '/')&&(i==lineStart)){
                                                        readState = 15; //end of this sequence
                                                        break;
                                                }else if(m_pFilter->IsValid(ch)){
                                                        if(gapChunk > 0){
                                                                if((gapChunk > 1 && seqChunkCount > 0) ||
                                                                  (gapChunk != 10 + m_newlineSize && seqChunkCount == 0)){
                                                                        if( !corruptWarning ){
                                                                                ErrorMsg("File is corrupt.  Proceed with caution.");
                                                                                corruptWarning = true;
                                                                        }
                                                                        curContig->SetRepeatSeqGap(false);
                                                                }
                                                                gapChunk = 0;
                                                        }
                                                        seqChunk++;
                                                        seqLength++;
                                                }else{
                                                        gapChunk++;
                                                        if(seqChunk == 10){
                                                                seqChunk = 0;
                                                                seqChunkCount++;
                                                                if(seqChunkCount == 6){
                                                                        //got a complete line.  start over
                                                                        seqChunkCount = 0;
                                                                }
                                                        }
                                                        if(ch == '\n')
                                                                lineStart = i + 1;
                                                }
                                                i++;
                                        }
                                        break;
                                case 15:
                                        if((ch == '\n')&&(buf[lineStart+1] == '/')){
                                                curContig->SetSectEnd(gnContigSequence, lineStart - 2 + streamPos);
                                                curContig->SetSeqLength(seqLength);
                                                m_contigList.push_back(curContig);
                                                curContig = 0;
                                                curSpec->SetLength(seqLength);
                                                curSpec = 0;
                                                seqLength = 0;
                                                lineStart = i + 1;
                                                sectionStart = i + 1;
                                                readState = 0;
                                        }
                                        break;
                        }
                }
                streamPos += bufReadLen;
        }
        if(curContig != 0){
                curContig->SetSectEnd(gnContigSequence, streamPos - 1);
                curContig->SetSeqLength(seqLength);
                m_contigList.push_back(curContig);
                curSpec->SetLength(seqLength);
        }
        if(curSpec != 0)
                if((curFrag->GetFeatureListLength() == 0) && (curFrag->GetHeaderListLength() == 0)
                        &&(curSpec->GetLength() == 0)){
                        m_spec->RemoveSpec(m_spec->GetSpecListLength() - 1);
                        delete curFrag;
                }
        m_ifstream.clear();
        return true;
}

---