Samples: B-form DNA, 6 Atracts and a region of increased flexibility

Features found in this sample CYCLIZE script:
  • uses the &read_sequence method of %seq generation!!
  • Creates a 156nt long B-form DNA.
  • Adds 6 in phase Atracts between positions 25-82.
  • Adds a regions of increased tilt and roll flexibility (9.2 versus 4.842) between positions 101-131.
  • Runs 1 simulation of 1x10^8 whole chains.
    • Notes:
  • This script is identical to the previous sample script, except it uses &read_sequence to generate the %seq hash, rather than by explicitly calculating the flex site positions.
  • If you want to only run multiple calculations, simply change the line that says foreach $repeat (1) to say foreach $repeat (1..10).
  • The value of $cyclize_parameters{whole_chains} can be adjusted to calculate more or less whole chains.
  • Notice the new parameter hash definition for %FLEX and the calls to the "&fill_name" and "&fill_params", which work to add this site of flexibility to a region of the DNA that was originally filled with standard B-form values. The list @flex_positions hold the nt numbers of where the flex site should be located.
    • To obtain this file, highlight the perl script below and copy/paste it into your favorite text editor, or click here for the raw text file (hold the shift key to save to disk).

    You will also need the text file containing the actual nucleotide sequence, get it here. 

    #! /usr/bin/perl -w
# 
require $ENV{'CYCLIZE_HOME'}."lib/cyclize_lib.pl";

# Build the cyclize parameters hash
%cyclize_parameters = (
    "whole_chains"     => 1e8,
    "nrad_stats"       => 1e7,
    "icalcs"           => 100,
    "radial_cutoff"    => 60,
    "axial_cutoff"     => 40,
    "torsional_cutoff" => 36,
    "nkeepers"         => 10,
);

# Set the default values for generating the DNA
$Temp       = 300;

# Build the parameter hashes defined below
&build_parameter_hashes; 

# Read in the sequence file
$file = "12a156.seq";
%seq = &read_sequence( $file );

# Assign B-form parameters to the whole thing
%seq = &fill_params( \%seq, \%BDNA );

# Find and set the A-tract parameters
%seq = &fill_Atract_params( \%seq, \%Atract, \%Atract_3, \%Atract_5 );

# Find the flex site
$flex_sequence = "ACGCCTATAAACGCCTATAAACGCCTTGCTC";
@flex_positions = &find_positions( \%seq, $flex_sequence );

# Set the flex site parameters
%seq = &fill_params( \%seq, \%FLEX, @flex_positions );

# &print_params( \%seq, "sequence.in" );
# exit;

for $repeat ( 1 ) {
    &Cyclize( \%seq, \%cyclize_parameters, "flex_example_$repeat" );
}

exit;

############################################################
# build_parameters_hashes
# 
# USAGE: 
# &build_parameter_hashes;
#
# -This subroutine initialize the defaults for the "parameter hashes" used
#  to construct the "sequence hash".  
# -A single parameter hash contains six keys: tilt, tilt_flex, roll,
#  roll_flex, twist and twist_flex.
############################################################
sub build_parameter_hashes {

    #####################################################
    # BDNA site parameter definitions
    $BDNA_tormod         = 2.4e-19; # torsional modulus (erg*cm)
    $BDNA_rpb            = 3.4e-8;  # rise/bp (cm)
    $BDNA_helical_repeat = 10.45;   # helical repeat (bp/360 degrees)

    %BDNA = (
        tilt       => 0,
        roll       => 0,
        twist      => 360/$BDNA_helical_repeat,
        tilt_flex  => 4.842,
        roll_flex  => 4.842,
        twist_flex => sqrt( $BDNA_rpb*$k*$Temp/$BDNA_tormod )*180/$PI,
        dz         => 3.4,
    );

    #####################################################
    # FLEX site parameter definitions
    $FLEX_tormod         = 2.4e-19; # torsional modulus (erg*cm)
    $FLEX_rpb            = 3.4e-8;  # rise/bp (cm)
    $FLEX_helical_repeat = 10.45;   # helical repeat (bp/360 degrees)

    %FLEX = (
        tilt       => 0,
        roll       => 0,
        twist      => 360/$FLEX_helical_repeat,
        tilt_flex  => 9.2,
        roll_flex  => 9.2,
        twist_flex => sqrt( $FLEX_rpb*$k*$Temp/$FLEX_tormod )*180/$PI,
        dz         => 3.4,
    );
    
    #####################################################
    # Atract site parameter definitions
    $Atract_tormod         = 2.4e-19; # torsional modulus (erg*cm)
    $Atract_rpb            = 3.4e-8;  # rise/bp (cm)
    $Atract_helical_repeat = 10.33;   # helical repeat (bp/360 degrees)

    # Default Atract values for middle of Atract region
    %Atract = (
        tilt       => 0,        
        roll       => 0,        
        twist      => 360/$Atract_helical_repeat,
        tilt_flex  => 4.842,
        roll_flex  => 4.842,
        twist_flex => sqrt( $Atract_rpb*$k*$Temp/$Atract_tormod )*180/$PI,
        dz         => 3.4,
    );
    $Atract{min} = 4;       # minimum number of A's in Atract
    $Atract{max} = 6;       # maximum number of A's in Atract

    # Default Atract values for 5' junction of Atract
    %Atract_5 = (
        tilt       => -7.7,
        roll       => 4.6,   
        twist      => 360/$Atract_helical_repeat,
        tilt_flex  => 4.842,
        roll_flex  => 4.842,
        twist_flex => sqrt( $Atract_rpb*$k*$Temp/$Atract_tormod )*180/$PI,
        dz         => 3.4,
    );

    # Default Atract values for 3' junction
    %Atract_3 = (
        tilt       => 9.7,
        roll       => 0,
        twist      => 360/$Atract_helical_repeat,
        tilt_flex  => 4.842,
        roll_flex  => 4.842,
        twist_flex => sqrt( $Atract_rpb*$k*$Temp/$Atract_tormod )*180/$PI,
        dz         => 3.4,
    );
}