libMems/Aligner.h

Go to the documentation of this file.

/*******************************************************************************
 * $Id: Aligner.h,v 1.23 2004/04/19 23:10:13 darling Exp $
 * This file is copyright 2002-2007 Aaron Darling and authors listed in the AUTHORS file.
 * This file is licensed under the GPL.
 * Please see the file called COPYING for licensing details.
 * **************
 ******************************************************************************/

#ifndef _Aligner_h_
#define _Aligner_h_

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include "libGenome/gnSequence.h"
#include "libMems/DNAMemorySML.h"
#include "libMems/GappedAligner.h"
#include "libMems/MatchList.h"
#include "libMems/Interval.h"
#include "libMems/IntervalList.h"
#include "libMems/MemHash.h"
#include "libMems/MaskedMemHash.h"
#include <map>
#include "libMems/NumericMatrix.h"
#include "libMems/GreedyBreakpointElimination.h"
#include <list>
#include "libMems/LCB.h"
#include "libMUSCLE/threadstorage.h"

namespace mems {

class LabeledMem{
public:
        Match* mem;
        uint label;
};

class LabeledMemComparator {
public:
        LabeledMemComparator( uint seq ){
                m_seq = seq;
        }
        LabeledMemComparator( LabeledMemComparator& lmc ){
                m_seq = lmc.m_seq;
        }
        boolean operator()(const LabeledMem& a, const LabeledMem& b) const{
                
                int64 a_start = a.mem->Start( m_seq ), b_start = b.mem->Start( m_seq );
                if( a_start == NO_MATCH || b_start == NO_MATCH ){
                        if( b_start != NO_MATCH )
                                return true;
                        return false;
                }
                if(a_start < 0)
                        a_start = -a_start;
//                      a_start = -a_start + a.mem->Length();
                if(b_start < 0)
                        b_start = -b_start;
//                      b_start = -b_start + b.mem->Length();
                int64 diff = a_start - b_start;
                return diff < 0;
        }
        ~LabeledMemComparator() {
        }
protected:
        uint m_seq;
private:
        LabeledMemComparator();
};

class PlacementMatch{
public:
        Match* mem;
        std::list< LabeledMem >::iterator iter;
};

class PlacementMatchComparator {
public:
        PlacementMatchComparator( uint seq ){
                m_seq = seq;
        }
        PlacementMatchComparator( PlacementMatchComparator& lmc ){
                m_seq = lmc.m_seq;
        }
        boolean operator()(const PlacementMatch& a, const PlacementMatch& b) const{
                
                int64 a_start = a.mem->Start( m_seq ), b_start = b.mem->Start( m_seq );
                if( a_start == NO_MATCH || b_start == NO_MATCH ){
                        if( b_start != NO_MATCH )
                                return true;
                        return false;
                }
                if(a_start < 0)
                        a_start = -a_start;
//                      a_start = -a_start + a.mem->Length();
                if(b_start < 0)
                        b_start = -b_start;
//                      b_start = -b_start + b.mem->Length();

                int64 diff = a_start - b_start;
                return diff < 0;
        }
        ~PlacementMatchComparator() {
        }
protected:
        uint m_seq;
private:
        PlacementMatchComparator();
};


typedef std::pair< mems::Match*, mems::Match* > search_cache_t;


class Aligner {
public:
        Aligner( uint seq_count );
        Aligner( const Aligner& al );
        Aligner& operator=( const Aligner& al );

        void align( MatchList& mlist, IntervalList& interval_list, double LCB_minimum_density, double LCB_minimum_range, boolean recursive, boolean extend_lcbs, boolean gapped_alignment, std::string tree_filename = "" );
        
        void Recursion( MatchList& r_list, Match* r_begin, Match* r_end, boolean nway_only = false );
        void GetBestLCB( MatchList& r_list, MatchList& best_lcb );
        void DoSomethingCool( MatchList& mlist, Interval& iv );
        
        void SetMinRecursionGapLength( gnSeqI min_r_gap );

        void SetMaxExtensionIterations( uint ext_iters ){ this->max_extension_iters = ext_iters; }

