alg
– Algorithms for Sequences¶
Submodules¶
Algorithms for Alignments¶
- MergePairwiseAlignments(pairwise_alns, ref_seq)¶
- Parameters:
pairwise_alns (
AlignmentList
) – A list of pairwise alignmentsref_seq (
SequenceHandle
) – The reference sequence
- Returns:
The merged alignment
- Return type:
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.
Example:
ref_seq = ost.seq.CreateSequence('ref', 'acdefghiklmn') seq_a1 = seq.CreateSequence('A1', 'acdefghikl-mn') seq_a2 = seq.CreateSequence('A2', 'atd-fghikllmn') seq_b1 = seq.CreateSequence('B1', 'acdefg-hiklmn') seq_b2 = seq.CreateSequence('B2', 'acd---qhirlmn') aln_a = seq.CreateAlignment() aln_a.AddSequence(seq_a1) aln_a.AddSequence(seq_a2) print(aln_a) # >>> A1 acdefghikl-mn # >>> A2 atd-fghikllmn aln_b = seq.CreateAlignment() aln_b.AddSequence(seq_b1) aln_b.AddSequence(seq_b2) print(aln_b) # >>> B1 acdefg-hiklmn # >>> B2 acd---qhirlmn aln_list = ost.seq.AlignmentList() aln_list.append(aln_a) aln_list.append(aln_b) merged_aln = ost.seq.alg.MergePairwiseAlignments(aln_list, ref_seq) print(merged_aln) # >>> ref acdefg-hikl-mn # >>> A2 atd-fg-hikllmn # >>> B2 acd---qhirl-mn
- ValidateSEQRESAlignment(aln, chain=None)¶
Checks if sequence in alignment has same connectivity as residues in chain. This looks for connected stretches in both the sequence and the chain and returns False if they don’t match. This uses the connectivity of the protein backbone.
- Parameters:
aln (
AlignmentHandle
) – Alignment of two sequences with second one expected to map to residues in chain.chain (
ChainHandle
) – Source of the sequence.
- 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.
- Parameters:
chain (
ChainHandle
) – Source of the sequenceseqres (
str
) – SEQRES sequencetry_resnum_first (
bool
) – If set to True, this first builds an alignment using residue numbers and checks if the one-letter-codes match. If they all match, this alignment is used (and possibly validated). Otherwise, it displays a warning and falls back to the connectivity-based alignment.validate (
bool
) – If set to True, the alignment is additionally checked byValidateSEQRESAlignment()
and raises a ValueError if the validation failed.
- Returns:
The alignment of the residues in the chain and the SEQRES entries.
- Return type:
- 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 chainhandle_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:
- 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 handleassign – 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).
- 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 alignmentignore_gaps (bool) – Whether to ignore gaps in the column.
- Returns:
List of column entropies
- SequenceIdentity(aln, ref_mode=seq.alg.RefMode.ALIGNMENT, seq_a=0, seq_b=1)¶
Calculates the sequence identity between two sequences at index seq_a and seq_b in a multiple sequence alignment.
- Parameters:
aln (
AlignmentHandle
) – multiple sequence alignmentref_mode (int) – influences the way the sequence identity is calculated. When set to seq.alg.RefMode.LONGER_SEQUENCE, the sequence identity is calculated as the number of matches divided by the length of the longer sequence. If set to seq.alg.RefMode.ALIGNMENT (the default), the sequence identity is calculated as the number of matches divided by the number of aligned residues (not including the gaps).
seq_a (int) – the index of the first sequence
seq_b (int) – the index of the second sequence
- Returns:
sequence identity in the range 0 to 100.
- Return type:
float
- SequenceSimilarity(aln, subst_weight, normalize=false, seq_a=0, seq_b=1)¶
Calculates the sequence similarity between two sequences at index seq_a and seq_b in a multiple sequence alignment.
- Parameters:
aln (
AlignmentHandle
) – Multiple sequence alignmentsubst_weight (
SubstWeightMatrix
) – the substitution weight matrix (see the BLOSUM Matrix section below)normalize (bool) – if set to True, normalize to the range of the substitution weight matrix
seq_a (int) – the index of the first sequence
seq_b (int) – the index of the second sequence
- Returns:
sequence similarity
- Return type:
float
Create pairwise alignments¶
OpenStructure provides naive implementations to create pairwise local, global and semi-global alignments between two sequences:
The use of parasail as a drop in replacement is optional and provides significant speedups. It must be enabled at compile time - see installation instructions.
Reference:
Jeff Daily. Parasail: SIMD C library for global, semi-global, and local pairwise sequence alignments. (2016) BMC Bioinformatics
Parasail allows to choose from various strategies but for the sake of
simplicity, this Python binding always calls
parasail_<mode>_trace_scan_sat
which seems reasonably fast across the
global, semi-global and local modes. See parasail documentation for more
information.
You can always check if the alignment algorithms use parasail or the naive implementations by calling:
- ParasailAvailable()¶
Returns True if OpenStructure has been compiled with parasail support, False if not.
- 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 sequenceseq2 (
ConstSequenceHandle
) – A valid sequencesubst_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 sequenceseq2 (
ConstSequenceHandle
) – A valid sequencesubst_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.
- 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--hi-----jk alns = seq.alg.SemiGlobalAlign(seq_a, seq_b, seq.alg.BLOSUM62) print(alns[0].ToString(80)) # >>> A abcdefghijklmnok # >>> B --cde--hijk-----
- Parameters:
seq1 (
ConstSequenceHandle
) – A valid sequenceseq2 (
ConstSequenceHandle
) – A valid sequencesubst_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.
Substitution Weight Matrices and BLOSUM Matrices¶
- class SubstWeightMatrix¶
Substitution weights for alignment algorithms
- GetWeight(olc_one, olc_two)¶
Get
int
weight for pair of characters- Parameters:
olc_one (
string
) – first characterolc_two (
string
) – second character
- SetWeight(olc_one, olc_two, weight)¶
Set
int
weight for pair of characters- Parameters:
olc_one (
string
) – first characterolc_two (
string
) – second characterweight (
int
) – the weight
- GetMinWeight()¶
Returns the minimal weight of the matrix
- GetMaxWeight()¶
Returns the maximum weight of the matrix
- GetName()¶
Getter for name (empty string if not set)
- SetName(name)¶
Setter for name
- Parameters:
name (
str
) – Name to be set
Four already preset BLOSUM (BLOcks SUbstitution Matrix) matrices are available at different levels of sequence identity:
BLOSUM45
BLOSUM62
BLOSUM80
BLOSUM100
Nucleotide substitution matrices:
NUC44: Nucleotide substitution matrix used in blastn that can deal with IUPAC ambiguity codes. ATTENTION: has been edited to explicitely encode T/U equivalence, i.e. you can just do m.GetWeight(‘G’, ‘U’) instead of first translating ‘U’ to ‘T’.
They can be directly accessed upon importing the sequence module:
from ost import seq
mat = seq.alg.BLOSUM62
print(mat.GetWeight('A', 'A'))
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 indexj (
int
) – Second index
- SetScore(i, j, score)¶
Sets matrix(i,j) to score
- Parameters:
i (
int
) – First indexj (
int
) – Second indexscore (
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 alignmentweights (:class`ContactWeightMatrix`) – The weight matrix
apc_correction (
bool
) – Whether to use the APC correctionzpx_transformation (
bool
) – Whether to transform the scores into Z-scoressmall_number_correction (
float
) – initial values for the probabilities of having a given pair of amino acids p(a,b).
- CalculateContactProbability(cpred_res, method)¶
Calculate the probability of a predicted contact to be correct. This simply transforms the score associated with a prediction into a probability.
- Parameters:
cpred_res (
ost.seq.alg.ContactPredictionScoreResult
) – A contact prediction resultmethod (
str
) – The method which was used for contact prediction. Should be one of “MI”, “MIp”, “MIpz”, “cevoMI”, “cevo”
- 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 alignmentweights (: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 alignmentweights (: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.
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?
- Returns:
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.
- Return type:
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.- Returns:
Distance map.
- Return type:
- Raises:
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.
- Return type:
- Parameters:
d_map (
DistanceMap
) – Distance map as created withCreateDistanceMap()
.sigma (
float
) – Used for weighting of variance measure (seeDistances.GetWeightedStdDev()
)
- Raises:
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.
- Return type:
- Parameters:
d_map (
DistanceMap
) – Distance map as created withCreateDistanceMap()
.- Raises:
Exception if d_map has no entries.
- CreateMeanlDDTHA(d_map)¶
- Returns:
lDDT calculation based on CA carbons of the structures with lddt distance threshold of 15 Angstrom and distance difference thresholds of [0.5, 1.0, 2.0, 4.0]. The reported values for a certain structure are the mean per-residue lDDT values given all other structures as reference. Structures are in the same order as passed when creating d_map.
- Return type:
- Parameters:
d_map (
DistanceMap
) – Distance map as created withCreateDistanceMap()
.- Raises:
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 (seeCreateDistanceMap()
).- GetDataSize()¶
- Returns:
Number of pairwise distances.
- Return type:
int
- GetAverage()¶
- Returns:
Average of all distances.
- Return type:
float
- Raises:
Exception if there are no distances.
- GetMin()¶
- GetMax()¶
- Returns:
Minimal/maximal distance.
- Return type:
tuple
(distance (float
), index (int
))- Raises:
Exception if there are no distances.
- GetDataElement(index)¶
- Returns:
Element at given index.
- Return type:
tuple
(distance (float
), index (int
))- Parameters:
index (
int
) – Index within list of distances (must be <GetDataSize()
).- Raises:
Exception if there are no distances or index out of bounds.
- GetStdDev()¶
- Returns:
Standard deviation of all distances.
- Return type:
float
- Raises:
Exception if there are no distances.
- GetWeightedStdDev(sigma)¶
- Returns:
Standard deviation of all distances multiplied by exp(
GetAverage()
/ (-2*sigma) ).- Return type:
float
- Parameters:
sigma (
float
) – Defines weight.- Raises:
Exception if there are no distances.
- GetNormStdDev()¶
- Returns:
Standard deviation of all distances divided by
GetAverage()
.- Return type:
float
- Raises:
Exception if there are no distances.
- class DistanceMap¶
Container returned by
CreateDistanceMap()
. Essentially a symmetricGetSize()
xGetSize()
matrix containing up toGetNumStructures()
distances (list stored asDistances
). Indexing of residues starts at 0 and corresponds to the positions in the originally used alignment (seeCreateDistanceMap()
).- GetDistances(i_res1, i_res2)¶
- Returns:
List of distances for given pair of residue indices.
- Return type:
- Parameters:
i_res1 (
int
) – Index of residue.i_res2 (
int
) – Index of residue.
- GetSize()¶
- Returns:
Number of residues in map.
- Return type:
int
- GetNumStructures()¶
- Returns:
Number of structures originally used when creating the map (see
CreateDistanceMap()
).- Return type:
int
- class VarianceMap¶
Container returned by
CreateVarianceMap()
. LikeDistanceMap
, it is a symmetricGetSize()
xGetSize()
matrix containing variance measures. Indexing of residues is as inDistanceMap
.- Get(i_res1, i_res2)¶
- Returns:
Variance measure for given pair of residue indices.
- Return type:
float
- Parameters:
i_res1 (
int
) – Index of residue.i_res2 (
int
) – Index of residue.
- GetSize()¶
- Returns:
Number of residues in map.
- Return type:
int
- 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.- Raises:
Exception if the file cannot be opened for writing.
- GetJsonString()¶
- GetData()¶
Gets all the data in this map at once. Note that this is much faster (10x speedup observed) than parsing
GetJsonString()
or usingGet()
on each element.
- class Dist2Mean¶
Container returned by
CreateDist2Mean()
. Stores distances to mean forGetNumResidues()
residues ofGetNumStructures()
structures. Indexing of residues is as inDistanceMap
. Indexing of structures goes from 0 toGetNumStructures()
- 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.
- Return type:
float
- Parameters:
i_res (
int
) – Index of residue.i_str (
int
) – Index of structure.
- GetNumResidues()¶
- Returns:
Number of residues.
- Return type:
int
- GetNumStructures()¶
- Returns:
Number of structures.
- Return type:
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.- Raises:
Exception if the file cannot be opened for writing.
- GetJsonString()¶
- Returns:
A JSON formatted list of
GetNumResidues()
lists withGetNumStructures()
distances.- Return type:
str
- GetData()¶
Gets all the data in this map at once. Note that this is much faster (10x speedup observed) than parsing
GetJsonString()
or usingGet()
on each element.- Returns:
A list of
GetNumResidues()
lists withGetNumStructures()
distances.- Return type:
list
oflist
offloat
- GetSubData(num_res_to_avg)¶
Gets subset of data in this map by averaging neighboring values for num_res_to_avg residues.
- Returns:
A list of ceil(
GetNumResidues()
/num_res_to_avg) lists withGetNumStructures()
distances.- Return type:
list
oflist
offloat
- class MeanlDDT¶
Container returned by
CreateMeanlDDTHA()
. Stores mean lDDT values forGetNumResidues()
residues ofGetNumStructures()
structures. Has the exact same functionality and behaviour asDist2Mean
HMM Algorithms¶
Openstructure implements basic HMM-related functionality that aims at
calculating an HMM-HMM alignment score as described in
Soding, Bioinformatics (2005) 21(7), 951-60. This is the score which is
optimized in the Viterbi algorithm of the hhalign tool.
As a prerequisite, OpenStructure also implements adding pseudo counts to
ost.seq.ProfileHandle
in order to avoid zero probabilities for
unobserved transitions/emissions. Given these requirements, all functions
in this section require HMM related data (transition probabilities, neff values,
etc.) to be set, which is the case if you load a file in hhm format.
- HMMScore(profile_0, profile_1, aln, s_0_idx, s_1_idx, match_score_offset=-0.03, correl_score_weight=0.1, del_start_penalty_factor=0.6, del_extend_penalty_factor=0.6, ins_start_penalty_factor=0.6, ins_extend_penalty_factor=0.6)¶
Scores an HMM-HMM alignment given in aln between profile_0 and profile_1. The score is described in Soding, Bioinformatics (2005) 21(7), 951-60 and consists of three components:
sum of column alignment scores of all aligned columns, the match_score_offset is applied to each of those scores
sum of transition probability scores, the prefactor of those scores can be controlled with penalty factors (del_start_penalty_factor etc.)
correlation score which rewards conserved columns occuring in clusters, correl_score_weight controls its contribution to the total score
You have to make sure that proper pseudo counts are already assigned before calling this function. You can find a usage example in this documentation. This score is not necessarily consistent with the output generated with hhalign, i.e. you take the hhalign output alignment and directly feed it into this function with the same profiles and expect an equal score. The reason is that by default, hhalign performs a re-alignment step but the output score actually relates to the initial alignment coming from the Viterbi alignment. To get consistent results, run hhalign with the -norealign flag.
- Parameters:
profile_0 (
ost.seq.ProfileHandle
) – First profile to be scoredprofile_1 (
ost.seq.ProfileHandle
) – Second profile to be scoredaln (
ost.seq.AlignmentHandle
) – Alignment connecting the two profiless_0_idx (
int
) – Idx of sequence in aln that describes profile_0s_1_idx (
int
) – Idx of sequence in aln that describes profile_1match_score_offset (
float
) – Offset which is applied to each column alignment scorecorrel_score_weight (
float
) – Prefactor to control contribution of correlation score to total scoredel_start_penalty_factor (
float
) – Factor which is applied for each transition score starting a deletiondel_extend_penalty_factor (
float
) – Factor which is applied for each transition score extending a deletionins_start_penalty_factor (
float
) – Factor which is applied for each transition score starting an insertionins_extend_penalty_factor (
float
) – Factor which is applied for each transition score extending an insertion
- Raises:
Exception if profiles don’t have HMM information assigned or specified sequences in aln don’t match with profile SEQRES. Potentially set sequence offsets are taken into account.
Example with pseudo count assignment:
from ost import io, seq
prof_query = io.LoadSequenceProfile("query.hhm")
prof_tpl = io.LoadSequenceProfile("tpl.hhm")
aln = io.LoadAlignment("aln.fasta")
# assign pseudo counts to transition probabilities
seq.alg.AddTransitionPseudoCounts(prof_query)
seq.alg.AddTransitionPseudoCounts(prof_tpl)
# hhblits/hhalign 3 assign different pseudo counts to
# query and template. The reason is computational efficiency.
# The more expensive Angermueller et al. pseudo counts
# are assigned to the query.
path_to_crf = "/path/to/hh-suite/data/context_data.crf"
lib = seq.alg.ContextProfileDB.FromCRF(path_to_crf)
seq.alg.AddAAPseudoCounts(prof_query, lib)
# templates are assigned the computationally cheaper pseudo
# counts derived from a Gonnet substitution matrix
seq.alg.AddAAPseudoCounts(prof_tpl)
# assign null model pseudo counts
# this should be done AFTER you assigned pseudo counts to emission
# probabilities as this affects the result
seq.alg.AddNullPseudoCounts(prof_query)
seq.alg.AddNullPseudoCounts(prof_tpl)
print("score:", seq.alg.HMMScore(prof_query, prof_tpl, aln, 0, 1))
- AddNullPseudoCounts(profile)¶
Adds pseudo counts to null model in profile as implemented in hhalign. Conceptually we’re mixing the original null model with the frequencies observed in the columns of profile. The weight of the original null model depends on the neff value of profile. This function should be called AFTER you already assigned pseudo counts to the emission probabilities as this affects the result.
- Parameters:
profile (
ost.seq.ProfileHandle
) – Profile to add pseudo counts- Raises:
Exception if profile doesn’t have HMM information assigned
- AddTransitionPseudoCounts(profile, gapb=1.0, gapd=0.15, gape=1.0)¶
Adds pseudo counts to the transition probabilities in profile as implemented in hhalign with equivalent parameter naming and default parameterization. The original transition probabilities are mixed with prior probabilities that are controlled by gapd and gape. Priors:
priorM2I = priorM2D = gapd * 0.0286
priorM2M = 1.0 - priorM2D - priorM2I
priorI2I = priorD2D = 1.0 * gape / (gape - 1.0 + 1.0/0.75)
priorI2M = priorD2M = 1.0 - priorI2I
Transition probabilities of column i starting from a match state are then estimated with pM2X = (neff[i] - 1) * pM2X + gape * priorM2X. Starting from an insertion/deletion state we have pI2X = neff_ins[i] * pI2X + gape * priorI2X. In the end, all probabilities are normalized such that (pM2M, pM2I, pM2D) sum up to one, (pI2M, pI2I) sum up to one and (pD2I, pD2D) sum up to one.
- Parameters:
profile (
ost.seq.ProfileHandle
) – Profile to add pseudo counts- Raises:
Exception if profile doesn’t have HMM information assigned
- AddAAPseudoCounts(profile, a=1.0, b=1.5, c=1.0)¶
Adds pseudo counts to the emission probabilities in profile by mixing in probabilities from the Gonnet matrix as implemented in hhalign with equivalent parameter naming and default parameterization. We only implement the diversity-dependent mode for the mixing factor tau (default in hhalign), which for column i depends on neff[i] , a , b and c .
- Parameters:
profile (
ost.seq.ProfileHandle
) – Profile to add pseudo counts- Raises:
Exception if profile doesn’t have HMM information assigned
- class ContextProfileDB¶
Database that contains context profiles which will be used to add pseudo counts as described by Angermueller et al., Bioinformatics (2012) 28, 3240-3247.
- FromCRF(filename)¶
Static load function which reads a crf file provided in an hh-suite installation. Default location: “path/to/hhsuite/data/context_data.crf”
- Parameters:
filename (
str
) – Filename of CRF file
- Save(filename)¶
Saves database in OST-internal binary format which can be loaded faster than a crf file.
- Parameters:
filename (
str
) – Filename to save db
- Load(filename)¶
Static load function that loads database in OST-internal binary format.
- Parameters:
filename (
str
) – Filename of db
- AddAAPseudoCounts(profile, db, a=0.9, b=4.0, c=1.0)
Adds pseudo counts to the emission probabilities in profile by utilizing context profiles as described in Angermueller et al., Bioinformatics (2012) 28, 3240-3247. We only implement the diversity-dependent mode for the mixing factor tau (default in hhalign), which for column i depends on neff[i] , a , b and c .
- Parameters:
profile (
ost.seq.ProfileHandle
) – Profile to add pseudo countsdb (
ContextProfileDB
) – Database of context profiles
- Raises:
Exception if profile doesn’t have HMM information assigned