Changeset 3410


Ignore:
Timestamp:
Apr 22, 2015, 1:05:43 AM (8 years ago)
Author:
Peter
Message:

add two new functions read_group and program_group to access data in BamHeader?. closes #817

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/bam_header.cc

    r3409 r3410  
    125125  }
    126126
     127  suite.err() << hdr.text() << "\n";
     128  try {
     129    std::string sample_name = hdr.read_group("foo", "SM");
     130    if (sample_name!="Tumour") {
     131      suite.add(false);
     132      suite.err() << "error: incorrect sample name: '" << sample_name << "'\n";
     133    }
     134
     135    std::string bwa_version = hdr.program_group("bwa", "VN");
     136    if (bwa_version!="0.6.1-r104") {
     137      suite.add(false);
     138      suite.err() << "error: incorrect bwa_version: '" << bwa_version << "'\n";
     139    }
     140  }
     141  catch (std::runtime_error& e) {
     142    suite.add(false);
     143    suite.err() << "error: exception: " << e.what() << "\n";
     144  }
    127145#endif
    128146}
  • trunk/test/data/foo.sam

    r2888 r3410  
    8383@SQ SN:GL000225.1 LN:211173
    8484@SQ SN:GL000192.1 LN:547496
     85@RG ID:foo  PL:ILLUMINA SM:Tumour
    8586@PG ID:bwa  PN:bwa  VN:0.6.1-r104
    86871403:1816-0 83  2 24819 60  100M  = 24506 -413  GCCTCTCATAGAATGCCAGGCGCCCCTGTGTTTGTCTCAGTGTCTCCAGCTTGAATTGCTTCACCCCGGCGAAGCAATGCTGCTCCTCTCTCCTGCATGG  DDD_DDbfi_hDehfihDDDfcgDeeeihD`cdegbige_Wiicihihdegiddhhiheaiia]ifh_hihicOhTidgIiiffhiff`aghiiiihihi  XT:A:U  NM:i:0  SM:i:37 AM:i:37 X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:100
  • trunk/yat/omic/BamHeader.cc

    r3409 r3410  
    2525
    2626#include "yat/utility/Exception.h"
     27#include "yat/utility/split.h"
    2728
    2829// we need to include 'sam_header.h' when compiling against libbam
     
    6768
    6869  BamHeader::BamHeader(const BamHeader& other)
     70    : read_group_(other.read_group_), program_group_(other.program_group_)
    6971  {
    7072#ifndef YAT_HAVE_HTSLIB
     
    7274#endif
    7375    header_ = bam_hdr_dup(other.header_);
     76  }
     77
     78
     79  const std::string& BamHeader::group(strMap2& table,const std::string& type,
     80                                      const std::string& id,
     81                                      const std::string& key) const
     82  {
     83    assert(type.size()==2);
     84    if (table.empty()) { // parse header
     85      std::istringstream ss(text());
     86      std::string line;
     87      while (getline(ss, line)) {
     88        if (line.substr(1,2)!=type)
     89          continue;
     90        // find ID value
     91        size_t pos = line.find("\tID:", 3);
     92        size_t end = line.find('\t', pos+4);
     93        std::string id_val = line.substr(pos+4, end-(pos+4));
     94        std::vector<std::string> vec;
     95        utility::split(vec, line, '\t');
     96        strMap& m = table[id_val];
     97        for (size_t i=1; i<vec.size(); ++i) {
     98          m[vec[i].substr(0,2)] = vec[i].substr(3);
     99        }
     100      }
     101      if (table.empty())
     102        table[""][""]; // to avoid re-parsing
     103    }
     104    return table[id][key];
    74105  }
    75106
     
    105136
    106137
     138  const std::string& BamHeader::program_group(const std::string& id,
     139                                              const std::string& key) const
     140  {
     141    return group(program_group_, "PG", id, key);
     142  }
     143
     144
     145  const std::string& BamHeader::read_group(const std::string& id,
     146                                           const std::string& key) const
     147  {
     148    return group(read_group_, "RG", id, key);
     149  }
     150
     151
    107152  void BamHeader::swap(BamHeader& other)
    108153  {
    109154    std::swap(header_, other.header_);
     155    std::swap(read_group_, other.read_group_);
     156    std::swap(program_group_, other.program_group_);
    110157  }
    111158
     
    134181  void BamHeader::text(const std::string& txt)
    135182  {
     183    read_group_.clear();
     184    program_group_.clear();
    136185#ifdef YAT_HAVE_HTSLIB
    137186    // destroy the old
  • trunk/yat/omic/BamHeader.h

    r3408 r3410  
    2929#include YAT_SAM_HEADER
    3030
     31#include <map>
    3132#include <string>
    3233
     
    8889
    8990    /**
     91       \brief Access value in \c \@PG lines.
     92
     93       A program group line in the header typically looks like
     94
     95       \code @PG  ID:bwa  PN:bwa  VN:0.6.1-r104 \endcode
     96
     97       and for this line \c program_group("bwa", "VN") returns
     98       \c "0.6.1-r104"
     99
     100       \return value for \a key for line with ID \a id.
     101
     102       \since New in yat 0.13
     103    */
     104    const std::string& program_group(const std::string& id,
     105                                     const std::string& key) const;
     106
     107    /**
     108       \brief Access value in \c \@RG lines.
     109
     110       A read group line in the header typically looks like
     111
     112       \code @RG  ID:foo  PL:ILLUMINA SM:Tumour \endcode
     113
     114       and for this line \c read_group("foo", "SM") returns \c "Tumour"
     115
     116       \return value for \a key for line with ID \a id.
     117
     118       \since New in yat 0.13
     119    */
     120    //
     121    const std::string& read_group(const std::string& id,
     122                                  const std::string& key) const;
     123
     124    /**
    90125       \brief Exchanges the content in \c *this and \a other
    91126
     
    146181#endif
    147182    bam_hdr_t* header_;
     183    typedef std::map<std::string, std::string> strMap;
     184    typedef std::map<std::string, strMap> strMap2;
     185    mutable strMap2 read_group_;
     186    mutable strMap2 program_group_;
    148187
    149188    friend class InBamFile;
    150189    friend class OutBamFile;
     190
     191    const std::string& group(strMap2& map, const std::string& type,
     192                             const std::string& id,
     193                             const std::string& key) const;
    151194
    152195    // using compiler generated copy and assignment
Note: See TracChangeset for help on using the changeset viewer.