origami_utils¶

Utility functions for DNA origamis designed in cadnano and simulated in oxDNA.

origami_utils.py includes the class: Origami

class origami_utils.Origami(system=False, cad2cuda_file=False, visibility=False)¶

Origami object.

Parameters:	system – base.System object. Can be obtained from readers.LorenzoReader.get_system cad2cuda_file – path to virt2nuc file. Can be obtained from cadnano_interface.py visibility – path to visibility file

base_to_vbase(vhelix, base)¶: Not used right now.

check_branch_migration_strict(vh1, vb1, vh2, vb2)¶

Check for branch migration using strict definition; returns true if there is branch migration, false otherwise.

Strict definition: at least one of the 4 junction base pairs is missing.

self.interaction_list must be filled using get_h_bond_list.

The two virtual helices and two virtual bases that define the junction’s position are given as the four arguments.

Parameters:	vh1 – first virtual helix index in cadnano vb1 – first virtual base index in cadnano vh2 – second virtual helix index in cadnano vb2 – second virtual base index in cadnano

compute_vh_midpoints()¶: Computes the base pair midpoints for each virtual helix just once to save time

corrug_angle(xpos, bbms1, bbms2, id1, id2, av_axis, ref_intra_in_plane)¶: Computes the corrugation angle as the dot product of the intra-strand vector at xpos with the reference vector, after both have been projected into the plane perpendicular to the average helix axis at the junction.

corrugation(hj, span=16)¶

Returns the corrugation a distance span around the junction hj. The corrugation pattern is quantified as the angle from the dot product between the average of the two intra-helix vectors at the junction and each of the other intra-helix vectors.

Run compute_vh_midpoints() first.

Parameters:	hj – should come from hjs[ii] where hjs=origami.get_holliday_junctions(). span – distance span around Holliday junction

get_1d_vecs(discard_unbonded=True)¶

Returns orthonormal vectors based on the orientation of the strand. The returned vectors are “invariant” co-ordinates x’, y’, z’. x’: along double helix; yy: average backbone-base vector for the 0th strand; z’: x’ cross yy; y’: z’ cross x’

Parameters:	discard_unbonded – boolean. If true, discard broken base pairs in the calculation. NB: it is not currently implemented in this function.

get_3d_twist(f_edge_vh, vhelix_extent=False, discard_unbonded=True)¶

Returns the 3d global twist (top-bottom and left-right). Run get_h_bond_list first.

Parameters:	f_edge_vh – path to vh_edge file (see the function parse_vh_edge) vhelix_extent – 2-element array including the beginning and end of the virtual helix indices that you want to use discard_unbonded – boolean. If true, discard broken base pairs in the calculation.

get_3d_vecs(vhelix_extent=False, discard_unbonded=True)¶

Returns orthonormal vectors based on the orientation of the origami.

vhelix_extent currently unused…. uses extent of double helices (assuming there is one duplex per virtual helix). NB: the author of this documentation believes that it is actually used. Perhaps the author of the code did not update the comments after vhelix_extent was implemented.

Assumes the skips are lined up across vhelices - in particular it will presumably break if there is a vb with a skip adjacent (in vhelix dimension) to a vb that represents the end of a double helix.

Parameters:	vhelix_extent – 2-element array including the beginning and end of the virtual helix indices that you want to use discard_unbonded – boolean. If true, discard broken base pairs in the calculation. Run get_h_bond_list first.

get_alignment(trim, period)¶

Returns the alignment of a pair of nucleotides, i.e. to see if ones that are supposed to line up do. Assumes no skip/loop / insertion/deletion

Parameters:	trim – number of base pairs to trim at the ends period – separation of pairs of base pairs

get_arm_vec(vhc, vbc, vbn)¶

This function was not documented when it was written. The author of this documentation believes that it returns a vector along an arm of a Holliday junction, or more generally, a vector pointing along a virtual helix.

Parameters:	vhc – virtual helix index in cadnano vbc – starting virtual base index in cadnano vbn – ending virtual base index in cadnano

get_av_helix_axis(vhelix_extent=False, discard_unbonded=True)¶

