wiggle.cpp

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