a probability of substituting any given position in the string with A C T G

So, in each clone, what's the error (mutation) limit?, does the mutation position hold any significance? do you need to filter your genomes afterwards to select only those still-functional ones that were not corrupted by the mutations but where the mutation may present phenotypic differences?

My initial thought was to make a mutlidimensional array, where each successive round of replication would be represented as values in the new array dimension

In light of the above, hidden markov models can use the transition probabilities while performing the mutations (maybe difficult when wanting to consider all the different cases out there), your best bet is to grab a statistician or to keep one handy. HMMs have been extensively used before in sequence alignment and detecting indels hence this scenario here is related.

A simpler approach is to use randomization to simulate mutations, check this discussion and its implementation or to consult the feasibility of using an already available solution.

UPDATE:Scanning the simulation program in the links I gave you; it seems to have a loose criteria on what base can be replaced with what so I followed another approach by introducing conditions on mutation types (favoring transitions rather than fatal transversions) and giving you control over the number of replication iterations then showing you the sequences aligned in proximity for good visuals

# a program that generates a 50 bp random DNA seq 
#and introduces 1 mutation per DNA replication
use strict;
use warnings;
use Data::Dumper;

#generating a 50 bp seq one base at a time.
#you can read sets of seqs from a file instead
my ($base,$virusGenome);
my @nucleotide = ('a','t','g','c');
for(1..50){
    $base = $nucleotide[rand(3)];    
    $virusGenome .=$base;
    }

my @seqArray; 
#a hash table of the permissible mutations
#no freqs rates, note the base case.
my %mutationBounds=('a'=>'T','c'=>'G','g'=>'C','t'=>'A'); 

print "how many replications?";
chomp(my $rounds=<>);
mutate(times=>$rounds, seq=>$virusGenome);
sub mutate{
    my %params=@_;
    print "original seq=> $params{seq}\n";
    for my $i (0..$params{times}-1){
        #convert seq from string to array 
        @seqArray = split("",$params{seq});

        #pick a base from a random position 
        my $randPos = int(rand(length($params{seq})));
        my $toReplace=$seqArray[$randPos]; 
        
        #mutate based on %mutationBounds
        # assign to a seq the altered $virusGenome. 
        my $seqMute = $params{seq};    
        substr ($seqMute, $randPos,1,$mutationBounds{$toReplace});    
+    
        print "mutant seq",$i+1,"=> ", $seqMute,"\n";

        #which rep cycle is done. What was replaced where.
        #print "Cycle " ,$i+1, " $toReplace replaced w $mutationBounds
+{$toReplace} @ $randPos\n";
        }
    };


#print Dumper(\%mutationBounds);
#print Dumper(\%params);
#print Dumper(\@seqArray);


##OUTPUT 
hisham@bioslayer:~/Desktop$ perl virusToy.pl
how many replications?3
original seq=> gaataggaataatagggagattgtaagagagtgggttagtatagttaata
mutant seq1=> gaataggaataatagggagattgtaagagagtgggttagtatagttaaAa
mutant seq2=> gaTtaggaataatagggagattgtaagagagtgggttagtatagttaata
mutant seq3=> gaataggaataataCggagattgtaagagagtgggttagtatagttaata
[download]

I just now saw this script that you added. Your knowledge of biology tells me you work in a related field. Most of the parameters such as rate, and a mutation probability table I will be taking directly from wet lab/ primary literature findings. I have some good estimates in mind. I am aware of the many issues with modeling biological systems, but in this case I have several natural test sets. This toy model is to be used only to assess the effect of uneven and perhaps biased sequence sampling. I will go over this script you sent tonight it looks much better than the one I made haha.

my script overlooked the mutation probability scores per base (unperlish laziness and lack of intuition at that moment) whereas yours has focused mainly on that, my approach was to provide room for the more common mutations rather than consulting with lesser common ones that may or may not be categorized as nonsense. While the user has a say at providing the number of replications to be conducted (a more prudent approach will have to consider that as the replication is repeated over and over there are more chances for errors to show, but my program and so is yours it seems were more obsessed with making mutations happen at the first level of replication and not based on the hierarchy that you indicated in your OP that each seq may generate more seq in a tree fashion)...

