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
  40         bam_header_t *header = bam_in->header;
  41         bool *used = new bool[header->n_targets];
  42         memset(used, 0, sizeof(bool) * header->n_targets);
  43
  44     int cur_tid = -1; //current tid;
  45         int cnt = 0;
  46     bam1_t *b = bam_init1();
  47     Wiggle wiggle;
  48         while (samread(bam_in, b) >= 0) {
  49                 if (b->core.flag & 0x0004) continue;
  50
  51                 if (b->core.tid != cur_tid) {
  52                         if (cur_tid >= 0) { used[cur_tid] = true; processor.process(wiggle); }
  53                         cur_tid = b->core.tid;
  54             wiggle.name = header->target_name[cur_tid];
  55             wiggle.length = header->target_len[cur_tid];
  56             wiggle.read_depth.assign(wiggle.length, 0.0);
  57                 }
  58         add_bam_record_to_wiggle(b, wiggle);
  59                 ++cnt;
  60                 if (cnt % 1000000 == 0) fprintf(stderr, "%d FIN\n", cnt);
  61         }
  62         if (cur_tid >= 0) { used[cur_tid] = true; processor.process(wiggle); }
  63
  64         for (int32_t i = 0; i < header->n_targets; i++)
  65                 if (!used[i]) {
  66                         wiggle.name = header->target_name[i];
  67                         wiggle.length = header->target_len[i];
  68                         wiggle.read_depth.clear();
  69                         processor.process(wiggle);
  70                 }
  71
  72         samclose(bam_in);
  73         bam_destroy1(b);
  74         delete[] used;
  75 }
  76
  77 UCSCWiggleTrackWriter::UCSCWiggleTrackWriter(const std::string& output_filename,
  78                                              const std::string& track_name) {
  79     fo = fopen(output_filename.c_str(), "w");
  80     fprintf(fo, "track type=wiggle_0 name=\"%s\" description=\"%s\" visibility=full\n",
  81             track_name.c_str(),
  82             track_name.c_str());
  83 }
  84
  85 UCSCWiggleTrackWriter::~UCSCWiggleTrackWriter() {
  86     fclose(fo);
  87 }
  88
  89 void UCSCWiggleTrackWriter::process(const Wiggle& wiggle) {
  90     int sp, ep;
  91
  92     if (wiggle.read_depth.empty()) return;
  93
  94     sp = ep = -1;
  95     for (size_t i = 0; i < wiggle.length; i++) {
  96         if (wiggle.read_depth[i] > 0) {
  97             ep = i;
  98         }
  99         else {
 100             if (sp < ep) {
 101                 ++sp;
 102                 fprintf(fo, "fixedStep chrom=%s start=%d step=1\n", wiggle.name.c_str(), sp + 1);
 103                 for (int j = sp; j <= ep; j++) fprintf(fo, "%.7g\n", wiggle.read_depth[j]);
 104             }
 105             sp = i;
 106         }
 107     }
 108     if (sp < ep) {
 109         ++sp;
 110         fprintf(fo, "fixedStep chrom=%s start=%d step=1\n", wiggle.name.c_str(), sp + 1);
 111         for (int j = sp; j <= ep; j++) fprintf(fo, "%.7g\n", wiggle.read_depth[j]);
 112     }
 113 }
 114
 115 ReadDepthWriter::ReadDepthWriter(std::ostream& stream)
 116     : stream_(stream) {
 117 }
 118
 119 void ReadDepthWriter::process(const Wiggle& wiggle) {
 120
 121     stream_ << wiggle.name << '\t'
 122             << wiggle.length << '\t';
 123
 124     if (wiggle.read_depth.empty()) { stream_ << "NA\n"; return; }
 125
 126     for (size_t i = 0; i < wiggle.length; ++i) {
 127         if (i > 0) stream_ << ' ';
 128         stream_ << wiggle.read_depth[i];
 129     }
 130     stream_ << '\n';
 131 }