1 from __future__
import print_function
7 from optparse
import OptionParser
10 usage =
"""%prog [options] <data_fn> <n_components> <out_txt_fn>
11 Create a GMM from either density file (.mrc), a pdb file (.pdb)
12 Will detect input format from extension.
13 Outputs as text and optionally as a density map
16 parser = OptionParser(usage)
17 parser.add_option(
"-t",
"--covar_type",dest=
"covar_type",default=
'full',
18 choices=[
'spherical',
'tied',
'diag',
'full'],
19 help=
"covariance type for the GMM")
20 parser.add_option(
"-m",
"--out_map",dest=
"out_map",default=
'',
21 help=
"write out the gmm to an mrc file")
22 parser.add_option(
"-a",
"--apix",dest=
"apix",default=1.0,type=
'float',
23 help=
"if you don't provide a map, set the voxel_size here (for sampling)")
24 parser.add_option(
"-r",
"--resolution",dest=
"resolution",default=1,type=
'float',
25 help=
"if you don't provide a map, set the resolution here (for sampling)")
26 parser.add_option(
"-n",
"--num_samples",dest=
"num_samples",default=1000000,type=
'int',
27 help=
"num samples to draw from the density map")
28 parser.add_option(
"-i",
"--num_iter",dest=
"num_iter",default=100,type=
'int',
29 help=
"num iterations of GMM")
30 parser.add_option(
"-s",
"--threshold",dest=
"threshold",default=0.0,type=
'float',
31 help=
"threshold for the map before sampling")
33 parser.add_option(
"-c",
"--input_anchors_fn",dest=
"input_anchors_fn",default=
'',
34 help=
"get initial centers from anchors file. ")
35 parser.add_option(
"-f",
"--force_radii",dest=
"force_radii",default=-1.0,
37 help=
"force radii to be this value (spherical) -1 means deactivated ")
38 parser.add_option(
"-w",
"--force_weight",dest=
"force_weight",default=-1.0,
40 help=
"force weight to be this value (spherical) -1 means deactivated ")
41 parser.add_option(
"-e",
"--force_weight_frac",dest=
"force_weight_frac",action=
"store_true",default=
False,
42 help=
"force weight to be 1.0/(num anchors). takes precedence over -w ")
43 parser.add_option(
"-o",
"--out_anchors_txt",dest=
"out_anchors_txt",default=
'',
44 help=
"write final GMM as anchor points (txt)")
45 parser.add_option(
"-q",
"--out_anchors_cmm",dest=
"out_anchors_cmm",default=
'',
46 help=
"write final GMM as anchor points (cmm)")
47 parser.add_option(
"-d",
"--use_dirichlet",dest=
"use_dirichlet",default=
False,
49 help=
"use dirichlet process for fit")
51 parser.add_option(
"-k",
"--multiply_by_mass",dest=
"multiply_by_mass",default=
False,
53 help=
"if set, will multiply all weights by the total mass of the particles (PDB ONLY)")
55 parser.add_option(
"-z",
"--use_cpp",dest=
"use_cpp",default=
False,
57 help=
"EXPERIMENTAL. Uses the IMP GMM code. Requires isd_emxl")
60 (options, args) = parser.parse_args()
62 parser.error(
"incorrect number of arguments")
66 options,args=parse_args()
67 data_fn,ncenters,out_txt_fn=args
68 ncenters=int(ncenters)
71 if not os.path.isfile(data_fn):
72 raise Exception(
"The data file you entered: "+data_fn+
" does not exist!")
75 ext=data_fn.split(
'.')[-1]
84 if options.multiply_by_mass:
89 print(
'ERROR: data_fn extension must be pdb, mrc, or npy')
92 if not options.use_cpp:
93 print(
'sampling points')
94 pts=IMP.isd.sample_points_from_density(dmap,options.num_samples,options.threshold)
99 if options.force_weight_frac:
100 force_weight=1.0/ncenters
102 force_weight=options.force_weight
104 print(
'weight forced to',force_weight)
105 if not options.use_dirichlet:
107 options.num_iter,options.covar_type,
108 force_radii=options.force_radii,
109 force_weight=options.force_weight,
110 mass_multiplier=mass_multiplier)
113 options.num_iter,options.covar_type,
114 mass_multiplier=mass_multiplier)
120 print(
"This option is experimental, only works if you have isd_emxl")
122 density_ps=IMP.isd_emxl.fit_gaussians_to_density(mdl,dmap,options.num_samples,
123 ncenters,options.num_iter,
128 if options.out_map!=
'':
133 if options.out_anchors_txt!=
'':
134 IMP.isd.gmm_tools.write_gmm_to_anchors(density_ps,options.out_anchors_txt,
135 options.out_anchors_cmm)
139 if __name__==
"__main__":
Select non water and non hydrogen atoms.
double get_mass(ResidueType c)
Get the mass from the residue type.
GenericHierarchies get_leaves(Hierarchy mhd)
Get all the leaves of the bit of hierarchy.
void read_pdb(TextInput input, int model, Hierarchy h)
Class for storing model, its restraints, constraints, and particles.
Class for sampling a density map from particles.
Basic utilities for handling cryo-electron microscopy 3D density maps.
DensityMap * read_map(std::string filename)
Read a density map from a file and return it.
algebra::BoundingBoxD< 3 > get_bounding_box(const DensityMap *m)
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...