variant.h

Go to the documentation of this file.

#ifndef gamgee__variant__guard
#define gamgee__variant__guard

#include "variant_header.h"
#include "individual_field.h"
#include "individual_field_value.h"
#include "shared_field.h"
#include "variant_filters.h"
#include "genotype.h"

#include "../utils/variant_utils.h"

#include "htslib/sam.h"
#include "boost/dynamic_bitset.hpp"

#include <string>
#include <memory>
#include <vector>

namespace gamgee {
enum class DiploidPLGenotype { HOM_REF = 0, HET = 1, HOM_VAR = 2};

class Variant {
 public:
  Variant() = default;                                                                                        
  explicit Variant(const std::shared_ptr<bcf_hdr_t>& header, const std::shared_ptr<bcf1_t>& body) noexcept;   
  Variant(const Variant& other);                                                                              
  Variant& operator=(const Variant& other);                                                                   
  Variant(Variant&& other) = default;                                                                         
  Variant& operator=(Variant&& other) = default;                                                              

  VariantHeader header() const {
    // Avoid a deep copy here by constructing using VariantHeader's internal shared header pointer
    return VariantHeader{m_header.m_header};
  }

  bool missing() const { return m_body == nullptr; }                 

  uint32_t chromosome()         const {return uint32_t(m_body->rid);}                                         
  std::string chromosome_name() const {return header().chromosomes()[chromosome()];}                          
  uint32_t alignment_start()    const {return uint32_t(m_body->pos+1);}                                       
  uint32_t alignment_stop()     const {return uint32_t(m_body->pos + m_body->rlen);}                          
  float    qual()               const {return m_body->qual;}                                                  
  uint32_t n_samples()          const {return uint32_t(m_body->n_sample);}                                    
  uint32_t n_alleles()          const {return uint32_t(m_body->n_allele);}                                    

  std::string id() const;                                                                                     
  std::string ref() const;                                                                                    
  std::vector<std::string> alt() const;                                                                       

  // filter field getter
  VariantFilters filters() const;                                                                             
  bool has_filter(const std::string& filter) const;                                                           

  // individual field getters (a.k.a "format fields")
  IndividualField<Genotype> genotypes() const;                                                                 
  IndividualField<IndividualFieldValue<int32_t>> integer_individual_field(const std::string& tag) const;       
  IndividualField<IndividualFieldValue<float>> float_individual_field(const std::string& tag) const;           
  IndividualField<IndividualFieldValue<std::string>> string_individual_field(const std::string& tag) const;    
  IndividualField<IndividualFieldValue<int32_t>> individual_field_as_integer(const std::string& tag) const;    
  IndividualField<IndividualFieldValue<float>> individual_field_as_float(const std::string& tag) const;        
  IndividualField<IndividualFieldValue<std::string>> individual_field_as_string(const std::string& tag) const; 
  IndividualField<IndividualFieldValue<int32_t>> integer_individual_field(const int32_t index) const;          
  IndividualField<IndividualFieldValue<float>> float_individual_field(const int32_t index) const;              
  IndividualField<IndividualFieldValue<std::string>> string_individual_field(const int32_t index) const;       
  IndividualField<IndividualFieldValue<int32_t>> individual_field_as_integer(const int32_t index) const;       
  IndividualField<IndividualFieldValue<float>> individual_field_as_float(const int32_t index) const;           
  IndividualField<IndividualFieldValue<std::string>> individual_field_as_string(const int32_t index) const;    

  // shared field getters (a.k.a "info fields")
  bool boolean_shared_field(const std::string& tag) const;                       
  SharedField<int32_t> integer_shared_field(const std::string& tag) const;       
  SharedField<float> float_shared_field(const std::string& tag) const;           
  SharedField<std::string> string_shared_field(const std::string& tag) const;    
  SharedField<int32_t> shared_field_as_integer(const std::string& tag) const;    
  SharedField<float> shared_field_as_float(const std::string& tag) const;        
  SharedField<std::string> shared_field_as_string(const std::string& tag) const; 
  bool boolean_shared_field(const int32_t index) const;                          
  SharedField<int32_t> integer_shared_field(const int32_t index) const;          
  SharedField<float> float_shared_field(const int32_t index) const;              
  SharedField<std::string> string_shared_field(const int32_t index) const;       
  SharedField<int32_t> shared_field_as_integer(const int32_t index) const;       
  SharedField<float> shared_field_as_float(const int32_t index) const;           
  SharedField<std::string> shared_field_as_string(const int32_t index) const;    

