alignment_processing.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_ALIGNMENT_PROCESSING_HPP__
#define SHORE_ALIGNMENT_PROCESSING_HPP__

#include <deque>
#include <iostream>
#include <list>
#include <map>
#include <set>
#include <utility>
#include <vector>

#include "shore/datatype/alignment_string_ops.hpp"
#include "shore/datatype/datatypes.hpp"
#include "shore/processing/feed.hpp"
#include "shore/processing/filter.hpp"
#include "shore/processing/extractor.hpp"
#include "shore/program/av_parser.hpp"
#include "shore/datatype/coor.hpp"
#include "shore/base/stringops.hpp"

namespace shore {

class serialization_core_access;

class alignment_queue
{
 private:

        std::deque<std::string> m_ids;
        std::deque<int> m_chromosomes;
        std::deque<long> m_positions;
        std::deque<size_t> m_sizes;
        std::deque<shore::Strand> m_strands;
        std::deque<int> m_hits;
        std::deque<std::string> m_alignments;
        std::deque<int> m_mismatches;
        std::deque<int> m_indexes;
        std::deque<shore::alignment::InsertState> m_pe_states;
        std::deque<std::string> m_quals1;
        std::deque<std::string> m_quals2;
        std::deque<std::string> m_quals3;


        std::deque<int> m_mapquals;
        std::deque<std::string> m_readgroups;
        std::deque<int> m_num_reads;
        std::deque<shore::read::PeOrientation> m_oris;
        std::deque<int> m_rnext_chromosomes;
        std::deque<long> m_rnext_positions;
        std::deque<shore::Strand> m_rnext_strands;

        size_t m_size;

        shore::alignment m_front;
        shore::alignment m_back;

        bool m_fifo;

 public:

        typedef shore::alignment current_type;
        typedef shore::alignment append_type;

        alignment_queue(const bool fifo=true);

        void push_back(const shore::alignment& f);

        void pop_front();
        void pop_back();

        void get(const size_t i,shore::alignment& res) const;

        size_t size() const;

        bool empty() const;

        const shore::alignment& front() const;
        const shore::alignment& back() const;
        shore::alignment& front();
        shore::alignment& back();

        void clear();

        const shore::alignment& current() const;
        bool has_data() const;
        void next();
        void append(const shore::alignment& f);
        void flush() {}
};

class alignment_3px
{
 private:

        typedef bool(*const FtOrd_t) (const shore::alignment&,const shore::alignment&);
        
        static bool cmp_coor(const shore::alignment& f0,const shore::alignment& f1);

        long m_xsize;

        alignment_queue m_fbuf;
        std::multiset<shore::alignment,FtOrd_t> m_rbuf;

        shore::refseq_coor m_flushpos;

        alignment_helper m_ah;

 public:

        typedef shore::alignment current_type;
        typedef shore::alignment append_type;

        alignment_3px(const long xsize);

        const shore::alignment& current() const;
        bool has_data() const;
        void next();
        void append(const shore::alignment& f);
        void flush();
};

class poissonifier
{
 public:
        
        typedef shore::alignment current_type;
        typedef shore::alignment append_type;

        class pf_core
        {
                public:

                        struct pos
                        {
                                shore::refseq_coor c;
                                double n;
                                double m1;
                                //          for right window.
                                double m2;
                                double v1;
                                double v2;
                                double d1;
                                double d2;
                                double d;
                                //double k1;
                                //double k2;
                                //double k;
                                double l;
                        };

                private:

                        std::deque<pos> m_limits;

                        size_t m_width;
                        size_t m_ws;

                        shore::refseq_coor m_centre;
                        std::deque<double> m_window;

                        bool m_verbose;


                        void move_window(const shore::refseq_coor& c);

                public:

                        pf_core(const size_t w=24);
                        virtual ~pf_core() {}

                        void set_verbose(const bool v)
                        {
                                m_verbose=v;
                        }

                        virtual void add(const shore::refseq_coor& c,const double w);
                        virtual void flush();

                        pos next_limit();
                        bool has_limit() const;
        };

 private:

        alignment_queue m_inbuff;
        alignment_queue m_outbuff;
        alignment_queue m_inbufr;
        alignment_queue m_outbufr;

        alignment_helper m_ah;

        pf_core m_pf;
        pf_core m_pr;
        alignment_3px m_3px;

        std::map<shore::refseq_coor,size_t> m_revlimits;

        bool m_isflushed;

        double m_lmax;

        shore::feed<shore::alignment> m_duplicates;

 public:

        poissonifier(const size_t w=24);

