---

gnSequence.cpp

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
// File:            gnSequence.h
// Purpose:         Source Seq, where sequence is stored in source file
// Description:     implements the baseSeq interface for dna source files.
// Changes:        
// Version:         libGenome 0.5.1 
// Author:          Aaron Darling 
// Modified by:     
// Copyright:       (c) Aaron Darling 
// Licenses:        See COPYING file for details
/////////////////////////////////////////////////////////////////////////////
#include "gn/gnSequence.h"
#include "gn/gnSourceFactory.h"
#include "gn/gnBaseSource.h"
#include "gn/gnGenomeSpec.h"
#include "gn/gnFragmentSpec.h"
#include "gn/gnSourceSpec.h"
#include "gn/gnStringSpec.h"
#include "gn/gnStringHeader.h"

// Constructor
gnSequence::gnSequence( )
{
        spec = new gnGenomeSpec();
        comparator = gnCompare::DNASeqCompare();
}
gnSequence::gnSequence( const gnSeqC* seq){
        spec = new gnGenomeSpec();
        if (seq[0] != 0){
                gnFragmentSpec *fragSpec = new gnFragmentSpec();
                spec->AddSpec(fragSpec);
                fragSpec->AddSpec(new gnStringSpec(seq));
        }
        comparator = gnCompare::DNASeqCompare();
}
gnSequence::gnSequence( const string& str){
        spec = new gnGenomeSpec();
        if (str.length() != 0){
                gnFragmentSpec *fragSpec = new gnFragmentSpec();
                spec->AddSpec(fragSpec);
                fragSpec->AddSpec(new gnStringSpec(str));
        }
        comparator = gnCompare::DNASeqCompare();
}
gnSequence::gnSequence( const gnGenomeSpec& gngs){
        spec = gngs.Clone();
        comparator = gnCompare::DNASeqCompare();
}
gnSequence::gnSequence( const gnFragmentSpec& gnfs){
        spec = new gnGenomeSpec();
        spec->AddSpec(gnfs.Clone());
        comparator = gnCompare::DNASeqCompare();
}
gnSequence::gnSequence( const gnContigSpec& gcs){
        spec = new gnGenomeSpec();
        gnFragmentSpec *fragSpec = new gnFragmentSpec();
        fragSpec->AddSpec(gcs.Clone());
        comparator = gnCompare::DNASeqCompare();
}
gnSequence::gnSequence( gnSeqC *bases, gnSeqI length){
        spec = new gnGenomeSpec();
        if (length > 0){
                gnFragmentSpec *fragSpec = new gnFragmentSpec();
                spec->AddSpec(fragSpec);
                fragSpec->AddSpec(new gnStringSpec(string(bases, length)));
        }
        comparator = gnCompare::DNASeqCompare();
}       
gnSequence::gnSequence( const gnSequence& seq )
{
        *this = seq;
}

// Destructor
gnSequence::~gnSequence()
{
        if (spec != NULL) 
                delete spec;
}

gnSequence& gnSequence::operator=( const gnSequence& seq){
        spec = seq.spec->Clone();
        filter_list = seq.filter_list;
        comparator = seq.comparator;
        return *this;
}

// Clone
gnSequence* gnSequence::Clone() const
{
        return new gnSequence( *this );
}
//IMPLEMENT ME!! Consider ambiguities!
int gnSequence::compare(const string& str) const{
        STACK_TRACE_START
                gnSeqI len = length();
                gnSeqI seq_len = str.length();
                gnSeqI comparelen = len < seq_len ? len : seq_len;
                gnSeqI compared = 0;
                while(comparelen > 0){
                        gnSeqI curlen = comparelen < BUFFER_SIZE ? comparelen : BUFFER_SIZE;
                        string bases = ToString(compared, curlen);
                        string seq_bases = str.substr(compared, curlen);
                        if(comparator->LessThan(bases, seq_bases))
                                return -1;
                        else if(comparator->LessThan(seq_bases, bases))
                                return 1;

                        compared += curlen;
                        comparelen -= curlen;
                }
                if(len < seq_len)
                        return -1;
                else if(len > seq_len)
                        return 1;
                return 0;
        STACK_TRACE_END
}