Returns the average helix axis. Optionally discard unbonded (but designed) ‘base-pairs’, and ignore the regions outside the ‘vhelix_extent’. Run get_h_bond_list first.

Parameters:	vhelix_extent – 2-element array including the beginning and end of the virtual helix indices that you want to use discard_unbonded – boolean. If true, discard broken base pairs in the calculation.

get_backback_midpoint(n_index)¶

Returns the midpoint of 2 backbone vectors corresponding to nucleotides that are hybridised according to cadnano scheme.

Parameters:	n_index – nucleotide index

get_backback_vec(n_index)¶

Returns the vector between 2 backbones corresponding to nucleotides that are hybridised according to cadnano scheme.

Parameters:	n_index – nucleotide index

get_bb_midpoint(n_index, pbc=True)¶

Returns the midpoint of 2 base vectors corresponding to nucleotides that are hybridised according to cadnano scheme.

Parameters:	n_index – nucleotide index pbc – boolean. If true, use the minimum image: i.e. make sure we get a sensible answer if the r1 and r2 correspond to nucleotides which were put at opposite (positions) – of the box due to pbc's. (ends) –

get_bb_vec(n_index)¶

Returns the vector between 2 bases corresponding to nucleotides that are hybridised according to cadnano scheme.

Parameters:	n_index – nucleotide index

get_bpwise_weave(vh_id, vvib_id)¶

Returns the weave for a pair of base pairs located at (vh_id, vvib_id) and (vh_id + 1, vvib_id). Assumes no skip/loop.

Parameters:	vh_id – virtual helix index in vhelix_indices vvib_id – virtual base index in vvib

get_cad2cudadna(infile, visibility=False)¶

This function was not documented when it was written. The author of this documentation believes that it is used for parsing the virt2nuc file when initialising an origami object. Users should not use this function directly.

Parameters:	infile – path to virt2nuc file visibility – path to visibility file

get_com()¶: Not currently used. Unsupported/may no longer work - e.g. take care in using self.width and get_bb_midpoint.

get_corners()¶: This function was not documented when it was written. The author of this documentation believes that it returns the indices of the nucleotides in the four corners of the origami.

get_flat_nucs(vh)¶

This function was not documented when it was written. The author of this documentation believes that it returns a list of nucleotides in a given virtual helix.

Parameters:	vh – virtual helix index in cadnano

get_flat_units(discard_unbonded)¶

Returns the width, height units of the origami if it had been laid out flat.

Parameters:	discard_unbonded – boolean. If true, discard broken base pairs in the calculation. Run get_h_bond_list first.

get_h_bond_list(infile)¶

Fill in self.interaction_list for the current system self._sys; self.interaction_list[nucid1] is a list of nucleotides that nucid1 has hydrogen bonds with

Parameters:	infile – path to input file for DNAnalysis

get_h_bond_list_output_bonds(infile, conffile, conf_num)¶

Deprecated function, use origami_utils.Origami.get_h_bond_list instead.

The old (slow) version using output bonds.

get_hj_alignment(start, period)¶

Find alignment of a pair of nucleotides, i.e. to see if ones that are supposed to line up do.

Assumes exactly 2 virtual helices, no skip/loop / insertion/deletion.

Parameters:	start – virtual base index of the first crossover (reading left to right) period – separation of the pair of base pairs

get_holliday_junctions(single=False, single_only=False)¶

Returns a list of hjs in format (vh, vb, vh_neighbour, vb_next). Assumes no insertions/deletions when using cad2cudadna. If somehow there were two holliday junctions between the same pairs of vbases (one for scaffold and one for staple strand), the 2nd would not be detected here.

Parameters:	single – boolean. If true, returns a list of all single crossovers rather than holliday junctions (i.e. pairs of crossovers). single_only – boolean. If true, returns a list of just single crossovers (double crossovers excluded). If single_only is true, the state of single should not matter.

get_junction_normal(hj, armlen=3)¶

Returns the normal vector of a Holliday junction. Assumes double crossover, no skip/loop (not sure if matters).

Parameters:	hj – should come from hjs[ii] where hjs=origami.get_holliday_junctions().

get_junction_normal2(hj)¶

