Import Potts epitope (amino acid sequence). Helper split() function on arbitrary...

diff --git a/pop_ss.cpp b/pop_ss.cpp
index dce29dd..2556ead 100644 (file)
--- a/pop_ss.cpp
+++ b/pop_ss.cpp
@@ -140,6 +140,16 @@ AARecognizer::AARecognizer(const std::string &seq, int start_index)
     }
 }
 
+AARecognizer::AARecognizer(const std::string &seq, const std::vector<unsigned> &s)
+: site(s), len(seq.length())
+{
+    ref.reserve(len);
+    for (size_t i=0; i<len; ++i) {
+        ref.push_back(str2aa.at(seq[i]));
+    }
+
+}
+
 // Did virus escape from immune pressure by mutating *any* of the sites away from ref?
 bool AARecognizer::escaped(const NTVirus &v) const {
     return AARecognizer::escaped(v.aa_seq);

diff --git a/pop_ss.h b/pop_ss.h
index 20273aa..221276f 100644 (file)
--- a/pop_ss.h
+++ b/pop_ss.h
@@ -85,6 +85,7 @@ public:
     std::vector<aa> ref;
     size_t len;
     AARecognizer(const std::string &seq, int start_index);
+    AARecognizer(const std::string &seq, const std::vector<unsigned> &s);
     bool escaped(const NTVirus &v) const;
     bool escaped(const std::vector<aa> &aa_seq) const;
 };

diff --git a/ss.cpp b/ss.cpp
index b3ea4ad..153b09c 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -386,15 +386,6 @@ void importState_Potts(RunParameters_SS &r)
 }
 
 
-
-// load epitope definitions from file
-
-void importEpitope_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, NTVirusSpecies *V, Species_parray &print_spec)
-{
-
-}
-
-
 // load epitope definitions from file
 
 void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V, Species_parray &print_spec)
@@ -554,6 +545,84 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea
 }
 
 
+// load epitope definitions from file
+
+void importEpitope_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, NTVirusSpecies *V, Species_parray &print_spec)
+{
+    if (r.num_CTL_clones != r.epitopefiles.size()) {
+        std::cerr << "ERROR in importEpitope_Potts: invalid number of CTL clones" << std::endl;
+        exit(1);
+    }
+
+    for (size_t i = 0; i != r.num_CTL_clones; ++i) {
+        CTLaaSpecies* M = new CTLaaSpecies("M",r.init_CTL_numM[i]);
+
+        std::ifstream input(r.epitopefiles[i].c_str());   // <infile>.ep
+        if (!input) { perror((std::string("ERROR in importEpitope: ") + r.epitopeInfile).c_str()); exit(1); }
+
+        double aff;
+        std::string aa_str_seq;
+
+        std::string readStr;
+        while (std::getline(input,readStr)) {
+            std::stringstream readStrStrm(readStr);
+
+            readStrStrm >> aff;
+            readStrStrm >> aa_str_seq;
+
+            std::string siteMaybe;      // either a single integer, or an A-B string (with A and B integers) denoting "site[1]-site[len]", inclusive
+            readStrStrm >> siteMaybe;
+            unsigned site;
+
+            std::vector<std::string> sites = split(siteMaybe,'-');
+
+            if (sites.size() == 2) {    // "A-B" string, not integer
+                site = strtouint(sites.at(0));
+                if (strtouint(sites.at(1)) + 1 != site + aa_str_seq.length()) {     // site[len] in [site[1],site[2],...,site[len]] must have value site[1] + len - 1
+                    std::cerr << "ERROR: based on length of epitope (" << aa_str_seq << "), "
+                              << "expected site range to be " << site << "-" << site+aa_str_seq.length()-1
+                              << ", not " << siteMaybe
+                              << std::endl;
+                    exit(1);
+                }
+                M->epitopes.push_back(AARecognizer(aa_str_seq,site));
+            }
+            else {  // list of indices
+                std::vector<unsigned> s;
+                s.push_back(site);
+                while (readStrStrm >> site) {
+                    s.push_back(site);
+                }
+                if (s.size() != aa_str_seq.size()) {
+                    std::cerr << "ERROR: length of epitope (" << aa_str_seq << ") does not match number of given sites"
+                              << std::endl;
+                    exit(1);
+                }
+                M->epitopes.push_back(AARecognizer(aa_str_seq,s));
+            }
+
+            M->affinity.push_back(aff);
+            M->num_ep += 1;
+
+        }
+
+        for (size_t e=0; e<M->num_ep; e++) {    // XXX
+            std::cout << M->affinity.back() << '\t';
+            for (size_t i=0; i<M->epitopes[e].len; i++) {
+                std::cout << aa2str.at(M->epitopes[e].ref[i]);
+            }
+            std::cout << '\t';
+            for (size_t i=0; i<M->epitopes[e].len; i++) {
+                std::cout << M->epitopes[e].site[i] << ' ';
+            }
+            std::cout << '\n';
+        }
+
+
+    }
+}
+
+
 void usage()
 {
     std::cout <<

diff --git a/ss.h b/ss.h
index b52d6f5..eb0c85b 100644 (file)
--- a/ss.h
+++ b/ss.h
@@ -48,7 +48,7 @@ inline int strtoint(const std::string &s) {
 
 // Converts a string to an unsigned integer
 
-inline int strtouint(const std::string &s) {
+inline unsigned int strtouint(const std::string &s) {
     
     std::istringstream i(s);
     unsigned int x;
@@ -71,6 +71,24 @@ inline double strtodouble(const std::string &s) {
 
 }
 
+// from http://stackoverflow.com/questions/236129/split-a-string-in-c
+inline std::vector<std::string> &split(const std::string &s, char delim, std::vector<std::string> &elems) {
+    std::stringstream ss(s);
+    std::string item;
+    while (std::getline(ss, item, delim)) {
+        elems.push_back(item);
+    }
+    return elems;
+}
+
+// from http://stackoverflow.com/questions/236129/split-a-string-in-c
+inline std::vector<std::string> split(const std::string &s, char delim) {
+    std::vector<std::string> elems;
+    split(s, delim, elems);
+    return elems;
+}
+
+
 // print a set of mutations to stream as a CSV-like string
 std::ostream& operator<<(std::ostream& os, std::vector<unsigned int> s)
 {