        void SearchWithinLCB( MatchList& mlist, std::vector< search_cache_t >& new_cache, bool leftmost = false, bool rightmost = false );
        void RecursiveAnchorSearch( MatchList& mlist, gnSeqI minimum_weight, std::vector< MatchList >& LCB_list, boolean entire_genome, std::ostream* status_out = NULL );

        void AlignLCB( MatchList& mlist, Interval& iv );
        void SetGappedAligner( GappedAligner& gal );
        void SetMaxGappedAlignmentLength( gnSeqI len );

        void SetPermutationOutput( std::string& permutation_filename, int64 permutation_weight );
        void WritePermutation( std::vector< LCB >& adjacencies, std::string out_filename );

        void SetRecursive( bool value ){ this->recursive = value; }
protected:
        TLS<MemHash> gap_mh;                    
        MaskedMemHash nway_mh;  
        uint32 seq_count;               
        boolean debug;                  
        double LCB_minimum_density;
        double LCB_minimum_range;

        uint max_extension_iters;       
        int64 cur_min_coverage; 
        gnSeqI min_recursive_gap_length;        
        void consistencyCheck( uint lcb_count, std::vector< LCB >& adjacencies, std::vector< MatchList >& lcb_list, std::vector< int64 >& weights );
        
        boolean recursive;              
        boolean extend_lcbs;    
        boolean gapped_alignment;       
        boolean currently_recursing;    
        boolean collinear_genomes;      
        GappedAligner* gal;

        std::string permutation_filename;
        int64 permutation_weight;

        std::vector< search_cache_t > search_cache; 
};

CREATE_EXCEPTION( AlignerError );

void transposeMatches( MatchList& mlist, uint seqI, const std::vector< int64 >& seq_regions );

void EliminateOverlaps( MatchList& ml );

void AaronsLCB( MatchList& mlist, std::set<uint>& breakpoints );


void ComputeLCBs( MatchList& meml, std::set<uint>& breakpoints, std::vector<MatchList>& lcb_list, std::vector<int64>& weights );
void computeLCBAdjacencies_v2( std::vector<MatchList>& lcb_list, std::vector< int64 >& weights, std::vector< LCB >& adjacencies );
void computeLCBAdjacencies_v2( IntervalList& iv_list, std::vector< int64 >& weights, std::vector< LCB >& adjacencies );
void scanLeft( int& left_recurseI, std::vector< LCB >& adjacencies, int min_weight, int seqI );
void scanRight( int& right_recurseI, std::vector< LCB >& adjacencies, int min_weight, int seqI );
void GetLCBCoverage( MatchList& lcb, uint64& coverage );

int64 greedyBreakpointElimination( gnSeqI minimum_weight, std::vector< LCB >& adjacencies, std::vector< int64 >& weights, std::ostream* status_out = NULL );
void filterMatches( std::vector< LCB >& adjacencies, std::vector< MatchList >& lcb_list, std::vector< int64 >& weights );

void CreateGapSearchList( std::vector< LCB >& adjacencies, const std::vector< genome::gnSequence* >& seq_table, std::vector< std::vector< int64 > >& iv_regions, boolean entire_genome );
void SearchLCBGaps( MatchList& new_matches, const std::vector< std::vector< int64 > >& iv_regions, MaskedMemHash& nway_mh );

static const uint MIN_ANCHOR_LENGTH = 9;


class SearchCacheComparator
{
public:
        SearchCacheComparator() : msc(0){};
        bool operator()( const search_cache_t& a, const search_cache_t& b ) const
        {
                bool lt = true;
                if( a.first == NULL )
                {
                        if( b.first == NULL )
                                lt = false;
                }else if( b.first == NULL )
                {
                        lt = false;
                }else if( !msc( a.first, b.first ) )
                        lt = false;
                else if( a.second == NULL )
                {
                        if( b.second == NULL )
                                lt = false;
                }else if( b.second == NULL )
                {
                        lt = false;
                }else if( !msc( a.second, b.second ) )
                        lt = false;

                return lt;              
        }
protected:
        mems::MatchStartComparator<mems::Match> msc;
};

static SearchCacheComparator cache_comparator;


}

#endif // _Aligner_h_

---