Returns the normal vector of a Holliday junction. Assumes double crossover, no skip/loop (not sure if matters).

A similar method to get_junction_normal.

Parameters:	hj – should come from hjs[ii] where hjs=origami.get_holliday_junctions().

get_local_twist(vh_id, vvib_id, type, conf=False)¶

Returns the local twist between given vhelix, vbase and the next one along. Assumes no skip/loop / insertion/deletion.

Parameters:	vh_id – virtual helix index in vhelix_indices vvib_id – virtual base index in vvib type – “base” for base or “back” for backbone

get_n_bp(vh)¶

Returns the number of base pairs in a virtual helix.

Parameters:	vh – virtual helix index in cadnano

get_n_nicks()¶: Returns the number of nicks (staple or scaffold) in an origami. Here a nick must be between adjacent strands (e.g. if a nucleotide is missing that isn’t counted as a nick at all).

get_nearest_native_bonded(vhi, vb, increment, bound)¶

Returns the nearest natively bonded base-pair. Run get_h_bond_list first.

Parameters:	vhi – virtual helix index in vhelix_indices vb – virtual base index in cadnano increment – must be +1 or -1, the amount to increment the virtual base by when searching for a virtual base with a native bond bound – if we get to this virtual base and we still haven’t found anything, we give up and throw and error

get_nicks()¶

Returns a list [scaffold_nicks, staple_nicks], each of which is a list of nicked neighbours in the origami. Here a nick must be between adjacent strands (e.g. if a nucleotide is missing that isn’t counted as a nick at all).

An element of each nick list has virtual helix number (NOT array index), virtual base 1, virtual base 2.

get_nucleotides(vhelix, vbase, type='default')¶

Given a cadnano virtual helix index and a virtual base index, returns the nucleotide index

Parameters:	vhelix – virtual helix index in cadnano vbase – virtual base index in cadnano type – “scaf” for scaffold, “stap” for staple, “double” for both

get_plane_vecs(discard_unbonded=True)¶

This function was not documented when it was written. The author of this documentation believes that it sets self.vec_long and self.vec_lat for an origami_utils.Origami object. These are believed to be two vectors defining a plane of an origami.

Parameters:	discard_unbonded – boolean. If true, discard broken base pairs in the calculation. Run get_h_bond_list first.

get_principal_axes(approxaxes, minnucs)¶: Unsupported function. Do not use.

get_radius(vh_id, vvib_id, new_method=False)¶

Returns the radius of a base pair (i.e. the radius of helix). Assumes no skip/loop.

Parameters:	vh_id – virtual helix index in vhelix_indices vvib_id – virtual base index in vvib new_method – boolean. If true, use a new method for getting the radius. It gives the same result for an ideal configuration but a greater radius for an equilibrated configuration.

get_rise(vh_id, vvib_id)¶

Returns the rise for a base pair. Assumes no skip/loop.

Parameters:	vh_id – virtual helix index in vhelix_indices vvib_id – virtual base index in vvib

get_vh_midpoints(vhi)¶

Returns a list of base pair midpoints for a virtual helix. Effectively ‘flattens’ a virtual helix, dealing with any skips or loops.

Parameters:	vhi – virtual helix index in vhelix_indices

get_vh_spline(vh, mode='midpoint', vhelix_extent=False, discard_unbonded=True, force_circular=False)¶

Returns a cartesian spline for the base-base midpoints of a virtual helix.

Parameters:

vh – virtual helix number in cadnano
mode – “midpoint”, which fits a spline through the midpoint between the two strands; or “both strands”, which fits two splines through each strand.
vhelix_extent – 2-element array including the beginning and end of the virtual helix indices that you want to use
discard_unbonded – boolean. If true, discard broken base pairs in the calculation. Run get_h_bond_list first.
force_circular – boolean. If true, force a closed curve; the ends of the midpoint curve (which is open) will be joined together by a straight line.

get_vhelix_ds_length(vhi)¶

Returns length of double helix. Requires one continuous double strand - otherwise how is double strand length defined?

Parameters:	vhi – virtual helix index in self.vhelix_indices

get_vhelix_neighbours(type)¶

This function was not documented when it was written. The author of this documentation believes that it returns a neighbour list of virtual helices.

