Phenix/CCTBX quick reference

This section summarizes some commonly-used functions in Phenix and cctbx.

Setup

Copy tutorial data

phenix.setup_tutorial tutorial_name=model-building-scripting

Set up a DataManager

from iotbx.data_manager import DataManager    # Load in the DataManager
dm = DataManager()             # Initialize the DataManager and call it dm
dm.set_overwrite(True)         # Overwrite files with the same name

Show all the keywords and summary documentation for a function

help(dm)
help(dm.get_map_model_manager)

Reading/writing files and creating managers

Get a map_model manager from a map and model

mmm_ligand = dm.get_map_model_manager(        # get map_model manager
  model_file="boxed_ligand_map.pdb",   # model file
  map_files="boxed_ligand_map.ccp4")   # map file

Read some models

model_for_morphing = dm.get_model("short_model_box_for_morph.pdb") # model
model_for_sequence = dm.get_model("short_model_main_chain_box.pdb") # model

Make a copy of a model so that when model is shifted we have the original model

original_model_for_morphing = model_for_morphing.deep_copy()  # copy of model

Read in a model of atp

model_atp = dm.get_model("atp.pdb")  # model file

Read in a restraints file

dm.process_restraint_file("atp.cif")  # Set up restraints
restraints_atp=dm.get_restraint("atp.cif") # get a restraints object

Read in map and get map_managers

mm = dm.get_real_map("boxed_ligand_map.ccp4")   # map file
mm_for_morphing = dm.get_real_map("short_model_box.ccp4")   # another map file
dm.remove_real_map("short_model_box.ccp4")   # remove map from data_manager

Get a map_model manager from a map only

mmm_map_only = dm.get_map_model_manager(        # get map_model manager
    map_files="short_model_box.ccp4")   # map file

Set defaults in a map_manager

mm.set_wrapping(False)   # do not wrap (do not extend outside supplied map)
mm.set_experiment_type("xray")   # alternatives are cryo_em or neutron
mm.set_scattering_table("n_gaussian")   # electron/neutron/wk1995/it1992

Read in reconstruction symmetry (“ncs”) and make an ncs object:

ncs_d7 = dm.get_ncs_spec("D7.ncs_spec")   # reconstruction symmetry file

Make a map_model_manager

map_manager = mm_for_morphing.deep_copy()  # use a separate copy
model = model_for_morphing.deep_copy()    # model, also a separate copy
from iotbx.map_model_manager import map_model_manager # import manager
mmm = map_model_manager(map_manager = map_manager, model = model)   #

Make a map_model_manager with reconstruction symmetry

map_manager = mm_for_morphing.deep_copy()  # use a separate copy
model = model_for_morphing.deep_copy()    # model, also a separate copy
from iotbx.map_model_manager import map_model_manager # import manager
mmm_ncs = map_model_manager(map_manager = map_manager, model = model, #
 ncs_object = ncs_d7)   # reconstruction symmetry

Shift a model to match origin shift in a map_manager

model_1 = model.deep_copy()  # copy of a model
map_manager.shift_model_to_match_map(model_1)  # shift model_1 to match map

Make a copy of a map_model_manager (deep copy where everything is separate from origin)

mmm_copy = mmm.deep_copy() # deep copy of a map-model manager

Set defaults in a map_model_manager

mmm.set_resolution(3)   # set default resolution to 3 A
mmm.set_experiment_type("xray")   # alternatives are cryo_em/neutron

Generate a model map using a map_model_manager

dm = DataManager()             # Initialize the DataManager and call it dm
dm.set_overwrite(True)
model = dm.get_model("structure_search_target.pdb")    # read in a model
model_mmm = dm.get_map_model_manager(  # getting a map_model_manager
    model_file = "structure_search_target.pdb",)   # model to read in
model_mmm.generate_map(d_min=3)        #   generate map to resolution of 3 A

Converting from one type of manager to another

