source: branches/lib/TREEGEN/rns.c

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