wiggle.cpp

   1 #include <cstring>
   2 #include <cstdlib>
   3 #include <cassert>
   4
   5 #include <stdint.h>
   6 #include "sam/bam.h"
   7 #include "sam/sam.h"
   8
   9 #include "wiggle.h"
  10
  11 void add_bam_record_to_wiggle(const bam1_t *b, Wiggle& wiggle) {
  12     uint8_t *p_tag = bam_aux_get(b, "ZW");
  13     float w = (p_tag != NULL ? bam_aux2f(p_tag) : 1.0);
  14     int pos = b->core.pos;
  15     uint32_t *p = bam1_cigar(b);
  16
  17     for (int i = 0; i < (int)b->core.n_cigar; i++, ++p) {
  18         int op = *p & BAM_CIGAR_MASK;
  19         int op_len = *p >> BAM_CIGAR_SHIFT;
  20
  21         switch (op) {
  22             //case BAM_CSOFT_CLIP : pos += op_len; break;
  23         case BAM_CINS : pos += op_len; break;
  24         case BAM_CMATCH :
  25             for (int j = 0; j < op_len; j++, ++pos) {
  26                 wiggle.read_depth[pos] += w;
  27             }
  28             break;
  29         case BAM_CREF_SKIP : pos += op_len; break;
  30         default : assert(false);
  31         }
  32     }
  33 }
  34
  35 void build_wiggles(const std::string& bam_filename,
  36                    WiggleProcessor& processor) {
  37     samfile_t *bam_in = samopen(bam_filename.c_str(), "rb", NULL);
  38         if (bam_in == 0) { fprintf(stderr, "Cannot open %s!\n", bam_filename.c_str()); exit(-1); }
  39         //assert(bam_in != 0);
  40
  41     int cur_tid = -1; //current tid;
  42         int cnt = 0;
  43     bam1_t *b = bam_init1();
  44     Wiggle wiggle;
  45         while (samread(bam_in, b) >= 0) {
  46                 if (b->core.flag & 0x0004) continue;
  47
  48                 if (b->core.tid != cur_tid) {
  49                         if (cur_tid >= 0) processor.process(wiggle);
  50                         cur_tid = b->core.tid;
  51             wiggle.name = bam_in->header->target_name[cur_tid];
  52             wiggle.read_depth.assign(bam_in->header->target_len[cur_tid], 0.0);
  53                 }
  54         add_bam_record_to_wiggle(b, wiggle);
  55                 ++cnt;
  56                 if (cnt % 1000000 == 0) fprintf(stderr, "%d FIN\n", cnt);
  57         }
  58         if (cur_tid >= 0) processor.process(wiggle);
  59
  60         samclose(bam_in);
  61         bam_destroy1(b);
  62 }
  63
  64 UCSCWiggleTrackWriter::UCSCWiggleTrackWriter(const std::string& output_filename,
  65                                              const std::string& track_name) {
  66     fo = fopen(output_filename.c_str(), "w");
  67     fprintf(fo, "track type=wiggle_0 name=\"%s\" description=\"%s\" visibility=full\n",
  68             track_name.c_str(),
  69             track_name.c_str());
  70 }
  71
  72 UCSCWiggleTrackWriter::~UCSCWiggleTrackWriter() {
  73     fclose(fo);
  74 }
  75
  76 void UCSCWiggleTrackWriter::process(const Wiggle& wiggle) {
  77     int sp, ep;
  78
  79     sp = ep = -1;
  80     for (size_t i = 0; i < wiggle.read_depth.size(); i++) {
  81         if (wiggle.read_depth[i] > 0) {
  82             ep = i;
  83         }
  84         else {
  85             if (sp < ep) {
  86                 ++sp;
  87                 fprintf(fo, "fixedStep chrom=%s start=%d step=1\n", wiggle.name.c_str(), sp + 1);
  88                 for (int j = sp; j <= ep; j++) fprintf(fo, "%.7g\n", wiggle.read_depth[j]);
  89             }
  90             sp = i;
  91         }
  92     }
  93     if (sp < ep) {
  94         ++sp;
  95         fprintf(fo, "fixedStep chrom=%s start=%d step=1\n", wiggle.name.c_str(), sp + 1);
  96         for (int j = sp; j <= ep; j++) fprintf(fo, "%.7g\n", wiggle.read_depth[j]);
  97     }
  98 }
  99
 100 ReadDepthWriter::ReadDepthWriter(std::ostream& stream)
 101     : stream_(stream) {
 102 }
 103
 104 void ReadDepthWriter::process(const Wiggle& wiggle) {
 105     stream_ << wiggle.name << '\t'
 106             << wiggle.read_depth.size() << '\t';
 107     for (size_t i = 0; i < wiggle.read_depth.size(); ++i) {
 108         if (i > 0) stream_ << ' ';
 109         stream_ << wiggle.read_depth[i];
 110     }
 111     stream_ << '\n';
 112 }