int gnSequence::compare(const gnSequence& seq) const{
        STACK_TRACE_START
                gnSeqI len = length();
                gnSeqI seq_len = seq.length();
                gnSeqI comparelen = len < seq_len ? len : seq_len;
                gnSeqI compared = 0;
                while(comparelen > 0){
                        gnSeqI curlen = comparelen < BUFFER_SIZE ? comparelen : BUFFER_SIZE;
                        string bases = ToString(curlen, compared+1);
                        string seq_bases = seq.ToString(curlen, compared+1);
                        if(comparator->LessThan(bases, seq_bases))
                                return -1;
                        else if(comparator->LessThan(seq_bases, bases))
                                return 1;

                        compared+= curlen;
                        comparelen -= curlen;
                }
                if(len < seq_len)
                        return -1;
                else if(len > seq_len)
                        return 1;
                return 0;
        STACK_TRACE_END
}

void gnSequence::insert( const gnSeqI offset, const gnSeqC *bases, const gnSeqI len){
        STACK_TRACE_START
                string str(bases, len);
                gnStringSpec gpbs(str);
                insert(offset, gpbs);
        STACK_TRACE_END
}
//Offset is the gene index, not a computer index, to insert before.
//This is the default insert.  It inserts entire contigs
void gnSequence::insert( const gnSeqI offset, const gnGenomeSpec& gpbs){
        STACK_TRACE_START
                if(offset == 0)
                        Throw_gnEx(SeqIndexOutOfBounds());
                if(offset == GNSEQI_END || spec->GetLength() < offset){ //simple append.
                        for(uint32 i=0; i < gpbs.GetSpecListLength(); i++)
                                spec->AddSpec(gpbs.GetSpec(i)->Clone());
                        return;
                }
                //clone this sequence
                gnGenomeSpec* tmpSpec = spec->Clone();
                //crop off the end of this sequence past the insert point
                //crop off the beginning of the cloned sequence up to the insert
                gnSeqI real_offset = offset - 1;
                gnSeqI endCrop = spec->GetLength() - real_offset;
                spec->CropEnd(endCrop);
                tmpSpec->CropStart(real_offset);
                //append the inserted sequence and the end of this sequence.
                insert(GNSEQI_END, gpbs);
                insert(GNSEQI_END, *tmpSpec);
                delete tmpSpec;
        STACK_TRACE_END
}
// Concatenation operators
gnSequence const gnSequence::operator+(const gnSequence& seq) const{
        STACK_TRACE_START
                gnSequence rval(*this);
                rval.append(seq);
                return rval;
        STACK_TRACE_END
}
// Subsequences and base deletion
gnSequence gnSequence::subseq(const gnSeqI offset, const gnSeqI length) const{
        STACK_TRACE_START
                if(length == 0)
                        return gnSequence();
                if(offset == 0)
                        Throw_gnEx(SeqIndexOutOfBounds());
                
                gnSequence tmpSeq;
                delete tmpSeq.spec;
                tmpSeq.spec = spec->CloneRange(offset - 1, length);
                return tmpSeq;
        STACK_TRACE_END
}
void gnSequence::erase( const gnSeqI offset, const gnSeqI len ){
        STACK_TRACE_START
                if(offset == 0)
                        Throw_gnEx(SeqIndexOutOfBounds());
                //range checking, etc.
                gnSeqI current_length = length();
                if(offset > current_length)
                        Throw_gnEx(SeqIndexOutOfBounds());
                gnSeqI endBase = offset + len - 1 < current_length ? offset + len - 1 : current_length;
                gnSeqI startBase = offset - 1;

                gnGenomeSpec* tmpSpec = spec->CloneRange(endBase, GNSEQI_END);
                uint32 lennard = tmpSpec->GetLength();
                spec->CropEnd(current_length - startBase);
                lennard = length();
                insert(GNSEQI_END, *tmpSpec);
                delete tmpSpec;
        STACK_TRACE_END
}
void gnSequence::splitContig(const gnSeqI splitI, const uint32 contigI){
        STACK_TRACE_START
                gnSeqI splitBase = splitI;
                gnSeqI current_length = length();
                if(splitI == 0)
                        Throw_gnEx(SeqIndexOutOfBounds());
                if(contigI == ALL_CONTIGS && splitI > current_length)
                                Throw_gnEx(SeqIndexOutOfBounds());
                else
                        localToGlobal(contigI, splitBase);
                gnGenomeSpec* tmpSpec = spec->Clone();
                tmpSpec->CropStart(splitBase);
                spec->CropEnd(current_length - splitBase);
                
                insert(GNSEQI_END, *tmpSpec);
                delete tmpSpec;
        STACK_TRACE_END
}

