ToolboxΒΆ
This tutorial shows the user what biceps.toolbox has to offer.
[1]:
import numpy as np
import biceps
BICePs - Bayesian Inference of Conformational Populations, Version 2.0
Warning on use of the timeseries module: If the inherent timescales of the system are long compared to those being analyzed, this statistical inefficiency may be an underestimate. The estimate presumes the use of many statistically independent samples. Tests should be performed to assess whether this condition is satisfied. Be cautious in the interpretation of the data.
Convert populations to energy
[2]:
pop_filename = 'population.txt'
# convert from population to energy(kT): E = -kT*log(P)
out = biceps.toolbox.convert_pop_to_energy(pop_filename)
print(out[::100]) # printing every 100
[5.583336321580334, 9.721165995742174, 9.115030192171858, 9.433483923290392, 6.495645628155304, 11.107460356862065, 6.094605306027495, 10.126631103850338, 8.8387768155437, 9.721165995742174]
Compute J coupling constants for natural amino acids
[3]:
# Compute J3_HN_HA for frames in trajectories
# for only 1 *gro/*pdb file
J = biceps.toolbox.get_J3_HN_HA('J3/example.gro',
model = "Habeck", outname='J3_out.npy')
# for 1 trajectory
J = biceps.toolbox.get_J3_HN_HA(traj='J3/example.xtc',top='J3/example.gro',
model = "Habeck", outname='J3_out.npy')
# for 1 trajectory with selected frames
frames = [0,1,3,5,7,9]
J = biceps.toolbox.get_J3_HN_HA(traj='J3/example.xtc',top='J3/example.gro',
frame = frames, model = "Habeck", outname='J3_out.npy')
saving output file...
Done!
saving output file...
Done!
saving output file...
Done!
/Users/RR/miniconda3/lib/python3.8/site-packages/numpy/lib/npyio.py:521: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
arr = np.asanyarray(arr)
Compute J coupling constants for non-natural amino acids
[5]:
biceps.toolbox.mkdir("J_coupling")
# compute J coupling constants of non-natural amino acids
# Karplus relations need to be determined
index = np.load('ind.npy') # atom indices
karplus_key = np.load('Karplus.npy').astype(str) # Karplus relation for each J coupling constant
for i in range(100):
J = biceps.toolbox.compute_nonaa_Jcoupling('top/%d.pdb'%i, indices=index, karplus_key=karplus_key)
if i<5:
print(J) # print the first 5
np.savetxt('J_coupling/%d.txt'%i,J)
[[11.59929633 11.5908202 3.77981085 3.06534659 11.46025735 0.61495801
10.80616333 2.94538664 1.24030298 2.77752424 1.3801523 10.88224604]]
[[11.59359379 11.58796156 11.5994007 10.86025507 4.40292064 12.30620068
3.43600049 10.7195684 1.07170021 2.98419957 10.8572357 1.18375748]]
[[11.58811557 11.58931135 10.49392263 4.02871242 11.33110466 12.59204923
10.9224119 2.64105364 1.43223646 2.05262705 10.99946341 1.99195757]]
[[11.59995984 11.59170756 11.34322617 10.98300376 2.75242845 0.56631515
1.14966803 10.75983707 3.06268945 2.78762916 1.35818234 10.85199877]]
[[11.58991007 11.58508194 10.75085679 3.20760985 11.53892434 0.03746871
10.9185836 2.76861831 1.46725946 2.17715562 10.99202389 1.80701938]]
Tools to be used after BICePs sampling
Convert ``biceps`` trajectory file to Pandas DataFrame
[6]:
file = "results/traj_lambda0.00.npz"
df = biceps.toolbox.npz_to_DataFrame(file)
print(df)
step E accept state para_index
0 0 135.228336 1 25 [[151], [116, 81]]
1 100 43.798846 1 58 [[170], [108, 92]]
2 200 35.797135 0 84 [[172], [110, 94]]
3 300 24.634310 0 79 [[176], [114, 95]]
4 400 19.433472 0 79 [[183], [113, 97]]
5 500 12.780269 0 87 [[187], [112, 98]]
6 600 8.845128 1 87 [[196], [117, 98]]
7 700 9.259623 0 87 [[194], [118, 100]]
8 800 10.022916 0 87 [[190], [124, 95]]
9 900 10.252859 1 87 [[193], [117, 96]]
# NOTE: The following cell is for pretty notebook rendering
[7]:
from IPython.core.display import HTML
def css_styling():
styles = open("../../../theme.css", "r").read()
return HTML(styles)
css_styling()
[7]: