This document is for OpenStructure version 1.5, the latest version is 2.7 !

seq.alg – Algorithms for Sequences

MergePairwiseAlignments(pairwise_alns, ref_seq)

Merge a list of pairwise alignments into a multiple sequence alignments. This function uses the reference sequence as the anchor and inserts gaps where needed. This is also known as the star method.

The resulting multiple sequence alignment provides a simple way to map between residues of pairwise alignments, e.g. to compare distances in two structural templates.

There are a few things to keep in mind when using this function:

  • The reference sequence mustn’t contain any gaps
  • The first sequence of each pairwise alignments corresponds to the reference sequence. Apart from the presence of gaps, these two sequences must be completely identical.
  • If the reference sequence has an offset, the first sequence of each pairwise alignment must have the same offset. This offset is inherited by the first sequence of the final output alignment.
  • The resulting multiple sequence alignment is by no means optimal. For better results, consider using a multiple-sequence alignment program such as MUSCLE or ClustalW.
  • Residues in columns where the reference sequence has gaps should not be considered as aligned. There is no information in the pairwise alignment to guide the merging, the result is undefined.
ValidateSEQRESAlignment(aln, chain=None)

Checks a sequence aligned to a SEQRES sequence to be free of strand breaks. Residues divided by gaps are not considered as breakage but may also not be connected.

Parameters:
Returns:

True if all residues (beside gapped ones) are connected, False otherwise.

AlignToSEQRES(chain, seqres, try_resnum_first=False, validate=True)

Aligns the residues of chain to the SEQRES sequence, inserting gaps where needed. The function uses the connectivity of the protein backbone to find consecutive peptide fragments. These fragments are then aligned to the SEQRES sequence.

All the non-ligand, peptide-linking residues of the chain must be listed in SEQRES. If there are any additional residues in the chain, the function raises a ValueError.

If ‘try_resnum_first’ is set, building the alignment following residue numbers is tried first.

If ‘validate’ is set (default), the alignment is checked using ValidateSEQRESAlignment().

Parameters:
  • chain (ChainHandle) – Source of the sequence
  • seqres (str) – SEQRES sequence
  • try_resnum_first (bool) – Try to align by residue number
  • validate (bool) – Validate alignment by ValidateSEQRESAlignment()
Returns:

The alignment of the residues in the chain and the SEQRES entries.

Return type:

AlignmentHandle

AlignmentFromChainView(chain, handle_seq_name='handle', view_seq_name='view')

Creates and returns the sequence alignment of the given chain view to the chain handle. The alignment contains two sequences, the first containing all non-ligand peptide-linking residues, the second containing all non-ligand peptide-linking residues that are part of the view.

Parameters:
  • chain (ChainView) – A valid chain
  • handle_seq_name – Name of the handle sequence in the output alignment
  • view_seq_name – Name of the view sequence in the output alignment
Returns:

The alignment

Return type:

AlignmentHandle

Conservation(aln, assign=true, prop_name="cons", ignore_gap=false)

Calculates conservation scores for each column in the alignment, according to the ConSurf method (Armon et al., J. Mol. Biol. (2001) 307, 447-463).

The conservation score is a value between 0 and 1. The bigger the number the more conserved the aligned residues are.

Parameters:
  • aln (AlignmentHandle) – An alignment handle
  • assign – If true, the conservation scores are assigned to attached residues. The name of the property can be changed with the prop_name parameter. Useful when coloring entities based on sequence conservation.
  • prop_name – The property name for assigning the conservation to attached residues. Defaults to ‘cons’.
  • ignore_gap – If true, the dissimilarity between two gaps is increased to 6.0 instead of 0.5 as defined in the original version. Without this, a stretch where in the alignment there is only one sequence which is aligned to only gaps, is considered highly conserved (depending on the number of gap sequences).
LocalAlign(seq1, seq2, subst_weight, gap_open=-5, gap_ext=-2)

Performs a Smith/Waterman local alignment of seq1 and seq2 and returns the best-scoring alignments as a list of pairwise alignments.

Example:

seq_a = seq.CreateSequence('A', 'acdefghiklmn')
seq_b = seq.CreateSequence('B', 'acdhiklmn')
alns = seq.alg.LocalAlign(seq_a, seq_b, seq.alg.BLOSUM62)
print alns[0].ToString(80)
# >>> A acdefghiklmn
# >>> B acd---hiklmn
Parameters:
  • seq1 (ConstSequenceHandle) – A valid sequence
  • seq2 (ConstSequenceHandle) – A valid sequence
  • subst_weigth – The substitution weights matrix
  • gap_open – The gap opening penalty. Must be a negative number
  • gap_ext – The gap extension penalty. Must be a negative number
Returns:

A list of best-scoring, non-overlapping alignments of seq1 and seq2. Since alignments always start with a replacement, the start is stored in the sequence offset of the two sequences.

GlobalAlign(seq1, seq2, subst_weight, gap_open=-5, gap_ext=-2)

Performs a Needleman/Wunsch global alignment of seq1 and seq2 and returns the best-scoring alignment.

Example:

seq_a = seq.CreateSequence('A', 'acdefghiklmn')
seq_b = seq.CreateSequence('B', 'acdhiklmn')
alns = seq.alg.GlobalAlign(seq_a, seq_b, seq.alg.BLOSUM62)
print alns[0].ToString(80)
# >>> A acdefghiklmn
# >>> B acd---hiklmn
Parameters:
  • seq1 (ConstSequenceHandle) – A valid sequence
  • seq2 (ConstSequenceHandle) – A valid sequence
  • subst_weigth – The substitution weights matrix
  • gap_open – The gap opening penalty. Must be a negative number
  • gap_ext – The gap extension penalty. Must be a negative number
Returns:

Best-scoring alignment of seq1 and seq2.

ShannonEntropy(aln, ignore_gaps=True)

Returns the per-column Shannon entropies of the alignment. The entropy describes how conserved a certain column in the alignment is. The higher the entropy is, the less conserved the column. For a column with no amino aids, the entropy value is set to NAN.

Parameters:
  • aln (AlignmentHandle) – Multiple sequence alignment
  • ignore_gaps (bool) – Whether to ignore gaps in the column.
Returns:

List of column entropies

SemiGlobalAlign(seq1, seq2, subst_weight, gap_open=-5, gap_ext=-2)

Performs a semi-global alignment of seq1 and seq2 and returns the best- scoring alignment. The algorithm is Needleman/Wunsch same as GlobalAlign, but without any gap penalty for starting or ending gaps. This is prefereble whenever one of the sequences is significantly shorted than the other. This make it also suitable for fragment assembly.

Example:

seq_a=seq.CreateSequence('A', 'abcdefghijklmnok')
seq_b=seq.CreateSequence('B', 'cdehijk')
alns=seq.alg.GlobalAlign(seq_a, seq_b, seq.alg.BLOSUM62)
print alns[0].ToString(80)
# >>> A abcdefghijklmnok
# >>> B --cde--hijk-----
Parameters:
  • seq1 (ConstSequenceHandle) – A valid sequence
  • seq2 (ConstSequenceHandle) – A valid sequence
  • subst_weigth – The substitution weights matrix
  • gap_open – The gap opening penalty. Must be a negative number
  • gap_ext – The gap extension penalty. Must be a negative number
Returns:

best-scoring alignment of seq1 and seq2.

Contact Prediction

This is a set of functions for predicting pairwise contacts from a multiple sequence alignment (MSA). The core method here is mutual information which uses coevolution to predict contacts. Mutual information is complemented by two other methods which score pairs of columns of a MSA from the likelyhood of certain amino acid pairs to form contacts (statistical potential) and the likelyhood of finding certain substitutions of aminio-acid pairs in columns of the MSA corresponding to interacting residues.

class ContactPredictionScoreResult

Object containing the results form a contact prediction.

matrix

An NxN FloatMatrix where N is the length of the alignment. The element i,j corresponds to the score of the corresponding columns of the MSA. High scores correspond to high likelyhood of a contact.

sorted_indices

List of all indices pairs i,j, containing (N*N-1)/2 elements, as the matrix is symmetrical and elements in the diagonal are ignored. The indices are sorted from the pair most likely to form a contact to the least likely one.

GetScore(i, j)

returns matrix(i,j)

Parameters:
  • i (int) – First index
  • j (int) – Second index
SetScore(i, j, score)

Sets matrix(i,j) to score

Parameters:
  • i (int) – First index
  • j (int) – Second index
  • score (float) – The score
PredictContacts(ali)

Predicts contacts from a multiple sequence alignment using a combination of Mutual Information (MI) and the Contact Substitution Score (CoEvoSc). MI is calculated with the APC and small number corrections as well as with a transformation into Z-scores. The CoEvoSc is calculated using the default PairSubstWeightMatrix (see seq.alg.LoadDefaultPairSubstWeightMatrix). The final score for a pair of columns (i,j) of ali is obtained from:

Sc(i,j)=MI(i,j)exp(CoEvoSc(i,j)) if (i,j) >=0

Sc(i,j)=MI(i,j)exp(1-CoEvoSc(i,j)) if (i,j) <0