// IO operators
istream& operator>>(istream& is, gnSequence& gps){      //read from source.
        string bases;
        is >> bases;
        gps.append(bases);
        return is;
}
ostream& operator<<(ostream& os, const gnSequence& gps){ //write to source.
        os << gps.ToString();
        return os;
}
string gnSequence::ToString( const gnSeqI len, const gnSeqI offset ) const
{
        STACK_TRACE_START
                string str;
                ToString(str, len, offset);
                return str;
        STACK_TRACE_END
}
boolean gnSequence::ToString( string& str, const gnSeqI len, const gnSeqI offset ) const
{
        STACK_TRACE_START
                gnSeqI real_offset = offset - 1;
                gnSeqI m_length = length();
                gnSeqI readSize = len > m_length - real_offset ? m_length - real_offset : len;
                Array<char> array_buf( readSize+1 );
                char *buf = array_buf.data;
                boolean success = spec->SeqRead( real_offset, buf, readSize, ALL_CONTIGS);
                buf[readSize] = '\0';
                str = buf;

                //now filter the string
                list<const gnBaseFilter*>::const_iterator iter = filter_list.begin();
                for(; iter != filter_list.end(); iter++)
                        (*iter)->Filter(str);

                if( success )
                        return true;
                return false;
        STACK_TRACE_END
}
boolean gnSequence::ToArray( gnSeqC* pSeqC, gnSeqI length, const gnSeqI offset ) const
{
        STACK_TRACE_START
                gnSeqI real_offset = offset - 1;
                if(offset == GNSEQI_END)
                        return false;
                gnSeqC* tmp = new gnSeqC[length];
                boolean success = spec->SeqRead( real_offset, tmp, length, ALL_CONTIGS);
                
                //now filter the array
                list<const gnBaseFilter*>::const_iterator iter = filter_list.begin();
                gnSeqC** tomp = &tmp;
                for(; iter != filter_list.end(); iter++)
                        (*iter)->Filter(tomp, length);
                memcpy(pSeqC, *tomp, length);
                delete[] *tomp;
                
                if( success )
                        return true;

                return false;
        STACK_TRACE_END
}
gnSeqC gnSequence::GetSeqC( const gnSeqI offset ) const
{
        STACK_TRACE_START
                char block;
                gnSeqI readLen = 1;
                gnSeqI real_offset = offset - 1;
                boolean success = spec->SeqRead( real_offset, &block, readLen, ALL_CONTIGS);

                //now filter the char
                list<const gnBaseFilter*>::const_iterator iter = filter_list.begin();
                for(; iter != filter_list.end(); iter++)
                        block = (*iter)->Filter(block);

                if( success )
                        return block;

                return GNSEQC_NULL;
        STACK_TRACE_END
}
gnSeqC gnSequence::operator[]( const gnSeqI offset ) const
{
        STACK_TRACE_START
                return GetSeqC( offset );
        STACK_TRACE_END
}

