Skip to content

Analysis Classes

Classes for groove width and torsion angle analysis.

mdna.analysis.GrooveAnalysis

Compute minor and major groove widths from a DNA trajectory.

Extracts phosphorus-atom positions from the two DNA strands, fits cubic splines through them, and measures inter-strand distances to derive groove widths for every frame.

Attributes:

Name Type Description
minor_widths ndarray

Minor-groove widths, shape (n_frames, n_points).
major_widths ndarray

Major-groove widths, shape (n_frames, n_points).
sequence list[str]

Sense-strand sequence letters.
base_pairs list[str]

Base-pair labels (e.g. ['G-C', 'A-T', ...]).
Example 
import mdtraj as md
traj = md.load('trajectory.xtc', top='topology.pdb')
grooves = GrooveAnalysis(traj)
fig, ax = plt.subplots()
grooves.plot_groove_widths(ax=ax)
Source code in mdna/analysis.py 
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
class GrooveAnalysis:
    """Compute minor and major groove widths from a DNA trajectory.

    Extracts phosphorus-atom positions from the two DNA strands, fits cubic
    splines through them, and measures inter-strand distances to derive
    groove widths for every frame.

    Attributes:
        minor_widths (numpy.ndarray): Minor-groove widths, shape ``(n_frames, n_points)``.
        major_widths (numpy.ndarray): Major-groove widths, shape ``(n_frames, n_points)``.
        sequence (list[str]): Sense-strand sequence letters.
        base_pairs (list[str]): Base-pair labels (e.g. ``['G-C', 'A-T', ...]``).

    Example:
        ```python
        import mdtraj as md
        traj = md.load('trajectory.xtc', top='topology.pdb')
        grooves = GrooveAnalysis(traj)
        fig, ax = plt.subplots()
        grooves.plot_groove_widths(ax=ax)
        ```
    """

    def __init__(self, raw_traj, points=None, parallel=joblib_available):
        """Initialize groove analysis.

        Args:
            raw_traj (md.Trajectory): Trajectory containing at least two DNA chains.
            points (int, optional): Number of spline interpolation points.
                Defaults to ``(n_bp - 1) * 4``.
            parallel (bool): Use joblib for parallel computation when available.
        """

        self.use_parallel = parallel
        self.points = points
        self.sequence = get_sequence_letters(raw_traj)
        self.base_pairs = get_base_pair_letters(raw_traj)
        self.get_phosphors_only_traj(raw_traj) 

        self.nbp = len(self.DNA_traj.top._residues)//2

        # Fit cubic spline curves based on phosphor atoms in strands
        if self.points:
            pass
        else:
            self.points = (self.nbp-1)*4

        self.fit_cubic_spline()
        self.compute_groove_widths()

    def describe(self):
        """Print a summary of the DNA trajectory (base-pair count and frames)."""
        print(f'Your DNA has {self.nbp} base pairs and contains {self.DNA_traj.n_frames} frames.')

    def get_phosphors_only_traj(self, raw_traj):
        """Extract phosphorus-only sub-trajectories for each strand.

        Populates :attr:`DNA_traj`, :attr:`phos_traj`, :attr:`strand_A`,
        and :attr:`strand_B`.

        Args:
            raw_traj (md.Trajectory): Full trajectory.
        """

        # Assuming first two chains are the DNA strands
        DNA_indices = [i.index for c in raw_traj.top._chains[:2] for i in c.atoms]
        self.DNA_traj = raw_traj.atom_slice(DNA_indices)

         # Select only phosphor atoms
        phos_indices = self.DNA_traj.top.select('name P')
        self.phos_traj = self.DNA_traj.atom_slice(phos_indices).center_coordinates()

        # Make new topology of chains
        phos_chain_A = self.phos_traj.topology._chains[0]
        phos_chain_B = self.phos_traj.topology._chains[1]

        # Split phos traj in two seperate traj for each strand
        self.strand_A = self.phos_traj.atom_slice([i.index for i in phos_chain_A.atoms])
        self.strand_B = self.phos_traj.atom_slice([i.index for i in phos_chain_B.atoms])

    def fit_cubic_spline(self):
        """Fit cubic splines through phosphorus positions and compute inter-strand distances.

        Populates :attr:`pairs`, :attr:`distances`, and :attr:`distance_matrices`.
        """

        # fit cubic spline curve
        curves_A = self.fit_curves(self.strand_A, self.points)
        curves_B = self.fit_curves(self.strand_B, self.points)

        # reshape back to trajectory xyz format
        a_curve_xyz = curves_A.T.swapaxes(1,2)
        b_curve_xyz = curves_B.T.swapaxes(1,2)

        # convert curves to xyz format
        curves = np.hstack((a_curve_xyz,b_curve_xyz))

        # Retrieve "predicted points/particles"
        n_points = curves.shape[1]
        points_strand_A = range(0,(int((n_points/2))))
        points_strand_B = range((int((n_points/2))),(int((n_points))))

        # Generate pairs between the two strands
        #self.pairs = list(itertools.product(points_strand_A, points_strand_B))
        self.pairs = [(point_A, point_B) for point_A in points_strand_A for point_B in points_strand_B]

        # Calculate pair distances
        self.distances = _compute_distance(curves, self.pairs)

        # Reshape pair distances to n x n matrix and account for vdWaals radii of P atoms
        s = self.distances.shape
        self.distance_matrices = self.distances.reshape(s[0],len(points_strand_A),len(points_strand_B))[::-1,:,::-1] - 0.58


    def fit_curves(self, traj, points):
        """Interpolate atom positions with a cubic spline.

        Args:
            traj (md.Trajectory): Single-strand phosphorus trajectory.
            points (int): Number of interpolation points.

        Returns:
            coordinates (numpy.ndarray): Interpolated coordinates.
        """

        # make curve with shape[n_particles, xyz, time]
        xyz = traj.xyz.T.swapaxes(0,1)
        n_particles = len(xyz)
        particle_list = list(range(0,n_particles))

        # initialize predictor based on # of actual particles
        predictor = CubicSpline(particle_list,xyz)

        # points of curve to interpolate the particles 
        return predictor(np.linspace(0,len(xyz),points))

    @staticmethod
    def find_first_local_minimum(arr):
        """Return the first local minimum in *arr*, or ``NaN`` if none exists.

        Args:
            arr (numpy.ndarray): 1-D array of values.

        Returns:
            minimum (float): First local minimum value.
        """
        for i in range(1, len(arr) - 1):
            if arr[i] < arr[i - 1] and arr[i] < arr[i + 1]:
                return arr[i]  # Return the first local minimum found
        return np.NaN  # Return None if no local minimum is found

    @staticmethod
    def split_array(array):
        """Split an anti-diagonal slice into minor and major halves.

        Args:
            array (numpy.ndarray): Anti-diagonal slice.

        Returns:
            halves (tuple[numpy.ndarray, numpy.ndarray]): ``(minor_half, major_half)``.
        """
        # Compute the midpoint of the array
        midpoint = len(array) // 2

        # If the length of the array is odd, exclude the middle element. 
        # If it's even, split the array into two equal halves.
        # Return the two halves that correspond to lower and upper triangles of the distance matrix
        # Reverse the first array because the order of the anti diagonal's in the opposite direction of backone/trailing diagonal
        return array[:midpoint][::-1], array[midpoint + len(array) % 2:]

    @staticmethod
    def get_anti_diagonal_slices(matrix):
        """Extract all anti-diagonal slices from a square matrix.

        Args:
            matrix (numpy.ndarray): Square distance matrix.

        Returns:
            slices (list[numpy.ndarray]): Anti-diagonal slices.
        """
        n = matrix.shape[0] # Get the size of the matrix and skip the first and last anti diagonal
        return [np.diagonal(np.flipud(matrix), offset) for offset in range(-(n - 2), n - 2)]


    def get_minor_major_widths(self, distance_matrix):
        """Compute minor and major groove widths from one distance matrix.

        Args:
            distance_matrix (numpy.ndarray): Inter-strand distance matrix for a single frame.

        Returns:
            widths (tuple[list[float], list[float]]): ``(minor_widths, major_widths)``.
        """

        # Split the distance matrix into anti diagonal slices
        diagonal_slices = self.get_anti_diagonal_slices(distance_matrix)

        minor_widths, major_widths = [],[]
        for slice_ in diagonal_slices:
            minor, major = self.split_array(slice_)
            major_widths.append(self.find_first_local_minimum(major))  # Find local minimum in the major slice and add it to the list
            minor_widths.append(self.find_first_local_minimum(minor))  # Find local minimum in the minor slice and add it to the list 

        return minor_widths, major_widths

    def compute_groove_widths(self):
        """Compute groove widths for all frames (parallel or sequential).

        Populates :attr:`minor_widths` and :attr:`major_widths`.
        """
        if self.use_parallel:
            # Parallelize the computation and subtract 0.58 to account for the vdw radius of the phosphorus atoms
            results = Parallel(n_jobs=-1)(delayed(self.get_minor_major_widths)(distance_matrix) for distance_matrix in self.distance_matrices)
        else:
            # Compute sequentially
            results = [self.get_minor_major_widths(distance_matrix) for distance_matrix in self.distance_matrices]

        minor_widths_batch, major_widths_batch = zip(*results)
        self.minor_widths = np.array(minor_widths_batch)
        self.major_widths = np.array(major_widths_batch)


    def plot_width(self, ax, groove, color=None, std=True, ls='-', lw=0.5):
        """Plot mean groove width with optional standard-deviation shading.

        Args:
            ax (matplotlib.axes.Axes): Target axes.
            groove (numpy.ndarray): Groove-width array, shape ``(n_frames, n_points)``.
            color (str, optional): Line colour.
            std (bool): If True, shade the ±1σ region.
            ls (str): Line style.
            lw (float): Line width.
        """
        # Suppress warnings for mean of empty slice
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", category=RuntimeWarning)    
            # Calculate the mean and standard deviation of major widths
            mean = np.nanmean(groove, axis=0)
            stds = np.nanstd(groove, axis=0)

        ax.plot(mean,color=color,ls=ls,lw=lw)
        if std:
            # Fill the area around the mean for widths
            ax.fill_between(range(len(mean)), mean - stds, mean + stds, alpha=0.25,color=color,ls='-',lw=lw)

    def plot_groove_widths(self, minor=True, major=True, std=True, color='k', c_minor=None, lw=0.5, c_major=None, ax=None, base_labels=True, ls='-'):
        """Plot minor and/or major groove widths.

        Args:
            minor (bool): Plot minor groove widths.
            major (bool): Plot major groove widths.
            std (bool): Show standard-deviation shading.
            color (str): Fallback colour when *c_minor*/*c_major* are not set.
            c_minor (str, optional): Minor-groove colour.
            c_major (str, optional): Major-groove colour.
            lw (float): Line width.
            ax (matplotlib.axes.Axes, optional): Target axes.  A new figure
                is created if ``None``.
            base_labels (bool): Label the x-axis with base-pair identifiers.
            ls (str): Line style.

        Returns:
            result (tuple | None): ``(fig, ax)`` when *ax* was ``None``.
        """

        # Create a figure and axes for plotting
        if ax is None:
            no_ax = True
            _, ax = plt.subplots()
        else:
            no_ax = False

        # If c_minor or c_major are not provided, use the general color
        if c_minor is None:
            c_minor = color
        if c_major is None:
            c_major = color
        if (minor and major) and (c_minor == c_major):
            c_minor = 'cornflowerblue'
            c_major = 'coral'

        if minor:
            self.plot_width(ax, self.minor_widths, std=std, color=c_minor,ls=ls,lw=lw)
        if major:
            self.plot_width(ax, self.major_widths, std=std, color=c_major,ls=ls,lw=lw)

        if base_labels:
            ax.set_xticks(np.linspace(0,len(self.major_widths[0]),len(self.base_pairs)).astype(float))
            ax.set_xticklabels(self.base_pairs,rotation=90)

        if no_ax:
            return _, ax

    def calculate_groove_depths(self):
        # Having defined a groove width by a minimal distance at some point along the nucleic acid fragment,
        # we have to calculate the corresponding groove depth. 

        # At a base pair level, this is defined as the distance from the centre of 
        # the backbone-to-backbone width vector to the mid-point of a vector defining the corresponding base pair. 
        # This vector is constructed using the C8 atom of purines and the C6 atom of pyrimidines 

        # For groove depths half-way between base pair levels,
        # we use the average of the corresponding base pair vector mid-points. 
        pass

    def calculate_groove_depths(self):

        # Calculate the C8 and C6 atom positions of purines and pyrimidines
        c8_indices = self.DNA_traj.top.select('name C8')
        c6_indices = self.DNA_traj.top.select('name C6')
        c8_positions = self.DNA_traj.xyz[:, c8_indices]
        c6_positions = self.DNA_traj.xyz[:, c6_indices]

        # Calculate mid-points of vectors defining base pairs
        base_pair_mid_points = 0.5 * (c8_positions + c6_positions)

        # Calculate backbone-to-backbone width vector centres
        backbone_width_vector_centres = 0.5 * (self.strand_A.xyz + self.strand_B.xyz)

        # Calculate groove depths at base pair level
        groove_depths_base_pair_level = np.linalg.norm(base_pair_mid_points - backbone_width_vector_centres, axis=-1)

        # For groove depths half-way between base pair levels, we use the average of the corresponding base pair vector mid-points
        groove_depths_half_way = 0.5 * (groove_depths_base_pair_level[:-1] + groove_depths_base_pair_level[1:])

        return groove_depths_base_pair_level, groove_depths_half_way

