<p>on OxDNA's documentation, i found a code for determining persistence length. However, i have ran simulations for different number of bases on OxDNA and have the following text file as results: Columns are: Time, X, Y, Z</p>
<p>Few lines:</p>
<pre><code>Time.    X.      Y.     Z
0.0       5.      5.     5
1.7e-13. 6.5.    7.5.   2.5
1.9e-13.  7.       8.    8.5
</code></pre>
<p>Basically, it gives me the 3D position of the double strand DNA molecule as a function of time.</p>
<p>I'm unsure how to use data from this text file to determine the persistence length using the OxDNA code.</p>
<p>Also, i'm wondering if there is a formula that can help to get persistence length from the data?</p>
<pre><code>from os import path
from typing import Dict, Tuple
from collections import namedtuple
import argparse
import numpy as np
import matplotlib.pyplot as plt
import oxpy
from oxDNA_analysis_tools.UTILS.RyeReader import describe
from oxDNA_analysis_tools.UTILS.oat_multiprocesser import oat_multiprocesser
from oxDNA_analysis_tools.UTILS.data_structures import TrajInfo
from oxDNA_analysis_tools.UTILS.logger import log, logger_settings

import time
start_time = time.time()

ComputeContext = namedtuple(&quot;ComputeContext&quot;,[&quot;traj_info&quot;,
                                              &quot;input_file&quot;,
                                              &quot;n1&quot;,
                                              &quot;n2&quot;])

def get_r(conf, nucid:int, pair_dict:Dict) -&gt; np.ndarray:
    &quot;&quot;&quot;
        Returns a normalized vector pointing from base midpoint of the `nucid`-th base pair to the midpoint of the next base pair.
        
        Parameters:
            conf (oxpy.config_info) : The current configuration
            nucid (int) : ID of the nucleotide in the first strand to compute the vector from

        Returns:
            (np.ndarray) : The vector pointing from the midpoint of nucid's base pair to the +1 base pair.
    &quot;&quot;&quot;
    box = np.array(conf.box.box_sides)
    pair = pair_dict[nucid]
    next_pair = pair_dict[nucid+1]

    firstA = conf.particles()[nucid].base_site()
    firstB = conf.particles()[pair].base_site()
    secondA = conf.particles()[nucid+1].base_site()
    secondB = conf.particles()[next_pair].base_site()

    first_midpos = (firstA + firstB) / 2
    second_midpos = (secondA + secondB) / 2 

    r = second_midpos - first_midpos        
    r -= box * np.rint(r / box)
    return r

def compute(ctx:ComputeContext, chunk_size:int, chunk_id:int):
    l0 = 0
    correlations = np.zeros(ctx.n2 - ctx.n1)
    correlations_counter = np.zeros_like(correlations)
    # Read in a conf and identify h-bonds
    with oxpy.Context():
        inp = oxpy.InputFile()
        inp.init_from_filename(ctx.input_file)
        inp[&quot;list_type&quot;] = &quot;cells&quot;
        inp[&quot;trajectory_file&quot;] = ctx.traj_info.path
        inp[&quot;analysis_bytes_to_skip&quot;] = str(ctx.traj_info.idxs[chunk_id*chunk_size].offset)
        inp[&quot;confs_to_analyse&quot;] = str(chunk_size)
        inp[&quot;analysis_data_output_1&quot;] = '{ \n name = stdout \n print_every = 1e10 \n col_1 = { \n id = my_obs \n type = hb_list \n } \n }'

        if (not inp[&quot;use_average_seq&quot;] or inp.get_bool(&quot;use_average_seq&quot;)) and &quot;RNA&quot; in inp[&quot;interaction_type&quot;]:
            log(&quot;Sequence dependence not set for RNA model, wobble base pairs will be ignored&quot;, level='warning')

        backend = oxpy.analysis.AnalysisBackend(inp)
        if backend.config_info().particles()[ctx.n1].strand_id != backend.config_info().particles()[ctx.n2].strand_id:
            raise RuntimeError(&quot;nucid_1 and nucid_2 must be on the same strand&quot;)

        while backend.read_next_configuration():
            pairs = backend.config_info().get_observable_by_id(&quot;my_obs&quot;).get_output_string(backend.config_info().current_step).strip().split('\n')[1:]

            # Extract paired nucleotides into a dict
            pair_dict = {}
            for p in pairs:
                p1 = int(p.split()[0])
                p2 = int(p.split()[1])
                if p1 &lt; ctx.n1 or p1 &gt; ctx.n2:
                    continue
                pair_dict[p1] = p2


            for j in range(ctx.n1, ctx.n2):
                # If there's no base pair, there's no midpoint
                if not j in pair_dict or not j+1 in pair_dict:
                    print(&quot;WARNING: Nucleotide {} or {} is unpaired.  Skipping...&quot;.format(j, j+1))
                    continue

                # Get the midpoint of the base pairs
                r0 = get_r(backend.config_info(), j, pair_dict)
                l0 += np.linalg.norm(r0)
                r0 = r0 / np.linalg.norm(r0)

                # Get the correlation between r0 and each subsequence base step
                for k in range(j, ctx.n2):
                    if not k in pair_dict or not k+1 in pair_dict:
                        continue

                    rk = get_r(backend.config_info(), k, pair_dict)
                    rk = rk / np.linalg.norm(rk)
                    correlations[k-j] += np.dot(r0, rk)
                    correlations_counter[k-j] += 1

    return(l0, correlations, correlations_counter)

