Decipher Potts encoding of sequence state/configuration.

diff --git a/Makefile b/Makefile
index 13ad714..98b38b5 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -19,7 +19,7 @@ LIBS = -lgsl -lgslcblas -lm
 # = `gsl-config --libs`
 
 ifeq ($(dbg),1)
-       DBGFLAG = -g -DDEBUG
+       DBGFLAG = -g -DDEBUG  # -D_GLIBCXX_DEBUG
 else
        DBGFLAG = -O3 -combine
 endif

diff --git a/seqTools.cpp b/seqTools.cpp
index 3a831de..44c9f0d 100644 (file)
--- a/seqTools.cpp
+++ b/seqTools.cpp
@@ -1,7 +1,14 @@
 #include <iostream>
+#include <vector>
+#include <array>
+#include <algorithm>
 #include <unordered_map>
 
-enum class aa {unknown, A,R,N,D,C, Q,E,G,H,I, L,K,M,F,P, S,T,W,Y,V, B,Z,X, star, gap,};
+typedef std::vector<std::vector<int> > IntVector;
+typedef std::vector<std::vector<char> > CharVector;
+
+//enum class aa {BEGIN, unknown=BEGIN, A,R,N,D,C, Q,E,G,H,I, L,K,M,F,P, S,T,W,Y,V, B,Z,X, star, gap, END};
+enum aa {aa_BEGIN, aa_unknown=aa_BEGIN, aa_A,aa_R,aa_N,aa_D,aa_C, aa_Q,aa_E,aa_G,aa_H,aa_I, aa_L,aa_K,aa_M,aa_F,aa_P, aa_S,aa_T,aa_W,aa_Y,aa_V, aa_B,aa_Z,aa_X, aa_star, aa_gap, aa_END};
 enum class nt {T,C,A,G};
 enum class codon {TTT,TTC,TTA,TTG, TCT,TCC,TCA,TCG, TAT,TAC,TAA,TAG, TGT,TGC,TGA,TGG,
                   CTT,CTC,CTA,CTG, CCT,CCC,CCA,CCG, CAT,CAC,CAA,CAG, CGT,CGC,CGA,CGG,
@@ -31,22 +38,22 @@ namespace std {
 }
 
 codon2aa_map codon2aa {
-    {codon::TTT, aa::F}, {codon::TTC, aa::F}, {codon::TTA, aa::L}, {codon::TTG, aa::L},
-    {codon::TCT, aa::S}, {codon::TCC, aa::S}, {codon::TCA, aa::S}, {codon::TCG, aa::S},
-    {codon::TAT, aa::Y}, {codon::TAC, aa::Y}, {codon::TAA, aa::star}, {codon::TAG, aa::star},
-    {codon::TGT, aa::C}, {codon::TGC, aa::C}, {codon::TGA, aa::star}, {codon::TGG, aa::W},
-    {codon::CTT, aa::L}, {codon::CTC, aa::L}, {codon::CTA, aa::L}, {codon::CTG, aa::L},
-    {codon::CCT, aa::P}, {codon::CCC, aa::P}, {codon::CCA, aa::P}, {codon::CCG, aa::P},
-    {codon::CAT, aa::H}, {codon::CAC, aa::H}, {codon::CAA, aa::Q}, {codon::CAG, aa::Q},
-    {codon::CGT, aa::R}, {codon::CGC, aa::R}, {codon::CGA, aa::R}, {codon::CGG, aa::R},
-    {codon::ATT, aa::I}, {codon::ATC, aa::I}, {codon::ATA, aa::I}, {codon::ATG, aa::M},
-    {codon::ACT, aa::T}, {codon::ACC, aa::T}, {codon::ACA, aa::T}, {codon::ACG, aa::T},
-    {codon::AAT, aa::N}, {codon::AAC, aa::N}, {codon::AAA, aa::K}, {codon::AAG, aa::K},
-    {codon::AGT, aa::S}, {codon::AGC, aa::S}, {codon::AGA, aa::R}, {codon::AGG, aa::R},
-    {codon::GTT, aa::V}, {codon::GTC, aa::V}, {codon::GTA, aa::V}, {codon::GTG, aa::V},
-    {codon::GCT, aa::A}, {codon::GCC, aa::A}, {codon::GCA, aa::A}, {codon::GCG, aa::A},
-    {codon::GAT, aa::D}, {codon::GAC, aa::D}, {codon::GAA, aa::E}, {codon::GAG, aa::E},
-    {codon::GGT, aa::G}, {codon::GGC, aa::G}, {codon::GGA, aa::G}, {codon::GGG, aa::G},
+    {codon::TTT, aa_F}, {codon::TTC, aa_F}, {codon::TTA, aa_L}, {codon::TTG, aa_L},
+    {codon::TCT, aa_S}, {codon::TCC, aa_S}, {codon::TCA, aa_S}, {codon::TCG, aa_S},
+    {codon::TAT, aa_Y}, {codon::TAC, aa_Y}, {codon::TAA, aa_star}, {codon::TAG, aa_star},
+    {codon::TGT, aa_C}, {codon::TGC, aa_C}, {codon::TGA, aa_star}, {codon::TGG, aa_W},
+    {codon::CTT, aa_L}, {codon::CTC, aa_L}, {codon::CTA, aa_L}, {codon::CTG, aa_L},
+    {codon::CCT, aa_P}, {codon::CCC, aa_P}, {codon::CCA, aa_P}, {codon::CCG, aa_P},
+    {codon::CAT, aa_H}, {codon::CAC, aa_H}, {codon::CAA, aa_Q}, {codon::CAG, aa_Q},
+    {codon::CGT, aa_R}, {codon::CGC, aa_R}, {codon::CGA, aa_R}, {codon::CGG, aa_R},
+    {codon::ATT, aa_I}, {codon::ATC, aa_I}, {codon::ATA, aa_I}, {codon::ATG, aa_M},
+    {codon::ACT, aa_T}, {codon::ACC, aa_T}, {codon::ACA, aa_T}, {codon::ACG, aa_T},
+    {codon::AAT, aa_N}, {codon::AAC, aa_N}, {codon::AAA, aa_K}, {codon::AAG, aa_K},
+    {codon::AGT, aa_S}, {codon::AGC, aa_S}, {codon::AGA, aa_R}, {codon::AGG, aa_R},
+    {codon::GTT, aa_V}, {codon::GTC, aa_V}, {codon::GTA, aa_V}, {codon::GTG, aa_V},
+    {codon::GCT, aa_A}, {codon::GCC, aa_A}, {codon::GCA, aa_A}, {codon::GCG, aa_A},
+    {codon::GAT, aa_D}, {codon::GAC, aa_D}, {codon::GAA, aa_E}, {codon::GAG, aa_E},
+    {codon::GGT, aa_G}, {codon::GGC, aa_G}, {codon::GGA, aa_G}, {codon::GGG, aa_G},
 };
 
 str2int_map aa2int {
@@ -66,19 +73,19 @@ int2str_map int2aa = {
 };
 
 aa2str_map aa2str = {
-    {aa::A,'A'}, {aa::R,'R'}, {aa::N,'N'}, {aa::D,'D'}, {aa::C,'C'},
-    {aa::Q,'Q'}, {aa::E,'E'}, {aa::G,'G'}, {aa::H,'H'}, {aa::I,'I'},
-    {aa::L,'L'}, {aa::K,'K'}, {aa::M,'M'}, {aa::F,'F'}, {aa::P,'P'},
-    {aa::S,'S'}, {aa::T,'T'}, {aa::W,'W'}, {aa::Y,'Y'}, {aa::V,'V'},
-    {aa::B,'B'}, {aa::Z,'Z'}, {aa::X,'X'}, {aa::star,'*'}, {aa::gap,'-'}, {aa::unknown,'?'}
+    {aa_A,'A'}, {aa_R,'R'}, {aa_N,'N'}, {aa_D,'D'}, {aa_C,'C'},
+    {aa_Q,'Q'}, {aa_E,'E'}, {aa_G,'G'}, {aa_H,'H'}, {aa_I,'I'},
+    {aa_L,'L'}, {aa_K,'K'}, {aa_M,'M'}, {aa_F,'F'}, {aa_P,'P'},
+    {aa_S,'S'}, {aa_T,'T'}, {aa_W,'W'}, {aa_Y,'Y'}, {aa_V,'V'},
+    {aa_B,'B'}, {aa_Z,'Z'}, {aa_X,'X'}, {aa_star,'*'}, {aa_gap,'-'}, {aa_unknown,'?'}
 };
 
 str2aa_map str2aa = {
-    {'A',aa::A}, {'R',aa::R}, {'N',aa::N}, {'D',aa::D}, {'C',aa::C},
-    {'Q',aa::Q}, {'E',aa::E}, {'G',aa::G}, {'H',aa::H}, {'I',aa::I},
-    {'L',aa::L}, {'K',aa::K}, {'M',aa::M}, {'F',aa::F}, {'P',aa::P},
-    {'S',aa::S}, {'T',aa::T}, {'W',aa::W}, {'Y',aa::Y}, {'V',aa::V},
-    {'B',aa::B}, {'Z',aa::Z}, {'X',aa::X}, {'*',aa::star}, {'-',aa::gap}, {'?',aa::unknown}
+    {'A',aa_A}, {'R',aa_R}, {'N',aa_N}, {'D',aa_D}, {'C',aa_C},
+    {'Q',aa_Q}, {'E',aa_E}, {'G',aa_G}, {'H',aa_H}, {'I',aa_I},
+    {'L',aa_L}, {'K',aa_K}, {'M',aa_M}, {'F',aa_F}, {'P',aa_P},
+    {'S',aa_S}, {'T',aa_T}, {'W',aa_W}, {'Y',aa_Y}, {'V',aa_V},
+    {'B',aa_B}, {'Z',aa_Z}, {'X',aa_X}, {'*',aa_star}, {'-',aa_gap}, {'?',aa_unknown}
 };
 
 codon2str_map codon2str = {
@@ -137,12 +144,97 @@ to x2y(from x)
 template <typename from, typename to>
 to x2y(from x) {return static_cast<to>(x);}
 
+
+void getSeq2State(FILE *input, CharVector &p)
+{
+    char c;
+    char o;
+    
+    while (fscanf(input,"%c",&c)==1) {
+    
+        p.push_back(std::vector<char>());
+        (p.back()).push_back(c);
+        
+        while (fscanf(input,"%c",&o)==1) {
+    
+            if (o=='\n' || o=='\r') break;
+            
+            fscanf(input,"%c",&c);
+            (p.back()).push_back(c);
+            
+        }
+        
+    }
+       
+}
+
 int main(int argc, char *argv[])
 {
+    FILE *f = fopen(argv[1],"r");
+    if (f==NULL) { perror("could not load file"); exit(1); }
+
+    CharVector Schar;
+    getSeq2State(f,Schar);
+    fclose(f);
+
+    for (int i=0; i<Schar.size(); i++) {
+        for (int j=0; j<Schar[i].size(); j++) {
+            std::cout << Schar[i][j] << ' ';
+        }
+        std::cout << '\n';
+    }
+
+    //IntVector aa2state;
+    std::vector<std::array<int,aa_END> > aa2state;
+    for (int i=0; i<Schar.size(); i++) {
+        aa2state.push_back(std::array<int,aa_END>());
+        // expect all lines to contain the 'X' character to represent all the other amino acids
+        std::vector<char>::iterator begin = Schar[i].begin(), end = Schar[i].end();
+        int Xidx = std::distance(begin, std::find(begin,end,'X'));
+        for (int aa=aa_BEGIN; aa != aa_END; ++aa) {
+            std::vector<char>::iterator it = std::find(begin,end,int2aa[aa]);
+            int idx;
+            if (it != end) {  // amino acid explicitly encoded
+                idx = std::distance(begin,it);
+            }
+            else {      // amino acid represented by 'X'
+                idx = Xidx;
+            }
+            // commonest amino acid goes to end
+            // maybe easier to subtract 1 and mod?
+            idx = idx==0 ? Schar[i].size()-1 : idx-1;
+            aa2state.back()[aa] = idx;
+        }
+    }
+
+    for (int i = 0; i < aa2state.size(); i++) {
+        for (int aa=aa_BEGIN; aa != aa_END; ++aa) {
+            std::cout << aa2state[i][aa] << ' ';
+        }
+        std::cout << '\n';
+    }
+
+    std::vector<char> configChar = {'A','Y'};
+    std::vector<int> config(configChar.size());
+    for (int i=0; i < config.size(); i++) {
+        std::cout << configChar[i] << ' ';
+        config[i] = aa2int[configChar[i]];
+    }
+    std::cout << '\n';
+
+    for (int i=0; i < config.size(); i++) {
+        std::cout << aa2state[i][config[i]] << ' ';
+    }
+    std::cout << '\n';
+
+
+
+    /*
     std::cout << x2y<int,float>(5) << std::endl;
     std::cout << (int)(codon2aa[codon::TTT]) << std::endl;
     std::cout << static_cast<int>(codon2aa[codon::TTC]) << std::endl;
     std::cout << aa2str[codon2aa[codon::TTC]] << std::endl;
+    */
     //std::cout << aa2str[str2aa[int2aa[aa2int[codon2aa[str2codon["TTC"]]]]]] << std::endl;
 
 #if 0

diff --git a/snip.cpp b/snip.cpp
index dde3cfd..a39157e 100644 (file)
--- a/snip.cpp
+++ b/snip.cpp
@@ -1,4 +1,6 @@
 #include <vector>\r
+#include <iostream>\r
+#include <assert.h>\r
 \r
 #include "snip.h"\r
 \r
@@ -39,6 +41,10 @@ void getCouplings(FILE *input, Vector &p) {
 \r
 \r
 // Set energies based on an initial configuration\r
+// Note: config[i] = 0 denotes the 2nd most common state,\r
+// whereas config[i] = J[i].size() is the top state\r
+// Also, "gauge fixing" sets field and couplings between\r
+// the commonest states to zero.\r
 \r
 double computeEnergy(const Vector &J, const std::vector<int> &config) {\r
     \r
@@ -47,6 +53,8 @@ double computeEnergy(const Vector &J, const std::vector<int> &config) {
     double E = 0;\r
     \r
     for (int i=0;i<config.size();i++) {\r
+\r
+        assert(config[i] <= J[i].size());\r
         \r
         if (config[i]!=J[i].size()) {\r
         \r
@@ -66,3 +74,22 @@ double computeEnergy(const Vector &J, const std::vector<int> &config) {
     \r
 }\r
 \r
+int main(int argc, char **argv)\r
+{\r
+    FILE* f = fopen(argv[1],"r");\r
+    Vector J;\r
+    getCouplings(f,J);\r
+    fclose(f);\r
+\r
+    for (int i=0; i<J.size(); i++) {\r
+        for (int j=0; j<J[i].size(); j++) {\r
+            std::cout << J[i][j] << ' ';\r
+        }\r
+        std::cout << '\n';\r
+    }\r
+\r
+    std::cout << computeEnergy(J,{3,2}) << '\n';\r
+\r
+    return 0;\r
+}\r
+\r