source: tags/arb_5.5/TREEGEN/rns.c

Last change on this file was 3257, checked in by westram, 20 years ago
  • added memory.h
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.0 KB
Line 
1#include <limits.h>
2#include <stdlib.h>
3#include <memory.h>
4
5#include "rns.h"
6#include "spreadin.h"
7
8
9/* /------------------------\ */
10/* |  Erzeugung der Ur-RNS  | */
11/* \------------------------/ */
12
13int        orgLen;         /* L„nge der Ur-RNS */
14double     orgHelixPart;   /* Anteil Helix-Bereich */
15static int rnsCreated; /* Anzahl bisher erzeugter RNS-Sequenzen */
16
17/* /------------\ */
18/* |  Mutation  | */
19/* \------------/ */
20
21int    timeSteps;        /* Anzahl Zeitschritte */
22Frand  mrpb_Init,        /* Initialisierungsfunktion fr 'mutationRatePerBase' */
23       l2hrpb_Init,      /* Initialisierungsfunktion fr 'loop2helixRatePerBase' */
24       pairPart,         /* Anteil paarender Helix-Bindungen */
25       mutationRate,     /* Mutationsrate */
26       splitRate,        /* Spaltungsrate */
27       helixGcDruck,     /* G-C-Druck       im Helix-Bereich */
28       helixGcRate,      /* Verh„ltnis G:C  im Helix-Bereich */
29       helixAtRate,      /* Verh„ltnis A:T  im Helix-Bereich */
30       loopGcDruck,      /* G-C-Druck       im Loop-Bereich */
31       loopGcRate,       /* Verh„ltnis G:C  im Loop-Bereich */
32       loopAtRate;       /* Verh„ltnis A:T  im Loop-Bereich */
33double transitionRate,   /* Transition-Rate */
34       transversionRate; /* Transversion-Rate */
35
36static double     *mutationRatePerBase,   /* positionsspez. Mutationsrate (wird nur einmal bestimmt und bleibt dann konstant) */
37                  *loop2helixRatePerBase; /* positionsspez. Rate fr Wechsel Loop-Base in Helix-Base und vv. (wird nur einmal bestimmt und bleibt dann konstant) */
38static int         mrpb_anz,              /* Anzahl Positionen */
39                   mrpb_allocated,        /* wirklich Gr”áe des Arrays */
40                   l2hrpb_anz,            /* Anzahl Positionen */
41                   l2hrpb_allocated;      /* wirklich Gr”áe des Arrays */
42static DoubleProb  helixMutationMatrix,   /* Mutationsmatrix fr Helix-Bereiche */
43                   loopMutationMatrix;    /* Mutationsmatrix fr Loop-Bereiche */
44
45/* /----------------------\ */
46/* |  Ausgabefilepointer  | */
47/* \----------------------/ */
48
49FILE *topo, /* Topologie */
50     *seq;  /* Sequenzen */
51
52/* /-------------\ */
53/* |  Sonstiges  | */
54/* \-------------/ */
55
56static int minDepth = INT_MAX, /* minimale Tiefe (Astanzahl) der Blattspitzen */
57           maxDepth = INT_MIN; /* maximale Tiefe der Blattspitzen */
58
59/* -------------------------------------------------------------------------- */
60/*      void dumpDepths(void) */
61/* ------------------------------------------------------ 24.05.95 22.27 ---- */
62void dumpDepths(void)
63{
64    printf("Minimale Baumtiefe = %i\n", minDepth);
65    printf("Maximale Baumtiefe = %i\n", maxDepth);
66}
67/* -------------------------------------------------------------------------- */
68/*      static void dumpRNS(RNS rns) */
69/* ------------------------------------------------------ 26.05.95 11.29 ---- */
70static void dumpRNS(RNS rns)
71{
72    int        b,
73               b1,
74               b2;
75    static int cleared,
76               h_cnt[BASETYPES+1][BASETYPES+1],
77               l_cnt[BASETYPES+1],
78               loop,
79               helix;
80
81    if (!cleared)
82    {
83        for (b1 = 0; b1<(BASETYPES+1); b1++)
84        {
85            for (b2 = 0; b2<(BASETYPES+1); b2++) h_cnt[b1][b2] = 0;
86            l_cnt[b1] = 0;
87        }
88
89        loop  = 0;
90        helix = 0;
91
92        cleared = 1;
93    }
94
95    if (rns)
96    {
97        for (b = 0; b<(rns->bases); b++)
98        {
99            char base = rns->base[b];
100
101            if (isHelical(base))
102            {
103                int bt1 = char2BaseType(base),
104                    bt2 = char2BaseType(rns->base[b+1]);
105
106                h_cnt[bt1][bt2]++;
107                helix++;
108                b++;
109            }
110            else
111            {
112                int bt = char2BaseType(base);
113
114                l_cnt[bt]++;
115                loop++;
116            }
117        }
118    }
119    else
120    {
121        printf("Helix-Basenpaare = %i\n"
122               "Loop-Basen       = %i\n"
123               "Helix:Loop       = %f\n",
124               helix,
125               loop,
126               (double)helix/(double)loop);
127
128        {
129            int gc      = h_cnt[BASE_C][BASE_G]+h_cnt[BASE_G][BASE_C],
130                at      = h_cnt[BASE_A][BASE_T]+h_cnt[BASE_T][BASE_A],
131                paarend = gc+at;
132
133            printf("GC-Paare              = %i\n"
134                   "AT-Paare              = %i\n"
135                   "Paare:Helix-Bindungen = %f\n"
136                   "GC-Paare:Paare        = %f\n",
137                   gc,
138                   at,
139                   (double)paarend/(double)helix,
140                   (double)gc/(double)paarend);
141        }
142
143        printf("\n");
144    }
145}
146/* -------------------------------------------------------------------------- */
147/*      static void initBaseSpecificProbs(int bases) */
148/* ------------------------------------------------------ 24.05.95 12.51 ---- */
149static void initBaseSpecificProbs(int bases)
150{
151    int b;
152
153    mrpb_anz            = bases;
154    mrpb_allocated      = bases;
155    mutationRatePerBase = malloc(bases*sizeof(double));
156
157    l2hrpb_anz            = bases;
158    l2hrpb_allocated      = bases;
159    loop2helixRatePerBase = malloc(bases*sizeof(double));
160
161    if (!mutationRatePerBase || !loop2helixRatePerBase) outOfMemory();
162
163    for (b = 0; b<bases; b++)
164    {
165        mutationRatePerBase[b]   = getFrand(mrpb_Init);
166        loop2helixRatePerBase[b] = getFrand(l2hrpb_Init);
167    }
168}
169/* -------------------------------------------------------------------------- */
170/*      static RNS allocRNS(int len) */
171/* ------------------------------------------------------ 20.05.95 16.04 ---- */
172static RNS allocRNS(int len)
173{
174    RNS rns = malloc(sizeof(*rns));
175
176    if (!rns) outOfMemory();
177
178/*     rns->bases = orgLen; */
179/*     rns->base  = malloc(sizeof(*(rns->base))*orgLen); */
180    rns->bases = len;
181    rns->base  = malloc(sizeof(*(rns->base))*len);
182
183    if (!rns->base) outOfMemory();
184
185    return rns;
186}
187/* -------------------------------------------------------------------------- */
188/*      RNS createOriginRNS(void) */
189/* ------------------------------------------------------ 14.05.95 14:54 ---- */
190/* */
191/*  Erzeugt eine Ur-RNS */
192/* */
193RNS createOriginRNS(void)
194{
195    RNS rns      = allocRNS(orgLen);
196    int helixLen = orgLen*orgHelixPart,
197        l;
198    str base     = rns->base;
199
200    printf("Generating origin species..\n");
201
202    initBaseSpecificProbs(orgLen);
203
204    rns->laufNr = rnsCreated++;
205
206    /* /------------------\ */
207    /* |  Helix erzeugen  |                                                             im Loop-Bereich */
208    /* \------------------/ */
209
210    if (helixLen%1) helixLen--;                            /* muá gerade L„nge haben, da nur Paare! */
211
212    assert(helixLen<=orgLen);
213
214    rns->helix   = helixLen/2;
215    rns->pairing = 0;
216
217    {
218        DoubleProb orgHelixProb;
219        Spreading  s;
220        int        b1,
221                   b2;
222        double     actPairPart     = getFrand(pairPart),
223                   actHelixGcDruck = getFrand(helixGcDruck),
224                   actHelixGcRate  = getFrand(helixGcRate),
225                   actHelixAtRate  = getFrand(helixAtRate),
226                   nonPairProb     = (1.0-actPairPart)/2.0;
227
228        for (b1 = 0; b1<BASETYPES; b1++)
229        {
230            for (b2 = 0; b2<BASETYPES; b2++)
231            {
232                if (isPairing(b1, b2))
233                {
234                    switch (b1)
235                    {
236                        case BASE_A:
237                        case BASE_T:
238                        {
239                            orgHelixProb[b1][b2] = (actPairPart*(1.0-actHelixGcDruck))/2.0;
240                            break;
241                        }
242                        case BASE_C:
243                        case BASE_G:
244                        {
245                            orgHelixProb[b1][b2] = (actPairPart*actHelixGcDruck)/2.0;
246                            break;
247                        }
248                    }
249                }
250                else
251                {
252                    double prob = nonPairProb;
253                    int    b    = b1;
254
255                    while (1)                              /* wird je einmal mit b1 und b2 ausgefhrt */
256                    {
257                        switch (b)
258                        {
259                            case BASE_A:
260                            {
261                                prob = prob*(1.0-actHelixGcDruck)*actHelixAtRate;
262                                break;
263                            }
264                            case BASE_C:
265                            {
266                                prob = prob*actHelixGcDruck*(1.0-actHelixGcRate);
267                                break;
268                            }
269                            case BASE_G:
270                            {
271                                prob = prob*actHelixGcDruck*actHelixGcRate;
272                                break;
273                            }
274                            case BASE_T:
275                            {
276                                prob = prob*(1.0-actHelixGcDruck)*(1.0-actHelixAtRate);
277                                break;
278                            }
279                        }
280
281                        if (b==b2) break;
282                        b = b2;
283                    }
284
285                    orgHelixProb[b1][b2] = prob;
286                }
287            }
288        }
289
290        s = newSpreading((double*)orgHelixProb, BASEQUAD);
291
292        for (l = 0; l<helixLen; l+=2)
293        {
294            int val = spreadRand(s),
295                B1  = val%BASETYPES,
296                B2  = val/BASETYPES;
297
298            base[l]   = helixBaseChar[B1];
299            base[l+1] = helixBaseChar[B2];
300
301            rns->pairing += isPairing(B1, B2);
302        }
303
304        freeSpreading(s);
305    }
306
307    /* /-----------------\ */
308    /* |  Loop erzeugen  | */
309    /* \-----------------/ */
310
311    {
312        SingleProb orgLoopProb;
313        Spreading  s;
314        double     actLoopGcDruck = getFrand(loopGcDruck),
315                   actLoopGcRate  = getFrand(loopGcRate),
316                   actLoopAtRate  = getFrand(loopAtRate);
317
318        orgLoopProb[BASE_A] = (1.0-actLoopGcDruck)*actLoopAtRate;
319        orgLoopProb[BASE_C] = actLoopGcDruck*(1.0-actLoopGcRate);
320        orgLoopProb[BASE_G] = actLoopGcDruck*actLoopGcRate;
321        orgLoopProb[BASE_T] = (1.0-actLoopGcDruck)*(1.0-actLoopAtRate);
322
323        s = newSpreading((double*)orgLoopProb, BASETYPES);
324        for (; l<orgLen; l++) base[l] = loopBaseChar[spreadRand(s)];
325        freeSpreading(s);
326    }
327
328    return rns;
329}
330/* -------------------------------------------------------------------------- */
331/*      void freeRNS(RNS rns) */
332/* ------------------------------------------------------ 20.05.95 19.45 ---- */
333void freeRNS(RNS rns)
334{
335    free(rns->base);
336    free(rns);
337}
338/* -------------------------------------------------------------------------- */
339/*      static RNS dupRNS(RNS rns) */
340/* ------------------------------------------------------ 20.05.95 20.32 ---- */
341static RNS dupRNS(RNS rns)
342{
343    RNS neu = malloc(sizeof(*rns));
344
345    if (!neu) outOfMemory();
346
347    memcpy(neu, rns, sizeof(*rns));
348
349    neu->base = malloc(rns->bases*sizeof(*(neu->base)));
350    memcpy(neu->base, rns->base, rns->bases);
351
352    neu->laufNr = rnsCreated++;
353
354    return neu;
355}
356/* -------------------------------------------------------------------------- */
357/*      static void dumpDoubleProb(double *d, int anz) */
358/* ------------------------------------------------------ 25.05.95 01.31 ---- */
359/*static void dumpDoubleProb(double *d, int anz) */
360/*{ */
361/*    while (anz--) printf("%-10f", *d++); */
362/*    printf("\n\n"); */
363/*} */
364/* -------------------------------------------------------------------------- */
365/*      static void calcMutationMatrix(DoubleProb mutationMatrix, double mu... */
366/* ------------------------------------------------------ 24.05.95 13.58 ---- */
367static void calcMutationMatrix(DoubleProb mutationMatrix, double muteRate, double gcDruck, double gcRate, double atRate, double *pairProb)
368{
369    double k   = transitionRate/transversionRate,
370           fa  = (1.0-gcDruck)*atRate,
371           fc  = gcDruck*(1.0-gcRate),
372           fg  = gcDruck*gcRate,
373           ft  = (1.0-gcDruck)*(1.0-atRate),
374           bfa = transversionRate*fa,
375           bfc = transversionRate*fc,
376           bfg = transversionRate*fg,
377           bft = transversionRate*ft,
378           kag = k/(fa+fg),
379           kct = k/(fc+ft);
380/*           sa  = (kag+3.0)*bfa,                            // Summe der "mutierenden" Positionen jeder Zeile */
381/*           sc  = (kct+3.0)*bfc, */
382/*           sg  = (kag+3.0)*bfg, */
383/*           st  = (kct+3.0)*bft; */
384
385    /* Auf aktuelle Mutationsrate normieren */
386
387/*    bfa = bfa*muteRate/sa; */
388/*    bfc = bfc*muteRate/sc; */
389/*    bfg = bfg*muteRate/sg; */
390/*    bft = bft*muteRate/st; */
391
392/*    printf("bfa=%f\n", bfa); */
393/*    printf("bfc=%f\n", bfc); */
394/*    printf("bfg=%f\n", bfg); */
395/*    printf("bft=%f\n", bft); */
396
397    /* Matrix besetzen */
398
399    mutationMatrix[BASE_A][BASE_A] = 1.0-(kag+3.0)*bfa;
400    mutationMatrix[BASE_C][BASE_A] = bfa;
401    mutationMatrix[BASE_G][BASE_A] = (kag+1.0)*bfa;
402    mutationMatrix[BASE_T][BASE_A] = bfa;
403
404    mutationMatrix[BASE_A][BASE_C] = bfc;
405    mutationMatrix[BASE_C][BASE_C] = 1.0-(kct+3.0)*bfc;
406    mutationMatrix[BASE_G][BASE_C] = bfc;
407    mutationMatrix[BASE_T][BASE_C] = (kct+1.0)*bfc;
408
409    mutationMatrix[BASE_A][BASE_G] = (kag+1.0)*bfg;
410    mutationMatrix[BASE_C][BASE_G] = bfg;
411    mutationMatrix[BASE_G][BASE_G] = 1.0-(kag+3.0)*bfg;
412    mutationMatrix[BASE_T][BASE_G] = bfg;
413
414    mutationMatrix[BASE_A][BASE_T] = bft;
415    mutationMatrix[BASE_C][BASE_T] = (kct+1.0)*bft;
416    mutationMatrix[BASE_G][BASE_T] = bft;
417    mutationMatrix[BASE_T][BASE_T] = 1.0-(kct+3.0)*bft;
418
419/*    { */
420/*        int von,                                         // Matrix ausgeben */
421/*            nach; */
422/* */
423/*        printf("       von %c     von %c     von %c     von %c     \n", */
424/*                helixBaseChar[0], */
425/*                helixBaseChar[1], */
426/*                helixBaseChar[2], */
427/*                helixBaseChar[3] ); */
428/* */
429/*        for (nach = BASE_A; nach<=BASE_T; nach++) */
430/*        { */
431/*            double sum = 0.0; */
432/* */
433/*            printf("nach %c ", helixBaseChar[nach]); */
434/* */
435/*            for (von = BASE_A; von<=BASE_T; von++) */
436/*            { */
437/*                printf("%-10f", mutationMatrix[von][nach]); */
438/*                sum += mutationMatrix[von][nach]; */
439/*            } */
440/* */
441/*            printf("    sum = %-10f\n", sum); */
442/*        } */
443/* */
444/*        printf("\n"); */
445/*    } */
446
447    if (pairProb)                                          /* soll pairProb berechnet werden? */
448    {
449        double mutatesTo[BASETYPES],
450               freq[BASETYPES];                            /* H„ufigkeit der einzelnen Basen */
451        int    von,
452               nach;
453
454        freq[BASE_A] = fa;
455        freq[BASE_C] = fc;
456        freq[BASE_G] = fg;
457        freq[BASE_T] = ft;
458
459        for (nach = 0; nach<BASETYPES; nach++)
460            mutatesTo[nach] = 0.0;
461
462        for (von = 0; von<BASETYPES; von++)
463            for (nach = 0; nach<BASETYPES; nach++)
464                mutatesTo[nach] += mutationMatrix[von][nach]*freq[von];
465
466        *pairProb = 2.0*mutatesTo[BASE_A]*mutatesTo[BASE_T] + 2.0*mutatesTo[BASE_C]*mutatesTo[BASE_G];
467    }
468}
469/* -------------------------------------------------------------------------- */
470/*      static int calcPairTrials(double pairProb, double actPairPart) */
471/* ------------------------------------------------------ 25.05.95 13.31 ---- */
472/* */
473/*  Berechnet die Anzahl Mutations-Wiederholungen, die notwendig sind, um */
474/*  mindestens 'actPairPart' Prozent paarende Bindungen zu erhalten, falls */
475/*  die Wahrscheinlichkeit eine paarende Bindung zu erzeugen gleich */
476/*  'pairProb' ist. */
477/* */
478static int calcPairTrials(double pairProb, double actPairPart)
479{
480    int    trials   = 1;
481    double failProb = 1.0-pairProb,
482           succProb = pairProb;
483
484    while (succProb<actPairPart)
485    {
486        pairProb *= failProb;
487        succProb += pairProb;
488        trials++;
489
490/*        printf("trials=%i succProb=%f actPairPart=%f\n", trials, succProb, actPairPart); */
491    }
492
493    return trials;
494}
495/* -------------------------------------------------------------------------- */
496/*      static void indent(int depth) */
497/* ------------------------------------------------------ 24.05.95 21.08 ---- */
498/*static void indent(int depth) */
499/*{ */
500/*    while (depth--) fputc(' ', topo); */
501/*} */
502/* -------------------------------------------------------------------------- */
503/*      static void mutateRNS(RNS rns, int steps, int depth) */
504/* ------------------------------------------------------ 20.05.95 19.50 ---- */
505/* */
506/*  Mutiert eine RNS bis zur n„chsten Spaltung */
507/* */
508/*  'steps'     Anzahl noch zu durchlaufender Zeitschritte */
509/* */
510static void mutateRNS(int no_of_father, RNS rns, int steps, int depth)
511{
512    int    splitInSteps,
513           s;
514    double mutationTime = 0.0;
515
516    /* /---------------------------------------\ */
517    /* |  Schritte bis zur Spaltung berechnen  | */
518    /* \---------------------------------------/ */
519
520    {
521        double actualSplitRate = getFrand(splitRate);
522
523        assert(actualSplitRate!=0);
524
525        splitInSteps = (int)(1.0/actualSplitRate);
526        if (splitInSteps>steps) splitInSteps = steps;
527
528        assert(splitInSteps>=1);
529    }
530
531    /* /----------------------------\ */
532    /* |  Zeitschritte durchlaufen  | */
533    /* \----------------------------/ */
534
535    for (s = 0; s<splitInSteps; s++)
536    {
537        int       b,
538                  pairTrials;                              /* Anzahl Versuche eine paarende Helixbindung herzustellen */
539        double    actMutationRate    = getFrand(mutationRate),
540                  actPairPart        = getFrand(pairPart);
541        Spreading s_helix[BASETYPES],
542                  s_loop[BASETYPES];
543
544        {
545            double pairProb;                               /* Wahrscheinlichkeit, daá ein Paar im helikalen Bereich entsteht */
546
547            calcMutationMatrix(helixMutationMatrix, 1.0, getFrand(helixGcDruck), getFrand(helixGcRate), getFrand(helixAtRate), &pairProb);
548            calcMutationMatrix(loopMutationMatrix, actMutationRate, getFrand(loopGcDruck), getFrand(loopGcRate), getFrand(loopAtRate), NULL);
549
550            pairTrials = calcPairTrials(pairProb, actPairPart);
551/*            printf("pairProb=%f pairTrials=%i\n", pairProb, pairTrials); */
552        }
553
554        for (b = 0; b<BASETYPES; b++)
555        {
556            s_helix[b] = newSpreading(&(helixMutationMatrix[b][0]), BASETYPES);
557            s_loop[b]  = newSpreading(&(loopMutationMatrix[b][0]), BASETYPES);
558        }
559
560        mutationTime += actMutationRate;                   /* Mutationszeit aufaddieren (Einheit ist Mutationsrate*Zeitschritte) */
561
562        /* /----------------------------------\ */
563        /* |  Alle Basen(-paare) durchlaufen  | */
564        /* \----------------------------------/ */
565
566        for (b = 0; b<(rns->bases); )
567        {
568            char base = rns->base[b];
569
570            if (!isDeleted(base))                          /* Deletes ignorieren */
571            {
572                /* /---------------------\ */
573                /* |  Helicale Bereiche  | */
574                /* \---------------------/ */
575
576                if (isHelical(base))
577                {
578                    int  trials = pairTrials,
579                         mut1   = randProb()<mutationRatePerBase[b]*actMutationRate,
580                         mut2   = randProb()<mutationRatePerBase[b+1]*actMutationRate;
581                    char base2  = rns->base[b+1];
582
583                    assert(isHelical(base2));
584
585                    if (mut1 || mut2)
586                    {
587                        int bt1 = char2BaseType(base),
588                            bt2 = char2BaseType(base2);
589
590                        if (isPairing(bt1, bt2))
591                        {
592                            rns->pairing--;
593                        }
594
595                        while (trials--)
596                        {
597                            if (mut1)
598                            {
599                                if (mut2)                  /* beide Basen mutieren */
600                                {
601                                    bt1 = spreadRand(s_helix[bt1]);
602                                    bt2 = spreadRand(s_helix[bt2]);
603                                }
604                                else                       /* nur 1.Base mutieren */
605                                {
606                                    bt1 = spreadRand(s_helix[bt1]);
607                                }
608                            }
609                            else                           /* nur 2.Base mutieren */
610                            {
611                                bt2 = spreadRand(s_helix[bt2]);
612                            }
613
614                            if (isPairing(bt1, bt2))       /* paarend? ja->abbrechen */
615                            {
616                                rns->pairing++;
617                                break;
618                            }
619                        }
620
621                        rns->base[b]   = helixBaseChar[bt1];
622                        rns->base[b+1] = helixBaseChar[bt2];
623                    }
624
625                    b++;
626                }
627
628                /* /-----------------\ */
629                /* |  Loop-Bereiche  | */
630                /* \-----------------/ */
631
632                else
633                {
634                    double mutationProb = actMutationRate*mutationRatePerBase[b];
635
636                    if (randProb()<mutationProb)
637                    {
638                        rns->base[b] = loopBaseChar[spreadRand(s_loop[char2BaseType(base)])];
639                    }
640                }
641            }
642
643            b++;
644        }
645
646        for (b = 0; b<BASETYPES; b++)
647        {
648            freeSpreading(s_helix[b]);
649            freeSpreading(s_loop[b]);
650        }
651    }
652
653    splitRNS(no_of_father, rns, mutationTime, steps-splitInSteps, depth+1);
654}
655/* -------------------------------------------------------------------------- */
656/*      void splitRNS(RNS origin, double age, int steps, int depth) */
657/* ------------------------------------------------------ 20.05.95 20.13 ---- */
658/* */
659/*  Spaltet eine RNS in zwei Species auf */
660/* */
661void splitRNS(int no_of_father, RNS origin, double age, int steps, int depth)
662{
663    int x;
664
665    dumpRNS(origin);
666
667    /* /---------------------\ */
668    /* |  Sequenz schreiben  | */
669    /* \---------------------/ */
670
671    if (no_of_father != -1) {
672        fprintf(seq, ">no%i son of no%i\n", origin->laufNr, no_of_father);
673    }
674    else {
675        fprintf(seq, ">no%i father of all species\n", origin->laufNr);
676    }
677    no_of_father = origin->laufNr; /* now i'm the father */
678    for (x = 0; x<(origin->bases); x++) fputc(origin->base[x], seq);
679    fputc('\n', seq);
680
681    if (steps)                                             /* Species splitten! */
682    {
683        double gcDruck_val      = helixGcDruck->val,       /* Frand-Werte merken */
684               pairPart_val     = pairPart->val,
685               mutationRate_val = mutationRate->val,
686               splitRate_val    = splitRate->val;
687
688/*        indent(depth); */
689        fprintf(topo, "(no%i:%f,\n", origin->laufNr, age);
690
691        {
692            RNS left = dupRNS(origin);                     /* linker Sohn */
693
694            mutateRNS(no_of_father, left, steps, depth);
695            freeRNS(left);
696        }
697
698        fputs(",\n", topo);
699
700        helixGcDruck->val = gcDruck_val;                   /* Frand-Werte wiederherstellen */
701        pairPart->val     = pairPart_val;
702        mutationRate->val = mutationRate_val;
703        splitRate->val    = splitRate_val;
704
705        {
706            RNS right = dupRNS(origin);                    /* rechter Sohn */
707
708            mutateRNS(no_of_father, right, steps, depth);
709            freeRNS(right);
710        }
711
712        fputc(')', topo);
713    }
714    else                                                   /* Baumspitze */
715    {
716        if      (depth>maxDepth) maxDepth = depth;
717        else if (depth<minDepth) minDepth = depth;
718
719/*        indent(depth); */
720        fprintf(topo, "no%i:%f", origin->laufNr, age);
721
722        if ((origin->laufNr%100) == 0) {
723            printf("generated Species: %i\n", origin->laufNr);
724        }
725    }
726
727    if (age==0.0) dumpRNS(NULL);
728}
729
730
Note: See TracBrowser for help on using the repository browser.