Analyzing a scattering profile

One of the most common tasks is analyzing a scattering profile. Here’s an example of how that might be done, from Guinier analysis through creating a 3D reconstruction.

Guinier fit and MW

import os
import shutil
import numpy as np
import bioxtasraw.RAWAPI as raw


#Load the settings
settings = raw.load_settings('./standards_data/SAXS.cfg')


#Load the profile of interest
profile_names = ['./reconstruction_data/glucose_isomerase.dat']
profiles = raw.load_profiles(profile_names)

gi_prof = profiles[0]


#Automatically calculate the Guinier range and fit
(rg, i0, rg_err, i0_err, qmin, qmax, qrg_min, qrg_max, idx_min, idx_max,
    r_sq) = raw.auto_guinier(gi_prof, settings=settings)


#Calculate M.W. using 6 different methods
mw_ref = raw.mw_ref(gi_prof, 0.47, settings=settings)
mw_abs = raw.mw_abs(gi_prof, 0.47, settings=settings)
mw_vp, pvol_cor, pvol, vp_qmax = raw.mw_vp(gi_prof, settings=settings)
mw_vc, vcor, mw_err, vc_qmax = raw.mw_vc(gi_prof, settings=settings)
mw_bayes, mw_prob, ci_lower, ci_upper, ci_prob = raw.mw_bayes(gi_prof)
mw_datclass, shape, dmax_datclass = raw.mw_datclass(gi_prof)

Calculate IFTs

#Calculate the IFT using BIFT
(gi_bift, gi_bift_dmax, gi_bift_rg, gi_bift_i0, gi_bift_dmax_err,
    gi_bift_rg_err, gi_bift_i0_err, gi_bift_chi_sq, gi_bift_log_alpha,
    gi_bift_log_alpha_err, gi_bift_evidence,
    gi_bift_evidence_err) = raw.bift(gi_prof, settings=settings)


#Calculate the IFT using GNOM
(gi_datgnom_ift, gi_datgnom_dmax, gi_datgnom_rg, gi_datgnom_i0,
    gi_datgnom_rg_err, gi_datgnom_i0_err, gi_datgnom_total_est,
    gi_datgnom_chi_sq, gi_datgnom_alpha, gi_datgnom_quality) = raw.datgnom(gi_prof)

Save the profile and IFTs

#Save the profile and IFTs
if not os.path.exists('./api_results'):
    os.mkdir('./api_results')

raw.save_profile(gi_prof, 'gi.dat', './api_results')
raw.save_ift(gi_datgnom_ift, 'gi.out', './api_results')
raw.save_ift(gi_bift, 'gi.ift', './api_results')


#Calculate the ambiguity using AMBIMETER
a_score, a_cats, a_eval = raw.ambimeter(gi_datgnom_ift)

Create a bead model reconstruction

#Create individual bead model reconstructions
if not os.path.exists('./api_results/gi_dammif'):
    os.mkdir('./api_results/gi_dammif')
else:
    files = os.listdir('./api_results/gi_dammif')
    for f in files:
        os.remove(os.path.join('./api_results/gi_dammif', f))

chi_sq_vals = []
rg_vals = []
dmax_vals = []
mw_vals = []
ev_vals = []

for i in range(5):
    chi_sq, rg, dmax, mw, ev = raw.dammif(gi_datgnom_ift,
        'gi_{:02d}'.format(i+1), './api_results/gi_dammif', mode='Fast')

    chi_sq_vals.append(chi_sq)
    rg_vals.append(rg)
    dmax_vals.append(dmax)
    mw_vals.append(mw)
    ev_vals.append(ev)


#Average the bead model reconstructions
damaver_files = ['gi_{:02d}-1.pdb'.format(i+1) for i in range(5)]

(mean_nsd, stdev_nsd, rep_model, result_dict, res, res_err,
    res_unit) = raw.damaver(damaver_files, 'gi',
    './api_results/gi_dammif')


#Cluster the bead model reconstructions
cluster_list, distance_list = raw.damclust(damaver_files, 'gi',
    './api_results/gi_dammif')


#Refine the bead model
chi_sq, rg, dmax, mw, ev = raw.dammin(gi_datgnom_ift, 'refine_gi',
    './api_results/gi_dammif', 'Refine',
    initial_dam='gi_damstart.pdb')


#Align the refined bead model with a high resolution structure
shutil.copy('./reconstruction_data/gi_complete/1XIB_4mer.pdb',
    './api_results/gi_dammif')

raw.supcomb('refine_gi-1.pdb', '1XIB_4mer.pdb', './api_results/gi_dammif')

Create an electron density reconstruction

#Create individual electron density reconstructions
if not os.path.exists('./api_results/gi_denss'):
    os.mkdir('./api_results/gi_denss')
else:
    files = os.listdir('./api_results/gi_denss')
    for f in files:
        os.remove(os.path.join('./api_results/gi_denss', f))

rhos = []
chi_vals = []
rg_vals = []
support_vol_vals = []
sides = []
fit_data = []

for i in range(5):
    (rho, chi_sq, rg, support_vol, side, q_fit, I_fit, I_extrap,
        err_extrap) = raw.denss(gi_datgnom_ift, 'gi',
        './api_results/gi_denss', mode='Fast')

    rhos.append(rho)
    chi_vals.append(chi_sq)
    rg_vals.append(rg)
    support_vol_vals.append(support_vol)
    sides.append(side)
    fit_data.append([q_fit, I_fit, I_extrap, err_extrap])


#Average the electron reconstructions
(average_rho, mean_cor, std_cor, threshold, res, scores,
    fsc) = raw.denss_average(np.array(rhos), side, 'gi_average',
    './api_results/gi_denss')


#Refine the electron density
(refined_rho, refined_chi_sq, refined_rg, refined_support_vol, refined_side,
    refined_q_fit, refined_I_fit, refined_I_extrap,
    refined_err_extrap) = raw.denss(gi_datgnom_ift, 'gi_refine',
    './api_results/gi_denss', mode='Fast',
    initial_model=average_rho)


#Align the electron density with a high resolution structure
shutil.copy('./reconstruction_data/gi_complete/1XIB_4mer.pdb',
    './api_results/gi_denss')

aligned_density, score = raw.denss_align(refined_rho, refined_side,
     '1XIB_4mer.pdb', './api_results/gi_denss',
     'gi_refined_aligned', './api_results/gi_denss')