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:	aln (AlignmentHandle) – Alignment 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.

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.