BACKGROUND
randomly shuffled sequences are routinely used in sequence analysis to evaluate the statistical significance of a biological sequence.
for example, a common method for assessing the thermodynamic stability of an rna sequence is to compare its folding free energy with those of a large sample of random sequences.
it is known that the stability of an rna secondary structure depends crucially on the stackings of adjacent base pairs; therefore the frequencies of distinct doublets in the random sequences are important considerations in such analysis  <cit> .
besides, natural biological sequences often manifest certain nearest-neighbor patterns: both eukaryotic and prokaryotic nucleic acid sequences show a consistent hierarchy in the doublet frequencies; in coding regions, the codon usage can also be markedly nonuniform.
in many cases, biologists need sophisticated shuffling tools that preserve not only the counts of distinct letters but also higher-order statistics such as doublet counts, triplet counts, and, in general, k-let counts.
methods for random sequence generation
several methods are commonly used to generate random sequences.
the basic permutation method works as follows: for a sequence s , pick a random number i between  <dig> and n, swap the two elements s  and s , then recurse on the subsequence s .
the random sequence generated by the basic permutation method preserves the exact count of each distinct letter in the alphabet, but does not preserve the higher-order statistics of k-let counts.
the markov method  <cit> , which is based on the markov chains, generates random sequences that preserve the k-let counts only on average: the counts of the individual sequences may deviate from the input distribution.
the swapping method  <cit> , a popular method which is now folklore, generates random sequences by repeatedly swapping disjoint subsequences flanked by the same -lets; it does preserve the k-let counts exactly, but produces random sequences that are only uniform asymptotically and may need a large number of swapping steps.
the euler algorithm preserves exact k-let counts
the euler algorithm is a less-known but very efficient algorithm for generating truly uniform random k-let-preserving sequences  <cit> .
we briefly review its history.
fitch  <cit>  first noticed that a doublet-preserving permutation is related to an eulerian walk of a directed multigraph; however, the algorithm he proposed does not generate all permutations with equal probability.
altschul and erickson  <cit>  presented the first algorithm  for generating truly uniform random sequences that preserve either the doublet counts or the triplet counts or both; however, a crucial step of their algorithm for generating random arborescences depends on a trial-and-error procedure, which is a potential bottleneck in performance.
this bottleneck was eliminated by kandel et al.
<cit> , who replaced the trial-and-error procedure with a simple and efficient procedure based on random walks in directed multigraphs.
they also generalized the euler algorithm to preserve the k-let counts for arbitrary k, and suggested a simple data structure for implementation.
this data structure is based on look-up tables and requires o space and time; it quickly become inefficient as the alphabet size σ and the let size k increase.
since the work by kandel et al.
<cit> , a better algorithm has been proposed by wilson  <cit>  for generating random arborescences, which is the crucial step of the euler algorithm that kandel et al.
<cit>  improves upon altschul and erickson  <cit> .
the superiority of wilson's arborescence generation algorithm to the two previous algorithms by altschul and erickson  <cit>  and by kandel et al.
<cit>  is both proved in the theoretical sense by wilson  <cit> , and demonstrated in the practical sense by a comparison of our implementation with a previous implementation .
implementations of the euler algorithm
we are aware of two previous implementations of earlier variants of the euler algorithm.
the dishuffle program by clote et al.
<cit>  implements the original version of the euler algorithm by altschul and erickson  <cit> .
the shufflet program by coward  <cit>  implements the improved version of the euler algorithm by kandel et al.
<cit> .
in this paper, we present a sequence analysis tool  for shuffling biological sequences while preserving the k-let counts.
the ushuffle program is based on the latest variant of the euler algorithm  <cit>  and uses wilson's algorithm  <cit>  in the crucial step of arborescence generation.
our goal is to provide a versatile tool that is as efficient and as flexible as possible:
arbitrary alphabet size and let size
in specific applications, the alphabet size σ and the let size k are often fixed: for biological sequences, typical alphabet sizes are  <dig>  and  <dig> , and typical let sizes are  <dig>  and  <dig> .
while it is tempting to implement the euler algorithm just for the fixed alphabet and let sizes at hand, we believe the flexibility of arbitrary alphabet and let sizes is useful.
the dishuffle program by clote et al.
<cit> , for example, is hard-coded for shuffling rna sequences preserving dinucleotide counts .
it is apparent that such an implementation cannot be used easily in other applications with different alphabet and let sizes.
efficiency
when the alphabet size and the let size are both small constants, the running time of the euler algorithm  is linear in the sequence length.
so it may appear that the efficiency of the shuffling program would not an issue since any conceivable downstream analysis of the randomized data would much slower than the shuffling.
however, we note that the linear running time has been proved only for the case that the alphabet and let sizes are constant  <cit> .
it is not at all clear whether the linear performance of the euler algorithm is scalable for arbitrary alphabet and let sizes.
as mentioned earlier, the "standard" data structure suggested by kandel et al.
<cit>  has time and space complexities o, which can become exponential when the alphabet size σ and the let size k become large, approaching the order of the sequence length.
indeed, as we will discuss later, we have reason to believe that this very data structure has been used in the shufflet program by coward  <cit> .
furthermore, the implementation of the euler algorithm  is non-trivial because of its heavy use of graph-theoretical concepts such as directed multigraphs and eulerian walks.
although wilson's celebrated algorithm  <cit>  dates back to  <dig>  and is well-known in the theoretical computer science community, coward's implementation of shufflet in  <dig>  <cit>  still uses the old arborescence algorithm by kandel et al.
<cit> .
we are not aware of any implementation of wilson's algorithm in bioinformatics applications.
by careful choices of algorithms and data structures, and by scrupulous algorithmic engineering, we strive for the most efficient implementation.
multiple forms and programming languages
the dishuffle program by clote et al.
<cit>  is written in python; the shufflet program by coward  <cit>  is a web application in c. to reach the widest audience, we have made our ushuffle program available in several forms.
it can be used as a command-line program, a web application, or a utility library.
source code in c, java, and c#, and integration instructions for perl and python are provided.
implementation
this section consists of four subsections.
in the first two subsections, we discuss, at a conceptual level, the euler algorithm and its crucial step of arborescence generation, in preparation of the discussion of implementation details.
in the third subsection, we present the algorithmic engineering details of our implementation.
in the fourth subsection, we describe the software organization and user interfaces of the ushuffle tool.
to justify our algorithm choices and to explain our optimization techniques, the discussions in the first three subsections are necessarily technical.
the readers who are not particularly interested in the theoretical discussion of graph algorithms or the technical details of algorithmic engineering can safely skip to the fourth subsection for the software organization and user interfaces.
the euler algorithm
in this subsection, we review some basic concepts of the euler algorithm.
directed multigraph
a k-let is a subsequence of k consecutive elements in a sequence.
let s be a sequence to be permuted.
let tk be a uniform random sequence that preserves the k-let counts of s.  to generate tk for k ≥  <dig>  the euler algorithm  <cit>  first constructs a directed multigraph g. we refer to figure  <dig> for an example.
for each distinct -let in s, g has a vertex.
for each k-let l in s, which contains two -lets l <dig> and l <dig> such that l <dig> precedes l <dig>  g has a directed edge from the vertex for l <dig> to the vertex for l <dig>  duplicates of k-lets may exist in s, so there may be multiple edges between the vertices.
correspondence between permutations and eulerian walks
as we scan the k-lets in s one by one, we also walk in the directed multigraph g from vertex to vertex.
when all the k-lets are scanned, each edge in g is visited exactly once: the walk is eulerian.
on the other hand, given an eulerian walk in g, we can recover a sequence by spelling out the -lets of the vertices along the walk .
since each k-let in s corresponds to an edge in g, every eulerian walk in g corresponds to a sequence with the same k-let counts as s. kandel et al.
<cit>  showed that, as long as an eulerian walk starts and ends at the same two vertices s and t that correspond to the starting and the ending -lets of s, the i-let counts for all  <dig> ≤ i ≤ k are preserved.
therefore, generating a uniform random sequence tk reduces to generating a uniform random eulerian walk in g from s to t.
correspondence between eulerian walks and arborescences
for an eulerian walk in g, each vertex v of g except the ending vertex t has a last edge ev that exits from v for the last time.
the set of last edges for all vertices except t forms an arborescence rooted at t: a directed spanning tree in which all vertices can reach t. given an arborescence a rooted at t, a random eulerian walk from s to t with the last edges conforming to a can be easily generated  <cit> :
<dig>  for each vertex v, collect the list of edges ev exiting from v. permute each edge list ev separately while keeping ev the last edge on the list.
<dig>  walk the graph g in accordance with the edge lists {ev}: start from s , take the first unmarked edge  from the list eu, mark the edge, then move to the next vertex v ; continue until all edges are marked and the walk ends at t.
in directed multigraphs, there is a nice correspondence between eulerian walks and arborescences: every arborescence rooted at t corresponds to exactly the same number of eulerian walks  <cit> .
therefore, generating a uniform random eulerian walk in g from s to t reduces to generating a uniform random arborescence in g rooted at t. in the next subsection, we discuss algorithms for generating random arborescences, some of which are based on, quite amusingly, random walks again.
generating random arborescences
in this subsection, we review the existing algorithms for arborescence generation, and explain our choice of wilson's algorithm  <cit> .
there are two major approaches to generating random arborescences and spanning trees: determinant algorithms and random-walk algorithms.
determinant algorithms
determinant algorithms are based on the matrix tree theorem .
for a graph g, the probability that a particular edge e appears in a uniform random spanning tree is the ratio of two numbers: the number of spanning trees that contain the edge e, and the total number of spanning trees.
the matrix tree theorem allows one to compute the exact number of spanning trees of a graph by evaluating the determinant of the combinatorial laplacian  of the graph.
a random spanning tree can be generated by repeatedly contracting or deleting edges according to their probabilities.
the first determinant algorithms were given by guénoche  <cit>  and kulkarni  <cit> : for a graph of n vertices and m edges, a random spanning tree can be generated in o time.
this running time was later improved to o  <cit> .
colbourn, myrvold, and neufeld  <cit>  simplified the o time algorithm and showed that the running time can be further reduced to o, the best upper bound for multiplying two n × n matrices  <cit> .
random-walk algorithms
random-walk algorithms use an entirely different approach to generating random spanning trees.
aldous  <cit>  and broder  <cit>   independently discovered an interesting connection between random spanning trees and random walks:
simulate a uniform random walk in a graph g starting at an arbitrary vertex s until all vertices are visited.
for each vertex v ≠ s, collect the edge {u, v} that corresponds to the first entrance to v. the collection t of edges is a uniform random spanning tree of g.
for a graph g and a vertex v in it, define the cover time cv as the expected number of steps a random walk starting from v takes to visit all vertices of g. the running time of the aldous-broder algorithm  <cit>  is clearly linear in the cover time.
in the context of shuffling biological sequences, kandel et al.
<cit>  extended aldous-broder algorithm  <cit>  to generate uniform random arborescences of eulerian directed graphs in the cover time.
wilson and propp  <cit>  then presented an algorithm for generating uniform random arborescences of general directed graphs in  <dig> cover times.
wilson's algorithm
wilson  <cit>  showed that random arborescences and spanning trees can be generated more quickly than the cover time by a cycle-popping algorithm which simulates loop-erased random walks.
for a graph g and two vertices u and v in it, define the hitting time hu,v as the expected number of steps a random walk takes from u to v. the running time of wilson's algorithm  <cit>  is linear in the maximum or mean hitting times of the corresponding stochastic graphs.
as wilson  <cit>  noted, the mean and maximum hitting times are always less than the cover time, and the differences can be quite significant in certain graphs.
therefore, for generating uniform random arborescences, wilson's algorithm  <cit>  is superior to kandel et al.
's algorithm  <cit> .
for completeness of presentation, we include in the following the pseudocode of wilson's algorithm  <cit> :
randomtreewithroot
1   for i ←  <dig> to n
2      intree  ← false
3   next  ← nil
4   intree  ← true
5   for i ←  <dig> to n
6      u ← i
7      while not intree
8         next  ← randomsuccessor
9          u ← next
10      u ← i
11      while not intree
12         intree  ← true
13         u ← next
14   return next
let eu be the set of directed edges exiting from the vertex u. the function randomsuccessor chooses a uniformly random edge  from eu, then returns the vertex v.
unlike the aldous-broder algorithm  <cit> , which simulates a single random walk from the root to visit all vertices, wilson's algorithm  <cit>  simulates multiple random walks: starting from each unvisited vertex, a random walk continues until it joins a growing arborescence which initially contains only the root.
a random walk follows the next  pointers; whenever a previously visited vertex is encountered again, a loop is formed and immediately erased because the next  pointer is overwritten .
as soon as a walk reaches the growing arborescence, all vertices in the walk join the arborescence as one more branch.
a comparison of the two approaches
we now give a comparison of the two approaches to generating random arborescences.
kandel et al.
<cit>  proved that the cover time of an eulerian directed multigraph of n vertices and m edges is o. from our preceding discussion on the cover time and the hitting time, it follows that the expected running time of wilson's algorithm  <cit>  on the same multigraph is at most o too, neglecting the log n factors.
for a multigraph, the number m of edges can be arbitrarily larger than the number n of vertices.
so it might appear that the determinant algorithm by colbourn et al.
<cit> , which runs in deterministic o time or even o time, would be a better alternative than the random walk algorithms  <cit> .
however, we note that when m is large the intermediate values of the determinant computation can be large too.
on the typical computer systems today, the arithmetic operations on floating-point numbers do not have enough precision to guarantee the accuracy and stability in the numerical computation of the determinant algorithms.
the random walk algorithms  <cit> , on the other hand, require only basic operations on small integers, and do not have these numerical problems.
therefore, we have decided to implement wilson's random-walk algorithm  <cit>  for arborescence generation.
implementation details
in this subsection, we describe the details of our implementation of the euler algorithm  <cit>  for generating k-let-preserving random sequences.
kandel et al.
's data structure
as suggested by kandel et al.
<cit> , a simple implementation of the euler algorithm  <cit>  may use a look-up table of size σk- <dig> for all possible -lets as vertices in the directed multigraph g, then build an adjacency matrix of size σk- <dig> × σk- <dig> for the edges in g. when both σ and k are small constants, the space requirement of this simple approach, σ2k- <dig>  may not look severe.
however, a calculation shows that, even for σ =  <dig>  and k =  <dig> , the space requirement amounts to
σ2k- <dig> =  <dig> =  <dig>   <dig>
on the other hand, the typical length of a protein sequence is below  <dig>  even though a sequence itself may be stored in only  <dig> kilo-bytes, the permutation algorithm still requires hundreds of times more space regardless.
the situation becomes even worse when k is further increased: even for the rather innocent-looking parameters σ =  <dig> and k =  <dig>  the space requirement
σ2k- <dig> =  <dig> >  <dig> =  <dig>
exceeds all  <dig> giga-bytes of memory that can be accommodated by a 32-bit computer!
we note that the two sets of parameters that coward  <cit>  used for experiments on his shufflet program were only
σ =  <dig>  k =  <dig>  σ2k- <dig> =  <dig>   <dig>   <dig>
and
σ =  <dig>  k =  <dig>  σ2k- <dig> =  <dig>   <dig>
we will discuss more about this in our comparison of ushuffle and shufflet in the results and discussion section.
representing directed multigraph in linear space
to make the ushuffle program scalable, it is clear that careful algorithmic engineering are necessary in the implementation.
as we discussed in the previous subsection on the euler algorithm, the directed multigraph g contains a vertex for each distinct -let in s. since the number of -lets in s is exactly l - k +  <dig>  g has at most l - k +  <dig> vertices, and hence exactly l - k +  <dig> directed edges between consecutive -lets.
this implies that the size of g is in fact linear in the length l of the sequence s to be permuted.
with suitable data structures, ushuffle needs only linear space.
in the following, we first explain the construction and representation of the directed multigraph g, then explain the random sequence generation after the graph construction.
the graph construction consists of two steps: determine the set of vertices, then add the directed edges.
determining vertices
we use a hashtable to determine the set of vertices.
the hashtable consists of a bucket array of size b = l - k +  <dig>  the number of  lets in s, and a linked list at each bucket to avoid collision by chaining  <cit> .
each -let x = x1x2⋯xk- <dig> has a polynomial hash code
h=x1ak−1+x2ak−2+⋯+xk−2a2+xk−1a=a+⋯+xk−2)a+xk−1)a,
where a = / <dig> is the reciprocal of the golden ratio; the index of x to the bucket array is
i = ⌊h·b⌋ mod b.
initialize the hashtable to be empty, then try to insert the -lets into the hashtable one by one.
if a -let is the first of its kind, it is assigned a new vertex number then inserted into the hashtable; its starting index to the sequence s is also recorded.
if a -let has been inserted before, it is not inserted to the hashtable: its vertex number and index to the sequence s are copied from those of the first -let of its kind.
after insertions, we can deduce the total number of vertices in the directed multigraph from the largest vertex number assigned.
the memory for vertices are then allocated.
adding directed edges
to add the directed edges, we use an adjacency-list representation to avoid the excessive memory requirement of an adjacency-matrix.
in an adjacency-list representation, two edge lists need to be maintained at each vertex: a list of incoming edges and a list of outgoing edges.
the outgoing edge lists are necessary for generating eulerian walks  <cit> .
the incoming edges lists are necessary for generating arborescences when kandel et al.
's algorithm  <cit>  is used .
we use wilson's algorithm  <cit>  for generating arborescences.
as we discussed in the previous section, wilson's algorithm  <cit>  is faster than kandel et al.
's algorithm  <cit> .
furthermore, we note here that wilson's algorithm  <cit>  has another advantage over kandel et al.
's algorithm  <cit>  in terms of the ease of implementation.
instead of one backward random walk from the ending vertex t to reach all other vertices as in kandel et al.
's algorithm  <cit> , wilson's algorithm  <cit>  uses multiple forward random walks from each unvisited vertex to join the arborescence rooted at t: the outgoing edge lists alone are sufficient for generating both the eulerian walks and the arborescences.
representing edge lists and managing memory
for maximum efficiency, we implement each edge list as an array of vertices.
the numbers of outgoing edges differ from vertex to vertex; if we allocate a fixed-size array for each vertex, then we would have to make each array large enough to hold all edges in the worst case, and the resulting space requirement would become quadratic in the length l of the sequence s. we could of course first count the number of outgoing edges for each vertex, then allocate a separate array just large enough for each vertex.
however, this would require us to call the relatively expensive memory allocation function once for each vertex.
in our implementation, we allocate one large array for all edges , then parcel out pieces to individual vertices.
to achieve this, we first scan the sequence s to count the number of outgoing edges for each vertex, then point the array  of each vertex to successive offsets of the large array.
with this optimization, the number of memory allocations is reduced to only 4: one for the hashtable bucket array, one for the array of -lets as hashtable entries, one of the array of vertices, and one for the array of edges.
the memory for the bucket array and the hashtable entries can be freed as soon as the directed multigraph is constructed.
sequence generation after graph construction
after the construction of the directed multigraph, we can generate a random sequence in three steps.
as discussed in the previous section, we need to first simulate the loop-erased random walks  <cit>  to generate an arborescence, next permute the individual edge lists while maintaining the last edges, then simulate an eulerian walk guided by the edge lists and output the sequence along the walk.
since each edge list is implemented as an array, the permutation can be executed very efficiently.
to output the random sequence along the walk is also easy, since each vertex keeps the starting index of its first occurrence in the input sequence.
software organization and user interfaces of the ushuffle tool
in this subsection, we describe the software organization and user interfaces of the ushuffle tool.
c library and command-line tool
our initial implementation of ushuffle is in the c programming language.
the c version of ushuffle consists of two components: a ushuffle library  and a command-line tool .
in a typical scenario, multiple k-let-preserving random sequences are generated for each input sequence.
the graph construction stage of the ushuffle program needs to be done only once for the multiple output sequences.
to give the users an option for optimization, we export three interface functions in the ushuffle library:
void shuffle;
void shuffle1;
void shuffle2;
the function shuffle accepts four parameters: s is the sequence to be permuted, t is the output random sequence, l is the length of s, and k is the let size k. the function shuffle simply calls shuffle <dig> first and shuffle <dig> next: shuffle <dig> implements the construction of the directed multigraph; shuffle <dig> implements the loop-erased random walks in the directed multigraph and the generation of the random sequence.
the statistical behavior of a random permutation depends heavily on the random number generator.
coward  <cit>  noted that the default implementations of random number generators on various platforms are often unsatisfying, so he implemented his own generator using an arguably better algorithm.
we note that there are numerous algorithms for random number generation, and new algorithms are continuously being proposed: whether one algorithm is superior to the other can be quite subjective.
instead of limiting the users to a particular implementation, we set the default generator to the random function from the standard c library, then export an interface function to allow sophisticated users to customize the generator:
typedef long ();
void set_randfunc;
the command-line ushuffle tool is a minimal front-end of the ushuffle library that illustrates a typical use of the library.
it has the following four options:
-s <string> specifies the input sequence,
-n <number> specifies the number of random sequences to generate,
-k <number> specifies the let size,
-seed <number> specifies the seed for random number generator.
java applet
the ushuffle program is ported to the java programming language.
beside having a library and a command-line tool, the java version of the ushuffle program can also run as an applet in a web browser.
we refer to figure  <dig> for a screenshot of the ushuffle java applet1: the interface of the applet is minimal and consists of three parts: an input text area at the top, an output text area at the bottom, and a control panel in the middle.
the control panel contains two text fields and a button.
the maximum let size k and the number n of output sequences can be set in the two text fields.
when the "shuffle" button is clicked, the applet takes the input sequence from the input text field, strips away the white spaces, generates n random sequences that preserve the k-let counts, then outputs the sequences in the output text area.
the output is in the fasta format when n > 1: each output sequence is preceded by a comment line containing a sequence number ranging from  <dig> to n.
the ushuffle java applet keeps all the output sequences in memory for display in the output text area.
when the number n of output sequences and the input sequence length l are exorbitantly large, for example, n =  <dig>   <dig>   <dig> and l =  <dig>  the total memory required to hold the output sequences may exceed the maximum heap size of the java virtual machine  and the applet may hang.
this is not a bug in our program but is due to the limit of jvm; nevertheless, we prepared a web page to instruct the users how to increase the maximum heap size of jvm.
c#/perl/python versions
the ushuffle program is also ported to the c# programming language.
perl and python are popular programming languages for bioinformatics; they allow easy integration with programs written in c. instead of porting the ushuffle program to perl and python at the source code level, we prepared two web pages to instruct the users how to extend the perl and python environments with the ushuffle library.
RESULTS
we have performed two sets of experiments to test the performance of two major forms of the ushuffle tool: we first benchmark the performance of the ushuffle c library, then compare the performance of the ushuffle java applet with the shufflet program by coward  <cit> .
performance of ushuffle c library
we tested the ushuffle c library on a desktop pc <dig> with test data consisting of both real biological sequences and artificially generated random sequences.
experiment on real biological sequences
the real biological sequences were acquired from two sources: first,  <dig> protein sequences  were sampled from the human protein reference database <dig>  one sequence from each of the  <dig> molecular classes; second,  <dig> micro rna precursor sequences  of mus.musculus  were extracted from the supplementary data <dig> of bonnet et al.
<cit> .
our experiments on these real biological sequences showed that the ushuffle library is extremely efficient: in just one second, it can generate either   <dig> doublet-preserving random sequences for each of  <dig> protein sequences, or   <dig> doublet-preserving random sequences for each of the  <dig> rna sequences.
experiments on artificially generated random sequences
in order to analyze the performance of ushuffle with various sets of parameters, we also performed a systematic test of ushuffle on artificially generated random sequences.
for simplicity, the sequence lengths were exact powers of two from  <dig> to  <dig>  that is, from around  <dig>   <dig> to around  <dig>   <dig>   <dig>  these numbers are somewhat arbitrary; nothing prevents a user from running ushuffle on very long sequences, even at the genome scale, as long as the computer has enough memory to store the input sequence and has some additional  memory required by our implementation.
for each sequence length,  <dig> uniform random sequences over the english alphabet  were generated as test sequences; for each test sequence,  <dig> k-let-preserving random sequences were then generated by ushuffle.
the total running time for ushuffle to generate the  <dig> ×  <dig> =  <dig> k-let-preserving random sequences was recorded for each sequence length.
two getrusage system calls were placed in the test program to sandwich the code region being benchmarked; the differences of the two timestamps were used to calculate the running times.
we refer to figure  <dig> for a log-log plot of the total running times of the ushuffle program for k =  <dig> and k =  <dig> at various sequence lengths.
the plot shows that the running time of the ushuffle program is essentially linear in the length of the sequence to be shuffled.
the absolute running times are not very effective in demonstrating the extreme efficiency of the ushuffle program.
we refer to figure  <dig> for a ratio plot that is more illustrative.
for k =  <dig> and k =  <dig>  and for each sequence length, the plot shows not the absolute running time of the ushuffle program but the ratio of two running times:
<dig>  the running time for the ushuffle program to generate the k-let-preserving sequences, and
<dig>  the running time for the simple permutation method  <cit>   to shuffle the same number of random sequences without preserving the k-let counts.
the ratio plot shows that the running time of the ushuffle program is on average only  <dig>  times that of the simple permutation method for k =  <dig>  and only  <dig> times for k =  <dig>  the simple permutation method is minimal: for each position of the input sequence it executes only one random function call plus one swap.
the ushuffle program, on the other hand, performs a lot more work; although the  <dig> k-let-preserving random sequences of each test sequence are generated by one shuffle <dig> and  <dig> shuffle <dig> function calls to avoid redundant multigraph construction, each shuffle <dig> function call still includes the generation of an arborescence by loop-erased random walks  <cit>  and the generation of an eulerian walk guided by the individual edge lists shuffled by simple permutations.
in light of the contrasting complexities of the ushuffle program and the simple permutation method, the small ratios of their running times are remarkable.
a careful reader will notice an interesting fact from figure  <dig>  when the sequence length increases to  <dig> , the running time of the ushuffle program for k =  <dig> is even less than the simple permutation method!
the "strange" phenomenon had kept us puzzling for a long time until we eventually convinced ourselves that this is not a bug but a feature.
we note that, in each step, the simple permutation method randomly swaps two elements scattered in a large array of  <dig> elements.
on the other hand, the ushuffle program performs random walks in small multigraphs  and permutes the individual edge lists  separately.
the memory references of the ushuffle program are much more local than those of the simple permutation method.
computers with modern memory architectures aggressively optimize code with local memory references by sophisticated caching schemes, which promotes the performance of the ushuffle program.
we refer to figure  <dig> for the running times of the ushuffle program at various values of the parameter k, where the test sequence length is fixed at  <dig>  the running time of the ushuffle program peaks at k =  <dig>  which is about three times its running time for k =  <dig>  then gradually decreases as k increases, and finally drops to zero at k =  <dig> because, with a sequence length of  <dig>  the only 1024-let-preserving random sequence is the input sequence itself.
this plot shows that the ushuffle program is efficient for all possible values of k.
comparison of ushuffle java applet with shufflet
there exist two other implementations of the euler algorithm.
the dishuffle program by clote et al.
<cit>  implements the original version of the euler algorithm by altschul and erickson  <cit> .
hard-coded for shuffling rna sequences preserving dinucleotide counts, dishuffle is not a general tool for arbitrary alphabet and let sizes.
another program, shufflet by coward  <cit> , implements the improved version of the euler algorithm by kandel et al.
<cit>  for arbitrary let size k. as we have explained in the implementation section, the arborescence generation algorithm by kandel et al.
<cit> , while superior to the algorithm by altschul and erickson  <cit> , is still inferior to wilson's algorithm  <cit> ; besides, its look-up table data structure is inefficient for large alphabet and let sizes.
in terms of functionality, the shufflet implementation  <cit>  is closer to our ushuffle implementation.
shufflet was written in the c programming language, and had been hosted as a web application .
we were unable to perform a comprehensive comparison of ushuffle and shufflet.
however, coward  <cit>  mentioned two experiments performed on a digital dec/alpha  <dig> web server:
<dig>   <dig> shufflings of a dna sequence of  <dig> nucleotides with k =  <dig> take about  <dig>  seconds;
<dig>   <dig> shufflings of a protein sequence of  <dig> amino acids with k =  <dig> take less than  <dig> second.
we performed similar experiments with the ushuffle java applet on an apple imac computer :
<dig>   <dig> shufflings of a dna sequence of  <dig> nucleotides with k =  <dig> take about  <dig>  seconds;
<dig>   <dig> shufflings of a protein sequence of  <dig> amino acids with k =  <dig> take less than  <dig> second.
assuming comparable performances of the two computers, we estimate that our ushuffle java applet is about 15– <dig> times faster than shufflet in the experiment on nucleotides , and about  <dig> times faster than shufflet in the experiment on amino acids .
we certainly understand the difficulty of such a comparison: a web server in  <dig> versus a desktop computer in 2005; a c program in a native machine versus a java applet in a virtual machine.
nevertheless the comparison illustrates the better scalability of our ushuffle java applet for large let sizes.
the difference between the two performance ratios, 15– <dig> versus  <dig>  suggests that ushuffle remains efficient even for large let size, while shufflet becomes more inefficient as the let size increases, due to  the use of the inefficient look-up table data structure by kandel et al.
<cit> .
CONCLUSIONS
the ushuffle tool is based on superior graph algorithms and is carefully engineered to be extremely efficient.
it achieves maximum flexibility by allowing arbitrary alphabet size and let size, and is available in many forms for different kinds of users.
we believe ushuffle is a useful tool for the bioinformatics community.
availability and requirements
project name: ushuffle.
project home page:
operating systems: platform independent.
programming languages: c, java, c#, perl, python.
other requirements: none.
licence: freebsd.
any restrictions to use by non-academics: none.
authors' contributions
mj designed the software and the experiments, implemented the c and java versions of ushuffle program, and wrote the technical report.
ja ported the java program to c#, investigated software licenses, and designed the ushuffle logo.
jg wrote the test program, performed the experiments, and wrote the instructions for perl integration.
mm wrote the java applet interface and the instructions for python integration.
all authors reviewed the source code and contributed to the home page construction.
note
1macos  <dig> .
<dig>  firefox  <dig> .
<dig> , java  <dig> .
<dig>
2dell xps m1710:  <dig>  ghz intel dual-core  <dig> processor,  <dig> gb dual in-line ram; microsoft windows vista business edition, cygwin  <dig> .
<dig>  gcc  <dig> .
<dig> with -o <dig> option.
<dig>
<dig>

 randomly shuffled sequences are routinely used in sequence analysis to evaluate the statistical significance of a biological sequence. for example, a common method for assessing the thermodynamic stability of an rna sequence is to compare its folding free energy with those of a large sample of random sequences. it is known that the stability of an rna secondary structure depends crucially on the stackings of adjacent base pairs; therefore the frequencies of distinct doublets in the random sequences are important considerations in such analysis. besides, natural biological sequences often manifest certain nearest-neighbor patterns: both eukaryotic and prokaryotic nucleic acid sequences show a consistent hierarchy in the doublet frequencies; in coding regions, the codon usage can also be markedly nonuniform. in many cases, biologists need sophisticated shuffling tools that preserve not only the counts of distinct letters but also higher-order statistics such as doublet counts, triplet counts, and, in general, k-let counts. methods for random sequence generation several methods are commonly used to generate random sequences. the basic permutation method works as follows: for a sequence s, pick a random number i between and n, swap the two elements s and s, then recurse on the subsequence s. the random sequence generated by the basic permutation method preserves the exact count of each distinct letter in the alphabet, but does not preserve the higher-order statistics of k-let counts. the markov method, which is based on the markov chains, generates random sequences that preserve the k-let counts only on average: the counts of the individual sequences may deviate from the input distribution. the swapping method, a popular method which is now folklore, generates random sequences by repeatedly swapping disjoint subsequences flanked by the same -lets; it does preserve the k-let counts exactly, but produces random sequences that are only uniform asymptotically and may need a large number of swapping steps. the euler algorithm preserves exact k-let counts the euler algorithm is a less-known but very efficient algorithm for generating truly uniform random k-let-preserving sequences. we briefly review its history. fitch first noticed that a doublet-preserving permutation is related to an eulerian walk of a directed multigraph; however, the algorithm he proposed does not generate all permutations with equal probability. altschul and erickson presented the first algorithm for generating truly uniform random sequences that preserve either the doublet counts or the triplet counts or both; however, a crucial step of their algorithm for generating random arborescences depends on a trial-and-error procedure, which is a potential bottleneck in performance. this bottleneck was eliminated by kandel et al., who replaced the trial-and-error procedure with a simple and efficient procedure based on random walks in directed multigraphs. they also generalized the euler algorithm to preserve the k-let counts for arbitrary k, and suggested a simple data structure for implementation. this data structure is based on look-up tables and requires o space and time; it quickly become inefficient as the alphabet size σ and the let size k increase. since the work by kandel et al., a better algorithm has been proposed by wilson for generating random arborescences, which is the crucial step of the euler algorithm that kandel et al. improves upon altschul and erickson. the superiority of wilson's arborescence generation algorithm to the two previous algorithms by altschul and erickson and by kandel et al. is both proved in the theoretical sense by wilson, and demonstrated in the practical sense by a comparison of our implementation with a previous implementation. implementations of the euler algorithm we are aware of two previous implementations of earlier variants of the euler algorithm. the dishuffle program by clote et al. implements the original version of the euler algorithm by altschul and erickson. the shufflet program by coward implements the improved version of the euler algorithm by kandel et al.. in this paper, we present a sequence analysis tool for shuffling biological sequences while preserving the k-let counts. the ushuffle program is based on the latest variant of the euler algorithm and uses wilson's algorithm in the crucial step of arborescence generation. our goal is to provide a versatile tool that is as efficient and as flexible as possible: arbitrary alphabet size and let size in specific applications, the alphabet size σ and the let size k are often fixed: for biological sequences, typical alphabet sizes are and, and typical let sizes are and. while it is tempting to implement the euler algorithm just for the fixed alphabet and let sizes at hand, we believe the flexibility of arbitrary alphabet and let sizes is useful. the dishuffle program by clote et al., for example, is hard-coded for shuffling rna sequences preserving dinucleotide counts. it is apparent that such an implementation cannot be used easily in other applications with different alphabet and let sizes. efficiency when the alphabet size and the let size are both small constants, the running time of the euler algorithm is linear in the sequence length. so it may appear that the efficiency of the shuffling program would not an issue since any conceivable downstream analysis of the randomized data would much slower than the shuffling. however, we note that the linear running time has been proved only for the case that the alphabet and let sizes are constant. it is not at all clear whether the linear performance of the euler algorithm is scalable for arbitrary alphabet and let sizes. as mentioned earlier, the "standard" data structure suggested by kandel et al. has time and space complexities o, which can become exponential when the alphabet size σ and the let size k become large, approaching the order of the sequence length. indeed, as we will discuss later, we have reason to believe that this very data structure has been used in the shufflet program by coward. furthermore, the implementation of the euler algorithm is non-trivial because of its heavy use of graph-theoretical concepts such as directed multigraphs and eulerian walks. although wilson's celebrated algorithm dates back to and is well-known in the theoretical computer science community, coward's implementation of shufflet in still uses the old arborescence algorithm by kandel et al.. we are not aware of any implementation of wilson's algorithm in bioinformatics applications. by careful choices of algorithms and data structures, and by scrupulous algorithmic engineering, we strive for the most efficient implementation. multiple forms and programming languages the dishuffle program by clote et al. is written in python; the shufflet program by coward is a web application in c. to reach the widest audience, we have made our ushuffle program available in several forms. it can be used as a command-line program, a web application, or a utility library. source code in c, java, and c#, and integration instructions for perl and python are provided. implementation this section consists of four subsections. in the first two subsections, we discuss, at a conceptual level, the euler algorithm and its crucial step of arborescence generation, in preparation of the discussion of implementation details. in the third subsection, we present the algorithmic engineering details of our implementation. in the fourth subsection, we describe the software organization and user interfaces of the ushuffle tool. to justify our algorithm choices and to explain our optimization techniques, the discussions in the first three subsections are necessarily technical. the readers who are not particularly interested in the theoretical discussion of graph algorithms or the technical details of algorithmic engineering can safely skip to the fourth subsection for the software organization and user interfaces. the euler algorithm in this subsection, we review some basic concepts of the euler algorithm. directed multigraph a k-let is a subsequence of k consecutive elements in a sequence. let s be a sequence to be permuted. let tk be a uniform random sequence that preserves the k-let counts of s. to generate tk for k ≥ the euler algorithm first constructs a directed multigraph g. we refer to figure for an example. for each distinct -let in s, g has a vertex. for each k-let l in s, which contains two -lets l and l such that l precedes l g has a directed edge from the vertex for l to the vertex for l duplicates of k-lets may exist in s, so there may be multiple edges between the vertices. correspondence between permutations and eulerian walks as we scan the k-lets in s one by one, we also walk in the directed multigraph g from vertex to vertex. when all the k-lets are scanned, each edge in g is visited exactly once: the walk is eulerian. on the other hand, given an eulerian walk in g, we can recover a sequence by spelling out the -lets of the vertices along the walk. since each k-let in s corresponds to an edge in g, every eulerian walk in g corresponds to a sequence with the same k-let counts as s. kandel et al. showed that, as long as an eulerian walk starts and ends at the same two vertices s and t that correspond to the starting and the ending -lets of s, the i-let counts for all ≤ i ≤ k are preserved. therefore, generating a uniform random sequence tk reduces to generating a uniform random eulerian walk in g from s to t. correspondence between eulerian walks and arborescences for an eulerian walk in g, each vertex v of g except the ending vertex t has a last edge ev that exits from v for the last time. the set of last edges for all vertices except t forms an arborescence rooted at t: a directed spanning tree in which all vertices can reach t. given an arborescence a rooted at t, a random eulerian walk from s to t with the last edges conforming to a can be easily generated: for each vertex v, collect the list of edges ev exiting from v. permute each edge list ev separately while keeping ev the last edge on the list. walk the graph g in accordance with the edge lists {ev}: start from s, take the first unmarked edge from the list eu, mark the edge, then move to the next vertex v; continue until all edges are marked and the walk ends at t. in directed multigraphs, there is a nice correspondence between eulerian walks and arborescences: every arborescence rooted at t corresponds to exactly the same number of eulerian walks. therefore, generating a uniform random eulerian walk in g from s to t reduces to generating a uniform random arborescence in g rooted at t. in the next subsection, we discuss algorithms for generating random arborescences, some of which are based on, quite amusingly, random walks again. generating random arborescences in this subsection, we review the existing algorithms for arborescence generation, and explain our choice of wilson's algorithm. there are two major approaches to generating random arborescences and spanning trees: determinant algorithms and random-walk algorithms. determinant algorithms determinant algorithms are based on the matrix tree theorem. for a graph g, the probability that a particular edge e appears in a uniform random spanning tree is the ratio of two numbers: the number of spanning trees that contain the edge e, and the total number of spanning trees. the matrix tree theorem allows one to compute the exact number of spanning trees of a graph by evaluating the determinant of the combinatorial laplacian of the graph. a random spanning tree can be generated by repeatedly contracting or deleting edges according to their probabilities. the first determinant algorithms were given by guénoche and kulkarni: for a graph of n vertices and m edges, a random spanning tree can be generated in o time. this running time was later improved to o. colbourn, myrvold, and neufeld simplified the o time algorithm and showed that the running time can be further reduced to o, the best upper bound for multiplying two n × n matrices. random-walk algorithms random-walk algorithms use an entirely different approach to generating random spanning trees. aldous and broder independently discovered an interesting connection between random spanning trees and random walks: simulate a uniform random walk in a graph g starting at an arbitrary vertex s until all vertices are visited. for each vertex v ≠ s, collect the edge {u, v} that corresponds to the first entrance to v. the collection t of edges is a uniform random spanning tree of g. for a graph g and a vertex v in it, define the cover time cv as the expected number of steps a random walk starting from v takes to visit all vertices of g. the running time of the aldous-broder algorithm is clearly linear in the cover time. in the context of shuffling biological sequences, kandel et al. extended aldous-broder algorithm to generate uniform random arborescences of eulerian directed graphs in the cover time. wilson and propp then presented an algorithm for generating uniform random arborescences of general directed graphs in cover times. wilson's algorithm wilson showed that random arborescences and spanning trees can be generated more quickly than the cover time by a cycle-popping algorithm which simulates loop-erased random walks. for a graph g and two vertices u and v in it, define the hitting time hu,v as the expected number of steps a random walk takes from u to v. the running time of wilson's algorithm is linear in the maximum or mean hitting times of the corresponding stochastic graphs. as wilson noted, the mean and maximum hitting times are always less than the cover time, and the differences can be quite significant in certain graphs. therefore, for generating uniform random arborescences, wilson's algorithm is superior to kandel et al.'s algorithm. for completeness of presentation, we include in the following the pseudocode of wilson's algorithm: randomtreewithroot 1   for i ← to n 2      intree ← false 3   next ← nil 4   intree ← true 5   for i ← to n 6      u ← i 7      while not intree 8         next ← randomsuccessor 9          u ← next 10      u ← i 11      while not intree 12         intree ← true 13         u ← next 14   return next let eu be the set of directed edges exiting from the vertex u. the function randomsuccessor chooses a uniformly random edge from eu, then returns the vertex v. unlike the aldous-broder algorithm, which simulates a single random walk from the root to visit all vertices, wilson's algorithm simulates multiple random walks: starting from each unvisited vertex, a random walk continues until it joins a growing arborescence which initially contains only the root. a random walk follows the next pointers; whenever a previously visited vertex is encountered again, a loop is formed and immediately erased because the next pointer is overwritten. as soon as a walk reaches the growing arborescence, all vertices in the walk join the arborescence as one more branch. a comparison of the two approaches we now give a comparison of the two approaches to generating random arborescences. kandel et al. proved that the cover time of an eulerian directed multigraph of n vertices and m edges is o. from our preceding discussion on the cover time and the hitting time, it follows that the expected running time of wilson's algorithm on the same multigraph is at most o too, neglecting the log n factors. for a multigraph, the number m of edges can be arbitrarily larger than the number n of vertices. so it might appear that the determinant algorithm by colbourn et al., which runs in deterministic o time or even o time, would be a better alternative than the random walk algorithms. however, we note that when m is large the intermediate values of the determinant computation can be large too. on the typical computer systems today, the arithmetic operations on floating-point numbers do not have enough precision to guarantee the accuracy and stability in the numerical computation of the determinant algorithms. the random walk algorithms, on the other hand, require only basic operations on small integers, and do not have these numerical problems. therefore, we have decided to implement wilson's random-walk algorithm for arborescence generation. implementation details in this subsection, we describe the details of our implementation of the euler algorithm for generating k-let-preserving random sequences. kandel et al.'s data structure as suggested by kandel et al., a simple implementation of the euler algorithm may use a look-up table of size σk- for all possible -lets as vertices in the directed multigraph g, then build an adjacency matrix of size σk- × σk- for the edges in g. when both σ and k are small constants, the space requirement of this simple approach, σ2k- may not look severe. however, a calculation shows that, even for σ = and k =, the space requirement amounts to σ2k- = = on the other hand, the typical length of a protein sequence is below even though a sequence itself may be stored in only kilo-bytes, the permutation algorithm still requires hundreds of times more space regardless. the situation becomes even worse when k is further increased: even for the rather innocent-looking parameters σ = and k = the space requirement σ2k- = > = exceeds all giga-bytes of memory that can be accommodated by a 32-bit computer! we note that the two sets of parameters that coward used for experiments on his shufflet program were only σ = k = σ2k- = and σ = k = σ2k- = we will discuss more about this in our comparison of ushuffle and shufflet in the results and discussion section. representing directed multigraph in linear space to make the ushuffle program scalable, it is clear that careful algorithmic engineering are necessary in the implementation. as we discussed in the previous subsection on the euler algorithm, the directed multigraph g contains a vertex for each distinct -let in s. since the number of -lets in s is exactly l - k + g has at most l - k + vertices, and hence exactly l - k + directed edges between consecutive -lets. this implies that the size of g is in fact linear in the length l of the sequence s to be permuted. with suitable data structures, ushuffle needs only linear space. in the following, we first explain the construction and representation of the directed multigraph g, then explain the random sequence generation after the graph construction. the graph construction consists of two steps: determine the set of vertices, then add the directed edges. determining vertices we use a hashtable to determine the set of vertices. the hashtable consists of a bucket array of size b = l - k + the number of lets in s, and a linked list at each bucket to avoid collision by chaining. each -let x = x1x2⋯xk- has a polynomial hash code h=x1ak−1+x2ak−2+⋯+xk−2a2+xk−1a=a+⋯+xk−2)a+xk−1)a, where a = / is the reciprocal of the golden ratio; the index of x to the bucket array is i = ⌊h·b⌋ mod b. initialize the hashtable to be empty, then try to insert the -lets into the hashtable one by one. if a -let is the first of its kind, it is assigned a new vertex number then inserted into the hashtable; its starting index to the sequence s is also recorded. if a -let has been inserted before, it is not inserted to the hashtable: its vertex number and index to the sequence s are copied from those of the first -let of its kind. after insertions, we can deduce the total number of vertices in the directed multigraph from the largest vertex number assigned. the memory for vertices are then allocated. adding directed edges to add the directed edges, we use an adjacency-list representation to avoid the excessive memory requirement of an adjacency-matrix. in an adjacency-list representation, two edge lists need to be maintained at each vertex: a list of incoming edges and a list of outgoing edges. the outgoing edge lists are necessary for generating eulerian walks. the incoming edges lists are necessary for generating arborescences when kandel et al.'s algorithm is used. we use wilson's algorithm for generating arborescences. as we discussed in the previous section, wilson's algorithm is faster than kandel et al.'s algorithm. furthermore, we note here that wilson's algorithm has another advantage over kandel et al.'s algorithm in terms of the ease of implementation. instead of one backward random walk from the ending vertex t to reach all other vertices as in kandel et al.'s algorithm, wilson's algorithm uses multiple forward random walks from each unvisited vertex to join the arborescence rooted at t: the outgoing edge lists alone are sufficient for generating both the eulerian walks and the arborescences. representing edge lists and managing memory for maximum efficiency, we implement each edge list as an array of vertices. the numbers of outgoing edges differ from vertex to vertex; if we allocate a fixed-size array for each vertex, then we would have to make each array large enough to hold all edges in the worst case, and the resulting space requirement would become quadratic in the length l of the sequence s. we could of course first count the number of outgoing edges for each vertex, then allocate a separate array just large enough for each vertex. however, this would require us to call the relatively expensive memory allocation function once for each vertex. in our implementation, we allocate one large array for all edges, then parcel out pieces to individual vertices. to achieve this, we first scan the sequence s to count the number of outgoing edges for each vertex, then point the array of each vertex to successive offsets of the large array. with this optimization, the number of memory allocations is reduced to only 4: one for the hashtable bucket array, one for the array of -lets as hashtable entries, one of the array of vertices, and one for the array of edges. the memory for the bucket array and the hashtable entries can be freed as soon as the directed multigraph is constructed. sequence generation after graph construction after the construction of the directed multigraph, we can generate a random sequence in three steps. as discussed in the previous section, we need to first simulate the loop-erased random walks to generate an arborescence, next permute the individual edge lists while maintaining the last edges, then simulate an eulerian walk guided by the edge lists and output the sequence along the walk. since each edge list is implemented as an array, the permutation can be executed very efficiently. to output the random sequence along the walk is also easy, since each vertex keeps the starting index of its first occurrence in the input sequence. software organization and user interfaces of the ushuffle tool in this subsection, we describe the software organization and user interfaces of the ushuffle tool. c library and command-line tool our initial implementation of ushuffle is in the c programming language. the c version of ushuffle consists of two components: a ushuffle library and a command-line tool. in a typical scenario, multiple k-let-preserving random sequences are generated for each input sequence. the graph construction stage of the ushuffle program needs to be done only once for the multiple output sequences. to give the users an option for optimization, we export three interface functions in the ushuffle library: void shuffle; void shuffle1; void shuffle2; the function shuffle accepts four parameters: s is the sequence to be permuted, t is the output random sequence, l is the length of s, and k is the let size k. the function shuffle simply calls shuffle first and shuffle next: shuffle implements the construction of the directed multigraph; shuffle implements the loop-erased random walks in the directed multigraph and the generation of the random sequence. the statistical behavior of a random permutation depends heavily on the random number generator. coward noted that the default implementations of random number generators on various platforms are often unsatisfying, so he implemented his own generator using an arguably better algorithm. we note that there are numerous algorithms for random number generation, and new algorithms are continuously being proposed: whether one algorithm is superior to the other can be quite subjective. instead of limiting the users to a particular implementation, we set the default generator to the random function from the standard c library, then export an interface function to allow sophisticated users to customize the generator: typedef long (); void set_randfunc; the command-line ushuffle tool is a minimal front-end of the ushuffle library that illustrates a typical use of the library. it has the following four options: -s <string> specifies the input sequence, -n <number> specifies the number of random sequences to generate, -k <number> specifies the let size, -seed <number> specifies the seed for random number generator. java applet the ushuffle program is ported to the java programming language. beside having a library and a command-line tool, the java version of the ushuffle program can also run as an applet in a web browser. we refer to figure for a screenshot of the ushuffle java applet1: the interface of the applet is minimal and consists of three parts: an input text area at the top, an output text area at the bottom, and a control panel in the middle. the control panel contains two text fields and a button. the maximum let size k and the number n of output sequences can be set in the two text fields. when the "shuffle" button is clicked, the applet takes the input sequence from the input text field, strips away the white spaces, generates n random sequences that preserve the k-let counts, then outputs the sequences in the output text area. the output is in the fasta format when n > 1: each output sequence is preceded by a comment line containing a sequence number ranging from to n. the ushuffle java applet keeps all the output sequences in memory for display in the output text area. when the number n of output sequences and the input sequence length l are exorbitantly large, for example, n = and l = the total memory required to hold the output sequences may exceed the maximum heap size of the java virtual machine and the applet may hang. this is not a bug in our program but is due to the limit of jvm; nevertheless, we prepared a web page to instruct the users how to increase the maximum heap size of jvm. c#/perl/python versions the ushuffle program is also ported to the c# programming language. perl and python are popular programming languages for bioinformatics; they allow easy integration with programs written in c. instead of porting the ushuffle program to perl and python at the source code level, we prepared two web pages to instruct the users how to extend the perl and python environments with the ushuffle library.  we have performed two sets of experiments to test the performance of two major forms of the ushuffle tool: we first benchmark the performance of the ushuffle c library, then compare the performance of the ushuffle java applet with the shufflet program by coward. performance of ushuffle c library we tested the ushuffle c library on a desktop pc with test data consisting of both real biological sequences and artificially generated random sequences. experiment on real biological sequences the real biological sequences were acquired from two sources: first, protein sequences were sampled from the human protein reference database one sequence from each of the molecular classes; second, micro rna precursor sequences of mus.musculus were extracted from the supplementary data of bonnet et al.. our experiments on these real biological sequences showed that the ushuffle library is extremely efficient: in just one second, it can generate either doublet-preserving random sequences for each of protein sequences, or doublet-preserving random sequences for each of the rna sequences. experiments on artificially generated random sequences in order to analyze the performance of ushuffle with various sets of parameters, we also performed a systematic test of ushuffle on artificially generated random sequences. for simplicity, the sequence lengths were exact powers of two from to that is, from around to around these numbers are somewhat arbitrary; nothing prevents a user from running ushuffle on very long sequences, even at the genome scale, as long as the computer has enough memory to store the input sequence and has some additional memory required by our implementation. for each sequence length, uniform random sequences over the english alphabet were generated as test sequences; for each test sequence, k-let-preserving random sequences were then generated by ushuffle. the total running time for ushuffle to generate the × = k-let-preserving random sequences was recorded for each sequence length. two getrusage system calls were placed in the test program to sandwich the code region being benchmarked; the differences of the two timestamps were used to calculate the running times. we refer to figure for a log-log plot of the total running times of the ushuffle program for k = and k = at various sequence lengths. the plot shows that the running time of the ushuffle program is essentially linear in the length of the sequence to be shuffled. the absolute running times are not very effective in demonstrating the extreme efficiency of the ushuffle program. we refer to figure for a ratio plot that is more illustrative. for k = and k = and for each sequence length, the plot shows not the absolute running time of the ushuffle program but the ratio of two running times: the running time for the ushuffle program to generate the k-let-preserving sequences, and the running time for the simple permutation method to shuffle the same number of random sequences without preserving the k-let counts. the ratio plot shows that the running time of the ushuffle program is on average only times that of the simple permutation method for k = and only times for k = the simple permutation method is minimal: for each position of the input sequence it executes only one random function call plus one swap. the ushuffle program, on the other hand, performs a lot more work; although the k-let-preserving random sequences of each test sequence are generated by one shuffle and shuffle function calls to avoid redundant multigraph construction, each shuffle function call still includes the generation of an arborescence by loop-erased random walks and the generation of an eulerian walk guided by the individual edge lists shuffled by simple permutations. in light of the contrasting complexities of the ushuffle program and the simple permutation method, the small ratios of their running times are remarkable. a careful reader will notice an interesting fact from figure when the sequence length increases to, the running time of the ushuffle program for k = is even less than the simple permutation method! the "strange" phenomenon had kept us puzzling for a long time until we eventually convinced ourselves that this is not a bug but a feature. we note that, in each step, the simple permutation method randomly swaps two elements scattered in a large array of elements. on the other hand, the ushuffle program performs random walks in small multigraphs and permutes the individual edge lists separately. the memory references of the ushuffle program are much more local than those of the simple permutation method. computers with modern memory architectures aggressively optimize code with local memory references by sophisticated caching schemes, which promotes the performance of the ushuffle program. we refer to figure for the running times of the ushuffle program at various values of the parameter k, where the test sequence length is fixed at the running time of the ushuffle program peaks at k = which is about three times its running time for k = then gradually decreases as k increases, and finally drops to zero at k = because, with a sequence length of the only 1024-let-preserving random sequence is the input sequence itself. this plot shows that the ushuffle program is efficient for all possible values of k. comparison of ushuffle java applet with shufflet there exist two other implementations of the euler algorithm. the dishuffle program by clote et al. implements the original version of the euler algorithm by altschul and erickson. hard-coded for shuffling rna sequences preserving dinucleotide counts, dishuffle is not a general tool for arbitrary alphabet and let sizes. another program, shufflet by coward, implements the improved version of the euler algorithm by kandel et al. for arbitrary let size k. as we have explained in the implementation section, the arborescence generation algorithm by kandel et al., while superior to the algorithm by altschul and erickson, is still inferior to wilson's algorithm; besides, its look-up table data structure is inefficient for large alphabet and let sizes. in terms of functionality, the shufflet implementation is closer to our ushuffle implementation. shufflet was written in the c programming language, and had been hosted as a web application. we were unable to perform a comprehensive comparison of ushuffle and shufflet. however, coward mentioned two experiments performed on a digital dec/alpha web server: shufflings of a dna sequence of nucleotides with k = take about seconds; shufflings of a protein sequence of amino acids with k = take less than second. we performed similar experiments with the ushuffle java applet on an apple imac computer: shufflings of a dna sequence of nucleotides with k = take about seconds; shufflings of a protein sequence of amino acids with k = take less than second. assuming comparable performances of the two computers, we estimate that our ushuffle java applet is about 15– times faster than shufflet in the experiment on nucleotides, and about times faster than shufflet in the experiment on amino acids. we certainly understand the difficulty of such a comparison: a web server in versus a desktop computer in 2005; a c program in a native machine versus a java applet in a virtual machine. nevertheless the comparison illustrates the better scalability of our ushuffle java applet for large let sizes. the difference between the two performance ratios, 15– versus suggests that ushuffle remains efficient even for large let size, while shufflet becomes more inefficient as the let size increases, due to the use of the inefficient look-up table data structure by kandel et al..  the ushuffle tool is based on superior graph algorithms and is carefully engineered to be extremely efficient. it achieves maximum flexibility by allowing arbitrary alphabet size and let size, and is available in many forms for different kinds of users. we believe ushuffle is a useful tool for the bioinformatics community. availability and requirements project name: ushuffle. project home page: operating systems: platform independent. programming languages: c, java, c#, perl, python. other requirements: none. licence: freebsd. any restrictions to use by non-academics: none. authors' contributions mj designed the software and the experiments, implemented the c and java versions of ushuffle program, and wrote the technical report. ja ported the java program to c#, investigated software licenses, and designed the ushuffle logo. jg wrote the test program, performed the experiments, and wrote the instructions for perl integration. mm wrote the java applet interface and the instructions for python integration. all authors reviewed the source code and contributed to the home page construction. note 1macos. firefox., java. 2dell xps m1710: ghz intel dual-core processor, gb dual in-line ram; microsoft windows vista business edition, cygwin. gcc. with -o option. 
