Changeset 138


Ignore:
Timestamp:
11/08/09 20:30:15 (10 months ago)
Author:
reed
Message:

all the logic associated with const_align appears to be working

Location:
current/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • current/src/align.cpp

    r133 r138  
    4444double aligner::align(alignment &aln) 
    4545{ 
    46         double d; 
    47         aln.data.clear();        
    48         if(aln.seqA.dna.size() >= aln.seqB.dna.size()) 
    49                 d = align_x(aln.seqA.dna, aln.seqB.dna, aln.data); 
    50         else 
    51         { 
    52                 d = align_x(aln.seqB.dna, aln.seqA.dna, aln.data); 
    53                 for(aln_data::iterator it = aln.data.begin(); it != aln.data.end(); ++it) 
    54                         *it = -*it; 
    55         } 
    56         return d; 
     46        aln.data.clear(); 
     47        return align_x(aln.seqA.dna, aln.seqB.dna, aln.data); 
    5748} 
    5849 

  • current/src/align.h

    r133 r138  
    4040                 
    4141        template<class OS> 
    42         inline void print(OS &os, const char *msg=NULL) const; 
     42        inline void print(OS &os, const char *msg=NULL, int dir=0, bool swapped=false) const; 
    4343         
    4444protected: 
    … …  
    129129 
    130130template<class OS> 
    131 inline void alignment::print(OS &os, const char *msg) const 
     131inline void alignment::print(OS &os, const char *msg, int dir, bool swapped) const 
    132132{ 
    133         std::string strA, strB, strC; 
     133        std::string strA, strB, strC, 
     134                nameA(seqA.name.substr(0, 14)), 
     135                nameB(seqB.name.substr(0, 14)); 
     136        //preallocate 
     137        strA.reserve(4048); 
     138        strB.reserve(4048); 
     139        strC.reserve(4048); 
     140 
    134141        std::string::const_iterator ait = seqA.dna.begin(), bit = seqB.dna.begin(); 
    135142        for(aln_data::const_iterator cit = data.begin(); cit != data.end(); ++cit) 
    … …  
    154161                } 
    155162        } 
    156                          
     163        if(swapped) { 
     164                swap(strA, strB); 
     165                swap(nameA, nameB); 
     166        } 
     167        seq_transform(strA, dir); 
     168        seq_transform(strB, dir); 
     169        seq_transform(strC, dir); 
     170 
    157171        os << "CLUSTAL multiple sequence alignment (Created by " << PACKAGE_STRING; 
    158172        if(msg != NULL) 
    … …  
    169183                ss = strA.substr(u, 60); 
    170184                a += ss.length() - std::count(ss.begin(), ss.end(), chGap); 
    171                 os << std::setw(14) << std::setiosflags(std::ios::left) << seqA.name.substr(0, 14) 
     185                os << std::setw(14) << std::setiosflags(std::ios::left) << nameA 
    172186                        << " " << ss << " " << a << std::endl; 
    173187                 
    174188                ss = strB.substr(u, 60); 
    175189                b += ss.length() - std::count(ss.begin(), ss.end(), chGap); 
    176                 os << std::setw(14) << std::setiosflags(std::ios::left) << seqA.name.substr(0, 14) 
     190                os << std::setw(14) << std::setiosflags(std::ios::left) << nameB 
    177191                         << " " << ss << " " << b << std::endl; 
    178192                 
    … …  
    183197 
    184198#endif 
    185  

  • current/src/ngila.cpp

    r137 r138  
    126126        } 
    127127         
    128         mydb.unique_sort(seq_db::DIR_ORI | (arg.const_align & 7)); 
    129                                          
     128        int direction = 0; 
     129        if(arg.const_align != 0) { 
     130                direction = mydb.unique_sort(seq_db::DIR_ORI | (arg.const_align & 7)); 
     131                mydb.transform(direction); 
     132                cerr << direction << endl; 
     133                if(!(arg.const_align & 8)) { // return to original order 
     134                        mydb.rearrange(mydb.db().get<id>().begin()); 
     135                } 
     136        } 
     137 
    130138        cost_model *pmod = NULL; 
    131139        string model_keys[] = { string("zeta"), string("geo"), string("cost") }; 
    … …  
    182190                if(cit != pvec.begin()) 
    183191                        cout << "//" << endl; 
    184                 alignment aln(mydb[cit->first], mydb[cit->second]); 
     192                // swap a and b so that a's hashed position is lower 
     193                size_t a = cit->first; 
     194                size_t b = cit->second; 
     195                bool swapped = false; 
     196                if(!(arg.const_align & 8) && mydb.db().project<hashid>(mydb.db().begin()+a) >  
     197                        mydb.db().project<hashid>(mydb.db().begin()+b) ) { 
     198                        swap(a,b); 
     199                        swapped = true; 
     200                } 
     201 
     202                alignment aln(mydb[a], mydb[b]); 
    185203                double dcost = alner.align(aln); 
    186                         dcost += pmod->offset(mydb[cit->first].dna, 
    187                                       mydb[cit->second].dna); 
     204                        dcost += pmod->offset(mydb[a].dna, 
     205                                      mydb[b].dna); 
    188206                ostringstream msg; 
    189207                msg << "Cost = " << setprecision(10) << dcost; 
    190208                 
    191                 aln.print(cout, msg.str().c_str()); 
     209                aln.print(cout, msg.str().c_str(), ((arg.const_align & 16) ? 0 : direction), swapped); 
    192210        } 
    193211         

  • current/src/seqdb.cpp

    r137 r138  
    115115                                complementer(sd.dna.rend())), DIR_RVC, sd)); 
    116116        } 
    117         foreach(const hash_data &h, hc) { 
    118                 cerr << h.ref.get().name << "\t" << h.dir << "\t" << h.hashed << endl; 
    119         } 
    120117        size_t hh = hc.begin()->hashed; 
    121118        int dd = hc.begin()->dir, mask = -1;  
    … …  
    142139        static const int ties[] = { 8,1,2,1,4,1,2,1,8,8,8,8,8,8,8,8 }; 
    143140        dd = ties[dd&15]; 
    144         //sort 
     141        // sort 
    145142        cont.rearrange(hc.get<1>().find(dd)); 
     143        // store this order 
     144        cont.get<hashid>().rearrange(cont.begin()); 
    146145        return dd; 
    147146} 
    148  

  • current/src/seqdb.h

    r137 r138  
    3535#include <boost/functional/hash.hpp> 
    3636 
     37 
     38inline char complement(char x) { 
     39        static char dnac[] = "\0\x01\x02\x03\x04\x05\x06\a\b\t\n\x0b\f\r\x0e\x0f" 
     40        "\x10\x11\x12\x13\x14\x15\x16\x17\x18\x19\x1a\x1b\x1c\x1d\x1e\x1f" 
     41        " !\"#$%&'()*+,-./0123456789:;<=>?@" 
     42        "TVGHEFCDIJMLKNOPQYSAABWXRZ" "[\\]^_`" 
     43        "tvghefcdijmlknopqysaabwxrz" "{|}~\x7f"; 
     44 
     45        return dnac[x]; 
     46} 
     47 
    3748template<class _Iterator> 
    3849class _complementer : public _Iterator { 
    … …  
    4253 
    4354        value_type operator*() const { 
    44                 static char dnac[] = "\0\x01\x02\x03\x04\x05\x06\a\b\t\n\x0b\f\r\x0e\x0f" 
    45                         "\x10\x11\x12\x13\x14\x15\x16\x17\x18\x19\x1a\x1b\x1c\x1d\x1e\x1f" 
    46                         " !\"#$%&'()*+,-./0123456789:;<=>?@" 
    47                         "TVGHEFCDIJMLKNOPQYSAABWXRZ" "[\\]^_`" 
    48                         "tvghefcdijmlknopqysaabwxrz" "{|}~\x7f"; 
    49                 return value_type(dnac[**static_cast<const base_type*>(this)]); 
     55                return value_type(complement(**static_cast<const base_type*>(this))); 
    5056        } 
    5157         
    … …  
    5763        return _complementer<_Iterator>(it); 
    5864} 
     65 
     66class seq_db; 
    5967 
    6068struct seq_data { 
    … …  
    7886                        boost::bind(gf, g, _1, 0) != std::string::npos), 
    7987                        dna.end()); 
     88                return *this; 
    8089        } 
    8190         
    … …  
    8594                        std::transform(dna.begin(), dna.end(), dna.begin(), std::ptr_fun(::toupper)); 
    8695        } 
     96 
     97        inline void transform(int dir); 
    8798}; 
    8899 
    … …  
    95106public: 
    96107        typedef boost::multi_index_container< seq_data, boost::multi_index::indexed_by< 
    97                 boost::multi_index::random_access<>, 
    98                 boost::multi_index::random_access<boost::multi_index::tag<id> >, 
     108                boost::multi_index::random_access<>, // current order 
     109                boost::multi_index::random_access<boost::multi_index::tag<id> >, // original order 
     110                boost::multi_index::random_access<boost::multi_index::tag<hashid> >, // hashed order 
    99111                boost::multi_index::ordered_unique< 
    100112                        boost::multi_index::tag<name>, 
    … …  
    137149                cont.rearrange(it); 
    138150        } 
    139          
    140         int unique_sort(int directions=DIR_ORI); 
    141          
    142 //      void create_hashes(int directions=seq_data::DIR_ORI) { 
    143 //              for(container::iterator it = cont.begin(); it != cont.end(); ++it) 
    144 //                      cont.modify(it, bind(&seq_data::create_hash, _1, directions)); 
    145 //      } 
     151                 
     152        void transform(int dir) { 
     153                for(container::iterator it = cont.begin(); it != cont.end(); ++it) 
     154                        cont.modify(it, bind(&seq_data::transform, _1, dir)); 
     155        } 
    146156         
    147157        bool parse_file(const char *csfile, bool bappend=false, bool bi=false); 
     158        int unique_sort(int directions=DIR_ORI); 
    148159         
    149160        seq_db() : ss_gaps("-") { } 
    … …  
    156167}; 
    157168 
     169inline void seq_transform(std::string& dna, int dir) { 
     170        switch(dir&7) { 
     171        default: 
     172        case 0: 
     173                break; 
     174        case seq_db::DIR_COM: 
     175                std::transform(dna.begin(), dna.end(), dna.begin(), complement); 
     176                break; 
     177        case seq_db::DIR_RVC: 
     178        case (seq_db::DIR_REV|seq_db::DIR_COM): 
     179                std::transform(dna.begin(), dna.end(), dna.begin(), complement); 
     180        case seq_db::DIR_REV: 
     181                std::reverse(dna.begin(), dna.end()); 
     182                break;           
     183        }; 
     184} 
     185 
     186inline void seq_data::transform(int dir) { 
     187        seq_transform(dna, dir); 
     188} 
    158189 
    159190#endif 
Note: See TracChangeset for help on using the changeset viewer.