Changeset 151

Timestamp:

02/19/10 19:08:58 (7 months ago)

Author:

reed

Message:

added several distance matrix output formats; removed identiy out in aln and fasta formats

Location:

current/src

Files:

• align.cpp (modified) (3 diffs)
• align.h (modified) (6 diffs)
• ngila.cpp (modified) (6 diffs)
• sort.h (modified) (1 prop)
• xm.h (modified) (1 prop)

current/src/align.cpp

@@ -27 +27 @@
 using namespace std;
 
-#ifdef min
-#       undef min
-#endif
-
 template<class T>
 inline T min3(T a, T b, T c)
 {
-        return min(a, min(b,c));
+        return (std::min)(a, (std::min)(b,c));
 }
 
@@ -67 +63 @@
         
         // Allocate travel table to maximum possible size
-        size_t szA = std::min(szMa,seqA.size())+1;
+        size_t szA = (std::min)(szMa,seqA.size())+1;
         tabTravel.resize(szA);
-        size_t szB = std::min(szMb,seqB.size())+1;
+        size_t szB = (std::min)(szMb,seqB.size())+1;
         for(travel_table::iterator it = tabTravel.begin();
             it != tabTravel.end(); ++it)
@@ -582 +578 @@
 }
 
+double alignment::identity() const {
+        aln_data::size_type matches = 0, mismatches = 0;
+        std::string::const_iterator ait = seqA.dna.begin(), bit = seqB.dna.begin();
+        for(aln_data::const_iterator cit = data.begin(); cit != data.end(); ++cit) {
+                if(*cit == 0) {
+                        if(*ait == *bit)
+                                ++matches;
+                        else
+                                ++mismatches;
+                        ++ait;
+                        ++bit;
+                } else if(*cit > 0)
+                        bit+=*cit;
+                else
+                        ait-=*cit;
+        }
+        return (1.0*matches)/(matches+mismatches);
+}
+

current/src/align.h

@@ -42 +42 @@
         inline void print(OS &os, int format=0, double cost=0.0, int dir=0, bool swapped=false) const;
         
+        double identity() const;
+        
 protected:
         const seq_data &seqA, &seqB;
@@ -50 +52 @@
 };
 
