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 | int orgLen; // Laenge der Ur-RNS |
---|
13 | double orgHelixPart; // Anteil Helix-Bereich |
---|
14 | static int rnsCreated; // Anzahl bisher erzeugter RNS-Sequenzen |
---|
15 | |
---|
16 | // ----------------- |
---|
17 | // Mutation |
---|
18 | |
---|
19 | int timeSteps; // Anzahl Zeitschritte |
---|
20 | Frand 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 |
---|
31 | double transitionRate, // Transition-Rate |
---|
32 | transversionRate; // Transversion-Rate |
---|
33 | |
---|
34 | static 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) |
---|
36 | static int mrpb_anz, // Anzahl Positionen |
---|
37 | mrpb_allocated, // wirklich Groesse des Arrays |
---|
38 | l2hrpb_anz, // Anzahl Positionen |
---|
39 | l2hrpb_allocated; // wirklich Groesse des Arrays |
---|
40 | static DoubleProb helixMutationMatrix, // Mutationsmatrix fuer Helix-Bereiche |
---|
41 | loopMutationMatrix; // Mutationsmatrix fuer Loop-Bereiche |
---|
42 | |
---|
43 | // --------------------------- |
---|
44 | // Ausgabefilepointer |
---|
45 | |
---|
46 | FILE *topo, // Topologie |
---|
47 | *seq; // Sequenzen |
---|
48 | |
---|
49 | // ------------------ |
---|
50 | // Sonstiges |
---|
51 | |
---|
52 | static int minDepth = INT_MAX, // minimale Tiefe (Astanzahl) der Blattspitzen |
---|
53 | maxDepth = INT_MIN; // maximale Tiefe der Blattspitzen |
---|
54 | |
---|
55 | void dumpDepths() |
---|
56 | { |
---|
57 | printf("Minimale Baumtiefe = %i\n", minDepth); |
---|
58 | printf("Maximale Baumtiefe = %i\n", maxDepth); |
---|
59 | } |
---|
60 | static 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 | } |
---|
136 | static 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 | } |
---|
156 | static 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 | } |
---|
169 | RNS 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 | } |
---|
305 | void freeRNS(RNS rns) |
---|
306 | { |
---|
307 | free(rns->base); |
---|
308 | free(rns); |
---|
309 | } |
---|
310 | static 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 | } |
---|
325 | static void calcMutationMatrix(DoubleProb mutationMatrix, 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 | } |
---|
383 | static 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 | } |
---|
402 | static 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 | } |
---|
543 | void 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 | |
---|