Get a model_building object from a map_model_manager

build = mmm.model_building()   # get model_building object

Get a map_model_manager from a model_building object

new_mmm = build.as_map_model_manager() # get map-model manager

Get a map_manager or model from a model_building object or map_model_manager

mm = build.map_manager()   # map_manager from model-building object
model = build.model()      # model from model-building object
mm = mmm.map_manager()     # map_manager from map-model manager
model = mmm.model()        # model from map-model manager

Working with map_managers

Tell DataManager to forget any previously-read map from a file, then read in map_manager from that file.

for file_name in dm.get_real_map_names():   # list of previously read maps
  dm.remove_real_map(file_name)   # forget previous reads
mm = dm.get_real_map("boxed_ligand_map.ccp4")   # read in map file

Deep-copy map_manager

copy_of_mm = mm.deep_copy()  #  Separate copy of map_data

Shift origin of map_manager object to place start of map data that is present at (0, 0, 0). Note: This is done automatically by map_model_manager if one is used.

mm.shift_origin() # shift origin to (0,0,0)

Shift origin of map_manager object back where it was originally

mm1 = copy_of_mm.deep_copy()   # map manager to work with
mm1.shift_origin_to_match_original()  # put origin back where it was

List all methods available for a map_manager

dir(mm)  # list all the methods available for mm

Set resolution, experiment_type, scattering_table

mm.set_resolution(3)   # nominal resolution is 3 A;
mm.set_experiment_type("xray")  # alternatives are cryo_em/neutron
mm.set_scattering_table("n_gaussian")   # electron/neutron/wk1995/it1992

Get origin shift in A (shift of position since map was read in)

shift_cart = mm.shift_cart() #   shift since start

Get position of original origin in grid units (same as how much to shift map to overlay original)

origin_shift_grid_units = mm.origin_shift_grid_units #  origin pos

Get gridding of working map (the part that is present)

n_real = mm.map_data().all()  #  gridding of working map