__init__(raw_traj, points=None, parallel=joblib_available)

Initialize groove analysis.

Parameters:

Name Type Description Default
raw_traj Trajectory

Trajectory containing at least two DNA chains.

 required
points int

Number of spline interpolation points. Defaults to (n_bp - 1) * 4.

 None
parallel bool

Use joblib for parallel computation when available.

 joblib_available
Source code in mdna/analysis.py 
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
def __init__(self, raw_traj, points=None, parallel=joblib_available):
    """Initialize groove analysis.

    Args:
        raw_traj (md.Trajectory): Trajectory containing at least two DNA chains.
        points (int, optional): Number of spline interpolation points.
            Defaults to ``(n_bp - 1) * 4``.
        parallel (bool): Use joblib for parallel computation when available.
    """

    self.use_parallel = parallel
    self.points = points
    self.sequence = get_sequence_letters(raw_traj)
    self.base_pairs = get_base_pair_letters(raw_traj)
    self.get_phosphors_only_traj(raw_traj) 

    self.nbp = len(self.DNA_traj.top._residues)//2

    # Fit cubic spline curves based on phosphor atoms in strands
    if self.points:
        pass
    else:
        self.points = (self.nbp-1)*4

    self.fit_cubic_spline()
    self.compute_groove_widths()

