BUILD_PROFILE -- Build a profile for a given sequence or alignment

RR_FILE = <string:1>	'$(LIB)/as1.sim.mat'	input residue-residue scoring file
DIRECTORY = <string:1>	''	directory list (e.g., 'dir1:dir2:dir3:./:/')
GAP_PENALTIES_1D = <real:2>	900 50	gap creation and extension penalties for sequence/sequence alignment
MATRIX_OFFSET = <real:1>	0.00	substitution matrix offset for local alignment
STOP_ON_ERROR = <integer:1>	1	whether to stop on error
N_PROF_ITERATIONS = <integer:1>	3	number of iterations in PROFILE_SEARCH
CHECK_PROFILE = <logical:1>	on	whether to monitor profile degenration
OUTPUT_SCORES = <logical:1>	off	whether to output individual scores in a build_profile scan
OUTPUT_SCORE_FILE = <string:1>	'default'	output file for writing out individual scores in seqfilter
MAX_ALN_EVALUE = <real:1>	0.1	Max. E-value of alignments to include in BUILD_PROFILE
GAPS_IN_TARGET = <logical:1>	off	whether to include gaps in target when using build_profile

Output:

MODELLER_STATUS = <integer:1>

Description:

This command iteratively scans a database of sequences to build a profile for the input sequence or alignment. The command calculates the score for a Smith-Waterman local alignment between the input sequence and each of the sequences in the database. The significance of the alignment scores (e-values) are calculated using a procedure similar to that described by Pearson (1998).

Alignments with e-values below MAX_ALN_EVALUE are then added to the current alignment. A position-specific scoring matrix is then calculated for the current alignment and is used to search the sequence database. This procedure is repeated for N_PROF_ITERATIONS or until there are are no significant alignments below the threshold, whichever occurs first.

The initial sequence or alignment can be read in either in the profile format, with READ_PROFILE, or as an alignment using READ_ALIGNMENT. In the latter case, the alignment has to be converted to the profile format using ALN_TO_PROF.

The output contains a multiple sequence alignment (assembled) of all the homologues of the input sequence found in the database. The output can be formatted as a profile with WRITE_PROFILE or converted into any of the standard alignment formats using PROF_TO_ALN. It can then be written out to a file with WRITE_ALIGNMENT.

The fit between the observed and theoretical distributions of the z-scores is calculated after each iteration and is reported in the log file. The fit is calculated using the Kolmogorov-Smirnov D-statistic. If the CHECK_PROFILE flag is set to 'on', then the command will not proceed if the fit deviates by more than 0.04 (D-statistic).

By default, regions of the alignment that introduce gaps in the target sequence are ignored (deleted) in the final multiple alignment. But if GAPS_IN_TARGET is set to 'on', then the gaps are retained. (See below for comments).

If the OUTPUT_SCORES flag is set to 'on', then the scores of each alignment between the input sequence and each database sequence, from all iterations, will be written out to the file specified in OUTPUT_SCORE_FILE.

Comments:

1. The procedure has been optimized only for the BLOSUM62 similarity matrix.

2. The dynamic programming algorithm has been optimized for performance on Intel Itanium2 architecture. Nevertheless, the calculation is sufficiently CPU intensive. It takes about 20 min for an iteration, using an input sequence of 250aa against a database containing 500,000 sequences on an Itanium2 machine. It could take much longer on any other machine.

3. It is advisable to have GAPS_IN_TARGET set to 'off', when scanning against large databases to avoid the local-alignments inserting a huge number of gaps in the final alignments.

4. The statistics will not be accurate (or may even fail) if the database does not have sequences that represent the entire range of lengths possible.

5. The method can be used for fold-assignment by first building a profile for the target sequence by scanning against a large non-redundant sequence database (like swissprot) and then using the resulting profile to scan once against a database of sequences extracted from PDB structures. GAPS_IN_TARGET can be set to 'on' in the second step to get the complete alignments that can then be used for modeling.

Example:

SET OUTPUT_CONTROL = 1 1 1 1 1

#-- Prepare the input files

#-- Read in the sequence database
SET MINMAX_DB_SEQ_LEN = 1 40000, CLEAN_SEQUENCES = on
READ_SEQUENCE_DB SEQ_DATABASE_FILE = 'pdb95.fsa', ;
                 SEQ_DATABASE_FORMAT = 'FASTA', ;
                 CHAINS_LIST = 'all'

#-- Write the sequence database in binary form
WRITE_SEQUENCE_DB SEQ_DATABASE_FILE = 'pdb95.bin', ;
                  SEQ_DATABASE_FORMAT = 'BINARY'

#-- Now, read in the binary database
READ_SEQUENCE_DB SEQ_DATABASE_FILE = 'pdb95.bin', ;
                 SEQ_DATABASE_FORMAT = 'BINARY', ;
                 CHAINS_LIST = 'all'

#-- Read in the target sequence/alignment
READ_ALIGNMENT FILE = 'toxin.ali', ALIGNMENT_FORMAT = 'PIR'

#-- Convert the input sequence/alignment into
#   profile format
ALN_TO_PROF

#-- Scan sequence database to pick up homologous sequences
SET MATRIX_OFFSET = -450
SET RR_FILE = '${LIB}/blosum62.sim.mat'
SET GAP_PENALTIES_1D = -500 -50

BUILD_PROFILE N_PROF_ITERATIONS = 5, ;
              CHECK_PROFILE = off, ;
              MAX_ALN_EVALUE = 0.01, ;
              GAPS_IN_TARGET = off

#-- Write out the profile
WRITE_PROFILE FILE = 'buildprofile.prf'

#-- Convert the profile back to alignment format
PROF_TO_ALN

#-- Write out the alignment file
WRITE_ALIGNMENT FILE = 'buildprofile.ali', ;
                ALIGNMENT_FORMAT = 'PIR'