Author(s)

• sequence_from_map: Tom Terwilliger

Purpose

The routine sequence_from_map will use a map to guess the sequence of a supplied segment or of segments built automatically.

Usage

How sequence_from_map works:

Sequence-from-map examines the features in a map at the position where side chains in a model are expected to be located and uses these features to guess the identity of the side chain at each position.

The procedure used by sequence-from-map is has several steps:

If map (and no model) is supplied, a model is created by searching for regular
secondary structure in the map.  First the unique part of the map is
identified and cut out to form a new map.
This is done first using the map_symmetry to find helical or point-group
symmetry in the map.  If any is found, a new map is created containing just
the unique part of the map. Then (if improper_symmetry is set), the new map
is examined to see if there is any local symmetry remaining. If so, the unique
part is cut out and a new map is created. Then the tool find_helices_strands
is used to analyze the new map to find regular secondary structure.

Once a map and model are available, each segment (a piece of chain without
breaks) is compared to the map and sequences compatible with the map
are identified.

The core of sequence-from-map is to generate a scoring matrix for each
segment that reflects the relative probability that each possible side chain
is located at each position in the segment.  This probability is estimated
from the map-model correlation for each side chain (after examining all
rotamers of the side chain and picking the one with the highest correlation).
A Z-score (value for this side chain minus mean for this side chain at all
positions, divided by the standard deviation of this quantity) is calculated
for each side chain at each position, where Z values less than zero are set
to zero.  These Z-scores are used as approximate minus-log-likelihood scores
for side chain probabilities and make up the scoring matrix for the segment.

A best-fitting sequence is generated for each segment by simply picking the
highest-scoring side chain at each position.  Optionally, if
score_by_residue_groups is set and residue groups are defined
(residue_groups = "VGASCTI P LDNEQM KR FHY W"), one representative of each
residue group is chosen to represent all the side chains in that group.

Optionally (if minimum_discrimination is defined and trim_models is set),
the scoring matrix is used to cut the model up into pieces, removing
any residues where the discrimination between the most probable and
least probable side chain is less than minimum_discrimination.  The
segments are then cut up at these points and smaller segments are created.

Also optionally (if trim_models is set to some value N), the N residues on
each end of each segment are removed before analysis.

If one or more potential sequences representing the molecule in the map are
provided, these sequences are scored, one at a time, by evaluating the
optimal alignment of the sequence to the segments in the model.

The alignment of a sequence to a set of segments is done by considering
all possible alignments for all segments, keeping the highest-scoring
segment-sequence alignment, then repeating the process with all remaining
segments.  Optionally, the alignments can be restricted to a set that
uses each residue in the sequence only once (allow_duplicates set to False).
Optionally, sequences that are shorter than the number of residues in the
model are rejected (if skip_if_too_short is set).

The raw score for one segment-sequence alignment is given by the sum of the
position-dependent minus-log-likelihood scores for all the side chains
in the alignment. This raw score is then adjusted by subtracting the mean
value of raw segment-sequence alignment scores for random sequences of the
same length as the sequence being considered.  The random sequences are
created using a residue frequency of eukaryotic amino acids (this can be
adjusted with the default_sequence parameter).  If positive_only is set,
negative segment-sequence alignment scores are ignored.

Once the best-fitting sequence is identified, the side chains of the
model are re-fit using this sequence and a model with the fitted side
chains is written out.

Examples

Standard run of sequence_from_map:

You can use sequence_from_map to guess the sequence of a fragment you have built:

phenix.sequence_from_map my_fragment.pdb my_map.mrc resolution=2.8

You can also use it to try and pick out which sequence file matches your map:

phenix.sequence_from_map my_map.mrc resolution=3 seq_dir=sequences/

where the directory sequences has a set of fasta or similar files representing possible sequences that might match your map.

Possible Problems

Specific limitations and problems:

Literature

• Automated side-chain model building and sequence assignment by template matching. T.C. Terwilliger. Acta Crystallogr D Biol Crystallogr 59, 45-9 (2003).

• Rapid model building of alpha-helices in electron-density maps. T.C. Terwilliger. Acta Crystallogr D Biol Crystallogr 66, 268-75 (2010).

• Rapid model building of beta-sheets in electron-density maps. T.C. Terwilliger. Acta Crystallogr D Biol Crystallogr 66, 276-84 (2010).

Additional information

List of all available keywords

• job_title = None Job title in PHENIX GUI, not used on command line
• input_files
  • map_file = None File with CCP4-style map. May have origin in any location.
  • multiple_seq_file = None Optional sequence file containing multiple potential sequences
  • seq_file = None Optional sequence file(s) containing potential sequences
  • seq_dir = None Optional directory containing sequence files
  • model_file = None Input PDB file with chains to be sequenced
  • symmetry_file = None Optional symmetry file (.ncs_spec format or MATRX records) with reconstruction symmetry. Used in identification of unique part of map.
  • probability_matrix_file = None Optional .pkl file with probability matrices. Note: if build_type is all then a prefix of either helices_strands or full_model will be added to the file name (if not already there)
  • model_pickle_file = None Optional .pkl file with final fragments. Note: if build_type is all then a prefix of either helices_strands or full_model will be added to the file name (if not already there)
