1 from __future__
import print_function
2 from IMP
import ArgumentParser
5 __doc__ =
"Perform analysis to determine sampling convergence."
15 parser = ArgumentParser(
16 description=
"First stages of analysis for assessing sampling "
18 parser.add_argument(
'--sysname',
'-n', dest=
"sysname",
19 help=
'name of the system', default=
"")
20 parser.add_argument(
'--path',
'-p', dest=
"path",
21 help=
'path to the good-scoring models', default=
"./")
22 parser.add_argument(
'--extension',
'-e', dest=
"extension",
23 help=
'extension of the file', choices=[
'rmf',
'pdb'],
25 parser.add_argument(
'--mode',
'-m', dest=
"mode", help=
'pyRMSD calculator',
26 choices=[
'cuda',
'cpu_omp',
'cpu_serial'], default=
"cuda")
27 parser.add_argument(
'--cores',
'-c', dest=
"cores", type=int,
28 help=
'number of cores for RMSD matrix calculations; '
29 'only for cpu_omp', default=1)
30 parser.add_argument(
'--subunit',
'-su', dest=
"subunit",
31 help=
'calculate RMSD/sampling and cluster precision/densities '
32 'etc over this subunit only', default=
None)
33 parser.add_argument(
'--align',
'-a', dest=
"align",
34 help=
'boolean flag to allow superposition of models',
35 default=
False, action=
'store_true')
36 parser.add_argument(
'--ambiguity',
'-amb', dest=
"symmetry_groups",
37 help=
'file containing symmetry groups', default=
None)
38 parser.add_argument(
'--scoreA',
'-sa', dest=
"scoreA",
39 help=
'name of the file having the good-scoring scores for sample A',
40 default=
"scoresA.txt")
41 parser.add_argument(
'--scoreB',
'-sb', dest=
"scoreB",
42 help=
'name of the file having the good-scoring scores for sample B',
43 default=
"scoresB.txt")
44 parser.add_argument(
'--rmfA',
'-ra', dest=
"rmf_A",
45 help=
'RMF file with conformations from Sample A', default=
None)
46 parser.add_argument(
'--rmfB',
'-rb', dest=
"rmf_B",
47 help=
'RMF file with conformations from Sample B', default=
None)
48 parser.add_argument(
'--gridsize',
'-g', dest=
"gridsize", type=float,
49 help=
'grid size for calculating sampling precision', default=10.0)
50 parser.add_argument(
'--skip',
'-s', dest=
"skip_sampling_precision",
51 help=
"This option will bypass the calculation of sampling "
52 "precision. This option needs to be used with the clustering "
53 "threshold option. Otherwise by default, sampling precision "
54 "is calculated and the clustering threshold is the "
55 "calculated sampling precision.", default=
False,
57 parser.add_argument(
'--cluster_threshold',
'-ct', dest=
"cluster_threshold",
59 help=
'final clustering threshold to visualize clusters. Assumes '
60 'that the user has previously calculated sampling precision '
61 'and wants clusters defined at a threshold higher than the '
62 'sampling precision for ease of analysis (lesser number of '
63 'clusters).', default=30.0)
64 parser.add_argument(
'--voxel',
'-v', dest=
"voxel", type=float,
65 help=
'voxel size for the localization densities', default=5.0)
66 parser.add_argument(
'--density_threshold',
'-dt', type=float,
67 dest=
"density_threshold",
68 help=
'threshold for localization densities', default=20.0)
69 parser.add_argument(
'--density',
'-d', dest=
"density",
70 help=
'file containing dictionary of density custom ranges',
72 parser.add_argument(
'--gnuplot',
'-gp', dest=
"gnuplot",
73 help=
"plotting automatically with gnuplot", default=
False,
75 return parser.parse_args()
78 def make_cluster_centroid(infname, frame, outfname, cluster_index,
79 cluster_size, precision, density, path):
83 if hasattr(RMF.NodeHandle,
'replace_child'):
84 print(infname, outfname)
85 inr = RMF.open_rmf_file_read_only(infname)
86 outr = RMF.create_rmf_file(outfname)
87 cpf = RMF.ClusterProvenanceFactory(outr)
88 RMF.clone_file_info(inr, outr)
89 RMF.clone_hierarchy(inr, outr)
90 RMF.clone_static_frame(inr, outr)
91 inr.set_current_frame(RMF.FrameID(frame))
93 RMF.clone_loaded_frame(inr, outr)
94 rn = outr.get_root_node()
95 children = rn.get_children()
96 if len(children) == 0:
99 prov = [c
for c
in rn.get_children()
if c.get_type() == RMF.PROVENANCE]
104 newp = rn.replace_child(prov,
"cluster.%d" % cluster_index,
107 cp.set_members(cluster_size)
108 cp.set_precision(precision)
109 cp.set_density(os.path.abspath(density))
113 print(infname, frame, outfname)
114 subprocess.call([
'rmf_slice', path + infname,
'-f', str(frame),
126 from scipy
import spatial
129 from IMP.sampcon import scores_convergence, clustering_rmsd
130 from IMP.sampcon import rmsd_calculation, precision_rmsd
134 idfile_A =
"Identities_A.txt"
135 idfile_B =
"Identities_B.txt"
141 with open(os.path.join(args.path, args.scoreA),
'r') as f:
143 score_A.append(float(line.strip(
"\n")))
145 with open(os.path.join(args.path, args.scoreB),
'r') as f:
147 score_B.append(float(line.strip(
"\n")))
149 scores = score_A + score_B
152 scores_convergence.get_top_scorings_statistics(scores, 0, args.sysname)
155 scores_convergence.get_scores_distributions_KS_Stats(
156 score_A, score_B, 100, args.sysname)
159 if args.extension ==
"pdb":
162 conforms, masses, radii, models_name = \
163 rmsd_calculation.get_pdbs_coordinates(
164 args.path, idfile_A, idfile_B)
166 args.extension =
"rmf3"
168 if args.rmf_A
is not None:
169 (ps_names, masses, radii, conforms, symm_groups, models_name,
170 n_models) = rmsd_calculation.get_rmfs_coordinates_one_rmf(
171 args.path, args.rmf_A, args.rmf_B, args.subunit,
172 args.symmetry_groups)
177 (ps_names, masses, radii, conforms,
178 models_name) = rmsd_calculation.get_rmfs_coordinates(
179 args.path, idfile_A, idfile_B, args.subunit)
181 print(
"Size of conformation matrix", conforms.shape)
183 if not args.skip_sampling_precision:
186 numpy.save(
"conforms", conforms)
187 inner_data = rmsd_calculation.get_rmsds_matrix(
188 conforms, args.mode, args.align, args.cores, symm_groups)
189 print(
"Size of RMSD matrix (flattened):", inner_data.shape)
191 conforms = numpy.load(
"conforms.npy")
192 os.unlink(
'conforms.npy')
194 import pyRMSD.RMSDCalculator
195 from pyRMSD.matrixHandler
import MatrixHandler
196 mHandler = MatrixHandler()
197 mHandler.loadMatrix(
"Distances_Matrix.data")
199 rmsd_matrix = mHandler.getMatrix()
200 distmat = rmsd_matrix.get_data()
202 distmat_full = sp.spatial.distance.squareform(distmat)
203 print(
"Size of RMSD matrix (unpacked, N x N):", distmat_full.shape)
206 if args.rmf_A
is not None:
207 sampleA_all_models = list(range(n_models[0]))
208 sampleB_all_models = list(range(n_models[0], n_models[1] + n_models[0]))
209 total_num_models = n_models[1] + n_models[0]
212 sampleB_all_models) = clustering_rmsd.get_sample_identity(
214 total_num_models = len(sampleA_all_models) + len(sampleB_all_models)
215 all_models = list(sampleA_all_models) + list(sampleB_all_models)
216 print(
"Size of Sample A:", len(sampleA_all_models),
217 " ; Size of Sample B: ", len(sampleB_all_models),
218 "; Total", total_num_models)
220 if not args.skip_sampling_precision:
222 print(
"Calculating sampling precision")
226 gridSize = args.gridsize
229 cutoffs_list = clustering_rmsd.get_cutoffs_list(distmat, gridSize)
230 print(
"Clustering at thresholds:", cutoffs_list)
233 pvals, cvs, percents = clustering_rmsd.get_clusters(
234 cutoffs_list, distmat_full, all_models, total_num_models,
235 sampleA_all_models, sampleB_all_models, args.sysname)
239 (sampling_precision, pval_converged, cramersv_converged,
240 percent_converged) = clustering_rmsd.get_sampling_precision(
241 cutoffs_list, pvals, cvs, percents)
244 with open(
"%s.Sampling_Precision_Stats.txt"
245 % args.sysname,
'w+')
as fpv:
246 print(
"The sampling precision is defined as the largest allowed "
247 "RMSD between the cluster centroid and a ", args.sysname,
248 "model within any cluster in the finest clustering for which "
249 "each sample contributes models proportionally to its size "
250 "(considering both significance and magnitude of the "
251 "difference) and for which a sufficient proportion of all "
252 "models occur in sufficiently large clusters. The sampling "
253 "precision for our ", args.sysname,
254 " modeling is %.3f" % (sampling_precision),
" A.", file=fpv)
256 print(
"Sampling precision, P-value, Cramer's V and percentage "
257 "of clustered models below:", file=fpv)
258 print(
"%.3f\t%.3f\t%.3f\t%.3f"
259 % (sampling_precision, pval_converged, cramersv_converged,
260 percent_converged), file=fpv)
263 final_clustering_threshold = sampling_precision
266 final_clustering_threshold = args.cluster_threshold
269 print(
"Clustering at threshold %.3f" % final_clustering_threshold)
270 (cluster_centers, cluster_members) = clustering_rmsd.precision_cluster(
271 distmat_full, total_num_models, final_clustering_threshold)
273 (ctable, retained_clusters) = clustering_rmsd.get_contingency_table(
274 len(cluster_centers), cluster_members, all_models,
275 sampleA_all_models, sampleB_all_models)
276 print(
"Contingency table:", ctable)
278 with open(
"%s.Cluster_Population.txt" % args.sysname,
'w+')
as fcp:
279 for rows
in range(len(ctable)):
280 print(rows, ctable[rows][0], ctable[rows][1], file=fcp)
283 density_custom_ranges = precision_rmsd.parse_custom_ranges(args.density)
286 fpc = open(
"%s.Cluster_Precision.txt" % args.sysname,
'w+')
290 for i
in range(len(retained_clusters)):
291 clus = retained_clusters[i]
295 conform_0 = conforms[all_models[cluster_members[clus][0]]]
298 if not os.path.exists(
"./cluster.%s" %i):
299 os.mkdir(
"./cluster.%s" %i)
300 os.mkdir(
"./cluster.%s/Sample_A/" % i)
301 os.mkdir(
"./cluster.%s/Sample_B/" % i)
303 shutil.rmtree(
"./cluster.%s" %i)
304 os.mkdir(
"./cluster.%s" %i)
305 os.mkdir(
"./cluster.%s/Sample_A/" % i)
306 os.mkdir(
"./cluster.%s/Sample_B/" % i)
310 gmd1 = precision_rmsd.GetModelDensity(
311 custom_ranges=density_custom_ranges,
312 resolution=args.density_threshold, voxel=args.voxel,
314 gmd2 = precision_rmsd.GetModelDensity(
315 custom_ranges=density_custom_ranges,
316 resolution=args.density_threshold, voxel=args.voxel,
318 gmdt = precision_rmsd.GetModelDensity(
319 custom_ranges=density_custom_ranges,
320 resolution=args.density_threshold, voxel=args.voxel,
324 both_file = open(
'cluster.'+str(i)+
'.all.txt',
'w')
325 sampleA_file = open(
'cluster.'+str(i)+
'.sample_A.txt',
'w')
326 sampleB_file = open(
'cluster.'+str(i)+
'.sample_B.txt',
'w')
331 for pi
in range(len(conform_0)):
340 cluster_precision = 0.0
345 for mem
in cluster_members[clus]:
347 model_index = all_models[mem]
351 if args.symmetry_groups:
352 rmsd, superposed_ps, trans = \
353 precision_rmsd.get_particles_from_superposed_amb(
354 conforms[model_index], conform_0, args.align, ps, trans,
357 rmsd, superposed_ps, trans = \
358 precision_rmsd.get_particles_from_superposed(
359 conforms[model_index], conform_0, args.align, ps, trans)
363 cluster_precision += rmsd
366 gmdt.add_subunits_density(superposed_ps)
367 print(model_index, file=both_file)
369 if model_index
in sampleA_all_models:
371 gmd1.add_subunits_density(superposed_ps)
372 print(model_index, file=sampleA_file)
375 gmd2.add_subunits_density(superposed_ps)
376 print(model_index, file=sampleB_file)
378 cluster_precision /= float(len(cluster_members[clus]) - 1.0)
380 print(
"Cluster precision (average distance to cluster centroid) "
381 "of cluster ", str(i),
" is %.3f" % cluster_precision,
"A",
389 density = gmdt.write_mrc(path=
"./cluster.%s" % i, file_prefix =
"LPD")
390 gmd1.write_mrc(path=
"./cluster.%s/Sample_A/" % i, file_prefix =
"LPD")
391 gmd2.write_mrc(path=
"./cluster.%s/Sample_B/" % i, file_prefix =
"LPD")
394 cluster_center_index = cluster_members[clus][0]
395 if args.rmf_A
is not None:
396 cluster_center_model_id = cluster_center_index
397 if cluster_center_index < n_models[0]:
398 make_cluster_centroid(os.path.join(args.path, args.rmf_A),
399 cluster_center_index,
400 os.path.join(
"cluster.%d" % i,
401 "cluster_center_model.rmf3"),
402 i, len(cluster_members[clus]),
403 cluster_precision, density, args.path)
405 make_cluster_centroid(os.path.join(args.path, args.rmf_B),
406 cluster_center_index - n_models[0],
407 os.path.join(
"cluster.%d" % i,
408 "cluster_center_model.rmf3"),
409 i, len(cluster_members[clus]),
410 cluster_precision, density, args.path)
413 cluster_center_model_id = all_models[cluster_center_index]
414 outfname = os.path.join(
"cluster.%d" % i,
415 "cluster_center_model." + args.extension)
416 if 'rmf' in args.extension:
417 make_cluster_centroid(
418 models_name[cluster_center_model_id], 0, outfname,
419 i, len(cluster_members[clus]),
420 cluster_precision, density, args.path)
422 shutil.copy(models_name[cluster_center_model_id], outfname)
431 thisdir = os.path.dirname(__file__)
433 for filename
in sorted(glob.glob(os.path.join(gnuplotdir,
"*.plt"))):
434 cmd = [
'gnuplot',
'-e',
'sysname="%s"' % args.sysname, filename]
436 subprocess.check_call(cmd)
438 if __name__ ==
'__main__':
def get_data_path
Return the full path to one of this module's data files.
static XYZR setup_particle(Model *m, ParticleIndex pi)
static XYZ setup_particle(Model *m, ParticleIndex pi)
Class for storing model, its restraints, constraints, and particles.
static Mass setup_particle(Model *m, ParticleIndex pi, Float mass)
Class to handle individual particles of a Model object.
Sampling exhaustiveness protocol.