1 |
/* |
2 |
* library_stats.cpp |
3 |
* TopHat |
4 |
* |
5 |
* Created by Cole Trapnell on 11/17/08. |
6 |
* Copyright 2008 Cole Trapnell. All rights reserved. |
7 |
* |
8 |
*/ |
9 |
|
10 |
#ifdef HAVE_CONFIG_H |
11 |
#include <config.h> |
12 |
#endif |
13 |
|
14 |
#include <cassert> |
15 |
#include <cstdio> |
16 |
#include <vector> |
17 |
#include <string> |
18 |
#include <map> |
19 |
#include <algorithm> |
20 |
#include <set> |
21 |
#include <stdexcept> |
22 |
#include <iostream> |
23 |
#include <fstream> |
24 |
#include <getopt.h> |
25 |
|
26 |
#include "common.h" |
27 |
#include "bwt_map.h" |
28 |
#include "inserts.h" |
29 |
#include "tokenize.h" |
30 |
|
31 |
using namespace std; |
32 |
|
33 |
struct LibraryStats |
34 |
{ |
35 |
LibraryStats() : |
36 |
m1_aligned_reads(0), |
37 |
m2_aligned_reads(0), |
38 |
m1_singletons(0), |
39 |
m2_singletons(0), |
40 |
spliced_reads(0), |
41 |
too_far_inserts(0), |
42 |
too_close_inserts(0), |
43 |
same_strand_inserts(0), |
44 |
happy_inserts(0) {} |
45 |
|
46 |
// A name for this statistics object |
47 |
string ref_region; |
48 |
|
49 |
// Read-level alignment counters |
50 |
uint32_t m1_aligned_reads; |
51 |
uint32_t m2_aligned_reads; |
52 |
uint32_t m1_singletons; |
53 |
uint32_t m2_singletons; |
54 |
uint32_t spliced_reads; |
55 |
|
56 |
// Insert category counts |
57 |
uint32_t too_far_inserts; |
58 |
uint32_t too_close_inserts; |
59 |
uint32_t same_strand_inserts; |
60 |
uint32_t happy_inserts; |
61 |
|
62 |
|
63 |
|
64 |
// Roll the stats for rhs into this object |
65 |
LibraryStats& operator+=(const LibraryStats& rhs) |
66 |
{ |
67 |
m1_aligned_reads += rhs.m1_aligned_reads; |
68 |
m2_aligned_reads += rhs.m2_aligned_reads; |
69 |
m1_singletons += rhs.m1_singletons; |
70 |
m2_singletons += rhs.m2_singletons; |
71 |
|
72 |
too_far_inserts += rhs.too_far_inserts; |
73 |
too_close_inserts += rhs.too_close_inserts; |
74 |
same_strand_inserts += rhs.same_strand_inserts; |
75 |
happy_inserts += rhs.happy_inserts; |
76 |
|
77 |
spliced_reads += rhs.spliced_reads; |
78 |
|
79 |
return *this; |
80 |
} |
81 |
|
82 |
uint32_t contiguously_aligned_reads() const |
83 |
{ |
84 |
return m2_aligned_reads + m1_aligned_reads - spliced_reads; |
85 |
} |
86 |
}; |
87 |
|
88 |
|
89 |
ostream& operator<<(ostream& os, const LibraryStats& L) |
90 |
{ |
91 |
os << L.ref_region << ":" << endl; |
92 |
os << "\t" << "Aligned reads:\t" << L.m2_aligned_reads + L.m1_aligned_reads << endl; |
93 |
os << "\t" << "Singletons:\t" << L.m2_singletons + L.m1_singletons << endl; |
94 |
os << "\t" << "Happy reads:\t" << 2 * L.happy_inserts << endl; |
95 |
os << "\t" << "Unhappy inserts, too far:\t" << 2 * L.too_far_inserts << endl; |
96 |
os << "\t" << "Unhappy inserts, too close:\t" << 2 * L.too_close_inserts << endl; |
97 |
os << "\t" << "Unhappy inserts, same strand:\t" << 2 * L.same_strand_inserts << endl; |
98 |
os << "\t" << "Spliced read alignments:\t" << L.spliced_reads << endl; |
99 |
os << "\t" << "Contiguous read alignments:\t" << L.contiguously_aligned_reads() << endl; |
100 |
|
101 |
return os; |
102 |
} |
103 |
|
104 |
|
105 |
// Print a simple listing of inner distances on the best insert mappings. |
106 |
// Useful if you want to build a histogram of insert sizes, or determine |
107 |
// the mean and std dev of your library insert size empirically. |
108 |
//void print_mate_distances(const string& dist_file_name, |
109 |
// const vector<InsertStatus>& best_status_for_inserts) |
110 |
//{ |
111 |
// ofstream dist_file(dist_file_name.c_str()); |
112 |
// for (size_t i = 0; i < best_status_for_inserts.size(); ++i) |
113 |
// { |
114 |
// if (best_status_for_inserts[i].mask < SINGLETON) |
115 |
// { |
116 |
// dist_file << best_status_for_inserts[i].inner_dist << endl; |
117 |
// } |
118 |
// } |
119 |
//} |
120 |
|
121 |
// Compute some statistics about the mapping for this library, and return them |
122 |
// in a LibraryStats object. |
123 |
|
124 |
|
125 |
|
126 |
LibraryStats compute_stats(const BestInsertAlignmentTable best_status_for_inserts) |
127 |
{ |
128 |
LibraryStats best_pairing_stats; |
129 |
//best_pairing_stats.ref_region = "Best alignment pairings"; |
130 |
|
131 |
// Analyze the table of best pairings for each insert and report the same |
132 |
// statistics |
133 |
for (size_t i = 0; i < best_status_for_inserts.size(); ++i) |
134 |
{ |
135 |
const pair<InsertAlignmentGrade, vector<InsertAlignment> >& st = best_status_for_inserts[i]; |
136 |
|
137 |
for (size_t j = 0; j < st.second.size(); ++j) |
138 |
{ |
139 |
|
140 |
InsertAlignment a = st.second[j]; |
141 |
if (a.left_alignment && a.right_alignment) |
142 |
{ |
143 |
pair<int, int> distances = pair_distances(*(a.left_alignment), |
144 |
*(a.right_alignment)); |
145 |
int inner_dist = distances.second; |
146 |
printf("%d\n", inner_dist); |
147 |
} |
148 |
} |
149 |
} |
150 |
return best_pairing_stats; |
151 |
} |
152 |
|
153 |
void load_hits(const vector<string>& filenames, |
154 |
RefSequenceTable& rt, |
155 |
ReadTable& it, |
156 |
HitTable& hits) |
157 |
{ |
158 |
for (size_t i = 0; i < filenames.size(); ++i) |
159 |
{ |
160 |
const string& filename = filenames[i]; |
161 |
HitFactory* hit_factory = NULL; |
162 |
|
163 |
if (filename.rfind(".sam") != string::npos) |
164 |
{ |
165 |
hit_factory = new SAMHitFactory(it,rt); |
166 |
} |
167 |
else |
168 |
{ |
169 |
hit_factory = new BowtieHitFactory(it,rt); |
170 |
} |
171 |
|
172 |
|
173 |
FILE* map = fopen(filename.c_str(), "r"); |
174 |
|
175 |
if (map == NULL) |
176 |
{ |
177 |
fprintf(stderr, "Error: could not open %s\n", filename.c_str()); |
178 |
exit(1); |
179 |
} |
180 |
fprintf(stderr, "Loading hits from %s\n", filename.c_str()); |
181 |
size_t num_hits_before_load = hits.total_hits(); |
182 |
get_mapped_reads(map, hits, *hit_factory, true); |
183 |
fprintf(stderr, "Loaded %d hits from %s\n", |
184 |
(int)hits.total_hits() - (int)num_hits_before_load, |
185 |
filename.c_str()); |
186 |
delete hit_factory; |
187 |
} |
188 |
} |
189 |
|
190 |
void load_hits(RefSequenceTable& rt, |
191 |
ReadTable& it, |
192 |
const string& left_read_maplist, |
193 |
HitTable& left_hits, |
194 |
const string& right_read_maplist, |
195 |
HitTable& right_hits) |
196 |
{ |
197 |
vector<string> left_filenames; |
198 |
tokenize(left_read_maplist, ",", left_filenames); |
199 |
load_hits(left_filenames, rt, it, left_hits); |
200 |
vector<string> right_filenames; |
201 |
tokenize(right_read_maplist, ",", right_filenames); |
202 |
load_hits(right_filenames, rt, it, right_hits); |
203 |
} |
204 |
|
205 |
void driver(const string& left_maps, |
206 |
const string& right_maps) |
207 |
{ |
208 |
// Load the set of left maps, and if provided, the set of right maps |
209 |
ReadTable it; |
210 |
RefSequenceTable rt(true); |
211 |
|
212 |
//HitFactory hit_factory(it,rt); |
213 |
|
214 |
HitTable left_hits; |
215 |
HitTable right_hits; |
216 |
|
217 |
load_hits(rt, it, left_maps, left_hits, right_maps, right_hits); |
218 |
|
219 |
BestInsertAlignmentTable best_pairings(it.size() + 1); |
220 |
|
221 |
insert_best_pairings(rt, it, left_hits, right_hits, best_pairings, true); |
222 |
|
223 |
// Print a simple listing of inner distances on the best insert mappings. |
224 |
// Useful if you want to build a histogram of insert sizes, or determine |
225 |
// the mean and std dev of your library insert size empirically. |
226 |
//print_mate_distances("mate_dist.txt", insert_best_pairings); |
227 |
|
228 |
// Compute and report some statistics about the mapping for this library. |
229 |
LibraryStats best_pairing_stats = compute_stats(best_pairings); |
230 |
//cout << best_pairing_stats; |
231 |
} |
232 |
|
233 |
|
234 |
void print_usage() |
235 |
{ |
236 |
fprintf(stderr, "Usage: library_stats <left_map1,...,left_mapN> <right_map1,...,right_mapN>\n"); |
237 |
|
238 |
// fprintf(stderr, "Usage: tophat_reports <coverage.wig> <junctions.bed> <accepted_hits.sam> <map1.bwtout> [splice_map1.sbwtout]\n"); |
239 |
} |
240 |
|
241 |
int main(int argc, char** argv) |
242 |
{ |
243 |
int parse_ret = parse_options(argc, argv, print_usage); |
244 |
if (parse_ret) |
245 |
return parse_ret; |
246 |
|
247 |
if(optind >= argc) |
248 |
{ |
249 |
print_usage(); |
250 |
return 1; |
251 |
} |
252 |
|
253 |
string map1_file_name = argv[optind++]; |
254 |
|
255 |
if(optind >= argc) |
256 |
{ |
257 |
print_usage(); |
258 |
return 1; |
259 |
} |
260 |
|
261 |
string map2_file_name = argv[optind++]; |
262 |
|
263 |
driver(map1_file_name, map2_file_name); |
264 |
|
265 |
return 0; |
266 |
} |