        shore::feed<shore::alignment> &duplicates();
        
        void set_verbose(const bool v,const bool w)
        {
                m_pf.set_verbose(v);
                m_pr.set_verbose(w);
        }

        const shore::alignment& current() const;
        bool has_data() const;
        void next();
        void append(const shore::alignment& f);
        void flush();

        double get_lmax() const { return m_lmax; }
};

inline std::ostream& operator<<(std::ostream& os,const poissonifier::pf_core::pos& p)
{
 os<<"chr="<<p.c.chromosome<<"; pos="<<p.c.position
        <<"; n="<<p.n
        <<"; m1="<<p.m1
        <<"; m2="<<p.m2
        <<"; v1="<<p.v1
        <<"; v2="<<p.v2
        <<"; d1="<<p.d1
        <<"; d2="<<p.d2
        <<"; d="<<p.d
        <<"; l="<<p.l;
 return os;
}

class seqid_translator
{
 private:

        std::vector<int> m_seqmap;
        shore::alignment m_current;
        bool m_hasdata;

 public:

        typedef shore::alignment append_type;
        typedef shore::alignment current_type;

        seqid_translator(const std::vector<int> &seqmap);

        const shore::alignment &current() const;

        bool has_data() const;

        void next();

        void append(const shore::alignment &f);

        void flush() {}
};

class besthit_filter
{
 private:

        typedef std::pair<shore::Strand,std::pair<shore::refseq_coor,std::string> > key_type;

        alignment_queue m_queue;
        shore::alignment m_fidwait;
        bool m_hasdata;
        bool m_idwait;
        std::string m_currentid;

        int m_current_ed1;
        int m_current_ed2;
        int m_current_hits1;
        int m_current_hits2;
        std::set<key_type> m_current_alignments1;
        std::set<key_type> m_current_alignments2;

        shore::signal<const shore::alignment&> m_sigdiscardhit;
        shore::signal<void> m_sigdiscardhit_flush;


        void clearbesthit1();
        void clearbesthit2();
        bool tweak_front();
        void prune_queue();

 public:

        typedef shore::alignment append_type;
        typedef shore::alignment current_type;

        besthit_filter();

        const shore::alignment &current() const;

        bool has_data() const;

        void next();

        void append(const shore::alignment &f);

        void flush();

        shore::signal<const shore::alignment&>& sigdiscardhit();
        shore::signal<void>& sigdiscardhit_flush();
};

class alignment_statistic
{
 private:

        struct element
        {
                int sequence;
                size_t size;
                int index;
                int edit_distance;
                int hits;
                char strand;


                element() {}
                element(const shore::alignment &f);
                element(const shore::read &f);


                template<class Archive>
                void serialize(Archive& ar,const unsigned int version)
                {
                        ar&sequence;
                        ar&size;
                        ar&index;
                        ar&edit_distance;
                        ar&hits;
                        ar&strand;
                }
        };

        typedef bool(*cmp_t)(const element &,const element &);

        std::map<element,size_t,cmp_t> m_stats;

        static bool cmp(const element &e1,const element &e2);

        friend class shore::serialization_core_access;

 public:

        typedef shore::alignment append_type;

        alignment_statistic();

        void append(const shore::alignment &f);
        void append(const shore::read &f);
        void append(const alignment_statistic &other);
        void flush() {}

        void save(std::ostream &out) const;
        void save(const std::string &outfile) const;

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

        void update(const std::string &statsfile,const std::string &mapfile,
                                const std::string &reffile=std::string());

        void clear();