compute_groove_widths()

Compute groove widths for all frames (parallel or sequential).

Populates :attr:minor_widths and :attr:major_widths.

Source code in mdna/analysis.py 
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
def compute_groove_widths(self):
    """Compute groove widths for all frames (parallel or sequential).

    Populates :attr:`minor_widths` and :attr:`major_widths`.
    """
    if self.use_parallel:
        # Parallelize the computation and subtract 0.58 to account for the vdw radius of the phosphorus atoms
        results = Parallel(n_jobs=-1)(delayed(self.get_minor_major_widths)(distance_matrix) for distance_matrix in self.distance_matrices)
    else:
        # Compute sequentially
        results = [self.get_minor_major_widths(distance_matrix) for distance_matrix in self.distance_matrices]

    minor_widths_batch, major_widths_batch = zip(*results)
    self.minor_widths = np.array(minor_widths_batch)
    self.major_widths = np.array(major_widths_batch)

describe()

Print a summary of the DNA trajectory (base-pair count and frames).

Source code in mdna/analysis.py 
79
80
81
def describe(self):
    """Print a summary of the DNA trajectory (base-pair count and frames)."""
    print(f'Your DNA has {self.nbp} base pairs and contains {self.DNA_traj.n_frames} frames.')

find_first_local_minimum(arr) staticmethod