I find your incorporation of the mutation table so marvelous yet I have some comments on the implementation however that I'll share after running your program. Combining these two programs though is a worthwhile learning opportunity.

p.s. Yes, I work as a bioinformatician

---

David R. Gergen said "We know that second terms have historically been marred by hubris and by scandal." and I am a two y.o. monk today :D, June,12th, 2011...

Even if every star in every galaxy in the universe had a duplicate earth, there wouldn't be enough disks & tape on all of them combined to store all the progeny of 1000 generations of a single 1000-codon sequence.

Taking a 1000-codon sequence through one possible lineage for 1000 generations takes maybe a millisecond, so you could fill a disk with possible end points very quickly.

But without some selection criteria, any single lineage of random mutations is as likely as any other, so there is nothing to be learnt from statistically sampling the final generation. All you'd be testing is the distribution of your source of randomness.

So the question becomes, what were you hoping to learn from the outcome of this experiment? If we knew what you were hoping to investigate, we might be able to suggest plausible approaches.

You want to do 1000 replications? If I understand you correctly that would amount to 5^1000 sequences, a number too big for my calculator. 5^100 is already a number with 69 digits

I hope you have a plan B ;-)

Earnestly, what do you want to do with those sequences eventually ? If you want some statistics out of them, maybe you can get away with math instead of storing all those sequences

Do you care about the individual nodes, or just each current genome sequence (not a geneticist, so pardon any terminology mess-ups) and the number of each in existence?

You are dealing with large numbers of values at the mutation rates and generations in your example. As was previously mentioned, 5**1000 is quite large. Is there any dieoff rate or fitness test at each generation that could be used to prune the dataset?

Perhaps a more manageable structure (although not at these rates and generations) would be a Hash or HoH (hash-of-hashes), where each key is the ACTG sequence, and the value is a hash with the current number of virii with that sequence, as well as some other bookkeeping information. If you don't need other bookkeeping information, perhaps a simple hash would do. Worst storage case on this (if I am thinking through this correctly) should be no worse (given a scaling factor) than an array of all of the virii (not that this value is small). If each mutation must work in groups of 3 ACGT values, perhaps you could even store your string as a bitstring of some sort (4**3 == 64).

Death of a virus would be a -- operation, and a successful mutation would be a ++ operation.

There is something wrong with the way I am passing the arrays to the hash .. it will output ie ARRAY(0x1f512580), not sure where exactly the problem is.

           my @last_gen=$allgen{$i-1};
[download]

           my @last_gen=@{$allgen{$i-1}};
[download]

Thankyou dearly. I spent way too much time trying to figure that out! Also why is it that %hash=(0=>A) must be followed by $hash{1}=B ( for next key value pair you set I mean. ) rather than using the first format twice? I can't find this anywhere in the camel book or online.

I guess you are asking why you can't do this:

%hash=(0=>A);
%hash=(1=>B);
[download]

Both these operations operate on the whole hash, initializing it with the data on the right side of the equation. The equation sign "=" means literally "the left side is made equal to the right side", not "the right side is added to or combined with the left side". So after the second equation any previous values in the hash will be lost

While when you use $hash{1}=B you only manipulate one value inside the hash, leaving the rest of the hash unchanged.

Thankyou dearly. I spent way too much time trying to figure that out! Also why is it that %hash=(0=>A) must be followed by $hash{1}=B ( for next key value pair you set I mean. ) rather than using the first format twice? I can't find this anywhere in the camel book or online.

It does not have to be followed. This is covered in perlintro under Perl variable types, or if you're new to programming, http://learn.perl.org/books/beginning-perl/

IIUC< Why not use AoA, array of arrays? Maybe an idea can be had in Making little array babies

ARRAY(0x1f512580) is a pointer in memory referring to the contents held in ARRAY, so to access these contents you need to dereference the reference References quick reference and Referencing in advanced data structures and of course the documentation perlref and perlreftut

Not really cuz not every biology problem is a BioPerl problem and not every BioPerl individual module can only be applied solely to biology problems. Biology in IT requires other tools from mother Perl such as PDL and Math::* modules for example.

Read Perl and Bioinformatics