ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/wiggles.cpp
Revision: 29
Committed: Tue Aug 2 21:24:54 2011 UTC (8 years, 2 months ago) by gpertea
File size: 3489 byte(s)
Log Message:
adding tophat source work

Line File contents
1 #ifdef HAVE_CONFIG_H
2 #include <config.h>
3 #endif
4
5 /*
6 * wiggles.cpp
7 * TopHat
8 *
9 * Created by Cole Trapnell on 12/12/08.
10 * Copyright 2008 Cole Trapnell. All rights reserved.
11 *
12 */
13 #include <cassert>
14 #include <stdio.h>
15 #include <vector>
16 #include <string>
17 #include "wiggles.h"
18
19 using namespace std;
20
21 void print_wiggle_header(FILE* coverage_out)
22 {
23 fprintf(coverage_out, "track type=bedGraph name=\"TopHat - read coverage\"\n");
24 }
25
26 void print_wiggle_for_ref(FILE* coverage_out,
27 const string& ref_name,
28 const vector<unsigned int>& DoC)
29 {
30 unsigned short last_doc = 0; // Last DoC value we wrote to the file
31 size_t last_pos = 0; // Postition where the last written DoC came from
32 //fprintf(coverage_out, "variableStep chrom=chr%s\n", name.c_str());
33 for (size_t i = 0; i < DoC.size(); ++i)
34 {
35 if (last_doc != DoC[i])
36 {
37 size_t j = last_pos;
38 while (i - j > 10000000)
39 {
40 fprintf(coverage_out,"%s\t%d\t%d\t%d\n",ref_name.c_str(),(int)j, (int)j + 10000000, last_doc);
41 j += 10000000;
42 }
43 if (i > 0)
44 {
45 fprintf(coverage_out,"%s\t%d\t%d\t%d\n",ref_name.c_str(),(int)j, (int)i, last_doc);
46 }
47 last_pos = i;
48 last_doc = DoC[i];
49 }
50 }
51 if (last_doc)
52 {
53 fprintf(coverage_out,"%s\t%d\t%d\t%d\n",ref_name.c_str(),(int)last_pos, (int)DoC.size() -1, last_doc);
54 }
55 }
56
57 void driver(FILE* map_file, FILE* coverage_file)
58 {
59 ReadTable it;
60 RefSequenceTable rt(true);
61
62 SAMHitFactory hit_factory(it,rt);
63
64 char bwt_buf[2048];
65 bwt_buf[0] = 0;
66
67 uint32_t last_ref = 0;
68
69 vector<unsigned int> DoC;
70
71 print_wiggle_header(coverage_file);
72
73 while (map_file && !feof(map_file))
74 {
75 fgets(bwt_buf, 2048, map_file);
76 // Chomp the newline
77 char* nl = strrchr(bwt_buf, '\n');
78 if (nl) *nl = 0;
79
80 // Get a new record from the tab-delimited Bowtie map
81 BowtieHit bh;
82
83 if (hit_factory.get_hit_from_buf(bwt_buf, bh, false))
84 {
85 if (bh.ref_id() != last_ref)
86 {
87 if (last_ref != 0)
88 {
89 // print wiggle
90
91 string ref_name = rt.get_name(last_ref);
92 print_wiggle_for_ref(coverage_file,
93 ref_name,
94 DoC);
95 }
96 DoC.clear();
97 }
98
99 last_ref = bh.ref_id();
100 add_hit_to_coverage(bh, DoC);
101 }
102 }
103
104 if (last_ref != 0)
105 {
106 string ref_name = rt.get_name(last_ref);
107 print_wiggle_for_ref(coverage_file,
108 ref_name,
109 DoC);
110
111 }
112 }
113
114 void print_usage()
115 {
116 fprintf(stderr, "Usage: wiggles <accepted_hits.sam> <coverage.wig>\n");
117 }
118
119
120 int main(int argc, char** argv)
121 {
122 fprintf(stderr, "wiggles v%s (%s)\n", PACKAGE_VERSION, SVN_REVISION);
123 fprintf(stderr, "---------------------------------------\n");
124
125 int parse_ret = parse_options(argc, argv, print_usage);
126 if (parse_ret)
127 return parse_ret;
128
129 if(optind >= argc)
130 {
131 print_usage();
132 return 1;
133 }
134
135 string map_filename = argv[optind++];
136
137 if(optind >= argc)
138 {
139 print_usage();
140 return 1;
141 }
142
143 string coverage_file_name = argv[optind++];
144
145 FILE* map_file = fopen(map_filename.c_str(), "r");
146 if (!map_file)
147 {
148 fprintf(stderr, "Error: cannot open map file %s for reading\n",
149 map_filename.c_str());
150 exit(1);
151 }
152
153 // Open the approppriate files
154
155 FILE* coverage_file = fopen((coverage_file_name).c_str(), "w");
156 if (coverage_file == NULL)
157 {
158 fprintf(stderr, "Error: cannot open WIG file %s for writing\n",
159 coverage_file_name.c_str());
160 exit(1);
161 }
162
163 driver(map_file, coverage_file);
164
165 return 0;
166 }