File: stream.h

package info (click to toggle)
anfo 0.98-5
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 2,256 kB
  • sloc: sh: 10,041; cpp: 8,548; ansic: 530; lisp: 337; makefile: 64
file content (898 lines) | stat: -rw-r--r-- 30,889 bytes parent folder | download | duplicates (5)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
//    Copyright 2009 Udo Stenzel
//    This file is part of ANFO
//
//    ANFO is free software: you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation, either version 3 of the License, or
//    (at your option) any later version.
//
//    Anfo is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.
//
//    You should have received a copy of the GNU General Public License
//    along with Anfo.  If not, see <http://www.gnu.org/licenses/>.

#ifndef INCLUDED_STREAM_H
#define INCLUDED_STREAM_H

#include "compress_stream.h"
#include "logdom.h"
#include "output.pb.h"
#include "sequence.h"
#include "util.h"

#include <google/protobuf/text_format.h>
#include <google/protobuf/io/coded_stream.h>
#include <google/protobuf/io/zero_copy_stream.h>
#include <google/protobuf/io/zero_copy_stream_impl.h>

#include <iostream>

#include <algorithm>
#include <deque>
#include <fstream>
#include <limits>
#include <memory>
#include <vector>

#include <sys/types.h>
#include <sys/wait.h>

// Including Elk defines some macros that collide with protobuf.  We
// undefine them (and hope they aren't needed...).

#include <elk/scheme.h>
#undef Print
#undef MAX_TYPE

