\documentstyle{article}
\begin{document}
\newcommand{\programtext}[1]{{\tt #1}}
\title{Dynamite Documentation, for version 1.5a}
\author{Ewan Birney\\
Sanger Centre\\
Wellcome Trust Genome Campus\\
Hinxton, Cambridge CB10 1SA,\\
England.\\
Email: birney@sanger.ac.uk}
\maketitle
\newpage
\tableofcontents
\newpage
\section{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
\begin{itemize}
\item 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.
\item 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.
\end{itemize}
\subsection{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.
\subsection{How to read this documentation}
Dynamite has quite a steep learning curve. This is due to three things
\begin{itemize}
\item Dynamite was written by a biochemist picking up Computer Science.
The syntax is a little weird
\item 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
\item 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
\end{itemize}
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
\subsubsection{Good C programmer, new to dynamic programming}
Read the first example, and then the ``Use your own types'' section, (\ref{own_objects}),
\subsubsection{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.
\subsubsection{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.
\section{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.
\subsection{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.
\subsection{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.
\subsection{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.
\section{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.
\section{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 assumme 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.
\section{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:
\begin{verbatim}
%zcat dyn1.0a.tar.Z | tar -xvf -
%cd dyn1.0a
%make
%cd examples
\end{verbatim}
I have not tested dynamite out at all on non UNIX systems. To make on linux
systems go
\begin{verbatim}
%cd dyc
%make linux
\end{verbatim}
\section{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
\subsection{proteinsw.dy}
proteinsw.dy is the actual dynamite source file to describe
smith waterman.
\begin{verbatim}
%{
#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
\end{verbatim}
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
\begin{verbatim}
%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 */
/*************************************************/
\end{verbatim}
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
\subsection{psw.c (driver file)}
psw.c contains the main function which is going to actually
use the function code produced from proteinsw.dy
\begin{verbatim}
/*
* 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;
}
\end{verbatim}
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
\begin{verbatim}
/*
* 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);
}
\end{verbatim}
There are two things going on here. The call
\begin{verbatim}
pal = PackAln_bestmemory_ProteinSW(query_cs,target_cs,comp,-gap,-ext);
\end{verbatim}
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
\begin{verbatim}
query_cs = new_ComplexSequence(query,evalfunc);
if( query_cs == NULL ) {
fatal("Unable to make a protein complex sequence from %s",query->name);
}
\end{verbatim}
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
\begin{itemize}
\item 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
\item The dynamite run-time libary contains default ComplexSequenceEvalSet's
for protein, cdna and genomic types
\end{itemize}
\subsection{Running the example}
Lets make and run the example
\begin{verbatim}
%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
//
\end{verbatim}
This gives back what you expect: the smith waterman alignment of the two proteins
\section{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
\begin{itemize}
\item Two sequences
\item A comparison matrix of 20x20 amino acids for the match of each amino acid
\item Gap open penalty
\item Gap extension penalty
\end{itemize}
To calculate the alignment one would take a mathematical recursion
something like
\begin{verbatim}
(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
\end{verbatim}
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.
\begin{verbatim}
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)
\end{verbatim}
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
\begin{verbatim}
matrix ProteinSW
endmatrix
\end{verbatim}
The 5 pieces of data are passed in at the top of the file
\begin{verbatim}
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"
\end{verbatim}
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.
\begin{verbatim}
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
\end{verbatim}
As this is the core of the dynamite defintion, lets review some
of the needs of state and source lines
\begin{itemize}
\item Each state name must be unique
\item Each source name must reference a state
\item 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.
\item Each source must have a calc line, even if that is calc="0"
\item Each source must have a query and target label. These can be defined as state
defaults. (We'll see more about labels /{labels} later/)
\item 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).
\end{itemize}
Because of the use of the state defaults issue, this means I could
have written the MATCH state as follows:
\begin{verbatim}
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
\end{verbatim}
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
\begin{verbatim}
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
\end{verbatim}
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
\begin{itemize}
\item Each matrix must have one start state, called anything, with
!start and !special defined
\item Each matrix must have one end state, called anything, with !end
and !special defined
\end{itemize}
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
\section{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:
\begin{itemize}
\item There are sequencing errors in the est sequence (? and maybe the genomic sequence)
\item There are introns in the genomic sequence which are not present in the est sequence
\end{itemize}
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
\begin{verbatim}
%{
#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
\end{verbatim}
Now, this file shows two different features of dynamite:
\begin{itemize}
\item Using different states (GENOMIC\_INTRON and GENOMIC\_INSERT) to represent
different gapping types
\item Using different types in dynamite
\end{itemize}
Let's just look at the driver C file first, to see how we will run the
program
\begin{verbatim}
/*
* 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;
}
\end{verbatim}
Basically this looks the same as the previous example program. The line
\begin{verbatim}
pal = PackAln_bestmemory_cDNA2Gen(query_cs,target_cs,dm,gap,-ext,-gap,-ext,0);
\end{verbatim}
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.
\begin{verbatim}
/*
* 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);
\end{verbatim}
This example can be made and tested by
\begin{verbatim}
%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
//
\end{verbatim}
\section{Example 3 - using your own Sequence objects}
\label{own_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 \emph{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
\begin{verbatim}
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')
\end{verbatim}
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
\begin{verbatim}
%{
#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
\end{verbatim}
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
\begin{verbatim}
calc="mat.comp[SEQ_POS(query,i)][SEQ_POS(target,j)]"
\end{verbatim}
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
\begin{verbatim}
query type="MySequence*" name="query" field:len="length_of_sequence"
\end{verbatim}
means dynamite generates the following code
\begin{verbatim}
/* to get out the length of the query dimension */
out->leni = query->length_of_sequence;
\end{verbatim}
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
\begin{verbatim}
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);
}
}
\end{verbatim}
The i,j numbers in PackAln are the C type numbers, starting
at 0, however they are best thought of as occuring \emph{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
\begin{verbatim}
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)
\end{verbatim}
Now we could write the Dynamite model like this
\begin{verbatim}
%{
#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
\end{verbatim}
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
\begin{verbatim}
%{
#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
\end{verbatim}
This model has a number of consequences
\begin{itemize}
\item 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
\item 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
\end{itemize}
It is possible to get around both these problems, but that is left as an exercise
for the reader.
\section{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
\begin{itemize}
\item The algorithms can be written without knowing too much
about the libraries they are linking to provide the actual C types
\item 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
\end{itemize}
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
\begin{verbatim}
type DNAMAT
real DnaMatrix*
endtype
method DNA_MAT_SCORE
map DnaMatrix_Score
arg DNAMAT
arg base
arg base
return int
endmethod
\end{verbatim}
This snippet does two things
\begin{itemize}
\item Defines DNAMAT as a dynamite type
\item Defines DNA\_MAT\_SCORE as 'method'.
\end{itemize}
The actual C-implementations are defined as being
\begin{itemize}
\item DnaMatrix* as the type (from the real definition). DnaMatrix can be found in the dnamatrix module in the library section.
\item The macro DnaMatrix\_Score which is defined at the top of the cdna2genomic.dy file
\end{itemize}
This illustrates two nice features in using the Dynamite types and methods
\begin{itemize}
\item You can use (in the C-implementation) a macro but still have it typechecked
\item 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.
\end{itemize}
This means that if we had inadvertently written:
\begin{verbatim}
calc="DNA_MAT_SCORE(dm,CDNA_BASE(query,i),GENOMIC_CODON(target,j))"
\end{verbatim}
The dynamite compiler produces the following warning:
\begin{verbatim}
%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
\end{verbatim}
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
\begin{verbatim}
calc="DNA_MAT_SCORE(dm,CDNA_BASE(query,j),GENOMIC_BASE(target,j))"
\end{verbatim}
Generates this warning:
\begin{verbatim}
%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
\end{verbatim}
\section{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:
\begin{verbatim}
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];
}
\end{verbatim}
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
\begin{verbatim}
/* 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));
}
\end{verbatim}
and then the methods file rather than looking like:
\begin{verbatim}
type COMPMAT
real CompMat*
endtype
method AAMATCH
map CompMat_AAMATCH
arg COMPMAT
arg aa
arg aa
return int
endmethod
\end{verbatim}
would look like
\begin{verbatim}
type COMPMAT
real CompMatrix*
endtype
method AAMATCH
map map_aa_number_CompMatrix
arg COMPMAT
arg aa
arg aa
return int
endmethod
\end{verbatim}
\section{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
\emph{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
\begin{verbatim}
#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;
}
\end{verbatim}
The main database searching call is at
\begin{verbatim}
if( search_ProteinSW(dbsi,out,querydb,prodb,mat,-12,-2) != SEARCH_OK )
fatal("Some sort of error in the database search. Dieing ungracefully");
\end{verbatim}
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
\begin{verbatim}
dbsi = new_DBSearchImpl_from_argv(&argc,argv);
\end{verbatim}
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
\begin{verbatim}
out = std_score_Hscore(10,-1);
\end{verbatim}
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
\begin{verbatim}
temp = get_Protein_from_ProteinDB(prodb,out->ds[i]->target);
\end{verbatim}
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.
\subsection{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.
\begin{verbatim}
../dyc/dyc -pthreads proteinsw.dy
\end{verbatim}
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.
\begin{verbatim}
cc -c -O -I../dynlibsrc/ -I../base -DPTHREAD proteinsw.c
\end{verbatim}
Finally in the link phase you must also specify the threads library as in
\begin{verbatim}
cc -o dbsearch dbsearch.o proteinsw.o seqaligndisplay.o
-L../dynlibsrc/ -L../base/ -ldyna -lwisebase -lm -lpthread
\end{verbatim}
For convience in the examples directory these steps are put in the makefile for
target dbthread. so {\tt 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
\newpage
\section{Dynamite language definition}
\subsection{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 {\tt \%\{ }
and {\tt \%\}} 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 \emph{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 {\tt \%\{} and {\tt \%\}}
and contains "marked-up" C. This is currently undocumented.
\subsection{Elements of the dynamite language}
The following language elements are defined within the
second, blueprints section.
\subsubsection{Comments}
Lines starting with \# are treated as comments.
\subsubsection{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:
\begin{itemize}
\item query
\item target
\item resource
\item extern
\item collapse
\item calcfunc
\item query\_label
\item target\_label
\item arg
\item real
\item dbtype
\item init
\item reload
\item close
\item addentry
\item name
\item map
\item return
\item func
\item des
\end{itemize}
Here is an example of a single-line statement:
\begin{verbatim}
resource type="COMPMAT" name="comp"
\end{verbatim}
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.
\subsubsection{Multi-line statements}
Multi-line statement blocks are delimited by pairs of keywords
of the form {\tt tag} and {\tt 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:
\begin{itemize}
\item matrix...endmatrix
\item state...endstate
\item source...endsource
\item method...endmethod
\item type...endtype
\item struct...endstruct
\item api...endapi
\item object...endobject
\end{itemize}
Here is an example of a multi-line statement:
\begin{verbatim}
source MATCH
calc="0"
endsource
\end{verbatim}
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.
\subsubsection{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:
\begin{itemize}
\item offi
\item offj
\item calc
\end{itemize}
The examples of single-line and multi-line statement blocks in the
preceding sections also contain name-value pairs. Here is another
example:
\begin{verbatim}
offi="3"
\end{verbatim}
\subsection{Definition of the file layout}
We now move to a more detailed description of the Dynamite file layout.
\subsubsection{Section 1 - C header}
\begin{verbatim}
%{
#include "dyna.h"
#include "my_file.h"
%}
\end{verbatim}
These lines are exported into the .h file verbatim.
All files using the standard dynamite run-time libraries
should use \#include "dyna.h"
The {\tt \%\{} and {\tt \%\}} delimiters must occur in the file.
If a second set of lines delimited by {\tt \%\{} and {\tt \%\}} are found,
these are exported into the .c file as ``marked-up'' C. This is only
really used by myself (ewan).
\subsubsection{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 {\tt \%\{} {\tt \%\}} 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:
\begin{itemize}
\item matrix...endmatrix - describes the central dynamic programming recursion routines
\item type...endtype - describes typing information for objects used by the generated code
\item method...endmethod - describes data-access methods for objects used by the generated code
\item struct...endstruct - describes the Dynamite object model (undocumented)
\item api...endapi - describes the Dynamite API model (undocumented)
\end{itemize}
The most commonly used blueprint is "matrix...endmatrix".
\subsection{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".
\subsubsection{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).
\begin{itemize}
\item query - describes the query sequence object. Only one of these is allowed.
\item target - describes the target sequence object. Only one of these is allowed.
\item resource - describes resource objects such as substitution matrices and other scoring parameters.
\item extern - describes resource objects with external linkage (i.e. linkage to \#define properties in C code).
\end{itemize}
In all of these lines, both of the following name-value pairs are MANDATORY:
\begin{itemize}
\item {\tt name="name"} name of the object
\item {\tt type="type"} type of the object, either a C-type or a logical dynamite type
\end{itemize}
\subsubsection{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
\begin{verbatim}
state NAME
endstate
\end{verbatim}
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:
\begin{verbatim}
state NAME !special
endstate
\end{verbatim}
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):
\begin{itemize}
\item {\tt offi="<int>"} - to specify the default i-offset (i.e. query sequence offset) for transitions
\item {\tt offj="<int>"} - to specify the default j-offset (i.e. target sequence offset) for transitions
\item {\tt calc="<expr>"} - to specify a transition-independent expression to be added to all transition scores
\end{itemize}
The following single-line statements can be used within the "state...endstate" block
(outside the scope of any enclosed "source...endsource" blocks):
\begin{itemize}
\item query\_label - to specify the default transition label for the query sequence
\item target\_label - to specify the default transition label for the target sequence
\end{itemize}
\subsubsection{Transition definitions}
Each "source...endsource" block must have the form
\begin{verbatim}
source NAME
endsource
\end{verbatim}
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:
\begin{itemize}
\item {\tt offi="int"} - offset in i for this transition
\item {\tt offj="int"} - offset in j for this transition
\end{itemize}
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:
\begin{itemize}
\item query\_label LABEL - label for query
\item target\_label LABEL - label for target
\end{itemize}
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.
\subsubsection{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:
\begin{itemize}
\item {\tt \+,\-,\*,\\} (N.B. floating point calculations are not recommended in calc lines)
\item {\tt \*} (pointer deference), . (struct deference), and -> (pointer + struct dereference)
\item {\tt \[\]} (array access)
\item standard function calls, including method calls to Dynamite objects
\item variables declared in query, target, resource and extern blocks
\end{itemize}
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.
\subsection{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
\begin{itemize}
\item PROTEIN
\item CDNA
\item GENOMIC
\item COMPMAT
\item aa
\item codon
\item base
\end{itemize}
\subsubsection{PROTEIN}
This represents a protein sequence. valid methods are
\begin{itemize}
\item aa AMINOACID(PROTEIN,int) - returns amino acid at this position
\end{itemize}
\subsubsection{CDNA}
This represents a cDNA sequence. valid methods are
\begin{itemize}
\item base CDNA\_BASE(CDNA,int) - returns the base at this position
\item codon CDNA\_CODON(CDNA,int) - returns the codon at this position
\end{itemize}
\subsubsection{GENOMIC}
This represents a Genomic sequence. valid methods are
\begin{itemize}
\item base GENOMIC\_BASE(GENOMIC,int) - returns the base at this position
\item codon GENOMIC\_CODON(GENOMIC,int) - returns the codon at this position
\item int GENOMIC\_5SS(GENOMIC,int) - returns the 5' splice potential for this position
\item int GENOMIC\_3SS(GENOMIC,int) - returns the 3' splice potential for this position
\item int GENOMIC\_REPEAT(GENOMIC,int) - returns non-cds potential (ie, definitely not a coding sequence) for this position
\item int GENOMIC\_CDSPOT(GENOMIC,int) - returns cds potential (ie, definitely a coding sequence) for this position
\end{itemize}
\subsubsection{COMPMAT}
This represents a comparison matrix for protein amino acids
\begin{itemize}
\item int AAMATCH(COMPMAT,aa,aa) returns the score for these two amino acids
\end{itemize}
\subsubsection{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.
\subsubsection{codon}
No functions associated with this. It is a number from 0-125 inclusive, where
this is calculated as \mbox{(base1*25)+(base2*5)+base3} (base defined below) of
a codon, at 125 means no possible codon at this position.
\subsubsection{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
\section{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
\begin{itemize}
\item Syntactic errors, due to writing the wrong commands in the dynamite file
\item Semantic errors in the Dynamite definition
\item Syntactic and semantic errors in the calc lines of the dynamite file
\end{itemize}
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...
\subsection{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.
\subsubsection{[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
\subsubsection{Unable to read GenericMatrix}
The syntax error when the matrix blueprint reading fails. The precise
error should have been posted beforehand
\subsubsection{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
\subsubsection{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
\subsubsection{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
{\tt modifer="something"} with no whitespace. This is an extremely bad thing I know.
\subsubsection{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
\subsection{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.
\begin{verbatim}
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
\end{verbatim}
\subsubsection{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.
\subsubsection{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.
\subsubsection{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.
\subsubsection{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
\subsubsection{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.
\subsubsection{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
\subsubsection{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.
\subsection{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 enviroment.
\subsubsection{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
\begin{verbatim}
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
\end{verbatim}
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.
\end{document}
distrib
>
Mageia
>
7
>
armv7hl
>
media
>
core-release
>
by-pkgid
>
af1a17339be7a0066c33c138028ddbfb
>
files
>
38
