Alignment.tcc

Go to the documentation of this file.

/*

Copyright (C) 2003-2009 Kevin Thornton, krthornt[]@[]uci.edu

Remove the brackets to email me.

This file is part of libsequence.

libsequence is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

libsequence is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
long with libsequence.  If not, see <http://www.gnu.org/licenses/>.

*/

// Code for the -*- C++ -*- namespace Sequence::Alignment
#include <Sequence/Alignment.hpp>
#include <Sequence/SeqConstants.hpp>
#include <Sequence/SeqProperties.hpp>
#include <boost/static_assert.hpp>
#include <boost/type_traits.hpp>
#include <iterator>
#include <algorithm>
namespace Sequence
{
  namespace Alignment
  {
    template < typename T >
    void GetData (std::vector < T >&seqarray,
                  const char *infilename) 
    
    {
      if(!infilename || infilename == NULL)
        return;

      std::ifstream infile(infilename);
      if (infile)
        {
          T temp;
          while (!(infile.eof ()))
            {
              infile >> temp;
              seqarray.push_back(temp);
            }
        }
    }

    template < typename T >
    std::istream & GetData (std::vector <
                            T >&seqarray,
                            std::istream & input_stream) 
    
    {
      T temp;
      if (input_stream)
        {
          while (!(input_stream.eof ()))
            {
              input_stream >> temp;
              seqarray.push_back(temp);
            }
        }
      return input_stream;
    }

    template < typename T >
    std::istream & ReadNObjects ( std::vector < T > &seqarray,
                                  unsigned n,
                                  std::istream & input_stream)
    {
      for(unsigned i = 0 ; i < n ; ++i)
        {
          if (!input_stream || input_stream.eof())
            return input_stream;
          T temp;
          input_stream >> temp;
          seqarray.push_back(temp);
        }
      return input_stream;
    }

    template < typename T >
    void EmptyVector (std::vector< T * > &seqarray)
    {
      for(unsigned i = 0 ; i < seqarray.size() ; ++i)
        delete seqarray[i];
      seqarray.resize(0);
    }

    template < typename T >
    bool Gapped (const std::vector < T >&data)
    {
      BOOST_STATIC_ASSERT((boost::is_base_and_derived<std::pair<std::string,std::string>,T>::value
                           || boost::is_same<std::pair<std::string,std::string>,T>::value));
      for (int i = 0; unsigned (i) < data.size (); ++i)
        //iterate over sequences
        {
          if( data[i].second.find('-') != std::string::npos )
            return true;
        }
      return false;
    }

    template < typename T >
    bool IsAlignment (const std::vector < T  >&data)
    {
       BOOST_STATIC_ASSERT((boost::is_base_and_derived<std::pair<std::string,std::string>,T>::value
                            || boost::is_same<std::pair<std::string,std::string>,T>::value));
      for (int i = 0; unsigned (i) < data.size (); ++i)
        if (data[i].second.length () != data[0].second.length ())
          return 0;

      return 1;
    }

    template<typename Iterator>
    bool validForPolyAnalysis( Iterator beg,
                               Iterator end )
    {
      //ensure that Iterator is an iterator type for a container
      //containing things in the inheritance hierarchy of Sequence::Seq
      //The assertion assumes that std::iterator_traits is specialized for Iterator
      BOOST_STATIC_ASSERT( (boost::is_base_and_derived< 
                            std::pair<std::string,std::string>,
                            typename std::iterator_traits<Iterator>::value_type
                            >::value || 
                            boost::is_same< 
                            std::pair<std::string,std::string>,
                            typename std::iterator_traits<Iterator>::value_type
                            >::value));
      while(beg < end)
        {
          if (std::find_if(beg->second.begin(),beg->second.end(),
                           Sequence::invalidPolyChar())
              != beg->second.end())
            {
              return false;
            }
          ++beg;
        }
      return true;
    }

