sam.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_IO_SAM_HPP__
#define SHORE_IO_SAM_HPP__

#include <map>
#include <string>

#include "shore/processing/pipeline.hpp"
#include "shore/datatype/alignment_string_ops.hpp"
#include "shore/datatype/alignment_processing.hpp"
#include "shore/datatype/datatypes.hpp"
#include "shore/stream/streams.hpp"
#include "shore/fmtio/text.hpp"

namespace shore {

enum SAMFlags
{
        SAM_PE_PAIRED_END=0x001,
        SAM_PE_PROPER_PAIR=0x002,
        SAM_UNMAPPED_SELF=0x004,
        SAM_PE_UNMAPPED_MATE=0x008,
        SAM_REVERSE_SELF=0x010,
        SAM_REVERSE_MATE=0x020,
        SAM_PE_READ1=0x040,
        SAM_PE_READ2=0x080,
        SAM_NO_PRIMARY_ALIGNMENT=0x100,
        SAM_QUALITY_FILTERED=0x200,
        SAM_DUPLICATE_READ=0x400
};

class sam_base
{
 public:

        class leftover_pipe
        :public shore::pipe_facade<leftover_pipe,shore::read,shore::read>
        {
                private:

                        friend class sam_base;

                        leftover_pipe();

                public:

                        void next();
                        void prepare(const shore::read &r);
                        void append(const shore::read &r);
                        void flush();
        };

        struct config
        {
                std::string refseq_fn;

                int max_ed;
                int max_gaps;
                int max_gapsize;
                int max_hits;

                int forced_readindex;
                shore::alignment::InsertState forced_pestate;

                bool discard_softclipped;

                bool ignore_xareads;
                
                bool bismark;

                config(const std::string &reffn=std::string());
        };

        struct cigarstats
        {
                int edit_distance;
                int ngaps;
                int nsoftclipped;
        };

 private:

        leftover_pipe m_leftovers;

 protected:

        config m_conf;

        static const int MAXTAGS=3224;

        std::pair<int,shore::alignment::InsertState> m_peconv[0x801];
        std::pair<int,shore::read::InsertState> m_peconv_left[0x801];

        int m_tid[256][256];

        struct samheader
        {
                struct sh_hd
                {
                        std::string vn;
                        std::string so;
                };

                struct sh_sq
                {
                        std::string sn;
                        size_t ln;
                        std::string as;
                        std::string m5;
                        std::string sp;
                        std::string ur;
                };

                sh_hd hd;

                std::vector<sh_sq> sq;
                std::map<std::string,size_t> m5_map;
                std::map<std::string,size_t> sn_map;


                void append(const std::string& s);

                void flush() {}
        };

        shore::sink<samheader,std::string> m_header;

        std::map<std::string,size_t> m_refmap;
        typedef std::map<std::string,size_t>::const_iterator ichr_t;
        std::vector<shore::sequence_record> m_refseq;
        std::vector<int> m_bamid2seqid;

        void load_reference();

        void xa_parse(shore::alignment &f,std::string &seq,std::string &cigar,
                      bool &isrev,cigarstats &cigsta,size_t &xaofs,
                      const std::string &xa);

 public:
 
        sam_base(const config &c);

        leftover_pipe &leftovers();

        void set_max_ed(const int ed);
        void set_max_gaps(const int g);
        void set_max_gapsize(const int g);
        void set_max_hits(const int h);

        void discard_softclipped(const bool d);