namespace streams {
	using namespace output ;
	using namespace std ;

inline uint32_t mk_msg_tag( uint32_t i ) { return i << 3 | 2 ; }

//! \defgroup outputfile Convenience functions to handle output files
//! The output file is a series of protocol buffer messages.  This makes
//! it compact and extensible.  Instead of writing a serial file, the
//! very same messages could be entered into a Berkeley-DB type
//! database.  Since conversion to all sorts of text files is a no
//! brainer, the binary format is not considered undue hardship on the
//! users of scripting languages.
//!
//! (If you are a user of a glorified logfile munching language and want
//! an ASCII-based table, define a sensible layout and ask for it.  Or
//! write the single function yourself, it isn't that hard.)
//!
//! The protobuf format itself is not self-delimiting, therefore we
//! prepend each message with the number of bytes it occupies, in
//! variable protobuf format.  The first message in a file
//! one is of type output::Header and should contain all sorts of
//! meta information, maybe even configuration of the mapper that was
//! used to create it (an option would be to include the full set of
//! policies).  Subsequent messages are of type output::Result.
//!
//! \note The output file is essentially a single, giant protobuf message.
//!       However, it is impossible to read it in one go (due to memory
//!       constraints), but it's also impossible to read it through a single
//!       CodedInputStream.  That's why methods in here tend to construct and
//!       destruct a fresh CodedInputStream on each invocation.
//!
//! @{

template< typename Msg >
void write_delimited_message( google::protobuf::io::CodedOutputStream& os, int tag, const Msg& m )
{
	os.WriteTag( mk_msg_tag( tag ) ) ;
	os.WriteVarint32( m.ByteSize() ) ;
	if( !m.SerializeToCodedStream( &os ) )
    {
        string s ;
        google::protobuf::TextFormat::PrintToString( m, &s ) ;
        cerr << s ;
        throw "error while serializing: " ;
    }
}

template< typename Msg >
void write_delimited_message( google::protobuf::io::ZeroCopyOutputStream *os, int tag, const Msg& m )
{
	google::protobuf::io::CodedOutputStream o( os ) ;
	write_delimited_message( o, tag, m ) ;
}

template< typename Msg >
bool read_delimited_message( google::protobuf::io::CodedInputStream& is, Msg &m )
{
	uint32_t size ;
	std::string code ;
	if( !is.ReadVarint32( &size ) ) return false ;
	int lim = is.PushLimit( size ) ;
	if( !m.ParseFromCodedStream( &is ) ) throw "error while deserializing" ;
	is.PopLimit( lim ) ;
	return true ;
}

void sanitize( Header& ) ;
void sanitize( Result& ) ;
void merge_sensibly( output::Header& lhs, const output::Header& rhs ) ;
void merge_sensibly( output::Footer& lhs, const output::Footer& rhs ) ;
void merge_sensibly( output::Result& lhs, const output::Result& rhs ) ;


//! \brief returns the hit to some genome
//! If an empty genome is asked for, returns the best hit.  Behaviour is
//! undefined if no suitable hit exists.
const output::Hit* hit_to( const output::Result& ) ;
const output::Hit* hit_to( const output::Result&, const string& ) ;
inline const output::Hit* hit_to( const output::Result& r, const char* g ) { return g ? hit_to( r ) : hit_to( r, string(g) ) ; }

template< typename I > const output::Hit* hit_to( const output::Result& r, I begin, I end )
{
	for( ; begin != end ; ++begin )
		if( const output::Hit* h = hit_to( r, *begin ) ) return h ;
	return 0 ;
}

//! \brief returns the mutable hit to some genome
//! If an empty genome is asked for, returns the best hit.  If no
//! suitable hit exists, a new one is created.
output::Hit* mutable_hit_to( output::Result* ) ;
output::Hit* mutable_hit_to( output::Result*, const string& ) ;

inline output::Hit::Operation cigar_op( uint32_t c ) { return (output::Hit::Operation)(c & 0xf) ; }
inline uint32_t cigar_len( uint32_t c ) { return c >> 4 ; }

//! \brief computes (trimmed) query length from CIGAR line
template< typename C > unsigned len_from_bin_cigar( const C& cig )
{
	unsigned l = 0 ;
	for( typename C::const_iterator i = cig.begin(), e = cig.end() ; i != e ; ++i )
	{
		switch( cigar_op( *i ) )
		{
			case output::Hit::Match:
			case output::Hit::Mismatch:
			case output::Hit::Insert:
			case output::Hit::SoftClip:
				l += cigar_len( *i ) ;
				break ;

			case output::Hit::Delete:
			case output::Hit::Skip:
			case output::Hit::HardClip:
			case output::Hit::Pad:
				break ;
		}
	}
	return l ;
}

inline uint32_t mk_cigar( output::Hit::Operation op, uint32_t len ) { return len << 4 | op ; }

inline void push_op( std::vector<unsigned>& s, unsigned m, output::Hit::Operation op )
{
	if( m && !s.empty() && (streams::cigar_op( s.back() ) == op) && streams::cigar_len( s.back() ) ) 
		s.back() = streams::mk_cigar( op, m + streams::cigar_len( s.back() ) ) ;
	else if( m ) s.push_back( streams::mk_cigar( op, m ) ) ;
}
inline void push_m( std::vector<unsigned>& s, unsigned m ) { push_op( s, m, output::Hit::Match ) ; }
inline void push_M( std::vector<unsigned>& s, unsigned m ) { push_op( s, m, output::Hit::Mismatch ) ; }
inline void push_i( std::vector<unsigned>& s, unsigned i ) { push_op( s, i, output::Hit::Insert ) ; }
inline void push_d( std::vector<unsigned>& s, unsigned d ) { push_op( s, d, output::Hit::Delete ) ; }

inline int effective_length( const Read& rd )
{ return (rd.has_trim_right() ? rd.trim_right() : rd.sequence().length()) - rd.trim_left() ; }

//! @}

//! \brief stream of result messages
//! Each stream is a header, followed by many results, followed by a
//! single footer.  The header will be cached internally, so it can be
//! asked for repeatedly.  Results are forgotten once read, the footer
//! is only available after all results have been read, and then it is
//! stored and can be read repeatedly.
//!
//! This class is both an input stream, an output stream and a stream
//! transducer.  Concrete implementations decide what sets of methods to
//! actually implement.  Streams behave as state machine, cycling
//! between need_input, have_output and end_of_stream.  It's undefined
//! behaviour to call methods in states where they don't make sense.

class Stream
{
	public:
		int refcount_ ;
		enum state { invalid, end_of_stream, need_input, have_output } ;

		static void cleanup( Stream* p ) { delete p ; }

	protected:
		// internal state---not strictly necessary here, but used almost
		// everywhere and therefore simply convenient
		Header hdr_ ;
		Result res_ ;
		Footer foot_ ;
		state state_ ;

