# Dynamite Documentation, for version 1.5a

### Ewan Birney
Sanger Centre
Wellcome Trust Genome Campus
Hinxton, Cambridge CB10 1SA,
England.
Email: birney@sanger.ac.uk

### 18/6/2001

## Contents

## 1 Introduction

Dynamite is a code generating language whose main purpose is to produce efficient code for dynamic programming. Dynamic programming is used throughout sequence analysis in a variety of different guises (sequence alignment, gene prediction etc). Dynamite is used extensively by myself (Ewan Birney), in particular for the Wise2 package.

Dynamite itself is probably of interest to a small number of people, who are quite skilled programmers that want to make new algorithms. However the programs developed by Dynamite are potentially useful for many different people.

The production of Dynamite has necessitated two other developments

- Firstly Dynamite has an internal “object” concept it uses. Dynamite is not a full OOP language, but rather can be classified as an Object based language. This is very useful for me internally, and may be useful for other people. In particular one development in this object model is to provide glue into Perl of the Dynamite objects.
- Secondly Dynamite requires a run-time library to link to. This run time library is designed inside Dynamite’s object model, which provides a consistent way of handling objects. The run-time library does contain things which are very dynamite specific but alot of the code is generic sequence handling, codon handling etc.

### 1.1 Thanks

Dynamite would not be possible without the help of a number of people. Richard Durbin, my supervisor has encouraged me and helped in the design of Dynamite. Ian Holmes has been a constant teacher of the theory behind this and a great help in many aspects of Dynamite. He also wrote the technical section at the end of this document. Lars Averstead helped by being an early tester and added some new functionality to the compiler. Finally I am very appreciative of the people who are currently struggling with Dynamite, such as Lisa Davies, Ramu Chenna and David Kulp.

### 1.2 How to read this documentation

Dynamite has quite a steep learning curve. This is due to three things

- Dynamite was written by a biochemist picking up Computer Science. The syntax is a little weird
- Dynamite provides a very flexible approach to dynamic programming. If you are not use to this class of algorithm then you will have to pick up some concepts on the way
- Dynamite is designed for library orientated C generation. You have to program ontop of the routines generated by Dynamite to provide an actual application. For people who aren’t used to C, there is the C learning curve, which is not as shallow as, say, Perl

This means the people who are most at ease with Dynamite are those who know C and dynamic programming. Unfortunately these are the people with the least need for Dynamite! Here is some help for reading this documentation

#### 1.2.1 Good C programmer, new to dynamic programming

Read the first example, and then the “Use your own types” section, (9section.9),

#### 1.2.2 new to C, familiar with dynamic programing/sequence alignment

Make sure your C set up is correct, by compiling (and running) a “Hello World” program. Read the examples, trying to compile and run them each time.

#### 1.2.3 new to C, new to dynamic programming

The curve is particularly steep for you. You need to make sure your C set up is good (find someone who compiles C on your system to help you) and then try to read around the topic, partly using this documentation and other dynamic programming references. I learnt dynamic programming from “Dynamic programming algorithms for biological sequence comparison.”, in Methods Enzymol 1992;210:575-601. I would highly recommend “Biological sequence analysis” by Messers Durbin, Eddy, Krogh and Mitchison.

## 2 An introduction to specific readers

Dynamite lies firmly between three different disciplines. The language has been specifically designed for bioinformatics problems, and therefore Dynamite could be of interest to people with a biochemical or molecular biology background. The algorithm basis however is the well known dynamic programming used in many aspects of computer science. In this case dynamite produces algorithms which are probably most similar to certain types of speech-recognition Hidden Markov Models, and in general considering these algorithms as probablistic finite state machines provides a neat theoretical way of setting parameters. Computer scientists may well be interested in Dynamite from this aspect. Finally Dynamite is interested in providing efficient, stable and extensible code, and there are aspects of Dynamite which are involved solely in programming issues.

### 2.1 Biological readers.

Dynamite is interested in making sequence alignment and database searching algorithms. It could be directly used, for example, to make an equivalent of SSEARCH in Bill Pearson’s FASTA package. The BLAST algorithms can be considered as an approximation to such algorithm types: In the more recent BLAST2 and WU-BLAST packages, a deliberate ’gapped alignment’ stage has been added which uses dynamic programming. A good introduction to Dynamic programming was written by Pearson and Miller in Methods in enzymology

However there are a number of possible dynamic programming algorithms to represent different types of biology. If one accepts certain restrictions (in particular, each position in a sequence must be treated independantly of any distant other location in a sequence), then one can model quite complicated biological concepts inside the database search. For example, one can use the fact that alpha helices have different residue biases than beta strands inside the comparison of two proteins.

### 2.2 Computer science readers.

Dynamic programming is used in so many different areas it is hard to know how best to describe dynamite’s role. The basic sequence alignment problem can be considered a shortest string edit problem, with somewhat arbitary costs for each edit coming from emperical estimates of biological mutations.

Another rich vein of similarity is with Hidden Markov Modeling of speech. In this case a word may be represented by a series of states that emit speech waveform at certain probabilities and transitions between states have defined probabilities. The dynamic programming algorithm is an efficient form of calculating probabilities of observing a given waveform given that it came from this word. This can be done over all possible paths (summed in probability space) or be used to find the most likely path through the state network (Viterbi algorithm). In molecular biology, the waveform is protein or DNA sequence, and the word is a model of a protein structure. An excellent package for this is the HMMer package by Sean Eddy. The summed matrix is very useful for determining parameters of such a model simply from data: the viterbi form of DP is often used for determining whether a new sequence does or does not contain this protein structure.

The final theoretical approach taken by many people has been a grammar based linguist approach to biological sequence analysis. This has many desirable features, in particular the heirarchy of phoneme, word, phrase, sentance etc seems very similar to the many biological heirarchies eg, base pair, exon, transcript, gene (to take one of many examples). Although the algorithms can be written as grammars, with classical grammar production rules, the traditional approach for parsing grammars, a depth first parse tree is rarely taken in biosequence analysis. Rather dynamic programming is employed which can be considered a breadth first approach. This is principly because in biosequence analysis there is no basic deterministic transitions, and all possible parses of a grammar are usually valid. Therefore depth first or complicated breadth first parsing does not help you in restricting the paths (or parses) that a “sequence” can take. In addition, it is likely that the memory is laid out arbitarily making implementations slow for unimportant reasons.

However, parse trees often allow much greater flexibility in adding new components in particular in heirarchies. In many ways, dynamite is an attempt to provide the same sort of flexibility in being able to change the grammar without scarificing the implmentation speed.

Searls and Murphy published a similar package to Dynamite focusing on some of the more theoretical aspects of Finite State Machines. Dynamite is more focused on practical aspects of the problem.

### 2.3 Programmers

I started out as a biologist. (I am not sure what I am now). Dynamite has really suffered from having to cope with some tough programming concepts in an pretty naive mind. Dynamite acts as a C-generating machine. You could think of it as a poorly written Cfront, or a program rather like yacc or a complex function template in C++. Because Dynamite has to produce reasonably efficient code it produces probably far too much C where I have by and large trusted the C optimiser to do a good job on it.

The dynamic programming algorithm has been implemented in many different ways. In perhaps its cleanest view, it can represented as an directed acyclic graph with different weights on each arc. The ’best’ (depending on your definition of best, max or min) path needs to be found (this is the viterbi algorithm).

This view however takes away a level of granularity of the DP problem as presented in biology. In essence there is one sub graph which is used to produce the entire graph, just that each sub graph is reparameterised on the basis of its position in the full graph. The inner loop of the algorithm therefore is usually written with the sub-graph explicitly laid out in the code, allowing for very good optimisation of memory look ups and also potential for micro-parallelisation (if the sub-graph allows it). The dynamite definition is precisely this sub-graph.

At a coarser level, the algorithms are used in a “embarrasingly parallel” fashion of it simply being repeated for a series of objects for a database search. There is no link between each object in terms of which ones need to be done first. Therefore very simple parallelisation strategies can be used.

So, in terms of implementation, there are many avenues to improve efficiency. These range from specialised hardware implementations which are flexible enough to take a sub-set of biologically sensible dynamic programming algorithms to simple pthreaded code for SMP boxees

The hope is that dynamite can hide the implementation development from the algorithm development, allowing people to focus on one or the other without having to worry that they will not be able to take advantage of the best algorithms or machines respectively.

## 3 I’m here to help

I am well aware that the learning curve to Dynamite is quite steep. Dynamite is my own language and suffers from the fact that it has been developed with only really one user. At the moment I am trying very hard to make it more useful to other people. (This documentation is a start).

If you do find something confusing, please email me (birney@sanger.ac.uk) and I’ll get a reply as soon as I can, usually within a day.

## 4 Overview of Dynamite (for programmers)

Programmers may appreciate a quick top-level overview of Dynamite. You can separate dynamite into 3 areas: the dynamite compiler which produces the C code, (acting much like yacc), the dynamite support libraries which are required for this code to run, ie hold some generic matrix and alignment structures, and the dynamite run-time library which provides predefined biologically relevant types. The run-time library and the dynamite compiler can communicate certain concepts through a file called ’methods’ which is bit like a dynamite compiler headerfile. When you are starting with Dynamite is probably not worth using this system, but if you want to get into the database searching aspects of dynamite, it is a necessity. This communication allows a stronger form of typing than offered by the C implementation: in particular macros can be typed using this method. In addition the communication is on a biologically ’logical’ level so theoretically one can drop out the dynamite run-time library, replace it with one’s own specific biological typed library (the ncbi toolkit for example) and recompile with a new ’methods’ file. This has to assume that certain biological concepts are represented in a particular way (for example, Base A is represented by 0, Base C 1 etc), but does allow the underlying data structures complete freedom.

If that’s all gibberish, then wait until I get around to writing an example of swapping dynamite run-time libraries.

## 5 Installation

The installation is a pretty straight forward on UNIX machines. Pick up the distribution from ’ftp://ftp.sanger.ac.uk/pub/birney/dynamite/dyn.x.tar.Z’ (where x is the latest release). Then uncompress and untar, cd into the directory made (something like dyn1.0a) and type make ie:

    %zcat dyn1.0a.tar.Z | tar -xvf -
    %cd dyn1.0a
    %make
    %cd examples

I have not tested dynamite out at all on non UNIX systems. To make on linux systems go

    %cd dyc
    %make linux

## 6 Example1 (psw)

Enough waffle! Lets see a real life example of a piece of code. This can be found in the examples directory

For our first example we’ll take the well known case of smith waterman protein alignment. This aligns to proteins using a comparison matrix and gap penalties

### 6.1 proteinsw.dy