Return the first local minimum in arr, or NaN if none exists.

Parameters:

Name Type Description Default
arr ndarray

1-D array of values.

 required

Returns:

Name Type Description
minimum float

First local minimum value.
Source code in mdna/analysis.py 
165
166
167
168
169
170
171
172
173
174
175
176
177
178
@staticmethod
def find_first_local_minimum(arr):
    """Return the first local minimum in *arr*, or ``NaN`` if none exists.

    Args:
        arr (numpy.ndarray): 1-D array of values.

    Returns:
        minimum (float): First local minimum value.
    """
    for i in range(1, len(arr) - 1):
        if arr[i] < arr[i - 1] and arr[i] < arr[i + 1]:
            return arr[i]  # Return the first local minimum found
    return np.NaN  # Return None if no local minimum is found

fit_cubic_spline()

Fit cubic splines through phosphorus positions and compute inter-strand distances.

Populates :attr:pairs, :attr:distances, and :attr:distance_matrices.

Source code in mdna/analysis.py 
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
def fit_cubic_spline(self):
    """Fit cubic splines through phosphorus positions and compute inter-strand distances.

    Populates :attr:`pairs`, :attr:`distances`, and :attr:`distance_matrices`.
    """

    # fit cubic spline curve
    curves_A = self.fit_curves(self.strand_A, self.points)
    curves_B = self.fit_curves(self.strand_B, self.points)

    # reshape back to trajectory xyz format
    a_curve_xyz = curves_A.T.swapaxes(1,2)
    b_curve_xyz = curves_B.T.swapaxes(1,2)

    # convert curves to xyz format
    curves = np.hstack((a_curve_xyz,b_curve_xyz))

    # Retrieve "predicted points/particles"
    n_points = curves.shape[1]
    points_strand_A = range(0,(int((n_points/2))))
    points_strand_B = range((int((n_points/2))),(int((n_points))))

    # Generate pairs between the two strands
    #self.pairs = list(itertools.product(points_strand_A, points_strand_B))
    self.pairs = [(point_A, point_B) for point_A in points_strand_A for point_B in points_strand_B]

    # Calculate pair distances
    self.distances = _compute_distance(curves, self.pairs)

    # Reshape pair distances to n x n matrix and account for vdWaals radii of P atoms
    s = self.distances.shape
    self.distance_matrices = self.distances.reshape(s[0],len(points_strand_A),len(points_strand_B))[::-1,:,::-1] - 0.58

fit_curves(traj, points)

Interpolate atom positions with a cubic spline.

Parameters:

Name Type Description Default
traj Trajectory

Single-strand phosphorus trajectory.

 required
points int

Number of interpolation points.

 required

Returns:

Name Type Description
coordinates ndarray

Interpolated coordinates.
Source code in mdna/analysis.py 
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def fit_curves(self, traj, points):
    """Interpolate atom positions with a cubic spline.

    Args:
        traj (md.Trajectory): Single-strand phosphorus trajectory.
        points (int): Number of interpolation points.

    Returns:
        coordinates (numpy.ndarray): Interpolated coordinates.
    """

    # make curve with shape[n_particles, xyz, time]
    xyz = traj.xyz.T.swapaxes(0,1)
    n_particles = len(xyz)
    particle_list = list(range(0,n_particles))

    # initialize predictor based on # of actual particles
    predictor = CubicSpline(particle_list,xyz)

    # points of curve to interpolate the particles 
    return predictor(np.linspace(0,len(xyz),points))

get_anti_diagonal_slices(matrix) staticmethod

Extract all anti-diagonal slices from a square matrix.

Parameters:

Name Type Description Default
matrix ndarray

Square distance matrix.

 required

Returns:

Name Type Description
slices list[ndarray]

Anti-diagonal slices.
Source code in mdna/analysis.py 
199
200
201
202
203
204
205
206
207
208
209
210
@staticmethod
def get_anti_diagonal_slices(matrix):
    """Extract all anti-diagonal slices from a square matrix.

    Args:
        matrix (numpy.ndarray): Square distance matrix.

    Returns:
        slices (list[numpy.ndarray]): Anti-diagonal slices.
    """
    n = matrix.shape[0] # Get the size of the matrix and skip the first and last anti diagonal
    return [np.diagonal(np.flipud(matrix), offset) for offset in range(-(n - 2), n - 2)]

get_minor_major_widths(distance_matrix)

Compute minor and major groove widths from one distance matrix.

Parameters:

Name Type Description Default
distance_matrix ndarray

Inter-strand distance matrix for a single frame.

 required

Returns:

Name Type Description
widths tuple[list[float], list[float]]

(minor_widths, major_widths).
Source code in mdna/analysis.py 
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
def get_minor_major_widths(self, distance_matrix):
    """Compute minor and major groove widths from one distance matrix.

    Args:
        distance_matrix (numpy.ndarray): Inter-strand distance matrix for a single frame.

    Returns:
        widths (tuple[list[float], list[float]]): ``(minor_widths, major_widths)``.
    """

    # Split the distance matrix into anti diagonal slices
    diagonal_slices = self.get_anti_diagonal_slices(distance_matrix)

    minor_widths, major_widths = [],[]
    for slice_ in diagonal_slices:
        minor, major = self.split_array(slice_)
        major_widths.append(self.find_first_local_minimum(major))  # Find local minimum in the major slice and add it to the list
        minor_widths.append(self.find_first_local_minimum(minor))  # Find local minimum in the minor slice and add it to the list 

    return minor_widths, major_widths