		virtual ~Stream() {} 				// want control over instantiation

	private:
		Stream( const Stream& ) ; 			// must not copy
		void operator = ( const Stream& ) ;	// must not copy

	public:
		Stream() : refcount_(0), state_( end_of_stream ) {}

		//! \brief returns stream state
		//! The state determines which methods may be called: get_header
		//! can be called unless the satet is \c invalid, get_result can
		//! be called in state have_output, put_result only in
		//! need_input, get_footer only in end_of_stream.  put_footer
		//! can only be called in need_input and it will signal the end
		//! of the input stream (if there is such a thing).  put_header
		//! can only be called in state invalid, it may change the state
		//! to either need_input or have_output.  Particular streams may
		//! need special initialization before the state bacomes
		//! something other than invalid.
		virtual state get_state() { return state_ ; }

		//! \brief returns the header
		//! The header can be requested any time, unless the stream is
		//! in the invalid state.
		virtual Header fetch_header() { return hdr_ ; }

		//! \brief reads the next result
		//! Every result can only be read once, internal iterator style.
		//! If the stream state is have_output, this method returns the
		//! next result.  Otherwise the behavior is undefined.
		virtual Result fetch_result() { state_ = need_input ; return res_ ; }

		//! \brief returns the footer
		//! Only after all results have been consumed is the footer
		//! available.  If anything goes wrong internally, the LSB
		//! should be set in \c exit_code.
		virtual Footer fetch_footer() { return foot_ ; }

		//! \brief sets the stream header
		//! Output streams and stream filters need the header to become
		//! valid streams.
		virtual void put_header( const Header& h ) { hdr_ = h ; state_ = need_input ; }

		//! \brief outputs one result
		//! This method can only be called in state need_input, and it
		//! will insert a result record into the stream.  Else the
		//! behaviour is undefined.
		virtual void put_result( const Result& r ) { res_ = r ; state_ = have_output ; }

		//! \brief sets the footer
		//! This method can only be called in state need_input, and it
		//! signals that no more input is available.  A filter will then
		//! flush internal buffers and signal end_of_stream.
		virtual void put_footer( const Footer& f ) { foot_ = f ; state_ = end_of_stream ; }

		//! \brief get the 'summary' of having processed this stream
		//! This functionality is dependent on Elk being present:  since
		//! the meaning of what the 'summary' is differs from stream to
		//! stream, the result is simply an Elk object.  The default is
		//! to return the exit code contained in the footer.
		virtual Object get_summary() const { return False ; }

		//! \brief some debugging info
		virtual string type_name() const { return "Stream" ; }

	protected:
		// doesn't belong here, but it's convenient
		void read_next_message( google::protobuf::io::CodedInputStream&, const std::string& ) ;
} ;

typedef ::Holder<Stream> StreamHolder ;

int transfer( Stream& in, Stream& out ) ;

//! \brief base class of streams that read from many streams
class StreamBundle : public Stream
{
	protected:
		std::deque< StreamHolder > streams_ ;
		typedef std::deque< StreamHolder >::const_iterator citer ;
		typedef std::deque< StreamHolder >::const_reverse_iterator criter ;

	public:
		void add_stream( StreamHolder s ) { streams_.push_back( s ) ; }
} ;

struct ParseError : public Exception {
	std::string msg_ ;
	ParseError( const std::string& m ) : msg_(m) {}
	virtual void print_to( std::ostream& s ) const { s << "AnfoReader: " << msg_ ; }
} ;


//! \brief presents ANFO files as series of messages
//! This class will read a stream of results formatted as a continuous
//! protobuf message.
class AnfoReader : public Stream
{
	private:
		std::auto_ptr< google::protobuf::io::ZeroCopyInputStream > is_ ;
		std::string name_ ;

	public: 
		AnfoReader( std::auto_ptr< google::protobuf::io::ZeroCopyInputStream > is, const std::string& name ) ;
		virtual Result fetch_result() ;

