yat/omic/BamHeader.cc

Code
Comments
Other
Rev Date Author Line
2883 03 Dec 12 peter 1 // $Id$
2883 03 Dec 12 peter 2
2993 03 Mar 13 peter 3 /*
4359 23 Aug 23 peter 4   Copyright (C) 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2023 Peter Johansson
2993 03 Mar 13 peter 5
2993 03 Mar 13 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
2993 03 Mar 13 peter 7
2993 03 Mar 13 peter 8   The yat library is free software; you can redistribute it and/or
2993 03 Mar 13 peter 9   modify it under the terms of the GNU General Public License as
2993 03 Mar 13 peter 10   published by the Free Software Foundation; either version 3 of the
2993 03 Mar 13 peter 11   License, or (at your option) any later version.
2993 03 Mar 13 peter 12
2993 03 Mar 13 peter 13   The yat library is distributed in the hope that it will be useful,
2993 03 Mar 13 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
2993 03 Mar 13 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2993 03 Mar 13 peter 16   General Public License for more details.
2993 03 Mar 13 peter 17
2993 03 Mar 13 peter 18   You should have received a copy of the GNU General Public License
3752 17 Oct 18 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
2993 03 Mar 13 peter 20 */
2993 03 Mar 13 peter 21
2883 03 Dec 12 peter 22 #include <config.h>
2883 03 Dec 12 peter 23
2883 03 Dec 12 peter 24 #include "BamHeader.h"
2883 03 Dec 12 peter 25
2999 14 Mar 13 peter 26 #include "yat/utility/Exception.h"
3410 21 Apr 15 peter 27 #include "yat/utility/split.h"
3412 22 Apr 15 peter 28 #include "yat/utility/stl_utility.h"
2999 14 Mar 13 peter 29
2883 03 Dec 12 peter 30 #include <cassert>
3408 16 Apr 15 peter 31 #include <cstdlib>
3408 16 Apr 15 peter 32 #include <cstring>
2999 14 Mar 13 peter 33 #include <string>
2999 14 Mar 13 peter 34 #include <sstream>
3408 16 Apr 15 peter 35 #include <utility>
2883 03 Dec 12 peter 36
2883 03 Dec 12 peter 37 namespace theplu {
2883 03 Dec 12 peter 38 namespace yat {
2883 03 Dec 12 peter 39 namespace omic {
2883 03 Dec 12 peter 40
2883 03 Dec 12 peter 41   BamHeader::BamHeader(void)
2883 03 Dec 12 peter 42     : header_(NULL)
2883 03 Dec 12 peter 43   {
2883 03 Dec 12 peter 44   }
2883 03 Dec 12 peter 45
2883 03 Dec 12 peter 46
3408 16 Apr 15 peter 47   BamHeader::~BamHeader(void)
3408 16 Apr 15 peter 48   {
3408 16 Apr 15 peter 49     bam_hdr_destroy(header_);
3408 16 Apr 15 peter 50   }
3408 16 Apr 15 peter 51
3408 16 Apr 15 peter 52
3408 16 Apr 15 peter 53   BamHeader::BamHeader(const BamHeader& other)
3410 21 Apr 15 peter 54     : read_group_(other.read_group_), program_group_(other.program_group_)
3408 16 Apr 15 peter 55   {
3408 16 Apr 15 peter 56     header_ = bam_hdr_dup(other.header_);
4265 13 Jan 23 peter 57     if (other.header_ && !header_)
4265 13 Jan 23 peter 58       throw utility::runtime_error("BamHeader: bam_hdr_dup failed");
3408 16 Apr 15 peter 59   }
3408 16 Apr 15 peter 60
3408 16 Apr 15 peter 61
3691 14 Sep 17 peter 62   BamHeader::BamHeader(BamHeader&& other) noexcept
3600 22 Jan 17 peter 63     : header_(NULL)
3600 22 Jan 17 peter 64   {
3600 22 Jan 17 peter 65     swap(other);
3600 22 Jan 17 peter 66   }
3600 22 Jan 17 peter 67
3600 22 Jan 17 peter 68
4290 03 Feb 23 peter 69   const std::map<std::string, std::string>&
4290 03 Feb 23 peter 70   BamHeader::group(strMap2& map, const std::string& type,
4290 03 Feb 23 peter 71                    const std::string& id) const
4290 03 Feb 23 peter 72   {
4290 03 Feb 23 peter 73     update_group(map, type);
4290 03 Feb 23 peter 74     strMap2::const_iterator line_table = map.find(id);
4290 03 Feb 23 peter 75     if (line_table == map.end()) {
4290 03 Feb 23 peter 76       std::string msg = "@" + type + ": unknown ID: '" + id + "'";
4290 03 Feb 23 peter 77       throw utility::runtime_error(msg);
4290 03 Feb 23 peter 78     }
4290 03 Feb 23 peter 79     return line_table->second;
4290 03 Feb 23 peter 80   }
4290 03 Feb 23 peter 81
4290 03 Feb 23 peter 82
4290 03 Feb 23 peter 83   BamHeader::read_group_iterator
4290 03 Feb 23 peter 84   BamHeader::group_begin(strMap2& group, const std::string& type) const
4290 03 Feb 23 peter 85   {
4290 03 Feb 23 peter 86     update_group(group, type);
4290 03 Feb 23 peter 87     BamHeader::strMap2::const_iterator it(group.begin());
4290 03 Feb 23 peter 88     return utility::pair_first_iterator(it);
4290 03 Feb 23 peter 89   }
4290 03 Feb 23 peter 90
4290 03 Feb 23 peter 91
4290 03 Feb 23 peter 92   BamHeader::read_group_iterator
4290 03 Feb 23 peter 93   BamHeader::group_end(strMap2& group, const std::string& type) const
4290 03 Feb 23 peter 94   {
4290 03 Feb 23 peter 95     update_group(group, type);
4290 03 Feb 23 peter 96     BamHeader::strMap2::const_iterator it(group.end());
4290 03 Feb 23 peter 97     return utility::pair_first_iterator(it);
4290 03 Feb 23 peter 98   }
4290 03 Feb 23 peter 99
4290 03 Feb 23 peter 100
3410 21 Apr 15 peter 101   const std::string& BamHeader::group(strMap2& table,const std::string& type,
3410 21 Apr 15 peter 102                                       const std::string& id,
3410 21 Apr 15 peter 103                                       const std::string& key) const
3410 21 Apr 15 peter 104   {
4290 03 Feb 23 peter 105     const std::map<std::string, std::string>& m = group(table, type, id);
4290 03 Feb 23 peter 106     strMap::const_iterator it = m.find(key);
4290 03 Feb 23 peter 107     if (it == m.end()) {
3477 09 Mar 16 peter 108       std::string msg = "@" + type + ": unknown KEY: '" + key + "'";
3477 09 Mar 16 peter 109       throw utility::runtime_error(msg);
3477 09 Mar 16 peter 110     }
3477 09 Mar 16 peter 111     return it->second;
3477 09 Mar 16 peter 112   }
3477 09 Mar 16 peter 113
3477 09 Mar 16 peter 114
3477 09 Mar 16 peter 115   void BamHeader::update_group(strMap2& table, const std::string& type) const
3477 09 Mar 16 peter 116   {
3410 21 Apr 15 peter 117     assert(type.size()==2);
3410 21 Apr 15 peter 118     if (table.empty()) { // parse header
3410 21 Apr 15 peter 119       std::istringstream ss(text());
3410 21 Apr 15 peter 120       std::string line;
3410 21 Apr 15 peter 121       while (getline(ss, line)) {
3410 21 Apr 15 peter 122         if (line.substr(1,2)!=type)
3410 21 Apr 15 peter 123           continue;
3410 21 Apr 15 peter 124         // find ID value
3410 21 Apr 15 peter 125         size_t pos = line.find("\tID:", 3);
3410 21 Apr 15 peter 126         size_t end = line.find('\t', pos+4);
3410 21 Apr 15 peter 127         std::string id_val = line.substr(pos+4, end-(pos+4));
3410 21 Apr 15 peter 128         std::vector<std::string> vec;
3410 21 Apr 15 peter 129         utility::split(vec, line, '\t');
3410 21 Apr 15 peter 130         strMap& m = table[id_val];
3410 21 Apr 15 peter 131         for (size_t i=1; i<vec.size(); ++i) {
3410 21 Apr 15 peter 132           m[vec[i].substr(0,2)] = vec[i].substr(3);
3410 21 Apr 15 peter 133         }
3410 21 Apr 15 peter 134       }
3410 21 Apr 15 peter 135     }
3410 21 Apr 15 peter 136   }
3410 21 Apr 15 peter 137
3410 21 Apr 15 peter 138
3358 23 Nov 14 peter 139   void BamHeader::parse_region(const std::string& reg, int& id, int& begin,
2999 14 Mar 13 peter 140                                int& end) const
2999 14 Mar 13 peter 141   {
2999 14 Mar 13 peter 142     assert(header_);
3358 23 Nov 14 peter 143     const char* b = reg.c_str();
3358 23 Nov 14 peter 144     const char* e = hts_parse_reg(b, &begin, &end);
3358 23 Nov 14 peter 145     // If begin > end suggests something is wrong and that is also how
3358 23 Nov 14 peter 146     // bam_parse_region (below in samtools impl) used to detect error
3358 23 Nov 14 peter 147     // and return non-zero.
3358 23 Nov 14 peter 148     if (begin<=end) {
3358 23 Nov 14 peter 149       std::string chr = reg.substr(0, e-b);
3358 23 Nov 14 peter 150       try {
3358 23 Nov 14 peter 151         id = tid(chr);
3358 23 Nov 14 peter 152         return;
3358 23 Nov 14 peter 153       }
3358 23 Nov 14 peter 154       catch (utility::runtime_error& e) {
3358 23 Nov 14 peter 155         ; // throw below instead
3358 23 Nov 14 peter 156       }
3358 23 Nov 14 peter 157     }
2999 14 Mar 13 peter 158     std::ostringstream ss;
2999 14 Mar 13 peter 159     ss << "invalid region: '" << reg << "'";
2999 14 Mar 13 peter 160     throw utility::runtime_error(ss.str());
2999 14 Mar 13 peter 161   }
2999 14 Mar 13 peter 162
2999 14 Mar 13 peter 163
3410 21 Apr 15 peter 164   const std::string& BamHeader::program_group(const std::string& id,
3410 21 Apr 15 peter 165                                               const std::string& key) const
3410 21 Apr 15 peter 166   {
3410 21 Apr 15 peter 167     return group(program_group_, "PG", id, key);
3410 21 Apr 15 peter 168   }
3410 21 Apr 15 peter 169
3410 21 Apr 15 peter 170
4290 03 Feb 23 peter 171   const std::map<std::string, std::string>&
4290 03 Feb 23 peter 172   BamHeader::program_group(const std::string& id) const
4290 03 Feb 23 peter 173   {
4290 03 Feb 23 peter 174     return group(program_group_, "PG", id);
4290 03 Feb 23 peter 175   }
4290 03 Feb 23 peter 176
4290 03 Feb 23 peter 177
4290 03 Feb 23 peter 178   BamHeader::program_group_iterator BamHeader::program_group_begin(void) const
4290 03 Feb 23 peter 179   {
4290 03 Feb 23 peter 180     return group_begin(program_group_, "PG");
4290 03 Feb 23 peter 181   }
4290 03 Feb 23 peter 182
4290 03 Feb 23 peter 183
4290 03 Feb 23 peter 184   BamHeader::program_group_iterator BamHeader::program_group_end(void) const
4290 03 Feb 23 peter 185   {
4290 03 Feb 23 peter 186     return group_end(program_group_, "PG");
4290 03 Feb 23 peter 187   }
4290 03 Feb 23 peter 188
4290 03 Feb 23 peter 189
3410 21 Apr 15 peter 190   const std::string& BamHeader::read_group(const std::string& id,
3410 21 Apr 15 peter 191                                            const std::string& key) const
3410 21 Apr 15 peter 192   {
3410 21 Apr 15 peter 193     return group(read_group_, "RG", id, key);
3410 21 Apr 15 peter 194   }
3410 21 Apr 15 peter 195
3410 21 Apr 15 peter 196
3477 09 Mar 16 peter 197   const std::map<std::string, std::string>&
3477 09 Mar 16 peter 198   BamHeader::read_group(const std::string& id) const
3477 09 Mar 16 peter 199   {
4290 03 Feb 23 peter 200     return group(read_group_, "RG", id);
3477 09 Mar 16 peter 201   }
3477 09 Mar 16 peter 202
3477 09 Mar 16 peter 203
3477 09 Mar 16 peter 204   BamHeader::read_group_iterator BamHeader::read_group_begin(void) const
3477 09 Mar 16 peter 205   {
4290 03 Feb 23 peter 206     return group_begin(read_group_, "RG");
3477 09 Mar 16 peter 207   }
3477 09 Mar 16 peter 208
3477 09 Mar 16 peter 209
3477 09 Mar 16 peter 210   BamHeader::read_group_iterator BamHeader::read_group_end(void) const
3477 09 Mar 16 peter 211   {
4290 03 Feb 23 peter 212     return group_end(read_group_, "RG");
3477 09 Mar 16 peter 213   }
3477 09 Mar 16 peter 214
3477 09 Mar 16 peter 215
3408 16 Apr 15 peter 216   void BamHeader::swap(BamHeader& other)
3408 16 Apr 15 peter 217   {
3408 16 Apr 15 peter 218     std::swap(header_, other.header_);
3410 21 Apr 15 peter 219     std::swap(read_group_, other.read_group_);
3410 21 Apr 15 peter 220     std::swap(program_group_, other.program_group_);
3408 16 Apr 15 peter 221   }
3408 16 Apr 15 peter 222
3408 16 Apr 15 peter 223
2960 17 Jan 13 peter 224   uint32_t BamHeader::target_length(size_t tid) const
2960 17 Jan 13 peter 225   {
2960 17 Jan 13 peter 226     assert(tid < static_cast<size_t>(n_targets()));
2960 17 Jan 13 peter 227     return header_->target_len[tid];
2960 17 Jan 13 peter 228   }
2960 17 Jan 13 peter 229
2960 17 Jan 13 peter 230
2883 03 Dec 12 peter 231   const char* BamHeader::target_name(size_t tid) const
2883 03 Dec 12 peter 232   {
2883 03 Dec 12 peter 233     assert(tid < static_cast<size_t>(n_targets()));
2883 03 Dec 12 peter 234     return header_->target_name[tid];
2883 03 Dec 12 peter 235   }
2883 03 Dec 12 peter 236
2883 03 Dec 12 peter 237
3408 16 Apr 15 peter 238   std::string BamHeader::text(void) const
3408 16 Apr 15 peter 239   {
3408 16 Apr 15 peter 240     assert(header_);
3408 16 Apr 15 peter 241     return std::string(header_->text, header_->l_text);
3408 16 Apr 15 peter 242   }
3408 16 Apr 15 peter 243
3409 17 Apr 15 peter 244
3408 16 Apr 15 peter 245   void BamHeader::text(const std::string& txt)
3408 16 Apr 15 peter 246   {
3573 12 Jan 17 peter 247     BamHeader tmp;
3573 12 Jan 17 peter 248     tmp.header_ = sam_hdr_parse(txt.size(), txt.c_str());
3573 12 Jan 17 peter 249
3573 12 Jan 17 peter 250     // throw if failed
3573 12 Jan 17 peter 251     if (!tmp.header_) {
3573 12 Jan 17 peter 252       std::ostringstream msg;
3573 12 Jan 17 peter 253       msg << "BamHeader::text: " << txt;
3573 12 Jan 17 peter 254       throw utility::runtime_error(msg.str());
3573 12 Jan 17 peter 255     }
3573 12 Jan 17 peter 256     tmp.header_->l_text = txt.size();
3573 12 Jan 17 peter 257     tmp.header_->text = (char*)realloc(tmp.header_->text, txt.size()+1);
3573 12 Jan 17 peter 258     if (!tmp.header_->text) {
3573 12 Jan 17 peter 259       std::ostringstream msg;
3573 12 Jan 17 peter 260       msg << "BamHeader::text(1): allocation failed\n" << txt;
3573 12 Jan 17 peter 261       throw utility::runtime_error(msg.str());
3573 12 Jan 17 peter 262     }
3573 12 Jan 17 peter 263     // copy to temporary
3573 12 Jan 17 peter 264     memcpy(tmp.header_->text, txt.c_str(), txt.size()+1);
3573 12 Jan 17 peter 265
3573 12 Jan 17 peter 266     // everything worked; swap data so this get the new stuff and
3573 12 Jan 17 peter 267     // destructor of tmp takes care of the old stuff.
3573 12 Jan 17 peter 268     swap(tmp);
3408 16 Apr 15 peter 269   }
3408 16 Apr 15 peter 270
3408 16 Apr 15 peter 271
2999 14 Mar 13 peter 272   int32_t BamHeader::tid(const std::string& name) const
2999 14 Mar 13 peter 273   {
3358 23 Nov 14 peter 274     int res = bam_name2id(header_, name.c_str());
3358 23 Nov 14 peter 275     if (res==-1) {
3358 23 Nov 14 peter 276       std::string msg("invalid segment name:");
3358 23 Nov 14 peter 277       msg += name;
3358 23 Nov 14 peter 278       throw utility::runtime_error(msg);
3358 23 Nov 14 peter 279     }
3358 23 Nov 14 peter 280     return res;
2999 14 Mar 13 peter 281   }
2999 14 Mar 13 peter 282
2999 14 Mar 13 peter 283
2883 03 Dec 12 peter 284   int32_t BamHeader::n_targets(void) const
2883 03 Dec 12 peter 285   {
2883 03 Dec 12 peter 286     assert(header_);
2883 03 Dec 12 peter 287     return header_->n_targets;
2883 03 Dec 12 peter 288   }
2883 03 Dec 12 peter 289
3408 16 Apr 15 peter 290
3408 16 Apr 15 peter 291   void swap(BamHeader& lhs, BamHeader& rhs)
3408 16 Apr 15 peter 292   {
3408 16 Apr 15 peter 293     lhs.swap(rhs);
3408 16 Apr 15 peter 294   }
3408 16 Apr 15 peter 295
3408 16 Apr 15 peter 296
3408 16 Apr 15 peter 297   BamHeader& BamHeader::operator=(const BamHeader& rhs)
3408 16 Apr 15 peter 298   {
3408 16 Apr 15 peter 299     BamHeader tmp(rhs);
3408 16 Apr 15 peter 300     swap(tmp);
3408 16 Apr 15 peter 301     return *this;
3408 16 Apr 15 peter 302   }
3408 16 Apr 15 peter 303
3600 22 Jan 17 peter 304
3600 22 Jan 17 peter 305   BamHeader& BamHeader::operator=(BamHeader&& rhs)
3600 22 Jan 17 peter 306   {
3600 22 Jan 17 peter 307     swap(rhs);
3600 22 Jan 17 peter 308     return *this;
3600 22 Jan 17 peter 309   }
3600 22 Jan 17 peter 310
2883 03 Dec 12 peter 311 }}}

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