proteinsw.dy is the actual dynamite source file to describe smith waterman.

    %{
    #include "dyna.h"

    %}

    matrix ProteinSW
    query  type="PROTEIN"  name="query"
    target type="PROTEIN"  name="target"
    resource type="COMPMAT" name="comp"
    resource type="int"    name="gap"
    resource type="int"    name="ext"
    state MATCH offi="1" offj="1"
      calc="AAMATCH(comp,AMINOACID(query,i),AMINOACID(target,j))"
      source MATCH
        calc="0"
        endsource
      source INSERT
        calc="0"
        endsource
      source DELETE
        calc="0"
        endsource
      source START
        calc="0"
        endsource
      query_label  SEQUENCE
      target_label SEQUENCE
    endstate
    state INSERT offi="0" offj="1"
      source MATCH
        calc="gap"
        endsource
      source INSERT
        calc="ext"
        endsource
      query_label  INSERT
      target_label SEQUENCE
    endstate
    state DELETE offi="1" offj="0"
      source MATCH
        calc="gap"
        endsource
      source DELETE
        calc="ext"
        endsource
      query_label  SEQUENCE
      target_label INSERT
    endstate
    state START !special !start
      query_label  START
      target_label START
    endstate
    state END !special !end
      source MATCH
        calc="0"
        endsource
      query_label  END
      target_label END
    endstate
    endmatrix

This file provides a single Dynamite defintion. Lets ignore the precise definitions and concentrate on how to use this file to make a workable application.

