## File: dynamite.tex

package info (click to toggle)
wise 2.4.1-21
 12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700170117021703170417051706170717081709171017111712171317141715171617171718171917201721172217231724172517261727172817291730173117321733173417351736173717381739174017411742174317441745174617471748174917501751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820182118221823182418251826182718281829183018311832183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185718581859186018611862186318641865186618671868186918701871187218731874187518761877187818791880188118821883188418851886188718881889189018911892189318941895189618971898189919001901190219031904190519061907190819091910191119121913191419151916191719181919192019211922192319241925192619271928192919301931193219331934193519361937193819391940194119421943194419451946194719481949195019511952195319541955195619571958195919601961196219631964196519661967196819691970197119721973197419751976197719781979198019811982198319841985198619871988198919901991199219931994199519961997199819992000200120022003200420052006200720082009201020112012201320142015201620172018201920202021202220232024202520262027202820292030203120322033203420352036203720382039204020412042204320442045204620472048204920502051205220532054205520562057205820592060206120622063206420652066206720682069207020712072207320742075207620772078207920802081208220832084208520862087208820892090209120922093209420952096209720982099210021012102210321042105210621072108210921102111211221132114211521162117211821192120212121222123212421252126212721282129213021312132213321342135213621372138213921402141214221432144214521462147214821492150215121522153215421552156215721582159216021612162216321642165216621672168216921702171217221732174217521762177217821792180218121822183218421852186218721882189219021912192219321942195219621972198219922002201220222032204220522062207220822092210221122122213221422152216221722182219222022212222222322242225222622272228222922302231223222332234223522362237223822392240224122422243224422452246224722482249225022512252225322542255225622572258225922602261226222632264226522662267226822692270227122722273227422752276227722782279228022812282228322842285228622872288228922902291229222932294229522962297229822992300230123022303230423052306230723082309231023112312231323142315231623172318231923202321232223232324232523262327232823292330233123322333233423352336233723382339234023412342234323442345234623472348234923502351235223532354235523562357235823592360236123622363236423652366236723682369237023712372237323742375237623772378237923802381238223832384238523862387238823892390239123922393239423952396239723982399240024012402240324042405240624072408240924102411241224132414241524162417241824192420242124222423242424252426242724282429243024312432243324342435243624372438243924402441244224432444244524462447244824492450245124522453245424552456245724582459246024612462246324642465246624672468246924702471247224732474247524762477247824792480248124822483248424852486248724882489249024912492249324942495249624972498249925002501250225032504250525062507250825092510251125122513251425152516251725182519252025212522252325242525252625272528252925302531253225332534253525362537253825392540254125422543254425452546254725482549255025512552255325542555255625572558255925602561256225632564256525662567256825692570257125722573257425752576257725782579258025812582258325842585258625872588258925902591259225932594259525962597259825992600260126022603260426052606260726082609261026112612261326142615261626172618261926202621262226232624262526262627262826292630263126322633263426352636263726382639264026412642264326442645264626472648264926502651265226532654265526562657265826592660266126622663266426652666266726682669267026712672267326742675267626772678267926802681268226832684268526862687268826892690269126922693269426952696269726982699270027012702270327042705270627072708270927102711271227132714271527162717271827192720272127222723272427252726272727282729273027312732273327342735273627372738273927402741274227432744274527462747274827492750275127522753275427552756275727582759276027612762276327642765276627672768276927702771277227732774277527762777277827792780278127822783278427852786278727882789279027912792279327942795279627972798279928002801280228032804280528062807280828092810281128122813281428152816281728182819282028212822282328242825282628272828282928302831283228332834283528362837283828392840284128422843284428452846284728482849285028512852285328542855285628572858285928602861286228632864286528662867286828692870287128722873287428752876287728782879288028812882288328842885288628872888 \documentstyle{article} \pdftrailerid{} \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} \date{18/6/2001} \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 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 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,) 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 ... 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 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 . 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 . 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;ilen;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] \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=""} - to specify the default i-offset (i.e. query sequence offset) for transitions \item {\tt offj=""} - to specify the default j-offset (i.e. target sequence offset) for transitions \item {\tt calc=""} - 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 environment. \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}