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, 11 months ago) by gpertea
File size: 3489 byte(s)
Log Message:
adding tophat source work

Line User Rev File contents
1 gpertea 29 #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     }