		//! \internal
		virtual string type_name() const { return "AnfoReader(" + name_ + ")" ; }
} ;

// \brief reader for all supported formats
// Here we take care not to open files before the header is requested.
// This is necessary to allow merging thousands of files without
// directly running into a filedescriptor limit.
// To this end, the UniversalReader can be initialized with or without a
// stream object.  If the stream exists, we take care not to read from
// it until the header is needed, and the name given serves just for
// informational purposes.  If no stream exists, we create a
// FileInputStream from the name (which must be a filename, obviously)
// when the header is requested.  At this point we also inspect the
// stream to determine its format and create the appropriate filters to
// decode it.
class UniversalReader : public Stream
{
	private:
		std::auto_ptr< google::protobuf::io::ZeroCopyInputStream > is_ ;
		std::string name_ ;
		StreamHolder str_ ;

		bool solexa_scores_ ;
		int origin_ ;
		string genome_ ;

	public: 
		UniversalReader(
				const std::string& name,
				google::protobuf::io::ZeroCopyInputStream* is = 0,
				bool solexa_scores = false,
				int origin = 33,
				const string& genome = ""
				) ;

		virtual state get_state() { return str_ ? str_->get_state() : invalid ; }
		virtual Header fetch_header() ;
		virtual Result fetch_result() { if( get_state() == have_output ) return str_->fetch_result() ; throw "calling sequence violated" ; }
		virtual Footer fetch_footer() { return str_->fetch_footer() ; }
		virtual string type_name() const { return "UniversalReader(" + name_ + ")" ; }
} ;


//! \brief new blocked native format
//! Writes in a format that can be read by stream::ChunkedReader.  The
//! file is made up of individually compressed blocks so that
//! near-random access is possible... in principle.
class ChunkedWriter : public Stream
{
	public:
		// sensible buffer size: big enough to make compression worthwhile,
		// small enough that BZip won't split it again
		enum { default_buffer_size = 850000 } ;

		// supported compression methods
		enum method { none, fastlz, gzip, bzip } ;

	private:
		std::auto_ptr< google::protobuf::io::ZeroCopyOutputStream > zos_ ;	// final output
		std::vector< char > buf_ ;											// in-memory buffer
		std::auto_ptr< google::protobuf::io::ArrayOutputStream > aos_ ;		// output to buffer
		Chan chan_ ;
		std::string name_ ;
		int64_t wrote_ ;
		uint8_t method_, level_ ;

		void flush_buffer( unsigned needed = 0 ) ;
		void init() ;

	public:
		static uint8_t method_of( int ) ;
		static uint8_t level_of( int ) ;

		ChunkedWriter( const pair< google::protobuf::io::ZeroCopyOutputStream*, string >&, int ) ;
		ChunkedWriter( int fd, int, const char* = "<pipe>" ) ;
		ChunkedWriter( const char* fname, int ) ;
		virtual ~ChunkedWriter() ;

		virtual void put_header( const Header& h ) ;
		virtual void put_result( const Result& r ) ;
		virtual void put_footer( const Footer& f ) ;
		virtual string type_name() const { return "ChunkedWriter(" + name_ + ")" ; }
} ;

class ChunkedReader : public Stream
{
	private:
		std::auto_ptr< google::protobuf::io::ZeroCopyInputStream > is_ ;
		std::vector< char > buf_ ;											// in-memory buffer
		std::auto_ptr< google::protobuf::io::ArrayInputStream > ais_ ;		// output to buffer
		std::string name_ ;

		bool get_next_chunk() ;

	public: 
		ChunkedReader( std::auto_ptr< google::protobuf::io::ZeroCopyInputStream > is, const std::string& name ) ;
		virtual Result fetch_result() ;

		//! \internal
		virtual string type_name() const { return "ChunkedReader(" + name_ + ")" ; }
} ;

//! \brief filters that drop or modify isolated records
//! Think "mapMaybe".
//! \todo Maybe some stats could be gathered into some sort of a result.

class Filter : public Stream
{
	public:
		virtual bool xform( Result& ) = 0 ;
		virtual void put_result( const Result& res ) { res_ = res ; if( xform( res_ ) ) state_ = have_output ; }
} ;

//! \brief filters that drop some alignments
//! Filtering only applies to hits to the specified genome(s), or to all
//! hits if no genomes are specified.  Other hits pass through.
//! \todo Maybe some stats could be gathered into some sort of a result.
class HitFilter : public Stream
{
	private:
		vector<string> gs_ ;