get_phosphors_only_traj(raw_traj)

Extract phosphorus-only sub-trajectories for each strand.

Populates :attr:DNA_traj, :attr:phos_traj, :attr:strand_A, and :attr:strand_B.

Parameters:

Name Type Description Default
raw_traj Trajectory

Full trajectory.

 required
Source code in mdna/analysis.py 
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
def get_phosphors_only_traj(self, raw_traj):
    """Extract phosphorus-only sub-trajectories for each strand.

    Populates :attr:`DNA_traj`, :attr:`phos_traj`, :attr:`strand_A`,
    and :attr:`strand_B`.

    Args:
        raw_traj (md.Trajectory): Full trajectory.
    """

    # Assuming first two chains are the DNA strands
    DNA_indices = [i.index for c in raw_traj.top._chains[:2] for i in c.atoms]
    self.DNA_traj = raw_traj.atom_slice(DNA_indices)

     # Select only phosphor atoms
    phos_indices = self.DNA_traj.top.select('name P')
    self.phos_traj = self.DNA_traj.atom_slice(phos_indices).center_coordinates()

    # Make new topology of chains
    phos_chain_A = self.phos_traj.topology._chains[0]
    phos_chain_B = self.phos_traj.topology._chains[1]

    # Split phos traj in two seperate traj for each strand
    self.strand_A = self.phos_traj.atom_slice([i.index for i in phos_chain_A.atoms])
    self.strand_B = self.phos_traj.atom_slice([i.index for i in phos_chain_B.atoms])

plot_groove_widths(minor=True, major=True, std=True, color='k', c_minor=None, lw=0.5, c_major=None, ax=None, base_labels=True, ls='-')

Plot minor and/or major groove widths.

Parameters:

Name Type Description Default
minor bool

Plot minor groove widths.

 True
major bool

Plot major groove widths.

 True
std bool

Show standard-deviation shading.

 True
color str

Fallback colour when c_minor/c_major are not set.

 'k'
c_minor str

Minor-groove colour.

 None
c_major str

Major-groove colour.

 None
lw float

Line width.

 0.5
ax Axes

Target axes. A new figure is created if None.

 None
base_labels bool

Label the x-axis with base-pair identifiers.

 True
ls str

Line style.

 '-'

Returns:

Name Type Description
result tuple | None

(fig, ax) when ax was None.
Source code in mdna/analysis.py 
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
def plot_groove_widths(self, minor=True, major=True, std=True, color='k', c_minor=None, lw=0.5, c_major=None, ax=None, base_labels=True, ls='-'):
    """Plot minor and/or major groove widths.

    Args:
        minor (bool): Plot minor groove widths.
        major (bool): Plot major groove widths.
        std (bool): Show standard-deviation shading.
        color (str): Fallback colour when *c_minor*/*c_major* are not set.
        c_minor (str, optional): Minor-groove colour.
        c_major (str, optional): Major-groove colour.
        lw (float): Line width.
        ax (matplotlib.axes.Axes, optional): Target axes.  A new figure
            is created if ``None``.
        base_labels (bool): Label the x-axis with base-pair identifiers.
        ls (str): Line style.

    Returns:
        result (tuple | None): ``(fig, ax)`` when *ax* was ``None``.
    """

    # Create a figure and axes for plotting
    if ax is None:
        no_ax = True
        _, ax = plt.subplots()
    else:
        no_ax = False

    # If c_minor or c_major are not provided, use the general color
    if c_minor is None:
        c_minor = color
    if c_major is None:
        c_major = color
    if (minor and major) and (c_minor == c_major):
        c_minor = 'cornflowerblue'
        c_major = 'coral'

    if minor:
        self.plot_width(ax, self.minor_widths, std=std, color=c_minor,ls=ls,lw=lw)
    if major:
        self.plot_width(ax, self.major_widths, std=std, color=c_major,ls=ls,lw=lw)

    if base_labels:
        ax.set_xticks(np.linspace(0,len(self.major_widths[0]),len(self.base_pairs)).astype(float))
        ax.set_xticklabels(self.base_pairs,rotation=90)

    if no_ax:
        return _, ax

plot_width(ax, groove, color=None, std=True, ls='-', lw=0.5)

Plot mean groove width with optional standard-deviation shading.

Parameters:

Name Type Description Default
ax Axes

Target axes.

 required
groove ndarray

Groove-width array, shape (n_frames, n_points).

 required
color str

Line colour.

 None
std bool

If True, shade the ±1σ region.

 True
ls str

Line style.

 '-'
lw float

Line width.

 0.5
Source code in mdna/analysis.py 
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
def plot_width(self, ax, groove, color=None, std=True, ls='-', lw=0.5):
    """Plot mean groove width with optional standard-deviation shading.

    Args:
        ax (matplotlib.axes.Axes): Target axes.
        groove (numpy.ndarray): Groove-width array, shape ``(n_frames, n_points)``.
        color (str, optional): Line colour.
        std (bool): If True, shade the ±1σ region.
        ls (str): Line style.
        lw (float): Line width.
    """
    # Suppress warnings for mean of empty slice
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", category=RuntimeWarning)    
        # Calculate the mean and standard deviation of major widths
        mean = np.nanmean(groove, axis=0)
        stds = np.nanstd(groove, axis=0)

    ax.plot(mean,color=color,ls=ls,lw=lw)
    if std:
        # Fill the area around the mean for widths
        ax.fill_between(range(len(mean)), mean - stds, mean + stds, alpha=0.25,color=color,ls='-',lw=lw)

