0
$\begingroup$

I'm a beginner mathematica programmer and I am tryin to create my first demonstration. This code should represent what happens if one changes the size of the reading frame, calculating the nucleotide percentage in random DNA.

When I run it, it stays in an infinite loop.

Manipulate[
 SeedRandom[1234];
 Dnasize = 1000;
 dnalist = RandomChoice[{A, T, G, C}, Dnasize];
 plist = Partition[dnalist, partition];
 iterationmax = Dnasize/partition;
 As = {};
 Ts = {};
 Gs = {};
 Cs = {};
 GCs = {};
 For[i = 1, i < iterationmax + 1, i++,
   AppendTo[As, Count[plist[[i]], A]];
   AppendTo[Gs, Count[plist[[i]], G]];
   AppendTo[Ts, Count[plist[[i]], T]];
   AppendTo[Cs, Count[plist[[i]], C]];
   AppendTo[
    GCs, (Gs[[i]] + Cs[[i]])/(As[[i]] + Gs[[i]] + Cs[[i]] + Ts[[i]])]
   ]

  ListLinePlot[GCs, 
   PlotLegends -> {"GC content across the molecule"}], {partition, 50,
   100}
 ]

Any help appreciated.

$\endgroup$
2
$\begingroup$

The first problem with your code is that you have not restricted partition to have integer values. This is easily fixed by specifying a step size of 1 for the control: {partition, 50, 100, 1}

The second problem is that when Dnasize isn't divisible by partitions you have too large a value for iterationmax. For example if partitions is 51 then plist will contain 19 sets of 51. So you need 19 iterations of the for loop, but iterationmax is 1000/51 = 19.6 (and therefore the loop will not terminate until i=20) You could instead use iterationmax = Length[plist] to ensure that the correct value is always used.

These two corrections will give you a working code, but perhaps you might be interested in some additional improvements. Instead of going through the elements of plist in a for loop, use Map:

SeedRandom[1234];
Dnasize = 1000;
dnalist = RandomChoice[{A, T, G, C}, Dnasize];

Manipulate[
 plist = Partition[dnalist, partition];
 As = Count[A] /@ plist;
 Gs = Count[G] /@ plist;
 Cs = Count[C] /@ plist;
 Ts = Count[T] /@ plist;
 GCs = (Gs + Cs)/(As + Gs + Cs + Ts);
 ListLinePlot[GCs, PlotLegends -> {"GC content across the molecule"}],
 {partition, 50, 100, 1}]

You could also go a step further and use Outer to do all the counts in one expression:

Manipulate[
 plist = Partition[dnalist, partition];
 {As, Gs, Ts, Cs} = Outer[Count[#2, #1] &, {A, G, T, C}, plist, 1];
 GCs = (Gs + Cs)/(As + Gs + Cs + Ts);
 ListLinePlot[GCs, PlotLegends -> {"GC content across the molecule"}],
 {partition, 50, 100, 1}]
$\endgroup$
  • $\begingroup$ Very good usage of "Outer"! Clever :-) $\endgroup$ – mgamer Dec 25 '14 at 22:57
  • $\begingroup$ thank you! This was very helpful also. $\endgroup$ – user3523464 Dec 27 '14 at 20:10
1
$\begingroup$

Welcome user3523464... There is a lot to do with the code. First you may notice that "c" in your for-loop ist in black color and not in blue. This is du to "C" is a reserved word in Mathematica. So you have to use characters here. Also you will run into trouble with the Partition command, since there will be lists of different length and you´ll get problems in the for-loop.. (I made a work-sound in the Dnasize). Also AppendTo is not the fastest Function, but in this application it does´t matter ;-)

Look at the list and string command of Mathematica thy are really powerful, besides this it is a good idea to begin all of "your" variable names with a lowercase character to avoid mixing up with reserved words.

Since you are a beginner in Mathematica I give you an example that works that is as close as possible to your code (the code is readable but far from being optimal). Put all the calculation outside the manipulate, that makes things easier.

calculate[GCs_, partition_] := Module[ {gtemp = GCs, ptemp = partition, Dnasize, As = {}, Ts = {}, Gs = {}, Cs = {}, dnalist, plist, iterationmax, i }, Dnasize = 50*ptemp; SeedRandom[1234]; dnalist = RandomChoice[{"A", "T", "G", "C"}, Dnasize]; plist = Partition[dnalist, ptemp]; iterationmax = Dnasize/ptemp; For[i = 1, i < iterationmax + 1, i++, AppendTo[As, Count[plist[[i]], "A"]]; AppendTo[Gs, Count[plist[[i]], "G"]]; AppendTo[Ts, Count[plist[[i]], "T"]]; AppendTo[Cs, Count[plist[[i]], "C"]]; AppendTo[ gtemp, (Gs[[i]] + Cs[[i]])/(As[[i]] + Gs[[i]] + Cs[[i]] + Ts[[i]])]] ; gtemp ]

Now the manipulation is easy

Manipulate[ ListLinePlot[calculate[GCs, partition], PlotLegends -> {"GC content across the molecule"}], {partition, 50, 100, 5}, (*Local variables*) {GCs, {}, None}]

$\endgroup$
  • $\begingroup$ wow. Thank you very much, I will study this code thorougly. $\endgroup$ – user3523464 Dec 25 '14 at 15:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.