ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/sam_juncs.cpp
Revision: 135
Committed: Mon Dec 12 22:28:38 2011 UTC (8 years, 10 months ago) by gpertea
File size: 3224 byte(s)
Log Message:
wip - SplicedSAMHitFactory() still not implemented

Line User Rev File contents
1 gpertea 29 /*
2     * sam_juncs.cpp
3     * TopHat
4     *
5     * Created by Cole Trapnell on 7/5/10.
6     * Copyright 2010 Cole Trapnell. All rights reserved.
7     *
8     */
9    
10     #ifdef HAVE_CONFIG_H
11     #include <config.h>
12     #else
13     #define PACKAGE_VERSION "INTERNAL"
14     #define SVN_REVISION "XXX"
15     #endif
16    
17     #include <getopt.h>
18    
19     #include "common.h"
20     #include "bwt_map.h"
21     #include "junctions.h"
22    
23    
24     void get_junctions_from_hitstream(HitStream& hitstream,
25     ReadTable& it,
26     JunctionSet& junctions)
27     {
28     HitsForRead curr_hit_group;
29    
30    
31     hitstream.next_read_hits(curr_hit_group);
32    
33    
34     uint32_t curr_obs_order = it.observation_order(curr_hit_group.insert_id);
35    
36    
37     // While we still have unreported hits...
38 gpertea 135 while(curr_obs_order != VMAXINT32)
39 gpertea 29 {
40     for (size_t i = 0; i < curr_hit_group.hits.size(); ++i)
41     {
42     const BowtieHit& bh = curr_hit_group.hits[i];
43     junctions_from_alignment(bh, junctions);
44     }
45    
46     //fprintf(stderr, "#Hits = %d\n", curr_hit_group.hits.size());
47    
48     //curr_hit_group = HitsForRead();
49     // Get next hit group
50    
51     hitstream.next_read_hits(curr_hit_group);
52     curr_obs_order = it.observation_order(curr_hit_group.insert_id);
53     }
54    
55     hitstream.reset();
56     }
57    
58    
59     void driver(FILE* hit_map)
60     {
61     ReadTable it;
62    
63     RefSequenceTable rt(sam_header, true);
64    
65     SAMHitFactory hit_factory(it,rt);
66    
67     //HitStream hitstream(hit_map, &hit_factory, false, true, true);
68    
69     JunctionSet junctions;
70    
71     while (hit_map && !feof(hit_map))
72     {
73     char bwt_buf[2048];
74     if (!fgets(bwt_buf, 2048, hit_map))
75     {
76     break;
77     }
78     // Chomp the newline
79     char* nl = strrchr(bwt_buf, '\n');
80     if (nl) *nl = 0;
81    
82     // Get a new record from the tab-delimited Bowtie map
83     BowtieHit bh;
84    
85     if (hit_factory.get_hit_from_buf(bwt_buf, bh, false))
86     {
87     junctions_from_alignment(bh, junctions);
88    
89     }
90     }
91    
92     for (JunctionSet::iterator itr = junctions.begin();
93     itr != junctions.end();
94     ++itr)
95     {
96     const char* ref_name = rt.get_name(itr->first.refid);
97    
98     fprintf(stdout,
99     "%s\t%d\t%d\t%c\n",
100     ref_name,
101     itr->first.left - 1,
102     itr->first.right,
103     itr->first.antisense ? '-' : '+');
104     }
105    
106     fprintf(stderr, "Extracted %lu junctions\n", junctions.size());
107     }
108    
109     void print_usage()
110     {
111     fprintf(stderr, "Usage: sam_juncs <hits.sam>\n");
112    
113     // fprintf(stderr, "Usage: tophat_reports <coverage.wig> <junctions.bed> <accepted_hits.sam> <map1.bwtout> [splice_map1.sbwtout]\n");
114     }
115    
116     int main(int argc, char** argv)
117     {
118     fprintf(stderr, "sam_juncs v%s (%s)\n", PACKAGE_VERSION, SVN_REVISION);
119     fprintf(stderr, "---------------------------------------\n");
120    
121     reads_format = FASTQ;
122    
123     int parse_ret = parse_options(argc, argv, print_usage);
124     if (parse_ret)
125     return parse_ret;
126    
127     if(optind >= argc)
128     {
129     print_usage();
130     return 1;
131     }
132    
133     string map_filename = argv[optind++];
134    
135     FILE* map_file = fopen(map_filename.c_str(), "r");
136     if (!map_file)
137     {
138     fprintf(stderr, "Error: cannot open map file %s for reading\n",
139     map_filename.c_str());
140     exit(1);
141     }
142    
143     driver(map_file);
144    
145     return 0;
146     }