gclib/tophat_cpp/library_stats.cpp

/*
 *  library_stats.cpp
 *  TopHat
 *
 *  Created by Cole Trapnell on 11/17/08.
 *  Copyright 2008 Cole Trapnell. All rights reserved.
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <cassert>
#include <cstdio>
#include <vector>
#include <string>
#include <map>
#include <algorithm>
#include <set>
#include <stdexcept>
#include <iostream>
#include <fstream>
#include <getopt.h>

#include "common.h"
#include "bwt_map.h"
#include "inserts.h"
#include "tokenize.h"

using namespace std;

struct LibraryStats
{
        LibraryStats() : 
                m1_aligned_reads(0), 
                m2_aligned_reads(0), 
                m1_singletons(0), 
                m2_singletons(0), 
                spliced_reads(0),
                too_far_inserts(0), 
                too_close_inserts(0), 
                same_strand_inserts(0), 
                happy_inserts(0) {}
        
        // A name for this statistics object
        string ref_region;
        
        // Read-level alignment counters
        uint32_t m1_aligned_reads;
        uint32_t m2_aligned_reads;
        uint32_t m1_singletons;
        uint32_t m2_singletons;
        uint32_t spliced_reads;
        
        // Insert category counts
        uint32_t too_far_inserts;
        uint32_t too_close_inserts;
        uint32_t same_strand_inserts;
        uint32_t happy_inserts;
        

        // Roll the stats for rhs into this object
        LibraryStats& operator+=(const LibraryStats& rhs)
        {
                m1_aligned_reads += rhs.m1_aligned_reads;
                m2_aligned_reads += rhs.m2_aligned_reads;
                m1_singletons += rhs.m1_singletons;
                m2_singletons += rhs.m2_singletons;
                
                too_far_inserts += rhs.too_far_inserts;
                too_close_inserts += rhs.too_close_inserts;
                same_strand_inserts += rhs.same_strand_inserts;
                happy_inserts += rhs.happy_inserts;
                
                spliced_reads += rhs.spliced_reads;
                
                return *this;
        }
                
        uint32_t contiguously_aligned_reads() const
        {
                return m2_aligned_reads + m1_aligned_reads - spliced_reads;
        }
};


ostream& operator<<(ostream& os, const LibraryStats& L)
{
        os << L.ref_region << ":" << endl;
        os << "\t" << "Aligned reads:\t" << L.m2_aligned_reads + L.m1_aligned_reads << endl;
        os << "\t" << "Singletons:\t" << L.m2_singletons + L.m1_singletons << endl;
        os << "\t" << "Happy reads:\t" << 2 * L.happy_inserts << endl;
        os << "\t" << "Unhappy inserts, too far:\t" << 2 * L.too_far_inserts << endl;
        os << "\t" << "Unhappy inserts, too close:\t" << 2 * L.too_close_inserts << endl;
        os << "\t" << "Unhappy inserts, same strand:\t" << 2 * L.same_strand_inserts << endl;
        os << "\t" << "Spliced read alignments:\t" << L.spliced_reads << endl;
        os << "\t" << "Contiguous read alignments:\t" << L.contiguously_aligned_reads() << endl;
        
        return os;
}


// Print a simple listing of inner distances on the best insert mappings.
// Useful if you want to build a histogram of insert sizes, or determine 
// the mean and std dev of your library insert size empirically.
//void print_mate_distances(const string& dist_file_name,
//                                                const vector<InsertStatus>& best_status_for_inserts)
//{
//      ofstream dist_file(dist_file_name.c_str());
//      for (size_t i = 0; i < best_status_for_inserts.size(); ++i)
//      {
//              if (best_status_for_inserts[i].mask < SINGLETON)
//              {
//                      dist_file << best_status_for_inserts[i].inner_dist << endl;
//              }
//      }
//}

// Compute some statistics about the mapping for this library, and return them
// in a LibraryStats object.


LibraryStats compute_stats(const BestInsertAlignmentTable best_status_for_inserts)
{
        LibraryStats best_pairing_stats;
        //best_pairing_stats.ref_region = "Best alignment pairings";

        // Analyze the table of best pairings for each insert and report the same 
        // statistics
        for (size_t i = 0; i < best_status_for_inserts.size(); ++i)
        {
                const pair<InsertAlignmentGrade, vector<InsertAlignment> >& st = best_status_for_inserts[i];
                
                for (size_t j = 0; j < st.second.size(); ++j)
                {
                        
                        InsertAlignment a = st.second[j];
                        if (a.left_alignment && a.right_alignment)
                        {
                                pair<int, int> distances = pair_distances(*(a.left_alignment),
                                                                                                                  *(a.right_alignment));
                                int inner_dist = distances.second;
                                printf("%d\n", inner_dist);
                        }
                }
        }
        return best_pairing_stats;
}

void load_hits(const vector<string>& filenames,
                           RefSequenceTable& rt,
                           ReadTable& it,
               HitTable& hits)
{
    for (size_t i = 0; i < filenames.size(); ++i)
    {
        const string& filename = filenames[i];
                HitFactory* hit_factory = NULL;
                
        if (filename.rfind(".sam") != string::npos)
                {
                        hit_factory = new SAMHitFactory(it,rt);
                }
                else
                {
                        hit_factory = new BowtieHitFactory(it,rt);
                }
                
                
        FILE* map = fopen(filename.c_str(), "r");
                
        if (map == NULL)
        {
            fprintf(stderr, "Error: could not open %s\n", filename.c_str());
            exit(1);
        }
        fprintf(stderr, "Loading hits from %s\n", filename.c_str());
        size_t num_hits_before_load = hits.total_hits();
        get_mapped_reads(map, hits, *hit_factory, true);
        fprintf(stderr, "Loaded %d hits from %s\n", 
                (int)hits.total_hits() - (int)num_hits_before_load, 
                filename.c_str());
                delete hit_factory;
    }
}

void load_hits(RefSequenceTable& rt,
                           ReadTable& it,
                           const string& left_read_maplist,
               HitTable& left_hits,
               const string& right_read_maplist,
               HitTable& right_hits)
{
    vector<string> left_filenames;
    tokenize(left_read_maplist, ",", left_filenames);
    load_hits(left_filenames, rt, it, left_hits);
        vector<string> right_filenames;
        tokenize(right_read_maplist, ",", right_filenames);
        load_hits(right_filenames, rt, it, right_hits);
}

void driver(const string& left_maps,
                        const string& right_maps)
{
        // Load the set of left maps, and if provided, the set of right maps
    ReadTable it;
    RefSequenceTable rt(true);
        
        //HitFactory hit_factory(it,rt);
        
    HitTable left_hits;
    HitTable right_hits;
        
        load_hits(rt, it, left_maps, left_hits, right_maps, right_hits);
        
        BestInsertAlignmentTable best_pairings(it.size() + 1);
        
        insert_best_pairings(rt, it, left_hits, right_hits, best_pairings, true);
        
        // Print a simple listing of inner distances on the best insert mappings.
        // Useful if you want to build a histogram of insert sizes, or determine 
        // the mean and std dev of your library insert size empirically.
        //print_mate_distances("mate_dist.txt", insert_best_pairings);
        
        // Compute and report some statistics about the mapping for this library.
        LibraryStats best_pairing_stats = compute_stats(best_pairings);
        //cout << best_pairing_stats;
}


void print_usage()
{
    fprintf(stderr, "Usage:   library_stats <left_map1,...,left_mapN> <right_map1,...,right_mapN>\n");
        
        //    fprintf(stderr, "Usage:   tophat_reports <coverage.wig> <junctions.bed> <accepted_hits.sam> <map1.bwtout> [splice_map1.sbwtout]\n");
}

int main(int argc, char** argv)
{       
        int parse_ret = parse_options(argc, argv, print_usage);
        if (parse_ret)
                return parse_ret;
        
    if(optind >= argc) 
    {
                print_usage();
                return 1;
        }
        
        string map1_file_name = argv[optind++];
        
        if(optind >= argc) 
        {
                print_usage();
                return 1;
        }
        
        string map2_file_name = argv[optind++];

    driver(map1_file_name, map2_file_name);
    
    return 0;
}
Revision:	29
Committed:	Tue Aug 2 21:24:54 2011 UTC (13 years, 1 month ago) by gpertea
File size:	7116 byte(s)
Log Message:	adding tophat source work
Line	File contents
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	}