source code in c, java, and c#, and integration instructions for perl and python are provided. it is a useful tool for the bioinformatics community. randomly shuffled sequences are routinely used in sequence analysis to evaluate the statistical significance of a biological sequence. we present a sequence analysis tool for generating uniform random permutations of biological sequences that preserve the exact k-let counts. the ushuffle tool implements the latest variant of the euler algorithm and uses wilson's algorithm in the crucial step of arborescence generation. the ushuffle tool achieves maximum flexibility by allowing arbitrary alphabet size and let size. the ushuffle tool surpasses existing implementation of the euler algorithm in both performance and flexibility. in many cases, biologists need sophisticated shuffling tools that preserve not only the counts of distinct letters but also higher-order statistics such as doublet counts, triplet counts, and, in general, k-let counts. it can be used as a command-line program, a web application, or a utility library.
 randomly shuffled sequences are routinely used in sequence analysis to evaluate the statistical significance of a biological sequence. in many cases, biologists need sophisticated shuffling tools that preserve not only the counts of distinct letters but also higher-order statistics such as doublet counts, triplet counts, and, in general, k-let counts.  we present a sequence analysis tool for generating uniform random permutations of biological sequences that preserve the exact k-let counts. the ushuffle tool implements the latest variant of the euler algorithm and uses wilson's algorithm in the crucial step of arborescence generation. it is carefully engineered and extremely efficient. the ushuffle tool achieves maximum flexibility by allowing arbitrary alphabet size and let size. it can be used as a command-line program, a web application, or a utility library. source code in c, java, and c#, and integration instructions for perl and python are provided.  the ushuffle tool surpasses existing implementation of the euler algorithm in both performance and flexibility. it is a useful tool for the bioinformatics community. 