CMakeLists.txt

@@ -35,6 +35,7 @@ set(disable-multi-promoters "" CACHE STRING "Forbid that a gene is translated by
 set(dna-factory-alg "L2G" CACHE STRING "Which memory allocation algorithm to use for managing the DnaFactory pool of DNAs")
 
 set(with-triangle OFF CACHE BOOL "Whether to enable triangle phenotypic target (else Gaussian)")
+set(with-eukaryote ON CACHE BOOL "Whether to enable eukaryotic version")
 set(with-detectclone ON CACHE BOOL "Whether to enable clones and not recompute them")
 set(with-floatconcentration OFF CACHE BOOL "Whether to enable the encoding of concentration has float (and not double)")
 set(with-perf-traces OFF CACHE BOOL "Whether to activate performance traces of (R-)Aevol")
@@ -150,6 +151,10 @@ if ( ${with-triangle} )
     add_definitions(-DPHENOTYPIC_TARGET_TRIANGLE)
 endif ()
 
+if ( ${with-eukaryote} )
+    add_definitions(-D__EUKARYOTE)
+endif ()
+
 if ( ${with-detectclone} )
     add_definitions(-D__DETECT_CLONE)
 endif ()

src/libaevol/7/DnaMutator.cpp

@@ -165,9 +165,14 @@ MutationEvent* DnaMutator::generate_next_mutation(int32_t length) {
         if (pos_1 < pos_2) {
           seqlen = pos_2 - pos_1;
         } else {
-          int32_t tmp1_len = length_ - pos_1;
-          int32_t tmp2_len = pos_2;
-          seqlen = tmp1_len + tmp2_len;
+          #ifdef __EUKARYOTE
+            Utils::exchange(pos_1, pos_2);
+            seqlen = pos_2 - pos_1;
+          #else
+            int32_t tmp1_len = length_ - pos_1;
+            int32_t tmp2_len = pos_2;
+            seqlen = tmp1_len + tmp2_len;
+          #endif
         }
 
       int32_t genome_size_after = length_ + seqlen;
@@ -206,8 +211,13 @@ MutationEvent* DnaMutator::generate_next_mutation(int32_t length) {
       if (pos_1 < pos_2) {
         genome_size_after = length_ - (pos_2 - pos_1);
       } else {
+        #ifdef __EUKARYOTE
+        Utils::exchange(pos_1, pos_2);
+        genome_size_after = length_ - (pos_2 - pos_1);
+        #else
         genome_size_after = length_ - (length_ - pos_1);
         genome_size_after = genome_size_after - pos_2;
+        #endif
       }
       
       // printf("Remove %d %d : LB %d LA %d\n",pos_1,pos_2,length_,genome_size_after);
@@ -226,6 +236,9 @@ MutationEvent* DnaMutator::generate_next_mutation(int32_t length) {
 
     }
     else if (random_value < nb_large_dupl_ + nb_large_del_ + nb_large_trans_) {
+      #ifdef __EUKARYOTE
+      std::cout << "Warning : performing a translocation with a eukaryotic individual"<<std::endl;
+      #endif
       nb_large_trans_--;
 
         if (length_ == 1)
@@ -292,6 +305,7 @@ MutationEvent* DnaMutator::generate_next_mutation(int32_t length) {
       // if (pos_1 == pos_2) return nullptr; // Invert everything <=> Invert nothing!
 
       // Invert the segment that don't contain OriC
+      // Convenient for prokaryotes, absolutely necessary for eukaryotes
       if (pos_1 > pos_2) Utils::exchange(pos_1, pos_2);
 
       mevent = new MutationEvent();
@@ -374,6 +388,12 @@ MutationEvent* DnaMutator::generate_next_mutation(int32_t length) {
 
       if (length_ - nb_pos_to_del < min_genome_length_)
         return nullptr;
+      
+      #ifdef __EUKARYOTE
+      if (pos + nb_pos_to_del > length_){
+        nb_pos_to_del = length_ - pos; // Doing the maximal possible deletion
+      }
+      #endif
 
       length_ = length_ - nb_pos_to_del;
 

src/libaevol/7/Dna_7.cpp

@@ -34,6 +34,8 @@
 #include "legacy/ExpManager.h"
 #include "ExpManager_7.h"
 
+#include <cassert>
+
 #define REDUCTION_FACTOR 16
 
 namespace aevol {
@@ -227,11 +229,19 @@ Mutation* Dna_7::do_switch(size_type pos
 #endif
 
   // Remove promoters containing the switched base
-  ann_chrsm_->promoter_list().remove_promoters_around(pos, Utils::mod(pos + 1, length_), length_);
+  #ifdef __EUKARYOTE
+    ann_chrsm_->promoter_list().remove_promoters_around(pos, pos + 1);
+  #else
+    ann_chrsm_->promoter_list().remove_promoters_around(pos, Utils::mod(pos + 1, length_), length_);
+  #endif
 
   // Look for potential new promoters containing the switched base
   if (length() >= PROM_SIZE)
+  #ifdef __EUKARYOTE
+    ann_chrsm_->promoter_list().look_for_new_promoters_around(*this, pos, pos + 1);
+  #else
     ann_chrsm_->promoter_list().look_for_new_promoters_around(*this, pos, Utils::mod(pos + 1, length()));
+  #endif
 
   #ifdef BASE_2
     return new PointMutation(pos);
@@ -241,6 +251,13 @@ Mutation* Dna_7::do_switch(size_type pos
 }
 
 Mutation* Dna_7::do_small_insertion(size_type pos, int16_t nb_insert, char* seq) {
+  size_type end_look;
+  #ifdef __EUKARYOTE
+    end_look = std::max(length_, pos + nb_insert);
+  #else
+    end_look = Utils::mod(pos + nb_insert, length_);
+  #endif
+
   // Remove the promoters that will be broken
   ann_chrsm_->promoter_list().remove_promoters_around(pos, length_);
 
@@ -258,7 +275,7 @@ Mutation* Dna_7::do_small_insertion(size_type pos, int16_t nb_insert, char* seq)
     }
     else {
       ann_chrsm_->promoter_list().move_all_promoters_after(pos, nb_insert, length_);
-      ann_chrsm_->promoter_list().look_for_new_promoters_around(*this, pos, Utils::mod(pos + nb_insert, length_));
+      ann_chrsm_->promoter_list().look_for_new_promoters_around(*this, pos, end_look);
     }
   }
 
@@ -269,7 +286,14 @@ Mutation* Dna_7::do_small_deletion(size_type pos, int16_t nb_del) {
   auto old_pos = pos;
 
   // Remove promoters containing at least one nucleotide from the sequence to  delete
-    ann_chrsm_->promoter_list().remove_promoters_around(pos, Utils::mod(pos + nb_del, length_), length_);
+  auto end_del = size_type{0};
+  #ifdef __EUKARYOTE
+    end_del = std::min(length(), pos + nb_del);
+  #else
+    end_del = Utils::mod(pos + nb_del, length());
+  #endif
+
+  ann_chrsm_->promoter_list().remove_promoters_around(pos, end_del, length());
 
   // Do the deletion and update promoter list
   if (pos + nb_del <= length()) { // the deletion does not contain the origin of
@@ -285,6 +309,10 @@ Mutation* Dna_7::do_small_deletion(size_type pos, int16_t nb_del) {
     }
   }
   else { // the deletion contains the origin of replication
+  #ifdef __EUKARYOTE
+    printf("Deleting across position 0 : not permitted in linear chromosomes");
+    assert(false);
+  #endif
     // Do the deletion
     int32_t nb_del_at_pos_0 = nb_del - length() + pos;
 
@@ -335,6 +363,12 @@ Mutation* Dna_7::do_duplication(size_type pos_1, size_type pos_2, size_type pos_
 
   }
   else { // if (pos_1 >= pos_2)
+    #ifdef __EUKARYOTE
+      std::cout << "*** Invalid duplicated segment for linear chromosomes *** \n";
+      std::cout << " pos1 : " << pos_1 << " and pos2: " << pos_2 << " should be inverted" << std::endl;
+      assert(false);
+    #endif
+
     // The segment to duplicate includes the origin of replication.
     // The copying process will be done in two steps.
     //
@@ -394,6 +428,11 @@ Mutation* Dna_7::do_duplication(size_type pos_1, size_type pos_2, size_type pos_
 }
 
 Mutation* Dna_7::do_translocation(size_type pos_1, size_type pos_2, size_type pos_3, size_type pos_4, bool invert) {
+  #ifdef __EUKARYOTE
+    printf("Code not tested with translocation in eukaryote. Are you sure ?");
+    assert(false);
+  #endif
+
   auto pos_min = std::min(pos_1, std::min(pos_2, std::min(pos_3, pos_4)));
 
   if (not invert) {
@@ -444,6 +483,10 @@ Mutation* Dna_7::do_inversion(size_type pos_1, size_type pos_2) {
   //                                             pos_2    <-'
   //
 
+  #ifdef __EUKARYOTE
+    assert(pos_1 <= pos_2);
+  #endif
+
         if (length_ == 1)
         {
             printf("*** genome of size 1 ; inversion not done *** \n");
@@ -490,6 +533,10 @@ Mutation* Dna_7::do_inversion(size_type pos_1, size_type pos_2) {
 
 
 Mutation* Dna_7::do_deletion(size_type pos_1, size_type pos_2) {
+  #ifdef __EUKARYOTE
+    assert(pos_1 <= pos_2);
+  #endif
+
   // Delete segment going from pos_1 (included) to pos_2 (excluded)
     if (length_ == 1)
     {

src/libaevol/7/Individual_7.cpp

@@ -142,6 +142,20 @@ void Individual_7::start_stop_RNA() {
         if (motif_id >= 26 && motif_id < 48) {
           // LAGGING
           int t_motif_id = motif_id - 26;
+          #ifdef __EUKARYOTE
+            if (dna_pos < PROM_SIZE) {
+              for (int i = 0 ; i < PROM_SIZE ; ++i) {
+                prom_dist_lagging[i] = 1;
+              }
+              motif_id = 47;
+            }
+            else {
+              prom_dist_lagging[t_motif_id] =
+                PROM_SEQ_LAG[t_motif_id] ==
+                annotated_chromosome_->dna_->data_[dna_pos - t_motif_id]
+                ? 0 : 1;
+            }
+          #else
           prom_dist_lagging[t_motif_id] =
               PROM_SEQ_LAG[t_motif_id] ==
                       annotated_chromosome_->dna_->data_[
@@ -150,8 +164,23 @@ void Individual_7::start_stop_RNA() {
                                                          t_motif_id :
                                                    dna_pos - t_motif_id]
                   ? 0 : 1;
+        #endif
         } else if (motif_id < 22) {
           // LEADING
+          #ifdef __EUKARYOTE
+            if (dna_pos > len - PROM_SIZE) {
+              for (int i = 0 ; i < PROM_SIZE ; ++i) {
+                prom_dist_leading[i] = 1;
+              }
+              motif_id = 21;
+            }
+            else {
+              prom_dist_leading[motif_id] =
+                PROM_SEQ_LEAD[motif_id] ==
+                annotated_chromosome_->dna_->data_[dna_pos + motif_id]
+                ? 0 : 1;
+            }
+          #else
           prom_dist_leading[motif_id] =
               PROM_SEQ_LEAD[motif_id] ==
                       annotated_chromosome_->dna_->data_[
@@ -162,6 +191,7 @@ void Individual_7::start_stop_RNA() {
                                                           motif_id]
                   ? 0
                   : 1;
+          #endif
         } else {
           // removed terminator search
         }
@@ -953,7 +983,11 @@ void Individual_7::translate_protein(const ExpSetup&, double w_max) {
 
         count_loop++;
         c_pos += CODON_SIZE;
-        c_pos = c_pos >= dna_length ? c_pos - dna_length : c_pos;
+        #ifdef __EUKARYOTE
+          assert(c_pos < dna_length);
+        #else
+          c_pos = c_pos >= dna_length ? c_pos - dna_length : c_pos;
+        #endif
       }
     } else { // LAGGING
       #ifdef BASE_2

src/libaevol/7/MutationEvent.cpp

@@ -26,6 +26,7 @@
 
 
 #include "MutationEvent.h"
+#include <cassert>
 
 namespace aevol {
 
@@ -55,6 +56,9 @@ void MutationEvent::small_deletion(int32_t pos, int32_t number) {
 }
 
 void MutationEvent::duplication(int32_t pos1, int32_t pos2, int32_t pos3) {
+  #ifdef __EUKARYOTE
+    assert(pos1 < pos2);
+  #endif
   type_ = MutationEventType::DUPLICATION;
   pos_1_ = pos1;
   pos_2_ = pos2;
@@ -63,6 +67,9 @@ void MutationEvent::duplication(int32_t pos1, int32_t pos2, int32_t pos3) {
 
 void MutationEvent::translocation(int32_t pos1, int32_t pos2, int32_t pos3,
                                     int32_t pos4, bool invert) {
+  #ifdef __EUKARYOTE
+    assert(pos1 < pos2);
+  #endif
   type_ = MutationEventType::TRANSLOCATION;
   pos_1_ = pos1;
   pos_2_ = pos2;
@@ -72,12 +79,18 @@ void MutationEvent::translocation(int32_t pos1, int32_t pos2, int32_t pos3,
 }
 
 void MutationEvent::inversion(int32_t pos1, int32_t pos2) {
+  #ifdef __EUKARYOTE
+    assert(pos1 < pos2);
+  #endif
   type_ = MutationEventType::INVERSION;
   pos_1_ = pos1;
   pos_2_ = pos2;
 }
 
 void MutationEvent::deletion(int32_t pos1, int32_t pos2) {
+  #ifdef __EUKARYOTE
+    assert(pos1 < pos2);
+  #endif
   type_ = MutationEventType::DELETION;
   pos_1_ = pos1;
   pos_2_ = pos2;

src/libaevol/biochemistry/GeneticUnit.cpp

@@ -828,9 +828,19 @@ void GeneticUnit::do_transcription() {
       rna->set_transcript_length(-1);
 
       int32_t i;
+      #ifdef __EUKARYOTE
+      for (i = (strand_id == Strand::LEADING ? transcript_start : 0) ;
+           i < (strand_id == Strand::LEADING ? genome_length : transcript_start);
+           ++i) {
+      #else
       for (i = 0; i < genome_length; ++i) {
+      #endif
         if (is_terminator(strand_id,
-                          transcript_start + (strand_id == Strand::LEADING ? i : -i))) {
+          #ifdef __EUKARYOTE
+            (strand_id == Strand::LEADING ? i : transcript_start-i))) {
+          #else
+            transcript_start + (strand_id == Strand::LEADING ? i : -i))) {
+          #endif
 
           // Found terminator => set transcript's length
           rna->set_transcript_length(i + 
@@ -1330,6 +1340,11 @@ bool GeneticUnit::is_promoter(Strand strand, int32_t pos, int8_t& dist) const {
 
   int8_t prom_dist[PROM_SIZE];
   if(strand == Strand::LEADING) {
+    #ifdef __EUKARYOTE
+      if (pos >= len - PROM_SIZE){
+        return false; // Too close to the border
+      }
+    #endif
 // #pragma omp parallel for
 // #pragma omp simd
       // if (indiv_->id() == 0)printf("Search for Promoter %d : ",pos);
@@ -1346,6 +1361,11 @@ bool GeneticUnit::is_promoter(Strand strand, int32_t pos, int8_t& dist) const {
       // if (indiv_->id() == 0)printf("\n");
   }
   else { // LAGGING
+    #ifdef __EUKARYOTE
+      if (pos < 0 + PROM_SIZE) {
+        return false; // Too close to the border
+      }
+    #endif
 // #pragma omp parallel for
 // #pragma omp simd
       for (int8_t i = 0; i < PROM_SIZE; i++) {
@@ -1374,6 +1394,11 @@ bool GeneticUnit::is_terminator(Strand strand, int32_t pos) const {
   const char* genome = dna_->data();
   int32_t len = dna_->length();
   if (strand == Strand::LEADING) {
+    #ifdef __EUKARYOTE
+      if (pos + TERM_SIZE >= len) {
+        return false;
+      }
+    #endif
     for (int8_t i = 0; i < TERM_STEM_SIZE; i++) {
       #ifdef BASE_2
       if (genome[Utils::mod(pos + i, len)] ==
@@ -1404,6 +1429,11 @@ bool GeneticUnit::is_terminator(Strand strand, int32_t pos) const {
   }
   else // (strand == LAGGING)
   {
+    #ifdef __EUKARYOTE
+      if (pos - TERM_SIZE < 0) {
+        return false;
+      }
+    #endif
     for (int8_t i = 0; i < TERM_STEM_SIZE; i++) {
       #ifdef BASE_2
       if (genome[Utils::mod(pos - i, len)] ==

test/gtest/MiniEukTest.cpp

+// ****************************************************************************
+//
+//          Aevol - An in silico experimental evolution platform
+//
+// ****************************************************************************
+//
+// Copyright: See the AUTHORS file provided with the package or <www.aevol.fr>
+// Web: http://www.aevol.fr/
+// E-mail: See <http://www.aevol.fr/contact/>
+// Original Authors : Guillaume Beslon, Carole Knibbe, David Parsons
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+//*****************************************************************************
+
+
+
+
+// =================================================================
+//                              Includes
+// =================================================================
+#include <inttypes.h>
+#include <cstring>
+
+#include <list>
+#include <vector>
+#include <memory>
+
+#include <gtest/gtest.h>
+
+#include "DnaFactory.h"
+#include "ExpManager_7.h"
+#include "ExpSetup.h"
+#include "fuzzy/FuzzyFactory_7.h"
+#include "Individual_7.h"
+#include "legacy/biochemistry/Rna.h"
+#include "legacy/biochemistry/Protein.h"
+#include "macros.h"
+#include "MutationParams.h"
+#include "Strand.h"
+
+using namespace aevol;
+
+//############################################################################
+//                                                                           #
+//                         Class IndividualTest                              #
+//                                                                           #
+//############################################################################
+class IndividualTest : public testing::Test
+{
+ protected:
+  virtual void SetUp(void);
+  virtual void TearDown(void);
+
+  // We have an version of each individual for each fuzzy set flavor
+  Individual_7* indiv1;
+};
+
+// ===========================================================================
+//                                 Public Methods
+// ===========================================================================
+void IndividualTest::SetUp(void)
+{
+  // Build ad-hoc genomes
+  // (and reverse to test the same things on the lagging strand.):
+  //
+  // indiv1: LEAD (prom-1 + AS + AG + AS + term)
+  //   size 22+5+28+6+11=72
+  // AS = Arbitrary Sequence
+  // AG = Arbitrary Gene
+  // Do not modify the sequences !
+
+  // Define a few arbitrary sequences
+  char as[5][10] = {
+    "0011",
+    "11101",
+    "110011",
+    "11000",
+    "000101"
+  };
+
+  const char* SHINE_DAL_SEQ = "011011";
+  // const char* SHINE_DAL_SEQ_LAG = "001001";
+  // Define an arbitrary gene
+  char gene[100];
+  sprintf(gene, "%s0011000100110110010001", SHINE_DAL_SEQ);
+
+  // char lag_gene[100];
+  // sprintf(lag_gene, "1100111011001001101110%s", SHINE_DAL_SEQ_LAG);
+
+  // Define an arbitrary terminator
+  char term[TERM_SIZE+1] = "01000001101";
+  // char bad_term_end[TERM_SIZE] = "0100000110";
+  // char bad_term_beg[TERM_SIZE] = "1000001101";
+
+  // Define a few arbitrary promoters
+  // char prom[2][23] = {
+  //   "0101010001110110010110", // dist from consensus: 2 => basal level: 0.6
+  //   "0101011001110010010010"  // dist from consensus: 1 => basal level: 0.8
+  // };
+
+  // char lag_prom[2][23] = {
+  //   "1001011001000111010101", // dist from consensus: 2 => basal level: 0.6
+  //   "1011011011000110010101"  // dist from consensus: 1 => basal level: 0.8
+  // };
+
+
+  char bad_prom[2][23] = {
+    "101010001110110010110", // removed first base : add 0 at end
+    "010101100111001001001"  // removed last base : add 0 at begin
+  };
+
+  //   char bad_lag_prom[2][23] = {
+  //   "001011001000111010101", // removed first base : add 1 at end
+  //   "101101101100011001010"  // removed last base : add 1 at begin
+  // };
+
+
+  // Initialize the experimental setup.
+  // These are needed in the GeneticUnit constructors.
+  ExpSetup* exp_s = new ExpSetup(nullptr);
+  MutationParams params_mut;
+
+
+  // Update fuzzy flavor and fuzzy factory // TODO<dpa> Still needed/useful?
+  FuzzyFlavor fuzzyFlavor = FuzzyFlavor::VECTOR;
+  exp_s->set_fuzzy_flavor(fuzzyFlavor);
+
+  DnaFactory* dna_factory = new DnaFactory(DnaFactory_Policy::LOCAL_GLOBAL_FIT, 48, 1000, 16);
+  FuzzyFactory_7* fuzzy_factory = new FuzzyFactory_7(fuzzyFlavor, 16, 300, 16);
+
+
+
+  // Build indiv1                         LEAD (prom-1 + AS + AG + AS + term)
+  // Construct a genome with these arbitrary sequences
+  char* genome = new char[100];
+  sprintf(genome, "%s%s%s%s%s0",bad_prom[0], as[1], gene, as[2], term);
+
+  indiv1 = new Individual_7(1.0, dna_factory, fuzzy_factory);
+  indiv1->annotated_chromosome_[A]->dna_ = dna_factory->get_dna(strlen(genome));
+  indiv1->annotated_chromosome_[A]->dna_->ann_chrsm_ = indiv1->annotated_chromosome_[A];
+  indiv1->annotated_chromosome_[A]->dna_->dna_factory_ = dna_factory;
+  indiv1->annotated_chromosome_[A]->dna_->set_indiv(genome, strlen(genome), indiv1);
+
+  indiv1->annotated_chromosome_[B]->dna_ = dna_factory->get_dna(strlen(genome));
+  indiv1->annotated_chromosome_[B]->dna_->ann_chrsm_ = indiv1->annotated_chromosome_[B];
+  indiv1->annotated_chromosome_[B]->dna_->dna_factory_ = dna_factory;
+  indiv1->annotated_chromosome_[B]->dna_->set_indiv(genome, strlen(genome), indiv1);
+  genome = nullptr;
+
+  // Do transcription and translation
+  indiv1->start_stop_RNA();
+  indiv1->prom_compute_RNA(*exp_s);
+  indiv1->start_protein();
+  indiv1->compute_protein();
+  indiv1->translate_protein(*exp_s, 1.0);
+
+}
+
+void IndividualTest::TearDown() {
+//  for (auto indiv : indivs1) {
+//    delete indiv;
+//  }
+}
+
+// For each version of each individual constructed with a different
+// fuzzy set implementation, we check that all values are correct.
+// We don't need to change the fuzzy set implementation in the
+// experimental setup again because it's only used at transcription and
+// translation time.
+
+TEST_F(IndividualTest, TestIndiv1) {
+  // Check that we have the right number of promoters, terminators etc
+  // and at the right positions
+  // "right" means those values we have computed by hand
+    for (auto chrsm: {A,B}){
+    // Check genome size
+    EXPECT_EQ(72, indiv1->annotated_chromosome_[chrsm]->dna_->length());
+    EXPECT_EQ(72, indiv1->annotated_chromosome_[chrsm]->length());
+
+    // Check Prom list
+    EXPECT_EQ(0, indiv1->annotated_chromosome_[chrsm]->promoter_list().promoter_count());
+
+    // Check RNA list
+    auto rna_list = indiv1->annotated_chromosome_[chrsm]->rnas();
+    EXPECT_EQ(0, rna_list.size());
+
+    // Check protein list
+    auto prot_list = indiv1->annotated_chromosome_[chrsm]->proteins();
+    EXPECT_EQ(0, prot_list.size());
+    }
+}
+
+
+
+// ===========================================================================
+//                                Protected Methods
+// ===========================================================================
+
+// ===========================================================================
+//                              Non inline accessors
+// ===========================================================================