-class aligner
-{
+class aligner {
 public:
-        struct indel
-        {
+        struct indel {
                 indel() { }
                 indel(size_t pp, size_t xx, double dd)
@@ -62 +62 @@
                 double d; // Score
         };
-        struct min_mid_cost
-        {
+        struct min_mid_cost {
                 double c; // Minimum cost
                 size_t s; // minimum position in s-vector
@@ -141 +140 @@
 
         std::string::const_iterator ait = seqA.dna.begin(), bit = seqB.dna.begin();
-        unsigned int matches = 0, mismatches = 0;
         for(aln_data::const_iterator cit = data.begin(); cit != data.end(); ++cit) {
                 if(*cit == 0) {
-                        if(*ait == *bit) {
+                        if(*ait == *bit)
                                 strC.append(1, '*');
-                                ++matches;
-                        } else {
+                        else
                                 strC.append(1, ' ');
-                                ++mismatches;
-                        }
                         strA.append(ait, ait+1); ++ait;
                         strB.append(bit, bit+1); ++bit;
@@ -173 +168 @@
         size_t sz = strA.size();        
         std::string ss;
-        double identity = double(matches)/(matches+mismatches);
         if(format == 0) {
                 os << "CLUSTAL multiple sequence alignment (Created by " << PACKAGE_STRING
-                   << "; Cost = " << std::setprecision(10) << cost  << " Identity = " << identity
+                   << "; Cost = " << std::setprecision(10) << cost
                    << ")" << std::endl << std::endl;
 
@@ -199 +193 @@
                 os << ">" << nameA << " " 
                    << "(Created by " << PACKAGE_STRING
-                   << "; Cost = " << std::setprecision(10) << cost  << " Identity = " << identity
+                   << "; Cost = " << std::setprecision(10) << cost
                    << ")" << std::endl;
                 for(size_t u = 0; u < sz; u += 80)

current/src/ngila.cpp

@@ -20 +20 @@
 #include <iomanip>
 #include <iostream>
+#include <limits>
 
 #include <boost/preprocessor.hpp>
 #include <boost/foreach.hpp>
 #include <boost/config.hpp>
+//BOOST_DISABLE_ASSERTS
 #include <boost/multi_array.hpp>
+#include <boost/algorithm/string/predicate.hpp>
 
 #include "ngila_app.h"
@@ -96 +99 @@
 }
 
-template<size_t _N>
-size_t key_switch(const std::string &ss, const std::string (&key)[_N]) {
-        for(size_t i=0;i<_N;++i) {
-                if(key[i].find(ss) == 0)
+template<std::size_t _N>
+std::size_t key_switch(const std::string &ss, const std::string (&key)[_N]) {
+        using boost::algorithm::starts_with;
+        for(std::size_t i=0;i<_N;++i) {
+                if(starts_with(key[i], ss))
                         return i;
         }
-        return (size_t)-1;
+        return (std::size_t)-1;
 }
 
@@ -169 +173 @@
         string pairs_keys[] = { string("first"), string("all"), string("each") };
         seq_db::size_type table_size = mydb.size();
-        switch(key_switch(arg.pairs, pairs_keys))
-        {
+        
+        switch(key_switch(arg.pairs, pairs_keys)) {
         case 0:
                 pvec.push_back(make_pair(0,1));
@@ -189 +193 @@
         };
         
-        string format_keys[] = { string("aln"), string("fasta"), string("smat") /*,string("imat")*/ };
+        string format_keys[] = {
+                string("aln"), string("fasta"),
+                string("dist"), string("dist-c"),
+                string("dist-d"), string("dist-i")
+        };
         int out_format = 0;
         if(!arg.output.empty()) {
@@ -223 +231 @@
                 }
         }
-        
-        typedef boost::multi_array<float, 2> tabmat;
-        tabmat cost_table(boost::extents[table_size][table_size]);
-
-        for(pair_vec::const_iterator cit = pvec.begin(); cit != pvec.end(); ++cit)
-        {
-                if(cit != pvec.begin())
-                        cout << "//" << endl;
+        ostream &myout = fout.is_open() ? static_cast<ostream&>(fout)
+                                        : static_cast<ostream&>(cout);
+        typedef boost::multi_array<float, 2> dist_mat;
+        dist_mat dist_table;
+        const bool do_dist = (out_format >= 2);
+        if(do_dist) {
+                dist_table.resize(boost::extents[table_size][table_size]);
+                fill(dist_table.data(), dist_table.data()+dist_table.num_elements(), numeric_limits<float>::quiet_NaN());
+                float diag = (out_format == 5) ? 1.0f : 0.0f;
+                for(seq_db::size_type i = 0; i < table_size; ++i)
+                        dist_table[i][i] = diag;
+        }
+        
+        for(pair_vec::const_iterator cit = pvec.begin(); cit != pvec.end(); ++cit) {
+                if(!do_dist && cit != pvec.begin())
+                        myout << "//" << endl;
                 // swap a and b so that a's hashed position is lower
                 size_t a = cit->first;
@@ -243 +259 @@
                 alignment aln(mydb[a], mydb[b]);
                 double dcost = alner.align(aln);
-                        dcost += pmod->offset(mydb[a].dna,
-                                      mydb[b].dna);
+                dcost += pmod->offset(mydb[a].dna, mydb[b].dna);
+                double dident = aln.identity();
                 
-                if(fout.is_open()) {
-                        aln.print(fout, out_format, dcost,
+                if(!do_dist) {
+                        aln.print(myout, out_format, dcost,
                                 ((arg.const_align & 16) ? 0 : direction), swapped);
                 } else {
-                        aln.print(cout, out_format, dcost,
-                                ((arg.const_align & 16) ? 0 : direction), swapped);
+                        double upper,lower;
+                        if(a > b)
+                                swap(a,b);
+                        switch(out_format) {
+                        default:
+                        case 2: //dist
+                                lower = 1.0-dident;
+                                upper = dcost;
+                                break;
+                        case 3: //dist-c
+                                upper = lower = dcost;
+                                break;
+                        case 4: //dist-d
+                                upper = lower = 1.0-dident;
+                                break;
+                        case 5: //dist-i
+                                upper = lower = dident;
+                                break;
+                        }
+                        dist_table[a][b] = (float)upper;
+                        dist_table[b][a] = (float)lower;
+                }
+        }
+        if(do_dist) {
+                myout << setprecision(10);
+                for(seq_db::size_type i = 0; i < table_size; ++i) {
+                        for(seq_db::size_type j = 0; j < table_size-1; ++j)
+                                myout << dist_table[i][j] << "\t";
+                        myout << dist_table[i][table_size-1] << endl;
                 }
         }