Parameters:	type – “sq” or “he”, corresponding to square or hexagonal lattice in cadnano

get_weave(verbose=False)¶

Returns the weave pattern as a list of lists. Each pair of virtual helices has a list associated with it. Each of those lists contains the interhelical distance as a function of base-pair index along the virtual helix.

This method copes with insertions/deletions by creating a list of base- base midpoints for each virtual helix, and assuming that each pair of virtual helices has the same number of base-base midpoints. This should mean that the base pairs are fairly well lined up and results in a reasonable definition of the weave pattern.

Parameters:	verbose – boolean. If true, prints skipped virtual helix pairs.

has_full_hbond(vh, vb)¶

Returns true if all cadnano base-pairs at this virtual base have their correct base-pairs in the configuration. Assumes self.interaction_list is already filled. Run get_h_bond_list first.

Parameters:	vh – virtual helix index in cadnano vb – virtual base index in cadnano

intra_helix_autocorr()¶

Returns the autocorrelation (i.e. the dot product) between the intra-helix vectors along each pair of virtual helices in the origami.

This method does assume that the base pairs on adjacent virtual helices ‘line up,’ in the sense that the list index of the closest base pair on one helix is the same as the index of the closest base pair on the neighbouring helix. A weak check for this is the condition that the lengths of the bbm lists (and so the number of base pairs in each virtual helix) are the same.

Note that the list this method returns is indexed by the actual virtual helix number given by cadnano, not the ‘virtual helix index’ which indexes the self.vhelix_indices list.

Run compute_vh_midpoints() first.

is_crossover(vh, vhn, vb, strand_type)¶

Returns a boolean of whether a crossover of type strand_type span the squares (vh, vb) and (vhn, vb).

Parameters:	vh – virtual helix index in cadnano vhn – another virtual helix index in cadnano vb – virtual base index in cadnano strand_type – “scaffold” or “staple”

is_junction(vh, vb1, vb2, hjs_double, hjs_single)¶

Returns ‘double’ if double crossover, ‘single’ if single crossover, False otherwise

Parameters:	vh – virtual helix index in cadnano vb1 – virtual base index in cadnano (left square) vb2 – virtual base index in cadnano (right square) hjs_double – array of holliday junctions (double crossovers). Should come from origami.get_holliday_junctions(). hjs_single – array of single crossovers. Should come from origami.get_holliday_junctions().

is_nicked(vh, vb1, vb2)¶

Returns true if there is a nick across vb1 and vb2, false otherwise.

Parameters:	vh – virtual helix index in cadnano vb1 – first virtual bsae index in cadnano vb2 – second virtual base index in cadnano

is_strand_end(vh, vb, strand_type)¶

Returns string “first” or “last” for first or last nucleotide of the strand, otherwise false. Assumes we don’t have any 1-nucleotide-long strands

Parameters:	vh – virtual helix index in cadnano vb – virtual base index in cadnano strand_type – “scaffold” or “staple”

map_base_to_vbase(vhelices_len, begin_bases, end_bases)¶

Not used right now.

Returns virtual base AND the index of the nucleotide at that virtual base (ie to deal with the possiblity of a loop).

min_distance(r1, r2)¶

Returns the minimum image distance in going from r1 to r2, in this system’s box.

Parameters:	r1 – first vector r2 – second vector

old_get_corrug_weave(vh1, vh2, vb, n, norm2, corrug_signed=False)¶: Returns distance between adjacent helices and decomposes into weave and corrugation for a pair of base pairs located at (vh1, vb) and (vh2, vb). Assumes no skip/loop.

parse_vh_edge(f_edge_vh)¶

Returns a tuple of lists, each element corresponds to a vhelix number that is either in the top (vhelices1) or bottom (vhelices2) of the origami.

Parameters:	f_edge_vh – path to vh_edge file

parse_vhelix_extent_file(vhelix_extentf)¶

Parses the vhelix_extent file and returns a list [vh_begin, vh_end], where vh_begin is a list of length <number of virtual helices> and each element being the left virtual base, and analogous for vh_end.

vhelix_extent file example: 10, 30