	public:
		HitFilter() {}
		HitFilter( const vector<string> &gs ) : gs_(gs) {}

		virtual void put_header( const Header& h ) {
			Stream::put_header( h ) ;
			hdr_.clear_is_sorted_by_coordinate() ;
			hdr_.clear_is_sorted_by_all_genomes() ;
		}

		virtual bool keep( const Hit& ) = 0 ;
		virtual void put_result( const Result& res ) ; 
} ;

namespace {
	template < typename C, typename V > bool contains( const C& c, const V& v )
	{ return find( c.begin(), c.end(), v ) != c.end() ; }
}

//! \brief deletes some alignments
//! A list of sequences and a list of genomes can be given.  An
//! alignment is deleted if a) the list of genomes is empty or the hit
//! genome is a member of that list and b) the list of sequences is
//! empty or the hit sequence is a member of that list.
class IgnoreHit : public HitFilter
{
	private:
		vector< string > ss_ ;

	public:
		IgnoreHit( const vector< string > &gs, const vector< string > &ss ) : HitFilter( gs ), ss_( ss ) {}
		virtual bool keep( const Hit& h ) { return !ss_.empty() && !contains( ss_, h.sequence() ) ; }
} ;

//! \brief deletes hits to uninteresting genomes
//! Hits to the specified genomes pass through, all others are dropped.
class OnlyGenome : public Filter
{
	private:
		vector< string > gs_ ;

	public:
		OnlyGenome( const vector< string > &gs ) : Filter(), gs_( gs ) {}
		virtual bool xform( Result& r ) ;
} ;

//! \brief stream that filters for a given score
//! Alignments that exceed the score are deleted, but the sequences are
//! kept.  If genomes are given, only alignments to that genome are
//! tested and removed, else all are considered individually and
//! selectively removed.  If no alignments remains, \c reason is set to
//! \c bad_alignment.  Most interesting downstream operations will
//! ignore sequences without alignments.  Scores are in Phred scale now,
//! experience tells that a typical butoff between good alignments and
//! junk is \f$ 7.5 \cdot ( \mbox{length} - 20 ) \f$

class ScoreFilter : public HitFilter
{
	private:
		double slope_ ;
		double intercept_ ;

	public:
		ScoreFilter( double s, double i, const vector<string> &gs ) : HitFilter(gs), slope_(s), intercept_(i) {}
		virtual bool keep( const Hit& ) ;
} ;

class TotalScoreFilter : public Filter
{
    private:
        double slope_ ;
        double intercept_ ;
        vector< string > gs_ ;

    public:
		TotalScoreFilter( double s, double i, const vector<string> &gs ) : slope_(s), intercept_(i), gs_(gs) {}
        virtual bool xform( Result& ) ;
} ;

        
//! \brief stream that filters for minimum mapping quality
//! All alignments of sequences where the difference to the next hit is
//! too low are deleted.  Filtering can be restricted to some genomes.
//! The sequence itself is always kept.
class MapqFilter : public HitFilter
{
	private:
		int minmapq_ ; ;

	public:
		MapqFilter( const vector<string> &gs, int q ) : HitFilter(gs), minmapq_(q) {}
		virtual bool keep( const Hit& ) ;
} ;

//! \brief filter for some average quality
class QualFilter : public Filter
{
	private:
		double minqual_ ; ;

	public:
		QualFilter( double q ) : minqual_(q) {}
		virtual bool xform( Result& ) ;
} ;

//! \brief stream that filters for minimum sequence length
//! All alignments of sequences that are too short are deleted, relying
//! on down stream filters to completely get rid of the sequences
//! themselves.
class LengthFilter : public Filter
{
	private:
		int minlength_, maxlength_ ;

	public:
		LengthFilter( int l = 0, int h = std::numeric_limits<int>::max() )
			: minlength_(l), maxlength_(h) {}
		virtual bool xform( Result& ) ;
} ;

class GcFilter : public Filter
{
	private:
		int mingc_, maxgc_ ;

