from __future__ import division from __future__ import with_statement from phaser import tbx_utils from phaser import output from libtbx import runtime_utils import libtbx.phil import os.path PROGRAM = "ensembler" VERSION = "0.2.0" MIN_SITE_COUNT = 3 MIN_MODEL_COUNT = 2 def find_chains(pdbs, logger): from phaser import mmt data = [] for f in pdbs: logger.info( msg = "Processing %s..." % f.name ) logger.indent() for info in f.get_chain_infos(): if not mmt.PROTEIN.recognize_chain_type( chain = info.chain ): logger.info( msg = "%s: not protein. Discarded" % info ) continue data.append( info ) logger.dedent() logger.info( msg = output.heading( text = "Chains found:" ) ) for ( index, info ) in enumerate( data, start = 1 ): message = "%d: %s" % ( index, info ) if info.identity is not None: message += ", identity wrt target: %.1f" % info.identity logger.info( msg = message ) return data def map_by_resid(infos, alignments, logger): logger.info( msg = "( resseq, icode ) will be used to map residues" ) logger.indent() # Do the mapping equiv_with = {} count = len( infos ) for ( index, info ) in enumerate( infos ): logger.info( "Mapping %s..." % info ) for rg in info.chain.residue_groups(): array = equiv_with.setdefault( ( rg.resseq_as_int(), rg.icode ), [ None ] * count ) array[ index ] = rg logger.dedent() # Convert it to array mapping mapping = [ equiv_with[ resid ] for resid in sorted( equiv_with ) ] return mapping def map_by_alignments(infos, alignments, logger): logger.info( msg = "Alignments will be used to map residues" ) assert infos assert all( alf.object.alignments for alf in alignments ) if not alignments: raise RuntimeError, "No alignments given" # Compare alignment first lines references = [ "".join( [ c for c in alf.object.alignments[0] if c != alf.object.gap ] ) for alf in alignments ] if not all( ali == references[0] for ali in references[1:] ): raise RuntimeError, "Unsuitable alignment set: different master sequences" logger.info( msg = "Master sequence:" ) logger.info( msg = references[0] ) # Do the mapping from phaser import mmt from phaser import rsam from iotbx import bioinformatics one_letter_for = dict( zip( mmt.PROTEIN.three_letter_codes, mmt.PROTEIN.one_letter_codes ) ) mapping = [] min_length = 2 for info in infos: logger.info( msg = "Mapping %s..." % info ) logger.indent() rgs = info.chain.residue_groups() for alf in alignments: ( index, rmap ) = rsam.map_to_alignment( alignment = alf.object, residues = rgs, consecutivity = mmt.PROTEIN.consecutivity, one_letter_for = one_letter_for, min_length = min_length ) if index < 0: continue master = alf.object.alignments[0] seq = alf.object.alignments[ index ] assert len( seq ) == len( rmap ), "len( %s ) != len( %s )" % ( seq, rmap ) logger.info( msg = "Matching sequence found in %s" % alf.name ) logger.info( msg = "Mapping to master:" ) logger.info( msg = bioinformatics.clustal_alignment( alignments = [ master, seq ], names = [ "MASTER", "CURRENT" ], ) ) log_residue_discard( present = rgs, mapping = rmap, logger = logger, min_length = min_length ) mapping.append( [ rg for ( c, rg ) in zip( master, rmap ) if c != alf.object.gap ] ) logger.dedent() break else: logger.dedent() raise RuntimeError, "No matching alignment for %s" % info # Reformatting mapping assert len( infos ) == len( mapping ) assert all( len( t ) == len( mapping[0] ) for t in mapping[1:] ) return zip( *mapping ) def map_by_ssm(infos, alignments, logger): logger.info( msg = "SSM will be used to map residues" ) assert infos logger.info( msg = "Reference: %s" % infos[0] ) reference = infos[0].chain reference_rgs = reference.residue_groups() mapping = [ reference_rgs ] indexer = lambda chain: dict( [ ( ( chain.id, rg.resseq_as_int(), rg.icode ), i ) for ( i, rg ) in enumerate( chain.residue_groups() ) ] ) import ccp4io_adaptbx logger.indent() for info in infos[1:]: logger.info( msg = "Mapping %s..." % info ) rgs = info.chain.residue_groups() try: ssm = ccp4io_adaptbx.SecondaryStructureMatching( reference = reference, moving = info.chain ) alignment = ccp4io_adaptbx.SSMAlignment( match = ssm, indexer = indexer ) except RuntimeError, e: raise RuntimeError, "SSM error: %s" % e m = [ None ] * len( reference_rgs ) for ( cur, ref ) in alignment.pairs: if ref is None or cur is None: continue assert ref < len( reference_rgs ) assert cur < len( rgs ) m[ ref ] = rgs[ cur ] mapping.append( m ) logger.dedent() # Reformatting mapping assert all( [ len( t ) == len( mapping[0] ) for t in mapping[1:] ] ) return zip( *mapping ) def map_by_multiple_alignment(infos, alignments, logger): logger.info( msg = "A multiple alignment will be used to map residues" ) if not alignments: raise RuntimeError, "No alignments given" if len( alignments ) != 1: logger.warning( msg = "Several alignment files specified; using first one only" ) return chain_mapping_via_alignment( infos = infos, alignment = alignments[0].object, logger = logger ) def map_by_muscle(infos, alignments, logger): logger.info( msg = "phenix.muscle will be used to map residues" ) sequences = [] from phaser import mmt from phaser import rsam from iotbx import bioinformatics one_letter_for = dict( zip( mmt.PROTEIN.three_letter_codes, mmt.PROTEIN.one_letter_codes ) ) for ( i, info ) in enumerate( infos, start = 1 ): seq = bioinformatics.sequence( name = "Chain_%d" % i, sequence = rsam.one_letter_sequence( rgs = info.chain.residue_groups(), one_letter_for = one_letter_for, unknown = mmt.PROTEIN.unknown_one_letter ) ) sequences.append( str( seq ) ) from mmtbx.msa import get_muscle_alignment alignment = get_muscle_alignment( fasta_sequences = "\n".join( sequences ) ) logger.info( msg = "Alignment:\n%s" % alignment ) logger.info( msg = "Chain assignment:" ) print_chain_infos( infos = infos, logger = logger ) return chain_mapping_via_alignment( infos = infos, alignment = alignment, logger = logger ) def chain_mapping_via_alignment(infos, alignment, logger): from phaser import mmt from phaser import rsam assert infos one_letter_for = dict( zip( mmt.PROTEIN.three_letter_codes, mmt.PROTEIN.one_letter_codes ) ) mapping = [] min_length = 2 logger.indent() for info in infos: logger.info( msg = "Mapping %s..." % info ) rgs = info.chain.residue_groups() ( index, rmap ) = rsam.map_to_alignment( alignment = alignment, residues = rgs, consecutivity = mmt.PROTEIN.consecutivity, one_letter_for = one_letter_for, min_length = min_length ) if index < 0: raise RuntimeError, "No matching alignment for %s" % info logger.info( msg = "Matching sequence found: index %d" % index ) log_residue_discard( present = rgs, mapping = rmap, logger = logger, min_length = min_length ) mapping.append( rmap ) logger.dedent() # Reformatting mapping assert all( [ len( t ) == len( mapping[0] ) for t in mapping[1:] ] ) return zip( *mapping ) def log_residue_discard(present, mapping, logger, min_length): substrahend = set( [ r for r in mapping if r ] ) unaccounted = [ rg for rg in present if rg not in substrahend ] if unaccounted: logger.info( msg = "%s residues in short (%s) segments discarded" % ( len( unaccounted ), "<%s residues" % min_length, ) ) MAPPING_PROCEDURE_FOR = { "resid": map_by_resid, "alignments": map_by_alignments, "multiple_alignment": map_by_multiple_alignment, "muscle": map_by_muscle, "ssm": map_by_ssm, } def print_chain_infos(infos, logger): logger.info( msg = output.blank() ) for ( index, info ) in enumerate( infos, start = 1 ): logger.info( "%d - %s" % ( index, info ) ) logger.info( msg = output.blank() ) def calculate_and_print_sequence_identity_statistics(infos, mapping, logger): from scitbx.array_family import flex count = len( infos ) assert 2 <= count overlaps = flex.int( [ 0 ] * ( count ** 2 ) ) matches = flex.int( [ 0 ] * ( count ** 2 ) ) for pos in mapping: assert len( pos ) == count for ( i_left, rg_left ) in enumerate( pos ): for ( i_right, rg_right ) in enumerate( pos ): if rg_left and rg_right: overlaps[ i_left * count + i_right ] += 1 if rg_left.unique_resnames()[0] == rg_right.unique_resnames()[0]: matches[ i_left * count + i_right ] += 1 # Checking that no zero overlap occurs for ( index, d ) in enumerate( overlaps ): if d == 0: left = int( index / count ) right = index - left * count logger.warning( msg = "No overlap between %s and %s" % ( infos[ left ], infos[ right ] ) ) assert matches[ index ] == 0 overlaps[ index ] = 1 identities = matches.as_double() / overlaps.as_double() logger.info( msg = "Pairwise sequence identities:" ) print_chain_infos( infos = infos, logger = logger ) reformatted = [ identities[ i_row * count : i_row * count + count ] for i_row in range( count ) ] logger.info( msg = tbx_utils.format_table( captions = [ " " ] + [ "%4d" % ( i + 1 ) for i in range( count ) ], formats = [ "%d" ] + [ "%4.2f" ] * count, lines = [ [ i ] + list( sl ) for ( i, sl ) in enumerate( reformatted, start = 1 ) ] ) ) logger.info( msg = output.blank() ) # Compensate for identity[i][i] == 1.0 averaged = [ ( flex.sum( sl ) - 1.0 ) / ( count - 1 ) for sl in reformatted ] logger.info( msg = "Averaged identities:" ) logger.info( msg = output.blank() ) for ( info, aver ) in zip( infos, averaged ): logger.info( msg = "%s: %5.2f" % ( info, aver ) ) return averaged def get_atom_selector(selection, logger): tokens = [ c.strip().upper() for c in selection.split(",") ] if "all" in tokens: logger.info( msg = "All atoms will be selected" ) selector = lambda atom: True else: logger.info( msg = "The following atoms will be selected: %s" % ", ".join( tokens ) ) selector = lambda atom: atom.name.strip() in tokens return selector def to_atom_mapping(mapping, selector, logger): assert mapping count = len( mapping[0] ) assert all( [ len( t ) == count for t in mapping[1:] ] ) atom_mapping = [] for ( position, equivalents ) in enumerate( mapping ): logger.debug( msg = "Mapping position %s..." % position ) atoms_named = {} for ( index, rg ) in enumerate( equivalents ): if not rg: continue known = set() for a in [ a for a in rg.atoms() if selector( atom = a ) ]: if a.name in known: logger.debug( msg = "Resid %s: multiple atoms %s" % ( rg.resid(), a.name ) ) continue atoms_named.setdefault( a.name, [ None ] * count )[ index ] = a # Append to mapping atom_mapping.extend( [ atoms_named[ n ] for n in sorted( atoms_named ) ] ) return atom_mapping def calculate_and_print_overlap_statistics(infos, mapping, logger): length = len( mapping ) count = len( infos ) assert 0 < count sites = [ 0 ] * count singles = [ 0 ] * count for pos in mapping: assert len( pos ) == count missing = pos.count( None ) assert missing < count is_single = ( missing == ( count - 1 ) ) for ( index, site ) in enumerate( pos ): if site: sites[ index ] += 1 if is_single: singles[ index ] += 1 overlaps = [ ( site - single ) / length for ( site, single ) in zip( sites, singles ) ] for ( info, site, single, overlap ) in zip( infos, sites, singles, overlaps ): logger.info( "%s: %s sites, %s singles (%.2f%% overlap)" % ( info, site, single, 100 * overlap, ) ) return overlaps def setup_unit_weighting_scheme(params): from phaser import multiple_superposition return multiple_superposition.UnitWeightScheme() def setup_robust_resistant_weighting_scheme(params): from phaser import multiple_superposition return multiple_superposition.RobustResistantWeightScheme( critical_value_square = params.robust_resistant.critical ) WEIGHTING_SCHEME_SETUP_FOR = { "unit": setup_unit_weighting_scheme, "robust_resistant": setup_robust_resistant_weighting_scheme, } def intersect_selection(mapping, logger): logger.info( msg = "Accepting positions that contain atoms for every chain" ) positions = [] for group in mapping: if None not in group: positions.append( group ) return positions def overlap_selection(mapping, logger): logger.info( msg = "Accepting positions that contain at least %s atoms" % MIN_MODEL_COUNT ) positions = [] for group in mapping: if MIN_MODEL_COUNT <= len( [ a for a in group if a ] ): positions.append( group ) return positions def diamond_superposition(sites, weights, logger): logger.info( msg = "Setting up Diamond superposition" ) from scitbx.array_family import flex from phaser import multiple_superposition return multiple_superposition.DiamondAlgorithm( site_sets = [ flex.vec3_double( [ a.xyz for a in m ] ) for m in zip( *sites ) ], weights = weights ) def wang_snoeyink_superposition(sites, weights, logger): logger.info( msg = "Setting up Wang-Snoeyink superposition" ) from scitbx.array_family import flex from phaser import multiple_superposition return multiple_superposition.WangSnoeyinkAlgorithm( site_sets = [ flex.vec3_double( [ a.xyz if a else ( 0, 0, 0 ) for a in m ] ) for m in zip( *sites ) ], selections = [ flex.bool( [ bool( a ) for a in m ] ) for m in zip( *sites ) ], weights = weights ) SUPERPOSITION_PROCEDURES_FOR = { "gapless": ( intersect_selection, diamond_superposition ), "gapped": ( overlap_selection, wang_snoeyink_superposition ), } def iterative_weighted_superposition( superposition, superposition_convergence, weighting, weights, weight_convergence, incremental_damping_factor, max_damping_factor, logger ): import math from scitbx.array_family import flex from phaser import multiple_superposition count = superposition.site_count while True: logger.info( msg = "Superposing..." ) try: ( rounds, residual ) = multiple_superposition.iterate( superposition = superposition, convergence = superposition_convergence ) except multiple_superposition.NoConvergenceException, e: logger.warning( msg = "Superposition: %s" % e.args[0] ) logger.warning( msg = "Continuing..." ) ( rounds, residual ) = e.args[1:3] logger.info( msg = "Residual: %s, %d iteration cycles" % ( residual, rounds ) ) diff_sq = superposition.distance_squares_from_average_structure() total_damping_factor = 1.0 while True: new_weights = weighting.for_difference_squares( squares = diff_sq ) try: superposition.prepare_new_weights( weights = new_weights ) except RuntimeError: logger.warning( msg = "Excessive weight shift, damping weight change" ) if max_damping_factor < total_damping_factor: raise RuntimeError, "Weight damping recovery exhausted" total_damping_factor *= incremental_damping_factor diff_sq = diff_sq / incremental_damping_factor continue break displacement = math.sqrt( flex.mean_sq( new_weights - weights ) ) logger.info( msg = "Weight change: %s" % displacement ) if displacement <= weight_convergence: break weights = new_weights superposition.set_new_weights() def calculate_rmsds(superposition, infos, logger): wrmsds = superposition.weighted_rmsds() unwrmsds = superposition.unweighted_rmsds() logger.info( msg = "Rmsds per chain:" ) print_chain_infos( infos = infos, logger = logger ) logger.info( msg = tbx_utils.format_table( captions = [ "Index", "Weighted rmsd", "Unweighted rmsd" ], formats = [ "%d", "%.3f", "%.3f" ], lines = zip( range( len( infos ) ), wrmsds, unwrmsds ) ) ) return ( wrmsds, unwrmsds ) def cluster_analysis(infos, superposition, critical, logger): count = superposition.set_count distances = [ [ 0.0 ] * count for i in range( count ) ] for l in range( count ): for r in range( l + 1, count ): distances[ l ][ r ] = superposition.rmsd_between( left = l, right = r ) distances[ r ][ l ] = distances[ l ][ r ] logger.info( msg = "Pairwise rms values" ) print_chain_infos( infos = infos, logger = logger ) logger.info( msg = tbx_utils.format_table( captions = [ " " ] + [ "%4d" % ( i + 1 ) for i in range( count ) ], formats = [ "%d" ] + [ "%4.2f" ] * count, lines = [ [ i ] + line for ( i, line ) in enumerate( distances, start = 1 ) ] ) ) logger.info( msg = output.heading( "Hierarchical clustering" ) ) from libtbx import cluster cl = cluster.HierarchicalClustering( data = range( count ), distance_function = lambda x, y: distances[ x ][ y ] ) levels = cl.getlevel( critical ) logger.info( msg = " %s clusters have been found" % len( levels ) ) print_chain_infos( infos = infos, logger = logger ) results = [ None ] * count for ( index, clust ) in enumerate( levels, start = 1 ): logger.info( " Cluster %s: %s" % ( index, ", ".join ( [ "%d" % ( s + 1 ) for s in sorted( clust ) ] ) ) ) for i in clust: assert results[ i ] is None results[ i ] = index return results def calculate_trimming_masks(infos, sites, deviations, threshold, logger): assert len( sites ) == len( deviations ) keeps = [ set() for i in infos ] logger.debug( msg = "Checking deviations of superposed positions..." ) for ( index, ( equivs, d ) ) in enumerate( zip( sites, deviations ) ): logger.debug( "Site %d: %.3f" % ( index, d ) ) if d <= threshold: logger.debug( "Accepting this position" ) assert len( equivs ) == len( keeps ) for ( a, k ) in zip( equivs, keeps ): if a: k.add( a.parent().parent() ) from scitbx.array_family import flex return [ flex.bool( [ rg in keep for rg in info.chain.residue_groups() ] ) for ( info, keep ) in zip( infos, keeps ) ] def trim_chains(infos, chains, logger): assert len( infos ) == len( chains ) from scitbx.array_family import flex for ( info, chain ) in zip( infos, chains ): mask = info.ensembler_trimming_mask assert len( chain.residue_groups() ) == len( mask ) logger.info( msg = "%s: removing %s residues" % ( info, flex.sum( ( ~mask ).as_int() ), ) ) for ( rg, to_remain ) in zip( chain.residue_groups(), mask ): if not to_remain: logger.debug( "Removing %s" % rg.resid() ) chain.remove_residue_group( rg ) SORTING_ON = { "input": lambda info: info.ensembler_input_index, "identity": lambda info: -info.ensembler_identity, "overlap": lambda info: -info.ensembler_overlap, "wrmsd": lambda info: info.ensembler_wrmsd, "unwrmsd": lambda info: info.ensembler_unwrmsd, } def write_merged_output(root, infos, chains, logger): logger.info( msg = output.heading( text = "Writing merged output file" ) ) logger.info( msg = output.blank() ) import iotbx.pdb import os.path pdb_root = iotbx.pdb.hierarchy.root() sequence_identity_remarks = [] identifier_remarks = [] logger.indent() for ( index, ( info, chain ) ) in enumerate( zip( infos, chains ), start = 1 ): model = iotbx.pdb.hierarchy.model() model.id = str( index ) pdb_root.append_model( model ) chain.id = "A" model.append_chain( chain ) logger.info( msg = "Model %s is %s" % ( index, info ) ) if info.identity is not None: remark = tbx_utils.get_phaser_model_identity_remark( mid = model.id, identity = info.identity ) else: remark = None sequence_identity_remarks.append( remark ) identifier_remarks.append( "REMARK MODEL %s: PDB %s, MODEL '%s', CHAIN %s" % ( index, os.path.basename( info.pdb.name ), info.model.id, info.chain.id, ) ) logger.dedent() output_file_name = "%s_merged.pdb" % root logger.info( msg = output.blank() ) logger.info( msg = "File name: %s" % output_file_name ) write_output_file( file_name = output_file_name, root = pdb_root, seqids = sequence_identity_remarks, identifiers = identifier_remarks, logger = logger ) logger.info( msg = output.section_close( text = "Output file written" ) ) return [ output_file_name ] def write_separate_output(root, infos, chains, logger): logger.info( msg = output.heading( text = "Writing separate output files" ) ) # Sort chains based on original file name and chainID output_data_for = {} for ( info, chain ) in zip( infos, chains ): output_data_for.setdefault( ( info.pdb.name, info.chain.id ), [] ).append( ( info, chain ) ) output_files = [] import iotbx.pdb import os.path for ( infile, cid ) in sorted( output_data_for ): output_file_name = "%s_%s%s.pdb" % ( root, os.path.basename( os.path.splitext( infile )[0] ), "_%s" % cid if cid else "", ) # Sort chains back according to modelID model_with = {} for ( info, chain ) in output_data_for[ ( infile, cid ) ]: data = model_with.setdefault( info.model.id, [] ) data.append( ( info, chain ) ) logger.info( "%s (model '%s') contains %s" % ( output_file_name, info.model.id, info, ) ) pdb_root = iotbx.pdb.hierarchy.root() sequence_identity_remarks = [] identifier_remarks = [] for ( index, mid ) in enumerate( sorted( model_with ) ): model = iotbx.pdb.hierarchy.model() model.id = mid pdb_root.append_model( model ) for ( info, chain ) in model_with[ mid ]: model.append_chain( chain ) if info.identity is not None: remark = tbx_utils.get_phaser_model_identity_remark( mid = mid, identity = info.identity ) else: remark = None sequence_identity_remarks.append( remark ) identifier_remarks.append( "REMARK MODEL %s: PDB %s, MODEL '%s', CHAIN %s" % ( index, os.path.basename( info.pdb.name ), mid, chain.id, ) ) write_output_file( file_name = output_file_name, root = pdb_root, seqids = sequence_identity_remarks, identifiers = identifier_remarks, logger = logger ) output_files.append( output_file_name ) logger.info( msg = output.section_close( text = "%d output files written" % len( output_files ) ) ) return output_files def write_chain_and_model_sorted_output_files(root, infos, chains, logger): logger.info( msg = output.heading( text = "Writing separate output file for each chain and model" ) ) logger.info( msg = output.blank() ) file_names = [] import iotbx.pdb import os.path logger.indent() for ( index, ( info, chain ) ) in enumerate( zip( infos, chains ), start = 1 ): pdb_root = iotbx.pdb.hierarchy.root() model = iotbx.pdb.hierarchy.model() pdb_root.append_model( model ) chain.id = "A" model.append_chain( chain ) logger.info( msg = "Model %s is %s" % ( index, info ) ) if info.identity is not None: seqremark = tbx_utils.get_phaser_model_identity_remark( mid = model.id, identity = info.identity ) else: seqremark = None identremark = "REMARK MODEL %s: PDB %s, MODEL '%s', CHAIN %s" % ( index, os.path.basename( info.pdb.name ), info.model.id, info.chain.id, ) output_file_name = "%s_%s%s.pdb" % ( root, os.path.basename( os.path.splitext( info.pdb.name )[0] ), "%s%s" % ( "_model%s" % info.model.id.strip() if info.model.id else "", "_chain%s" % info.chain.id if info.chain.id else "", ) ) write_output_file( file_name = output_file_name, root = pdb_root, seqids = [ seqremark ], identifiers = [ identremark ], logger = logger ) logger.info( msg = "File name: %s" % output_file_name ) file_names.append( output_file_name ) logger.dedent() logger.info( msg = output.section_close( text = "%d output files written" % len( file_names ) ) ) return file_names def write_output_file(file_name, root, seqids, identifiers, logger): with open( file_name, "w" ) as fout: if None not in seqids: fout.write( "%s\n" % "\n".join( seqids ) ) elif seqids.count( None ) != len( seqids ): logger.warning( msg = "Inconsistent sequence identity remarks, discarded" ) fout.write( "%s\n" % "\n".join( identifiers ) ) fout.write( "%s\n" % root.as_pdb_string() ) WRITE_OUTPUT_AS = { "merged": write_merged_output, "separate": write_separate_output, } PHIL_MASTER = libtbx.phil.parse( """ input .help = "Input files" { model = None .help = "Input model file" .type = path .optional = True .multiple = True alignment = None .help = "Input alignment file" .type = path .optional = True .multiple = True } output .help = "Output file(s)" { location = "" .type = path .short_caption = Output directory .style = hidden gui_output_dir = None .type = path .short_caption = Output directory .help = Sets base output directory for Phenix GUI - not used when \ run from the command line. .style = output_dir bold root = ensemble .help = "Output file root" .type = str .optional = False .short_caption = Output file root style = %(style)s .help = "Output options for ensemble" .type = choice .optional = False .short_caption = Output ensemble as sort = %(sort)s .help = "Sort ensemble components" .type = choice .optional = False .short_caption = Sort models by include scope libtbx.phil.interface.tracking_params } configuration .help = "General parameters for ensemble generation" .short_caption = Ensemble generation settings .style = menu_item { superposition .help = "Superposition setup" { method = %(superposition)s .help = "Superposition algorithm" .type = choice .optional = False .short_caption = Superposition algorithm convergence = 1.0E-4 .help = "Convergence criterion for superposition" .type = float .optional = False .short_caption = Superposition RMS covergence threshold } mapping = %(mapping)s .help = "Residue correspondence mapping" .type = choice .optional = False .short_caption = Residue mapping atoms = CA .help = "Include atom in superpositions" .type = str .optional = False .short_caption = Atoms to include in superposition clustering = 0.5 .help = "Cutoff distance for cluster analysis" .type = float .optional = False .short_caption = Cutoff distance for cluster analysis weighting .help = "Parameters for weighting scheme" .style = box { scheme = %(weighting)s .help = "Weighting scheme" .type = choice .optional = False .short_caption = Weighting scheme convergence = 1.0E-3 .help = "Convergence criterion for weight iteration" .type = float .optional = False .short_caption = Convergence criterion for weight iteration incremental_damping_factor = 1.5 .help = "Damping factor in recovery cycle" .type = float .optional = False max_damping_factor = 3.34 .help = "Quit recovery if cumulative damping factor is above" .type = float .optional = False robust_resistant .help = "Setting for robust-resistance scheme" { critical = 9 .help = "tolerance" .type = float .optional = False .short_caption = Robust resistant rms weighting term } } trim = False .help = "Trim differing loops" .type = bool .optional = False .short_caption = Trim residues deviating more than threshold trimming .help = "Trimming parameters" { threshold = 3.0 .help = "Max position RMSD for inclusion" .type = float .optional = False .short_caption = Trimming threshold } } """ % { "mapping": tbx_utils.choice_string( possible = MAPPING_PROCEDURE_FOR, default = "ssm" ), "superposition": tbx_utils.choice_string( possible = SUPERPOSITION_PROCEDURES_FOR, default = "gapless" ), "weighting": tbx_utils.choice_string( possible = WEIGHTING_SCHEME_SETUP_FOR, default = "robust_resistant" ), "sort": tbx_utils.choice_string( possible = SORTING_ON, default = "input" ), "style": tbx_utils.choice_string( possible = WRITE_OUTPUT_AS, default = "merged" ) }, process_includes=True) master_params = PHIL_MASTER # for documentation def process(pdbs, alignments, params, logger): # Assign functions (to check whether everything exists) matching_method = MAPPING_PROCEDURE_FOR[ params.configuration.mapping ] ( site_selection_method, setup_superposition_method ) = SUPERPOSITION_PROCEDURES_FOR[ params.configuration.superposition.method ] setup_weighting_method = WEIGHTING_SCHEME_SETUP_FOR[ params.configuration.weighting.scheme ] sorting_method = SORTING_ON[ params.output.sort ] # Collates chains in pdb files and discard unknowns logger.info( msg = output.subtitle( text = "Chain selection" ) ) infos = find_chains( pdbs = pdbs, logger = logger ) for ( index, info ) in enumerate( infos, start = 1 ): info.ensembler_input_index = index if len( infos ) < MIN_MODEL_COUNT: raise RuntimeError, "Less than %s chains for superposition" % MIN_MODEL_COUNT # Match residue groups - atom selection applied at this step logger.info( msg = output.subtitle( text = "Residue mapping" ) ) mapping = matching_method( infos = infos, alignments = alignments, logger = logger ) logger.info( msg = "Residue mapping: %s aligned positions" % len( mapping ) ) logger.info( msg = output.underlined( text = "Sequence identity statistics:" ) ) identities = calculate_and_print_sequence_identity_statistics( infos = infos, mapping = mapping, logger = logger ) assert len( identities ) == len( infos ) for ( info, value ) in zip( infos, identities ): info.ensembler_identity = value # Convert to atom mapping logger.info( msg = output.subtitle( text = "Atom mapping" ) ) selector = get_atom_selector( selection = params.configuration.atoms, logger = logger ) mapping = to_atom_mapping( mapping = mapping, selector = selector, logger = logger ) logger.info( msg = "Atom mapping: %s aligned positions" % len( mapping ) ) logger.info( msg = output.underlined( text = "Overlap statistics:" ) ) overlaps = calculate_and_print_overlap_statistics( infos = infos, mapping = mapping, logger = logger ) assert len( overlaps ) == len( infos ) for ( info, value ) in zip( infos, overlaps ): info.ensembler_overlap = value # Convert to sites and setup superposition logger.info( msg = output.subtitle( text = "Superposition" ) ) sites = site_selection_method( mapping = mapping, logger = logger ) if len( sites ) < MIN_SITE_COUNT: raise RuntimeError, "Less than %s sites for superposition" % MIN_SITE_COUNT logger.info( msg = "Site selection: %s sites for superposition" % len( sites ) ) from scitbx.array_family import flex weights = flex.double( [ 1 ] * len( sites ) ) superposition = setup_superposition_method( sites = sites, weights = weights, logger = logger ) # Superpose structures logger.info( msg = "Weighting scheme: %s" % params.configuration.weighting.scheme ) weighting = setup_weighting_method( params = params.configuration.weighting ) logger.info( msg = "Iterating weighting with superposition..." ) iterative_weighted_superposition( superposition = superposition, superposition_convergence = params.configuration.superposition.convergence, weighting = weighting, weights = weights, weight_convergence = params.configuration.weighting.convergence, incremental_damping_factor = params.configuration.weighting.incremental_damping_factor, max_damping_factor = params.configuration.weighting.max_damping_factor, logger = logger ) logger.info( msg = "Superposition complete" ) logger.info( msg = "Rmsd between all chains: %.3f" % superposition.rmsd() ) # Calculate and save results and statistics logger.info( msg = output.subtitle( text = "Superposition results" ) ) logger.info( msg = output.underlined( text = "Transformations" ) ) ( rotations, translations ) = superposition.transformations() assert len( infos ) == len( translations ) assert len( infos ) == len( rotations ) for ( info, rot, tra ) in zip( infos, rotations, translations ): logger.info( msg = "%s:" % info ) logger.info( msg = " Rotation:\n " + ( "%7.3f" * 9 ) % rot.elems ) logger.info( msg = " Translation: " + ( "%9.4f" * 3 ) % tra.elems ) info.ensembler_rotation = rot info.ensembler_translation = tra logger.info( msg = output.subtitle( text = "Rmsds" ) ) ( wrmsds, unwrmsds ) = calculate_rmsds( superposition = superposition, infos = infos, logger = logger ) assert len( wrmsds ) == len( infos ) assert len( unwrmsds ) == len( infos ) for ( info, wrmsd, unwrmsd ) in zip( infos, wrmsds, unwrmsds ): info.ensembler_wrmsd = wrmsd info.ensembler_unwrmsd = unwrmsd # Cluster analysis logger.info( msg = output.subtitle( text = "Cluster analysis" ) ) clusters = cluster_analysis( superposition = superposition, infos = infos, critical = params.configuration.clustering, logger = logger ) assert len( clusters ) == len( infos ) for ( info, value ) in zip( infos, clusters ): info.ensembler_cluster = value # Calculate trimming data (on demand) if params.configuration.trim: logger.info( msg = output.subtitle( text = "Chain trimming" ) ) logger.info( msg = "Calculating data..." ) trimming_masks = calculate_trimming_masks( infos = infos, sites = sites, deviations = superposition.distance_squares_from_average_structure(), threshold = params.configuration.trimming.threshold, logger = logger ) assert len( infos ) == len( trimming_masks ) for ( info, mask ) in zip( infos, trimming_masks ): info.ensembler_trimming_mask = mask logger.info( msg = output.subtitle( text = "Summary" ) ) logger.info( msg = "Sort order: %s" % params.output.sort ) infos.sort( key = sorting_method ) print_chain_infos( infos = infos, logger = logger ) logger.info( msg = tbx_utils.format_table( captions = [ "Index", "Input index", "Identity", "Overlap", "Weighted rmsds", "Unweighted rmsds", ], formats = [ "%d", "%d", "%.3f", "%.3f", "%.3f", "%.3f" ], lines = [ [ i, info.ensembler_input_index, info.ensembler_identity, info.ensembler_overlap, info.ensembler_wrmsd, info.ensembler_unwrmsd, ] for ( i, info ) in enumerate( infos, start = 1 ) ] ) ) return infos def run(args, out, verbosity = output.SingleStream.INFO): logger = output.SingleStream( stream = out, level = verbosity ) logger.info( msg = output.banner( text = "%s version %s" % ( PROGRAM, VERSION ) ) ) accumulator = tbx_utils.PhilAccumulator( master_phil = PHIL_MASTER ) from iotbx import bioinformatics accumulator.register_file_handler( handler = tbx_utils.ExtensionFileHandler( extensions = bioinformatics.known_alignment_formats(), template = "input.alignment=%s" ) ) accumulator.register_file_handler( handler = tbx_utils.ExtensionFileHandler( extensions = [ ".pdb", ".ent" ], template = "input.model=%s" ) ) for argument in args: argument.process( accumulator = accumulator ) from libtbx.utils import Sorry try: merged_phil = accumulator.merge() except RuntimeError, e: raise Sorry, e logger.info( msg = output.underlined( text = "All configuration options:" ) ) logger.info( msg = merged_phil.as_str() ) try: params = merged_phil.extract() except RuntimeError, e: raise Sorry, e # Sanity checks if not params.input.model: raise Sorry, "No input pdb files" # Main functionality logger.info( msg = output.underlined( text = "Reading pdb files:" ) ) logger.indent() logger.info( msg = "\n".join( "%d. %s" % ( i, f ) for ( i, f ) in enumerate( params.input.model, start = 1 ) ) ) logger.dedent() pdbs = tbx_utils.read_in_files( file_names = params.input.model, object_type = tbx_utils.PDBObject, ) logger.info( msg = output.section_close( text = "%d files read" % len( pdbs ) ) ) if params.input.alignment: logger.info( msg = output.underlined( text = "Reading alignment files:" ) ) logger.indent() logger.info( msg = "\n".join( "%d. %s" % ( i, f ) for ( i, f ) in enumerate( params.input.alignment ) ) ) logger.dedent() alignments = tbx_utils.read_in_files( file_names = params.input.alignment, object_type = tbx_utils.AlignmentObject ) logger.info( msg = output.section_close( text = "%d files read" % len( pdbs ) ) ) else: alignments = [] try: infos = process( pdbs = pdbs, alignments = alignments, params = params, logger = logger ) except RuntimeError, e: raise Sorry, e logger.info( msg = output.subtitle( text = "Prepare output" ) ) # Create copies chains = [ info.chain.detached_copy() for info in infos ] for ( info, chain ) in zip( infos, chains ): sites = ( info.ensembler_rotation.elems * chain.atoms().extract_xyz() + info.ensembler_translation ) chain.atoms().set_xyz( sites ) if params.configuration.trim: assert all( hasattr( info, "ensembler_trimming_mask" ) for info in infos ) logger.info( msg = "Chain trimming" ) logger.info( msg = output.blank() ) trim_chains( infos = infos, chains = chains, logger = logger ) # Write output import os.path try: file_names = WRITE_OUTPUT_AS[ params.output.style ]( root = os.path.join( params.output.location, params.output.root ), infos = infos, chains = chains, logger = logger ) except IOError, e: raise Sorry, e logger.info( msg = "%s finished" % PROGRAM ) return file_names # XXX: adaptor for Phenix GUI, which handles exceptions separately. class launcher (runtime_utils.target_with_save_result) : def run (self) : os.makedirs(self.output_dir) os.chdir(self.output_dir) import sys from libtbx.utils import Sorry factory = tbx_utils.PhilArgumentFactory( master_phil = PHIL_MASTER ) try : output_files = run( args= [ factory( argument = arg ) for arg in self.args ], out=sys.stdout) #output_files = process(merged_phil, command_name="phenix.ensembler") except Exception, e : raise # raise Sorry("Error encountered in processing: %s" % str(e)) else : if (output_files is None) or (len(output_files) == 0) : raise Sorry("Empty result set; no output file(s) written.") else : if isinstance(output_files, list) : return [ os.path.abspath(fn) for fn in output_files ] return output_files def finish_job (result) : files = [] if isinstance(result, str) : if os.path.isfile(result) : pdb_out = os.path.abspath(result) files.append((pdb_out, "Multi-model PDB file")) elif isinstance(result, list) : for file_name in result : pdb_out = os.path.abspath(file_name) files.append((pdb_out, "Individual model")) return (files, [])