VCF.h

#ifndef theplu_yat_omic_vcf
#define theplu_yat_omic_vcf

// $Id: VCF.h 3999 2020-10-08 23:22:32Z peter $

/*
  Copyright (C) 2018, 2019, 2020 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library 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.

  The yat library 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 yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include "yat/utility/config_public.h"

#include "yat/utility/Exception.h"
#include "yat/utility/split.h"
#include "yat/utility/utility.h"
#include "yat/utility/yat_assert.h"

#include <boost/cstdint.hpp>

#include <algorithm>
#include <vector>
#include <string>

namespace theplu {
namespace yat {
namespace omic {

  class VCF
  {
  public:
    VCF(void);

    explicit VCF(const std::string&);

    explicit VCF(std::istream& is);

    void chromosome(const std::string& chr);

    const std::string& chromosome(void) const;

    void pos(int32_t p);

    const int32_t& pos(void) const;

    void id(const std::string& id);

    const std::string& id(void) const;

    void ref(const std::string& r);

    void ref(char r);

    const std::string& ref(void) const;

    void alt(const std::string& a);

    void alt(char a);

    const std::string& alt(void) const;

    const std::vector<std::string>& alts(void) const;

    void alts(const std::vector<std::string>& a);

    void alts(std::vector<std::string>&& a);

    unsigned int n_alleles(void) const;

    const std::string& qual(void) const;

    void qual(const std::string& q);

    void qual(unsigned int q);

    const std::string& filter(void) const;

    const std::vector<std::string>& filters(void) const;

    void add_filter(const std::string&);

    void filter(const std::string&);

    class Info
    {
    public:
      void add(const std::string& key);

      template<typename T>
      void add(const std::string& key, const T& value);

      void clear(void);

      size_t count(const std::string& key) const;

      template<typename T>
      void get(T& result, const std::string& key) const;

      void remove(const std::string& key);

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

      void set(std::string&& s);

      const std::string& str(void) const;

    private:
      friend std::ostream& operator<<(std::ostream&, const VCF::Info&);
      friend class VCF;

      // find range [begin, end) such that
      // 1) [begin, end) == key
      // 2) begin is equal to str_.begin() or begin[-1] is ';'
      // 3) end is equal to str_.end() or end[0] == '=' or end[0] == ';'
      std::string::iterator find(const std::string& key);
      // same as find (above) but returns a const_iterator
      std::string::const_iterator cfind(const std::string& key) const;

      void init(std::string& str);
      void init(std::string&& str);

      // For now we store data as raw string (as it's stored in
      // file). This has no overhead and is thus the fastest in case
      // when the info field is not accessed nor modified. For the
      // case when heavily accessed and modified, it might be
      // preferable to also store data in a map<string, string> with a
      // lazy implementation switching between storages depending on
      // usage (i.e. avoid buoilding the map until needed).
      std::string str_;

      // only call if initialised
      void validate(void) const;
    }; // end class Info


    VCF::Info& info(void);

    const VCF::Info& info(void) const;

    class Data
    {
    public:
      Data(void);

      template<typename T>
      void add(const std::string& key, const std::vector<T>& data);

      void clear(void);

      size_t count(const std::string& key) const;

      const std::vector<std::string>& format(void) const;

      template<typename T>
      void get(const std::string& key, std::vector<T>& data) const;

    private:
      friend std::ostream& operator<<(std::ostream&, const VCF::Data&);
      friend class VCF;

      // return true if class has been initialised
      bool initialised(void) const;
      void init(std::vector<std::string>::const_iterator,
                std::vector<std::string>::const_iterator);
      size_t index(const std::string& key) const;
      bool initialised_;
      std::vector<std::string> format_;
      std::vector<std::vector<std::string> > data_;

      // only call if initialised
      void validate(void) const;
    }; // end class Data

    const Data& data(void) const;

    Data& data(void);

    unsigned int n_samples(void) const;

    void n_samples(int n);

  private:
    friend std::ostream& operator<<(std::ostream&, const VCF&);
    std::vector<std::string> vec_;
    int32_t position_;

    Info info_;
    Data data_;

    void init(std::vector<std::string>& vec);
    // update vec_[4] from alts_
    void update_alt(void);
    void print(std::ostream& os) const;

    // if NDEBUG no-op
    // otherwise throw if not valid
    void validate(void) const;

    // private lazy variables
    mutable std::vector<std::string> alts_;
    mutable std::vector<std::string> filters_;

    // using compiler generated copy
    // VCF(const VCF&)
    // VCF& operator=(const VCF&)
    // using compiler generated move
    // VCF(VCF&&)
    // VCF& operator=(VCF&&)
  };

  std::ostream& operator<<(std::ostream& os, const VCF& vcf);

  std::istream& operator>>(std::istream& is, VCF& vcf);

  std::ostream& operator<<(std::ostream& os, const VCF::Info& info);

  std::ostream& operator<<(std::ostream& os, const VCF::Data& data);

  bool is_indel(const VCF&);

  bool is_snv(const VCF&);

  bool is_dnv(const VCF&);

  bool is_mnv(const VCF&);


  // We would prefer to hide this privately with VCF, but explicit
  // template specialization must be in namespace scope.
  //
  namespace detail {

    // class that append a string from T; T might be string, numerical
    // or a vector of above.
    template <typename T>
    struct Appender
    {
      Appender(char delim) : delim_(delim) {}
      void operator()(std::string& val, const T& x) const
      {
        val += utility::convert(x);
      }
    private:
      char delim_;
    };


    template <>
    struct Appender<std::string>
    {
      Appender(char delim) : delim_(delim) {}
      void operator()(std::string& val, const std::string& x) const
      {
        val += x;
      }
    private:
      char delim_;
    };


    template <>
    struct Appender<std::vector<std::string> >
    {
      Appender(char delim) : delim_(delim) {}
      void operator()(std::string& val,const std::vector<std::string>& x) const
      {
        for (size_t i=0; i<x.size(); ++i) {
          if (i)
            val += delim_;
          val += x[i];
        }
      }
    private:
      char delim_;
    };


    template <typename T>
    struct Appender<std::vector<T> >
    {
      Appender(char delim) : delim_(delim) {}
      void operator()(std::string& val, const std::vector<T>& x) const
      {
        for (size_t i=0; i<x.size(); ++i) {
          if (i)
            val += delim_;
          val += utility::convert(x[i]);
        }
      }
    private:
      char delim_;
    };


    // Class that converts a string to T
    template<typename T>
    struct Converter
    {
      void operator()(const std::string& x, T& result) const
      {
        result = utility::convert<T>(x);
      }
    };


    template<>
    struct Converter<std::string>
    {
      void operator()(const std::string& x, std::string& result) const
      {
        result = x;
      }
    };


    template<typename T>
    struct Converter<std::vector<T> >
    {
      void operator()(const std::string& x, std::vector<T>& result) const
      {
        std::vector<std::string> vec;
        utility::split(vec, x, ',');
        result.clear();
        result.reserve(vec.size());
        for (size_t i=0; i<vec.size(); ++i)
          result.push_back(utility::convert<T>(vec[i]));
      }
    };


    template<>
    struct Converter<std::vector<std::string> >
    {
      void
      operator()(const std::string& x, std::vector<std::string>& result) const
      {
        result.clear();
        utility::split(result, x, ',');
      }
    };

  }



  // template implementations

  // VCF::Info

  template<typename T>
  void VCF::Info::add(const std::string& key, const T& value)
  {
    add(key);
    str_ += '=';
    detail::Appender<T> append(',');
    append(str_, value);
  }


  template<typename T>
  void VCF::Info::get(T& result, const std::string& key) const
  {
    detail::Converter<T> converter;
    std::string value;
    std::string k = key + "=";
    size_t begin = str_.find(k);

    if (begin == std::string::npos) {
      converter(value, result);
      return;
    }

    if (begin && str_[begin-1]!=';') {
      k = ";" + key + "=";
      begin = str_.find(k, begin+key.size());
      if (begin == std::string::npos) {
        converter(value, result);
        return;
      }
    }

    begin += k.size();
    int end = str_.find(";", begin);
    value = str_.substr(begin, end-begin);
    converter(value, result);
  }


  // VCF::Data

  template<typename T>
  void VCF::Data::add(const std::string& key, const std::vector<T>& input)
  {
    YAT_ASSERT(key != "");
    format_.push_back(key);
    detail::Appender<T> append(',');
    for (size_t i=0; i<input.size(); ++i) {
      data_[i].push_back("");
      append(data_[i].back(), input[i]);
    }
  }


  template<typename T>
  void VCF::Data::get(const std::string& key, std::vector<T>& data) const
  {
    size_t idx = index(key);
    detail::Converter<T> converter;
    data.resize(data_.size());
    for (size_t i=0; i<data.size(); ++i)
      converter(data_[i][idx], data[i]);
  }


  /*
  template<typename T>
  void VCF::add_data(const std::string& key, const std::vector<T>& val)
  {
    if (format() != "")
      format() += ":";
    format() += key;


    YAT_ASSERT(val.size()+9 == vec_.size());
    for (size_t i=0; i<val.size(); ++i) {
      std::string& data = vec_[i+9];
      if (data != "")
        data += ":";
      append(data, val[i]);
    }
  }
  */
}}}
#endif