	public:
		GcFilter( int l = 0, int h = 100 ) : mingc_(l), maxgc_(h) {}
		virtual bool xform( Result& ) ;
} ;

//! \brief a stream that removes sequences without a specific hit
//! This is intended to shrink a file by removing junk that didn't
//! align at all.
class RequireHit : public Filter
{
	private:
		vector< string > gs_, ss_ ;

	public:
		RequireHit( const vector<string> &gs, const vector<string> &ss ) : gs_(gs), ss_(ss) {}
		virtual bool xform( Result& ) ;
} ;

//! \brief a stream that removes sequences without a specific best hit
//! This is intended to shrink a file by removing junk that didn't
//! align to the expected genome.
class RequireBestHit : public Filter
{
	private:
		vector< string > gs_, ss_ ;

	public:
		RequireBestHit( const vector<string> &gs, const vector<string> &ss ) : gs_(gs), ss_(ss) {}
		virtual bool xform( Result& ) ;
} ;

class Subsample : public Filter
{
	private:
		const float f_ ;

	public:
		Subsample( float f ) : f_(f) {}
		virtual bool xform( Result& ) ;
} ;
//! \brief filters for minimum multiplicity
//! Only results that stem from duplicate removal with a minimum
//! multiplicity are retained.  Intended for the analysis of libraries
//! sequenced with high redundancy to lower sequencing error.
class MultiFilter : public Filter
{
	private:
		unsigned n_ ;

	public:
		MultiFilter( unsigned n ) : n_(n) {}
		virtual bool xform( Result& r ) { return r.members_size() + r.nmembers() >= n_ ; }
} ;

//! \brief stream that filters out low quality bases
//! Bases with insufficient quality are replaced by N symbols.
//! (Originally I used gap symbols, which doesn't make sense and
//! actually confused the legacy tool downstream.  Ns should be fine and
//! are actually the correct symbol in a certain sense.)
//! This suppresses the counting of low quality bases in whatever
//! follows downstream.
class QualMasker : public Filter
{
	private:
		int q_ ;

	public:
		QualMasker( int q ) : q_(q) {}
		virtual bool xform( Result& ) ;
} ;

/*! \brief stream with PCR duplicates removed

    A set of duplicates is (more or less by definition) a set of reads
    that map to the same position.
   
	Any set of duplicates is merged (retaining the original reads in an
	auxilliary structure), and a consensus is called with new quality
	scores.  A quality score (for base A) is defined as \f[ Q := -10
	\log_{10} P( \bar{A} | \omega ), \f] that is, a representation of
	the probability of the base being wrong given some observation \f$
	\omega. \f$  It follows that \f[ Q = -10 \log_{10} \frac{
	P(\omega|\bar{A}) P(\bar{A}) }{ P(\omega) } = -10 \log_{10} P( \omega
	| \bar{A} ) \f] by assuming a uniform base composition in the sample
	(\f$ P(A)=\frac{1}{4} \f$) and after base calling (\f$
	P(\omega)=\frac{1}{4} \f$).

	When combining observations \f$ \omega_1, \ldots, \omega_n \f$
	conditional on the same base, they are all independent, hence the
	probabilities conditioned on the actual base can be multiplied.  To
	get the final quality when calling an A:
	\f[ P(A|\bigcap_n \omega_n) = \frac{ P(\bigcap_n \omega_n | A) P(A) }{ P(\bigcap_n \omega_n) }
	                            = \frac{ P(A) \prod_n P(\omega_n | A) }{ P(\bigcap_n \omega_n) } \f]

	The unconditional probabilities in the denominator are not
	independent (for they collectively depend on the true base), we have
	to replace them by total probabilites:
	\f[ = \frac{ P(A) \prod_n P(\omega_n | A) }{ \sum_N P(\bigcap_n \omega_n | N) P(N) } \f]

	Again, assuming a uniform base composition, the priors are all equal
	and cancel, and now the conditional probabilities are again
	independent.
	\f[ = \frac{ \prod_n P(\omega_n | A) }{ \sum_N \prod_n P(\omega_n | N) } \f]

	Tracking this incrementally is easy, we need to store the four
	products.  Computation needs to be done in the log domain to avoid
	loss of precision.  Likewise, computation of a quality must not involve \f$ P(A|\omega) \f$,
	since that is too close to one to be representable, so we calculate the quality as
	\f[ Q := -10 \log_{10} P( \bar{A} | \omega ) = -10 \log_{10} ( P(C | \omega ) + P(G | \omega ) + P(T | \omega ) ) \f]

	What remains is how to get an estimate of \f$
	P(\omega | N) \f$.  We easily get \f$ P(\omega | A) \f$ (see above),
	which will normally be very close to one, and the sum of the other
	three.  To estimate the summands, we set
	\f[ P(\omega|N) = \frac{P(N|\omega) * P(\omega)}{P(N)} =
	P(N|\omega) = P(N|\bar{A}) * P(\bar{A}|\omega) \f]
	and estimate \f$ P(N|\bar{A}) \f$ from misclassfication statistics
	of the base caller.
 */

class RmdupStream : public Stream
{
	private:
		output::Result cur_ ;
		std::vector< Logdom > quals_[4] ;
		double slope_ ;
		double intercept_ ;
		vector<string> gs_ ;
		int maxq_ ;
		bool have_foot_ ;