• output_files
  • pdb_out = sequence_from_map.pdb Model with sequence determined from map. Note: if build_type is all then a prefix of either helices_strands or full_model will be added to the file name (if not already there)
  • output_probability_matrix_file = None Optional output .pkl file with probability matrices
  • output_model_pickle_file = None Optional output .pkl file with final fragments
  • temp_dir = temp_dir Temporary directory
• crystal_info
  • resolution = None High-resolution limit for map analysis.
  • scattering_table = n_gaussian wk1995 it1992 *electron neutron Choice of scattering table for structure factor calculations. Standard for X-ray is n_gaussian, for cryoEM is electron.
  • sequences = None Sequences to be matched as text strings
  • default_sequence = None Sequence reflecting relative usage of residues/nucleotides. Default is amino acid composition for eukaryotes.
  • chain_type = *PROTEIN DNA RNA Chain type
  • molecular_mass = None Molecular mass of the entire structure in the map (not the mass of one chain). Does not have to be accurate, but helpful if using symmetry to extract unique part of map.
• sequencing
  • build_type = *all helices_strands full_model You can choose to build just helices/strands (quick) or a full model (takes more time). Note you can define how many regions to analyze with regions_to_keep (set it to 1 to run quickly). You can choose all to run both helices_strands and full_model.
  • regions_to_keep = None You can specify a limit to the number of regions to keep when generating the asymmetric unit of density. Default is 1 for build_type=full_model and None for build_type=helices_strands.
  • minimum_final_z_score = 4 Skip additional builds if a solution with this Z-score or better is found.
  • minimum_length = 10 Minimum length of a segment to consider it for sequencing
  • maximum_length = 50 If segment is longer than this, try to cut it into pieces, splitting at residues that are poorly-defined. Assumption is very long segments probably have register shifts.
  • maximum_segments = 10 Maximum number of segments to consider (remainder ignored)
  • maximum_sequence_length = 2000 Maximum sequence length to consider.
  • subtract_expected = None *single_try multiple_tries Subtract expected score based on random sequence for a single try or for the number of tries carried out.
  • positive_only = False Just sum scores that are positive and skip others
  • allow_duplicates = True Allow duplication of use of sequence.
  • skip_if_too_short = False Skip sequences that are shorter than the number of residues in the model.
  • refinement = True Run refinement after finding model
  • macro_cycles = 5 Refinement macro-cycles
  • trim_from_ends = None You can take final fragments and trim them N residues in from each end.
  • residue_groups = "VGASCTI P LDNEQM KR FHY W" Optional groups of residues to score together
  • score_by_residue_groups = False Use residue groups in final scoring of side chains
  • trim_models = False Trim models to break at residues with low discrimination. Default is False.
  • vary_minimum_cc = False Vary the minimum_cc in optimizing percent to reject. Alternative is vary_minimum_z_score.
  • percent_to_reject = 7.5 Adjust minimum_cc and minimum_z_score to reject this percentage of residues.
  • minimum_z_score = 0.5 Minimum residue Z-score for a residue (or residue groups if score_by_residue_groups is set) to keep the choice. Otherwise report residue as X.
  • minimum_cc = None Minimum map-model correlation for a residue (or residue groups if score_by_residue_groups is set) to keep the choice. Otherwise report residue as X.
  • minimum_discrimination = 0.01 Minimum discrimination between lowest and highest probability residue (or residue groups if score_by_residue_groups is set) to keep the choice. Otherwise report residue as X.
  • random_sequences = 10 Number of random sequences of each length to use as baseline
  • use_default_sequence = False Base the random sequences on the default sequence. Alternative is to base them on the actual sequence
  • symmetry = ANY Optional symmetry used in reconstruction. For example D7, C3, C2 I (icosahedral),T (tetrahedral), or ANY (try everything and use the highest symmetry found). Not needed if symmetry is supplied.
  • include_helical_symmetry = False Include search for helical symmetry if symmetry=ALL
  • improper_symmetry = False After any proper symmetry is found, take new unique part of map and search for any symmetry at all with find_ncs_from_density. Useful if your molecule has both proper reconstruction symmetry such as C3 and within the reconstruction unit there is local symmetry which can be of any kind. Final result is just the unique part of the map.
  • separate_chains = True Write each segment as a separate chain with unique chain ID
  • null_score = -2 Score for a segment that is not placed
  • n_direct = 10 Extreme values for N<= n_direct calculated directly
• control
  • nproc = 1 Number of processors to use
  • random_seed = 171731 Random seed
  • verbose = False Verbose output
  • quick = False Quick run
• guiGUI-specific parameter required for output directory
  • output_dir = None