This function could be extended to parse a file with different vh_begin and vh_end for each virtual helix; the machinery that uses this function should be able to deal with that without modification.

Parameters:	vhelix_extentf – path to vhelix_extent file

prepare_principal_axes_calc()¶: Unsupported function. Do not use.

square_contains_exactly_one_nucleotide(vh, vb, strand_type)¶

Returns a boolean of whether the square vh, vb contains exactly one nucleotide.

Parameters:	vh – virtual helix index in cadnano vb – virtual base index in cadnano strand_type – “scaffold” or “staple”

update_system(system)¶

This function was not documented when it was written. The author of this documentation believes that it updates the system of the origami object.

Parameters:	system – base.System object. Can be obtained from readers.LorenzoReader.get_system

vb2nuci(vhi, vb0)¶

This function was not documented when it was written. The author of this documentation believes that it returns the number of nucleotides corresponding to a given virtual helix and virtual base. It is not always 1 because of insertions/deletions.

Parameters:	vhi – virtual helix index in self.vhelix_indices vb0 – virtual base index in cadnano (maybe, the author is not sure)

vb2vhelix_nucid(vh, vb, vb_nuci, mode='default')¶

Given a virtual helix, a virtual base, and a nucleotide id for that virtual base, returns the index of that nucleotide in a scheme where the leftmost (lowest virtual base index) nucleotide has index 0 and all nucleotides are given a unique integer index.

Parameters:

vh – virtual helix index in cadnano
vb – virtual base index in cadnano
vb_nuci – index of particular nucleotide in the virtual base (if the virtual base does not have a loop there is always exactly one nucleotide, which has index 0)
mode – gets given to get_nucleotides as type; ‘default’ will count a vb that has either a scaffold or staple strand, while ‘double’ will only count vbs that have both a scaffold and staple strand.

vbase_to_base(vhelix, vbase)¶

Not used right now.

Returns base corresponding to first base at a particular virtual base i.e. taking into account skips and loops (counting from the left of a cadnano diagram).

origami_utils.angle_sense(v1, v2, axis)¶

Returns the angle between two vectors, using a 3rd vector to determine the sense of the angle.

Parameters:	v1 – first vector v2 – second vector axis – third vector

origami_utils.array_index(mylist, myarray)¶

Returns the index in a list mylist of a numpy array myarray.

Parameters:	mylist – list myarray – numpy array

origami_utils.dihedral_angle_sense(v1, v2, axis, sanity_check=False)¶

Returns the dihedral angle between two vectors, using a 3rd vector to determine the sense of the angle and to define the axis normal to the plane we want the vectors in

Parameters:	v1 – first vector v2 – second vector axis – third vector sanity_check – boolean. If true, check whether v1 and v2 are not too parallel to the axis, and returns false if either of them is.

origami_utils.get_base_spline(strand, reverse=False)¶

Returns a cartesian spline that represents a fit through the bases for the strand ‘strand’.

Parameters:	strand – base.Strand object reverse – boolean. If true, strand is reversed

origami_utils.get_bb_midpoint(system, strand, n_index, interaction_list)¶

Deprecated function, use origami_utils.Origami.get_bb_midpoint instead.

Returns the midpoint vector between 2 hybridised bases.

origami_utils.get_nucleotide(vhelix, vbase, vhelix_indices)¶

Deprecated function, use origami_utils.Origami.get_nucleotides instead.

Returns the system nucleotide index of a nucleotide given a position on the origami.

origami_utils.get_pos_midpoint(r1, r2, box)¶

Returns the midpoint of two vectors r1 and r2.

Uses the minimum image: i.e. make sure we get a sensible answer if the positions r1 and r2 correspond to nucleotides which were put at opposite ends of the box due to periodic boundary conditions.

Parameters:	r1 – first vector r2 – second vector box – box dimensions in a numpy array

origami_utils.get_sayar_twist(s1, s2, smin, smax, npoints=1000, circular=False, integral_type='simple')¶

Returns the twist for a given pair of spline fits, one through the bases of each strand.

From Sayar et al. Phys. Rev. E, 81, 041916 (2010)