		// XXX double err_prob_[4][4] ; // get this from config or
		// something?

		bool is_duplicate( const Result& , const Result& ) const ;
		void add_read( const Result& ) ;
		void call_consensus() ;
		bool good_score( const Hit* ) const ;

	public:
		//! \brief sets parameters
		//! \param s Slope of score function, bad alignments are disregarded.
		//! \param i Intercept of score function.
		//! \param q (Assumed) quality of the polymerase.
		RmdupStream( double s, double i, int q ) :
			slope_(s), intercept_(i), maxq_( std::min(q,127) ), have_foot_(false) {}

		virtual void put_header( const Header& h )
		{
			if( !h.has_is_sorted_by_all_genomes() && !h.is_sorted_by_coordinate_size() )
				throw "RmdupStream: need sorted stream to remove duplicates" ;
			gs_.assign( h.is_sorted_by_coordinate().begin(), h.is_sorted_by_coordinate().end() ) ;
			Stream::put_header( h ) ;
		}

		virtual void put_result( const Result& ) ;
		virtual void put_footer( const Footer& ) ;
		virtual Result fetch_result() ;
} ;

//! \brief a stream that concatenates its input streams
//! The headers and footers are merged sensibly (plain merge with removal of
//! redundant information), then the results are simply concatenated.
class ConcatStream : public StreamBundle
{
	public:
		virtual Header fetch_header() ;
		virtual Result fetch_result() ;
} ;

class FastqReader : public Stream
{
	private:
		std::auto_ptr< google::protobuf::io::ZeroCopyInputStream > is_ ;
		bool sol_scores_ ;
		char origin_ ;

		void read_next_message() {
            if( read_fastq( is_.get(), *res_.mutable_read(), sol_scores_, origin_ ) ) {
                state_ = have_output ;
                sanitize( res_ ) ;
            } else {
                state_ = end_of_stream ;
				is_.reset( 0 ) ;
            }
		}

	public: 
		FastqReader( std::auto_ptr< google::protobuf::io::ZeroCopyInputStream > is, bool solexa_scores, char origin ) :
			is_( is ), sol_scores_(solexa_scores), origin_(origin) { read_next_message() ; }

		virtual Result fetch_result() { Result r ; std::swap( r, res_ ) ; read_next_message() ; return r ; }
		virtual string type_name() const { return "FastqReader" ; }
} ;

class SamReader : public Stream
{
	private:
		std::auto_ptr< google::protobuf::io::ZeroCopyInputStream > is_ ;
		string name_, genome_ ;
		Chan progress_ ;
		int nmsg_ ;

		void read_next_message() {
            if( read_sam( is_.get(), genome_, res_ ) ) {
                state_ = have_output ;
                sanitize( res_ ) ;
				if( (++nmsg_ & 0xffff) == 0 ) {
					stringstream ss ;
					ss << name_ << ": " << nmsg_ << " records" ;
					progress_( Console::info, ss.str() ) ;
				}
            } else {
                state_ = end_of_stream ;
				is_.reset( 0 ) ;
				progress_.close() ;
				foot_.set_exit_code( 0 ) ;
            }
		}