        static void cigar2aln(std::string &aln,cigarstats &cs,
                              const std::string &read,const std::string &cigar,
                              const shore::sequence_record &chr,
                              const size_t start,const int max_gapsize=1);
};

class sam_reader
:public sam_base
{
 public:

        typedef shore::alignment current_type;

 private:

        static const int MAXCOLS=MAXTAGS+11;

        enum SAMCol
        {                       // SAM doc column names
                SC_ID=0,            // QNAME
                SC_FLAG=1,          // FLAG
                SC_CHRNAME=2,       // RNAME
                SC_POS=3,           // POS
                SC_MAPPINGQUALITY=4,// MAPQ
                SC_CIGAR=5,         // CIGAR
                SC_MATE_SEQNAME=6,  // MRNM
                SC_MATE_POS=7,      // MPOS
                SC_INSERTSIZE=8,    // ISIZE
                SC_READSEQ=9,       // SEQ
                SC_QUALITY=10       // QUAL
        };

        typedef boost::function<void (const std::string*const)> tagfun_t;
        tagfun_t m_tagfuns[MAXTAGS+1];

        shore::line_reader m_in;

        shore::alignment m_buf1;
        shore::alignment m_buf2;
        shore::alignment *m_current;
        shore::alignment *m_last;
        shore::read m_current_left;
        
        std::string m_linebuf;

        std::string m_splitbuf[MAXCOLS];
        int m_ntok;
        cigarstats m_cigsta;

        bool m_hasdata;
        bool m_needinit;

        const std::string* m_xastring;
        size_t m_xapos;
        shore::refseq_coor m_xa_parentcoor;
        
        std::string m_xrtag;
        std::string m_xgtag;
        std::string m_xmtag;
        void convert_alignment_string();

        bool m_emit_as_leftover;
        int m_currentsamflag;
        bool m_current_is_revcomp;

        std::string m_rangeindexfile;

        shore::alignment::cmp_type m_cmp;

        void tagfun_invalid(const std::string*const ptag);
        void tagfun_MD(const std::string*const ptag);
        void tagfun_X0(const std::string*const ptag);
        void tagfun_XA(const std::string*const ptag);
        void tagfun_NH(const std::string*const ptag);
        void tagfun_RG(const std::string*const ptag);
        void tagfun_XR(const std::string*const ptag);
        void tagfun_XG(const std::string*const ptag);
        void tagfun_XM(const std::string*const ptag);

 public:

        sam_reader(const std::string& fn,const config &c,
                   shore::alignment::cmp_type cmp=shore::alignment::cmp_null);

        ~sam_reader();

        bool has_data();

        const shore::alignment& current() const;

        void next();

        const std::string &current_line() const;

        void set_range(const shore::refseq_coor& s,const shore::refseq_coor& e,
                       const long maxrl=0l);
};

typedef shore::source<sam_reader> sam_source;

class bam_reader
:public sam_base
{
 public:

        typedef shore::alignment current_type;

 private:

        typedef boost::function<void (const char**)> tagfun_t;
        tagfun_t m_tagfuns[MAXTAGS+1];

        static const char CIGARCONV[17];
        static const char BASECONV[17];

        shore::istreams m_in;

        shore::alignment m_buf1;
        shore::alignment m_buf2;
        shore::alignment *m_current;
        shore::alignment *m_last;
        shore::read m_current_left;
        
        cigarstats m_cigsta;

        bool m_hasdata;
        bool m_needinit;

        bool m_emit_as_leftover;
        int m_currentsamflag;
        bool m_current_is_revcomp;

        shore::alignment::cmp_type m_cmp;

        size_t m_readbufsize;
        char *m_readbuf;

        size_t m_nextblocklen;

        std::string m_cigarbuf;
        std::string m_seqbuf;
        std::string m_xabuf;
        size_t m_xaofs;
        shore::refseq_coor m_xa_parentcoor;

        void ensure_buf(size_t size);
 
        void read_header();

        void tagfun_invalid(const char **ptag);
        void tagfun_default(const char **ptag);
        void tagfun_MD(const char **ptag);
        void tagfun_X0(const char **ptag);
        void tagfun_XA(const char **ptag);
        void tagfun_NH(const char **ptag);
        void tagfun_RG(const char **ptag);

 public:

        bam_reader(const std::string& fn,const config &c,
                   shore::alignment::cmp_type cmp=shore::alignment::cmp_null);
        ~bam_reader();

        bool has_data();

        const shore::alignment& current() const;

        void next();

        void set_range(const shore::refseq_coor& s,const shore::refseq_coor& e,
                       const long maxrl=0l);
};

typedef shore::source<bam_reader> bam_source;

class sam_writer
{
 private:

        shore::ostreams m_os;

        std::string m_mdbuf;
        std::string m_cgbuf;
        std::string m_sqbuf;
        std::string m_qualbuf;

        int m_peconv[9];

        int m_seqmax;

        shore::alignment_tokenizer m_alntok;

        std::set<std::string> m_allowed_readgroups;

        bool m_correct_leftover_flags;
        bool m_orphan_partners_only;

        std::map<std::string,std::pair<shore::refseq_coor,shore::Strand> > m_orphan_coors;

 public:

        enum SamSorting
        {
                SAM_SO_UNKNOWN,
                SAM_SO_UNSORTED,
                SAM_SO_COORDINATE,
                SAM_SO_QUERYNAME
        };

        typedef shore::alignment append_type;

        sam_writer(const std::string& fn,const std::vector<std::string> &allowed_readgroups);
        sam_writer(std::ostream*const os,const std::vector<std::string> &allowed_readgroups);

        void write_header(const std::string& fasta_fn=std::string(),
                          const SamSorting sort=SAM_SO_UNKNOWN);

        void set_correct_leftover_flags(const bool b);

        void set_orphan_leftovers_only(const bool b);

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

        void flush();
};

S_MAKE_ENUM_TRAITS(sam_writer::SamSorting)

typedef shore::sink<sam_writer> sam_sink;

typedef shore::sink<sam_writer,shore::read> sam_read_sink;

} // namespace

#endif // SHORE_IO_SAM_HPP__