Added fill_sequence parameter to SegmentModel.set_sequence and to the mrdna...

Added fill_sequence parameter to SegmentModel.set_sequence and to the mrdna and enrgmd scripts, with default value 'T', since that is most common for origami. Thanks to Elija Feigl for recommending the change.

bin/enrgmd

@@ -14,6 +14,16 @@ parser.add_argument('-o','--output-prefix', type=str, default=None,
                     help="Name for your job's output")
 parser.add_argument('-d','--directory', type=str,  default=None,
                     help='Directory for simulation; does not need to exist yet')
+
+parser.add_argument('--sequence-file', type=str, default=None,
+                    help='Sequence of longest strand.')
+
+parser.add_argument('--fill-sequence', choices=['A','T','C','G','random'], default='T',
+                    help='Fill sequence where otherwise unset.')
+
+parser.add_argument('--crossover-to-intrahelical-cutoff', type=float, default=-1,
+                    help='Set the distance (in Angstroms) beyond which crossovers between helix ends are converted to intrahelical connections; a negative value means no crossovers will be converted')
+
 parser.add_argument('--backbone-scale', type=float, default=1.0,
                     help='Factor to scale DNA backbone in atomic model; try 0.25 to avoid clashes for atomistic simulations')
 
@@ -55,6 +65,15 @@ if __name__ == '__main__':
                     
     model = read_model( str(infile) )
 
+    if args.crossover_to_intrahelical_cutoff >= 0:
+        model.convert_crossovers_at_ends_beyond_cutoff_to_intrahelical(
+            args.crossover_to_intrahelical_cutoff)
+
+    if args.sequence_file is not None:
+        from mrdna.model.dna_sequence import read_sequence_file
+        seq = read_sequence_file(args.sequence_file)
+        model.set_sequence( seq, args.fill_sequence )
+
     d_orig = os.getcwd()
     try:
         if directory is None:

bin/mrdna

@@ -36,6 +36,9 @@ parser.add_argument('--dimensions', type=str, default=None,
 parser.add_argument('--sequence-file', type=str, default=None,
                     help='Sequence of longest strand.')
 
+parser.add_argument('--fill-sequence', choices=['A','T','C','G','random'], default='T',
+                    help='Fill sequence where otherwise unset.')
+
 parser.add_argument('--output-period', type=float, default=1e4,
                     help='Simulation steps between DCD frames')
 # parser.add_argument('--minimization-steps', type=float, default=0,
@@ -132,7 +135,7 @@ def main():
     if args.sequence_file is not None:
         from mrdna.model.dna_sequence import read_sequence_file
         seq = read_sequence_file(args.sequence_file)
-        model.set_sequence( seq )
+        model.set_sequence( seq, args.fill_sequence )
 
     if args.output_prefix is not None:
         prefix = args.output_prefix

mrdna/resources/oxDNA1_sequence_dependent_parameters.txt

+STCK_FACT_EPS = 0.18f
+
+STCK_G_C =  1.684
+
+STCK_C_G =  1.737
+
+STCK_G_G =  1.604
+STCK_C_C =  1.604
+
+STCK_G_A =  1.590
+STCK_T_C =  1.590
+
+STCK_A_G =  1.610
+STCK_C_T =  1.610
+
+STCK_T_G =  1.654
+STCK_C_A =  1.654
+
+STCK_G_T =  1.671
+STCK_A_C =  1.671
+
+STCK_A_T =  1.553
+
+STCK_T_A =  1.634
+
+STCK_A_A =  1.709
+STCK_T_T =  1.709
+
+HYDR_A_T =  0.893
+HYDR_T_A =  0.893
+
+HYDR_C_G =  1.243
+HYDR_G_C =  1.243

mrdna/resources/oxDNA2_sequence_dependent_parameters.txt

+STCK_FACT_EPS = 0.18
+STCK_G_C = 1.69339
+STCK_C_G = 1.74669
+STCK_G_G = 1.61295
+STCK_C_C = 1.61295
+STCK_G_A = 1.59887
+STCK_T_C = 1.59887
+STCK_A_G = 1.61898
+STCK_C_T = 1.61898
+STCK_T_G = 1.66322
+STCK_C_A = 1.66322
+STCK_G_T = 1.68032
+STCK_A_C = 1.68032
+STCK_A_T = 1.56166
+STCK_T_A = 1.64311
+STCK_A_A = 1.84642
+STCK_T_T = 1.58952
+HYDR_A_T = 0.88537
+HYDR_T_A = 0.88537
+HYDR_C_G = 1.23238
+HYDR_G_C = 1.23238

mrdna/segmentmodel.py

@@ -772,6 +772,21 @@ class Segment(ConnectableElement, Group):
                 if self.sequence[i] is None:
                     self.sequence[i] = random.choice(bases)
 
+    def assign_unset_sequence(self, fill_sequence='T'):
+        if fill_sequence == 'random':
+            self.randomize_unset_sequence()
+        elif fill_sequence in 'A T C G'.split():
+            if self.sequence is None:
+                self.sequence = [fill_sequence]*self.num_nt
+            else:
+                assert(len(self.sequence) == self.num_nt) # TODO move
+                for i in range(len(self.sequence)):
+                    if self.sequence[i] is None:
+                        self.sequence[i] = 'T'
+        else:
+            raise ValueError("'fill_sequence' parameter must be one of 'random' 'A' 'T' 'C' or 'G'.")
+
+
     def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead ):
         raise NotImplementedError
 
@@ -3208,7 +3223,7 @@ class SegmentModel(ArbdModel):
             crossovers.append(tmp)
         return crossovers
 
-    def set_sequence(self, sequence, force=True):
+    def set_sequence(self, sequence, force=True, fill_sequence='random'):
         if force:
             self.strands[0].set_sequence(sequence)
         else:
@@ -3217,9 +3232,10 @@ class SegmentModel(ArbdModel):
             except:
                 ...
         for s in self.segments:
-            s.randomize_unset_sequence()
+            s.assign_unset_sequence(fill_sequence)
+
 
-    def _set_cadnano2_sequence(self, cadnano_sequence_csv_file, replace_unset=None):
+    def _set_cadnano2_sequence(self, cadnano_sequence_csv_file, sequence_fill='T'):
         try:
             self.segments[0]._cadnano_helix
         except:
@@ -3240,8 +3256,8 @@ class SegmentModel(ArbdModel):
                 if values[0] == "Start": continue
                 start,end,seq = values[:3]
                 strand = starts[start]
-                if replace_unset in base_set:
-                    seq = [replace_unset if s == '?' else s for s in seq]
+                if sequence_fill in base_set:
+                    seq = [sequence_fill if s == '?' else s for s in seq]
                 else:
                     seq = [random.choice(base_set) if s == '?' else s for s in seq]
                 if len([s for s in seq if s not in base_set]) > 0: