You are reading the documentation for version 1.6 of OpenStructure. You may also want to read the documentation for: 1.1 1.2 1.3 1.4 1.5 devel

# `seq.alg` – Algorithms for Sequences¶

## Algorithms for Alignments¶

`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: aln (`AlignmentHandle`) – Alignment chain (`ChainHandle`) – Source of the sequence 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()` The alignment of the residues in the chain and the SEQRES entries. `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 The alignment `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 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 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. 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 best-scoring alignment of seq1 and seq2.
`Renumber`(seq_handle, sequence_number_with_attached_view=1)

Function to renumber an entity according to an alignment between the model sequence and the full-length target sequence. The aligned model sequence or the alignment itself with an attached view needs to be provided. Upon succcess, the renumbered entity is returned. If an alignment is given, the sequence must

```from ost.seq.alg import renumber
from ost.bindings.clustalw import *
pdb_seq = seq.SequenceFromChain("model", ent.chains[0])
aln = ClustalW(s, pdb_seq)
aln.AttachView(1, ent.chains[0].Select(""))
e = Renumber(aln.sequences[1])
io.SavePDB(e, "renum.pdb")
```
Parameters: seq_handle (`SequenceHandle` / `AlignmentHandle`) – Sequence or alignment handle with attached view. sequence_number_with_attached_view (`int`) – Sequence number for the aln. handle (not used if seq. handle given) `RuntimeError` if unknown type of seq_handle or if attached view is missing or if the given alignment sequence is inconsistent.

## 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.

## Get and analyze distance matrices from alignments¶

Given a multiple sequence alignment between a reference sequence (first sequence in alignment) and a list of structures (remaining sequences in alignment with an attached view to the structure), this set of functions can be used to analyze differences between the structures.

Example:

```# SETUP: aln is multiple sequence alignment, where first sequence is the
#        reference sequence and all others have a structure attached

# clip alignment to only have parts with at least 3 sequences (incl. ref.)
# -> aln will be cut and clip_start is 1st column of aln that was kept
clip_start = seq.alg.ClipAlignment(aln, 3)

# get variance measure and distance to mean for each residue pair
d_map = seq.alg.CreateDistanceMap(aln)
var_map = seq.alg.CreateVarianceMap(d_map)
dist_to_mean = seq.alg.CreateDist2Mean(d_map)

# report min. and max. variances
print "MIN-MAX:", var_map.Min(), "-", var_map.Max()
# get data and json-strings for further processing
var_map_data = var_map.GetData()
var_map_json = var_map.GetJsonString()
dist_to_mean_data = dist_to_mean.GetData()
dist_to_mean_json = dist_to_mean.GetJsonString()
```
`ClipAlignment`(aln, n_seq_thresh=2, set_offset=true, remove_empty=true)

Clips alignment so that first and last column have at least the desired number of structures.

Parameters: aln (`AlignmentHandle`) – Multiple sequence alignment. Will be cut! n_seq_thresh (`int`) – Minimal number of sequences desired. set_offset (`bool`) – Shall we update offsets for attached views? remove_empty (`bool`) – Shall we remove sequences with only gaps in cut aln? Starting column (0-indexed), where cut region starts (w.r.t. original aln). -1, if there is no region in the alignment with at least the desired number of structures. `int`
`CreateDistanceMap`(aln)

Create distance map from a multiple sequence alignment.

The algorithm requires that the sequence alignment consists of at least two sequences. The sequence at index 0 serves as a frame of reference. All the other sequences must have an attached view and a properly set sequence offset (see `SetSequenceOffset()`).

For each of the attached views, the C-alpha distance pairs are extracted and mapped onto the corresponding C-alpha distances in the reference sequence.

Parameters: aln (`AlignmentHandle`) – Multiple sequence alignment. Distance map. `DistanceMap` Exception if aln has less than 2 sequences or any sequence (apart from index 0) is lacking an attached view.
`CreateVarianceMap`(d_map, sigma=25)
Returns: Variance measure for each entry in d_map. `VarianceMap` d_map (`DistanceMap`) – Distance map as created with `CreateDistanceMap()`. sigma (`float`) – Used for weighting of variance measure (see `Distances.GetWeightedStdDev()`) Exception if d_map has no entries.
`CreateDist2Mean`(d_map)
Returns: Distances to mean for each structure in d_map. Structures are in the same order as passed when creating d_map. `Dist2Mean` d_map (`DistanceMap`) – Distance map as created with `CreateDistanceMap()`. Exception if d_map has no entries.
class `Distances`

Container used by `DistanceMap` to store a pair wise distance for each structure. Each structure is identified by its index in the originally used alignment (see `CreateDistanceMap()`).

`GetDataSize`()
Returns: Number of pairwise distances. `int`
`GetAverage`()
Returns: Average of all distances. `float` Exception if there are no distances.
`GetMin`()
`GetMax`()
Returns: Minimal/maximal distance. `tuple` (distance (`float`), index (`int`)) Exception if there are no distances.
`GetDataElement`(index)
Returns: Element at given index. `tuple` (distance (`float`), index (`int`)) index (`int`) – Index within list of distances (must be < `GetDataSize()`). Exception if there are no distances or index out of bounds.
`GetStdDev`()
Returns: Standard deviation of all distances. `float` Exception if there are no distances.
`GetWeightedStdDev`(sigma)
Returns: Standard deviation of all distances multiplied by exp( `GetAverage()` / (-2*sigma) ). `float` sigma (`float`) – Defines weight. Exception if there are no distances.
`GetNormStdDev`()
Returns: Standard deviation of all distances divided by `GetAverage()`. `float` Exception if there are no distances.
class `DistanceMap`

Container returned by `CreateDistanceMap()`. Essentially a symmetric `GetSize()` x `GetSize()` matrix containing up to `GetNumStructures()` distances (list stored as `Distances`). Indexing of residues starts at 0 and corresponds to the positions in the originally used alignment (see `CreateDistanceMap()`).

`GetDistances`(i_res1, i_res2)
Returns: List of distances for given pair of residue indices. `Distances` i_res1 (`int`) – Index of residue. i_res2 (`int`) – Index of residue.
`GetSize`()
Returns: Number of residues in map. `int`
`GetNumStructures`()
Returns: Number of structures originally used when creating the map (see `CreateDistanceMap()`). `int`
class `VarianceMap`

Container returned by `CreateVarianceMap()`. Like `DistanceMap`, it is a symmetric `GetSize()` x `GetSize()` matrix containing variance measures. Indexing of residues is as in `DistanceMap`.

`Get`(i_res1, i_res2)
Returns: Variance measure for given pair of residue indices. `float` i_res1 (`int`) – Index of residue. i_res2 (`int`) – Index of residue.
`GetSize`()
Returns: Number of residues in map. `int`
`Min`()
`Max`()
Returns: Minimal/maximal variance in the map. `float`
`ExportDat`(file_name)
`ExportCsv`(file_name)
`ExportJson`(file_name)

Write all variance measures into a file. The possible formats are:

• “dat” file: a list of “i_res1+1 i_res2+1 variance” lines
• “csv” file: a list of ”;” separated variances (one line for each i_res1)
• “json” file: a JSON formatted file (see `GetJsonString()`)
Parameters: file_name (`str`) – Path to file to be created. Exception if the file cannot be opened for writing.
`GetJsonString`()
Returns: A JSON formatted list of `GetSize()` lists with `GetSize()` variances `str`
`GetData`()

Gets all the data in this map at once. Note that this is much faster (10x speedup observed) than parsing `GetJsonString()` or using `Get()` on each element.

Returns: A list of `GetSize()` lists with `GetSize()` variances. `list` of `list` of `float`
class `Dist2Mean`

Container returned by `CreateDist2Mean()`. Stores distances to mean for `GetNumResidues()` residues of `GetNumStructures()` structures. Indexing of residues is as in `DistanceMap`. Indexing of structures goes from 0 to `GetNumStructures()` - 1 and is in the same order as the structures in the originally used alignment.

`Get`(i_res, i_str)
Returns: Distance to mean for given residue and structure indices. `float` i_res (`int`) – Index of residue. i_str (`int`) – Index of structure.
`GetNumResidues`()
Returns: Number of residues. `int`
`GetNumStructures`()
Returns: Number of structures. `int`
`ExportDat`(file_name)
`ExportCsv`(file_name)
`ExportJson`(file_name)

Write all distance measures into a file. The possible formats are:

• “dat” file: a list of “i_res+1 distances” lines (distances are space separated)
• “csv” file: a list of ”;” separated distances (one line for each i_res)
• “json” file: a JSON formatted file (see `GetJsonString()`)
Parameters: file_name (`str`) – Path to file to be created. Exception if the file cannot be opened for writing.
`GetJsonString`()
Returns: A JSON formatted list of `GetNumResidues()` lists with `GetNumStructures()` distances. `str`
`GetData`()

Gets all the data in this map at once. Note that this is much faster (10x speedup observed) than parsing `GetJsonString()` or using `Get()` on each element.

Returns: A list of `GetNumResidues()` lists with `GetNumStructures()` distances. `list` of `list` of `float`

## Search

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

## Previous topic

`seq` – Sequences and Alignments

## Next topic

`bindings` – Interfacing external programs