gnSeqI gnSequence::contigListSize() const{
        STACK_TRACE_START
                return spec->GetSpecListLength();
        STACK_TRACE_END
}
gnSeqI gnSequence::contigListLength() const{
        STACK_TRACE_START
                return spec->GetSpecListLength();
        STACK_TRACE_END
}

//find the contig associated with base
uint32 gnSequence::contigIndexByBase( const gnSeqI baseI) const{
        STACK_TRACE_START
                return spec->GetSpecIndexByBase(baseI-1);
        STACK_TRACE_END
}
gnSequence gnSequence::contig( const uint32 contigI) const{
        STACK_TRACE_START
                if(contigI == ALL_CONTIGS)
                        return *this;
                return gnSequence(*spec->GetSpec(contigI));
        STACK_TRACE_END
}
//returns a gnSequence pointer containing the specified contig.
gnSequence gnSequence::contigByBase( const gnSeqI baseI) const{
        STACK_TRACE_START
                return gnSequence(*spec->GetSpecByBase(baseI-1));
        STACK_TRACE_END
}
uint32 gnSequence::contigIndexByName( string& contigName) const{
        STACK_TRACE_START
                return spec->GetSpecIndexByName(contigName);
        STACK_TRACE_END
}

gnSeqI gnSequence::contigStart( const uint32 contigI) const{
        STACK_TRACE_START
                int64 leafI = contigI;
                // add one for the stupid geneticists whose indices start at 1
                return spec->GetSpecStartBase( leafI ) + 1;
        STACK_TRACE_END
}

gnSeqI gnSequence::contigLength( const uint32 contigI) const{
        STACK_TRACE_START
                return spec->GetSpec( contigI )->GetLength();
        STACK_TRACE_END
}

string gnSequence::contigName( const uint32 contigI) const{
        STACK_TRACE_START
                return spec->GetSpec(contigI)->GetName();
        STACK_TRACE_END
}

void gnSequence::globalToLocal(uint32& contigI, gnSeqI& baseI) const{
        STACK_TRACE_START
                contigI = contigIndexByBase(baseI);
                baseI -= (contigStart(contigI) - 1);
        STACK_TRACE_END
}

void gnSequence::localToGlobal(const uint32 contigI, gnSeqI& baseI) const{
        STACK_TRACE_START
                if(baseI > contigLength(contigI))
                        Throw_gnEx(SeqIndexOutOfBounds());
                baseI += contigStart(contigI) - 1;
        STACK_TRACE_END
}

void gnSequence::globalToSource(uint32& contigI, gnSeqI& baseI) const{
        STACK_TRACE_START
                baseI--;        //convert from geneticist coordinates
                gnSeqI seq_contig = spec->GetSpecIndexByBase(baseI);
                gnSeqI new_base = baseI - spec->GetSpecStartBase(seq_contig);
                gnFragmentSpec *fragment_spec = spec->GetSpec(seq_contig);
                seq_contig = fragment_spec->GetSpecIndexByBase(new_base);
                new_base = new_base - fragment_spec->GetSpecStartBase(seq_contig);
                gnContigSpec* contig_spec = fragment_spec->GetSpec(seq_contig);
                contigI = contig_spec->GetSourceContigIndex();
                gnSeqI contig_start = contig_spec->GetStart();
                if(contig_spec->IsReverseComplement()){
                        gnSeqI source_len = contig_spec->GetSourceLength();
                        //it may seem counter intuitive, but we can add the source_len and mod by source_len
                        //to deal with circularity
        //              contig_start = (contig_start - contig_spec->GetLength() + source_len) % source_len;
                        //if we are reverse complement then the contig_start will be the ending base pair
                        //we want to subtract new_base from it. and we can add/mod by source_len to deal
                        //with circles.
                        baseI = (contig_start - new_base - 1 + source_len) % source_len;
                }else
                        baseI =  contig_start + new_base + 1;
        STACK_TRACE_END
}

void gnSequence::localToSource(uint32& contigI, gnSeqI& baseI) const{
        STACK_TRACE_START
                localToGlobal(contigI, baseI);
                globalToSource(contigI, baseI);
        STACK_TRACE_END
}

