find_shared_impl.hpp

Go to the documentation of this file.

#ifndef FIND_SHARED_IMPL_HPP_
#define FIND_SHARED_IMPL_HPP_
#include <algorithm>
#include <functional>
#include <ostream>
#include <string>
#include <vector>
#include "boost/algorithm/string/classification.hpp" //operator!
#include "boost/algorithm/string/predicate.hpp" //all
#include "boost/assert.hpp"
#include "boost/container/flat_map.hpp"
#include "boost/foreach.hpp"
#include "boost/multi_index_container.hpp"
#include "boost/multi_index/hashed_index.hpp"
#include "boost/multi_index/member.hpp"
#include "boost/multi_index/ordered_index.hpp"
#include "boost/range.hpp"
#include "vdj_pipe/detail/identity_predicate.hpp"
#include "vdj_pipe/exception.hpp"
#include "vdj_pipe/file_stream.hpp"
#include "vdj_pipe/gdst_search.hpp"
#include "vdj_pipe/min_match_length.hpp"
#include "vdj_pipe/object_ids.hpp"

namespace vdj_pipe{

struct Read_ginfo {
   std::string id_str_;
   Read_id rid_;
   Mapped_id mid_;
};

class Group_count {

public:
   Group_count() {}

   explicit Group_count(const std::size_t n_groups)
   : v_(n_groups, 0)
   {}

   unsigned group_count() const {
      return std::count_if(v_.begin(), v_.end(), Identity<unsigned>());
   }

   void operator()(const Mapped_id mid) {++v_[mid()];}
   unsigned operator[](const Mapped_id mid) const {return v_[mid()];}

   unsigned copy_count() const {
      return accumulate(v_.begin(), v_.end(), 0U);
   }

   void operator()(Group_count const& gc) {
      transform(
               v_.begin(),
               v_.end(),
               gc.v_.begin(),
               v_.begin(),
               std::plus<unsigned>()
      );

   }

   std::vector<unsigned> v_;
};

class Read_ginfo_store {
   typedef boost::multi_index_container<
            Read_ginfo,
            boost::multi_index::indexed_by<
               boost::multi_index::hashed_unique<
                  boost::multi_index::tag<struct id_tag>,
                  boost::multi_index::member<
                  Read_ginfo, Read_id, &Read_ginfo::rid_
                  >
               >,
               boost::multi_index::hashed_unique<
                  boost::multi_index::tag<struct id_str_tag>,
                  boost::multi_index::member<
                  Read_ginfo, std::string, &Read_ginfo::id_str_
                  >
               >,
               boost::multi_index::ordered_non_unique<
                  boost::multi_index::tag<struct mapped_tag>,
                  boost::multi_index::member<
                  Read_ginfo, Mapped_id, &Read_ginfo::mid_
                  >
               >
            >
   > store_t;

   typedef store_t::index<id_tag>::type id_index;
   typedef id_index::const_iterator id_iterator;

   typedef store_t::index<id_str_tag>::type id_str_index;
   typedef id_str_index::const_iterator id_str_iterator;

   typedef store_t::index<mapped_tag>::type mid_index;
   typedef mid_index::const_iterator mid_iterator;
   typedef boost::iterator_range<mid_iterator> mid_range;

public:

   struct Err : public base_exception {};

   std::size_t size() const {return m_.size();}

   Read_id insert(std::string const& id, const Mapped_id mid) {
      const Read_id rid(m_.size() + 1);
      const Read_ginfo rg = {id, rid, mid};
      m_.insert(rg);
      return rid;
   }

   mid_range by_mid(const Mapped_id mid) const {
      return boost::make_iterator_range(m_.get<mapped_tag>().equal_range(mid));
   }

   Read_ginfo const& operator[](const Read_id rid) const {
      return *m_.get<id_tag>().find(rid);
   }

private:
   store_t m_;
};

struct Sequence_gstats {
   Sequence_gstats()
   : gc_(), gem_(0), gme_(0), tem_(0), tme_(0), mm5_(0), mm3_(0)
   {}

   Sequence_gstats(Group_count const& gc, const Mapped_id mid)
   : gc_(gc), gem_(gc[mid]), gme_(0), tem_(gc.copy_count()), tme_(0),
     mm5_(0), mm3_(0)
   {}