  template <class VALUE, template<class> class ITER>                                                                                         
  static boost::dynamic_bitset<> select_if(                                                                   
      const ITER<VALUE>& first,         
      const ITER<VALUE>& last,          
      const std::function<bool (const decltype(*first)& value)> pred)  // decltype(*first) will always yield an l-value reference to VALUE (VALUE&). I just added the & here to be explicit about my intention of capturing the parameter by reference.
  {
    const auto n_samples = last - first; 
    auto selected_samples = boost::dynamic_bitset<>(n_samples);
    auto it = first;
    for (auto i = 0; i < n_samples; i++) {
      selected_samples[i] = pred(*it++);
    }
    return selected_samples;
  }

  AlleleMask allele_mask() const;

 private:
  VariantHeader m_header;                                                                        
  std::shared_ptr<bcf1_t> m_body;                                                                

  bcf_fmt_t*  find_individual_field(const std::string& tag) const { return bcf_get_fmt(m_header.m_header.get(), m_body.get(), tag.c_str());  }
  bcf_info_t* find_shared_field(const std::string& tag)     const { return bcf_get_info(m_header.m_header.get(), m_body.get(), tag.c_str()); }
  bcf_fmt_t*  find_individual_field(const uint32_t index) const { return bcf_get_fmt_id(m_body.get(), index); }
  bcf_info_t* find_shared_field(const uint32_t index)     const { return bcf_get_info_id(m_body.get(), index); }
  bool check_field(const int32_t type_field, const int32_t type_value, const int32_t index) const;
  inline AlleleType allele_type_from_difference(const int diff) const;

  template<class FIELD_TYPE, class INDEX_OR_TAG> SharedField<FIELD_TYPE> shared_field_as(const INDEX_OR_TAG& p) const;
  template<class FIELD_TYPE, class INDEX_OR_TAG> IndividualField<IndividualFieldValue<FIELD_TYPE>> individual_field_as(const INDEX_OR_TAG& p) const;

  friend class VariantWriter;
  friend class VariantBuilder; 

  // TODO: remove this friendship and these mutators after Issue #320 is resolved

  friend class ReferenceBlockSplittingVariantIterator;

  inline void set_alignment_start(const int32_t start) { m_body->pos = start - 1; }
  inline void set_alignment_stop(const int32_t end) { m_body->rlen = end - m_body->pos; }

  inline void set_reference_allele(const char* ref, const int32_t ref_length)
  {
    bcf_unpack(m_body.get(), BCF_UN_STR);
    //Try to avoid calling bcf_update_alleles as it causes significant overhead in
    //terms of memory re-allocation etc.
    //If enough space is available, over-write the current reference allele value
    //Seems like an un-necessary hack, but trust me every single such overhead counts
    if(m_body->rlen >= ref_length)
    {
      memcpy(m_body->d.allele[0], ref, ref_length); //memcpy is slightly faster than strcpy
      m_body->d.allele[0][ref_length] = '\0';   //append null - why not include in memcpy? see set_reference_allele with char option
      m_body->d.shared_dirty |= BCF1_DIRTY_ALS;
    }
    else
    {
      //expensive - causes significant reallocs within htslib
      //hacks to interpret const ptrs to non-const. Again, saves some overhead
      //No need to free d.allele[0] because it points to shared string within bcf1_1.shared
      m_body->d.allele[0] = const_cast<char*>(ref);
      //Re-use same d.allele char** as update_alleles writes to different region of memory anyway
      bcf_update_alleles(const_cast<const bcf_hdr_t*>(m_header.m_header.get()), m_body.get(), const_cast<const char**>(m_body->d.allele), m_body->n_allele);
    }
  }
  inline void set_reference_allele(const char* ref)  { set_reference_allele(ref, static_cast<int32_t>(strlen(ref))); }
  inline void set_reference_allele(const char ref_base)  { set_reference_allele(&ref_base, 1); }
};

}  // end of namespace

#endif /* defined(gamgee__variant__guard) */