	public: 
		SamReader( std::auto_ptr< google::protobuf::io::ZeroCopyInputStream > is, const string& n, const string& g ) :
			is_( is ), name_( n ), genome_( g ), nmsg_(0) { read_next_message() ; }

		virtual Header fetch_header() { hdr_.add_is_sorted_by_coordinate( "" ) ; return Stream::fetch_header() ; }
		virtual Result fetch_result() { Result r ; std::swap( r, res_ ) ; read_next_message() ; return r ; }
		virtual string type_name() const { return "SamReader" ; }
} ;

class SffReader : public Stream
{
	private:
		std::auto_ptr< google::protobuf::io::ZeroCopyInputStream > is_ ;
		string name_ ;
		unsigned remaining_ ;
		unsigned number_of_flows_ ;

		const void* buf_ ;
		int buf_size_ ;

		uint8_t read_uint8() ;
		uint16_t read_uint16() ;
		uint32_t read_uint32() ;
		void read_string( unsigned, string* ) ;
		void skip( int ) ;

	public:
		SffReader( auto_ptr< google::protobuf::io::ZeroCopyInputStream > is, const string& name ) : 
			is_( is ), name_( name ), buf_(0), buf_size_(0) {}

		virtual Header fetch_header() ;
		virtual Result fetch_result() ;
		virtual string type_name() const { return "SffReader(" + name_ + ")" ; }
} ;

class BamReader : public Stream
{
	private:
		std::auto_ptr< google::protobuf::io::ZeroCopyInputStream > is_ ;
		string name_ ;
		string genome_ ;
		vector< string > refseqs_ ;

		const void* buf_ ;
		int buf_size_ ;

		uint8_t read_uint8() ;
		uint16_t read_uint16() ;
		uint32_t read_uint32() ;
		void read_string( unsigned, string* ) ;

	public:
		BamReader( auto_ptr< google::protobuf::io::ZeroCopyInputStream > is, const string& name, const string& genome ) ;
		virtual Result fetch_result() ;
		virtual string type_name() const { return "BamReader(" + name_ + ")" ; }
} ;

} // namespace streams

class PipeInputStream : public google::protobuf::io::FileInputStream 
{
	private:
		pid_t cpid_ ;

	public:
		PipeInputStream( int fd, pid_t cpid ) : google::protobuf::io::FileInputStream( fd ), cpid_( cpid ) {} 
		virtual ~PipeInputStream() { Close() ; throw_errno_if_minus1( waitpid( cpid_, 0, 0 ), "waiting for pipe process" ) ; }
} ;

std::pair< PipeInputStream*, std::string > make_PipeInputStream( const std::string& ) ;

class PipeOutputStream : public google::protobuf::io::FileOutputStream 
{
	private:
		pid_t cpid_ ;

	public:
		PipeOutputStream( int fd, pid_t cpid ) : google::protobuf::io::FileOutputStream( fd ), cpid_( cpid ) {} 
		virtual ~PipeOutputStream() { Close() ; throw_errno_if_minus1( waitpid( cpid_, 0, 0 ), "waiting for pipe process" ) ; }
} ;

std::pair< PipeOutputStream*, std::string > make_PipeOutputStream( const std::string& ) ;

//! \brief adapts a ZeroCopyOutputStream into a streambuf
class zero_copy_output_buf : public std::streambuf {
	private:
		std::auto_ptr< google::protobuf::io::ZeroCopyOutputStream > os_ ;

	public:
		zero_copy_output_buf( google::protobuf::io::ZeroCopyOutputStream* os ) : os_(os) {}
		virtual ~zero_copy_output_buf() ;
		virtual int sync() ;
		virtual int_type overflow( int_type ) ;
} ;

class zero_copy_ostream : public std::ostream {
	private:
		zero_copy_output_buf b_ ;

	public:
		zero_copy_ostream( google::protobuf::io::ZeroCopyOutputStream* os ) 
			: std::ostream( &b_ ), b_( os ) {}
} ;

#endif