   static std::ostream& print_header(std::ostream& os, const char delim) {
      os
    /* ISSUE-41:
      << "group_copies" << delim
      << "group_exact_copies" << delim
      << "group_mismatched_ends" << delim
      << "5'-mismatch" << delim
      << "3'-mismatch" << delim
      << "total_copies" << delim
      << "total_exact_copies" << delim */
      << "exact_copies"
      ;
      return os;
   }

   std::ostream& print(
            std::ostream& os,
            const char delim,
            const Mapped_id mid
   ) const {
      os
    /* ISSUE-41:
      << gc_[mid] << delim
      << gem_ << delim
      << gme_ << delim
      << (int)mm5_ << delim
      << (int)mm3_ << delim
      << gc_.copy_count() << delim
      << tem_ << delim
      << tme_ */
      << gc_[mid]
      ;
      return os;
   }

   Group_count gc_;
   unsigned gem_; 
   unsigned gme_; 
   unsigned tem_; 
   unsigned tme_; 
   unsigned char mm5_; 
   unsigned char mm3_; 
};

class Group_unique {
   typedef boost::container::flat_map<Seq_id, Sequence_gstats> seq_map;
   typedef seq_map::value_type pair;
   typedef std::vector<unsigned> xstats;

public:
   Group_unique(
            Seq_store& ss,
            Read_ginfo_store const& rgs,
            const unsigned n_mids,
            const Mapped_id mid,
            Get_match_length const& gml
   )
   : st_(ss.sequence_map()),
     m_(),
     ss_(ss),
     rgs_(rgs),
     n_mids_(n_mids),
     mid_(mid),
     gml_(gml)
   {
      if( ss.empty() ) return;
      const unsigned shortest = ss_[ss_.by_size().front()].size();
      for(
               unsigned len = ss_[ss_.by_size().back()].size();
               len >= shortest;
               --len
      ) {
         BOOST_FOREACH(const Seq_id sid, ss_.by_size(len)) {
            Group_count gc(n_mids_);
            BOOST_FOREACH(const sub_seq sp, ss_.maps_from(sid)) {
               gc(rgs_[sp.id_].mid_);
            }
            if( gc[mid_] ) find_duplicates(sid, gc);
         }
      }
   }

   unsigned n_reads() const {return m_.size();}

   void write_group_unique(
            File_ostream& fos,
            const unsigned min_duplicates,
            const bool consensus_trim,
        const std::string group_name
   ) const {
      BOOST_FOREACH(pair const& p, m_) {
         const Read_id rid = find_id(p.first);
         if( ! min_duplicates ||  min_duplicates < p.second.gc_[mid_] ) {
            boost::string_ref seq = ss_[p.first].sequence();
            if( consensus_trim ) {
               seq.remove_prefix(p.second.mm5_);
               seq.remove_suffix(p.second.mm3_);
            }
        //std::cout << min_duplicates << ' ' << rgs_[rid].id_str_ << ' ' << p.second.gc_[mid_] << ' ' << p.second.gem_ << '\n';
        std::ostringstream stm;
        stm << p.second.gc_[mid_];
        std::string descr = rgs_[rid].id_str_ + "|MID=" + group_name + "|DUPCOUNT=" + stm.str();
            fos.write(descr, seq);
         }
      }
   }

   void write_unique(
            File_ostream& fos,
            const unsigned min_duplicates
   ) const {
      BOOST_FOREACH(pair const& p, m_) {
         BOOST_ASSERT(p.second.gc_.group_count());
         BOOST_ASSERT(p.second.gc_[mid_]);
         if(
                  p.second.gc_.group_count() > 1 ||
                  ( min_duplicates && p.second.gc_[mid_] <= min_duplicates )
         ) continue;
         const Read_id rid = find_id(p.first);
     std::ostringstream stm;
     stm << p.second.gc_[mid_];
     std::string descr = rgs_[rid].id_str_ + "|DUPCOUNT=" + stm.str();
         fos.write(descr, ss_[p.first].sequence());
      }
   }

   void write_redundancy(std::ostream& os, const char delim) const {
      os << "read_id" << delim;
      Sequence_gstats::print_header(os, delim) << '\n';

      BOOST_FOREACH(pair const& p, m_) {
         const Read_id rid = find_id(p.first);
         os << rgs_[rid].id_str_ << delim;
         p.second.print(os, delim, mid_) << '\n';
      }
   }