Get wrapping of map (wrapped means infinite, repeating with gridding of map. Map must be full size (not boxed)

wrapping = mm.wrapping()  #  wrapping of working map

Determine if map is full size, if map is zero-based, if map is consistent with wrapping

is_full_size = mm.is_full_size()  #  is map full size (not boxed)
is_zero_based = mm.origin_is_zero()  #  is zero-based
is_consistent_with_wrapping = mm.is_consistent_with_wrapping() # ok to wrap

Determine if map_manager is similar to another one (same gridding, origin)

is_similar = mm.is_similar(mm1)  # map_managers are similar

Get map_data (flex array) from map_manager object

map_data = mm.map_data()  #  map_data

Customized copy of map_manager object with new data

some_map_data = mm.map_data() + 5 #  some map_data
new_mm = mm.customized_copy(map_data = some_map_data)  # new map manager

Set values of data in existing map_manager or zero out a map_manager

new_mm.set_map_data(map_data = map_data)  # new_mm data is now map_data
                # Note:  change new_mm and you are also changing map_data
new_mm.initialize_map_data()  # set values of map_data in new_mm to zero

Get full size version of map (zeros outside existing region)

full_size = mm.as_full_size_map()  # full size version padded with zero

Resample map on new grid (new grid must be multiple of existing if origin_shift_grid_units is not (0,0,0)

n_real = mm.map_data().all()  # current gridding
new_n_real = [n_real[0],n_real[1],n_real[2] * 2]  # double gridding along c
resampled = mm.resample_on_different_grid(new_n_real)  # resample

Get correlation to another map

cc = mm.cc_to_other_map_manager(new_mm) # map correlation to other map

Resolution filter, gaussian filter, binary filter:

low_res_map =  mm.resolution_filter(10)  # 10 A map
gaussian_blur = mm.gaussian_filter(1) # blur with 1 A Gaussian
binary_filtered = mm.binary_filter(threshold = 0.5)  # set value to 1 if
                                # average of 27 in box around this point > 0.5

Create binary mask around density, around edges, around atoms (does not apply mask, just creates it)

mm.create_mask_around_density(resolution = 3) # mask around density
mm.create_mask_around_edges(soft_mask_radius = 3) # mask around edges
model = dm.get_model("structure_search_target.pdb")  # get a model
mm.create_mask_around_atoms(model = model)  # radius=max(resolution,3)

Make the mask a soft mask

mm.soft_mask()  #  make the mask a soft mask. Default radius is resolution

Apply mask to working map or get mask as a map_manager

mm.apply_mask()  # apply the mask to working map (changes working map)
mask_as_map_manager = mm.get_mask_as_map_manager()   # mask as a map manager

Get map values at coordinates

sites_cart = model.get_sites_cart()  #  some coordinates
density_values = mm.density_at_sites_cart(sites_cart)  # density at coords

Normalize the map (set mean to zero and sd to 1)

mm.set_mean_zero_sd_one()  #  normalize map

Find reconstruction symmetry in map

mm.find_map_symmetry() #  find map symmetry
ncs_obj = mm.ncs_object()  # get the map symmetry as an ncs object

Find highest grid points in map and return as coordinates. Note: sites are relative to origin of map

sites_cart = mm.find_n_highest_grid_points_as_sites_cart() # highest points

Find n_atoms grid points in map in high density separated by dist_min

sites_cart = mm.trace_atoms_in_map(dist_min = 2, n_atoms = 31) # separated
                             # sites in density

Get map as Fourier coefficients and back. Note: origin stays at (0,0,0) throughout.

map_coeffs = mm.map_as_fourier_coefficients(d_min = 3)  # map coeffs to 3 A
new_mm = mm.fourier_coefficients_as_map_manager(map_coeffs) # new map_manager

Working with map_data arrays (3D arrays)

Deep-copy map_data

copy_of_map_data = map_data.deep_copy()  #  Separate copy of map_data

Size of map_data

print ("Size of almost anything is anything.size(): %s" %(map_data.size()))

Pointer-only copy of map_data

pointer_to_map_data = map_data  #  just an alias to map_data

Get or set value of map_data at indices i,j,k

value = map_data[3,4,5] # get value of map_data at indices 3,4,5
map_data[3,4,5] = value # set value of map_data at indices 3,4,5

Select all map_data > 0. Note: selection array has same shape as map_data

sel = (map_data > 0)  #  array with True at all grid points with map_data > 0

Set all selected values to value. Note: if map_data is a 3D array, sel must also be a 3D array

map_data.set_selected(sel, value) # set values where sel is True to value

Work with 3D map as a 1D array (Note: the data in the arrays are identical; change one and the other changes too

map_data_1d = map_data.as_1d()  # 1-D array containing the data in map_data

Empty 1D array

from scitbx.array_family import flex   # import flex
new_array = flex.double(map_data_1d.size(), 0) # new array of zeroes

Convert a 1D array to 3D

acc = map_data.accessor()  # get the accessor for our 3d map data
new_array.reshape(acc)     # change new_array shape to match map_data

Working with 1-D versions of map_data arrays (flex.double)

Duplicate (deep copy, new numbers)

dc =  map_data_1d.deep_copy()    # copy of map_data_1d

Multiply, divide, add, check result

dc = dc * 3   # multiply all values in dc by 3 and make it the new dc
dc *= 3   # multiply all values in dc by 3 and make it the new dc
new_array = dc + (0.5 * dc ) * (dc)  # each value in new_array is
                                     #  value in dc + 0.5 *value in dc squared
new_array_2 = dc + 0.5 * flex.pow2(dc) # same result
assert new_array == new_array_2      # assert all values in 2 arrays match

Get the minimum, maximum, and mean

m = map_data_1d.min_max_mean()
print ("Min, max, mean: (%.3f, %.3f, %.3f" %(m.min, m.max, m.mean))

Get the standard deviation

sd = map_data_1d.standard_deviation_of_the_sample()

Working with model objects

Get coordinates (orthogonal Angstrom units, flex array) from model object

sites_cart = model.get_sites_cart()  # coordinates in model. Order as PDB file

Get coordinates (fractional, flex array) from model object

sites_frac = model.get_sites_frac()  # fractional coordinates

Convert fractional to orthogonal A and back

sites_frac = model.crystal_symmetry().unit_cell().fractionalize(sites_cart)  #
sites_cart = model.crystal_symmetry().unit_cell().orthogonalize(sites_frac)  #

Print out first 10 coordinates

for x,y,z in sites_cart[:10]:                 # go through sites_cart
  print ("SITE: (%.3f, %.3f, %.3f)" %(x,y,z))   # print them out

Select part of a model

model = dm.get_model("structure_search_target.pdb")   # read in model
part_of_model = model.apply_selection_string('resseq 6:25')  # select part

Separately create a selection (a flex.bool array) and then apply it

sel = model.selection('resseq 6:25')    # get a selection array
part_of_model = model.select(sel)       # apply the selection

Working with coordinate arrays (3-D, flex.vec3_double)

Get the minimum x,y,z or maximum x,y,z

min_xyz = sites_cart.min() # lowest value of x, of y, and of z (3 numbers)
max_xyz = sites_cart.max() # highest value of x, of y, and of z (3 numbers)

Get the distances of a set of coordinates from the origin

distances = sites_cart.norms()  # distances from origin

Get the distances between 2 sets of coordinates

sites_cart_2 = sites_cart.reversed()  # just something different
diffs = (sites_cart - sites_cart_2)  # difference vectors
distances = diffs.norms()  # distances between sites_cart and sites_cart_2

Get rms length of a set of coordinates

rms_length = sites_cart.rms_length()  # rms length

Get rms difference between two sets of coordinates

rms_difference = sites_cart.rms_difference(sites_cart_2)  # rms diff

Select all the sites where x > 10 and add 3 to their y values.

x = sites_cart.parts()[0]     # x values for each coordinate
y = sites_cart.parts()[1]     # y values
z = sites_cart.parts()[2]     # z values
sel = (x > 10)   #   all sites where x value is > 10

sel_y = y.select(sel)  # selected y values
new_y = y.set_selected(sel,sel_y + 3) # add 3 to selected y values
                              # Note: sel_y.size() == sel.count(True)
from scitbx.array_family import flex
new_sites_cart = flex.vec3_double(x, new_y, z) #  new sites with new y

Working with ncs (reconstruction symmetry) objects

Read in reconstruction symmetry (“ncs”) and make an ncs object:

ncs_d7 = dm.get_ncs_spec("D7.ncs_spec")   # reconstruction symmetry file

Print out ncs information

ncs_d7.display_all()  #  print out ncs

Get ncs from a model

ncs_d7_model = dm.get_model("D7.pdb")   # model with D7 symmetry
ncs_d7_model.search_for_ncs()  # find NCS in this model
ncs_spec_object = ncs_d7_model.get_ncs_obj().get_ncs_info_as_spec() # ncs
print (ncs_spec_object.display_all()) #   print out summary

Find ncs from a map and apply to a model using the model_building tool

mmm = dm.get_map_model_manager(model_file="D7.pdb",   # model
   map_files="D7.ccp4")    #map
build=mmm.model_building()    # get a model-building object
build.map_manager().find_map_symmetry(symmetry = "D7")  # find symmetry in map
build.map_manager().ncs_object().display_all() # summary
model_no_ncs = build.model().apply_selection_string("chain A") #
print (model_no_ncs)#  make sure we got one chain
model_ncs = build.apply_ncs(model = model_no_ncs) # apply ncs to model_no_ncs
print (model_ncs)  # now we have 14 chains