Just need to integrate along the contour parameter s that is common to both splines. We need the normalised tangent vector to the spline formed by the midpoint of the two input splines t(s), the normalised normal vector formed by the vectors between the splines u(s), and the derivative of the normalised normal vector between the splines d/ds (u(s)). NB, the normal u(s) vector should be orthogonal to the tangent vector t(s); we ensure this by using only the component orthogonal to t(s).

Using integral_type = ‘simple’ and npoints = 200, it will give a correct twist, or at least one that gives a conserved linking number when combined with get_sayar_writhe.

Parameters:

s1 – list of 3 splines corresponding to 3-D spline through strand 1’s bases (e.g. use get_base_spline())
s2 – list of 3 splines corresponding to 3-D spline through strand 2’s bases – NB the splines should run in the same direction, i.e. one must reverse one of the splines if they come from get_base_spline (e.g. use get_base_spline(reverse = True))
smin – minimum value for s, which parameterises the splines
smax – maximum value for s, which parameterises the splines
npoints – number of points for the discrete integration
circular – boolean. If true, spline is circular.
integral_type – “simple”

origami_utils.get_sayar_writhe(splines1, smin, smax, splines2=False, npoints=1000, debug=False, circular=False, integral_type='simple')¶

Returns the writhe for a 3D spline fit through a set of duplex midpoints.

From Sayar et al. Phys. Rev. E, 81, 041916 (2010).

Using integral_type = ‘simple’ and npoints = 200, it will give a correct writhe, or at least one that gives a conserved linking number when combined with get_sayar_twist.

Parameters:

splines1 – list of 3 splines corresponding to either (if not splines2) a 3D spline through the duplex or (if splines2) strand 1’s bases
smin – minimum value for s, which parameterises the splines
smax – maximum value for s, which parameterises the splines
splines2 – optionally, (see splines1) list of 3 splines corresponding to a 3D spline through strand2’s bases
npoints – number of points for the discrete integration
debug – print a load of debugging information
circular – boolean. If true, spline is circular.
integral_type – “simple”

origami_utils.get_scaffold_index(system)¶

Returns the index of the scaffold strand in the system.

Parameters:	system – base.System object. Can be obtained from readers.LorenzoReader.get_system.

origami_utils.get_vhelix_vis(fname)¶

Reads the ‘virtual helix visibility’ from a text file and ignore any virtual helices set to invisible.

The syntax is the same as for the base.py visibility.

Modified from the base.py strand visibility parsing.

Parameters:	fname – path to visibility file

origami_utils.min_distance(r1, r2, box)¶

Returns the minimum image distance in going from r1 to r2, in a box of size box.

Same as base.Nucleotide.distance().

Parameters:	r1 – first vector r2 – second vector box – box dimensions in a numpy array

origami_utils.norm(vec)¶

Returns a normalised vector.

Parameters:	vec – vector

origami_utils.open_arrow_debug_file(filename, type='w')¶: This function was not documented when it was written. The author of this documentation believes that it is a debugging function for VMD.

origami_utils.open_box_debug_file(filename)¶: This function was not documented when it was written. The author of this documentation believes that it is a debugging function for VMD.

origami_utils.parse_scanfile(infile)¶

This function was not documented when it was written. The author of this documentation believes that it parses a “scan file”, the format and purpose of which is unclear.

Parameters:	infile – path to scan file

origami_utils.parse_vh_data(filename, origami)¶

Gets vhelices data from file - format is either <auto,[number of vhelices]>, or, if a region of an origami is to be analysed, <[region width],[list of starting nucleotide index for the region for each vhelix]>.

Parameters:	filename – path to vhelix data file origami – origami_utils.Origami object

origami_utils.print_arrow_debug_line(begin, end, file)¶: This function was not documented when it was written. The author of this documentation believes that it is a debugging function for VMD.

origami_utils.print_box_debug_line(vecs, file)¶: This function was not documented when it was written. The author of this documentation believes that it is a debugging function for VMD.

origami_utils.vecs2spline(vecs, per)¶

This function was not documented when it was written. The author of this documentation believes that it returns a spline interpolating the given vectors.

Parameters:	vecs – list of vectors (arrays of length 3) per – boolean. If true, forces spline to be periodic.