source: tags/ms_r16q2/TREEGEN/rns.c

Last change on this file was 7811, checked in by westram, 14 years ago

merge from dev [7748] [7749] [7750]

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