split_array(array) staticmethod

Split an anti-diagonal slice into minor and major halves.

Parameters:

Name Type Description Default
array ndarray

Anti-diagonal slice.

 required

Returns:

Name Type Description
halves tuple[ndarray, ndarray]

(minor_half, major_half).
Source code in mdna/analysis.py 
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
@staticmethod
def split_array(array):
    """Split an anti-diagonal slice into minor and major halves.

    Args:
        array (numpy.ndarray): Anti-diagonal slice.

    Returns:
        halves (tuple[numpy.ndarray, numpy.ndarray]): ``(minor_half, major_half)``.
    """
    # Compute the midpoint of the array
    midpoint = len(array) // 2

    # If the length of the array is odd, exclude the middle element. 
    # If it's even, split the array into two equal halves.
    # Return the two halves that correspond to lower and upper triangles of the distance matrix
    # Reverse the first array because the order of the anti diagonal's in the opposite direction of backone/trailing diagonal
    return array[:midpoint][::-1], array[midpoint + len(array) % 2:]

mdna.analysis.TorsionAnalysis

Compute backbone torsion angles and BI/BII conformational states.

Calculates the epsilon and zeta dihedral angles of the DNA backbone and classifies each base step as BI or BII based on the sign of epsilon − zeta.

Attributes:

Name Type Description
epsilon ndarray

Epsilon torsion angles, shape (n_frames, n_steps).
zeta ndarray

Zeta torsion angles, shape (n_frames, n_steps).
B_state ndarray

Fraction of BII per base step.
Example 
torsions = TorsionAnalysis(traj)
torsions.B_state  # BII fraction per base step
Source code in mdna/analysis.py 
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
class TorsionAnalysis:
    """Compute backbone torsion angles and BI/BII conformational states.

    Calculates the epsilon and zeta dihedral angles of the DNA backbone
    and classifies each base step as BI or BII based on the sign of
    epsilon − zeta.

    Attributes:
        epsilon (numpy.ndarray): Epsilon torsion angles, shape ``(n_frames, n_steps)``.
        zeta (numpy.ndarray): Zeta torsion angles, shape ``(n_frames, n_steps)``.
        B_state (numpy.ndarray): Fraction of BII per base step.

    Example:
        ```python
        torsions = TorsionAnalysis(traj)
        torsions.B_state  # BII fraction per base step
        ```
    """

    def __init__(self, traj, degrees=True, chain=0):
        """Initialize torsion analysis.

        Args:
            traj (md.Trajectory): DNA-containing trajectory.
            degrees (bool): Report angles in degrees (default) or radians.
            chain (int): Chain index to analyse (0 = sense, 1 = anti-sense).
        """
        self.chain = chain
        self.dna = self.load_trajectory_and_slice_dna(traj)
        self.epsilon, self.zeta = self.compute_BI_BII(degrees=degrees)
        self.B_state = self.get_B_state(self.epsilon - self.zeta)

    def load_trajectory_and_slice_dna(self, traj):
        """Slice a trajectory to keep only canonical DNA residues.

        Args:
            traj (md.Trajectory): Input trajectory.

        Returns:
            dna (md.Trajectory): Sub-trajectory containing DG, DC, DA, and DT residues.
        """
        dna = traj.atom_slice(traj.top.select('resname DG DC DA DT'))
        return dna

    def get_backbone_indices(self, chainid, ref_atoms):
        """Collect backbone atom objects matching *ref_atoms* from a chain.

        Args:
            chainid (int): Chain index.
            ref_atoms (list[str]): Atom names to select (e.g. ``["C3'", "O3'"]``).

        Returns:
            indices (list): Atom objects.
        """
        indices = []
        # find torsions based on the epsilon and zeta atoms
        # finally map the torsions for all base steps 
        if chainid == 0:
            residues = self.dna.top._chains[chainid].residues
        else:
            residues = self.dna.top._chains[chainid]._residues

        for res in residues:
            for at in res.atoms:
                if at.name in ref_atoms:
                    indices.append(at)
        return indices

    def get_torsions(self, indices, ref_atoms):
        """Group sequential atoms into torsion-angle quartets.

        Args:
            indices (list): Ordered backbone atom objects.
            ref_atoms (list[str]): Reference atom-name pattern for one torsion.

        Returns:
            torsions (list[list]): Each element is a four-atom group defining one torsion.
        """
        # Find the chunks based on ref_atoms
        torsions = []
        i = 0
        while i < len(indices):
            ref = [at.name for at in indices[i:i+len(ref_atoms)]]
            if ref == ref_atoms:
                torsions.append(indices[i:i+len(ref_atoms)])
                i += len(ref_atoms)
            else:
                i += 1
        return torsions

    def get_torsion_indices(self, chainid, ref_atoms):
        """Get torsion-angle atom groups for a chain.

        Args:
            chainid (int): Chain index.
            ref_atoms (list[str]): Atom-name pattern.

        Returns:
            torsions (list[list]): Torsion atom groups.
        """
        indices = self.get_backbone_indices(chainid, ref_atoms)
        torsions = self.get_torsions(indices, ref_atoms)
        return torsions

    def convert_torsion_indices_to_atom_indices(self, torsion_indices):
        """Convert atom objects to integer indices for MDTraj.

        Args:
            torsion_indices (list[list]): Torsion atom-object groups.

        Returns:
            atom_indices (list[list[int]]): Integer atom indices.
        """
        atom_indices = []
        for torsion in torsion_indices:
            atom_indices.append([at.index for at in torsion])
        return atom_indices

    def compute_BI_BII(self, degrees=True):
        """Compute epsilon and zeta backbone torsion angles.

        Args:
            degrees (bool): Return angles in degrees.

        Returns:
            angles (tuple[numpy.ndarray, numpy.ndarray]): ``(epsilon, zeta)`` arrays.
        """

        epsilon_atoms = ["C4'","C3'","O3'","P"] 
        zeta_atoms = ["C3'","O3'","P","O5'"]

        epsi_0 = self.get_torsion_indices(0, epsilon_atoms)
        epsi_1 = self.get_torsion_indices(1, epsilon_atoms)
        zeta_0 = self.get_torsion_indices(0, zeta_atoms)
        zeta_1 = self.get_torsion_indices(1, zeta_atoms)

        print(len(epsi_0), len(epsi_1), len(zeta_0), len(zeta_1))

        # From here only the antisense strand is used
        if self.chain == 1:
            e_torsion_indices = self.convert_torsion_indices_to_atom_indices(epsi_1)
            z_torsion_indices = self.convert_torsion_indices_to_atom_indices(zeta_1)
        else:
            e_torsion_indices = self.convert_torsion_indices_to_atom_indices(epsi_0)
            z_torsion_indices = self.convert_torsion_indices_to_atom_indices(zeta_0)

        epsi = md.compute_dihedrals(self.dna, e_torsion_indices)
        zeta = md.compute_dihedrals(self.dna, z_torsion_indices)

        if degrees:
            epsi = np.degrees(epsi)
            zeta = np.degrees(zeta)

        print(epsi.shape, zeta.shape)
        return epsi, zeta

    def get_B_state(self, diff):
        """Classify each base step as BI (0) or BII (1) from epsilon−zeta.

        Args:
            diff (numpy.ndarray): ``epsilon - zeta`` array.

        Returns:
            state (numpy.ndarray): BII fraction per base step.
        """
        state = np.zeros_like(diff)
        state[diff < 0] = 0  # BI
        state[diff > 0] = 1  # BII
        return np.round(np.sum(state,axis=0)/state.shape[0],2)