def persistence_length(traj_info:TrajInfo, inp_file:str, n1:int, n2:int, ncpus:int=1) -&gt; Tuple[float, np.ndarray]:
    &quot;&quot;&quot;
        Computes the persistence length of a bonded sequence of nucleotides.

        Note that while n1 and n2 must be on the same strand, there can be multiple distinct duplexes within the strand.

        Only paired nucleotides will be considered in the persistence length calculation

        Parameters:
            traj_info (TrajInfo): TrajInfo object for the trajectory you want to analyze
            inp_file (str): The path to the input file used to run the simulation
            n1 (int): ID of the particle to start the analysis
            n2 (int): ID of the particle to end the analysis. 
        
        Returns:
            (Tuple[float, np.ndarray]) : Tuple containing the average contour length and correlations between each pair step
    &quot;&quot;&quot;

    ctx = ComputeContext(traj_info, inp_file, n1, n2)
    l0 = 0
    correlations = np.zeros(n2 - n1)
    correlations_counter = np.zeros_like(correlations)

    def callback(i, r):
        nonlocal l0, correlations, correlations_counter
        l0 += r[0]
        correlations += r[1]
        correlations_counter += r[2]

    oat_multiprocesser(traj_info.nconfs, ncpus, compute, callback, ctx)

    l0 /= traj_info.nconfs
    correlations /= correlations_counter
    return(l0, correlations)

def fit_PL(correlations:np.ndarray, plt_name:str) -&gt; float:
    &quot;&quot;&quot;
        Fits persistence length from tangent vector correlations

        Parameters:
            correlations (np.ndarray) : Array of offsets vs correlation
            plt_name (str) : Name to save the resulting plot to

        Returns:
            (float) : Persistence length in nucleotides
    &quot;&quot;&quot;
    # Fit the PL to the correlations
    x = np.arange(0, len(correlations))
    log_corr = np.log(correlations)
    A, B = np.polyfit(x, log_corr, 1)
    pl = -1/A

    # Make a plot
    fig, ax = plt.subplots()
    ax.scatter(x, log_corr, alpha=0.5)
    trend = np.poly1d([A, B])
    ax.plot(x, trend(x), label=f'Fit y = (-x/{pl:.1f})+{B:.3f}')
    ax.spines[['top', 'right']].set_visible(False)
    ax.set_xlabel('Offset')
    ax.set_ylabel('ln(correlation)')
    plt.legend()
    plt.tight_layout()
    plt.savefig(plt_name, dpi=300)
    log(f&quot;Saving figure to {plt_name}&quot;)

    return pl

