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 (7 years, 10 months ago) by gpertea
File size: 3224 byte(s)
Log Message:
wip - SplicedSAMHitFactory() still not implemented

Line File contents
1 /*
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 while(curr_obs_order != VMAXINT32)
39 {
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 }