13 bool no_fractional_weight = false;
15 void add_bam_record_to_wiggle(const bam1_t *b, Wiggle& wiggle) {
18 if (no_fractional_weight) w = 1.0;
20 uint8_t *p_tag = bam_aux_get(b, "ZW");
21 if (p_tag == NULL) return;
25 int pos = b->core.pos;
26 uint32_t *p = bam1_cigar(b);
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;
33 //case BAM_CSOFT_CLIP : pos += op_len; break;
34 case BAM_CINS : pos += op_len; break;
36 for (int j = 0; j < op_len; j++, ++pos) {
37 wiggle.read_depth[pos] += w;
40 case BAM_CREF_SKIP : pos += op_len; break;
41 default : assert(false);
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); }
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);
55 int cur_tid = -1; //current tid;
57 bam1_t *b = bam_init1();
59 while (samread(bam_in, b) >= 0) {
60 if (b->core.flag & 0x0004) continue;
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);
69 add_bam_record_to_wiggle(b, wiggle);
71 if (cnt % 1000000 == 0) std::cout<< cnt<< std::endl;
73 if (cur_tid >= 0) { used[cur_tid] = true; processor.process(wiggle); }
75 for (int32_t i = 0; i < header->n_targets; i++)
77 wiggle.name = header->target_name[i];
78 wiggle.length = header->target_len[i];
79 wiggle.read_depth.clear();
80 processor.process(wiggle);
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",
96 UCSCWiggleTrackWriter::~UCSCWiggleTrackWriter() {
100 void UCSCWiggleTrackWriter::process(const Wiggle& wiggle) {
103 if (wiggle.read_depth.empty()) return;
106 for (size_t i = 0; i < wiggle.length; i++) {
107 if (wiggle.read_depth[i] > 0) {
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, "%.7g\n", wiggle.read_depth[j]);
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, "%.7g\n", wiggle.read_depth[j]);
126 ReadDepthWriter::ReadDepthWriter(std::ostream& stream)
130 void ReadDepthWriter::process(const Wiggle& wiggle) {
132 stream_ << wiggle.name << '\t'
133 << wiggle.length << '\t';
135 if (wiggle.read_depth.empty()) { stream_ << "NA\n"; return; }
137 for (size_t i = 0; i < wiggle.length; ++i) {
138 if (i > 0) stream_ << ' ';
139 stream_ << wiggle.read_depth[i];