2dex.hpp

Go to the documentation of this file.

/*
 * Copyright 2008,2009,2010,2011,2012 Stephan Ossowski, Korbinian Schneeberger,
 * Felix Ott, Joerg Hagmann, Alf Scotland, Sebastian Bender
 * 
 * This file is part of SHORE.
 * 
 * SHORE 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.
 * 
 * SHORE 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
 * along with SHORE.  If not, see <http://www.gnu.org/licenses/>.
 */


#ifndef SHORE_ALGO_TWODEX_HPP__
#define SHORE_ALGO_TWODEX_HPP__

#include <iterator>
#include <set>
#include <string>
#include <vector>

#include <boost/iterator/iterator_facade.hpp>

#include "shore/algo/d2tree.hpp"
#include "shore/datatype/datatypes.hpp"
#include "shore/container/intpack.hpp"
#include "shore/datatype/alignment_string_ops.hpp"
#include "shore/base/memops.hpp"
#include "shore/datatype/coor.hpp"
#include "shore/base/stringops.hpp"

namespace shore {

class size_iterator
:public boost::iterator_facade<
  size_iterator,size_t const,boost::random_access_traversal_tag>
{
 private:

        friend class boost::iterator_core_access;

        size_t m_value;

        void increment()
        {
                ++m_value;
        }

        void decrement()
        {
                --m_value;
        }

        void advance(difference_type n)
        {
                m_value+=n;
        }

        difference_type distance_to(const size_iterator& other) const
        {
                return other.m_value-m_value;
        }

        bool equal(const size_iterator& other) const
        {
                return other.m_value==m_value;
        }

        const reference dereference() const
        {
                return m_value;
        }

 public:
        
        size_iterator(const size_t value=0)
        :m_value(value)
        {}
        
        size_iterator(const size_iterator& other)
        :m_value(other.m_value)
        {}

        size_iterator& operator=(const size_iterator& other)
        {
                m_value=other.m_value;
                return *this;
        }

        const size_t operator[](difference_type n) const
        {
                return m_value+n;
        }
};

struct blockcmp_pos
{
        const intpack::const_iterator& pos;

        blockcmp_pos(intpack::const_iterator& p)
        :pos(p)
        {}

        blockcmp_pos(const intpack& p)
        :pos(p.begin())
        {}

        bool operator()(const size_t i,const size_t j) const
        {
                return pos[i]<pos[j];
        }

        bool operator()(const size_t i,const std::pair<size_t,size_t>& posend) const
        {
                return pos[i]<posend.first;
        }

        bool operator()(const std::pair<long,long>& posend,const size_t i) const
        {
                return posend.first<long(pos[i]);
        }
};

struct blockcmp_end
{
        const intpack::const_iterator& pos;
        const intpack::const_iterator& sizes;

        blockcmp_end(intpack::const_iterator& p,intpack::const_iterator& s)
        :pos(p),sizes(s)
        {}

        blockcmp_end(const intpack& p,const intpack& s)
        :pos(p.begin()),sizes(s.begin())
        {}

        bool operator()(const size_t i,const size_t j) const
        {
                return (pos[i]+sizes[i])<(pos[j]+sizes[j]);
        }

        bool operator()(const size_t i,const std::pair<size_t,size_t>& posend) const
        {
                return (pos[i]+sizes[i])<posend.second;
        }

        bool operator()(const std::pair<long,long>& posend,const size_t i) const
        {
                return posend.second<long(pos[i]+sizes[i]);
        }
};

struct blockcmp_posend
{
        const intpack::const_iterator& pos;
        const intpack::const_iterator& sizes;

        blockcmp_posend(intpack::const_iterator& p,intpack::const_iterator& s)
        :pos(p),sizes(s)
        {}

        blockcmp_posend(const intpack& p,const intpack& s)
        :pos(p.begin()),sizes(s.begin())
        {}

        bool operator()(const size_t i,const size_t j) const
        {
                return pos[i]<(pos[j]+sizes[j]);
        }

        bool operator()(const size_t i,const std::pair<size_t,size_t>& posend) const
        {
                return pos[i]<posend.second;
        }

        bool operator()(const std::pair<long,long>& posend,const size_t i) const
        {
                return posend.first<long(pos[i]+sizes[i]);
        }
};

class twodex
{
 public:

        enum D2Type
        {
                D2_SIZE,
                D2_INCLUSIVE,
                D2_EXCLUSIVE,
                D2_MAPLIST,
                D2_CIGAR
        };