__init__(traj, degrees=True, chain=0)

Initialize torsion analysis.

Parameters:

Name Type Description Default
traj Trajectory

DNA-containing trajectory.

 required
degrees bool

Report angles in degrees (default) or radians.

 True
chain int

Chain index to analyse (0 = sense, 1 = anti-sense).

 0
Source code in mdna/analysis.py 
375
376
377
378
379
380
381
382
383
384
385
386
def __init__(self, traj, degrees=True, chain=0):
    """Initialize torsion analysis.

    Args:
        traj (md.Trajectory): DNA-containing trajectory.
        degrees (bool): Report angles in degrees (default) or radians.
        chain (int): Chain index to analyse (0 = sense, 1 = anti-sense).
    """
    self.chain = chain
    self.dna = self.load_trajectory_and_slice_dna(traj)
    self.epsilon, self.zeta = self.compute_BI_BII(degrees=degrees)
    self.B_state = self.get_B_state(self.epsilon - self.zeta)

compute_BI_BII(degrees=True)

Compute epsilon and zeta backbone torsion angles.

Parameters:

Name Type Description Default
degrees bool

Return angles in degrees.

 True

Returns:

Name Type Description
angles tuple[ndarray, ndarray]

(epsilon, zeta) arrays.
Source code in mdna/analysis.py 
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
def compute_BI_BII(self, degrees=True):
    """Compute epsilon and zeta backbone torsion angles.

    Args:
        degrees (bool): Return angles in degrees.

    Returns:
        angles (tuple[numpy.ndarray, numpy.ndarray]): ``(epsilon, zeta)`` arrays.
    """

    epsilon_atoms = ["C4'","C3'","O3'","P"] 
    zeta_atoms = ["C3'","O3'","P","O5'"]

    epsi_0 = self.get_torsion_indices(0, epsilon_atoms)
    epsi_1 = self.get_torsion_indices(1, epsilon_atoms)
    zeta_0 = self.get_torsion_indices(0, zeta_atoms)
    zeta_1 = self.get_torsion_indices(1, zeta_atoms)

    print(len(epsi_0), len(epsi_1), len(zeta_0), len(zeta_1))

    # From here only the antisense strand is used
    if self.chain == 1:
        e_torsion_indices = self.convert_torsion_indices_to_atom_indices(epsi_1)
        z_torsion_indices = self.convert_torsion_indices_to_atom_indices(zeta_1)
    else:
        e_torsion_indices = self.convert_torsion_indices_to_atom_indices(epsi_0)
        z_torsion_indices = self.convert_torsion_indices_to_atom_indices(zeta_0)

    epsi = md.compute_dihedrals(self.dna, e_torsion_indices)
    zeta = md.compute_dihedrals(self.dna, z_torsion_indices)

    if degrees:
        epsi = np.degrees(epsi)
        zeta = np.degrees(zeta)

    print(epsi.shape, zeta.shape)
    return epsi, zeta

