from phaser import chisel from phaser import mmt from phaser import tbx_utils from phaser import output from iotbx import bioinformatics import libtbx.phil from libtbx import runtime_utils from libtbx.utils import Sorry import operator import sys # Global constants VERSION = "0.3.3" PROGRAM = "sculptor" KNOWN = set( [ mmt.PROTEIN ] ) # Reusable PHIL extracts SIMILARITY_PHIL = """ calculation .help = "Normalized values (~1: perfect match; ~0: random sequence: ~-1: gap)" .short_caption = Sequence similarity calculation .help = "Configure sequence similarity calculation" .style = box auto_align { matrix = %(matrix)s .help = "Similarity matrix" .short_caption = Similarity matrix .type = choice .optional = False window = %(window)s .help = "Averaging window width" .type = int .optional = False weighting = %(weighting)s .help = "Weighting scheme" .type = choice .optional = False }""" TRANSFORMATION_PHIL = """ factor = %s .help = "Transform values by multiplying with a factor" .type = float .optional = False .short_caption = Scale current B-factors by .style = bold""" SS_MATRICES = mmt.PROTEIN.similarity_matrix_suite.names() # # Factory functions # # Infrastructure def get_unit_weights(window): from scitbx.array_family import flex return flex.double( [ 1 ] * ( 2 * window + 1 ) ) def get_triangular_weights(window): from scitbx.array_family import flex return flex.double( range( 1, window + 1 ) + [ window + 1 ] + range( window, 0, -1 ) ) WEIGHTING_SCHEME_NAMED = { "uniform": get_unit_weights, "triangular": get_triangular_weights, } # Mainchain step def gap_mainchain_algorithm(params): return chisel.MDAGap() gap_mainchain_algorithm.PHIL = \ """ .short_caption = Filter by residue presence .help = "Delete residue if aligned with gap" .style = box auto_align { }""" def similarity_mainchain_algorithm(params): w = params.calculation.window weights = WEIGHTING_SCHEME_NAMED[ params.calculation.weighting ]( window = w ) return chisel.MDASimilarity( ss_calc = chisel.LinearAveragedSequenceSimilarity( matrix = params.calculation.matrix, unknown = 0, window = w, weights = weights ), threshold = params.threshold ) similarity_mainchain_algorithm.PHIL = \ """ .short_caption = Filter by sequence similarity, truncate by specified threshold .help = "Delete residue if sequence similarity is low" .style = box auto_align { threshold = -0.2 .help = "Threshold to accept a residue" .short_caption = Threshold to accept a residue .type = float .optional = False %s }""" % ( SIMILARITY_PHIL % { "matrix": tbx_utils.choice_string( possible = SS_MATRICES, default = SS_MATRICES[0] ), "window": 5, "weighting": tbx_utils.choice_string( possible = WEIGHTING_SCHEME_NAMED, default = "triangular" ), } ) def threshold_adjust_similarity_mainchain_algorithm(params): w = params.calculation.window weights = WEIGHTING_SCHEME_NAMED[ params.calculation.weighting ]( window = w ) return chisel.MDAThresholdAdjustSimilarity( ss_calc = chisel.LinearAveragedSequenceSimilarity( matrix = params.calculation.matrix, unknown = 0, window = w, weights = weights ), offset = params.offset ) threshold_adjust_similarity_mainchain_algorithm.PHIL = \ """ .short_caption = Filter by sequence similarity, truncate to target completeness .help = "Delete residues based on sequence similarity to get same number of gaps \ as the Schwarzenbacher algorithm" .style = box auto_align { offset = 0.0 .help = "Completeness in fraction of model length \ (0.0 = completeness from Schwarzenbacher algorithm, useful range: +/-0.05)" .short_caption = Model completeness wrt Schwarzenbacher-algorithm .type = float .optional = False %s }""" % ( SIMILARITY_PHIL % { "matrix": tbx_utils.choice_string( possible = SS_MATRICES, default = SS_MATRICES[0] ), "window": 5, "weighting": tbx_utils.choice_string( possible = WEIGHTING_SCHEME_NAMED, default = "triangular" ), } ) def remove_long_mainchain_algorithm(params): return chisel.MDARemoveLong( min_length = params.min_length) remove_long_mainchain_algorithm.PHIL = \ """ .short_caption = Filter by residue presence .help = "Delete residue if aligned with gap" .style = box auto_align { min_length = 3 .help = "Minimum length for mainchain segment to remove" .type = int .optional = False }""" MAINCHAIN_ALGORITHM_NAMED = { "gap": gap_mainchain_algorithm, "threshold_based_similarity": similarity_mainchain_algorithm, "completeness_based_similarity": threshold_adjust_similarity_mainchain_algorithm, "remove_long": remove_long_mainchain_algorithm, } def remove_short_polishing_algorithm(params): return chisel.MPARemoveShort( min_length = params.minimum_length ) remove_short_polishing_algorithm.PHIL = \ """ .help = "Delete short unconnected segments" { minimum_length = 3 .help = "Minimum length" .type = int .optional = False .short_caption = Minimum chain length to keep }""" def keep_regular_polishing_algorithm(params): return chisel.MPAKeepRegular( max_length = params.maximum_length ) keep_regular_polishing_algorithm.PHIL = \ """ .help = "Keep residues in secondary structure" { maximum_length = 1 .help = "Maximum length" .type = int .optional = False .short_caption = Maximum chain length to keep }""" POLISHING_ALGORITHM_NAMED = { "remove_short": remove_short_polishing_algorithm, "keep_regular": keep_regular_polishing_algorithm, } # Sidechain step def schwarzenbacher_sidechain_pruning_algorithm(params): return chisel.SPASchwarzenbacher( level = params.pruning_level ) schwarzenbacher_sidechain_pruning_algorithm.PHIL = \ """ .help = "Truncate atoms if target residue != source residue" { pruning_level = 2 .help = "Level of truncation" .type = int .optional = False }""" def similarity_sidechain_pruning_algorithm(params): w = params.calculation.window weights = WEIGHTING_SCHEME_NAMED[ params.calculation.weighting ]( window = w ) return chisel.SPASimilarity( ss_calc = chisel.LinearAveragedSequenceSimilarity( matrix = params.calculation.matrix, unknown = 0, window = w, weights = weights ), lower = params.full_truncation_limit, upper = params.full_length_limit, level = params.pruning_level ) similarity_sidechain_pruning_algorithm.PHIL = \ """ .help = "Truncate atoms based on sequence similarity" .short_caption = Similarity-based truncation .style = box auto_align { pruning_level = 2 .help = "Level of intermediate truncation" .type = int .optional = False full_length_limit = 0.2 .help = "Limit of no truncation" .type = float .optional = False full_truncation_limit = -0.2 .help = "Limit for full truncation" .type = float .optional = False %s }""" % ( SIMILARITY_PHIL % { "matrix": tbx_utils.choice_string( possible = SS_MATRICES, default = SS_MATRICES[0] ), "window": 5, "weighting": tbx_utils.choice_string( possible = WEIGHTING_SCHEME_NAMED, default = "triangular" ), } ) SIDECHAIN_PRUNING_ALGORITHM_NAMED = { "schwarzenbacher": schwarzenbacher_sidechain_pruning_algorithm, "similarity": similarity_sidechain_pruning_algorithm, } # Bfactor step def original_bfactor_algorithm(params): return chisel.BPAOriginalBfactor( factor = params.factor ) original_bfactor_algorithm.PHIL = \ """ .help = "Use original bfactors to predict new B-values" .short_caption = Original B-factors .style = auto_align box { %s }""" % ( TRANSFORMATION_PHIL % 1 ) def similarity_bfactor_algorithm(params): w = params.calculation.window weights = WEIGHTING_SCHEME_NAMED[ params.calculation.weighting ]( window = w ) return chisel.BPASequenceSimilarity( ss_calc = chisel.LinearAveragedSequenceSimilarity( matrix = params.calculation.matrix, unknown = 0, window = w, weights = weights ), factor = params.factor ) similarity_bfactor_algorithm.PHIL = \ """ .help = "Use sequence similarity to predict new B-values" .short_caption = Sequence similarity-based .style = box auto_align { %s %s }""" % ( TRANSFORMATION_PHIL % -100, SIMILARITY_PHIL % { "matrix": tbx_utils.choice_string( possible = SS_MATRICES, default = SS_MATRICES[0] ), "window": 5, "weighting": tbx_utils.choice_string( possible = WEIGHTING_SCHEME_NAMED, default = "triangular" ), }, ) def asa_bfactor_algorithm(params): if params.precision < 1: raise Sorry, "Invalid precision: >0 required" if params.probe_radius < 0: raise Sorry, "Invalid probe_radius: >=0 required" return chisel.BPAAccessibleSurfaceArea( probe = params.probe_radius, precision = params.precision, factor = params.factor ) asa_bfactor_algorithm.PHIL = \ """ .help = "Use accessible surface area to predict new B-values" .short_caption = Surface area-based .style = box auto_align { %s precision = 960 .help = "Number of points per atom" .type = int .optional = False probe_radius = 1.4 .help = "Radius for probing surface accessibility" .type = float .optional = False }""" % ( TRANSFORMATION_PHIL % 2 ) BFACTOR_ALGORITHM_NAMED = { "original": original_bfactor_algorithm, "similarity": similarity_bfactor_algorithm, "asa": asa_bfactor_algorithm, } # Renumbering step def target_renumber_step(params): return chisel.MRASequence.Target( start = params.start ) def model_renumber_step(params): return chisel.MRASequence.Model( start = params.start ) def no_renumber_step(params): return chisel.MRAOriginal() RENUMBER_STEP_BASED_ON = { "target": target_renumber_step, "model": model_renumber_step, "original": no_renumber_step, } # Rename step def target_sequence_rename_step(keep_ptm): return chisel.SRSingleLetterTarget() if keep_ptm else chisel.SRSequence() def model_sequence_rename_step(): return chisel.SROriginal() # Output step def regularize_hierarchy(root): for a in root.atoms(): if a.occ < 0.005: a.parent().remove_atom( a ) for model in root.models(): if not model.atoms(): root.remove_model( model ) continue for chain in model.chains(): if not chain.atoms(): model.remove_chain( chain ) continue for rg in chain.residue_groups(): if not rg.atoms(): chain.remove_residue_group( rg ) continue for ag in rg.atom_groups(): if not ag.atoms(): rg.remove_atom_group( ag ) root.atoms_reset_serial() def write_pdb_output(root, blueprint, identities, file_root, suffix, logger): logger.info( msg = "PDB output requested" ) logger.debug( msg = "Associate modifications with hierarchy objects..." ) copy_root = root.deep_copy() indexer = chisel.ReverseIndexer( root = copy_root ) instructions = blueprint.instructions( r_indexer = indexer ) logger.info( msg = "Applying modifications...") for ( i, e ) in instructions: m = i.log( entity = e ) if e.is_connected(): logger.debug( m ) i.apply_to( entity = e ) else: logger.debug( msg = "%s - execution omitted because entity not attached" % m ) logger.info( msg = "Regularizing hierarchy..." ) regularize_hierarchy( root = copy_root ) logger.info( msg = "Generating sequence identity records..." ) remarks = [] for ( d, i ) in identities: if not i: continue h = list( d.as_entity( indexer = indexer ).hierarchy() ) assert h if h[0].level() == 0: remarks.append( tbx_utils.get_phaser_model_identity_remark( mid = h[1].object().id, identity = i * 100.0 ) ) logger.info( msg = "Preparing output file..." ) file_name = "%s_%s.pdb" % ( file_root, suffix ) fout = file( file_name, "w" ) for r in remarks: fout.write( "%s\n" % r ) fout.write( "%s\n" % copy_root.as_pdb_string() ) fout.close() return file_name def write_coot_output(root, blueprint, identities, file_root, suffix, logger): import iotbx.pdb logger.info( msg = "Coot output requested" ) logger.debug( msg = "Associate modifications with hierarchy objects..." ) copy_root = root.deep_copy() indexer = chisel.ReverseIndexer( root = copy_root ) instructions = blueprint.instructions( r_indexer = indexer ) logger.debug( msg = "Applying modifications (but not deletions)..." ) moribund_residues = [] moribund_atoms = [] attached_atoms = [] level_rg = chisel.Entity.LEVEL_FOR[ iotbx.pdb.hierarchy.residue_group ] level_ag = chisel.Entity.LEVEL_FOR[ iotbx.pdb.hierarchy.atom_group ] level_a = chisel.Entity.LEVEL_FOR[ iotbx.pdb.hierarchy.atom ] for ( i, e ) in instructions: if isinstance( i, chisel.DeleteInstruction ): logger.debug( msg = "%s deferred" % i.log( entity = e ) ) if e.level() == level_rg: moribund_residues.append( e ) elif e.level() == level_a: moribund_atoms.append( e ) elif isinstance( i, chisel.AttachInstruction ): logger.debug( msg = "%s applied and logged" % i.log( entity = e ) ) if e.level() == level_ag: attached_atoms.append( ( e, i.data.name ) ) i.apply_to( entity = e ) else: logger.debug( msg = "%s applied" % i.log( entity = e ) ) i.apply_to( entity = e ) logger.info( msg = "Regularizing hierarchy..." ) regularize_hierarchy( root = copy_root ) logger.info( msg = "Preparing output file..." ) file_name = "%s_%s.coot" % ( file_root, suffix ) fout = file( file_name, "w" ) fout.write( "(\n" ) fout.write( "%s,\n" % str( chisel.to_coot_s_expression( root = copy_root ) ) ) fout.write( "(\n" ) for e in moribund_residues: fout.write( " %s,\n" % str( e.coot_selector() ) ) fout.write( " ),\n" ) fout.write( "(\n" ) for e in moribund_atoms: fout.write( " %s,\n" % str( e.coot_selector() ) ) fout.write( " ),\n" ) fout.write( "(\n" ) for ( e, n ) in attached_atoms: fout.write( " %s,\n" % str( e.coot_selector() + ( n, e.object().altloc.strip() ) ) ) fout.write( " ),\n" ) fout.write( ")\n" ) fout.close() return file_name OUTPUT_METHOD_FOR = { "pdb": write_pdb_output, "coot": write_coot_output, } # # Factory functions sculptor # def get_algorithms(collection, params): return [ collection[ n ]( params = getattr( params, n ) if hasattr( params, n ) else None ) for n in params.use ] def get_sculptor_object(deletion, polishing, renumber, pruning, completion, rename, bfactor, keep_ptm = False): return chisel.Sculptor( deletions = get_algorithms( collection = MAINCHAIN_ALGORITHM_NAMED, params = deletion ), polishings = get_algorithms( collection = POLISHING_ALGORITHM_NAMED, params = polishing ), renumber = RENUMBER_STEP_BASED_ON[ renumber.use ]( params = renumber ), prunings = get_algorithms( collection = SIDECHAIN_PRUNING_ALGORITHM_NAMED, params = pruning ), completions = set( completion ), rename = target_sequence_rename_step( keep_ptm = keep_ptm ) if rename else model_sequence_rename_step(), bfactors = get_algorithms( collection = BFACTOR_ALGORITHM_NAMED, params = bfactor ), minimum_b = bfactor.minimum ) def get_discarder_object(params): if not params: known = [] else: known = reduce( operator.add, [ [ i.strip() for i in ind.split(",") ] for ind in params ] ) return chisel.Discarder( keep = known ) def process_multiprotocol_phil(params, logger): logger.info( msg = "Predefined protocols will be used" ) protocol_phil = libtbx.phil.parse( PHIL_PROTOCOL_MP ) protocol_keys = set() for key in params.macromolecule.protocols: if key in PROTOCOL_COMBINATIONS: protocol_keys.update( PROTOCOL_COMBINATIONS[ key ] ) else: protocol_keys.add( key ) protocols = [] for key in protocol_keys: logger.info( msg = output.heading( text = "Protocol %s:" % key ) ) mpp = protocol_phil.fetch( sources = [ get_protocol_phil( key = key ) ] ) logger.info( msg = mpp.as_str() ) prot = mpp.extract() protocols.append( ( "protocol %s" % key, "_%s" % key, get_sculptor_object( deletion = prot.macromolecule.deletion, polishing = prot.macromolecule.polishing, pruning = prot.macromolecule.pruning, completion = prot.macromolecule.completion, bfactor = prot.macromolecule.bfactor, renumber = params.macromolecule.renumber, rename = params.macromolecule.rename, keep_ptm = params.macromolecule.keep_ptm_if_base_residues_agree, ), ) ) return protocols def process_user_defined_protocol_phil(params, logger): logger.info( msg = "User-specified protocol will be used" ) logger.info( msg = output.blank() ) protocol_phil = libtbx.phil.parse( PHIL_PROTOCOL_MP ) mpp = protocol_phil.fetch( sources = [ protocol_phil.format( python_object = params ) ] ) logger.info( msg = mpp.as_str() ) return [ ( "user defined protocol", "", get_sculptor_object( deletion = params.macromolecule.deletion, polishing = params.macromolecule.polishing, renumber = params.macromolecule.renumber, pruning = params.macromolecule.pruning, completion = params.macromolecule.completion, rename = params.macromolecule.rename, bfactor = params.macromolecule.bfactor, keep_ptm = params.macromolecule.keep_ptm_if_base_residues_agree, ) ) ] # Predefined protocols # This has to be first as later PHILs use it to know the keys PROTOCOLS = { "1": """ macromolecule { deletion { use = gap } polishing { use = None } pruning { use = schwarzenbacher schwarzenbacher { pruning_level = 2 } } bfactor { use = original minimum = 10 } completion = cbeta } """, #1 "2": """ macromolecule { deletion { use = completeness_based_similarity completeness_based_similarity { offset = 0.0 calculation { matrix = blosum62 window = 5 weighting = triangular } } } polishing { use = None } pruning { use = schwarzenbacher schwarzenbacher { pruning_level = 2 } } bfactor { use = original minimum = 10 } completion = cbeta } """, #2 "3": """ macromolecule { deletion { use = gap } polishing { use = None } pruning { use = similarity similarity { pruning_level = 2 full_length_limit = 0.2 full_truncation_limit = -0.2 calculation { matrix = blosum62 window = 1 weighting = triangular } } } bfactor { use = original minimum = 10 } completion = cbeta } """, #3 "4": """ macromolecule { deletion { use = completeness_based_similarity completeness_based_similarity { offset = 0.0 calculation { matrix = blosum62 window = 5 weighting = triangular } } } polishing { use = None } pruning { use = similarity similarity { pruning_level = 2 full_length_limit = 0.2 full_truncation_limit = -0.2 calculation { matrix = blosum62 window = 1 weighting = triangular } } } bfactor { use = original minimum = 10 } completion = cbeta } """, #4 "5": """ macromolecule { deletion { use = gap } polishing { use = None } pruning { use = schwarzenbacher schwarzenbacher { pruning_level = 2 } } bfactor { use = similarity minimum = 10 similarity { factor = -80 calculation { matrix = blosum62 window = 5 weighting = triangular } } } completion = cbeta } """, #5 "6": """ macromolecule { deletion { use = completeness_based_similarity completeness_based_similarity { offset = 0.0 calculation { matrix = blosum62 window = 5 weighting = triangular } } } polishing { use = None } pruning { use = schwarzenbacher schwarzenbacher { pruning_level = 2 } } bfactor { use = similarity minimum = 10 similarity { factor = -80 calculation { matrix = blosum62 window = 5 weighting = triangular } } } completion = cbeta } """, #6 "7": """ macromolecule { deletion { use = gap } polishing { use = None } pruning { use = schwarzenbacher schwarzenbacher { pruning_level = 2 } } bfactor { use = asa minimum = 10 asa { factor = 12 precision = 960 probe_radius = 1.4 } } completion = cbeta } """, #7 "8": """ macromolecule { deletion { use = completeness_based_similarity completeness_based_similarity { offset = 0.0 calculation { matrix = blosum62 window = 5 weighting = triangular } } } polishing { use = None } pruning { use = schwarzenbacher schwarzenbacher { pruning_level = 2 } } bfactor { use = asa minimum = 10 asa { factor = 12 precision = 960 probe_radius = 1.4 } } completion = cbeta } """, #8 "9": """ macromolecule { deletion { use = gap } polishing { use = None } pruning { use = schwarzenbacher schwarzenbacher { pruning_level = 2 } } bfactor { use = similarity+asa minimum = 10 asa { factor = 8 precision = 960 probe_radius = 1.4 } similarity { factor = -60 calculation { matrix = blosum62 window = 5 weighting = triangular } } } completion = cbeta } """, #9 "10": """ macromolecule { deletion { use = completeness_based_similarity completeness_based_similarity { offset = 0.0 calculation { matrix = blosum62 window = 5 weighting = triangular } } } polishing { use = None } pruning { use = schwarzenbacher schwarzenbacher { pruning_level = 2 } } bfactor { use = similarity+asa minimum = 10 asa { factor = 8 precision = 960 probe_radius = 1.4 } similarity { factor = -60 calculation { matrix = blosum62 window = 5 weighting = triangular } } } completion = cbeta } """, #10 "11": """ macromolecule { deletion { use = gap } polishing { use = None } pruning { use = similarity similarity { pruning_level = 2 full_length_limit = 0.2 full_truncation_limit = -0.2 calculation { matrix = blosum62 window = 1 weighting = triangular } } } bfactor { use = similarity+asa minimum = 10 asa { factor = 8 precision = 960 probe_radius = 1.4 } similarity { factor = -60 calculation { matrix = blosum62 window = 5 weighting = triangular } } } completion = cbeta } """, #11 "12": """ macromolecule { deletion { use = completeness_based_similarity completeness_based_similarity { offset = 0.0 calculation { matrix = blosum62 window = 5 weighting = triangular } } } polishing { use = None } pruning { use = similarity similarity { pruning_level = 2 full_length_limit = 0.2 full_truncation_limit = -0.2 calculation { matrix = blosum62 window = 1 weighting = triangular } } } bfactor { use = similarity+asa minimum = 10 asa { factor = 8 precision = 960 probe_radius = 1.4 } similarity { factor = -60 calculation { matrix = blosum62 window = 5 weighting = triangular } } } completion = cbeta } """, #12 "13": """ macromolecule { deletion { use = threshold_based_similarity threshold_based_similarity { threshold = -0.2 calculation { matrix = blosum62 window = 3 weighting = triangular } } } polishing { use = None } pruning { use = schwarzenbacher schwarzenbacher { pruning_level = 2 } } bfactor { use = original minimum = 10 } completion = cbeta } """, #13 } PROTOCOL_COMBINATIONS = { "all": PROTOCOLS.keys(), "minimal": [ "1", "8", "12" ], } PROTOCOL_PHILS = {} def get_protocol_phil(key): if key not in PROTOCOL_PHILS: if key not in PROTOCOLS: raise Sorry, "Unknown protocol '%s': not in %s" % ( key, ",".join( sorted( PROTOCOLS ) ), ) PROTOCOL_PHILS[ key ] = libtbx.phil.parse( PROTOCOLS[ key ] ) return PROTOCOL_PHILS[ key ] def get_known_protocols(): return PROTOCOLS.keys() + PROTOCOL_COMBINATIONS.keys() # Phil setup def configuration_phil(phil_for): return "\n".join( [ "%s\n%s" % ( k, v.PHIL ) for ( k, v ) in phil_for.items() ] ) PHIL_INPUT = """ input .help = "Input files" { model .help = "Input pdb file" .optional = True .short_caption = PDB file { file_name = None .help = "PDB file name" .short_caption = PDB file .optional = False .type = path .style = bold file_type:pdb selection = all .help = "Selection string" .optional = False .short_caption = Atom selection .input_size = 400 .type = str .style = bold remove_alternate_conformations = False .help = "Remove alternate conformations" .type = bool .optional = False .style = bold noauto sanitize_occupancies = False .help = "Sets occupancies > 1.0 to 1.0" .type = bool .optional = False .style = noauto } alignment .help = "Input alignment file" .multiple = True .optional = True .short_caption = Sequence alignment file .style = auto_align { file_name = None .optional = False .type = path target_index = 1 .help = "Index of target sequence in alignment" .type = int .optional = False } sequence .help = "Input sequence file" .multiple = True .optional = True .style = auto_align { file_name = None .help = "Sequence file" .short_caption = Sequence file .type = path chain_ids = None .optional = False .short_caption = Chain IDs .type = strings } } """ PHIL_OUTPUT = """ output .help = "Output options" { include scope libtbx.phil.interface.tracking_params folder = . .help = "Output file folder" .type = path .optional = False .short_caption = "Output file folder" .input_size = 400 .style = directory root = sculpt .help = "Output file root" .type = str .optional = False .short_caption = Output file base .input_size = 400 .style = bold format = %(output_format)s .help = "Output file format" .type = choice .optional = False .style = hidden } """ % { "output_format": tbx_utils.choice_string( possible = OUTPUT_METHOD_FOR, default = "pdb" ), } PHIL_COMMON_OPERATIONS = """ renumber .short_caption = Numbering { use = %(renumber_use)s .help = "Mainchain numbering; (* = selected; None: disable)" .type = choice .optional = False .short_caption = Mainchain numbering method start = 1 .help = "Number for first residue" .type = int .optional = False .short_caption = First residue number for sequence-based numbering } rename = True .help = "True: enable; False: disable" .type = bool .optional = False .short_caption = Rename residues based on target sequence .style = bold keep_ptm_if_base_residues_agree = False .help = "Keep post-translational modification if residues agree" .type = bool .optional = False """ % { "renumber_use": tbx_utils.choice_string( possible = RENUMBER_STEP_BASED_ON, default = "target" ), } PHIL_PROTOCOL_OPERATIONS = """ deletion .help = "Configure mainchain deletion" .short_caption = Main-chain deletion .style = menu_item auto_align box { use = %(mainchain_use)s .help = "Available algorithms (* = active)" .short_caption = Mainchain deletion method .type = choice ( multi = True ) .optional = True .style = bold noauto cols:2 %(mainchain_configure)s } polishing .help = "Configure mainchain polishing" .short_caption = Main-chain polishing .style = menu_item auto_align box { use = %(polishing_use)s .help = "Available algorithms (* = active)" .type = choice ( multi = True ) .short_caption = Mainchain polishing method .optional = True .style = bold noauto %(polishing_configure)s } pruning .help = "Configure sidechain pruning" .short_caption = Sidechain pruning .style = menu_item auto_align box { use = %(pruning_use)s .help = "Available algorithms (* = active)" .short_caption = Sidechain pruning method .type = choice ( multi = True ) .optional = True .style = bold noauto %(pruning_configure)s } bfactor .help = "Configure bfactor prediction" .short_caption = B-factor modifications .style = menu_item auto_align box { use = %(bfactor_use)s .help = "Available algorithms (* = active)" .type = choice ( multi = True ) .optional = True .short_caption = B-factor predication algorithm .style = bold noauto minimum = 10 .help = "Minimum allowed value (a constant is added if any B-factors would fall below this value)" .type = float .optional = False .short_caption = Add constant to current B-factors .style = bold %(bfactor_configure)s } completion = %(completion_use)s .help = "Sidechain completion algorithms (* = active)" .type = choice ( multi = True ) .optional = True .short_caption = Sidechain completion method .style = menu_item """ % { "mainchain_use": tbx_utils.choice_string( possible = MAINCHAIN_ALGORITHM_NAMED, default = "gap" ), "mainchain_configure": configuration_phil( phil_for = MAINCHAIN_ALGORITHM_NAMED ), "polishing_use": tbx_utils.choice_string( possible = POLISHING_ALGORITHM_NAMED, default = None ), "polishing_configure": configuration_phil( phil_for = POLISHING_ALGORITHM_NAMED ), "pruning_use": tbx_utils.choice_string( possible = SIDECHAIN_PRUNING_ALGORITHM_NAMED, default = "schwarzenbacher" ), "pruning_configure": configuration_phil( phil_for = SIDECHAIN_PRUNING_ALGORITHM_NAMED ), "bfactor_use": tbx_utils.choice_string( possible = BFACTOR_ALGORITHM_NAMED, default = "original" ), "bfactor_configure": configuration_phil( phil_for = BFACTOR_ALGORITHM_NAMED ), "completion_use": tbx_utils.choice_string( possible = reduce( operator.add, [ [ a.NAME for a in mm_type.sidechain_completion ] for mm_type in KNOWN ] ), default = "cbeta" ), } PHIL_SCULPTOR = """ macromolecule .help = "Workflow step configuration" { %(protocol_operations)s %(common_operations)s } """ % { "protocol_operations": PHIL_PROTOCOL_OPERATIONS, "common_operations": PHIL_COMMON_OPERATIONS, } PHIL_SCULPTOR_MP = """ macromolecule .help = "Workflow step configuration" { protocols = %(protocols)s .help = "Available protocols (* = active)" .type = choice ( multi = True ) .optional = False %(common_operations)s } """ % { "protocols": tbx_utils.choice_string( possible = get_known_protocols(), default = "1" ), "common_operations": PHIL_COMMON_OPERATIONS, } PHIL_PROTOCOL_MP = """ macromolecule .help = "Workflow step configuration" { %(protocol_operations)s } """ % { "protocol_operations": PHIL_PROTOCOL_OPERATIONS, } PHIL_OTHER = """ hetero = None .help = "Keep named hetero residues" .type = strings .multiple = False min_hssp_length = 6 .help = "Length of residue segment that indicates a reliable match" .type = int .multiple = False .optional = False min_matching_fraction = 0.4 .help = "Minimum matching fraction in residue-to-alignment matching" .type = float .multiple = False .optional = False """ PHIL_MASTER = libtbx.phil.parse( """ %(input)s %(output)s %(sculptor)s %(other)s """ % { "input": PHIL_INPUT, "output": PHIL_OUTPUT, "sculptor": PHIL_SCULPTOR, "other": PHIL_OTHER, }, process_includes=True ) SETUP_DEFAULT = ( PHIL_MASTER, process_user_defined_protocol_phil ) PHIL_MULTIPROTOCOL_MASTER = libtbx.phil.parse( """ %(input)s %(output)s %(sculptor)s %(other)s """ % { "input": PHIL_INPUT, "output": PHIL_OUTPUT, "sculptor": PHIL_SCULPTOR_MP, "other": PHIL_OTHER, }, process_includes=True ) SETUP_MULTIPROTOCOL = ( PHIL_MULTIPROTOCOL_MASTER, process_multiprotocol_phil ) # For documentation master_params = PHIL_MASTER def transfer_solvent_to_separate_chains(root): import iotbx.pdb resnames = set( mmt.SOLVENT.three_letter_codes ) for c in root.chains(): is_solvent = [ i for ( i, rg ) in enumerate( c.residue_groups() ) if rg.atom_groups()[0].resname in resnames ] if not is_solvent or is_solvent[0] == 0: continue solvent_chain = iotbx.pdb.hierarchy.chain() solvent_chain.id = c.id for rg in c.residue_groups()[ is_solvent[0] : ]: c.remove_residue_group( rg ) solvent_chain.append_residue_group( rg ) c.parent().append_chain( solvent_chain ) def remove_alternate_conformations(root): from scitbx.array_family import flex for rg in root.residue_groups(): ags = rg.atom_groups() assert ags if ags[0].altloc == "": if 2 <= len( ags ): ags[1].altloc = "" rg.merge_atom_groups( ags[0], ags[1] ) else: ags[0].altloc = "" for ag in ags[1:]: rg.remove_atom_group( ag ) atoms = root.atoms() atoms.set_occ( flex.double( [ 1.0 ] * len( atoms ) ) ) def sanitize_occupancies(root): from scitbx.array_family import flex atoms = root.atoms() occ = atoms.extract_occ() sanitized = [ min( o, 1.0 ) for o in occ ] atoms.set_occ( flex.double( sanitized ) ) def restrict_model_count(root, max_count): for m in root.models()[max_count:]: root.remove_model( m ) if len( root.models() ) == 1: root.models()[0].id = str( 0 ) def chain_samples(factory, chain_infos, alignments, sequence_for, min_length, min_fraction, logger): from phaser import rsam collection = rsam.SequenceDB() for ali in alignments: for ( index, seq ) in enumerate( ali.object.sequence_strings() ): collection.add( sequence = seq, data = ( index, ali ) ) known = [] unknown = [] for cinfo in chain_infos: logger.info( msg = "Processing: %s" % cinfo ) mm_type = mmt.determine( chain = cinfo.chain ) logger.info( "Chain type: %s" % mm_type.name ) if mm_type in KNOWN: try: cs = factory( chain = cinfo.chain, mmt = mm_type, collection = collection, sequence = sequence_for.get( cinfo.chain.id ), min_length = min_length, min_fraction = min_fraction, logger = logger ) known.append( cs ) except RuntimeError, e: logger.warning( msg = e ) logger.warning( msg = "Treating as UNKNONW type" ) cs = chisel.ChainSample( chain = cinfo.chain, mmt = mm_type ) unknown.append( cs ) else: cs = chisel.ChainSample( chain = cinfo.chain, mmt = mm_type ) unknown.append( cs ) cs.store[ "info" ] = cinfo.identifier[1:] return ( known, unknown ) def unaligned_chain_sample(chain, mmt, collection, sequence, min_length, min_fraction, logger): return chisel.ChainSample( chain = chain, mmt = mm_type ) def aligned_chain_sample(chain, mmt, collection, sequence, min_length, min_fraction, logger): from phaser import rsam from phaser import application_utils logger.info( msg = "Chain-sequence matching:" ) one_letter_for = mmt.one_letter_for() seq = rsam.one_letter_sequence( rgs = chain.residue_groups(), one_letter_for = one_letter_for, unknown = mmt.unknown_one_letter ) logger.debug( msg = "Chain sequence: %s" % seq ) if sequence: logger.info( msg = "Target sequence provided" ) ali_obj = application_utils.create_alignment( seq_a = sequence.sequence, seq_b = seq, name_a = "MASTER", name_b = "CHAIN_%s" % chain.id ) assert ali_obj.multiplicity() == 2 logger.info( "Alignment with target:" ) logger.info( ali_obj ) ali = tbx_utils.AlignmentObject( name = "phenix.alignment", alignment = ali_obj ) index = 1 chain_ali = bioinformatics.clustal_alignment( names = [ "CHAIN" ], alignments = [ seq ] ) else: logger.info( msg = "Select alignment from predefined set" ) if collection.count() == 0: raise RuntimeError, "No alignments are available" logger.debug( msg = "Pre-screening alignment sequences..." ) logger.debug( msg = "Scoring sequences based on exact overlap..." ) matchings = collection.matching( sequence = seq ) assert matchings longest = max( matchings ) if longest < min_length: logger.warning( msg = "The longest exact overlap is shorter than %s residues (min_hssp_length)" % min_length ) raise RuntimeError, "No sufficiently similar alignment sequences have been found" logger.debug( msg = "Longest overlap is %d residues" % longest ) bests = [ d for ( l, d ) in zip( matchings, collection.store ) if l == longest ] if 1 < len( bests ): logger.warning( msg = "There are %d sequences with longest overlap" % len( bests ) ) logger.warning( msg = "The first one will be selected" ) ( index, ali ) = bests[0][1] logger.info( msg = "Best matching: name = %s, index = %s" % ( ali.name, index ) ) logger.debug( msg = ali.object ) logger.info( msg = output.heading( text = "Residue mapping" ) ) cs = chisel.AlignedChainSample( chain = chain, alignment = ali.object, index = index, mmt = mmt, min_hssp_length = min_length, min_matching_fraction = min_fraction ) if cs.unaligned: logger.warning( msg = "Unaligned residues:" ) unknowns = set( cs.unaligned ) midline = "".join( [ "*" if rg in unknowns else " " for rg in chain.residue_groups() ] ) show_ali = bioinformatics.clustal_alignment( alignments = [ seq ], names = [ "PDB" ] ) logger.warning( msg = show_ali.format( aln_width = 60, caption_width = 15, middle_line = midline ) ) mapping = [ ( c, rg ) for ( c, rg ) in zip( cs.model_sequence(), cs.aligned ) if c != cs.alignment.gap ] chain_ali = bioinformatics.clustal_alignment( alignments = [ "".join( p[0] for p in mapping ), "".join( one_letter_for.get( p[1].atom_groups()[0].resname, mmt.unknown_one_letter ) if p[1] else cs.alignment.gap for p in mapping ) ], names = [ "ALI", "PDB" ] ) identity_fraction = chain_ali.identity_fraction() if identity_fraction < 1.0: logger.warning( msg = "Sequence mismatches" ) logger.warning( msg = "Matching fraction: %.2f" % identity_fraction ) logger.warning( chain_ali ) if identity_fraction < min_fraction: raise RuntimeError, "Unable to align: matching fraction < min_matching_fraction" return cs def chain_sample_factory(processors, logger): if all( p[-1].processable( sample_type = chisel.ChainSample ) for p in processors ): logger.info( msg = "No chain-alignment matching necessary" ) factory = unaligned_chain_sample else: assert all( p[-1].processable( sample_type = chisel.AlignedChainSample ) for p in processors ) logger.info( msg = "Chain-alignment matching will be performed" ) factory = aligned_chain_sample return factory def sample_process(indexer, samples, processor, logger): instructions = [] for cs in samples: logger.info( msg = "Processing '%s'..." % str( cs.store[ "info" ] ) ) assert processor.processable( sample_type = cs.__class__ ) logger.info( msg = " Finding modifications..." ) instr = processor.run( sample = cs, logger = logger ) logger.info( msg = " (%d modifications found)" % len( instr ) ) instructions.append( instr ) bp = chisel.Blueprint.from_instructions( instructions = reduce( operator.add, instructions, [] ), indexer = indexer ) return bp def identity_process(indexer, known, logger): identities = [] for cs in known: i = ( chisel.Entity.Chain( cs.chain ).as_descriptor( indexer = indexer ), cs.identity() ) identities.append( i ) return identities def process_single_pdb( pdb, factory, alignments, sequence_for, sculptors, discarder, min_length, min_fraction, prefix, suffix, output_method, logger ): logger.info( msg = output.underlined( text = "Find chains" ) ) ( known, unknown ) = chain_samples( factory = factory, chain_infos = pdb.get_chain_infos(), alignments = alignments, sequence_for = sequence_for, min_length = min_length, min_fraction = min_fraction, logger = logger ) indexer = chisel.ForwardIndexer( root = pdb.object ) logger.info( msg = "Calculating identity records..." ) logger.info( msg = output.blank() ) identities = identity_process( indexer = indexer, known = known, logger = logger ) if unknown: logger.info( msg = output.underlined( text = "Processing hetero chains" ) ) bp_unknown = sample_process( indexer = indexer, samples = unknown, processor = discarder, logger = logger ) if known: logger.info( msg = output.underlined( text = "Processing protein chains" ) ) outfiles = [] for ( name, tag, s ) in sculptors: logger.info( msg = "Running: %s" % name ) logger.info( msg = output.blank() ) bp_known = sample_process( indexer = indexer, samples = known, processor = s, logger = logger ) for d in bp_unknown.data: bp_known.add( datum = d ) outfile = output_method( root = pdb.object, blueprint = bp_known, identities = identities, file_root = prefix, suffix = "%s%s" % ( suffix, tag ), logger = logger ) logger.info( msg = "%s written\n" % outfile ) outfiles.append( outfile ) return outfiles def run(args, out = None, setup = SETUP_DEFAULT, verbosity = output.SingleStream.INFO): logger = output.SingleStream( stream = out, level = verbosity ) logger.info( msg = output.banner( text = "%s version %s" % ( PROGRAM, VERSION ) ) ) ( phil_master, phil_processor ) = setup 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.file_name=%s" ) ) accumulator.register_file_handler( handler = tbx_utils.ExtensionFileHandler( extensions = [ ".pdb", ".ent" ], template = "input.model.file_name=%s" ) ) for argument in args: argument.process( accumulator = accumulator ) 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 logger.info( msg = output.underlined( text = "Requested protocols" ) ) sculptors = phil_processor( params = params, logger = logger ) discarder = get_discarder_object( params = params.hetero ) # Sanity checks if not params.input.model.file_name: raise Sorry, "No pdb files specified" # Main functionality logger.info( msg = output.underlined( "Reading PDB file" ) ) logger.info( msg = "File name: %s" % params.input.model.file_name ) pdb = tbx_utils.PDBObject.from_file( file_name = params.input.model.file_name ) if params.input.model.selection != "all": if not params.input.model.selection: params.input.model.selection = "none" logger.info( msg = "Atom selection: %s" % params.input.model.selection ) asc = pdb.object.atom_selection_cache().selection( params.input.model.selection ) root = pdb.object.select( asc, True ) if not root.atoms(): raise Sorry, "No atoms left after atom selection" pdb = tbx_utils.PDBObject( root = root, name = pdb.name ) if params.input.model.remove_alternate_conformations: logger.info( msg = "Removing alternate conformations..." ) remove_alternate_conformations( root = pdb.object ) if params.input.model.sanitize_occupancies: logger.info( msg = "Removing occupancies > 1.0..." ) sanitize_occupancies( root = pdb.object ) logger.info( msg = "Separating solvent from macromolecules..." ) transfer_solvent_to_separate_chains( root = pdb.object ) # Read in alignments alignments = [] if params.input.alignment: logger.info( msg = output.underlined( text = "Reading alignment files" ) ) for ( index, inp ) in enumerate( params.input.alignment, start = 1 ): logger.info( msg = "%d. %s" % ( index, inp.file_name ) ) p = tbx_utils.AlignmentObject.from_file( file_name = inp.file_name ) if inp.target_index != 1: logger.info( msg = "Assigning sequence #%d as target..." % inp.target_index ) try: p.object.assign_as_target( index = inp.target_index - 1 ) except IndexError, e: raise Sorry, "%s: %d" % ( e, inp.target_index ) alignments.append( p ) logger.verbose( msg = output.heading( text = "Alignment read" ) ) logger.verbose( msg = p.object ) logger.info( msg = output.blank() ) # Read in sequence files sequence_for = {} if params.input.sequence: logger.info( msg = output.underlined( text = "Reading sequence files" ) ) for ( index, inp ) in enumerate(params.input.sequence): logger.info( msg = "%d. %s" % ( index, inp.file_name ) ) p = tbx_utils.SequenceObject.from_file( file_name = inp.file_name ) logger.debug( msg = "Sequences read:" ) for o in p.object: logger.debug( msg = str( o ) ) if 1 < len( p.object ): logger.warning( msg = "File contains multiple sequences, only the first will be used" ) if inp.chain_ids is None: raise Sorry, "Empty chain_id for %s" % inp.file_name applicable = reduce( operator.add, [ s.split( "," ) for s in inp.chain_ids ], [] ) logger.info( "Applicable chains: %s" % ", ".join( [ str( i ) for i in applicable ] ) ) sequence_for.update( dict( [ ( k, p.object[0] ) for k in applicable ] ) ) logger.info( msg = output.blank() ) # Setup values corresponding to all protocols import os.path factory = chain_sample_factory( processors = sculptors, logger = logger ) output_method = OUTPUT_METHOD_FOR[ params.output.format ] logger.info( msg = output.subtitle( text = "Processing PDB file" ) ) logger.info( msg = "File name: %s" % pdb.name ) ( head, tail ) = os.path.split( pdb.name ) try: outfiles = process_single_pdb( pdb = pdb, factory = factory, alignments = alignments, sequence_for = sequence_for, sculptors = sculptors, discarder = discarder, min_length = params.min_hssp_length, min_fraction = params.min_matching_fraction, prefix = os.path.join( params.output.folder, os.path.basename( params.output.root ) ), suffix = os.path.splitext( tail )[0], output_method = output_method, logger = logger ) except RuntimeError, e: raise Sorry, e return outfiles class launcher (runtime_utils.target_with_save_result) : def run (self) : return run_phil(args=self.args, out=sys.stdout) def run_phil (args, out=sys.stdout) : factory = tbx_utils.PhilArgumentFactory( master_phil = PHIL_MASTER ) return run( args=[ factory( argument = arg ) for arg in args ], out=out) def finish_job (result) : import os.path files = [] for file_name in result : if os.path.isfile(file_name) : pdb_out = os.path.abspath(file_name) files.append((pdb_out, "Final model")) return (files, []) def validate_params (params) : if (None in [params.input.model, params.input.alignment]) : raise Sorry("Missing the PDB file or alignment file.") return True