1 /*
2 * reads_of_interest.cpp
3 * cufflinks
4 *
5 * Created by Cole Trapnell on 11/29/09.
6 * Copyright 2009 Cole Trapnell. All rights reserved.
7 *
8 */
9
10 #ifdef HAVE_CONFIG_H
11 #include <config.h>
12 #endif
13
14 #include <stdlib.h>
15 #include <getopt.h>
16 #include <string>
17 #include <numeric>
18 #include <cfloat>
19 #include <iostream>
20
21 #include "common.h"
22 #include "hits.h"
23 #include "bundles.h"
24 #include "abundances.h"
25
26 #if ENABLE_THREADS
27 const char *short_options = "";
28 #else
29 const char *short_options = "";
30 #endif
31
32 static struct option long_options[] = {
33 {0, 0, 0, 0} // terminator
34 };
35
print_usage()36 void print_usage()
37 {
38 //NOTE: SPACES ONLY, bozo
39 fprintf(stderr, "Usage: cuffdiff <transcripts.gtf> <sample_hits.sam>\n");
40 }
41
parse_options(int argc,char ** argv)42 int parse_options(int argc, char** argv)
43 {
44 int option_index = 0;
45 int next_option;
46 do {
47 next_option = getopt_long(argc, argv, short_options, long_options, &option_index);
48 switch (next_option) {
49 case -1: /* Done with options. */
50 break;
51
52 default:
53 print_usage();
54 return 1;
55 }
56 } while(next_option != -1);
57
58 return 0;
59 }
60
driver(FILE * ref_gtf,FILE * sam_hit_file)61 void driver(FILE* ref_gtf, FILE* sam_hit_file)
62 {
63 ReadTable it;
64 RefSequenceTable rt(true, false);
65
66 vector<Scaffold> ref_mRNAs;
67
68 SAMHitFactory hs(it, rt);
69
70 boost::shared_ptr<HitFactory> hit_factory(createSamHitFactory(hit_file_name, it, rt));
71 BundleFactory& bundle_factory = *(new BundleFactory(hit_factory, bundle_mode));
72
73 boost::crc_32_type ref_gtf_crc_result;
74 vector<boost::shared_ptr<Scaffold> > ref_mRNAs;
75 if (ref_gtf)
76 {
77 ::load_ref_rnas(ref_gtf, bundle_factory.ref_table(), ref_mRNAs, ref_gtf_crc_result, false, false);
78 bundle_factory.set_ref_rnas(ref_mRNAs);
79 }
80 rt.print_rec_ordering();
81 vector<boost::shared_ptr<Scaffold> > mask_rnas;
82 boost::crc_32_type mask_gtf_crc_result;
83 if (mask_gtf)
84 {
85 ::load_ref_rnas(mask_gtf, bundle_factory.ref_table(), mask_rnas, mask_gtf_crc_result, false, false);
86 bundle_factory.set_mask_rnas(mask_rnas);
87 }
88
89 long num_fragments = 0;
90 long num_reads = 0;
91
92 const vector<MateHit>& hits = bundle.hits();
93 for (size_t i = 0; i < bundle.hits().size(); ++i)
94 {
95 bool compatible = false;
96 for (vector<Scaffold>::iterator ri = ref_mrnas.begin();
97 ri != ref_mrnas.end();
98 ++ri)
99 {
100 if (Scaffold::overlap_in_genome(*ri, hits[i], 0) &&
101 Scaffold::compatible(*ri,hits[i]))
102 {
103 compatible = true;
104 break;
105 }
106 }
107 if (!compatible)
108 continue;
109
110 if (hits[i].left_alignment() || hits[i].right_alignment())
111 num_fragments++;
112 if (hits[i].left_alignment())
113 {
114 printf("%s\n", hits[i].left_alignment()->hitfile_rec().c_str());
115 num_reads++;
116 }
117 if (hits[i].right_alignment())
118 {
119 printf("%s\n", hits[i].right_alignment()->hitfile_rec().c_str());
120 num_reads++;
121 }
122 }
123
124 fprintf(stderr, "Extracted %ld fragments, %ld reads\n", num_fragments, num_reads);
125 }
126
main(int argc,char ** argv)127 int main(int argc, char** argv)
128 {
129 int parse_ret = parse_options(argc,argv);
130 if (parse_ret)
131 return parse_ret;
132
133 if(optind >= argc)
134 {
135 print_usage();
136 return 1;
137 }
138
139 string ref_gtf_filename = argv[optind++];
140
141 string sam_hits_file_name = argv[optind++];
142 // Open the approppriate files
143 FILE* sam_hits_file = fopen(sam_hits_file_name.c_str(), "r");
144 if (sam_hits_file == NULL)
145 {
146 fprintf(stderr, "Error: cannot open SAM file %s for reading\n",
147 sam_hits_file_name.c_str());
148 exit(1);
149 }
150
151 FILE* ref_gtf = NULL;
152 if (ref_gtf_filename != "")
153 {
154 ref_gtf = fopen(ref_gtf_filename.c_str(), "r");
155 if (!ref_gtf)
156 {
157 fprintf(stderr, "Error: cannot open GTF file %s for reading\n",
158 ref_gtf_filename.c_str());
159 exit(1);
160 }
161 }
162
163 driver(ref_gtf, sam_hits_file);
164
165 return 0;
166 }
167