        double calc_total() const;
};

class hits_correction_single_pipe
:public shore::pipe_facade<shore::hits_correction_single_pipe,
                           shore::alignment,shore::alignment>
{
 private:

        friend class pipeline_core_access;

        shore::alignment m_current;
        alignment_queue m_queue;
        long m_current_hits;


        void next();
        void prepare(const shore::alignment& f);
        void append(const shore::alignment& f);

        void flush();
};

class flatten_alignment_pipe
:public shore::pipe_facade<shore::flatten_alignment_pipe,
                           shore::alignment,shore::read>
{
 private:

        friend class pipeline_core_access;

        shore::alignment_helper m_ah;

        std::vector<std::string> m_current_ids;
        shore::read m_current;

        bool m_remove_redundant;


        void next();
        void prepare(const shore::alignment &f);
        void append(const shore::alignment &f);
        void flush();

 public:

        flatten_alignment_pipe(const bool remove_redundant=true);
};

class mapping_pair_finder_pipe
:public shore::pipe_facade<mapping_pair_finder_pipe,
                           shore::alignment,std::vector<shore::alignment> >
{
 private:

        friend class pipeline_core_access;

        typedef std::pair<shore::Strand,std::string> pairkey_t;

        typedef std::pair<std::vector<shore::alignment>,bool> paircollect_t;

        pairkey_t m_key;
        std::deque<paircollect_t> m_pairwait_q;
        uint64_t m_queueofs;
        std::map<shore::refseq_coor,std::map<pairkey_t,size_t> > m_nextones;

        bool m_find_all_concordant;
        long m_find_unique_distance;


        void emit_next();

        void next();

        void prepare(const shore::alignment &a);
        void append(const shore::alignment &a);
        void flush();

 public:

        mapping_pair_finder_pipe(bool conc,long any);
};

class mapping_pair_joiner_pipe
:public shore::pipe_facade<mapping_pair_joiner_pipe,
                           std::vector<shore::alignment>,shore::alignment>
{
 private:

        friend class pipeline_core_access;

        shore::alignment_helper m_ah;
        shore::alignment m_current;
        std::string m_buf;

        void next();

        void prepare(const std::vector<shore::alignment> &a);
        void append(const std::vector<shore::alignment> &a);
        void flush();
};

class alignment_trimmer_pipe
:public shore::pipe_facade<alignment_trimmer_pipe,shore::alignment,shore::alignment>
{
 private:

        friend class pipeline_core_access;

        shore::alignment_helper m_ah;
        shore::alignment_tokenizer m_at;
        shore::alignment m_current;
        const shore::refseq_range m_range;

        void append(const shore::alignment &a);

 public:

        alignment_trimmer_pipe(const shore::refseq_range &r);
};

class alignment_filter_pipe
:public shore::pipe_facade<alignment_filter_pipe,
                           shore::alignment,shore::alignment>
{
 public:

        class outofrange_exception
        :public std::exception
        {};

        struct config
        {
                int p3x;

                bool join_happy;
                long join_unique;

                std::pair<size_t,size_t> hitrange;
                std::pair<size_t,size_t> edrange;
                shore::OptionalStrand strand;
                std::set<size_t> lengths;
                std::set<std::pair<shore::alignment::InsertState,int> > peflags;

                double duplicates;

                shore::refseq_range posrange;

                int wpoiss;

                bool position_sorted;

                config();

                void publish_options(shore::av_parser& p,
                  const char*const D="strand,T",
                  const char*const F="duplicate-estimation,F",
                  const char*const N="read-lengths,N",
                  const char*const X="p3fix,X",
                  const char*const B="duplicates,B",
                  const char*const R="region,R",
                  const char*const M="mm-range,M",
                  const char*const H="hits-range,H",
                  const char*const pe="peflags",
                  const char*const jh="join-happy",
                  const char*const ju="join-unique");

                void sanity_check(shore::av_parser& p);

                bool operator==(const config& c2) const
                {
                        return
                                (p3x==c2.p3x)
                                &&(join_happy==c2.join_happy)
                                &&(join_unique==c2.join_unique)
                                &&(hitrange==c2.hitrange)
                                &&(edrange==c2.edrange)
                                &&(strand==c2.strand)
                                &&(lengths==c2.lengths)
                                &&(peflags==c2.peflags)
                                &&(duplicates==c2.duplicates)
                                &&(posrange==c2.posrange)
                                &&(wpoiss==c2.wpoiss);
                }

                bool is_default() const
                {
                        return (*this)==(config());
                }
        };

 private:

        friend class pipeline_core_access;

        const config m_defaults;
        const config m_conf;

        shore::alignment_helper m_ah;

        shore::feed<shore::alignment> m_feed;
        shore::extractor<shore::alignment> m_extractor;

        shore::splitter<shore::alignment> m_rangefilter;
        mapping_pair_finder_pipe m_pairfinder;
        mapping_pair_joiner_pipe m_pairjoiner;
        shore::splitter<shore::alignment> m_firstfilter;
        shore::pipe<shore::alignment_3px> m_3primex;
        shore::pipe<shore::poissonifier> m_dupfilter;
        shore::splitter<shore::alignment> m_lastfilter;

        shore::thru<shore::alignment> m_failures;

        bool flt_hits(const shore::alignment& f);
        bool flt_editdist(const shore::alignment& f);
        bool flt_strand(const shore::alignment& f);
        bool flt_length(const shore::alignment& f);
        bool flt_PE(const shore::alignment& f);
        bool flt_position_sorted(const shore::alignment& f);
        bool flt_position_unsorted(const shore::alignment& f);

        struct flt_duplicates
        {
                config conf;

                shore::alignment_helper m_ah;

                double m_duplicacy_f;
                shore::refseq_coor m_lastpos_f;
                 std::map<shore::refseq_coor,double> m_duplicacy_r;
                
                flt_duplicates(const config& c);

                bool operator()(const shore::alignment& f);
        };

        void init_filters();

        template<typename T>
        void connect_rangefilter(T &source)
        {
                if(!m_rangefilter.empty())
                {
                        source|m_rangefilter;
                        connect_pejoin(m_rangefilter);
                }
                else connect_pejoin(source);
        }

        template<typename T>
        void connect_pejoin(T &source)
        {
                if(m_conf.join_happy||(m_conf.join_unique>0))
                {
                        source|m_pairfinder|m_pairjoiner;
                        connect_filter1(m_pairjoiner);
                }
                else connect_filter1(source);
        }

        template<typename T>
        void connect_filter1(T &source)
        {
                if(!m_firstfilter.empty())
                {
                        source|m_firstfilter;
                        m_firstfilter.failures()|m_failures;
                        connect_3primex(m_firstfilter);
                }
                else connect_3primex(source);
        }

        template<typename T>
        void connect_3primex(T &source)
        {
                if(m_conf.p3x!=0)
                {
                        source|m_3primex;
                        connect_dup(m_3primex);
                }
                else connect_dup(source);
        }

        template<typename T>
        void connect_dup(T &source)
        {
                if(m_conf.wpoiss!=0)
                {
                        source|m_dupfilter;
                        m_dupfilter->duplicates()|m_failures;
                        connect_last(m_dupfilter);
                }
                else connect_last(source);
        }

        template<typename T>
        void connect_last(T &source)
        {
                if(!m_lastfilter.empty())
                {
                        source|m_lastfilter|m_extractor;
                        m_lastfilter.failures()|m_failures;
                }
                else source|m_extractor;
        }

        void next();

        void prepare(const shore::alignment &a);
        void append(const shore::alignment &a);
        void flush();

 public:

        alignment_filter_pipe(const config &c);

        shore::thru<shore::alignment> &failures();
};

class alignment_condenser_pipe
:public shore::pipe_facade<alignment_condenser_pipe,
                           shore::alignment,shore::alignment>
{
 public:

        enum IDStyle
        {
                ID_PRESERVE,
                ID_PREFIXED,
                ID_PLAIN
        };

 private:

        friend class pipeline_core_access;

        int m_compact_aln;
        int m_qgrid;
        bool m_discard_qual_other;

        IDStyle m_readid_style;

        shore::alignment m_buf;

        std::vector<std::pair<bool,int> > m_tokbuf;

        shore::alignment_tokenizer m_tok;

        std::map<std::string,int> *m_readid_map;


        template<typename Iter>
        void condense_quality(Iter beg,Iter end) const
        {
                for(size_t i=0;(beg!=end)&&(i<m_tokbuf.size());++i)
                {
                        // Match
                        if(m_tokbuf[i].first)
                        {
                                Iter save=beg;

                                int qmin=uint8_t(*beg);
                                int qmax=uint8_t(*beg);

                                for(int j=0;(beg!=end)&&(j<m_tokbuf[i].second);
                                        ++j,++beg)
                                {
                                        const int qcur=uint8_t(*beg);

                                        if((qcur>(qmin+m_qgrid))
                                           ||(qcur<(qmax-m_qgrid)))
                                        {
                                                while(save!=beg)
                                                {
                                                        (*save)=char(qmin);
                                                        ++save;
                                                }

                                                qmin=qmax=qcur;
                                        }
                                        else
                                        {
                                                qmin=std::min(qmin,qcur);
                                                qmax=std::max(qmax,qcur);
                                        }
                                }
                                while(save!=beg)
                                {
                                        (*save)=char(qmin);
                                        ++save;
                                }
                        }
                        // Not a match
                        else
                        {
                                for(int j=0;(beg!=end)&&(j<m_tokbuf[i].second);
                                        ++j,++beg)
                                {}
                        }
                }
        }

        void append(const shore::alignment &a);

 public:

        alignment_condenser_pipe(const int compact_aln=1,const int qgrid=0,
                            const bool discard_qual_other=false,
                            const IDStyle readid_style=ID_PRESERVE,
                            std::map<std::string,int> *const readid_map=0);
};

S_MAKE_ENUM_TRAITS(alignment_condenser_pipe::IDStyle)

} // namespace

#endif // SHORE_ALIGNMENT_PROCESSING_HPP__