Alignment.align2d() — align sequences with structures

align2d(overhang=0, align_block=0, rr_file='$(LIB)/as1.sim.mat', align_what='BLOCK', off_diagonal=100, max_gap_length=999999, local_alignment=False, matrix_offset=0.0, gap_penalties_1d=(-100.0, 0.0), gap_penalties_2d=(3.5, 3.5, 3.5, 0.2, 4.0, 6.5, 2.0, 0.0, 0.0), surftyp=1, fit=True, fix_offsets=(0.0, -1.0, -2.0, -3.0, -4.0), input_weights_file=None, output_weights_file=None, n_subopt=0, subopt_offset=0.0, input_profile_file=None, output_profile_file=None, weigh_sequences=False, smooth_prof_weight=10, weights_type='SIMILAR', break_break_bonus=10000.0, io=None)

IMPORTANT NOTE: This command is obsolete, and is no longer maintained. It is strongly recommended that you use Alignment.salign() instead.

This command aligns a block of sequences (second block) with a block of structures (first block). It is the same as the Alignment.align() command except that a variable gap opening penalty is used. This gap penalty depends on the 3D structure of all sequences in block 1. The variable gap penalty can favor gaps in exposed regions, avoid gaps within secondary structure elements, favor gaps in curved parts of the mainchain, and minimize the distance between the two C_α positions spanning a gap. The Alignment.align2d() command is preferred for aligning a sequence with structure(s) in comparative modeling because it tends to place gaps in a better structural context. See Section A.1.2 for the dynamic programming algorithm that implements the variable gap penalty. gap_penalties_2d specifies parameters $\omega_H$, $\omega_S$, $\omega_B$, $\omega_C$, $\omega_D$, $d_o$, $\gamma$, $t$ and $\omega_SC$. (Section A.1.2). The default gap penalties gap_penalties_1d and gap_penalties_2d as well as the rr_file substitution matrix were found to be optimal in pairwise alignments of structures and sequences sharing from 30% to 45% sequence identity [Madhusudhan et al., 2006].

The linear gap penalty function for inserting a gap in block 1 of structures is: $g = f_1(H, S, B, C, SC) u + l v$ where $u$ and $v$ are the usual gap opening and extension penalties, $l$ is gap length, and $f_1$ is a function that is at least 1, but can be larger to make gap opening more difficult in the following circumstances: between two consecutive (i.e., $i, i+1$) helical positions, two consecutive β-strand positions, two consecutive buried positions, or two consecutive positions where the mainchain is locally straight. This function is $f_1 = 1 + [\omega_H H_i H_{i+1} + \omega_S S_i S_{i+1} + \omega_B B_i B_{i+1} + \omega_C C_i C_{i+1} + \omega_SC SC_i SC_{i+1}]$, $H_i$ is the fraction of helical residues at position $i$ in block 1, $S_i$ is the fraction of β-strand residues at position $i$ in block 1, $B_i$ is the average relative sidechain buriedness of residues at position $i$ in block 1, $C_i$ is the average straightness of residues at position $i$ in block 1, and $SC_i$ is the structural conservedness at position $i$ in block 1. See Section 6.6.32 for the definition of these features. The original straightness is modified here by assigning maximal straightness of 1 to all residues in a helix or a β-strand. The structural conservedness of the residues in block 1 are imported from an external source input_profile_file. The structural conservedness at a particular position gives the likelihood of the occurrence of a gap when structurally similar regions from all known protein structures are aligned structurally.

The linear gap penalty function for opening a gap in block 2 of sequences is: $g = f_2(H, S, B, C, D, SC) u + l v$ where $f_2$ is a function that is at least 1, but can be larger to make the gap opening in block 2 more difficult in the following circumstances: when the first gap position is aligned with a helical residue, a β-strand residue, a buried residue, extended mainchain, or when the whole gap in block 2 is spanned by two residues in block 1 that are far apart in space. This function is $f_2 = 1 + [\omega_H H_i + \omega_S S_i + \omega_B B_i + \omega_C C_i + \omega_D \sqrt{d-d_o} + \omega_SC SC_i]$. $d$ is the distance between the two C_α atoms spanning the gap, averaged over all structures in block 1 and $d_o$ is the distance that is small enough to correspond to no increase in the opening gap penalty (e.g., 8.6Å).

To find the best alignment, this method backtracks through the dynamic programming matrix, effectively adding gaps up to max_gap_length. Thus, for optimum performance you may want to reduce this parameter from its default value.

Other parameters are described in Alignment.align().

When fit is False, no alignment is done and the routine returns only the average structural information, which can be written out by the Alignment.write() command.

Example: examples/commands/align2d.py

# Demonstrating ALIGN2D, aligning with variable gap penalty

from modeller import *
log.verbose()
env = Environ()

env.libs.topology.read('$(LIB)/top_heav.lib')
env.io.atom_files_directory = ['../atom_files']

# Read aligned structure(s):
aln = Alignment(env)
aln.append(file='toxin.ali', align_codes='2ctx')
aln_block = len(aln)

# Read aligned sequence(s):
aln.append(file='toxin.ali', align_codes='2nbt')

# Structure sensitive variable gap penalty sequence-sequence alignment:
aln.align2d(overhang=0, gap_penalties_1d=(-100, 0),
            gap_penalties_2d=(3.5, 3.5, 3.5, 0.2, 4.0, 6.5, 2.0, 0., 0.),
            align_block=aln_block)

aln.write(file='align2d.ali', alignment_format='PIR')
aln.write(file='align2d.pap', alignment_format='PAP',
          alignment_features='INDICES HELIX BETA STRAIGHTNESS ' + \
                             'ACCESSIBILITY CONSERVATION')
aln.check()

# Color the first template structure according to gaps in alignment:
aln = Alignment(env)
aln.append(file='align2d.ali', align_codes=('2ctx', '2nbt'),
           alignment_format='PIR', remove_gaps=True)
mdl = Model(env)
mdl.read(file=aln['2ctx'].atom_file,
         model_segment=aln['2ctx'].range)
mdl.color(aln=aln)
mdl.write(file='2ctx.aln.pdb')

# Color the first template structure according to secondary structure:
mdl.write_data(file='2ctx', output='SSM')
mdl.write(file='2ctx.ssm.pdb')

# Superpose the target structure onto the first template:
mdl2 = Model(env)
mdl2.read(file=aln['2nbt'].atom_file,
          model_segment=aln['2nbt'].range)
sel = Selection(mdl).only_atom_types('CA')
sel.superpose(mdl2, aln)
mdl2.write(file='2nbt.fit.pdb')