def cli_parser(prog=&quot;persistence_length.py&quot;):
    parser = argparse.ArgumentParser(prog=prog, description=&quot;Calculates persistence length and contour length of a paired sequence of DNA.&quot;)
    parser.add_argument('traj_file', type=str, nargs=1, help=&quot;The trajectory file to analyze&quot;)
    parser.add_argument('input', type=str, nargs=1, help=&quot;The input file associated with the trajectory&quot;)
    parser.add_argument('nucid_1', type=int, nargs=1, help=&quot;Nucleotide ID to start the calculation at&quot;)
    parser.add_argument('nucid_2', type=int, nargs=1, help=&quot;Nucleotide ID to stop the calculation at&quot;)
    parser.add_argument('-p', '--parallel', metavar='num_cpus', nargs=1, type=int, dest='parallel', help=&quot;(optional) How many cores to use&quot;)
    parser.add_argument('-d', '--data', metavar='data_file', nargs=1, help='If set, the correlations will be written to a txt file in the format `offset correlation`')
    parser.add_argument('-n', '--plot_name', nargs=1, help='Name to save the plot showing the fit of persistence length to correlations.  Defaults to persistence_length.png')
    parser.add_argument('-q', metavar='quiet', dest='quiet', action='store_const', const=True, default=False, help=&quot;Don't print 'INFO' messages to stderr&quot;)
    return parser

def main():
    parser = cli_parser(path.basename(__file__))
    args = parser.parse_args()

    logger_settings.set_quiet(args.quiet)
    from oxDNA_analysis_tools.config import check
    check([&quot;python&quot;, &quot;matplotlib&quot;, &quot;numpy&quot;])

    traj = args.traj_file[0]
    inp_file = args.input[0]
    n1 = args.nucid_1[0]
    n2 = args.nucid_2[0]

    # -p sets the number of cpus to use
    if args.parallel:
        ncpus = args.parallel[0]
    else:
        ncpus = 1

    ti, di = describe(None, traj)

    l0, correlations = persistence_length(di, inp_file, n1, n2, ncpus)

    if args.data:
        log(f&quot;Writing correlations to {args.data[0]}&quot;)
        with open(args.data[0], 'w+') as f:
            for i, c in enumerate(correlations):
                f.write(f'{i} {c}\n')

    if args.plot_name:
        plot_name = args.plot_name[0]
    else:
        plot_name = 'persistence_length.png'

    pl = fit_PL(correlations, plot_name)
    print(&quot;Persistence length: {:.1f} nucleotides&quot;.format(pl))
    print(&quot;Overall bonded contour length between n1 and n2 is:&quot;, l0*0.8518, &quot;nm&quot;)

    print(&quot;--- %s seconds ---&quot; % (time.time() - start_time))
                
if __name__ == '__main__':
    main()
</code></pre>
<p>Link: <a href="https://github.com/lorenzo-rovigatti/oxDNA/blob/master/analysis/src/oxDNA_analysis_tools/persistence_length.py" rel="nofollow noreferrer">https://github.com/lorenzo-rovigatti/oxDNA/blob/master/analysis/src/oxDNA_analysis_tools/persistence_length.py</a></p>

CC BY-SA 4.0
2024-08-01T19:45:59.853Z
Thanks. Can you please show the full command line you're using to run this code?
gringer
CC BY-SA 4.0
2024-08-01T19:48:12.447Z
The file example you have written here has dots separating fields, but the previous "columns are" line suggest it should be comma-separated. Which matches the actual file you are using?
gringer
<python><modelling><dna>
Finding persistence length using OxDNA on python
# Finding persistence length using OxDNA on python

(Asked by: username_0 on 2024-07-29)

on OxDNA's documentation, i found a code for determining persistence length. However, i have ran simulations for different number of bases on OxDNA and have the following text file as results: Columns are: Time, X, Y, Z
Few lines:
Time.    X.      Y.     Z
0.0       5.      5.     5
1.7e-13. 6.5.    7.5.   2.5
1.9e-13.  7.       8.    8.5

Basically, it gives me the 3D position of the double strand DNA molecule as a function of time.
I'm unsure how to use data from this text file to determine the persistence length using the OxDNA code.
Also, i'm wondering if there is a formula that can help to get persistence length from the data?
from os import path
from typing import Dict, Tuple
from collections import namedtuple
import argparse
import numpy as np
import matplotlib.pyplot as plt
import oxpy
from oxDNA_analysis_tools.UTILS.RyeReader import describe
from oxDNA_analysis_tools.UTILS.oat_multiprocesser import oat_multiprocesser
from oxDNA_analysis_tools.UTILS.data_structures import TrajInfo
from oxDNA_analysis_tools.UTILS.logger import log, logger_settings

import time
start_time = time.time()

ComputeContext = namedtuple("ComputeContext",["traj_info",
                                              "input_file",
                                              "n1",
                                              "n2"])

def get_r(conf, nucid:int, pair_dict:Dict) -> np.ndarray:
    """
        Returns a normalized vector pointing from base midpoint of the `nucid`-th base pair to the midpoint of the next base pair.
        
        Parameters:
            conf (oxpy.config_info) : The current configuration
            nucid (int) : ID of the nucleotide in the first strand to compute the vector from

        Returns:
            (np.ndarray) : The vector pointing from the midpoint of nucid's base pair to the +1 base pair.
    """
    box = np.array(conf.box.box_sides)
    pair = pair_dict[nucid]
    next_pair = pair_dict[nucid+1]

    firstA = conf.particles()[nucid].base_site()
    firstB = conf.particles()[pair].base_site()
    secondA = conf.particles()[nucid+1].base_site()
    secondB = conf.particles()[next_pair].base_site()

    first_midpos = (firstA + firstB) / 2
    second_midpos = (secondA + secondB) / 2 

    r = second_midpos - first_midpos        
    r -= box * np.rint(r / box)
    return r

def compute(ctx:ComputeContext, chunk_size:int, chunk_id:int):
    l0 = 0
    correlations = np.zeros(ctx.n2 - ctx.n1)
    correlations_counter = np.zeros_like(correlations)
    # Read in a conf and identify h-bonds
    with oxpy.Context():
        inp = oxpy.InputFile()
        inp.init_from_filename(ctx.input_file)
        inp["list_type"] = "cells"
        inp["trajectory_file"] = ctx.traj_info.path
        inp["analysis_bytes_to_skip"] = str(ctx.traj_info.idxs[chunk_id*chunk_size].offset)
        inp["confs_to_analyse"] = str(chunk_size)
        inp["analysis_data_output_1"] = '{ \n name = stdout \n print_every = 1e10 \n col_1 = { \n id = my_obs \n type = hb_list \n } \n }'

        if (not inp["use_average_seq"] or inp.get_bool("use_average_seq")) and "RNA" in inp["interaction_type"]:
            log("Sequence dependence not set for RNA model, wobble base pairs will be ignored", level='warning')

        backend = oxpy.analysis.AnalysisBackend(inp)
        if backend.config_info().particles()[ctx.n1].strand_id != backend.config_info().particles()[ctx.n2].strand_id:
            raise RuntimeError("nucid_1 and nucid_2 must be on the same strand")

        while backend.read_next_configuration():
            pairs = backend.config_info().get_observable_by_id("my_obs").get_output_string(backend.config_info().current_step).strip().split('\n')[1:]

            # Extract paired nucleotides into a dict
            pair_dict = {}
            for p in pairs:
                p1 = int(p.split()[0])
                p2 = int(p.split()[1])
                if p1 < ctx.n1 or p1 > ctx.n2:
                    continue
                pair_dict[p1] = p2


            for j in range(ctx.n1, ctx.n2):
                # If there's no base pair, there's no midpoint
                if not j in pair_dict or not j+1 in pair_dict:
                    print("WARNING: Nucleotide {} or {} is unpaired.  Skipping...".format(j, j+1))
                    continue

                # Get the midpoint of the base pairs
                r0 = get_r(backend.config_info(), j, pair_dict)
                l0 += np.linalg.norm(r0)
                r0 = r0 / np.linalg.norm(r0)

                # Get the correlation between r0 and each subsequence base step
                for k in range(j, ctx.n2):
                    if not k in pair_dict or not k+1 in pair_dict:
                        continue

                    rk = get_r(backend.config_info(), k, pair_dict)
                    rk = rk / np.linalg.norm(rk)
                    correlations[k-j] += np.dot(r0, rk)
                    correlations_counter[k-j] += 1

    return(l0, correlations, correlations_counter)

def persistence_length(traj_info:TrajInfo, inp_file:str, n1:int, n2:int, ncpus:int=1) -> Tuple[float, np.ndarray]:
    """
        Computes the persistence length of a bonded sequence of nucleotides.

        Note that while n1 and n2 must be on the same strand, there can be multiple distinct duplexes within the strand.

        Only paired nucleotides will be considered in the persistence length calculation

        Parameters:
            traj_info (TrajInfo): TrajInfo object for the trajectory you want to analyze
            inp_file (str): The path to the input file used to run the simulation
            n1 (int): ID of the particle to start the analysis
            n2 (int): ID of the particle to end the analysis. 
        
        Returns:
            (Tuple[float, np.ndarray]) : Tuple containing the average contour length and correlations between each pair step
    """

    ctx = ComputeContext(traj_info, inp_file, n1, n2)
    l0 = 0
    correlations = np.zeros(n2 - n1)
    correlations_counter = np.zeros_like(correlations)

    def callback(i, r):
        nonlocal l0, correlations, correlations_counter
        l0 += r[0]
        correlations += r[1]
        correlations_counter += r[2]

    oat_multiprocesser(traj_info.nconfs, ncpus, compute, callback, ctx)

    l0 /= traj_info.nconfs
    correlations /= correlations_counter
    return(l0, correlations)

def fit_PL(correlations:np.ndarray, plt_name:str) -> float:
    """
        Fits persistence length from tangent vector correlations

        Parameters:
            correlations (np.ndarray) : Array of offsets vs correlation
            plt_name (str) : Name to save the resulting plot to

        Returns:
            (float) : Persistence length in nucleotides
    """
    # Fit the PL to the correlations
    x = np.arange(0, len(correlations))
    log_corr = np.log(correlations)
    A, B = np.polyfit(x, log_corr, 1)
    pl = -1/A

    # Make a plot
    fig, ax = plt.subplots()
    ax.scatter(x, log_corr, alpha=0.5)
    trend = np.poly1d([A, B])
    ax.plot(x, trend(x), label=f'Fit y = (-x/{pl:.1f})+{B:.3f}')
    ax.spines[['top', 'right']].set_visible(False)
    ax.set_xlabel('Offset')
    ax.set_ylabel('ln(correlation)')
    plt.legend()
    plt.tight_layout()
    plt.savefig(plt_name, dpi=300)
    log(f"Saving figure to {plt_name}")

    return pl

def cli_parser(prog="persistence_length.py"):
    parser = argparse.ArgumentParser(prog=prog, description="Calculates persistence length and contour length of a paired sequence of DNA.")
    parser.add_argument('traj_file', type=str, nargs=1, help="The trajectory file to analyze")
    parser.add_argument('input', type=str, nargs=1, help="The input file associated with the trajectory")
    parser.add_argument('nucid_1', type=int, nargs=1, help="Nucleotide ID to start the calculation at")
    parser.add_argument('nucid_2', type=int, nargs=1, help="Nucleotide ID to stop the calculation at")
    parser.add_argument('-p', '--parallel', metavar='num_cpus', nargs=1, type=int, dest='parallel', help="(optional) How many cores to use")
    parser.add_argument('-d', '--data', metavar='data_file', nargs=1, help='If set, the correlations will be written to a txt file in the format `offset correlation`')
    parser.add_argument('-n', '--plot_name', nargs=1, help='Name to save the plot showing the fit of persistence length to correlations.  Defaults to persistence_length.png')
    parser.add_argument('-q', metavar='quiet', dest='quiet', action='store_const', const=True, default=False, help="Don't print 'INFO' messages to stderr")
    return parser

def main():
    parser = cli_parser(path.basename(__file__))
    args = parser.parse_args()

    logger_settings.set_quiet(args.quiet)
    from oxDNA_analysis_tools.config import check
    check(["python", "matplotlib", "numpy"])

    traj = args.traj_file[0]
    inp_file = args.input[0]
    n1 = args.nucid_1[0]
    n2 = args.nucid_2[0]

    # -p sets the number of cpus to use
    if args.parallel:
        ncpus = args.parallel[0]
    else:
        ncpus = 1

    ti, di = describe(None, traj)

    l0, correlations = persistence_length(di, inp_file, n1, n2, ncpus)

    if args.data:
        log(f"Writing correlations to {args.data[0]}")
        with open(args.data[0], 'w+') as f:
            for i, c in enumerate(correlations):
                f.write(f'{i} {c}\n')

    if args.plot_name:
        plot_name = args.plot_name[0]
    else:
        plot_name = 'persistence_length.png'

    pl = fit_PL(correlations, plot_name)
    print("Persistence length: {:.1f} nucleotides".format(pl))
    print("Overall bonded contour length between n1 and n2 is:", l0*0.8518, "nm")

    print("--- %s seconds ---" % (time.time() - start_time))
                
if __name__ == '__main__':
    main()

Link: https://github.com/lorenzo-rovigatti/oxDNA/blob/master/analysis/src/oxDNA_analysis_tools/persistence_length.py

**Tags:** python, modelling, dna
> **Comments:**
> **username_1**: Thanks. Can you please show the full command line you're using to run this code?  
> **username_1**: The file example you have written here has dots separating fields, but the previous "columns are" line suggest it should be comma-separated. Which matches the actual file you are using?  

---

## 0 answers