At the moment there is no option for the dynamite compiler to build all the necessary C code including the main function and ’pretty’ output displays. To make the actual dynamite code one needs to compile the file proteinsw.dy with dyc (the dynamite compiler). This makes a .c and a .h file of ANSI C code

    %pwd
    /nfs/disk21/birney/prog/dyn_release/dyn1.0a/examples
    %../bin/dyc proteinsw.dy
    %more proteinsw.c
    #ifdef _cplusplus
    extern "C" {
    #endif
    #include "proteinsw.h"

    # line 5 "proteinsw.c"

      /*****************   C functions  ****************/
      /*             Written using dynamite            */
      /*            Wed Jan 21 15:04:45 1998           */
      /*            email birney@sanger.ac.uk          */
      /* http://www.sanger.ac.uk/Software/Dynamite     */
      /*************************************************/

Then this file can be compiled to a .o file and linked to an appropiate “driver” main C file. In this example it is in psw.c

### 6.2 psw.c (driver file)

psw.c contains the main function which is going to actually use the function code produced from proteinsw.dy

    /*
     * include proteinsw.h - will include the dynamite
     * produced declarations provided
     */

    #include "proteinsw.h"

    /*
     * seqaligndisplay - fancy display
     *
     */
    #include "seqaligndisplay.h"

    void show_help(FILE * ofp)
    {
      fprintf(ofp,"\npsw <options> seq1 seq2\nBoth sequences in fasta format\n"
        "\tOPTIONS\n"
        "\t-g gap penalty (default 12)\n"
        "\t-e ext penatly (default 2)\n"
        "\t-m comp matrix (default blosum62.bla)\n"
        "\t-r show raw output\n"
        "\t-l show label output\n"
        "\t-f show fancy output\n"
        "\t   (default, -f, all outputs can be shown together\n"
        );
    }

    int main(int argc,char ** argv)
    {
      Sequence * query;
      Sequence * target;
      ComplexSequence * query_cs;
      ComplexSequence * target_cs;
      ComplexSequenceEvalSet  * evalfunc;
      CompMat * comp;
      char * comp_file;
      int gap = (12);
      int ext = (2);

      boolean show_raw_output = FALSE;
      boolean show_label_output = FALSE;
      boolean show_fancy_output = FALSE;
      boolean has_outputted = FALSE;

      PackAln * pal;
      AlnBlock * alb;

      /*
       * Process command line options
       * -h or -help gives us help
       * -g for gap value (an int) - rely on commandline error processing
       * -e for ext value (an int) - rely on commandline error processing
       * -m for matrix (a char)
       * -r - raw matrix output
       * -l - label output
       * -f - fancy output
       *
       *
       * Use calls to commandline.h functions
       *
       */

      if( strip_out_boolean_argument(&argc,argv,"h") == TRUE || strip_out_boolean_argument(&argc,argv,"-help") == TRUE) {
        show_help(stdout);
        exit(1);
      }

      show_raw_output = strip_out_boolean_argument(&argc,argv,"r");
      show_label_output = strip_out_boolean_argument(&argc,argv,"l");
      show_fancy_output = strip_out_boolean_argument(&argc,argv,"f");

      /** if all FALSE, set fancy to TRUE **/

      if( show_raw_output == FALSE && show_label_output == FALSE )
        show_fancy_output = TRUE;

      (void) strip_out_integer_argument(&argc,argv,"g",&gap);
      (void) strip_out_integer_argument(&argc,argv,"e",&ext);

      comp_file = strip_out_assigned_argument(&argc,argv,"m");
      if( comp_file == NULL)
        comp_file = "blosum62.bla";

      if( argc != 3 ) {
        warn("Must have two arguments for sequence 1 and sequence 2 %d",argc);
        show_help(stdout);
        exit(1);
      }

      /*
       * Read in two sequences
       */

      if( (query=read_fasta_file_Sequence(argv[1])) == NULL ) {
        fatal("Unable to read the sequence in file %s",argv[1]);
      }

      if( (target=read_fasta_file_Sequence(argv[2])) == NULL ) {
        fatal("Unable to read the sequence in file %s",argv[2]);
      }

      /*
       * Open a blosum matrix. This will be opened from WISECONFIGDIR
       * or WISEPERSONALDIR if it is not present in the current directory.
       */

      comp = read_Blast_file_CompMat(comp_file);

      if( comp == NULL ) {
        fatal("unable to read file %s",comp_file);
      }

      /*
       * Convert sequences to ComplexSequences:
       * To do this we need an protein ComplexSequenceEvalSet
       *
       */

      evalfunc = default_aminoacid_ComplexSequenceEvalSet();

      query_cs = new_ComplexSequence(query,evalfunc);
      if( query_cs == NULL ) {
        fatal("Unable to make a protein complex sequence from %s",query->name);
      }

      target_cs = new_ComplexSequence(target,evalfunc);
      if( target_cs == NULL ) {
        fatal("Unable to make a protein complex sequence from %s",target->name);
      }

      /*
       * Make an alignment. I don't care about the implementation:
       * If the sequences are small enough then it should use explicit memory.
       * Long sequences should use divide and conquor methods.
       *
       * Calling PackAln_bestmemory_ProteinSW is the answer
       * This function decides on the best method considering the
       * memory and changes accordingly. It frees the matrix memory
       * at the end as well.
       *
       */

      pal = PackAln_bestmemory_ProteinSW(query_cs,target_cs,comp,-gap,-ext,NULL);

      if( pal == NULL ) {
        fatal("Unable to make an alignment from %s and %s",query->name,target->name);
      }

      /*
       * ok, make other alignment forms, and be ready to show
       */

      alb = convert_PackAln_to_AlnBlock_ProteinSW(pal);

      /*
       * show output. If multiple outputs, divide using //
       */

      if( show_raw_output == TRUE ) {
        show_simple_PackAln(pal,stdout);
        puts("//\n");
      }

      if( show_label_output == TRUE ) {
        show_flat_AlnBlock(alb,stdout);
      }

      if( show_fancy_output == TRUE ) {
        write_pretty_seq_align(alb,query,target,15,50,stdout);
        puts("//\n");
      }

      /*
       * Destroy the memory.
       */

      free_Sequence(query);
      free_Sequence(target);
      free_CompMat(comp);
      free_ComplexSequence(query_cs);
      free_ComplexSequence(target_cs);
      free_PackAln(pal);
      free_AlnBlock(alb);

      return 0;
    }

This seems to have alot of code to do a simple thing. However, if you notice, most of the code in interested in I/O. The core algorithm is in this section

      /*
       * Convert sequences to ComplexSequences:
       * To do this we need an protein ComplexSequenceEvalSet
       *
       */

      evalfunc = default_aminoacid_ComplexSequenceEvalSet();

      query_cs = new_ComplexSequence(query,evalfunc);
      if( query_cs == NULL ) {
        fatal("Unable to make a protein complex sequence from %s",query->name);
      }

      target_cs = new_ComplexSequence(target,evalfunc);
      if( target_cs == NULL ) {
        fatal("Unable to make a protein complex sequence from %s",target->name);
      }

      /*
       * Make an alignment. I don't care about the implementation:
       * If the sequences are small enough then it should use explicit memory.
       * Long sequences should use divide and conquor methods.
       *
       * Calling PackAln_bestmemory_ProteinSW is the answer
       * This function decides on the best method considering the
       * memory and changes accordingly. It frees the matrix memory
       * at the end as well.
       *
       */

      pal = PackAln_bestmemory_ProteinSW(query_cs,target_cs,comp,-gap,-ext,NULL);

      if( pal == NULL ) {
        fatal("Unable to make an alignment from %s and %s",query->name,target->name);
      }

There are two things going on here. The call

      pal = PackAln_bestmemory_ProteinSW(query_cs,target_cs,comp,-gap,-ext);

is there the actual call to the dynamite produced algorithm. It is only one function - dynamite has done all the decision making about what implementation of the algorithm to use, and also the algorithm itself and free’d all memory usage.

Before that there was the rather counter intuitive production of ’ComplexSequence’ objects in lines like

      query_cs = new_ComplexSequence(query,evalfunc);
      if( query_cs == NULL ) {
        fatal("Unable to make a protein complex sequence from %s",query->name);
      }

Why do this? Well the answer is that Dynamite really needs a slightly more extended concept of a ’sequence’ to allow to cope with things like splice sites or repeats (to give two obvious example) which somehow exsist in DNA sequence without being easily derivable from the sequence.

(A deeper answer is that the dynamite configuration setup which has been shipped with the dynamite compiler has provided this need for ComplexSequences. In fact this configuration can change, allowing dynamite to produce code for nearly any type of data layout that can be got at via C syntax. But that’s beyond the scope of this level of explanation)

The only thing at the moment you need to realise is that

- Each ’type’ of biological sequence is going to represented by a specific ComplexSequenceEvalSet. This ComplexSequenceEvalSet maps a sequence to the correct data layout for the algorithm
- The dynamite run-time libary contains default ComplexSequenceEvalSet’s for protein, cdna and genomic types

### 6.3 Running the example

Lets make and run the example

    %make psw
    %psw
    Warning Error
            Must have two arguments for sequence 1 and sequence 2 1

    psw <options> seq1 seq2
    Both sequences in fasta format
            OPTIONS
            -g gap penalty (default 12)
            -e ext penatly (default 2)
            -m comp matrix (default blosum62.bla)
            -r show raw output
            -l show label output
            -f show fancy output
               (default, -f, all outputs can be shown together
    %psw roac.pep roah.pep

    Q22037               EPENLRKIFVGGLTSNTTDDLMREFYSQFGEITDIIVMRDPTTKRSRGF
                         EPE LRK+F+GGL+  TTD+ +R  + Q+G +TD +VMRDP TKRSRGF
    ROA1_HUMAN           EPEQLRKLFIGGLSFETTDESLRSHFEQWGTLTDCVVMRDPNTKRSRGF

    Q22037               GFVTFSGKTEVDAAMKQRPHIIDGKTVDPKRAVPRDDKNRSESNVSTKR
                         GFVT++   EVDAAM  RPH +DG+ V+PKRAV R+D  R  ++++ K+
    ROA1_HUMAN           GFVTYATVEEVDAAMNARPHKVDGRVVEPKRAVSREDSQRPGAHLTVKK

    Q22037               LYVSGVREDHTEDMLTEYFTKYGTVTKSEIILDKATQKPRGFGFVTFDD
                         ++V G++ED  E  L +YF +YG +   EI+ D+ + K RGF FVTFDD
    ROA1_HUMAN           IFVGGIKEDTEEHHLRDYFEQYGKIEVIEIMTDRGSGKKRGFAFVTFDD

    Q22037               HDSVDQCVLQKSHMVNGHRCDVRKGLSKDEMSKAQMNRDRETRGGRSRD
                         HDSVD+ V+QK H VNGH C+VRK LSK EM+ A       ++ GRS
    ROA1_HUMAN           HDSVDKIVIQKYHTVNGHNCEVRKALSKQEMASAS-----SSQRGRSGS

    Q22037               GQRGGYNGGG-GGGGGWGGPAQRGGPGAYGGP-GGGGQGGYGGDYGG--
                         G  GG  GGG GG   +G      G G +GG  GGGG GG G  Y G
    ROA1_HUMAN           GNFGGGRGGGFGGNDNFGRGGNFSGRGGFGGSRGGGGYGGSGDGYNGFG

    Q22037               GWGQQGGGGQGGWGGPQQQQGGG-GWGQQGGGGQGGWGGPQQQQQGGWG
                           G  GGGG G  GG +    GG G+G QG G  GG G       GG
    ROA1_HUMAN           NDGGYGGGGPGYSGGSRGYGSGGQGYGNQGSG-YGGSGSYDSYNNGGGR

    Q22037               GPQQGGGGGGWGGQGQQQGGWGGQSGAQQWAHAQGGN
                         G   GG G  +GG G         + +  +   +GGN
    ROA1_HUMAN           G-FGGGSGSNFGGGGSYNDFGNYNNQSSNFGPMKGGN

    //

This gives back what you expect: the smith waterman alignment of the two proteins

## 7 Explaining proteinsw.dy

Lets go back to proteinsw.dy to find out how this file represents the smith waterman algorithm.

To calculate the algorithm one needs 5 pieecs of data

- Two sequences
- A comparison matrix of 20x20 amino acids for the match of each amino acid
- Gap open penalty
- Gap extension penalty

To calculate the alignment one would take a mathematical recursion something like

      (Pseudo-code)

      N = length of sequence 1
      M = length of sequence 2
      best         = -infinity
      Match[N][M]  = -infinity
      Insert[N][M] = -infinity
      Delete[N][M] = -infinity

      for( i goes 0 to N-1)
        for( j goes 0 to M-1)
          Match(i,j) = max {
              0,
              Match(i-1,j-1),
              Insert(i-1,j-1),
              Delete(i-1,j-1)
              } + MatchScore(seq[i],seq[j])
          Insert(i,j) = max {
              Match(i-1,j) - gap_open,
              Insert(i-1,j) - gap_ext
              }
          Delete(i,j) = max {
              Match(i,j-1) - gap_open,
              Delete(i,j-1) - gap_ext
              }
          best = max(best,Match(i,j))

      return best

The alignment is the set of (i,j,<state>) triples, where state is one of (Match,Insert,Delete) which provided the best score.

Sometimes this algorithm is written in a different way with a concept of a ’cell’. There is absolutely no difference of this form compared to the upper form.

      for( i goes 0 to N-1)
        for( j goes 0 to M-1)
          Cell(i,j).Match = max {
              0,
              Match(i-1,j-1),
              Insert(i-1,j-1),
              Delete(i-1,j-1)
              } + MatchScore(seq[i],seq[j])
          Cell(i,j).Insert= max {
              Match(i-1,j) - gap_open,
              Insert(i-1,j) - gap_ext
              }
          Cell(i,j).Delete= max {
              Match(i,j-1) - gap_open,
              Delete(i,j-1) - gap_ext
              }
          best = max(best,Cell(i,j).Match)

Now lets look at the Dynamite definition file:

The dynamite definition is started by a idenitifer for this code and finished with an endmatrix line

    matrix ProteinSW
    endmatrix

The 5 pieces of data are passed in at the top of the file

    query  type="PROTEIN"  name="query"
    target type="PROTEIN"  name="target"
    resource type="COMPMAT" name="comp"
    resource type="int"    name="gap"
    resource type="int"    name="ext"

The two sequences are given special ’resource’ types called query and target, meaning this are the two datatypes which are laid along the axis of the matrix. This is not so important for the calculation of a single matrix, but is important for the database searching modes.

The rest of the dynamite definition is interested in the actual recursion defintions. You can see the similarity of the recursion and the state defintions: Each ’cell unit’ or matrix (like Match, Insert or Delete) are provided with a state <name> ... endstate set of lines. Inside each state lines are the ’source’ lines which indicate the calculations done to define each number. Most times these calculations are done inside the max line of the different possible ways of making the numbers, but some of these calculations (for example, MatchScore(seq1[i],seq2[j]) ) can occur outside of the max system.

    state MATCH offi="1" offj="1"
      calc="AAMATCH(comp,AMINOACID(query,i),AMINOACID(target,j))"
      source MATCH
        calc="0"
        endsource
      source INSERT
        calc="0"
        endsource
      source DELETE
        calc="0"
        endsource
      source START
        calc="0"
        endsource
      query_label  SEQUENCE
      target_label SEQUENCE
    endstate
    state INSERT offi="0" offj="1"
      source MATCH
        calc="gap"
        endsource
      source INSERT
        calc="ext"
        endsource
      query_label  INSERT
      target_label SEQUENCE
    endstate
    state DELETE offi="1" offj="0"
      source MATCH
        calc="gap"
        endsource
      source DELETE
        calc="ext"
        endsource
      query_label  SEQUENCE
      target_label INSERT
    endstate

As this is the core of the dynamite defintion, lets review some of the needs of state and source lines

- Each state name must be unique
- Each source name must reference a state
- Each source must have a valid offi (offset in i) and offj (offset in j) line. These lines can be defined as defaults for the entire states.
- Each source must have a calc line, even if that is calc="0"
- Each source must have a query and target label. These can be defined as state defaults. (We’ll see more about labels /labels later/)
- Each state can have an optional calc line which is notionally ’added’ to each source’s calc line. (the implementation does make efficient use of this, and adds it after the max calculation if you are interested).

Because of the use of the state defaults issue, this means I could have written the MATCH state as follows:

    state MATCH
      source MATCH offi="1" offj="1"
              calc="AAMATCH(comp,AMINOACID(query,i),AMINOACID(target,j))"
            query_label  SEQUENCE
        target_label SEQUENCE
        endsource
      source INSERT offi="1" offj="1"
              calc="AAMATCH(comp,AMINOACID(query,i),AMINOACID(target,j))"
        query_label  SEQUENCE
        target_label SEQUENCE
        endsource
      source DELETE offi="1" offj="1"
              calc="AAMATCH(comp,AMINOACID(query,i),AMINOACID(target,j))"
        query_label  SEQUENCE
        target_label SEQUENCE
        endsource
      source START offi="1" offj="1"
              calc="AAMATCH(comp,AMINOACID(query,i),AMINOACID(target,j))"
        query_label  SEQUENCE
        target_label SEQUENCE
        endsource
    endstate

But as you can see this is a much longer definition, and if you want to change something then you have to edit quite a few lines. In addition the dynamite compiler is not yet clever enough to figure out that the calc lines are identical and therefore can be optimised. (your C optimiser may however manage this!).

The final point of the dynamite defintion is to define start and end points. In smith waterman you ’start at any point’ and ’end at any point’ (actually not true. You can start in any match, and end in any match in one interpretation).

This is represented in Dynamite as two states, both are *special* states, which are used in a broader context later, and one has the identifier !start, and one !end. I tend to actually call the states START and END, but this is not required

    state START !special !start
      query_label  START
      target_label START
    endstate
    state END !special !end
      source MATCH
        calc="0"
        endsource
      query_label  END
      target_label END
    endstate

The end state here is updated for all MATCH states, and provides the best score over the alignment.

Just to review the start/end initiation of the matrix, these are the definitions which are requires

- Each matrix must have one start state, called anything, with !start and !special defined
- Each matrix must have one end state, called anything, with !end and !special defined

The start state is the only state which is initialised to zero. This means that all alignments must start in the start state.

The end of alignment is considered to be the end state with the highest score.

As you might see, dynamite can code for many other algorithms than just smith waterman. Lets see another example, est2gen

## 8 Example2 (est2gen)

est2gen compares an est sequence to genomic sequence. In this case, assuming that the est sequence came from the same gene we have two different processes occuring:

- There are sequencing errors in the est sequence (? and maybe the genomic sequence)
- There are introns in the genomic sequence which are not present in the est sequence

The intron is going to look like a ’long gap’ in the est sequence. We are going to have ’short gaps’ as well in the est sequence compared to the genomic sequence because of sequencing error. In addition introns are going to be started with GT (perhaps a better description would be good) and finish with AG.

The solution is to write a dynamite definition file which has two types of gaps in the genomic direction, indicated by the different state, GENOMIC_INSERT and GENOMIC_INTRON. Notice how the source lines to and from GENOMIC_INSERT and GENOMIC_INTRON are different

    %{
    #include "dyna.h"

    #define DnaMatrix_Score(dnamat,base1,base2) (dnamat->score[base1][base2])

    %}

    type DNAMAT
    real    DnaMatrix*
    endtype

    method DNA_MAT_SCORE
    map    DnaMatrix_Score
    arg    DNAMAT
    arg    base
    arg    base
    return int
    endmethod

    matrix   cDNA2Gen
    query    type="CDNA" name="query"
    target   type="GENOMIC" name="target"
    resource type="DNAMAT" name="dm"
    resource type="Score" name="cdna_open"
    resource type="Score" name="cdna_ext"
    resource type="Score" name="gen_open"
    resource type="Score" name="gen_ext"
    resource type="Score" name="intron_open"
    state MATCH offi="1" offj="1"
      calc="DNA_MAT_SCORE(dm,CDNA_BASE(query,i),GENOMIC_BASE(target,j))"
      source MATCH
        calc="0"
        endsource
      source CDNA_INSERT
        calc="0"
        endsource
      source START
        calc="0"
        endsource
      source GENOMIC_INSERT
        calc="0"
        endsource
      source GENOMIC_INTRON
        calc="GENOMIC_3SS(target,j-1)"
        target_label 3SS
        endsource
      query_label SEQUENCE
      target_label SEQUENCE
    endstate
    state CDNA_INSERT offi="1" offj="0"
      source MATCH
        calc="cdna_open"
        endsource
      source CDNA_INSERT
        calc="cdna_ext"
        endsource
      query_label  SEQUENCE
      target_label INSERT
    endstate
    state GENOMIC_INSERT offi="0" offj="1"
      source MATCH
        calc="gen_open"
        endsource
      source GENOMIC_INSERT
        calc="gen_ext"
        endsource
      query_label INSERT
      target_label SEQUENCE
    endstate
    state GENOMIC_INTRON offi="0" offj="1"
      source MATCH offj="1"
        calc="GENOMIC_5SS(target,j) + intron_open"
        target_label 5SS
        endsource
      source GENOMIC_INTRON
        calc="0"
        endsource
      query_label INSERT
      target_label CENTRAL_INTRON
    endstate
    state START !special !start
    endstate
    state END !special !end
      source MATCH
        calc="0"
        endsource
      query_label END
      target_label END
    endstate
    collapse INSERT CENTRAL_INTRON
    endmatrix

Now, this file shows two different features of dynamite:

- Using different states (GENOMIC_INTRON and GENOMIC_INSERT) to represent different gapping types
- Using different types in dynamite

Let’s just look at the driver C file first, to see how we will run the program

    /*
     * include cdna2genomic.h - will include the dynamite
     * produced declarations provided
     */

    #include "cdna2genomic.h"

    /*
     * fancy display
     */
    #include "estgendisplay.h"

    void show_help(FILE * ofp)
    {
      fprintf(ofp,"\nest2gen <options> est-seq genomic-seq\nBoth sequences in fasta format\n"
        "\tOPTIONS\n"
        "\t-g gap penalty (default 2)\n"
        "\t-e ext penatly (default 1)\n"
        "\t-m match score (default 4)\n"
        "\t-n mismatch score (default -3)\n"
        "\t-r show raw output\n"
        "\t-l show label output\n"
        "\t-f show fancy output\n"
        "\t   (default, -f, all outputs can be shown together\n)"
        );
    }

    int main(int argc,char ** argv)
    {
      Sequence * query;
      Sequence * target;
      ComplexSequence * query_cs;
      ComplexSequence * target_cs;

      int gap = (2);
      int ext = (1);
      int match = (4);
      int mismatch = (-3);

      Score cdna_open,cdna_ext,gen_open,gen_ext,intron_open;

      DnaMatrix * dm;
      ComplexSequenceEvalSet * cses_g;
      ComplexSequenceEvalSet * cses_c;

      boolean show_raw_output = FALSE;
      boolean show_label_output = FALSE;
      boolean show_fancy_output = FALSE;
      boolean has_outputted = FALSE;

      PackAln * pal;
      AlnBlock * alb;

      /*
       * Process command line options
       *
       * Use calls to commandline.h functions
       *
       */

      if( strip_out_boolean_argument(&argc,argv,"h") == TRUE || strip_out_boolean_argument(&argc,argv,"-help") == TRUE) {
        show_help(stdout);
        exit(1);
      }

      show_raw_output = strip_out_boolean_argument(&argc,argv,"r");
      show_label_output = strip_out_boolean_argument(&argc,argv,"l");
      show_fancy_output = strip_out_boolean_argument(&argc,argv,"f");

      /** if all FALSE, set fancy to TRUE **/

      if( show_raw_output == FALSE && show_label_output == FALSE )
        show_fancy_output = TRUE;

      (void) strip_out_integer_argument(&argc,argv,"g",&gap);
      (void) strip_out_integer_argument(&argc,argv,"e",&ext);
      (void) strip_out_integer_argument(&argc,argv,"m",&match);
      (void) strip_out_integer_argument(&argc,argv,"n",&mismatch);

      if( argc != 3 ) {
        warn("Must have two arguments for sequence 1 and sequence 2 %d",argc);
        show_help(stdout);
        exit(1);
      }

      /*
       * Read in two sequences
       */

      if( (query=read_fasta_file_Sequence(argv[1])) == NULL ) {
        fatal("Unable to read the sequence in file %s",argv[1]);
      }

      if( (target=read_fasta_file_Sequence(argv[2])) == NULL ) {
        fatal("Unable to read the sequence in file %s",argv[2]);
      }

      /*
       * build dna matrix
       */

      dm = identity_DnaMatrix(match,mismatch);

      /*
       * Convert sequences to ComplexSequences:
       * To do this we need a cdna ComplexSequenceEvalSet and a genomic one
       *
       * Really our genomic model should be alot more complex. The 'default' one
       * has 0 at GT----AG for 5' and 3' splice sites, and NEGI (= -infinity)
       * elsewhere.
       *
       * We could build up something much better, using complexconsensi and
       * other machinery, but not for now <grin>. See the genewise code if you
       * want to get scared.
       */

      cses_g = default_genomic_ComplexSequenceEvalSet();
      cses_c = default_cDNA_ComplexSequenceEvalSet();

      query_cs = new_ComplexSequence(query,cses_c);
      if( query_cs == NULL ) {
        fatal("Unable to make a protein complex sequence from %s",query->name);
      }

      target_cs = new_ComplexSequence(target,cses_g);
      if( target_cs == NULL ) {
        fatal("Unable to make a protein complex sequence from %s",target->name);
      }

      free_ComplexSequenceEvalSet(cses_g);
      free_ComplexSequenceEvalSet(cses_c);

      /*
       * Make an alignment using best memory
       *
       * cDNA2Gen has alot more parameter space than the parameters to this
       * program. Firstly we are treating errors similarly on each side of the
       * sequences (? correct).
       *
       * Secondly there is a rather complex interaction between the gap/extension
       * of what is thought to be sequencing error and the introns. Here we have
       * one more parameter, and intron open penalty, which can be set, to prevent
       * the more permissive use of introns to 'cheat' gaps.
       *
       * One good way to parameterise all this would be to have a probabilistic
       * model of the processes, derive probabilities and then map them to ints
       * (probability.h has got these mappings, such as Probability2Score).
       *
       * but this is not done here...
       */

      pal = PackAln_bestmemory_cDNA2Gen(query_cs,target_cs,dm,gap,-ext,-gap,-ext,0,NULL);

      if( pal == NULL ) {
        fatal("Unable to make an alignment from %s and %s",query->name,target->name);
      }

      /*
       * ok, make other alignment forms, and be ready to show
       */

      alb = convert_PackAln_to_AlnBlock_cDNA2Gen(pal);

      /*
       * show output. If multiple outputs, divide using //
       */

      if( show_raw_output == TRUE ) {
        show_simple_PackAln(pal,stdout);
        puts("//\n");
      }

      if( show_label_output == TRUE ) {
        show_flat_AlnBlock(alb,stdout);
        puts("//\n");
      }

      if( show_fancy_output == TRUE ) {
        write_pretty_estgen_seq_align(alb,query,target,15,50,stdout);
        puts("//\n");
      }

      /*
       * Destroy the memory.
       */

      free_Sequence(query);
      free_Sequence(target);
      free_DnaMatrix(dm);
      free_ComplexSequence(query_cs);
      free_ComplexSequence(target_cs);
      free_PackAln(pal);
      free_AlnBlock(alb);

      return 0;
    }

Basically this looks the same as the previous example program. The line

      pal = PackAln_bestmemory_cDNA2Gen(query_cs,target_cs,dm,gap,-ext,-gap,-ext,0);

is the major algorithm call.

Notice the rather different ’ComplexSequence’ building. Here we need two different complexsequence builds, one for the genomic dna and one for the cdna.

      /*
       * Convert sequences to ComplexSequences:
       * To do this we need a cdna ComplexSequenceEvalSet and a genomic one
       *
       * Really our genomic model should be alot more complex. The 'default' one
       * has 0 at GT----AG for 5' and 3' splice sites, and NEGI (= -infinity)
       * elsewhere.
       *
       * We could build up something much better, using complexconsensi and
       * other machinery, but not for now <grin>. See the genewise code if you
       * want to get scared.
       */

      cses_g = default_genomic_ComplexSequenceEvalSet();
      cses_c = default_cDNA_ComplexSequenceEvalSet();

      query_cs = new_ComplexSequence(query,cses_c);
      if( query_cs == NULL ) {
        fatal("Unable to make a protein complex sequence from %s",query->name);
      }

      target_cs = new_ComplexSequence(target,cses_g);
      if( target_cs == NULL ) {
        fatal("Unable to make a protein complex sequence from %s",target->name);
      }

      free_ComplexSequenceEvalSet(cses_g);
      free_ComplexSequenceEvalSet(cses_c);

This example can be made and tested by

    %make est2gen
    %est2gen hn.est hn.gen

    HSHNCPA1        1    GGTGGCATTAAAGAAGACACTGAAGAACATCACCTAAGAGATGGT-TAT
                         GGTGGCATTAAAGAAGACACTGAAGAACATCACCTAAGAGAT  T T T
    HSHNRNPA        2132 GGTGGCATTAAAGAAGACACTGAAGAACATCACCTAAGAGAT--TAT-T

    HSHNCPA1        49   TTGAACAGTATGGAAAAATTGAAGTGATTGAAATCATGACTGACCGAGG
                         TTGAACAGTATGGAAAAATTGAAGTGATTGAAATCATGACTGACCGAGG
    HSHNRNPA        2178 TTGAACAGTATGGAAAAATTGAAGTGATTGAAATCATGACTGACCGAGG

    HSHNCPA1        98   CAGTGGCAAGAAAAGGGGCTTTGCCT-TAGTAACCTTTGACGACCATGA
                         CAGTGGCAAGAAAAGGGGCTTTGCCT T GTAACCTTTGACGACCATGA
    HSHNRNPA        2227 CAGTGGCAAGAAAAGGGGCTTTGCCTTT-GTAACCTTTGACGACCATGA

    HSHNCPA1        146  CTCCGTGGATAAGATTGTCA             TTCAGAAATACCATAC
                         CTCCGTGGATAAGATTGTCA             TTCAGAAATACCATAC
    HSHNRNPA        2275 CTCCGTGGATAAGATTGTCA<-2295:2387->TTCAGAAATACCATAC

    HSHNCPA1        182  TGTGAATGGCCACAACTGTGAAGTTAGAAAAGCCCTGTCAAAGCAAGAG
                         TGTGAATGGCCACAACTGTGAAGTTAGAAAAGCCCTGTCAAAGCAAGAG
    HSHNRNPA        2404 TGTGAATGGCCACAACTGTGAAGTTAGAAAAGCCCTGTCAAAGCAAGAG

    HSHNCPA1        231  ATGGCTAGTGCTTCATCCAGCCAAAGAG             GTCGAAGT
                         ATGGCTAGTGCTTCATCCAGCCAAAGAG             GTCGAAGT
    HSHNRNPA        2453 ATGGCTAGTGCTTCATCCAGCCAAAGAG<-2481:2566->GTCGAAGT

    HSHNCPA1        267  GGTTCTGGAAACTTTGGTGGTGGTCGTGGAGGTGGTTTCGGTGGGAATG
                         GGTTCTGGAAACTTTGGTGGTGGTCGTGGAGGTGGTTTCGGTGGGAATG
    HSHNRNPA        2575 GGTTCTGGAAACTTTGGTGGTGGTCGTGGAGGTGGTTTCGGTGGGAATG

    HSHNCPA1        316  ACAACTTCGGTCGTGGAGGAAACTTCAGTGGTCGTG
                         ACAACTTCGGTCGTGGAGGAAACTTCAGTGGTCGTG
    HSHNRNPA        2624 ACAACTTCGGTCGTGGAGGAAACTTCAGTGGTCGTG<-2660:2793->

    HSHNCPA1        352  GTNG-CTTTGGTGGCAGCCGTGGTGGTGGTGGATATGGTGGCAGTGGGG
                         GT G CTTTGGTGGCAGCCGTGGTGGTGGTGGATATGGTGGCAGTGGGG
    HSHNRNPA        2794 GT-GGCTTTGGTGGCAGCCGTGGTGGTGGTGGATATGGTGGCAGTGGGG

    HSHNCPA1        400  ATGGCTATAATGGATTTGGCAATGATG             GAAGCAATT
                         ATGGCTATAATGGATTTGGCAATGATG             GAAGCAATT
    HSHNRNPA        2842 ATGGCTATAATGGATTTGGCAATGATG<-2869:3805->GAAGCAATT

    HSHNCPA1        436  TTGGAGGTGGTGGAAGCTACAATGATTTTGGG--tttatgcA-CAATCA
                         TTGGAGGTGGTGGAAGCTACAATGATTTTGGG  ++ +  +A CAATCA
    HSHNRNPA        3815 TTGGAGGTGGTGGAAGCTACAATGATTTTGGGAATT-A--CAACAATCA

    HSHNCPA1        482  GTCTTCAAATTTTGGACCCATGAAGGGAGGAAATTTTGGAGGCAGAAGC
                         GTCTTCAAATTTTGGACCCATGAAGGGAGGAAATTTTGGAGGCAGAAGC
    HSHNRNPA        3861 GTCTTCAAATTTTGGACCCATGAAGGGAGGAAATTTTGGAGGCAGAAGC

    HSHNCPA1        531  TCTGGCCCCTATGGCGGTGGAGGCCAATACTTTGCAAAACCACGAAACC
                         TCTGGCCCCTATGGCGGTGGAGGCCAATACTTTGCAAAACCACGAAACC
    HSHNRNPA        3910 TCTGGCCCCTATGGCGGTGGAGGCCAATACTTTGCAAAACCACGAAACC

    HSHNCPA1        580  AAG             GTGGCTATGGCGGTTCCAGCAGCAGCAGTAGCT
                         AAG             GTGGCTATGGCGGTTCCAGCAGCAGCAGTAGCT
    HSHNRNPA        3959 AAG<-3962:4251->GTGGCTATGGCGGTTCCAGCAGCAGCAGTAGCT

    HSHNCPA1        616  ATGGCAGTGGCAGAAGATTT
                         ATGGCAGTGGCAGAAGATTT
    HSHNRNPA        4285 ATGGCAGTGGCAGAAGATTT

    //

## 9 Example 3 - using your own Sequence objects

If you are a seasoned programmer, you might want to use Dynamite’s dynamic programming engine, in particular the divide-and-conquor (linear memory) code. In this case, you have your own Sequence structures and your own comparison matrix code: you just want to use Dynamites dynamic programming engine.

First off, Dynamite is very flexible, which was one of the design goals of the program. For the basic class of dynamic programming (compare two iterated objects using a regular grammar) you can write most types in Dynamite. I have seen Dynamite be used for aspects as diverse as biological sequences, musical copyright databases and financial data.

For the first example, let’s recode smith waterman using our own sequence objects. The objects will be of the following type

    struct MySequence {
       char * seq;
       int  length_of_sequence;
    }

       and

    struct MyComparisonMatrix {
       int comp[26][26];
    }

      and the following macro

    /* macro gets out the sequence at i as 0-26 number */
    #define SEQ_POS(obj,i) (toupper(obj->seq[i])-'A')

How do we use these objects in Dynamite? Quite simply in fact. Dynamite can take raw C types and the calc lines are parsed as a strict subset of C. This means we can provide any calculation function in the middle of the dynamic programming loop. Here is smith waterman written with these objects

    %{
    #include "dyna.h"

    %}

    matrix MyProteinSW
    query  type="MySequence*"  name="query" field:len="length_of_sequence"
    target type="MySequence*"  name="target" field:len="length_of_sequence"
    resource type="MyComparisonMatrix" name="mat"
    resource type="int"    name="gap"
    resource type="int"    name="ext"
    state MATCH offi="1" offj="1"
      calc="mat.comp[SEQ_POS(query,i)][SEQ_POS(target,j)]"
      source MATCH
        calc="0"
        endsource
      source INSERT
        calc="0"
        endsource
      source DELETE
        calc="0"
        endsource
      source START
        calc="0"
        endsource
      query_label  SEQUENCE
      target_label SEQUENCE
    endstate
    state INSERT offi="0" offj="1"
      source MATCH
        calc="gap"
        endsource
      source INSERT
        calc="ext"
        endsource
      query_label  INSERT
      target_label SEQUENCE
    endstate
    state DELETE offi="1" offj="0"
      source MATCH
        calc="gap"
        endsource
      source DELETE
        calc="ext"
        endsource
      query_label  SEQUENCE
      target_label INSERT
    endstate
    state START !special !start
      query_label  START
      target_label START
    endstate
    state END !special !end
      source MATCH
        calc="0"
        endsource
      query_label  END
      target_label END
    endstate
    endmatrix

When the current dynamite compiles this code, it will generate alot of warnings that things are ’out of scope’ or untypeable. Ignore these and/or pass the -nocwarn flag to dyc to cut some of them out. I need to configure the dynamite better with respect to warnings.

Notice how the calc line

      calc="mat.comp[SEQ_POS(query,i)][SEQ_POS(target,j)]"

is precisely what you would write in normal C to get the score of the ith amino acid from the query sequence matched to the jth amino acid from the target sequence. The dynamite compiler expands this expression to be scoped properely (so that mat, query and target are held in the ’matrix’ scope and can be refered between different parts of the dynamite generated code sensibly).

Notice we passed in the Sequence objects as pointers and the Matrix object as a struct. In theory you should be able to pass in objects any way you wish, except I made a mistake and hard coded the concept of the length being from a pointer construct. In other words the line

    query  type="MySequence*"  name="query" field:len="length_of_sequence"

means dynamite generates the following code

      /* to get out the length of the query dimension */
      out->leni = query->length_of_sequence;

In other words, dynamite can only cope with query and target objects being pointers to structs with a field being their length. I consider this a bug which I should remove, with the length being a more general calc like function (ie, it could be any piece of valid C in the correct scope). Apologies.

To use this piece of Dynamite you would only have call

     PackAln * pal;
     PackAlnUnit * unit;

     /* do your own IO for your objects */

     pal = PackAln_bestmemory_MyProteinSW(seqone,seqtwo,mat,-12,-2,NULL);

     /* to investigate PackAln, loop over it */
     for(i=0;i<pal->len;i++) {
       unit = pal->pau[i];
       if( unit->state == 0 ) {
          fprintf(stdout,"%d,%d are matched with score %d\n",unit->i+1,unit->j+1,unit->score);
       }
      }

The i,j numbers in PackAln are the C type numbers, starting at 0, however they are best thought of as occuring between the letters of the actual sequence, starting at -1 for each sequence (so the number 0 sits between the first and second resiudes of the sequence). For single amino acid sequences this is usually what you expect, but for multi residue matches (eg, codons), you need to be a little careful of what is going on.

The state number starts at 0 as the first state of the Dynamite file, and increasing down the file. When I get around to writing the documentation about the AlnBlock structure sensibly I will explain how to use that sensibly.

Lets look at extending the dynamite model for raising the gap penalty in secondary structure regions of the protein. Say our structure was now

    struct MySequence {
       char * seq;
       int  length_of_sequence;
       char * ss; /* string of H,E or C */
    }

    #define IS_ALPHA_SEQ(obj,i) (obj->ss[i] == 'H' ? 1 : 0)

Now we could write the Dynamite model like this

    %{
    #include "dyna.h"

    %}

    matrix MyProteinSW2
    query  type="MySequence*"  name="query" field:len="length_of_sequence"
    target type="MySequence*"  name="target" field:len="length_of_sequence"
    resource type="MyComparisonMatrix" name="mat"
    resource type="int"    name="gap"
    resource type="int"    name="gap_alpha"
    resource type="int"    name="ext"
    state MATCH offi="1" offj="1"
      calc="mat.comp[SEQ_POS(query,i)][SEQ_POS(target,j)]"
      source MATCH
        calc="0"
        endsource
      source INSERT
        calc="0"
        endsource
      source DELETE
        calc="0"
        endsource
      source START
        calc="0"
        endsource
      query_label  SEQUENCE
      target_label SEQUENCE
    endstate
    state INSERT offi="0" offj="1"
      source MATCH
        calc="(IS_ALPHA_SEQ(target,j) == 1 ? gap_alpha : gap)"
        endsource
      source INSERT
        calc="ext"
        endsource
      query_label  INSERT
      target_label SEQUENCE
    endstate
    state DELETE offi="1" offj="0"
      source MATCH
        calc="(IS_ALPHA_SEQ(query,i) == 1 ? gap_alpha : gap)"
        endsource
      source DELETE
        calc="ext"
        endsource
      query_label  SEQUENCE
      target_label INSERT
    endstate
    state START !special !start
      query_label  START
      target_label START
    endstate
    state END !special !end
      source MATCH
        calc="0"
        endsource
      query_label  END
      target_label END
    endstate
    endmatrix

Notice that we are editing the calc lines for MATCH to INSERT and MATCH to DELETE. These edits are linked to the correct sequence (INSERT extends in the j direction, and so we look up on the j line).

This is cute but it actually one can do much better, by having a separate model of evolution, including comparison matrix and gap penalties in the model. In this case we will partion the sequences only into Alpha regions and non alpha regions

    %{
    #include "dyna.h"

    #define ALPHA_SEQ_SCORE(obj,pos) (obj->ss[pos] == 'H' ? 0 : NEGI)
    %}

    matrix MyProteinAlphaSW
    query  type="MySequence*"  name="query" field:len="length_of_sequence"
    target type="MySequence*"  name="target" field:len="length_of_sequence"
    resource type="MyComparisonMatrix" name="mat_alpha"
    resource type="int"    name="gap_alpha"
    resource type="int"    name="ext_alpha"
    resource type="MyComparisonMatrix" name="mat"
    resource type="int"    name="gap"
    resource type="int"    name="ext"
    state MATCH offi="1" offj="1"
      calc="mat.comp[SEQ_POS(query,i)][SEQ_POS(target,j)]"
      source MATCH
        calc="0"
        endsource
      source MATCH_ALPHA
        calc="0"
        endsource
      source INSERT
        calc="0"
        endsource
      source DELETE
        calc="0"
        endsource
      source START
        calc="0"
        endsource
      query_label  SEQUENCE
      target_label SEQUENCE
    endstate
    state INSERT offi="0" offj="1"
      source MATCH
        calc="gap"
        endsource
      source INSERT
        calc="ext"
        endsource
      query_label  INSERT
      target_label SEQUENCE
    endstate
    state DELETE offi="1" offj="0"
      source MATCH
        calc="ext"
        endsource
      source DELETE
        calc="ext"
        endsource
      query_label  SEQUENCE
      target_label INSERT
    endstate
    state MATCH_ALPHA offi="1" offj="1"
      calc="ALPHA_SEQ_SCORE(query,i) + ALPHA_SEQ_SCORE(target,j) + mat_alpha.comp[SEQ_POS(query,i)][SEQ_POS(target,j)]"
      source MATCH
        calc="0"
        endsource
      source MATCH_ALPHA
        calc="0"
        endsource
      source INSERT_ALPHA
        calc="0"
        endsource
      source DELETE_ALPHA
        calc="0"
        endsource
      source START
        calc="0"
        endsource
      query_label  SEQUENCE
      target_label SEQUENCE
    endstate
    state INSERT_ALPHA offi="0" offj="1"
      calc="ALPHA_SEQ_SCORE(target,j)"
      source MATCH_ALPHA
        calc="gap_alpha"
        endsource
      source INSERT_ALPHA
        calc="ext_alpha"
        endsource
      query_label  INSERT
      target_label SEQUENCE
    endstate
    state DELETE_ALPHA offi="1" offj="0"
      calc="ALPHA_SEQ_SCORE(query,i)"
      source MATCH_ALPHA
        calc="ext_alpha"
        endsource
      source DELETE_ALPHA
        calc="ext_alpha"
        endsource
      query_label  SEQUENCE
      target_label INSERT
    endstate
    state START !special !start
      query_label  START
      target_label START
    endstate
    state END !special !end
      source MATCH
        calc="0"
        endsource
      query_label  END
      target_label END
    endstate
    endmatrix

This model has a number of consequences

- Both the query and target sequences must be predicted alpha for the model to enter the _ALPHA states, and no non alpha regions are tolerated in the _ALPHA states
- The model is only allowed to switch from alpha to non alpha states and visa versa at matched positions. If a switch from alpha to non alpha or visa versa occurs “naturally” at a gap position, the algorithm will be forced to make a match

It is possible to get around both these problems, but that is left as an exercise for the reader.

## 10 Extending the Dynamite types

The previous example showed how to use your own C types in the Dynamite generated code. These types can be connected up to the Dynamite typing system. This gives a number of benefits

- The algorithms can be written without knowing too much about the libraries they are linking to provide the actual C types
- You can type-safe macros. The macros provide efficient decoding of your data structures but don’t provide good type saftey. By using dynamite types, the type safety can be enforced by the Dynamite compiler

Because of the way I wrote this documentation, this section has been written to refer to the est2gen example. The est2gen example showed how the Dynamite types can be extended to accommodate your own C structures.

The region in the dynamite blueprint area does the extension

    type DNAMAT
    real    DnaMatrix*
    endtype

    method DNA_MAT_SCORE
    map    DnaMatrix_Score
    arg    DNAMAT
    arg    base
    arg    base
    return int
    endmethod

This snippet does two things

- Defines DNAMAT as a dynamite type
- Defines DNA_MAT_SCORE as ’method’.

The actual C-implementations are defined as being

- DnaMatrix* as the type (from the real definition). DnaMatrix can be found in the dnamatrix module in the library section.
- The macro DnaMatrix_Score which is defined at the top of the cdna2genomic.dy file

This illustrates two nice features in using the Dynamite types and methods

- You can use (in the C-implementation) a macro but still have it typechecked
- The types are far extended from the C-implementation. For example, ’bases’, ’codons’ and ’amino-acids’ are all implemented as integers with appropiate ranges. Most C-compilers have to allow conversions between typedefs. The dynamite compiler however treats each type independently.

This means that if we had inadvertently written:

      calc="DNA_MAT_SCORE(dm,CDNA_BASE(query,i),GENOMIC_CODON(target,j))"

The dynamite compiler produces the following warning:

    %dyc cdna2genomic.dy
    Warning Error
    ->In preparing matrix cDNA2Gen
    -->In parsing calc line for state [MATCH] (source ind.)
            Mis-type in argument 3 of DNA\_MAT\_SCORE: wanted [base] got [codon]
    Warning Error
    ->In preparing matrix cDNA2Gen
            Failed to parse calc lines

    The following parser errors were considered fatal:
    Mistyped arguments

In addition to this type-checking, Dynamite can check some of the semantics in the calc-line usage. In particular the index i should only refer to the query dimension, whereas the index j should only refer to the target dimension.

Thus the code

      calc="DNA_MAT_SCORE(dm,CDNA_BASE(query,j),GENOMIC_BASE(target,j))"

Generates this warning:

    %dyc cdna2genomic.dy
    Warning Error
    ->In preparing matrix cDNA2Gen
    -->In parsing calc line for state [MATCH] (source ind.)
            For function CDNA\_BASE, you have arguments j and query, which do not
    expect to paired directly in a function. This is just a warning that
    you can ignore

## 11 Replacing Dynamite types

The Dynamite typing system really comes into its own when you want to use someone else’s libraries to provide the basic ’sequence’ objects to use.

All the dynamite types go through the same typing system as outlined above. The dynamite compiler picks up the file methods from either the current directory or $WISEPERSONALDIR or $WISECONFIGDIR. These provide the methods associated with the dynamite run-time library which is distributed with the dynamite compiler.

Say in proteinsw.dy you wanted to use your own implementation of the scoring matrix, which was a more generic system, able to cope with completely arbitary alphabets. The C code you had was:

    typedef struct {
      int ** score;    /* a len-by-len array */
      char * alphabet; /* the alphabet of the matrix */
      int len;         /* the length of the matrix */
    } CompMatrix;

    int score_CompMatrix(CompMatrix * mat,char one,char two)
    {
       int a,b;
       char * pos;

       if( (pos = strchr(mat->alphabet,one)) == NULL ) {
      throw_warning("Letter %c does not exist in comparison matrix",one);
            return 0;
       }
       a = pos - mat->alphabet;

       if( (pos = strchr(mat->alphabet,two)) == NULL ) {
      throw_warning("Letter %c does not exist in comparison matrix",two);
            return 0;
       }
       b = pos - mat->alphabet;

       return mat->score[a][b];
    }

How can we adapt this so that proteinsw.dy could use it. Well, the answer is a very simple change in the methods. We need to write one piece of C - like

    /* this could be a macro */

    int map_aa_number_CompMatrix(CompMatrix * mat,int a,int b)
    {
       return score_CompMatrix(mat,(char)('A'+a),(char)('A'+b));
    }

and then the methods file rather than looking like:

    type COMPMAT
    real CompMat*
    endtype

    method AAMATCH
    map CompMat_AAMATCH
    arg COMPMAT
    arg aa
    arg aa
    return int
    endmethod

would look like

    type COMPMAT
    real CompMatrix*
    endtype

    method AAMATCH
    map map_aa_number_CompMatrix
    arg COMPMAT
    arg aa
    arg aa
    return int
    endmethod

## 12 Database Searching code

Database searching code is probably one of the first things you would like to do with the Dynamite generated code. Dynamite builds code for a number of different database implementations - at the time of writing, single threaded and multi threaded implementations are provided for.

Dynamite will build a single database searching function which will compare (potentially) a database of query objects to a database of target objects. In each case, to be able to build the code, Dynamite must have a specification of the type as a logical type with additional attributes indicating how to initialise, reload and close a database associated with a type. For most uses this means sticking to provided types in the methods file (ie things like PROTEIN or GENOMIC_DNA). To provide a pthreads port yet more information is required. However, it is possible to write your own database objects and provide the necessary information in the methods file to indicate how dynamite should loop over the database. For the moment you should contact Ewan about writing your own database objects.

For many cases you will be using a database-to-database search. The more common, single object vs a database of objects is provided by ways of making a database object which simply has a single entry. These generally are cache’d sensible so it is efficient (whereas a normal database object will free the memory of each entry object once it has been used).

The implementation of the database search is provided in the DBSeachImpl object. This object can be constructed from the command line, allowing the driving program to be somewhat immune to additional implementations being provided by the dynamite compiler. Each implementation is only guarenteed to provide a single score for each pair of comparisons. For most cases as well it will provide some sort of on-the-fly indexing of the database and ways to retrieve the sequences. At the moment, you need to manually write the alignment code if you want it.

Below a simple example of a database search of protein smith waterman, comparing one sequence vs a database a protein sequences. It is dbsearch.c in the examples directory

    #include "proteinsw.h"

    void show_help (void)
    {
      fprintf(stdout,"dbsearch [options] <protein-seq> <protein-fasta-database>\n");
      fprintf(stdout,"Valid options are\n");
      /** add more options here sometime, eg comp matrix and gap penalty*/

      /** print out dbsearch options. We don't know here what implementations are
          either possible or how they are specified. Of course, there is the problem
          that we could clash our options with the dbsearchimpl options, but that
          is not too likely, and this makes this program future proof wrt to new
          implementations
      */

      show_help_DBSearchImpl(stdout);
    }

    int main ( int argc, char ** argv)
    {
      Hscore * out;
      DBSearchImpl * dbsi;
      Protein * temp;
      Sequence * query;
      ProteinDB * querydb;
      ProteinDB * prodb;
      CompMat * mat;

      ComplexSequence * query_cs;
      ComplexSequence * target_cs;
      ComplexSequenceEvalSet  * evalfunc;

      PackAln * pal;
      AlnBlock * alb;

      int i;

      /*
       * processes the command line, removing options
       * that it wants to in order to make the new DBSearchImpl
       *
       * The great thing about this is that this programs does not
       * care about which implementation is used, and does not know either (!)
       *
       */

      dbsi = new_DBSearchImpl_from_argv(&argc,argv);

      if( argc != 3 ) {
        show_help();
        exit(1);
      }

      /*
       * first argument is a single sequence. Read it in and make it
       * into a database
       */

      query = read_fasta_file_Sequence(argv[1]);
      if( query == NULL )
        fatal("Cannot read sequence in %s\n",argv[1]);

      querydb = new_ProteinDB_from_single_seq(query);

      /*
       * Second argument is a real database. This call is
       * a nice short cut for doing this.
       */

      prodb = single_fasta_ProteinDB(argv[2]);

      if( prodb == NULL )
        fatal("Cannot read protein database in %s\n",argv[2]);

      /*
       * This is where all the results are stored. It also
       * on-the-fly stores distribution information ready
       * to be fitted by a extreme value distribution
       */

      /* 10 means a score cutoff of 10, -1 means don't report on stderr search progress */

      out = std_score_Hscore(10,-1);

      mat = read_Blast_file_CompMat("blosum62.bla");

      if( search_ProteinSW(dbsi,out,querydb,prodb,mat,-12,-2) != SEARCH_OK )
        fatal("Some sort of error in the database search. Dieing ungracefully");

      sort_Hscore_by_score(out);

      evalfunc = default_aminoacid_ComplexSequenceEvalSet();

      query_cs = new_ComplexSequence(query,evalfunc);
      if( query_cs == NULL ) {
        fatal("Unable to make a protein complex sequence from %s",query->name);
      }

      for(i=0;i<10 && i < out->len;i++) {
          fprintf(stdout,"Comparison to %s was %d score\n",out->ds[i]->target->name,out->ds[i]->score);

          /*
           * Retrieve the protein from the database
           */

          temp = get_Protein_from_ProteinDB(prodb,out->ds[i]->target);

          /*
           * Make a complex sequence of it - see psw for info on this
           */

          target_cs = new_ComplexSequence(temp->baseseq,evalfunc);
          if( target_cs == NULL ) {
     fatal("Unable to make a protein complex sequence from %s",temp->baseseq->name);
          }

          /*
           * Actually align it
           */

          pal = PackAln_bestmemory_ProteinSW(query_cs,target_cs,mat,-12,-2,NULL);

          if( pal == NULL ) {
     fatal("Unable to make an alignment from %s and %s",query->name,temp->baseseq->name);
          }

          alb = convert_PackAln_to_AlnBlock_ProteinSW(pal);

          write_pretty_seq_align(alb,query,temp->baseseq,15,50,stdout);
          puts("//\n");

          free_Protein(temp);
          free_ComplexSequence(target_cs);
      }

      return 0;
    }

The main database searching call is at

      if( search_ProteinSW(dbsi,out,querydb,prodb,mat,-12,-2) != SEARCH_OK )
        fatal("Some sort of error in the database search. Dieing ungracefully");

It takes the two database objects (querydb and prodb) instead of the normal query and target sequences, and the resources required for this comparison.

The dbsi variable is an instance of the DBSearchImpl object made from the command line with the function call

      dbsi = new_DBSearchImpl_from_argv(&argc,argv);

The out variable is an instance of the Hscore object. You can find out more about this object by reading the hscore.h file in the dynlibsrc section, with additional functionality in the dynlibcross.dy section. It is made in this case by

      out = std_score_Hscore(10,-1);

This indicates that it should store every score greater than 10 (not that sensible) and not report at all about progress during the search.

Later on in the program it goes back to pull the sequences out of the database, using the call

          temp = get_Protein_from_ProteinDB(prodb,out->ds[i]->target);

This call is not an inherent part of Dynamite, but is provided for convience in the proteindb database. The protein object is then used to actually align the sequences.

### 12.1 Compiling with pthreads

As it stands the dynamite compiler is not compiling with pthread support. To do so you need to pass the -pthreads switch to the dyc compiler.

     ../dyc/dyc -pthreads proteinsw.dy

This generates the pthreading code in the corresponding .c file, but guards it with preprocessor defines so it need not be compiled if so wished. To compile the C program you need to provide the PTHREAD symbol in the preprocessor.

     cc -c -O -I../dynlibsrc/ -I../base -DPTHREAD proteinsw.c

Finally in the link phase you must also specify the threads library as in

       cc -o dbsearch dbsearch.o proteinsw.o seqaligndisplay.o
       -L../dynlibsrc/ -L../base/  -ldyna -lwisebase -lm -lpthread

For convience in the examples directory these steps are put in the makefile for target dbthread. so make dbthread should produce a pthreaded piece of code.

To actually run the exectuable with threads you need to pass in the -pthread switch to dbsearch

## 13 Dynamite language definition

### 13.1 Overview

Dynamite files consist of three sections. The first two of these are mandatory, the last is optional. The first section is a yacc-like prefix block enclosed by %{ and %} delimiters. This block is copied directly to the output and will typically contain #include statements and other header information. The second section contains of a series of blueprint definitions written in the Dynamite language. These include type and access information for objects such as query and target sequences, as well as the main recursion information for the eventual generated DP code. This section is interpreted by the Dynamite parser. The third (optional) section is also delimited by %{ and %} and contains "marked-up" C. This is currently undocumented.

### 13.2 Elements of the dynamite language

The following language elements are defined within the second, blueprints section.

#### 13.2.1 Comments

Lines starting with # are treated as comments.

#### 13.2.2 Single-line statements

Single-line statements are identified by a single keyword at the start of the line whose scope extends to the end of that line. Following this keyword will typically be additional information pertaining to the keyword.

The following single-line keywords are reserved:

- query
- target
- resource
- extern
- collapse
- calcfunc
- query_label
- target_label
- arg
- real
- dbtype
- init
- reload
- close
- addentry
- name
- map
- return
- func
- des

Here is an example of a single-line statement:

    resource type="COMPMAT" name="comp"

Many of these single-line keywords are only valid within certain scopes. The scoping rules for individual keywords are described within the relevant sections below.

#### 13.2.3 Multi-line statements

Multi-line statement blocks are delimited by pairs of keywords of the form tag and endtag. Within the block delimited by these keywords will typically be additional information and constructs pertaining to that keyword. The following multi-line keywords are reserved:

- matrix...endmatrix
- state...endstate
- source...endsource
- method...endmethod
- type...endtype
- struct...endstruct
- api...endapi
- object...endobject

Here is an example of a multi-line statement:

    source MATCH
      calc="0"
    endsource

Many of these multi-line delimiter keywords are only valid within certain scopes. The scoping rules for individual keywords are described within the relevant sections below.

#### 13.2.4 Name-value pairs

Name-value pairs can occur within both single-line and multi-line statement blocks and are used to specify parameters in a flexible way. They take the form of a name, followed by an "=" sign, followed by a quotation-mark-enclosed value. There must be no whitespace between the name and the "=" sign, or between the "=" sign and the value.

The following name keywords are reserved:

- offi
- offj
- calc

The examples of single-line and multi-line statement blocks in the preceding sections also contain name-value pairs. Here is another example:

    offi="3"

### 13.3 Definition of the file layout

We now move to a more detailed description of the Dynamite file layout.

#### 13.3.1 Section 1 - C header

    %{
    #include "dyna.h"
    #include "my_file.h"

    %}

These lines are exported into the .h file verbatim. All files using the standard dynamite run-time libraries should use #include "dyna.h"

The %{ and %} delimiters must occur in the file.

If a second set of lines delimited by %{ and %} are found, these are exported into the .c file as “marked-up” C. This is only really used by myself (ewan).

#### 13.3.2 Section 2 - Dynamite blueprint definitions

Dynamite has a series of ’blueprints’ which are sort of the top level block which the dynamite compiler acts on. For example, each DP matrix definition is a single blueprint.

Following the first %{ %} block can be any number of dynamite blueprint definitions. These are multi-line statement blocks as described above. A complete list of allowed blueprints follows:

- matrix...endmatrix - describes the central dynamic programming recursion routines
- type...endtype - describes typing information for objects used by the generated code
- method...endmethod - describes data-access methods for objects used by the generated code
- struct...endstruct - describes the Dynamite object model (undocumented)
- api...endapi - describes the Dynamite API model (undocumented)

The most commonly used blueprint is "matrix...endmatrix".

### 13.4 Definition of the matrix blueprint

We now consider the matrix blueprint in more detail.

Essentially the Dynamite compiler needs to know two things: (1) what kinds of sequence-like objects are to be compared, and (2) what recursion relations are to comprise the main loop of the dynamic programming routines. The Dynamite code must therefore specify (1) type and access information for the sequence objects, (2) information about the transitions between the states of the alignment automaton, and (3) how to calculate the scores of these transitions. Requirements (1) and (3) are essentially "data" requirements and we consider these first, in the section entitled "Data objects". Requirement (2) is to do with the model definition and is considered in the section entitled "Model definition".

#### 13.4.1 Data objects

The following single-line statements are used to specify resources. These resources will be visible to the single-line "calc" statements that specify how to calculate the scores of transitions between states (see "Model definition", below).

- query - describes the query sequence object. Only one of these is allowed.
- target - describes the target sequence object. Only one of these is allowed.
- resource - describes resource objects such as substitution matrices and other scoring parameters.
- extern - describes resource objects with external linkage (i.e. linkage to #define properties in C code).

In all of these lines, both of the following name-value pairs are MANDATORY:

- name="name" name of the object
- type="type" type of the object, either a C-type or a logical dynamite type

#### 13.4.2 Model definition

The main recursion is defined in terms of a finite-state automaton. Each state is described by a "state...endstate" multi-line block. Within each "state...endstate" block, one or more "source...endsource" blocks describe the various possible transitions into that state. States do not need to be declared in any particular order.

Each "state...endstate" block has the following format

    state NAME
    endstate

where NAME is the name of the state.

Each "state...endstate" block must contain at least one "source...endsource" block.

A state can also be declared "special" as follows:

    state NAME !special
    endstate

DP calculations for special states occur "outside" the dynamic programming matrix, in that only one cell score for special states is maintained for each residue of the target sequence.

The following name-value pairs can be used within the "state...endstate" block (outside the scope of any enclosed "source...endsource" blocks):

- offi="<int>" - to specify the default i-offset (i.e. query sequence offset) for transitions
- offj="<int>" - to specify the default j-offset (i.e. target sequence offset) for transitions
- calc="<expr>" - to specify a transition-independent expression to be added to all transition scores

The following single-line statements can be used within the "state...endstate" block (outside the scope of any enclosed "source...endsource" blocks):

- query_label - to specify the default transition label for the query sequence
- target_label - to specify the default transition label for the target sequence

#### 13.4.3 Transition definitions

Each "source...endsource" block must have the form

    source NAME
    endsource

where the NAME is the name of the source state for this transition.

Within each "source...endsource" block, the following name-value pairs may override the defaults specified in the "state...endstate" block:

- offi="int" - offset in i for this transition
- offj="int" - offset in j for this transition

These attributes must be specified somewhere, so if no default is specified in the "state...endstate" block, then they must be specified here.

Within each "source...endsource" block, the following single-line statements may override the defaults specified in the "state...endstate" block:

- query_label LABEL - label for query
- target_label LABEL - label for target

These attributes must be specified somewhere, so if no default is specified in the "state...endstate" block, then they must be specified here.

Each "source...endsource" block muys

Each "source...endsource" block must contain a single calc definition. The format for these is described in the following section.

#### 13.4.4 Calc lines

Calc lines describe how to calculate the scores of transitions between states. They provide the major interface between Dynamite-generated code and your own C routines. The calc lines can use the following C expression syntax:

- ,,,
  (N.B. floating point calculations are not recommended in calc lines)
- (pointer deference), . (struct deference), and -> (pointer + struct dereference)
- (array access)
- standard function calls, including method calls to Dynamite objects
- variables declared in query, target, resource and extern blocks

Dynamite will automatically type-check method calls to pre-defined Dynamite objects. This type safety is stronger than that offered by the C compiler and is recommended. Where C functions unknown to Dynamite are invoked, Dynamite will issue a warning that it cannot type-check the function call.

### 13.5 Definition of the run-time set up.

The Dynamite runtime set up is found in a file called methods, picked up from .,$WISEPERSONALDIR or $WISECONFIGDIR. This file is shipped to be compatible with the dynamite run-time libraries which come with the system.

The following types are defined "logical" types

- PROTEIN
- CDNA
- GENOMIC
- COMPMAT
- aa
- codon
- base

#### 13.5.1 PROTEIN

This represents a protein sequence. valid methods are

- aa AMINOACID(PROTEIN,int) - returns amino acid at this position

#### 13.5.2 CDNA

This represents a cDNA sequence. valid methods are

- base CDNA_BASE(CDNA,int) - returns the base at this position
- codon CDNA_CODON(CDNA,int) - returns the codon at this position

#### 13.5.3 GENOMIC

This represents a Genomic sequence. valid methods are

- base GENOMIC_BASE(GENOMIC,int) - returns the base at this position
- codon GENOMIC_CODON(GENOMIC,int) - returns the codon at this position
- int GENOMIC_5SS(GENOMIC,int) - returns the 5’ splice potential for this position
- int GENOMIC_3SS(GENOMIC,int) - returns the 3’ splice potential for this position
- int GENOMIC_REPEAT(GENOMIC,int) - returns non-cds potential (ie, definitely not a coding sequence) for this position
- int GENOMIC_CDSPOT(GENOMIC,int) - returns cds potential (ie, definitely a coding sequence) for this position

#### 13.5.4 COMPMAT

This represents a comparison matrix for protein amino acids

- int AAMATCH(COMPMAT,aa,aa) returns the score for these two amino acids

#### 13.5.5 aa

No functions associated with this. It is a number from 0-25 inclusive, where ’A’ (alanine) is 0 etc... This means we have more numbers than amino acids, but that is fine.

#### 13.5.6 codon

No functions associated with this. It is a number from 0-125 inclusive, where this is calculated as (base1*25)+(base2*5)+base3 (base defined below) of a codon, at 125 means no possible codon at this position.

#### 13.5.7 base

No functions associated with this. It is a number from 0-4 inclusive where A = 0, G = 1, C = 2, T = 3 and N = 4

## 14 Errors Reported by the Dynamite Compiler

A number of different errors can be reported by the Dynamite compiler. Some of them are warnings but the dynamite compiler still generates valid code whereas others cause the dynamite compiler to stop. The errors can be divided into three types

- Syntactic errors, due to writing the wrong commands in the dynamite file
- Semantic errors in the Dynamite definition
- Syntactic and semantic errors in the calc lines of the dynamite file

Syntatic errors in the blueprint file will sometime cause a run-away set of parsing errors. These are terminated with a fatal call to prevent your stderr filling up with boring messages. Otherwise dynamite tries to gather as many errors as possible before exiting.

At the moment the dynamite compiler does not remove the .c file when it encounters an error (it should do!). This means that if an error is encountered you may well have to remove the .c file yourself for the makefile to remake it etc...

### 14.1 Syntax errors in the matrix blueprint parsing

Generally the syntax errors in the matrix blueprint parsing are quite informative. They come with a line indicating where they lie in the file. This line is the currently read line in the parser, and so might be not be precisely where the error is. The hardest error to understand is the run-away parsing error in which a “endsource” or similar multi line ending section is not closed.

Below details some of the more common errors found in parsing.

#### 14.1.1 [Dynamite Level] Did not understand line [ source MATCH]. Probably a run-away parsing error, so failing now

This error indicates that the parser has left the GenericMatrix parsing and returned to the top level dynamite level. It is usually preceded by the next error

#### 14.1.2 Unable to read GenericMatrix

The syntax error when the matrix blueprint reading fails. The precise error should have been posted beforehand

#### 14.1.3 Got the line [state INSERT offi="0" offj="1"], a state start line inside a state. Expect you forgot an endstate

The error message says it all. This is when you have forgotten to close a multi-line brace, like state...endstate or source...endsource

#### 14.1.4 got an end tag [endxxx] but expecting a [endyyy]

xxx might be matrix whereas yyy might be source. This means that you haven’t balanced your state...endstate or source...endsource blocks correctly. You have probably missed out a endsource or a endstate line somewhere

#### 14.1.5 You have specified a modifier [name] to XXX but it has either no ’=’ sign or no quoted argument. The ’=’ character should be flush to both the tag and the quoted (using ) argument

The dynamite parser is particularly bad about the modifer lines. They have to be of the form modifer="something" with no whitespace. This is an extremely bad thing I know.

#### 14.1.6 Could not understand line in YYY parse

YYY might be “matrix” or “state” or “source”. This means that the parser got a line it couldn’t handle. It is quite likely that you have missed out a “endxxxx” line

### 14.2 Semantic errors

Semantic errors are not reported with a line number. This is because the Dynamite compiler does not associate a line number with the parts of the datastructure which is carries around (it should do).

However semantic errors are labeled by where (logically) they come in the Dynamite file: for example the following error was made by mistyping “MATCH” in a source line, and therefore making a source line that had no state. The first error is quite descriptive about what is going on.

    adnah:[/wise2/dyc]<369>: dyc proteinsw.dy
    Warning Error
    In preparing matrix ProteinSW
            In matrix ProteinSW - State MATCH asks for source MATCH2 but there is
    no State of that name
    Warning Error
    In preparing matrix ProteinSW
            Unable to cross reference state and source
    Warning Error
    In preparing matrix ProteinSW
            Failing simple cross-checks, aborting before calc-line parsing
    Fatal Error
            A Dynamite blueprint fails semantic checks.
    Please refer to previous errors for the precise problem

#### 14.2.1 A Dynamite blueprint fails semantic checks. Please refer to previous errors for the precise problem

This error indicates that a the top level a semantic problem was found, even though the syntax was parsed ok.

#### 14.2.2 Failing simple cross-checks, aborting before calc-line parsing

This is the error message likely to be got just before the last error. This means that the semantics of the matrix blueprint failed, and the Dynamite parser did not attempt to parse the calc lines.

#### 14.2.3 Start/End points are faulty

This indicates that the start or end points where not there or not special states. There is individual error messages for each of those.

#### 14.2.4 You have not got a !start special, but you do have a START special state. Presuming that you wanted to make that the START ;)

This is warning. You haven’t used the proper dynamite syntax of !start to indicate the start state. You do however have a special state called START and there is all likelhood that this is what you want as the start. There is a similar error message for end

#### 14.2.5 Unable to prepare labels

This means that you have missed out a “query_label” or a “target_label” definition for one source line. Even if you don’t use label alignments, the labels are still required. Remember that you can place them as defaults for a particular state.

#### 14.2.6 Unable to resolve all the cell offset refs into correct offsets

This means that you have probably forgotten an offi and or an offj modifier for a source line. Remember that you can set them as defaults at the state line if you so wish. It is probably preceeded by the next error

#### 14.2.7 In state INSERT (source MATCH), both offi and offj are zero: dynamite cannot currently handle cell internal references

The most likely explanation for this is that you have not specified the offsets for this source at all. If you have specified them as 0,0 this is illegal in Dynamite. In needs at least one offset to be non zero.

### 14.3 Syntax and Semantics of calc line parsing

Once all the semantics of the Dynamite blueprint is ok, the compiler then turns its attention to the calc lines. This is like a mini-parser operating inside the dynamite parser, but unlike the dynamite parser, this once was written in yacc/lex and is a more vanilla parser environment.

#### 14.3.1 Parser Syntax error on calc line

This is the error which indicates that the yacc parser cannot get through the calc line. Above it should be precisely which calc line gives the error and where on it the parser halted. For example

    adnah:[/wise2/dyc]<377>: dyc proteinsw.dy
    Warning Error
    In preparing matrix ProteinSW
    In parsing calc line for state [MATCH] source [START]
            Calc line parser error: [syntax error]
    Occured at:
    0 += 2
    ---^
    Warning Error
    In preparing matrix ProteinSW
            Failed to parse calc lines

    The following parser errors were considered fatal:
    Parser Syntax error on calc line

    Fatal Error
            A Dynamite blueprint fails semantic checks.
    Please refer to previous errors for the precise problem

This will cause a fatal error which you will need to solve before proceeding. A number of other errors can be generated by the parsing of calc lines, but most of them do not cause a fatal error at the end of the day.

> This document was translated from LATEX by H E V E A.