        enum QueryMode
        {
                QRY_OVERLAPPING,
                QRY_INCLUDING,
                QRY_INCLUDED
        };

        struct build_spec
        {
                std::string colspec;
                size_t blocksize;
                size_t maxgap;
                D2Type d2type;
                char commentchar;

                build_spec();
        };

        struct seq_info
        {
                size_t id;
                long offset;
                long firstpos;

                seq_info(const size_t i)
                :id(i),offset(std::numeric_limits<long>::min()),
                 firstpos(std::numeric_limits<long>::max())
                {}
        };

        struct header
        {
                size_t blocksize;

                int chr_col;
                int pos_col;
                int d2_col;
                D2Type d2type;
                char commentchar;

                bool sorted;

                std::streamsize seq_total;
                std::streamsize file_total;
                size_t nnzblocks;

                std::map<std::string,seq_info> sequences;
                std::vector<std::map<std::string,seq_info>::iterator> seqv;

                int bitsperid;
                int bitsperpos;
                int bitspersize;

                size_t id_bytes;
                size_t pos_bytes;
                size_t size_bytes;

                header();

                size_t parse(std::istream& in);
                void cleardata();

                void print(std::ostream& os) const;
        };

        struct bytecount
        {
                std::streamsize count;

                bytecount();
                void operator()(const std::string& s);
        };

 private:

        struct block
        {
                typedef std::pair<shore::refseq_coor,shore::refseq_coor> range_t;

                std::vector<range_t> ranges;
                std::vector<range_t> rangebuf;

                void add_range(const shore::refseq_coor& p,const long e,const long maxgap)
                {
                        range_t orange(p,shore::refseq_coor(p.chromosome,e));

                        for(size_t i=0;i<ranges.size();++i)
                        {
                                if(maxgap<0)
                                {
                                        // maxgap not set, blocks may cross chromosome borders
                                        if(ranges[i].first<orange.first)
                                                orange.first=ranges[i].first;
                                        if(orange.second<ranges[i].second)
                                                orange.second=ranges[i].second;
                                }
                                else if((ranges[i].first.chromosome!=p.chromosome)
                                        ||(ranges[i].second.position<=(p.position-maxgap))
                                        ||((e+maxgap)<=ranges[i].first.position))
                                {
                                        rangebuf.push_back(ranges[i]);
                                }
                                else
                                {
                                        if(ranges[i].first.position<orange.first.position)
                                                orange.first.position=ranges[i].first.position;

                                        if(orange.second.position<ranges[i].second.position)
                                                orange.second.position=ranges[i].second.position;
                                }
                        }

                        ranges.clear();
                        ranges.swap(rangebuf);
                        ranges.push_back(orange);
                }
        };

        header m_header;

        shore::mmapping m_mmapping;

        const char *m_data;
        intpack::const_iterator m_idbeg;
        intpack::const_iterator m_idend;
        intpack::const_iterator m_posbeg;
        intpack::const_iterator m_posend;
        intpack::const_iterator m_sizebeg;
        intpack::const_iterator m_sizeend;

        d2tree<blockcmp_pos,blockcmp_end,blockcmp_posend> m_tree;        

        QueryMode m_querymode;


        void load(const std::string &indexfile);

 public:

        static void build(const std::string &fn,const std::string &colspec,
                          const size_t blocksize,const size_t maxgap,
                          const D2Type d2type,const char commentchar,
                          std::ostream *const log=0);

        static void build(const std::string &fn,const build_spec &bs,
                          std::ostream *const log=0);

        static void parse_xend(long *const xend,const long pos,
                               const D2Type d2type,const std::string& d2,
                               shore::alignment_helper &ah);

        twodex(const std::string& indexfile,const QueryMode qm=QRY_OVERLAPPING);

        void set_querymode(const QueryMode qm);

        const header& get_header() const;

        void query(std::set<size_t>& res,const std::pair<long,long>& r) const;
        void query(std::set<size_t>& res,const shore::genomic_range& r) const;

        std::pair<long,long> range_enc(const shore::genomic_range& r) const;

        int test(const shore::genomic_range& r1,const shore::genomic_range& r2) const;
        int test(const shore::genomic_range& r1,const std::pair<long,long>& r2) const;

        static int test(const std::pair<long,long>& r1,
                        const std::pair<long,long>& r2,
                        const QueryMode querymode);
};

S_MAKE_ENUM_TRAITS(shore::twodex::QueryMode)
S_MAKE_ENUM_TRAITS(shore::twodex::D2Type)

} // namespace shore

#endif // SHORE_ALGO_TWODEX_HPP__