Changeset 3464


Ignore:
Timestamp:
Feb 9, 2016, 1:10:28 AM (6 years ago)
Author:
Peter
Message:

implementing idxstats functions against libhts. refs #852

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/bam.cc

    r3463 r3464  
    326326  InBamFile in;
    327327  in.open(file);
     328#ifdef YAT_HAVE_HTSLIB
     329  // test idxstats
     330  if (in.n_mapped(1)!=795) {
     331    suite.add(false);
     332    suite.err() << "error: InBamFile::n_mapped: "
     333                << in.n_mapped(1) << " expected: 795\n";
     334  }
     335  if (in.n_unmapped(1)!=5) {
     336    suite.add(false);
     337    suite.err() << "error: InBamFile::n_unmapped: "
     338                << in.n_unmapped(1) << " expected: 5\n";
     339  }
     340  if (in.n_no_coordinate()!=0) {
     341    suite.add(false);
     342    suite.err() << "error: InBamFile::n_no_coordinate: "
     343                << in.n_no_coordinate() << " expected: 0\n";
     344  }
     345#endif
     346
    328347  BamRead bam;
    329348  in.read(bam);
  • trunk/yat/omic/BamFile.cc

    r3463 r3464  
    5858    bam_index_destroy(index_);
    5959#endif
     60  }
     61
     62
     63  uint64_t InBamFile::get_idx_stat(int tid, bool return_mapped) const
     64  {
     65    uint64_t mapped;
     66    uint64_t unmapped;
     67#if YAT_HAVE_HTSLIB
     68    if (hts_idx_get_stat(index(), tid, &mapped, &unmapped)) {
     69      throw utility::runtime_error("InBamFile::get_idx_stat");
     70    }
     71#else
     72    assert(0 && "function not implemented");
     73#endif
     74    if (return_mapped)
     75      return mapped;
     76    return unmapped;
    6077  }
    6178
     
    8299    }
    83100    return index_;
     101  }
     102
     103
     104  uint64_t InBamFile::n_mapped(int tid) const
     105  {
     106    return get_idx_stat(tid, true);
     107  }
     108
     109
     110  uint64_t InBamFile::n_unmapped(int tid) const
     111  {
     112    return get_idx_stat(tid, false);
     113  }
     114
     115
     116  uint64_t InBamFile::n_no_coordinate(void) const
     117  {
     118#if YAT_HAVE_HTSLIB
     119    return hts_idx_get_n_no_coor(index());
     120#else
     121    assert(0 && "function not implemented");
     122    return 0;
     123#endif
    84124  }
    85125
  • trunk/yat/omic/BamFile.h

    r3363 r3464  
    171171
    172172    /**
     173       \return number of mapped reads in segment \a tid
     174     */
     175    uint64_t n_mapped(int tid) const;
     176
     177    /**
     178       Typically unmapped reads have same coordinate as its mate.
     179
     180       \return number of unmapped reads with segment \a tid
     181     */
     182    uint64_t n_unmapped(int tid) const;
     183
     184    /**
     185       \return number of reads with no coordinate
     186     */
     187    uint64_t n_no_coordinate(void) const;
     188
     189    /**
    173190       \brief Open an input bam file.
    174191
     
    201218#endif
    202219  private:
     220    uint64_t get_idx_stat(int tid, bool return_mapped) const;
    203221    BamHeader header_;
     222    // always access index_ via function index(), so index is loaded
     223    // if needed
    204224    mutable index_type* index_;
    205225  };
Note: See TracChangeset for help on using the changeset viewer.