gnSequence gnSequence::contigByName( string& contigName) const{
        STACK_TRACE_START
                uint32 contigIndex = spec->GetSpecIndexByName(contigName);
                return gnSequence(*spec->GetSpec(contigIndex));
        STACK_TRACE_END
}

void gnSequence::setContigName( const uint32 contigI, const string& contig_name) {
        STACK_TRACE_START
                if(contigI == ALL_CONTIGS)
                        spec->SetName(contig_name);
                else
                        spec->GetSpec(contigI)->SetName(contig_name);
        STACK_TRACE_END
}

void gnSequence::setReverseComplement( const boolean revComp, const uint32 contigI){
        STACK_TRACE_START
                if(contigI == ALL_CONTIGS)
                        spec->SetReverseComplement(revComp);
                else{
                        gnFragmentSpec* contig = spec->GetSpec(contigI);
                        contig->SetReverseComplement(revComp);
                }
        STACK_TRACE_END
}

boolean gnSequence::isReverseComplement( const uint32 contigI ){
        STACK_TRACE_START
                if(contigI == ALL_CONTIGS)
                        return spec->IsReverseComplement();
                gnFragmentSpec* contig = spec->GetSpec(contigI);
                return contig->IsReverseComplement();
        STACK_TRACE_END
}

uint32 gnSequence::getHeaderListLength(const uint32 contigI) const{
        STACK_TRACE_START
                if(contigI == ALL_CONTIGS)
                        return spec->GetHeaderListLength();
                else{
                        return spec->GetSpec(contigI)->GetHeaderListLength();
                }
        STACK_TRACE_END
}

gnBaseHeader* gnSequence::getHeader(const uint32 contigI, const uint32 headerI) const{
        STACK_TRACE_START
                if(contigI == ALL_CONTIGS)
                        return spec->GetHeader(headerI);
                else{
                        return spec->GetSpec(contigI)->GetHeader(headerI);
                }
        STACK_TRACE_END
}

void gnSequence::addHeader(const uint32 contigI, gnBaseHeader* header, const uint32 headerI){
        STACK_TRACE_START
                if(contigI == ALL_CONTIGS)
                        spec->AddHeader(header, headerI);
                else{
                        spec->GetSpec(contigI)->AddHeader(header, headerI);
                }
        STACK_TRACE_END
}

void gnSequence::removeHeader(const uint32 contigI, const uint32 headerI){
        STACK_TRACE_START
                if(contigI == ALL_CONTIGS)
                        spec->RemoveHeader(headerI);
                else
                        spec->GetSpec(contigI)->RemoveHeader(headerI);
        STACK_TRACE_END
}

void gnSequence::merge(const gnSeqI startBase, const gnSeqI endBase){
        STACK_TRACE_START
        STACK_TRACE_END
}

void gnSequence::mergeContigs(const uint32 startC, const uint32 endC){
        STACK_TRACE_START
                spec->MergeFragments(startC, endC);
        STACK_TRACE_END
}

bool gnSequence::LoadSource(const string sourcename){
        STACK_TRACE_START
                gnSourceFactory* m_sourceFactory = gnSourceFactory::GetSourceFactory();
                gnBaseSource *m_pSource = m_sourceFactory->AddSource(sourcename, true);
                if (m_pSource!=NULL)
                {
                        if(spec != NULL)
                                delete spec;
                        spec = m_pSource->GetSpec();
                        return true;
                }
                return false;
        STACK_TRACE_END
}

gnSeqI gnSequence::find(const gnSequence& search, const gnSeqI offset)const{
        STACK_TRACE_START
                //this is really adHoc... should probably be switched
                string searchIn=ToString();
                string find= search.ToString();
                string::size_type pos = searchIn.find(find, offset);
                if (pos == string::npos)
                        return GNSEQI_ERROR;
                else{
                        return pos++; 
                }
        STACK_TRACE_END
}

---