Parameters:ali (AlignmentHandle) – The multiple sequence alignment
CalculateMutualInformation(aln,weights=LoadConstantContactWeightMatrix(),
apc_correction=true,zpx_transformation=true,small_number_correction=0.05)

Calculates the mutual information (MI) from a multiple sequence alignemnt. Contributions of each pair of amino-acids are weighted using the matrix weights (weighted mutual information). The average product correction (apc_correction) correction and transformation into Z-scores (zpx_transofrmation) increase prediciton accuracy by reducing the effect of phylogeny and other noise sources. The small number correction reduces noise for alignments with small number of sequences of low diversity.

Parameters:
  • aln (AlignmentHandle) – The multiple sequences alignment
  • weights (:class`ContactWeightMatrix`) – The weight matrix
  • apc_correction (bool) – Whether to use the APC correction
  • zpx_transformation (bool) – Whether to transform the scores into Z-scores
  • small_number_correction (float) – initial values for the probabilities of having a given pair of amino acids p(a,b).
CalculateContactScore(aln, weights=LoadDefaultContactWeightMatrix())

Calculates the Contact Score (CoSc) from a multiple sequence alignment. For each pair of residues (i,j) (pair of columns in the MSA), CoSc(i,j) is the average over the values of the weights corresponding to the amino acid pairs in the columns.

Parameters:
  • aln (AlignmentHandle) – The multiple sequences alignment
  • weights (:class`ContactWeightMatrix`) – The contact weight matrix
CalculateContactSubstitutionScore(aln,ref_seq_index=0,
weights=LoadDefaultPairSubstWeightMatrix())

Calculates the Contact Substitution Score (CoEvoSc) from a multiple sequence alignment. For each pair of residues (i,j) (pair of columns in the MSA), CoEvoSc(i,j) is the average over the values of the weights corresponding to substituting the amino acid pair in the reference sequence (given by ref_seq_index) with all other pairs in columns (i,j) of the aln.

Parameters:
  • aln (AlignmentHandle) – The multiple sequences alignment
  • weights (:class`ContactWeightMatrix`) – The pair substitution weight matrix
LoadDefaultContactWeightMatrix()
Returns:CPE, a ContactWeightMatrix that was calculated from a large (>15000) set of high quality crystal structures as CPE=log(CF(a,b)/NCF(a,b)) and then normalised so that all its elements are comprised between 0 and 1. CF(a,b) is the frequency of amino acids a and b for pairs of contacting residues and NCF(a,b) is the frequency of amino acids a and b for pairs of non-contacting residues. Apart from weights for the standard amino acids, this matrix gives a weight of 0 to all pairs for which at least one amino-acid is a gap.
LoadConstantContactWeightMatrix()
Returns:A ContactWeightMatrix. This matrix gives a weight of one to all pairs of standard amino-acids and a weight of 0 to pairs for which at least one amino-acid is a gap.
LoadDefaultPairSubstWeightMatrix()
Returns:CRPE, a PairSubstWeightMatrix that was calculated from a large (>15000) set of high quality crystal structures as CRPE=log(CRF(ab->cd)/NCRF(ab->cd)) and then normalised so that all its elements are comprised between 0 and 1. CRF(ab->cd) is the frequency of replacement of a pair of amino acids a and b by a pair c and d in columns of the MSA corresponding to contacting residues and NCRF(ab->cd) is the frequency of replacement of a pair of amino acids a and b by a pair c and d in columns of the MSA corresponding to non-contacting residues. Apart from weights for the standard amino acids, this matrix gives a weight of 0 to all pair substitutions for which at least one amino-acid is a gap.
class PairSubstWeightMatrix(weights, aa_list)

This class is used to associate a weight to any substitution from one amino-acid pair (a,b) to any other pair (c,d).

weights

A FloatMatrix4 of size NxNxNxN, where N=len(aa_list)

aa_list

A CharList of one letter codes of the amino acids for which weights are found in the weights matrix.

class ContactWeightMatrix(weights, aa_list)

This class is used to associate a weight to any pair of amino-acids.

weights

A FloatMatrix of size NxN, where N=len(aa_list)

aa_list

A CharList of one letter codes of the amino acids for which weights are found in the weights matrix.

Search

Enter search terms or a module, class or function name.

Contents

Documentation is available for the following OpenStructure versions:

dev / 2.7 / 2.6 / 2.5 / 2.4 / 2.3.1 / 2.3 / 2.2 / 2.1 / 2.0 / 1.9 / 1.8 / 1.7.1 / 1.7 / 1.6 / (Currently viewing 1.5) / 1.4 / 1.3 / 1.2 / 1.11 / 1.10 / 1.1

This documentation is still under heavy development!
If something is missing or if you need the C++ API description in doxygen style, check our old documentation for further information.