How to plot a gene graph for a DNA sequence say ATGCCGCTGCGC?

Question

I need to generate a random walk based on the DNA sequence of a virus, given its base pair sequence of 2k base pairs. The sequence looks like "ATGCGTCGTAACGT". The path should turn right for an A, left for a T, go upwards for a G and downwards for a C. How can I use either Matlab, Mathematica or SPSS for this purpose?

Answer 1

I did not previously know of Mark McClure's blog about Chaos Game representation of gene sequences, but it reminded me of an article by Jose Manuel Gutiérrez (The Mathematica Journal Vol 9 Issue 2), which also gives a chaos game algorithm for an IFS using (the four bases of) DNA sequences. A detailed description may be found here (the original article).

The method may be used to produce plots such as the following. Just for the hell of it, I've included (in the RHS panels) the plots generated with the corresponding complementary DNA strand (cDNA).

• Mouse Mitochondrial DNA (LHS) and its complementary strand (cDNA) (RHS).

These plots were generated from GenBank Identifier gi|342520. The sequence contains 16295 bases.

(One of the examples used by Jose Manuel Gutiérrez. If anyone is interested, plots for the human equivalent may be generated from gi|1262342).

• Human Beta Globin Region (LHS) and its cDNA (RHS)

Generated from gi|455025| (the example used my Mark McClure). The sequence contains 73308 bases

There are pretty interesting plots! The (sometimes) fractal nature of such plots is known, but the symmetry obvious in the LHS vs RHS (cDNA) versions was very surprising (at least to me).

The nice thing is that such plots for any DNA sequence may be very easily generated by directly importing the sequence (from, say, Genbank), and then using the power of Mma.
All you need it the accession number! ('Unknown' nucleotides such as "R" may need to be zapped) (I am using Mma v7).

The Original Implimenation (slightly modified) (by Jose Manuel Gutiérrez)

Important Update

On the advise of Mark McClure, I have changed Point/@Orbit[s, Union[s]] to Point@Orbit[s, Union[s]].

This speeds things up very considerably. See Mark's comment below.

Orbit[s_List, {a_, b_, c_, d_}] :=    OrbitMap[s /. {a -> {0, 0}, b -> {0, 1}, c -> {1, 0},       d -> {1, 1}}]; OrbitMap =    Compile[{{m, _Real, 2}}, FoldList[(#1 + #2)/2 &, {0, 0}, m]]; IFSPlot[s_List] :=   Show[Graphics[{Hue[{2/3, 1, 1, .5}], AbsolutePointSize[2.5],      Point @ Orbit[s, Union[s]]}], AspectRatio -> Automatic,    PlotRange -> {{0, 1}, {0, 1}},    GridLines -> {Range[0, 1, 1/2^3], Range[0, 1, 1/2^3]}] 

This gives a blue plot. For green, change Hue[] to Hue[{1/3,1,1,.5}]

The following code now generates the first plot (for mouse mitochondrial DNA)

 IFSPlot[Flatten@       Characters@        Rest@Import[          "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=\     nucleotide&id=342520&rettype=fasta&retmode=text", "Data"]] 

To get a cDNA plot I used the follow transformation rules (and also changed the Hue setting)

IFSPlot[    ....   "Data"] /. {"A" -> "T", "T" -> "A", "G" -> "C",     "C" -> "G"}] 

Thanks to Sjoerd C. de Vries and telefunkenvf14 for help in directly importing sequences from the NCBI site.

Splitting things up a bit, for the sake of clarity.

Import a Sequence

mouseMitoFasta=Import["http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&id=342520&rettype=fasta&retmode=text","Data"]; 

The method given for importing sequences in the original Mathematica J. article is dated.

A nice check

First@mouseMitoFasta

Output:

{>gi|342520|gb|J01420.1|MUSMTCG Mouse mitochondrion, complete genome} 

Generation of the list of bases

mouseMitoBases=Flatten@Characters@Rest@mouseMitoFasta 

Some more checks

{Length@mouseMitoBases, Union@mouseMitoBases,Tally@mouseMitoBases} 

Output:

{16295,{A,C,G,T},{{G,2011},{T,4680},{A,5628},{C,3976}}} 

The second set of plots was generated in a similar manner from gi|455025. Note that the sequence is long!

{73308,{A,C,G,T},{{G,14785},{A,22068},{T,22309},{C,14146}}} 

One final example (containing 265922 bp), also showing fascinating 'fractal' symmetry. (These were generated with AbsolutePointSize[1] in IFSPlot).

The first line of the fasta file:

{>gi|328530803|gb|AFBL01000008.1| Actinomyces sp. oral taxon 170 str. F0386 A_spOraltaxon170F0386-1.0_Cont9.1, whole genome shotgun sequence}

The corresponding cDNA plot is again shown in blue on RHS

Finally, Mark's method also gives very beautiful plots (for example with gi|328530803), and may be downloaded as a notebook.