convert_torsion_indices_to_atom_indices(torsion_indices)

Convert atom objects to integer indices for MDTraj.

Parameters:

Name Type Description Default
torsion_indices list[list]

Torsion atom-object groups.

 required

Returns:

Name Type Description
atom_indices list[list[int]]

Integer atom indices.
Source code in mdna/analysis.py 
460
461
462
463
464
465
466
467
468
469
470
471
472
def convert_torsion_indices_to_atom_indices(self, torsion_indices):
    """Convert atom objects to integer indices for MDTraj.

    Args:
        torsion_indices (list[list]): Torsion atom-object groups.

    Returns:
        atom_indices (list[list[int]]): Integer atom indices.
    """
    atom_indices = []
    for torsion in torsion_indices:
        atom_indices.append([at.index for at in torsion])
    return atom_indices

get_B_state(diff)

Classify each base step as BI (0) or BII (1) from epsilon−zeta.

Parameters:

Name Type Description Default
diff ndarray

epsilon - zeta array.

 required

Returns:

Name Type Description
state ndarray

BII fraction per base step.
Source code in mdna/analysis.py 
512
513
514
515
516
517
518
519
520
521
522
523
524
def get_B_state(self, diff):
    """Classify each base step as BI (0) or BII (1) from epsilon−zeta.

    Args:
        diff (numpy.ndarray): ``epsilon - zeta`` array.

    Returns:
        state (numpy.ndarray): BII fraction per base step.
    """
    state = np.zeros_like(diff)
    state[diff < 0] = 0  # BI
    state[diff > 0] = 1  # BII
    return np.round(np.sum(state,axis=0)/state.shape[0],2)

get_backbone_indices(chainid, ref_atoms)

Collect backbone atom objects matching ref_atoms from a chain.

Parameters:

Name Type Description Default
chainid int

Chain index.

 required
ref_atoms list[str]

Atom names to select (e.g. ["C3'", "O3'"]).

 required

Returns:

Name Type Description
indices list

Atom objects.
Source code in mdna/analysis.py 
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
def get_backbone_indices(self, chainid, ref_atoms):
    """Collect backbone atom objects matching *ref_atoms* from a chain.

    Args:
        chainid (int): Chain index.
        ref_atoms (list[str]): Atom names to select (e.g. ``["C3'", "O3'"]``).

    Returns:
        indices (list): Atom objects.
    """
    indices = []
    # find torsions based on the epsilon and zeta atoms
    # finally map the torsions for all base steps 
    if chainid == 0:
        residues = self.dna.top._chains[chainid].residues
    else:
        residues = self.dna.top._chains[chainid]._residues

    for res in residues:
        for at in res.atoms:
            if at.name in ref_atoms:
                indices.append(at)
    return indices

get_torsion_indices(chainid, ref_atoms)

Get torsion-angle atom groups for a chain.

Parameters:

Name Type Description Default
chainid int

Chain index.

 required
ref_atoms list[str]

Atom-name pattern.

 required

Returns:

Name Type Description
torsions list[list]

Torsion atom groups.
Source code in mdna/analysis.py 
446
447
448
449
450
451
452
453
454
455
456
457
458
def get_torsion_indices(self, chainid, ref_atoms):
    """Get torsion-angle atom groups for a chain.

    Args:
        chainid (int): Chain index.
        ref_atoms (list[str]): Atom-name pattern.

    Returns:
        torsions (list[list]): Torsion atom groups.
    """
    indices = self.get_backbone_indices(chainid, ref_atoms)
    torsions = self.get_torsions(indices, ref_atoms)
    return torsions

get_torsions(indices, ref_atoms)

Group sequential atoms into torsion-angle quartets.

Parameters:

Name Type Description Default
indices list

Ordered backbone atom objects.

 required
ref_atoms list[str]

Reference atom-name pattern for one torsion.

 required

Returns:

Name Type Description
torsions list[list]

Each element is a four-atom group defining one torsion.
Source code in mdna/analysis.py 
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
def get_torsions(self, indices, ref_atoms):
    """Group sequential atoms into torsion-angle quartets.

    Args:
        indices (list): Ordered backbone atom objects.
        ref_atoms (list[str]): Reference atom-name pattern for one torsion.

    Returns:
        torsions (list[list]): Each element is a four-atom group defining one torsion.
    """
    # Find the chunks based on ref_atoms
    torsions = []
    i = 0
    while i < len(indices):
        ref = [at.name for at in indices[i:i+len(ref_atoms)]]
        if ref == ref_atoms:
            torsions.append(indices[i:i+len(ref_atoms)])
            i += len(ref_atoms)
        else:
            i += 1
    return torsions

load_trajectory_and_slice_dna(traj)

Slice a trajectory to keep only canonical DNA residues.

Parameters:

Name Type Description Default
traj Trajectory

Input trajectory.

 required

Returns:

Name Type Description
dna Trajectory

Sub-trajectory containing DG, DC, DA, and DT residues.
Source code in mdna/analysis.py 
388
389
390
391
392
393
394
395
396
397
398
def load_trajectory_and_slice_dna(self, traj):
    """Slice a trajectory to keep only canonical DNA residues.

    Args:
        traj (md.Trajectory): Input trajectory.

    Returns:
        dna (md.Trajectory): Sub-trajectory containing DG, DC, DA, and DT residues.
    """
    dna = traj.atom_slice(traj.top.select('resname DG DC DA DT'))
    return dna