    template < typename T >
    unsigned UnGappedLength (const std::vector <T>&data)
    {
      BOOST_STATIC_ASSERT((boost::is_base_and_derived<std::pair<std::string,std::string>,T>::value
                           || boost::is_same<std::pair<std::string,std::string>,T>::value));
      bool site_gapped = 0;
      int len = 0;
      if (!IsAlignment(data))
        return Sequence::SEQMAXUNSIGNED;

      for (size_t j = 0; j < data[0].second.length (); ++j)

        {
          site_gapped = 0;
          for (size_t i = 0;  i < data.size ();  ++i)
            {
              if (data[i][j] == '-')
                {
                  site_gapped = 1;
                  i = data.size();
                }
            }
          if (!(site_gapped))
            ++len;
        }
      return len;
    }

    template <typename T>
    void RemoveGaps (std::vector <T> &data)
    {

       BOOST_STATIC_ASSERT((boost::is_base_and_derived<std::pair<std::string,std::string>,T>::value
                            || boost::is_same<std::pair<std::string,std::string>,T>::value));
       size_t i, j,datasize=data.size();
       size_t length = data[0].second.length ();
      std::vector < std::string > ungapped_sequences(data.size());
      bool site_is_gapped;

      for (i = 0; i < length; ++i)
        {       //iterate over sites
          for ( j = 0, site_is_gapped = 0;
                j < datasize;  ++j)
            {
              if (data[j].second[i] == '-')
                {
                  site_is_gapped = 1;
                  j = datasize;
                }
            }
          if (!(site_is_gapped))
            {
              for ( j = 0 ; j != data.size();  ++j)
                ungapped_sequences[j] += data[j].second[i];
            }
        }

      //redo the data
      for (j = 0; j != datasize ; ++j)
        {
          data[j] = T(data[j].first,ungapped_sequences[j]);
        }
    }

    //only remove gaps from the beginning
    //and end of an alignment
    template < typename T >
    void RemoveTerminalGaps (std::vector <T>&data)
    {
      BOOST_STATIC_ASSERT((boost::is_base_and_derived<std::pair<std::string,std::string>,T>::value
                           || boost::is_same<std::pair<std::string,std::string>,T>::value));
      size_t i,j,length = data[0].second.length ();     //how much we have to iterate over
      std::vector < std::string > trimmed_sequences;
      size_t leftmost, rightmost, numUngapped;
      size_t offset;
      size_t size = data.size ();

      leftmost = SEQMAXUNSIGNED;
      rightmost = length + 1;   //offset by one b/c its an array...

      //find the leftmost site where all sites in the alignment are ungapped
      for (i = 0; i < length; ++i)
        {       //iterate over sites
          for (numUngapped = 0, j = 0; j != data.size (); ++j)
            {
              if (data[j].second[i] != '-')
                ++numUngapped;
            }
          if (numUngapped == size)
            {
              leftmost = i;
              i = length + 1;
            }
        }

      //find the rigthmost site where all sites in the alignment are ungapped
      bool exit_condition = false;
      for (i = length - 1;  i < data[0].second.length() && exit_condition == false; --i)
        {
          for (numUngapped = 0, j = 0; j != data.size (); ++j)
            {
              if (data[j].second[i] != '-')
                ++numUngapped;
            }
          if (numUngapped == size)
            {
              rightmost = i;
              exit_condition = true;
            }
        }

      //now, fill the array of trimmed sequences
      offset = rightmost - leftmost;
      for (j = 0; j != data.size (); ++j)
        trimmed_sequences.push_back (data[j].second.substr  (leftmost, rightmost-leftmost+1));

      //now, redo the seq array for the current object
      for ( j = 0; j != data.size (); ++j)
        {
          data[j] =T(data[j].first,trimmed_sequences[j]);
        }
    }

    template < typename T >
    void RemoveFixedOutgroupInsertions( std::vector<T> & data,
                                        unsigned site,
                                        const unsigned & ref )
    {
      BOOST_STATIC_ASSERT((boost::is_base_and_derived<std::pair<std::string,std::string>,T>::value
                           || boost::is_same<std::pair<std::string,std::string>,T>::value));
      size_t nsam = data.size()-1;
      size_t s = data[0].second.length()-1,max=s+1;
      for(size_t i=0 ; i<max ; --s,++i)
        {
          unsigned ngap=0;
          for(size_t ind=0;ind<data.size();++ind)
            {
              if (ind != ref && data[ind].second[s] == '-')
                {
                  ngap++;
                }
            }
          if(ngap==nsam)
            {
              for(unsigned ind=0;ind<data.size();++ind)
                {
                  data[ind].second.erase(s,1);
                }
            }
        }
    }

