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