wiggle.cpp

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