    template < typename T >
    std::vector < T >Trim (const std::vector < T >&data,
                           const std::vector <int> &sites) throw (Sequence::SeqException)
    {
      BOOST_STATIC_ASSERT((boost::is_base_and_derived<std::pair<std::string,std::string>,T>::value
                           || boost::is_same<std::pair<std::string,std::string>,T>::value));
      size_t i, j, numseqs = data.size (), numIntervals =  sites.size ();
      unsigned start, stop;
      std::vector < T >trimmedData(numseqs);
      std::vector < std::string > trimmedTemp(numseqs);
      if (sites.empty ())
        {
          throw SeqException ("Sequence::Alignment::Trim(): empty vector of positions passed to function");
        }
      if (numIntervals % 2 != 0)
        {
          throw SeqException ("Sequence::Alignment::Trim(): odd number of positions passed");
        }

      for (i = 0; i < numIntervals; i += 2)
        {
          start = sites[i];
          stop = sites[i + 1];
          for (j = 0; j < numseqs; ++j)
            {
              trimmedTemp[j] += data[j].second.substr (start, stop - start + 1);
            }
        }

      for (i = 0; i < numseqs; ++i)
        {
          trimmedData[i]=T(data[i].first,trimmedTemp[i]);
        }

      return trimmedData;
    }

    template < typename T >
    std::vector < T >TrimComplement (const std::vector < T >&data,
                                     const std::vector < int > &sites) 
      throw (Sequence::SeqException)
    {
      BOOST_STATIC_ASSERT((boost::is_base_and_derived<std::pair<std::string,std::string>,T>::value
                           || boost::is_same<std::pair<std::string,std::string>,T>::value));
      std::vector < int >newSites;
      size_t i, j, start, stop, numseqs = data.size (), numIntervals = sites.size (), lastval;

      if (sites.empty ())
        {
          throw SeqException ("Sequence::Alignment::TrimComplement(): empty vector of positions passed to function");
        }
      if (numIntervals % 2 != 0)
        {
          throw SeqException ("Sequence::Alignment::TrimComplement(): odd numer of positions passed to function");
        }

      std::vector < T >trimmedData(numseqs);
      std::vector < std::string > trimmedTemp(numseqs);

      size_t odd_even;
      if (sites[0] == 0)
        {
          for (i = 1; i < numIntervals; ++i)
            {
              odd_even = i+1;
              if (odd_even%2==0.)
                {
                  newSites.push_back (sites[i] +  1);
                }
              else if (odd_even%2!=0.)
                {
                  newSites.push_back (sites[i] -  1);
                }
            }
        }
      else if (sites[0] > 0)
        {
          newSites.push_back (0);
          for (i = 0; i < numIntervals; ++i)
            {
              odd_even = i+1;
              if (odd_even%2==0.)
                {
                  newSites.push_back (sites[i] +  1);
                }
              else if (odd_even%2!=0.)
                {
                  newSites.push_back (sites[i] -  1);
                }
            }
        }

      lastval = newSites[newSites.size () - 1];
      newSites.pop_back ();
      numIntervals = newSites.size ();
      for (i = 0; i < numIntervals; i += 2)
        {
          start = newSites[i];
          stop = newSites[i + 1];
          for (j = 0; j < numseqs; ++j)
            {
              trimmedTemp[j] +=
                data[j].second.substr (start,  stop - start + 1);
            }
        }

      for (j = 0; j < numseqs; ++j)
        {
          trimmedTemp[j] += data[j].second.substr (lastval);
        }

      for (i = 0; i < numseqs; ++i)
        {
          trimmedData[i] = T(data[i].first,trimmedTemp[i]);
        }

      return trimmedData;
    }

  }
}

---