yat/omic/VCF.h

Code
Comments
Other
Rev Date Author Line
3747 14 Aug 18 peter 1 #ifndef theplu_yat_omic_vcf
3747 14 Aug 18 peter 2 #define theplu_yat_omic_vcf
3747 14 Aug 18 peter 3
3747 14 Aug 18 peter 4 // $Id$
3747 14 Aug 18 peter 5
3747 14 Aug 18 peter 6 /*
4316 17 Feb 23 peter 7   Copyright (C) 2018, 2019, 2020, 2022, 2023 Peter Johansson
3747 14 Aug 18 peter 8
3747 14 Aug 18 peter 9   This file is part of the yat library, http://dev.thep.lu.se/yat
3747 14 Aug 18 peter 10
3747 14 Aug 18 peter 11   The yat library is free software; you can redistribute it and/or
3747 14 Aug 18 peter 12   modify it under the terms of the GNU General Public License as
3747 14 Aug 18 peter 13   published by the Free Software Foundation; either version 3 of the
3747 14 Aug 18 peter 14   License, or (at your option) any later version.
3747 14 Aug 18 peter 15
3747 14 Aug 18 peter 16   The yat library is distributed in the hope that it will be useful,
3747 14 Aug 18 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
3747 14 Aug 18 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3747 14 Aug 18 peter 19   General Public License for more details.
3747 14 Aug 18 peter 20
3747 14 Aug 18 peter 21   You should have received a copy of the GNU General Public License
3755 17 Oct 18 peter 22   along with yat. If not, see <http://www.gnu.org/licenses/>.
3747 14 Aug 18 peter 23 */
3747 14 Aug 18 peter 24
3747 14 Aug 18 peter 25 #include "yat/utility/config_public.h"
3747 14 Aug 18 peter 26
3747 14 Aug 18 peter 27 #include "yat/utility/Exception.h"
3747 14 Aug 18 peter 28 #include "yat/utility/split.h"
3747 14 Aug 18 peter 29 #include "yat/utility/utility.h"
3747 14 Aug 18 peter 30 #include "yat/utility/yat_assert.h"
3747 14 Aug 18 peter 31
3747 14 Aug 18 peter 32 #include <boost/cstdint.hpp>
3747 14 Aug 18 peter 33
3747 14 Aug 18 peter 34 #include <algorithm>
4145 24 Feb 22 peter 35 #include <map>
3747 14 Aug 18 peter 36 #include <vector>
3747 14 Aug 18 peter 37 #include <string>
3747 14 Aug 18 peter 38
3747 14 Aug 18 peter 39 namespace theplu {
3747 14 Aug 18 peter 40 namespace yat {
3747 14 Aug 18 peter 41 namespace omic {
3747 14 Aug 18 peter 42
3747 14 Aug 18 peter 43   /**
3747 14 Aug 18 peter 44      Class to parse and modify an entry in a VcfFile (Variant Call
3747 14 Aug 18 peter 45      Format) as defined in https://samtools.github.io/hts-specs/VCFv4.2.pdf
3747 14 Aug 18 peter 46
3747 14 Aug 18 peter 47      \since new in yat 0.16
3747 14 Aug 18 peter 48    */
3747 14 Aug 18 peter 49   class VCF
3747 14 Aug 18 peter 50   {
3747 14 Aug 18 peter 51   public:
3747 14 Aug 18 peter 52     /**
3747 14 Aug 18 peter 53        \brief default constructor
3747 14 Aug 18 peter 54      */
3747 14 Aug 18 peter 55     VCF(void);
3747 14 Aug 18 peter 56
3747 14 Aug 18 peter 57     /**
4144 23 Feb 22 peter 58        Create a VCF representing \c n samples.
4144 23 Feb 22 peter 59
4144 23 Feb 22 peter 60        \since New in yat 0.20
4144 23 Feb 22 peter 61      */
4144 23 Feb 22 peter 62     explicit VCF(size_t n);
4144 23 Feb 22 peter 63
4144 23 Feb 22 peter 64     /**
3747 14 Aug 18 peter 65        \brief Construct entry from a std::string
3747 14 Aug 18 peter 66      */
3747 14 Aug 18 peter 67     explicit VCF(const std::string&);
3747 14 Aug 18 peter 68
3747 14 Aug 18 peter 69     /**
3747 14 Aug 18 peter 70        Equivalent with
3747 14 Aug 18 peter 71        \code
3747 14 Aug 18 peter 72        std::string line;
3747 14 Aug 18 peter 73        std::getline(is, line);
3747 14 Aug 18 peter 74        VCF(line);
3747 14 Aug 18 peter 75        \endcode
3747 14 Aug 18 peter 76      */
3747 14 Aug 18 peter 77     explicit VCF(std::istream& is);
3747 14 Aug 18 peter 78
3747 14 Aug 18 peter 79     /**
3747 14 Aug 18 peter 80        \brief set chromosome name
3747 14 Aug 18 peter 81      */
3747 14 Aug 18 peter 82     void chromosome(const std::string& chr);
3747 14 Aug 18 peter 83
3747 14 Aug 18 peter 84     /**
3747 14 Aug 18 peter 85        \return Chromosome name
3747 14 Aug 18 peter 86      */
3747 14 Aug 18 peter 87     const std::string& chromosome(void) const;
3747 14 Aug 18 peter 88
3747 14 Aug 18 peter 89     /**
3747 14 Aug 18 peter 90        \brief set position
3747 14 Aug 18 peter 91      */
3747 14 Aug 18 peter 92     void pos(int32_t p);
3747 14 Aug 18 peter 93
3747 14 Aug 18 peter 94     /**
3747 14 Aug 18 peter 95        \return position
3747 14 Aug 18 peter 96      */
3747 14 Aug 18 peter 97     const int32_t& pos(void) const;
3747 14 Aug 18 peter 98
3747 14 Aug 18 peter 99     /**
3747 14 Aug 18 peter 100        \brief Variant identifier
3747 14 Aug 18 peter 101
3747 14 Aug 18 peter 102        Usually the dbSNP rsID
3747 14 Aug 18 peter 103      */
3747 14 Aug 18 peter 104     void id(const std::string& id);
3747 14 Aug 18 peter 105
3747 14 Aug 18 peter 106     /**
3747 14 Aug 18 peter 107        \return Variant identifier
3747 14 Aug 18 peter 108      */
3747 14 Aug 18 peter 109     const std::string& id(void) const;
3747 14 Aug 18 peter 110
3747 14 Aug 18 peter 111     /**
3747 14 Aug 18 peter 112        \brief set refererence sequence at position involved in variant
3747 14 Aug 18 peter 113      */
3747 14 Aug 18 peter 114     void ref(const std::string& r);
3747 14 Aug 18 peter 115
3747 14 Aug 18 peter 116     /**
3747 14 Aug 18 peter 117        \brief set refererence sequence at position involved in variant
3747 14 Aug 18 peter 118      */
3747 14 Aug 18 peter 119     void ref(char r);
3747 14 Aug 18 peter 120
3747 14 Aug 18 peter 121     /**
3747 14 Aug 18 peter 122        \return reference sequence at position
3747 14 Aug 18 peter 123      */
3747 14 Aug 18 peter 124     const std::string& ref(void) const;
3747 14 Aug 18 peter 125
3747 14 Aug 18 peter 126     /**
3747 14 Aug 18 peter 127        \brief set alternative sequence(s)
3747 14 Aug 18 peter 128
3747 14 Aug 18 peter 129        \param a is a comma-delimited list of sequences
3747 14 Aug 18 peter 130      */
3747 14 Aug 18 peter 131     void alt(const std::string& a);
3747 14 Aug 18 peter 132
3747 14 Aug 18 peter 133     /**
3747 14 Aug 18 peter 134        \brief set alternative sequence(s)
3747 14 Aug 18 peter 135
3747 14 Aug 18 peter 136        \param a is a comma-delimited list of sequences
3747 14 Aug 18 peter 137      */
3747 14 Aug 18 peter 138     void alt(char a);
3747 14 Aug 18 peter 139
3747 14 Aug 18 peter 140     /**
3747 14 Aug 18 peter 141        \return alternative sequence(s)
3747 14 Aug 18 peter 142      */
3747 14 Aug 18 peter 143     const std::string& alt(void) const;
3747 14 Aug 18 peter 144
3747 14 Aug 18 peter 145     /**
3747 14 Aug 18 peter 146        \return alternative sequences
3747 14 Aug 18 peter 147      */
3747 14 Aug 18 peter 148     const std::vector<std::string>& alts(void) const;
3747 14 Aug 18 peter 149
3747 14 Aug 18 peter 150     /**
3747 14 Aug 18 peter 151        \brief set alternative alleles
3747 14 Aug 18 peter 152     */
3747 14 Aug 18 peter 153     void alts(const std::vector<std::string>& a);
3747 14 Aug 18 peter 154
3747 14 Aug 18 peter 155     /**
3747 14 Aug 18 peter 156        \brief set alternative alleles
3747 14 Aug 18 peter 157     */
3747 14 Aug 18 peter 158     void alts(std::vector<std::string>&& a);
3747 14 Aug 18 peter 159
3747 14 Aug 18 peter 160     /**
3747 14 Aug 18 peter 161        \return number of alternative sequences
3747 14 Aug 18 peter 162      */
3747 14 Aug 18 peter 163     unsigned int n_alleles(void) const;
3747 14 Aug 18 peter 164
3747 14 Aug 18 peter 165     /**
3747 14 Aug 18 peter 166        \brief phred-scaled probability that all samples are wildtype
3747 14 Aug 18 peter 167      */
3747 14 Aug 18 peter 168     const std::string& qual(void) const;
3747 14 Aug 18 peter 169
3747 14 Aug 18 peter 170     /**
3747 14 Aug 18 peter 171        \brief set quality
3747 14 Aug 18 peter 172      */
3747 14 Aug 18 peter 173     void qual(const std::string& q);
3747 14 Aug 18 peter 174
3747 14 Aug 18 peter 175     /**
3747 14 Aug 18 peter 176        \brief set quality
3747 14 Aug 18 peter 177      */
3747 14 Aug 18 peter 178     void qual(unsigned int q);
3747 14 Aug 18 peter 179
3747 14 Aug 18 peter 180     /**
3747 14 Aug 18 peter 181        \return semicolon delimited list of filters that the variant fails
3747 14 Aug 18 peter 182      */
3747 14 Aug 18 peter 183     const std::string& filter(void) const;
3747 14 Aug 18 peter 184
3747 14 Aug 18 peter 185     /**
3747 14 Aug 18 peter 186        \return filters
3747 14 Aug 18 peter 187      */
3747 14 Aug 18 peter 188     const std::vector<std::string>& filters(void) const;
3747 14 Aug 18 peter 189
3747 14 Aug 18 peter 190     /**
3747 14 Aug 18 peter 191        \brief append list of filters
3747 14 Aug 18 peter 192      */
3747 14 Aug 18 peter 193     void add_filter(const std::string&);
3747 14 Aug 18 peter 194
3747 14 Aug 18 peter 195     /**
3747 14 Aug 18 peter 196        \brief set filters
3747 14 Aug 18 peter 197      */
3747 14 Aug 18 peter 198     void filter(const std::string&);
3747 14 Aug 18 peter 199
3747 14 Aug 18 peter 200     /**
4144 23 Feb 22 peter 201        Create a VCF with subset of samples in \c other defined by \c index.
4144 23 Feb 22 peter 202
4144 23 Feb 22 peter 203        If present, INFO fields AC and AN are recalculated to match the
4144 23 Feb 22 peter 204        subset.
4144 23 Feb 22 peter 205
4144 23 Feb 22 peter 206        \see VCF::Data.subset(const std::vector<size_t>&)
4144 23 Feb 22 peter 207
4144 23 Feb 22 peter 208        \since New in yat 0.20
4144 23 Feb 22 peter 209      */
4144 23 Feb 22 peter 210     void subset(const std::vector<size_t>& index);
4144 23 Feb 22 peter 211
4144 23 Feb 22 peter 212     /**
3747 14 Aug 18 peter 213        \brief class handling the INFO field in an VCF entry
3747 14 Aug 18 peter 214      */
3747 14 Aug 18 peter 215     class Info
3747 14 Aug 18 peter 216     {
3747 14 Aug 18 peter 217     public:
4145 24 Feb 22 peter 218       /// default constructor
4145 24 Feb 22 peter 219       Info(void) = default;
4145 24 Feb 22 peter 220
3747 14 Aug 18 peter 221       /**
3747 14 Aug 18 peter 222          Add entry \a key to semicolon-delimited list of entries.
4316 17 Feb 23 peter 223
4316 17 Feb 23 peter 224          If \a key is already present, the behaviour is undefined.
3747 14 Aug 18 peter 225        */
3747 14 Aug 18 peter 226       void add(const std::string& key);
3747 14 Aug 18 peter 227
3747 14 Aug 18 peter 228       /**
3747 14 Aug 18 peter 229          T is either string or numeric or a std::vector with strings
3747 14 Aug 18 peter 230          or numerics.
4316 17 Feb 23 peter 231
4316 17 Feb 23 peter 232          If \a key is already present, the behaviour is undefined.
3747 14 Aug 18 peter 233        */
3747 14 Aug 18 peter 234       template<typename T>
3747 14 Aug 18 peter 235       void add(const std::string& key, const T& value);
3747 14 Aug 18 peter 236
3747 14 Aug 18 peter 237       /**
3747 14 Aug 18 peter 238          Clear Info field.
3747 14 Aug 18 peter 239        */
3747 14 Aug 18 peter 240       void clear(void);
3747 14 Aug 18 peter 241
3747 14 Aug 18 peter 242       /**
3790 04 Apr 19 peter 243          \return 1 if key exists
3790 04 Apr 19 peter 244        */
3790 04 Apr 19 peter 245       size_t count(const std::string& key) const;
3790 04 Apr 19 peter 246
3790 04 Apr 19 peter 247       /**
4316 17 Feb 23 peter 248         get value corresponding to \a key. If \a key is not present,
4316 17 Feb 23 peter 249         result is as value is empty string.
3823 16 Jul 19 peter 250
3823 16 Jul 19 peter 251         T is string, numeric, vector<string> or vector<numeric>
3790 04 Apr 19 peter 252        */
3790 04 Apr 19 peter 253       template<typename T>
3790 04 Apr 19 peter 254       void get(T& result, const std::string& key) const;
3790 04 Apr 19 peter 255
3790 04 Apr 19 peter 256       /**
3747 14 Aug 18 peter 257          \brief remove Info entry with key \a key
3747 14 Aug 18 peter 258        */
3747 14 Aug 18 peter 259       void remove(const std::string& key);
3747 14 Aug 18 peter 260
3747 14 Aug 18 peter 261       /**
3747 14 Aug 18 peter 262          Set info to \a s
3747 14 Aug 18 peter 263        */
3747 14 Aug 18 peter 264       void set(const std::string& s);
3747 14 Aug 18 peter 265
3747 14 Aug 18 peter 266       /**
3747 14 Aug 18 peter 267          Set info to \a s
3747 14 Aug 18 peter 268        */
3747 14 Aug 18 peter 269       void set(std::string&& s);
3747 14 Aug 18 peter 270
3747 14 Aug 18 peter 271       /**
3747 14 Aug 18 peter 272          \return info as a string
3747 14 Aug 18 peter 273        */
3747 14 Aug 18 peter 274       const std::string& str(void) const;
3747 14 Aug 18 peter 275
4145 24 Feb 22 peter 276       /**
4316 17 Feb 23 peter 277          If \a key is absent, equivalent to call add(key, value). If
4316 17 Feb 23 peter 278          \a key present, erase its value and set to \a value.
4316 17 Feb 23 peter 279
4145 24 Feb 22 peter 280          T is either string or numeric or a std::vector with strings
4145 24 Feb 22 peter 281          or numerics.
4171 10 May 22 peter 282
4171 10 May 22 peter 283          \since New in yat 0.20
4145 24 Feb 22 peter 284        */
4145 24 Feb 22 peter 285       template<typename T>
4145 24 Feb 22 peter 286       void set(const std::string& key, const T& value);
4145 24 Feb 22 peter 287
3747 14 Aug 18 peter 288     private:
3747 14 Aug 18 peter 289       friend std::ostream& operator<<(std::ostream&, const VCF::Info&);
3747 14 Aug 18 peter 290       friend class VCF;
3747 14 Aug 18 peter 291
4145 24 Feb 22 peter 292       explicit Info(std::string&& str);
4145 24 Feb 22 peter 293       // This is the raw string as passed in the constructor. keys_
4145 24 Feb 22 peter 294       // and map_ are populated lazily and when those containers are
4145 24 Feb 22 peter 295       // alteredm str_ is cleared to indicate that it's not uptodate
4145 24 Feb 22 peter 296       std::string str_;
3790 04 Apr 19 peter 297
4145 24 Feb 22 peter 298       bool indexed(void) const;
4145 24 Feb 22 peter 299       // This function populates map_ and keys_. This is done lazily,
4145 24 Feb 22 peter 300       // i.e., it's not done untile map and vector are needed, thus
4145 24 Feb 22 peter 301       // index is called within const functions as well and index() is
4145 24 Feb 22 peter 302       // declared const, which can be thought of it doesn't change the
4145 24 Feb 22 peter 303       // object per se.
4145 24 Feb 22 peter 304       void index(void) const;
4145 24 Feb 22 peter 305       std::map<std::string, std::string> map_;
4146 25 Feb 22 peter 306       // We store the keys in a vector to keep track of the order
4145 24 Feb 22 peter 307       std::vector<std::string> keys_;
3747 14 Aug 18 peter 308
3747 14 Aug 18 peter 309       // only call if initialised
3747 14 Aug 18 peter 310       void validate(void) const;
3747 14 Aug 18 peter 311     }; // end class Info
3747 14 Aug 18 peter 312
3747 14 Aug 18 peter 313
3747 14 Aug 18 peter 314     /**
3747 14 Aug 18 peter 315        \brief set info
3747 14 Aug 18 peter 316      */
3747 14 Aug 18 peter 317     VCF::Info& info(void);
3747 14 Aug 18 peter 318
3747 14 Aug 18 peter 319     /**
3747 14 Aug 18 peter 320        \return container holiding information in INFO field (column 8)
3747 14 Aug 18 peter 321      */
3747 14 Aug 18 peter 322     const VCF::Info& info(void) const;
3747 14 Aug 18 peter 323
3747 14 Aug 18 peter 324     /**
3747 14 Aug 18 peter 325        \brief class holiding data stored in FORMAT and data fields
3747 14 Aug 18 peter 326      */
3747 14 Aug 18 peter 327     class Data
3747 14 Aug 18 peter 328     {
3747 14 Aug 18 peter 329     public:
3747 14 Aug 18 peter 330       /// Default constructor
3747 14 Aug 18 peter 331       Data(void);
3747 14 Aug 18 peter 332
4144 23 Feb 22 peter 333       /**
4144 23 Feb 22 peter 334          Create Data representing \c n samples
4144 23 Feb 22 peter 335
4144 23 Feb 22 peter 336          \since New in yat 0.20
4144 23 Feb 22 peter 337        */
4144 23 Feb 22 peter 338       explicit Data(size_t n);
4144 23 Feb 22 peter 339
3747 14 Aug 18 peter 340       /// T is string, numeric or vector string or numeric
3747 14 Aug 18 peter 341       template<typename T>
3747 14 Aug 18 peter 342       void add(const std::string& key, const std::vector<T>& data);
3747 14 Aug 18 peter 343
3747 14 Aug 18 peter 344       /**
3747 14 Aug 18 peter 345          Set FORMAT and all data fields to empty string. Number of
3747 14 Aug 18 peter 346          samples remains unmodified.
3747 14 Aug 18 peter 347        */
3747 14 Aug 18 peter 348       void clear(void);
3747 14 Aug 18 peter 349
3747 14 Aug 18 peter 350       /**
3790 04 Apr 19 peter 351          \return 1 if key exists
3790 04 Apr 19 peter 352        */
3790 04 Apr 19 peter 353       size_t count(const std::string& key) const;
3790 04 Apr 19 peter 354
3790 04 Apr 19 peter 355       /**
3747 14 Aug 18 peter 356          \c FORMAT
3747 14 Aug 18 peter 357
3747 14 Aug 18 peter 358          \return vector of format in the data fields
3747 14 Aug 18 peter 359       */
3747 14 Aug 18 peter 360       const std::vector<std::string>& format(void) const;
3747 14 Aug 18 peter 361
3747 14 Aug 18 peter 362       /// T is string, numeric or vector string or numeric
3747 14 Aug 18 peter 363       template<typename T>
3747 14 Aug 18 peter 364       void get(const std::string& key, std::vector<T>& data) const;
3747 14 Aug 18 peter 365
4144 23 Feb 22 peter 366       /**
4145 24 Feb 22 peter 367          \return number of samples
4145 24 Feb 22 peter 368
4145 24 Feb 22 peter 369          \since New in yat 0.20
4145 24 Feb 22 peter 370        */
4145 24 Feb 22 peter 371       size_t n_samples(void) const;
4145 24 Feb 22 peter 372
4145 24 Feb 22 peter 373       /**
4170 09 May 22 peter 374          T is string, numeric or vector string or numeric
4170 09 May 22 peter 375
4170 09 May 22 peter 376          If key exists in FORMAT, replace data with \c data. Otherwise
4170 09 May 22 peter 377          add data with key; see data(key, data).
4170 09 May 22 peter 378
4170 09 May 22 peter 379          \since New in yat 0.20
4170 09 May 22 peter 380       */
4170 09 May 22 peter 381       template<typename T>
4170 09 May 22 peter 382       void set(const std::string& key, const std::vector<T>& data);
4170 09 May 22 peter 383
4170 09 May 22 peter 384       /**
4144 23 Feb 22 peter 385          Construct Data representing a subset of samples. The object
4144 23 Feb 22 peter 386          will have index.size() samples; if index[0] is k, the first
4144 23 Feb 22 peter 387          sample correspond to the (k+1)th sample in the old object.
4144 23 Feb 22 peter 388
4144 23 Feb 22 peter 389          \since New in yat 0.20
4144 23 Feb 22 peter 390        */
4144 23 Feb 22 peter 391       void subset(const std::vector<size_t>& index);
4144 23 Feb 22 peter 392
3747 14 Aug 18 peter 393     private:
3747 14 Aug 18 peter 394       friend std::ostream& operator<<(std::ostream&, const VCF::Data&);
3747 14 Aug 18 peter 395       friend class VCF;
3747 14 Aug 18 peter 396
3747 14 Aug 18 peter 397       // return true if class has been initialised
3747 14 Aug 18 peter 398       bool initialised(void) const;
3747 14 Aug 18 peter 399       void init(std::vector<std::string>::const_iterator,
3747 14 Aug 18 peter 400                 std::vector<std::string>::const_iterator);
3747 14 Aug 18 peter 401       size_t index(const std::string& key) const;
3747 14 Aug 18 peter 402       bool initialised_;
3747 14 Aug 18 peter 403       std::vector<std::string> format_;
3747 14 Aug 18 peter 404       std::vector<std::vector<std::string> > data_;
3747 14 Aug 18 peter 405
3747 14 Aug 18 peter 406       // only call if initialised
3747 14 Aug 18 peter 407       void validate(void) const;
3747 14 Aug 18 peter 408     }; // end class Data
3747 14 Aug 18 peter 409
3747 14 Aug 18 peter 410     /**
3747 14 Aug 18 peter 411        \brief access data
3747 14 Aug 18 peter 412     */
3747 14 Aug 18 peter 413     const Data& data(void) const;
3747 14 Aug 18 peter 414
3747 14 Aug 18 peter 415     /**
3747 14 Aug 18 peter 416        \brief access data
3747 14 Aug 18 peter 417     */
3747 14 Aug 18 peter 418     Data& data(void);
3747 14 Aug 18 peter 419
3747 14 Aug 18 peter 420     /**
3747 14 Aug 18 peter 421        \return number of samples
3747 14 Aug 18 peter 422      */
3747 14 Aug 18 peter 423     unsigned int n_samples(void) const;
3747 14 Aug 18 peter 424
3747 14 Aug 18 peter 425     /**
3747 14 Aug 18 peter 426        \brief set number of samples to \a n
3747 14 Aug 18 peter 427      */
3747 14 Aug 18 peter 428     void n_samples(int n);
3747 14 Aug 18 peter 429
3747 14 Aug 18 peter 430   private:
3747 14 Aug 18 peter 431     friend std::ostream& operator<<(std::ostream&, const VCF&);
3747 14 Aug 18 peter 432     std::vector<std::string> vec_;
3747 14 Aug 18 peter 433     int32_t position_;
3747 14 Aug 18 peter 434
3747 14 Aug 18 peter 435     Info info_;
3747 14 Aug 18 peter 436     Data data_;
3747 14 Aug 18 peter 437
3747 14 Aug 18 peter 438     void init(std::vector<std::string>& vec);
3747 14 Aug 18 peter 439     // update vec_[4] from alts_
3747 14 Aug 18 peter 440     void update_alt(void);
3747 14 Aug 18 peter 441     void print(std::ostream& os) const;
3747 14 Aug 18 peter 442
4145 24 Feb 22 peter 443     void update_AC_AN(void);
4145 24 Feb 22 peter 444
3747 14 Aug 18 peter 445     // if NDEBUG no-op
3747 14 Aug 18 peter 446     // otherwise throw if not valid
3747 14 Aug 18 peter 447     void validate(void) const;
3747 14 Aug 18 peter 448
3747 14 Aug 18 peter 449     // private lazy variables
3747 14 Aug 18 peter 450     mutable std::vector<std::string> alts_;
3747 14 Aug 18 peter 451     mutable std::vector<std::string> filters_;
3747 14 Aug 18 peter 452
3772 26 Oct 18 peter 453     // using compiler generated copy
3772 26 Oct 18 peter 454     // VCF(const VCF&)
3772 26 Oct 18 peter 455     // VCF& operator=(const VCF&)
3772 26 Oct 18 peter 456     // using compiler generated move
3772 26 Oct 18 peter 457     // VCF(VCF&&)
3772 26 Oct 18 peter 458     // VCF& operator=(VCF&&)
3747 14 Aug 18 peter 459   };
3747 14 Aug 18 peter 460
3747 14 Aug 18 peter 461   /**
3747 14 Aug 18 peter 462      \brief output operator for VCF
3747 14 Aug 18 peter 463      \relates VCF
3747 14 Aug 18 peter 464
3747 14 Aug 18 peter 465      \return stream
3747 14 Aug 18 peter 466    */
3747 14 Aug 18 peter 467   std::ostream& operator<<(std::ostream& os, const VCF& vcf);
3747 14 Aug 18 peter 468
3747 14 Aug 18 peter 469   /**
3747 14 Aug 18 peter 470      \brief input operator for VCF
3772 26 Oct 18 peter 471
3772 26 Oct 18 peter 472      If \a is is good and next character is not EOF, a VCF is read
3772 26 Oct 18 peter 473      into \a vcf as if it was created with the istream constructor.
3772 26 Oct 18 peter 474
3747 14 Aug 18 peter 475      \relates VCF
3747 14 Aug 18 peter 476
3747 14 Aug 18 peter 477      \return stream
3747 14 Aug 18 peter 478    */
3747 14 Aug 18 peter 479   std::istream& operator>>(std::istream& is, VCF& vcf);
3747 14 Aug 18 peter 480
3747 14 Aug 18 peter 481   /**
3747 14 Aug 18 peter 482      \brief output operator for VCF::Info
3747 14 Aug 18 peter 483      \relates VCF::Info
3747 14 Aug 18 peter 484
3747 14 Aug 18 peter 485      \return stream
3747 14 Aug 18 peter 486    */
3747 14 Aug 18 peter 487   std::ostream& operator<<(std::ostream& os, const VCF::Info& info);
3747 14 Aug 18 peter 488
3747 14 Aug 18 peter 489   /**
3747 14 Aug 18 peter 490      \brief output operator for VCF::Data
3747 14 Aug 18 peter 491      \relates VCF::Data
3747 14 Aug 18 peter 492
3747 14 Aug 18 peter 493      \return stream
3747 14 Aug 18 peter 494    */
3747 14 Aug 18 peter 495   std::ostream& operator<<(std::ostream& os, const VCF::Data& data);
3747 14 Aug 18 peter 496
3747 14 Aug 18 peter 497   /**
3747 14 Aug 18 peter 498      \return true if at least one variant is an indel, i.e., if at
3747 14 Aug 18 peter 499      least one alt allele is not the same size as ref.
3920 31 May 20 peter 500
3920 31 May 20 peter 501      \note ALT being an angle-bracketed ID string is ignored
3747 14 Aug 18 peter 502    */
3747 14 Aug 18 peter 503   bool is_indel(const VCF&);
3747 14 Aug 18 peter 504
3747 14 Aug 18 peter 505   /**
3747 14 Aug 18 peter 506      \return true if at least one variant is a single-nucleotide
3747 14 Aug 18 peter 507      variant (SNV), i.e., ref is 1bp and at least one alt allele is
3920 31 May 20 peter 508      1bp (but not '*').
3920 31 May 20 peter 509
3920 31 May 20 peter 510      \note ALT being an angle-bracketed ID string is ignored
3747 14 Aug 18 peter 511    */
3747 14 Aug 18 peter 512   bool is_snv(const VCF&);
3747 14 Aug 18 peter 513
3747 14 Aug 18 peter 514   /**
3747 14 Aug 18 peter 515      \return true if at least one variant is a double-nucleotide
3747 14 Aug 18 peter 516      variant (DNV), i.e., ref is 2bp and at least one alt allele is
3747 14 Aug 18 peter 517      2bp.
3920 31 May 20 peter 518
3920 31 May 20 peter 519      \note ALT being an angle-bracketed ID string is ignored
3747 14 Aug 18 peter 520    */
3747 14 Aug 18 peter 521   bool is_dnv(const VCF&);
3747 14 Aug 18 peter 522
3747 14 Aug 18 peter 523   /**
3747 14 Aug 18 peter 524      \return true if at least one variant is a multi-nucleotide
3747 14 Aug 18 peter 525      variant (MNV), i.e., ref is >2bp and at least one alt allele is
3747 14 Aug 18 peter 526      same length as ref.
3920 31 May 20 peter 527
3920 31 May 20 peter 528      \note ALT being an angle-bracketed ID string is ignored
3747 14 Aug 18 peter 529    */
3747 14 Aug 18 peter 530   bool is_mnv(const VCF&);
3747 14 Aug 18 peter 531
3747 14 Aug 18 peter 532
3747 14 Aug 18 peter 533   // We would prefer to hide this privately with VCF, but explicit
3747 14 Aug 18 peter 534   // template specialization must be in namespace scope.
3747 14 Aug 18 peter 535   //
3747 14 Aug 18 peter 536   /// \cond IGNORE_DOXYGEN
3747 14 Aug 18 peter 537   namespace detail {
3747 14 Aug 18 peter 538
3747 14 Aug 18 peter 539     // class that append a string from T; T might be string, numerical
3747 14 Aug 18 peter 540     // or a vector of above.
3747 14 Aug 18 peter 541     template <typename T>
3747 14 Aug 18 peter 542     struct Appender
3747 14 Aug 18 peter 543     {
3747 14 Aug 18 peter 544       Appender(char delim) : delim_(delim) {}
3747 14 Aug 18 peter 545       void operator()(std::string& val, const T& x) const
3747 14 Aug 18 peter 546       {
3747 14 Aug 18 peter 547         val += utility::convert(x);
3747 14 Aug 18 peter 548       }
3747 14 Aug 18 peter 549     private:
3747 14 Aug 18 peter 550       char delim_;
3747 14 Aug 18 peter 551     };
3747 14 Aug 18 peter 552
3747 14 Aug 18 peter 553
3747 14 Aug 18 peter 554     template <>
3747 14 Aug 18 peter 555     struct Appender<std::string>
3747 14 Aug 18 peter 556     {
3747 14 Aug 18 peter 557       Appender(char delim) : delim_(delim) {}
3747 14 Aug 18 peter 558       void operator()(std::string& val, const std::string& x) const
3747 14 Aug 18 peter 559       {
3747 14 Aug 18 peter 560         val += x;
3747 14 Aug 18 peter 561       }
3747 14 Aug 18 peter 562     private:
3747 14 Aug 18 peter 563       char delim_;
3747 14 Aug 18 peter 564     };
3747 14 Aug 18 peter 565
3747 14 Aug 18 peter 566
3747 14 Aug 18 peter 567     template <>
3747 14 Aug 18 peter 568     struct Appender<std::vector<std::string> >
3747 14 Aug 18 peter 569     {
3747 14 Aug 18 peter 570       Appender(char delim) : delim_(delim) {}
3747 14 Aug 18 peter 571       void operator()(std::string& val,const std::vector<std::string>& x) const
3747 14 Aug 18 peter 572       {
3747 14 Aug 18 peter 573         for (size_t i=0; i<x.size(); ++i) {
3747 14 Aug 18 peter 574           if (i)
3747 14 Aug 18 peter 575             val += delim_;
3747 14 Aug 18 peter 576           val += x[i];
3747 14 Aug 18 peter 577         }
3747 14 Aug 18 peter 578       }
3747 14 Aug 18 peter 579     private:
3747 14 Aug 18 peter 580       char delim_;
3747 14 Aug 18 peter 581     };
3747 14 Aug 18 peter 582
3747 14 Aug 18 peter 583
3747 14 Aug 18 peter 584     template <typename T>
3747 14 Aug 18 peter 585     struct Appender<std::vector<T> >
3747 14 Aug 18 peter 586     {
3747 14 Aug 18 peter 587       Appender(char delim) : delim_(delim) {}
3747 14 Aug 18 peter 588       void operator()(std::string& val, const std::vector<T>& x) const
3747 14 Aug 18 peter 589       {
3747 14 Aug 18 peter 590         for (size_t i=0; i<x.size(); ++i) {
3747 14 Aug 18 peter 591           if (i)
3747 14 Aug 18 peter 592             val += delim_;
3747 14 Aug 18 peter 593           val += utility::convert(x[i]);
3747 14 Aug 18 peter 594         }
3747 14 Aug 18 peter 595       }
3747 14 Aug 18 peter 596     private:
3747 14 Aug 18 peter 597       char delim_;
3747 14 Aug 18 peter 598     };
3747 14 Aug 18 peter 599
3747 14 Aug 18 peter 600
3747 14 Aug 18 peter 601     // Class that converts a string to T
3747 14 Aug 18 peter 602     template<typename T>
3747 14 Aug 18 peter 603     struct Converter
3747 14 Aug 18 peter 604     {
3747 14 Aug 18 peter 605       void operator()(const std::string& x, T& result) const
3747 14 Aug 18 peter 606       {
3747 14 Aug 18 peter 607         result = utility::convert<T>(x);
3747 14 Aug 18 peter 608       }
3747 14 Aug 18 peter 609     };
3747 14 Aug 18 peter 610
3747 14 Aug 18 peter 611
3747 14 Aug 18 peter 612     template<>
3747 14 Aug 18 peter 613     struct Converter<std::string>
3747 14 Aug 18 peter 614     {
3747 14 Aug 18 peter 615       void operator()(const std::string& x, std::string& result) const
3747 14 Aug 18 peter 616       {
3747 14 Aug 18 peter 617         result = x;
3747 14 Aug 18 peter 618       }
3747 14 Aug 18 peter 619     };
3747 14 Aug 18 peter 620
3747 14 Aug 18 peter 621
3747 14 Aug 18 peter 622     template<typename T>
3747 14 Aug 18 peter 623     struct Converter<std::vector<T> >
3747 14 Aug 18 peter 624     {
3747 14 Aug 18 peter 625       void operator()(const std::string& x, std::vector<T>& result) const
3747 14 Aug 18 peter 626       {
3747 14 Aug 18 peter 627         std::vector<std::string> vec;
3747 14 Aug 18 peter 628         utility::split(vec, x, ',');
3747 14 Aug 18 peter 629         result.clear();
3747 14 Aug 18 peter 630         result.reserve(vec.size());
3747 14 Aug 18 peter 631         for (size_t i=0; i<vec.size(); ++i)
3747 14 Aug 18 peter 632           result.push_back(utility::convert<T>(vec[i]));
3747 14 Aug 18 peter 633       }
3747 14 Aug 18 peter 634     };
3747 14 Aug 18 peter 635
3747 14 Aug 18 peter 636
3747 14 Aug 18 peter 637     template<>
3747 14 Aug 18 peter 638     struct Converter<std::vector<std::string> >
3747 14 Aug 18 peter 639     {
3747 14 Aug 18 peter 640       void
3747 14 Aug 18 peter 641       operator()(const std::string& x, std::vector<std::string>& result) const
3747 14 Aug 18 peter 642       {
3747 14 Aug 18 peter 643         result.clear();
3747 14 Aug 18 peter 644         utility::split(result, x, ',');
3747 14 Aug 18 peter 645       }
3747 14 Aug 18 peter 646     };
3747 14 Aug 18 peter 647
3747 14 Aug 18 peter 648   }
3747 14 Aug 18 peter 649
3747 14 Aug 18 peter 650   /// \endcond
3747 14 Aug 18 peter 651
3747 14 Aug 18 peter 652
3747 14 Aug 18 peter 653   // template implementations
3747 14 Aug 18 peter 654
3747 14 Aug 18 peter 655   // VCF::Info
3747 14 Aug 18 peter 656
3747 14 Aug 18 peter 657   template<typename T>
3747 14 Aug 18 peter 658   void VCF::Info::add(const std::string& key, const T& value)
3747 14 Aug 18 peter 659   {
4145 24 Feb 22 peter 660     index();
4145 24 Feb 22 peter 661     str_ = "";
4145 24 Feb 22 peter 662     keys_.push_back(key);
3747 14 Aug 18 peter 663     detail::Appender<T> append(',');
4145 24 Feb 22 peter 664     append(map_[key], value);
3747 14 Aug 18 peter 665   }
3747 14 Aug 18 peter 666
3747 14 Aug 18 peter 667
3790 04 Apr 19 peter 668   template<typename T>
3790 04 Apr 19 peter 669   void VCF::Info::get(T& result, const std::string& key) const
3790 04 Apr 19 peter 670   {
4145 24 Feb 22 peter 671     index();
3790 04 Apr 19 peter 672     detail::Converter<T> converter;
4145 24 Feb 22 peter 673     auto it = map_.find(key);
3790 04 Apr 19 peter 674
4021 26 Nov 20 peter 675     try {
4145 24 Feb 22 peter 676       if (it == map_.end())
4145 24 Feb 22 peter 677         converter("", result);
4145 24 Feb 22 peter 678       else
4145 24 Feb 22 peter 679         converter(it->second, result);
3790 04 Apr 19 peter 680     }
4021 26 Nov 20 peter 681     catch (utility::runtime_error& e) {
4021 26 Nov 20 peter 682       std::ostringstream msg;
4021 26 Nov 20 peter 683       msg << "invalid VCF::Info format for key '" << key << "'";
4021 26 Nov 20 peter 684       std::throw_with_nested(utility::runtime_error(msg.str()));
4021 26 Nov 20 peter 685     }
3790 04 Apr 19 peter 686   }
3790 04 Apr 19 peter 687
3790 04 Apr 19 peter 688
4145 24 Feb 22 peter 689   template<typename T>
4145 24 Feb 22 peter 690   void VCF::Info::set(const std::string& key, const T& value)
4145 24 Feb 22 peter 691   {
4145 24 Feb 22 peter 692     index();
4145 24 Feb 22 peter 693     str_ = "";
4145 24 Feb 22 peter 694     auto it = map_.lower_bound(key);
4316 17 Feb 23 peter 695     if (it==map_.end() || it->first!=key) {
4145 24 Feb 22 peter 696       keys_.push_back(key);
4145 24 Feb 22 peter 697       it = map_.insert(it, std::make_pair(key, ""));
4145 24 Feb 22 peter 698     }
4316 17 Feb 23 peter 699     else
4316 17 Feb 23 peter 700       it->second.clear();
4145 24 Feb 22 peter 701     detail::Appender<T> append(',');
4145 24 Feb 22 peter 702     append(it->second, value);
4145 24 Feb 22 peter 703   }
4145 24 Feb 22 peter 704
3747 14 Aug 18 peter 705   // VCF::Data
3747 14 Aug 18 peter 706
3747 14 Aug 18 peter 707   template<typename T>
3747 14 Aug 18 peter 708   void VCF::Data::add(const std::string& key, const std::vector<T>& input)
3747 14 Aug 18 peter 709   {
3747 14 Aug 18 peter 710     YAT_ASSERT(key != "");
3747 14 Aug 18 peter 711     format_.push_back(key);
3747 14 Aug 18 peter 712     detail::Appender<T> append(',');
3747 14 Aug 18 peter 713     for (size_t i=0; i<input.size(); ++i) {
3747 14 Aug 18 peter 714       data_[i].push_back("");
3747 14 Aug 18 peter 715       append(data_[i].back(), input[i]);
3747 14 Aug 18 peter 716     }
3747 14 Aug 18 peter 717   }
3747 14 Aug 18 peter 718
3747 14 Aug 18 peter 719
3747 14 Aug 18 peter 720   template<typename T>
3747 14 Aug 18 peter 721   void VCF::Data::get(const std::string& key, std::vector<T>& data) const
3747 14 Aug 18 peter 722   {
3747 14 Aug 18 peter 723     size_t idx = index(key);
3747 14 Aug 18 peter 724     detail::Converter<T> converter;
3747 14 Aug 18 peter 725     data.resize(data_.size());
4021 26 Nov 20 peter 726     try {
4021 26 Nov 20 peter 727       for (size_t i=0; i<data.size(); ++i)
4021 26 Nov 20 peter 728         converter(data_[i][idx], data[i]);
4021 26 Nov 20 peter 729     }
4021 26 Nov 20 peter 730     catch (utility::runtime_error& e) {
4021 26 Nov 20 peter 731       std::ostringstream msg;
4021 26 Nov 20 peter 732       msg << "invalid VCF::Data format for key '" << key << "'";
4021 26 Nov 20 peter 733       std::throw_with_nested(utility::runtime_error(msg.str()));
4021 26 Nov 20 peter 734     }
3747 14 Aug 18 peter 735   }
3747 14 Aug 18 peter 736
3747 14 Aug 18 peter 737
4170 09 May 22 peter 738   template<typename T>
4170 09 May 22 peter 739   void VCF::Data::set(const std::string& key, const std::vector<T>& input)
4170 09 May 22 peter 740   {
4170 09 May 22 peter 741     YAT_ASSERT(key != "");
4170 09 May 22 peter 742     auto it = std::find(format_.begin(), format_.end(), key);
4170 09 May 22 peter 743     if (it == format_.end())
4170 09 May 22 peter 744       add(key, input);
4170 09 May 22 peter 745
4170 09 May 22 peter 746     size_t idx = it - format_.begin();
4170 09 May 22 peter 747     detail::Appender<T> append(',');
4170 09 May 22 peter 748     for (size_t i=0; i<input.size(); ++i) {
4170 09 May 22 peter 749       data_[i][idx] = "";
4170 09 May 22 peter 750       append(data_[i][idx], input[i]);
4170 09 May 22 peter 751     }
4170 09 May 22 peter 752   }
4170 09 May 22 peter 753
4170 09 May 22 peter 754
3747 14 Aug 18 peter 755 }}}
3747 14 Aug 18 peter 756 #endif

Generated on Sat Oct 14 02:22:54 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)