   unsigned find_shared(xstats& xs1, xstats& xs2) {
      xs1.resize(n_mids_, 0);
      xs2.resize(n_mids_, 0);
      BOOST_FOREACH(const Seq_id sid, ss_) {
         Group_count gc(n_mids_);
         BOOST_FOREACH(const sub_seq sp, ss_.maps_from(sid)) {
            gc(rgs_[sp.id_].mid_);
         }

         if( gc[mid_] || find_shared(sid, gc) ) {
            transform(
                     xs2.begin(),
                     xs2.end(),
                     gc.v_.begin(),
                     xs2.begin(),
                     std::plus<unsigned>()
            );
         }
      }

      unsigned unique = 0;
      BOOST_FOREACH(pair const& p, m_) {
         bool is_unique = true;
         for(std::size_t n = 0; n != n_mids_; ++n) {
            if( p.second.gc_.v_[n] ) ++xs1[n];
            if( is_unique && n != mid_() && p.second.gc_.v_[n] ) {
               is_unique = false;
            }
         }
         unique += is_unique;
      }
      return unique;
   }

   void group_redundancy_histogram(xstats& h) const {
      BOOST_FOREACH(pair const& p, m_) {
         const std::size_t n = p.second.gc_[mid_];
         if( h.size() <= n ) h.resize(n + 1, 0);
         ++h[n];
      }
   }

   void redundancy_histogram(xstats& h) const {
      BOOST_FOREACH(pair const& p, m_) {
         const std::size_t n = p.second.gc_.copy_count();
         if( h.size() <= n ) h.resize(n + 1, 0);
         ++h[n];
      }
   }

private:
   gdst::Gdst st_;
   seq_map m_;
   Seq_store const& ss_;
   Read_ginfo_store const& rgs_;
   unsigned n_mids_;
   Mapped_id mid_;
   Get_match_length const& gml_;

   void store_stats(
            const unsigned seq_len,
            const unsigned match_pos,
            const unsigned match_len,
            super_seq const& super,
            Group_count const& gc
   ) {
      Sequence_gstats& stats = m_[super.id_];
      stats.gc_(gc);
      if(
               gdst::mismatching_end(
                        match_pos,
                        seq_len,
                        super.pos_,
                        ss_[super.id_].size(),
                        match_len
               )
      ) {
         stats.gme_ += gc[mid_];
         stats.tme_ += gc.copy_count();
      }
      if( stats.mm5_ < match_pos ) stats.mm5_ = match_pos;
      const unsigned match_end = seq_len - match_pos - match_len;
      if( stats.mm3_ < match_end ) stats.mm3_ = match_end;
   }

   bool check_duplicates(
            std::string const& seq,
            Group_count const& gc
   ) {
      BOOST_ASSERT( all(seq, ! Is_ambiguous()) );
      const gdst::Common_subseq cs = st_.find_longest(seq, gml_(seq.size()));
      if( cs.seq_.empty() ) return false;
      BOOST_FOREACH(super_seq const& super, cs.seq_) {
         BOOST_ASSERT(m_.find(super.id_) != m_.end());
         store_stats(seq.size(), cs.start_, cs.len_, super, gc);
      }
      return true;
   }

   bool find_shared(
            const Seq_id sid,
            Group_count const& gc
   ) {
      std::string const& seq = ss_[sid].sequence();
      if( ! all(seq, ! Is_ambiguous()) ) return false;
      return check_duplicates(seq, gc);
   }

   Read_id find_id(const Seq_id sid) const {
      BOOST_FOREACH(const sub_seq sp, ss_.maps_from(sid)) {
         const Read_id rid = sp.id_;
         if( rgs_[rid].mid_ == mid_ ) return rid;
      }
      BOOST_ASSERT(false);
      return Read_id();
   }

   bool is_only_group(const Seq_id sid, const Mapped_id mid) const {
      bool found = false;
      BOOST_FOREACH(const sub_seq sp, ss_.maps_from(sid)) {
         const Read_id rid = sp.id_;
         if( rgs_[rid].mid_ == mid ) found = true;
         if( rgs_[rid].mid_ == mid_ ) return false;
      }
      return found;
   }

   void find_duplicates(
            const Seq_id sid,
            Group_count const& gc
   ) {
      std::string const& seq = ss_[sid].sequence();
      if( all(seq, ! Is_ambiguous()) ) {
         if( check_duplicates(seq, gc) ) return;
         st_.insert(sid);
      }
      BOOST_ASSERT(m_.find(sid) == m_.end());
      m_.emplace(sid, Sequence_gstats(gc, mid_));
   }
};

}//namespace vdj_pipe
#endif /* FIND_SHARED_IMPL_HPP_ */