yat/omic/VcfHeader.cc

Code
Comments
Other
Rev Date Author Line
3750 17 Oct 18 peter 1 // $Id$
3750 17 Oct 18 peter 2
3750 17 Oct 18 peter 3 /*
4308 10 Feb 23 peter 4   Copyright (C) 2018, 2020, 2021, 2023 Peter Johansson
3750 17 Oct 18 peter 5
3750 17 Oct 18 peter 6   The yat library is free software; you can redistribute it and/or
3750 17 Oct 18 peter 7   modify it under the terms of the GNU General Public License as
3750 17 Oct 18 peter 8   published by the Free Software Foundation; either version 3 of the
3750 17 Oct 18 peter 9   License, or (at your option) any later version.
3750 17 Oct 18 peter 10
3750 17 Oct 18 peter 11   The yat library is distributed in the hope that it will be useful,
3750 17 Oct 18 peter 12   but WITHOUT ANY WARRANTY; without even the implied warranty of
3750 17 Oct 18 peter 13   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
3750 17 Oct 18 peter 14   General Public License for more details.
3750 17 Oct 18 peter 15
3750 17 Oct 18 peter 16   You should have received a copy of the GNU General Public License
3754 17 Oct 18 peter 17   along with yat. If not, see <http://www.gnu.org/licenses/>.
3750 17 Oct 18 peter 18 */
3750 17 Oct 18 peter 19
3750 17 Oct 18 peter 20 #include <config.h>
3750 17 Oct 18 peter 21
3750 17 Oct 18 peter 22 #include "VcfHeader.h"
3750 17 Oct 18 peter 23
3762 19 Oct 18 peter 24 #include <yat/utility/Exception.h>
3750 17 Oct 18 peter 25 #include <yat/utility/split.h>
3750 17 Oct 18 peter 26 #include <yat/utility/utility.h>
3750 17 Oct 18 peter 27
3750 17 Oct 18 peter 28 #include <algorithm>
3750 17 Oct 18 peter 29 #include <cassert>
4308 10 Feb 23 peter 30 #include <ctime>
3750 17 Oct 18 peter 31 #include <istream>
3750 17 Oct 18 peter 32 #include <iterator>
3750 17 Oct 18 peter 33 #include <ostream>
3750 17 Oct 18 peter 34 #include <sstream>
4040 15 Feb 21 peter 35 #include <utility>
3750 17 Oct 18 peter 36
3750 17 Oct 18 peter 37 namespace theplu {
3750 17 Oct 18 peter 38 namespace yat {
3750 17 Oct 18 peter 39 namespace omic {
3750 17 Oct 18 peter 40
3750 17 Oct 18 peter 41   VcfHeader::VcfHeader(const std::string& version)
3750 17 Oct 18 peter 42   {
3750 17 Oct 18 peter 43     add("fileformat", "VCFv" + version);
3750 17 Oct 18 peter 44     std::vector<std::string> tmp;
3750 17 Oct 18 peter 45     samples(tmp);
3750 17 Oct 18 peter 46   }
3750 17 Oct 18 peter 47
3750 17 Oct 18 peter 48
3750 17 Oct 18 peter 49   VcfHeader::VcfHeader(std::istream& is)
3750 17 Oct 18 peter 50   {
3762 19 Oct 18 peter 51     if (!is.good())
3762 19 Oct 18 peter 52       throw utility::IO_error("VcfHeader: istream is not good");
3762 19 Oct 18 peter 53
3750 17 Oct 18 peter 54     assert(is.good());
3750 17 Oct 18 peter 55     if (is.peek()!='#') {
3762 19 Oct 18 peter 56       std::ostringstream ss;
3762 19 Oct 18 peter 57       char buffer[16];
3762 19 Oct 18 peter 58       is.get(buffer, 16, '\n');
3762 19 Oct 18 peter 59       ss << "VcfHeader: invalid format: " << buffer;
3762 19 Oct 18 peter 60       throw utility::IO_error(ss.str());
3750 17 Oct 18 peter 61     }
3750 17 Oct 18 peter 62     std::string line;
3750 17 Oct 18 peter 63     while (getline(is, line)) {
4040 15 Feb 21 peter 64       lines_.push_back(std::move(line));
3750 17 Oct 18 peter 65       if (is.peek()!='#') // next line does not start with '#'
3750 17 Oct 18 peter 66         break;
3750 17 Oct 18 peter 67     }
3750 17 Oct 18 peter 68   }
3750 17 Oct 18 peter 69
3750 17 Oct 18 peter 70
3750 17 Oct 18 peter 71   void VcfHeader::add(const std::string& key, const std::string& value)
3750 17 Oct 18 peter 72   {
3750 17 Oct 18 peter 73     std::string line("##");
3750 17 Oct 18 peter 74     line += key;
3750 17 Oct 18 peter 75     line += "=";
3750 17 Oct 18 peter 76     line += value;
3750 17 Oct 18 peter 77     if (lines_.empty()) {
4040 15 Feb 21 peter 78       lines_.push_back(std::move(line));
3750 17 Oct 18 peter 79       return;
3750 17 Oct 18 peter 80     }
3750 17 Oct 18 peter 81     typedef std::vector<std::string>::iterator iterator;
3750 17 Oct 18 peter 82     bool found_section = false;
3750 17 Oct 18 peter 83     for (iterator i = lines_.begin(); i!=lines_.end(); ++i) {
3750 17 Oct 18 peter 84       if (i->substr(0,2)!="##") {
3750 17 Oct 18 peter 85         lines_.insert(i, line);
3750 17 Oct 18 peter 86         return;
3750 17 Oct 18 peter 87       }
3750 17 Oct 18 peter 88       if (i->size() >= key.size()+2 && i->substr(2, key.size())==key)
3750 17 Oct 18 peter 89         found_section = true;
3750 17 Oct 18 peter 90       else if (found_section) {
3750 17 Oct 18 peter 91         lines_.insert(i, line);
3750 17 Oct 18 peter 92         return;
3750 17 Oct 18 peter 93       }
3750 17 Oct 18 peter 94     }
3750 17 Oct 18 peter 95   }
3750 17 Oct 18 peter 96
3750 17 Oct 18 peter 97
3750 17 Oct 18 peter 98   void VcfHeader::add(const std::string& key, const std::string& ID, int number,
3750 17 Oct 18 peter 99                       const std::string& type, const std::string& description)
3750 17 Oct 18 peter 100   {
3750 17 Oct 18 peter 101     remove(key, ID);
3750 17 Oct 18 peter 102     std::stringstream value;
3750 17 Oct 18 peter 103     value << "<ID=" << ID << ",Number=";
3750 17 Oct 18 peter 104     if (number)
3750 17 Oct 18 peter 105       value << number;
3750 17 Oct 18 peter 106     else
3750 17 Oct 18 peter 107       value << ".";
3750 17 Oct 18 peter 108     value << ",Type=" << type << ",Description=\"" << description << "\">";
3750 17 Oct 18 peter 109     add(key, value.str());
3750 17 Oct 18 peter 110   }
3750 17 Oct 18 peter 111
3750 17 Oct 18 peter 112
3750 17 Oct 18 peter 113   void VcfHeader::add(const std::string& key,
3750 17 Oct 18 peter 114                       const std::string& ID, const std::string& number,
3750 17 Oct 18 peter 115                       const std::string& type,
3750 17 Oct 18 peter 116                       const std::string& description)
3750 17 Oct 18 peter 117   {
3750 17 Oct 18 peter 118     remove(key, ID);
3750 17 Oct 18 peter 119     std::stringstream value;
3750 17 Oct 18 peter 120     value << "<ID=" << ID << ",Number=" << number
3750 17 Oct 18 peter 121           << ",Type=" << type << ",Description=\"" << description << "\">";
3750 17 Oct 18 peter 122     add(key, value.str());
3750 17 Oct 18 peter 123   }
3750 17 Oct 18 peter 124
3750 17 Oct 18 peter 125
3750 17 Oct 18 peter 126   void VcfHeader::add_info(const std::string& ID, const std::string& number,
3750 17 Oct 18 peter 127                            const std::string& type,
3750 17 Oct 18 peter 128                            const std::string& description)
3750 17 Oct 18 peter 129   {
3750 17 Oct 18 peter 130     add("INFO", ID, number, type, description);
3750 17 Oct 18 peter 131   }
3750 17 Oct 18 peter 132
3750 17 Oct 18 peter 133
3750 17 Oct 18 peter 134   void VcfHeader::add_format(const std::string& ID, const std::string& number,
3750 17 Oct 18 peter 135                              const std::string& type,
3750 17 Oct 18 peter 136                              const std::string& description)
3750 17 Oct 18 peter 137   {
3750 17 Oct 18 peter 138     add("FORMAT", ID, number, type, description);
3750 17 Oct 18 peter 139   }
3750 17 Oct 18 peter 140
3750 17 Oct 18 peter 141
3750 17 Oct 18 peter 142   void VcfHeader::add_program_command(int argc, char* argv[])
3750 17 Oct 18 peter 143   {
3750 17 Oct 18 peter 144     std::string key = utility::basename(argv[0]) + "Command";
4308 10 Feb 23 peter 145     std::stringstream ss;
4308 10 Feb 23 peter 146
3750 17 Oct 18 peter 147     for (int i=1; i<argc; ++i) {
3750 17 Oct 18 peter 148       if (i!=1)
4308 10 Feb 23 peter 149         ss << ' ';
4308 10 Feb 23 peter 150       ss << argv[i];
3750 17 Oct 18 peter 151     }
4308 10 Feb 23 peter 152     ss << "; Date=";
4308 10 Feb 23 peter 153     std::time_t time = std::time(nullptr);
4308 10 Feb 23 peter 154     ss << std::asctime(std::localtime(&time));
4308 10 Feb 23 peter 155
4308 10 Feb 23 peter 156     std::string val;
4308 10 Feb 23 peter 157     getline(ss, val);
4308 10 Feb 23 peter 158
3750 17 Oct 18 peter 159     add(key, val);
3750 17 Oct 18 peter 160   }
3750 17 Oct 18 peter 161
3750 17 Oct 18 peter 162
3750 17 Oct 18 peter 163   void VcfHeader::add_program_version(const std::string& prog,
3750 17 Oct 18 peter 164                                       const std::string& version)
3750 17 Oct 18 peter 165   {
3750 17 Oct 18 peter 166     add(prog + "Version", version);
3750 17 Oct 18 peter 167   }
3750 17 Oct 18 peter 168
3750 17 Oct 18 peter 169
3891 26 Mar 20 peter 170   void VcfHeader::erase(size_t i)
3891 26 Mar 20 peter 171   {
3891 26 Mar 20 peter 172     assert(i+1 < lines_.size());
3891 26 Mar 20 peter 173     lines_.erase(lines_.begin()+i);
3891 26 Mar 20 peter 174   }
3891 26 Mar 20 peter 175
3891 26 Mar 20 peter 176
3750 17 Oct 18 peter 177   const std::string& VcfHeader::header(void) const
3750 17 Oct 18 peter 178   {
3750 17 Oct 18 peter 179     assert(lines_.size());
3750 17 Oct 18 peter 180     return lines_.back();
3750 17 Oct 18 peter 181   }
3750 17 Oct 18 peter 182
3750 17 Oct 18 peter 183
3750 17 Oct 18 peter 184   void VcfHeader::header(const std::string& line)
3750 17 Oct 18 peter 185   {
3750 17 Oct 18 peter 186     assert(lines_.size());
3750 17 Oct 18 peter 187     lines_.back() = line;
3750 17 Oct 18 peter 188   }
3750 17 Oct 18 peter 189
3750 17 Oct 18 peter 190
3750 17 Oct 18 peter 191   const std::string& VcfHeader::line(size_t i) const
3750 17 Oct 18 peter 192   {
3750 17 Oct 18 peter 193     assert(i<lines_.size());
3750 17 Oct 18 peter 194     return lines_[i];
3750 17 Oct 18 peter 195   }
3750 17 Oct 18 peter 196
3750 17 Oct 18 peter 197
3750 17 Oct 18 peter 198   void VcfHeader::key_value(size_t line, std::string& key,
3750 17 Oct 18 peter 199                             std::string& value) const
3750 17 Oct 18 peter 200   {
3750 17 Oct 18 peter 201     assert(line<size());
3750 17 Oct 18 peter 202     const std::string& str = lines_[line];
3750 17 Oct 18 peter 203     if (str.size()<3 || str[0]!='#' || str[1]!='#')
3750 17 Oct 18 peter 204       return;
3750 17 Oct 18 peter 205     size_t end = str.find_first_of('=', 2);
3750 17 Oct 18 peter 206     if (end==std::string::npos)
3750 17 Oct 18 peter 207       return;
3750 17 Oct 18 peter 208     key = str.substr(2, end-2);
3750 17 Oct 18 peter 209     value = str.substr(end+1);
3750 17 Oct 18 peter 210   }
3750 17 Oct 18 peter 211
3750 17 Oct 18 peter 212
3750 17 Oct 18 peter 213   std::string VcfHeader::key(size_t line) const
3750 17 Oct 18 peter 214   {
3750 17 Oct 18 peter 215     std::string key;
3750 17 Oct 18 peter 216     std::string val;
3750 17 Oct 18 peter 217     key_value(line, key, val);
3750 17 Oct 18 peter 218     return key;
3750 17 Oct 18 peter 219   }
3750 17 Oct 18 peter 220
3750 17 Oct 18 peter 221
3750 17 Oct 18 peter 222   void VcfHeader::nicify(void)
3750 17 Oct 18 peter 223   {
3750 17 Oct 18 peter 224     std::stable_sort(lines_.begin(), lines_.end(), Compare());
3750 17 Oct 18 peter 225   }
3750 17 Oct 18 peter 226
3750 17 Oct 18 peter 227
3750 17 Oct 18 peter 228   void VcfHeader::remove(const std::string& KEY, const std::string& ID)
3750 17 Oct 18 peter 229   {
3750 17 Oct 18 peter 230     for (size_t i=0; i<lines_.size(); ++i) {
3750 17 Oct 18 peter 231       if (key(i) == KEY && value(i, "ID") == ID) {
3750 17 Oct 18 peter 232         lines_.erase(lines_.begin()+i);
3750 17 Oct 18 peter 233         --i;
3750 17 Oct 18 peter 234       }
3750 17 Oct 18 peter 235     }
3750 17 Oct 18 peter 236   }
3750 17 Oct 18 peter 237
3750 17 Oct 18 peter 238
3750 17 Oct 18 peter 239   std::string VcfHeader::value(size_t line) const
3750 17 Oct 18 peter 240   {
3750 17 Oct 18 peter 241     std::string key;
3750 17 Oct 18 peter 242     std::string val;
3750 17 Oct 18 peter 243     key_value(line, key, val);
3750 17 Oct 18 peter 244     return val;
3750 17 Oct 18 peter 245   }
3750 17 Oct 18 peter 246
3750 17 Oct 18 peter 247
3750 17 Oct 18 peter 248   std::string VcfHeader::value(size_t line, const std::string& key) const
3750 17 Oct 18 peter 249   {
3750 17 Oct 18 peter 250     const std::string& str = lines_[line];
3750 17 Oct 18 peter 251     std::string::const_iterator begin = str.begin();
3750 17 Oct 18 peter 252     begin = std::find(begin, str.end(), '=');
3750 17 Oct 18 peter 253     begin = std::find(begin, str.end(), '<');
3750 17 Oct 18 peter 254     ++begin; // '<'
3750 17 Oct 18 peter 255
3750 17 Oct 18 peter 256     while (begin < str.end()) {
3750 17 Oct 18 peter 257       std::string::const_iterator middle = std::find(begin, str.end(), '=');
3750 17 Oct 18 peter 258       std::string::const_iterator end(middle+1);
3750 17 Oct 18 peter 259       bool quoted = false;
3750 17 Oct 18 peter 260       while (end < str.end()) {
3750 17 Oct 18 peter 261         if (*end == '"')
3750 17 Oct 18 peter 262           quoted = !quoted;
3750 17 Oct 18 peter 263         else if (!quoted && *end == ',')
3750 17 Oct 18 peter 264           break;
3750 17 Oct 18 peter 265         ++end;
3750 17 Oct 18 peter 266       }
3750 17 Oct 18 peter 267       if (*(end-1) == '>')
3750 17 Oct 18 peter 268         --end;
3750 17 Oct 18 peter 269       if (static_cast<size_t>(middle - begin) == key.size() &&
3750 17 Oct 18 peter 270           std::equal(begin,middle,key.begin()))
3750 17 Oct 18 peter 271         return std::string(middle+1,end);
3750 17 Oct 18 peter 272       begin = end + 1;
3750 17 Oct 18 peter 273     }
3750 17 Oct 18 peter 274     return "";
3750 17 Oct 18 peter 275   }
3750 17 Oct 18 peter 276
3750 17 Oct 18 peter 277
3750 17 Oct 18 peter 278   const std::vector<std::string>& VcfHeader::samples(void) const
3750 17 Oct 18 peter 279   {
3750 17 Oct 18 peter 280     if (samples_.empty()) {
3750 17 Oct 18 peter 281       utility::split(samples_, header(), '\t');
3750 17 Oct 18 peter 282       if (samples_.size()>9)
3750 17 Oct 18 peter 283         samples_.erase(samples_.begin(), samples_.begin()+9);
3750 17 Oct 18 peter 284       else
3750 17 Oct 18 peter 285         samples_.clear();
3750 17 Oct 18 peter 286     }
3750 17 Oct 18 peter 287     return samples_;
3750 17 Oct 18 peter 288   }
3750 17 Oct 18 peter 289
3750 17 Oct 18 peter 290
3750 17 Oct 18 peter 291   void VcfHeader::samples(const std::vector<std::string>& s)
3750 17 Oct 18 peter 292   {
3750 17 Oct 18 peter 293     samples_ = s;
3750 17 Oct 18 peter 294     std::string hdr("#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT");
3750 17 Oct 18 peter 295     if (lines_.size() &&  lines_.back().size()>1 && lines_.back()[1]!='#')
4040 15 Feb 21 peter 296       lines_.back() = std::move(hdr);
3750 17 Oct 18 peter 297     else
4040 15 Feb 21 peter 298       lines_.push_back(std::move(hdr));
3750 17 Oct 18 peter 299
3750 17 Oct 18 peter 300     for (size_t i=0; i<s.size(); ++i) {
3750 17 Oct 18 peter 301       lines_.back() += "\t";
3750 17 Oct 18 peter 302       lines_.back() += s[i];
3750 17 Oct 18 peter 303     }
3750 17 Oct 18 peter 304   }
3750 17 Oct 18 peter 305
3750 17 Oct 18 peter 306
3750 17 Oct 18 peter 307   size_t VcfHeader::size(void) const
3750 17 Oct 18 peter 308   {
3750 17 Oct 18 peter 309     return lines_.size();
3750 17 Oct 18 peter 310   }
3750 17 Oct 18 peter 311
3750 17 Oct 18 peter 312
3750 17 Oct 18 peter 313   std::ostream& operator<<(std::ostream& os, const VcfHeader& header)
3750 17 Oct 18 peter 314   {
3750 17 Oct 18 peter 315     std::copy(header.lines_.begin(), header.lines_.end(),
3750 17 Oct 18 peter 316               std::ostream_iterator<std::string>(os, "\n"));
3750 17 Oct 18 peter 317     return os;
3750 17 Oct 18 peter 318   }
3750 17 Oct 18 peter 319
3760 18 Oct 18 peter 320
3760 18 Oct 18 peter 321   bool VcfHeader::Compare::operator()(const std::string& lhs,
3760 18 Oct 18 peter 322                                       const std::string& rhs) const
3760 18 Oct 18 peter 323   {
3760 18 Oct 18 peter 324     return index(lhs) < index(rhs);
3760 18 Oct 18 peter 325   }
3760 18 Oct 18 peter 326
3760 18 Oct 18 peter 327
3760 18 Oct 18 peter 328   int VcfHeader::Compare::index(const std::string& s) const
3760 18 Oct 18 peter 329   {
3760 18 Oct 18 peter 330     if (match(s, "##fileformat="))
3760 18 Oct 18 peter 331       return 0;
3760 18 Oct 18 peter 332     if (match(s, "#CHROM"))
3760 18 Oct 18 peter 333       return 99;
3760 18 Oct 18 peter 334
3760 18 Oct 18 peter 335     if (match(s, "##reference="))
3760 18 Oct 18 peter 336       return 4;
3760 18 Oct 18 peter 337     if (match(s, "##contig="))
3760 18 Oct 18 peter 338       return 4;
3760 18 Oct 18 peter 339     if (match(s, "##INFO="))
3760 18 Oct 18 peter 340       return 5;
3760 18 Oct 18 peter 341     if (match(s, "##FILTER="))
3760 18 Oct 18 peter 342       return 6;
3760 18 Oct 18 peter 343     if (match(s, "##FORMAT="))
3760 18 Oct 18 peter 344       return 7;
3760 18 Oct 18 peter 345
3760 18 Oct 18 peter 346     return 9;
3760 18 Oct 18 peter 347   }
3760 18 Oct 18 peter 348
3760 18 Oct 18 peter 349
3760 18 Oct 18 peter 350   bool VcfHeader::Compare::match(const std::string& s,
3760 18 Oct 18 peter 351                                  const std::string& prefix) const
3760 18 Oct 18 peter 352   {
4041 15 Feb 21 peter 353     return s.size() >= prefix.size() &&
4041 15 Feb 21 peter 354       std::equal(prefix.begin(), prefix.end(), s.begin());
3760 18 Oct 18 peter 355   }
3760 18 Oct 18 peter 356
3750 17 Oct 18 peter 357 }}}

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