Using Mathematica to Teach Bioinformatic Skills

Brian G. Higgins
Department of Chemical Engineering & Materials Science
University of California, Davis

Mathematica Developer Conference, April 2003

Introduction

In this talk I will discuss how Mathematica can be used to manipulate and process biological data using built-in string manipulation functions. These capabilities do not seem to be well appreciated by the bioinformatics community who traditionally have used Perl as their programming environment. Recent developments have made Mathematica an attractive programming environment for managing and analyzing biological data. For example, with the  J/Link toolkit it is now possible to integrate Mathematica and Java so that modules from  BioJava  can be accessed directly from Mathematica, thereby extending the Mathematica computing environment for bioinformatics. One of Mathematica's newest features is its ability to import/export XML data, which can be parsed and manipulated using Mathematica's  SymbolicXML function.

In this presentation I will give a brief overview  of several bioinformatics problems and then I will show  how Mathematica  can be used to teach programming skills  for bioinformatics. I will discuss the parsing and management of biological data from various sources (GenBank, PDP, BLAST). I will also demonstrate how Mathematica can be used to manipulate DNA and protein sequences.

This notebook makes extensive use of Symbolic XML function; thus users should use Mathematica Version 4.2 or later. The XML files can be downloaded from my website. The notebook also accsses several text files. These can also be found at 
		http://www.higgins.ucdavis.edu/biomath.php

Molecular Biology Data

Bioinformatics

"Study of how information is represented and transmitted in biological systems" (Bergeron, 2003)

Involves the use of computer technology for data analysis, data storage, data mining, knowledge discovery, and biological modeling.

Applications involving molecular biology data:

- Pattern matching: sequence alignment (e.g., Blast)
	- Restriction maps (e.g., clone libraries, genome maps)
	- Evolutionary relationships (e.g., trees and sequences, parsimony, sequence assembly)
	- RNA secondary structure
	- Data analysis (e.g., microarrays, cluster analysis)
	- Development of biological databases  for data mining (e.g., web based tools)
	- Modeling and simulation (e.g., metabolic pathway analysis, protein structure)
					etc.,

Impact

The genome project has resulted in an enormous amount of molecular biology data that can be retrieved from online sources.

Researchers and students now have access to real-world data sets that can be used for discovery, manipulation, aggregation and analysis.

Academic programs and courses in bioinformatics are now offered by most major research universities.

Tools for Accessing Molecular Biology Data

Computing Environment

Computer programming is central for manipulating archival molecular biology data.
Perl is the preferred programming language. 
Other programming environments are Java, Python, C, C++
Unix is the preferred OS.

Positive attributes of Perl

-  Strong string manipulation/processing capabilities
	-  Suitable for CGI-programming ( e.g linking a web server with a web application)
	-  Easy to learn (?)
	- Modules for common bioinformatic operations

Concerns with using Perl

-  Not adequate for complex projects in computational biology
	- More efficient to do complex calculations/graphics with other  languages
	- Non-computer science students often are not familiar with Unix/Linux OS

Is there a role for Mathematica?

Role for Mathematica

Mathematica has powerful built-in string manipulation functions. 
Mathematica is an integrated development environment for doing symbolic, numerics and data visualization.
With Mathematica's J/Link toolkit it is now possible to integrate Mathematica and Java so that existing modules  written in Java can be accessed directly from Mathematica , thereby extending the Mathematica computing environment for bioinformatics. 
Mathematica's newest features is its ability to import/export XML data, which can be parsed and manipulated using Mathematica's  SymbolicXML function.

Central Dogma in Biology

Initialization Cells

Needs["Graphics`Legend`"];

DNABases := 
      Module[{BaseNmaes, col1, leg, pp1}, 
        baseNames = {"Cytosine(C)", "Thymine(T)", "Guanine(G)", "\
Adenine(A)"}; 
col1[n_] := 
          Switch[n, 0, RGBColor[1, 0, 0], 1, RGBColor[1, 1, 0], 2, 
            RGBColor[0, 1, 0], 3, 
            RGBColor[0, 1, 1]]; 
leg = 
          Table[{col1[Mod[n, 4]], Part[baseNames, n + 1]}, {n, 0, 
              3}]; 
pp1 = 
          Show[Graphics[{RGBColor[1, 0, 0], 
                Disk[{0, 0}, 3, {-π/4, π/4}], RGBColor[1, 1, 0], 
                Disk[{0, 0}, 3, {π/4, 3  π/4}], RGBColor[0, 1, 0], 
                Disk[{0, 0}, 3, {3  π/4, 5  π/4}], RGBColor[0, 1, 1], 
                Disk[{0, 0}, 3, {5  π/4, 7  π/4}]}], 
            DisplayFunction -> Identity, 
            AspectRatio -> 1]; 
ShowLegend[
          pp1, {leg, LegendShadow -> None, 
            LegendSize -> {1.4, 1.4}, 
            LegendPosition -> {1.3, - .5}}, 
          DisplayFunction -> $DisplayFunction]]

\(\(AminoAcids = {{"Alanine", "Ala", "A"}, {"Arginine", "Arg", "R"}, {"Aspartic Acid", "Asp", "D"}, {"Asparginine", "Asn", "N"}, {"Cysteine", "Cys", "C"}, {"Glutamic Acid", "Glu", "E"}, {"Glutamine", "Gln", "Q"}, {"Glycine", "Gly", "G"}, {"Histine", "His", "H"}, {"Isoleucine", "Ile", "I"}, {"Leucine", "Leu", "L"}, {"Lysine", "Lys", "K"}, {"Methionine", "Met", "M"}, {"Phenylalanine", "Phe", "F"}, {"Proline", "Pro", "P"}, {"Serine", "Ser", "S"}, {"Threonine", "Thr", "T"}, {"Tryptophan", "Trp", "W"}, {"Tyrosine", "Tyr", "Y"}, {"Valine", "Val", "V"}};\)\)

\(AATable := 
      DisplayForm[
        GridBox[Prepend[
            AminoAcids, {"Amino Acid", \ "3-Letter Code", \ "1-Letter Code"}], ColumnLines -> {True, True}, 
          RowLines -> {True, False}, ColumnSpacings -> {1, 1}, 
          RowSpacings -> {2, 1}, ColumnAlignments -> {Left, Center, Center}, 
          GridFrame -> True]]\)

DNA Alphabet

DNA is  a polymer made up of small molecules called nucleotides. 
There are 4 nucleotides, identified by the bases:
		 adenine (a), cytosine (c), guanine (g), and thymine (t)

We can view DNA as a word defined over the four letter alphabet  (Subscript[𝒜, DNA])={a,c,g,t}

RNA is also a word defined over the alphabet of ribonucleotides (Subscript[𝒜, RNA])={a,c,g,u}

The central dogma (Crick, 1958) ( now viewed as a simplified dogma)

DNA⟶RNA⟶Protein

Both RNA and proteins are made from instructions in the DNA.  
	Replication- DNA is copied
	Transcription - DNA templates RNA
	Translation - RNA templates proteins

DNABases;

Amino acids are encoded by triplets of nucleotides of DNA, called codons. There are 64 possible codons, but 22  encoded amino acids( 20 common amino acids, plus stop and start codons). Here is a list of the common amino acids

AATable

RNA that is translated into protein is known as messenger RNA (mRNA)

The phase of codon reading is called the reading frame. There is a total of 6 reading frames, including reading frames associated with the complementary DNA strand.

For convenience we will not represent the codons in terms of the base Uracil. During transcription RNA converts Thymine into Uracil. Our codons will use Thymine.

There are 20 common amino acids (22 genetically encoded amino acids) and they are defined by the codons according to the 3 letter base sequences (codons)

The codon-amino acid transcription rules can be represented visually as

The transcription wheel shows codon redundancy directly.

We can also represent the codon-amino acids transcriptions in terms of Mathematica rules

\(\(CodonRules = {"tca" -> "S", "tcc" -> "S", "tcg" -> "S", "tct" -> "S", "ttc" -> "F", \ \ \ "ttt" -> "F", "tta" -> "L", "ttg" -> "L", "tac" -> "Y", "tat" -> "Y", "taa" -> "_", "tag" -> "_", "tgc" -> "C", "tgt" -> "C", "tga" -> "_", "tgg" -> "W", "cta" -> "L", "ctc" -> "L", "ctg" -> "L", "ctt" -> "L", "cca" -> "P", "ccc" -> "P", "ccg" -> "P", "cct" -> "P", "cac" -> "H", "cat" -> "H", "caa" -> "Q", "cag" -> "Q", "cga" -> "R", "cgc" -> "R", "cgg" -> "R", "cgt" -> "R", "ata" -> "I", "att" -> "I", "atc" -> "I", "atg" -> "M", "aca" -> "T", "acc" -> "T", "acg" -> "T", "act" -> "T", "aac" -> "N", "aat" -> "N", "aaa" -> "K", "aag" -> "K", "agc" -> "S", "agt" -> "S", "aga" -> "R", "agg" -> "R", "gta" -> "V", "gtc" -> "V", "gtg" -> "V", "gtt" -> "V", "gca" -> "A", "gcc" -> "A", "gcg" -> "A", "gct" -> "A", "gac" -> "D", "gat" -> "D", "gaa" -> "E", "gag" -> "E", "gga" -> "G", "ggc" -> "G", "ggg" -> "G", "ggt" -> "G"};\)\)

Importing Biological Data into Mathematica

Importing Genbank file into Mathematica

In this first example we are going to show how we can write a simple parser that will allow us to import GenBank (Genetic Sequence Data Bank) files into Mathematica. The database can be found at the following URL:
				http://www.ncbi.nlm.nih.gov/Genbank/index.html

The file ID we are going to examine is U49845.  This file gives the amino acid sequence for several genes from the bacterium Saccharomyces cerevisiae, which is commonly known as baker's or budding yeast. The DNA sequence for these genes comes from  chromosome IX.

We use of the ReadList function to read in the text file from my directory: (The file U49845.txt can be downloaded from http://www.higgins.ucdavis.edu/biomath.php )

myGenBankfile = 
        ReadList[Experimental`FileBrowse[False], Word,  
          WordSeparators -> {"\r"}];

Here is what the first 10 lines of the file look like

Table[myGenBankfile[[i]], {i, 1, 10}] // TableForm

Note that the text at the beginning of each line is not always defined by a record type. In the above output we have a LOCUS, DEFINITION, ACCESSION, etc. as record types. But there are also lines without a record type, e.g. the line containing the string

		"Eukaryota; Fungi; Ascomycota;....."

Note also that each record of the file is a string with lots of white space.

The number of records in this file can be determined using the Length function

Length[myGenBankfile]

Suppose we are interested in getting the DNA sequence data for this record. Here is a portion of the DNA sequence

Table[myGenBankfile[[i]], {i, 81, 90}] // TableForm

We see it is not possible to select the DNA rows from this list by record type as the DNA rows have no record type. However,  we note that the DNA sequence begins after the row with the record type  ORIGIN. Thus we will define a parser that determines the row position containing the word ORIGIN, and then  use that information to extract the string data for the DNA sequence. We convert each row into a Stream, use ReadList, with the Null string as a record separator, drop the first element which is the position number of the base in that record, and then use StringJoin to obtain the final  base sequence for our mRNA. Here is our parser

\(genBankParser[file_] := 
      Module[{pat1, OriginRowNumber, DNAsequence, DNAseqString, s}, 
        pat1 = First[
            Select[file, 
              StringMatchQ[#, "ORIGIN*"] &]]; 
OriginRowNumber = 
          First[Flatten[
              Position[file, 
                x_String /; 
                  x == pat1]]]; 
DNAsequence = 
          Table[file[\([i]\)], {i, OriginRowNumber + 1, 
              Length[file] - 1}]; 
For[
          DNAseqString = ""; i = 1, i < Length[DNAsequence] + 1, 
          s = StringToStream[DNAsequence[\([i]\)]]; 
          DNAseqString = 
            DNAseqString <> 
              StringJoin[
                Drop[ReadList[s, Word, RecordSeparators -> {" "}], 
                  1]]; \(i++\)]; 
DNAseqString]\)

DNAseq1 = genBankParser[myGenBankfile]

StringLength[DNAseq1]

Importing XML Data from the EBI Database

Several public databases (DDBJ/EMBL/GenBank) of molecular biology data on the internet now support XML format for presenting results from a query search. In this section we are going to explore how to manipulate molecular biology data from a query search using the EMBL-EBI nucleotide database. Details of the XML project at EMBL can be found at the following link

				http://www.ebi.ac.uk/xembl/

The particular XML format we are going to use is called BSML. With Mathematica 4.2 it is now possible to import XML files into Mathematica and then convert them to SymbolicXML, There are several ways to import XML data into Mathematica. One can use the Import function or XMLGet function. In this notebook we will use the latter. The XMLGet function allows one to query a URL directly from the notebook. In the next section we will perform the URL query to obtain the XML data from

	"http://www.higgins.ucdavis.edu/xmldata/xmldata/SC49845.xml" .
 
 Before we do, let's look at the file's format. The original XML document has 147 lines of code; here are the first 15 lines of code:

<?xml version='1.0' encoding='UTF-8'?>
<?format DECIMAL="."?>
<Bsml>
 <Definitions>
  <Sequences>
   <Sequence id='SC49845'
       ic-acckey='U49845'
       title='SC49845'
       comment='Saccharomyces cerevisiae TCP1-beta gene, partial cds; and Axl2p (AXL2) and Rev7p (REV7) genes, complete cds. '
       length='5028'
       topology='linear'
       molecule='dna'>
    <Attribute name='version'
        content='U49845.1'/>

The XML document can been parsed by Mathematica's XMLGet function to obtain a SymbolicXML expression:

xmlbiodata = 
        XML`Parser`XMLGet["http://www.higgins.ucdavis.edu/xmldata/SC49845.\
xml"] ;

Shown below are the first 30 lines of the SymbolicXML expression

Short[xmlbiodata, 30]

Here is an annotated version of the output

XMLObject["Document"][\({\(XMLObject["Declaration"]\)["Version" -> "1.0", "Encoding" -> "UTF-8"], \(XMLObject["ProcessingInstruction"]\)["format", "DECIMAL=".""]}\),XMLElement["Bsml",{},{XMLElement["Definitions",{},{XMLElement["Sequences",{},{XMLElement["Sequence",\({"id" -> "SC49845", "ic-acckey" -> "U49845", "title" -> "SC49845", "comment" -> "Saccharomyces cerevisiae TCP1-beta  gene, partial cds; and Axl2p (AXL2) and Rev7p (REV7) genes, complete cds. ", "length" -> "5028", "topology" -> "linear", "molecule" -> "dna"}\),\({XMLElement["Attribute", {"name" -> "version", "content" -> "U49845.1"}, {}], «12»}\)]}]}]}],{}]

If we inspect the data file after it has been parsed into SymbolicXML we see that it is a highly nested function, consisting of XMLElements. We can extract data from the parsed file using the built-in functions in Mathematica. In particular, we will make use of the Cases function with appropriately designed  patterns/rules. Here is a simple example that extracts the data shown above in red:

Cases[xmlbiodata, 
      XMLElement ["Sequence", {x__}, __] -> 
        x, ∞]

To extract the nucleotide sequence data from the XML file, we must search within the Seq-data tag, which is defined by the following XMLElement:
		XMLElement["Seq-data", {}, {seq}]
where seq is the nucleotide sequence. We can extract the data using the following pattern

DNAseq2 = 
        First[Cases[xmlbiodata, 
            XMLElement["Seq-data", {}, {x_}] :> 
              StringReplace[
                x, " " -> ""], ∞]];

Note we make use of StringReplace to remove any null spaces in the sequence due to spurious formatting in the XML file. The number of bases in this sequence is

StringLength[DNAseq2]

We can check that our XML parser gives the same result found earlier using the GenBank file

DNAseq2 == DNAseq1

DNAseq2

Manipulating DNA Sequences

We can use the concept of base pairing to determine complement of a DNA sequence using the following function

\(DNAComplement[x_String] := 
      StringReplace[
        x, {"a" -> "t", "t" -> "a", "g" -> "c", "c" -> "g"}]\)

Let's apply this function to our DNA sequence

DNAseq2c = DNAComplement[DNAseq2];

We can visualize the DNA using the following function (thanks to David Parks, Mathgroup for assistance)

\(SequenceView[DNA1_, DNA2_, start_, rowlength_, rows_] := 
      Module[{subdata1, subdata2, subdata, \ 
          rowgrid1}, 
rowgrid1[
            subdata_]\  := 
\ \ \ \ \ \ GridBox[{subdata[\([1]\)], \ 
              MapAt[StyleBox[#, FontColor -> RGBColor[1, 0, 0]] &, \ 
                Table["|", \ {rowlength}], 
                Position[subdata[\([1]\)], 
                  StyleBox[x__]]], 
\ \ \ \ \ \ \ \ \ subdata[\([2]\)]}, \ 
            GridBaseline\  -> \ 
              Bottom]; 
\ \ subdata1\  = 
\ \ \ \ \ \ Take[
            Characters[DNA1], \ {start, \ 
              Min[start\  + \ rows\ rowlength\  - \ 1, \ 
                Length[Characters[
                    DNA1]]]}]; 
\ \ \ \ subdata1\  = \ 
          Partition[subdata1, \ rowlength, \ 
            rowlength, \ {1, \ 
              1}, \ " "]; 
subdata2 = 
\ \ \ \ \ \ Take[Characters[DNA2], \ {start, \ 
              Min[start\  + \ rows\ rowlength\  - \ 1, \ 
                Length[Characters[
                    DNA2]]]}]; 
\ \ \ \ subdata2\  = \ 
          Partition[subdata2, \ rowlength, \ 
            rowlength, \ {1, \ 1}, \ " "]; 
\ \ \ \ subdata\  = 
          Map[\ {rowgrid1[#]}\  &\ , 
            Transpose[{subdata1, subdata2}]]; 
\ \ \ \ DisplayForm[
          FrameBox[GridBox[subdata, \ RowLines\  -> \ True]]]]\)

The SequenceView function takes 5 arguments. The first two arguments are strings that represent the complementary DNA strands. The third argument is the position of the base pair where viewing begins; the fourth argument is the number of base pairs in a row; and the fifth argument is the number of rows. Here is a view of the first 125 base pairs of the sequence.

SequenceView[DNAseq2, DNAseq2c, 1, 25, 10]

We can readily determine the number of individual bases in our sequence

Map[
      Length[Cases[
            Characters[
              DNAseq2], #]] &, {"a", "c", "g", "t"}]

Comparing DNA Sequences of Equal Length

We can use the results from the previous section to write a function that compares two sequences and then displays the sequence in a cartoon and shows the bases that are aligned in a different color. First, we have the following function that generates a random DNA sequence

randomnucl := 
      Switch[Random[], _?(# < 0.25 &), "a", _?(# < 
                0.5 &), "t", _?(# < 
                0.75 &), "c", _, "g"]

Next, we generate two sequences at random each with 500 bases

DNASeqRandom1 = Table[randomnucl, {500}];
DNASeqRandom2 = Table[randomnucl, {500}];

The following function does the comparison and displays the results

\(CompareSequences[seq1_, seq2_, start_, rowlength_, rows_] := 
      Module[{pos, seq11, seq21, subdata1, subdata2, subdata, \ 
          rowgrid1}, 
pos = 
          Position[\((seq1 - seq2)\), 
            x_ /; x == 0]; 
seq11 = 
          MapAt[ToUpperCase, seq1, pos] /. 
            x_String?UpperCaseQ :> 
              StyleBox[ToLowerCase[x], 
                FontColor -> 
                  RGBColor[1, 0, 0]]; 
seq21 = 
          MapAt[ToUpperCase, seq2, pos] /. 
            x_String?UpperCaseQ :> 
              StyleBox[ToLowerCase[x], 
                FontColor -> 
                  RGBColor[1, 0, 0]]; 
rowgrid1[
            subdata_]\  := 
\ \ \ \ \ \ GridBox[{subdata[\([1]\)], \ 
              MapAt[StyleBox[#, FontColor -> RGBColor[1, 0, 0]] &, \ 
                Table["|", \ {rowlength}], 
                Position[subdata[\([1]\)], 
                  StyleBox[x__]]], 
\ \ \ \ \ \ \ \ \ subdata[\([2]\)]}, \ 
            GridBaseline\  -> \ 
              Bottom]; 
\ \ subdata1\  = 
\ \ \ \ \ \ Take[
            seq11, \ {start, \ 
              Min[start\  + \ rows\ rowlength\  - \ 1, \ 
                Length[seq11]]}]; 
\ \ \ \ subdata1\  = \ 
          Partition[subdata1, \ rowlength, \ 
            rowlength, \ {1, \ 
              1}, \ " "]; 
subdata2 = 
\ \ \ \ \ \ Take[seq21, \ {start, \ 
              Min[start\  + \ rows\ rowlength\  - \ 1, \ 
                Length[seq21]]}]; 
\ \ \ \ subdata2\  = \ 
          Partition[subdata2, \ rowlength, \ 
            rowlength, \ {1, \ 1}, \ " "]; 
\ \ \ \ subdata\  = 
          Map[\ {rowgrid1[#]}\  &\ , 
            Transpose[{subdata1, subdata2}]]; 
\ \ \ \ DisplayForm[
          FrameBox[GridBox[subdata, \ RowLines\  -> \ True]]]]\)

Here is an example

CompareSequences[DNASeqRandom1, DNASeqRandom2, 10, 25, 10]

Comparing DNA Sequences using Dot Plots

A dot plot is simple method for analyzing and comparing two DNA sequences. In its simplest form a matrix is constructed where the rows are the character positions in sequence A and the columns are the character positions in sequence B. A dot is placed in location (i,j) if the (i^th)character of sequence A matches the (j^th)character of sequence B.  Regions that are similar are revealed as a set of diagonal dots. Random matches are given by isolated dots.  Because two DNA sequences will have a high probability of random matches, the noise can be filtered by using a sliding window that compares blocks of DNA bases. Two blocks are said to match if a mismatch criterion(level) is satisfied.  If the mismatch level is zero then all bases in the two windows must match. In this example we will show how to construct dot plots when the mismatch level is zero.

In this example we will examine DNA sequences that produce human hemoglobins, which is a protein used by red blood cells to distribute oxygen to other tissues and cells in the body. In a normal adult hemoglobin A is the most prevalent form  comprising about 95% of all hemoglobin. Two other forms also occur : hemoglobin (Subscript[A, 2]) and hemoglobin F. Hemoglobin A is composed of two alpha and two beta globins, while hemoglobin F consists of two alpha and two gamma globins.  The human alpha globin gene cluster is located on chromosome 16 and spans about 30000 bases.

\(DNADotPlot[DNA1_, DNA2_, BlockSize_, offset_] := 
      Module[{dnablocks1, dnablocks2, blocklists, orderedlists, grouplists, 
          altPatterns, matchlists}, 
dnablocks1 = 
          MapIndexed[List[#, \((First[#2] - 1)\) offset + 1] &, 
            Map[StringJoin[#] &, 
              Partition[Characters[DNA1], BlockSize, offset], 
              1]]; 
dnablocks2 = 
          MapIndexed[List[#, \((First[#2] - 1)\) offset + 1] &, 
            Map[StringJoin[#] &, 
              Partition[Characters[DNA2], BlockSize, offset], 
              1]]; 
blocklists = {dnablocks1, 
            dnablocks2}; 
orderedlists = \(#[\([Ordering[#[\([All, 1]\)]]]\)] &\) /@ 
            blocklists; 
grouplists = \(Split[#, #1[\([1]\)] === #2[\([1]\)] &] &\) /@ orderedlists; 
        altPatterns = \((Alternatives @@ \((Intersection @@ 
                  orderedlists[\([All, All, 
                      1]\)])\))\); 
matchlists = 
          Transpose[\(Cases[#, {{altPatterns, _}, ___}] &\) /@ 
              grouplists]; 
ListPlot[
          Flatten[Apply[Outer[List, ##, 1] &, matchlists, {1}], 
              2] /. {{x1_String, y1_}, {x2_String, y2_}} -> {y1, y2}, 
          AspectRatio -> 1, 
          PlotStyle -> {RGBColor[0, 0, 1], PointSize[0.01]}, 
          Frame -> True, FrameLabel -> {"DNA1", "DNA2"}, 
          RotateLabel -> False, Axes -> False]]\)

The algorithm used in the above module was suggested by Allan Hayes from the Mathgroup.

We will import the DNA sequences that encode for various hemoglobin forms from the local web site, and then use our DNADotPlot function to see if there are similarities in the sequences. The files we will be importing are in XML format and were downloaded originally from GenBank. We will use the XMLGet function to query a URL directly from the notebook. We will import 4 files.

HBalpha1XML = 
        XML`Parser`XMLGet["http://www.higgins.ucdavis.edu/xmldata/\
hemaglobin_alpha1.xml"] ;

HBalpha2XML = 
        XML`Parser`XMLGet["http://www.higgins.ucdavis.edu/xmldata/\
hemaglobin_alpha2.xml"] ;

HBbetaXML = 
        XML`Parser`XMLGet["http://www.higgins.ucdavis.edu/xmldata/\
hemaglobin_beta.xml"] ;

HBgammaAXML = 
        XML`Parser`XMLGet["http://www.higgins.ucdavis.edu/xmldata/\
hemaglobin_gammaA.xml"] ;

We can extract the DNA sequences from these files by selecting the data from the appropriate XMLElement

{HBalpha1, HBalpha2, HBbeta, HBgammaA} = 
        Map[First[
              Flatten[Cases[#, 
                  XMLElement["IUPACna", _, x_] -> 
                    x, ∞]]] &, {HBalpha1XML, HBalpha2XML, 
            HBbetaXML, HBgammaAXML}];

In our first plot we will compare the two alpha sequences

DNADotPlot[HBalpha1, HBalpha2, 6, 1];

We see that there is considerable overlap except near the end of each sequence. In our next plot we compare the beta encoding sequence with the alpha encoding sequence. In this case there is no meaningful overlap

DNADotPlot[HBbeta, HBalpha1, 4, 1];

Finally we compare the gamma encoding sequence  with the beta encoding sequence. In this case there is some overlap.

DNADotPlot[HBgammaA, HBbeta, 4, 1];

Working with Amino Acid Sequences

Information about the amino acid sequences is contained in the following XMLElement:

				XMLElement["Feature", ... ....]

We can use the following pattern to extract out all "Feature" data from the XML file

Cases[xmlbiodata, 
      XMLElement["Feature", ___], ∞]

Suppose we want to extract out the data for a particular gene. The amino acid sequence for the AXL2 gene can be found in the following  record

\(Cases[xmlbiodata, 
      XMLElement["Feature", {__, "title" -> "AXL2", ___}, ___], ∞]\)

If we inspect the above record we can define the following function that extracts the amino acid sequence from the record for a particular gene as well as the start and end positions for coding, which is defined in the last XMLElement of the "Feature" XMLElement. Since not all coding begins by reading the DNA sequence of bases from left to right we must also  add in a third argument for this XMLElement that determines whether the reading is from right to left (i.e.using the complement DNA strand). This argument is given as a rule: Complement -> int , where int = 1, 2,  or 3. We define an optional  pattern for this argument since if the complement strand is not used  it is omitted.

\(geneData[geneName_String] := 
      Flatten[Cases[xmlbiodata, 
          XMLElement["Feature", {__, "title" -> 
                  geneName, ___}, {___, 
                XMLElement["Qualifier", {"value-type" -> "translation", "value" -> x_}, {}], ___, 
                XMLElement["Interval-loc", {"startpos" -> 
                      n1_, "endpos" -> n2_\ , 
                    y_: ""}, {}]}] :> {x, {n1, n2}, 
              y}, ∞], 1]\)

Here is the data for the AXL2 gene:

{AminoAcidSeq1, loc1, comp1} = geneData["AXL2"]

Let us test out whether this amino acid sequence can be reconstructed from the DNA sequence. We proceed as follows. First, we select the bases from our sequence in positions 687-3158 using the StringTake function

baseSeq1 = StringTake[DNAseq, ToExpression[loc1]];

We can convert this base sequence to an amino acid sequence using the following codon rules

\(\(CodonRules = {"tca" -> "S", "tcc" -> "S", "tcg" -> "S", "tct" -> "S", "ttc" -> "F", \ \ \ "ttt" -> "F", "tta" -> "L", "ttg" -> "L", "tac" -> "Y", "tat" -> "Y", "taa" -> "_", "tag" -> "_", "tgc" -> "C", "tgt" -> "C", "tga" -> "_", "tgg" -> "W", "cta" -> "L", "ctc" -> "L", "ctg" -> "L", "ctt" -> "L", "cca" -> "P", "ccc" -> "P", "ccg" -> "P", "cct" -> "P", "cac" -> "H", "cat" -> "H", "caa" -> "Q", "cag" -> "Q", "cga" -> "R", "cgc" -> "R", "cgg" -> "R", "cgt" -> "R", "ata" -> "I", "att" -> "I", "atc" -> "I", "atg" -> "M", "aca" -> "T", "acc" -> "T", "acg" -> "T", "act" -> "T", "aac" -> "N", "aat" -> "N", "aaa" -> "K", "aag" -> "K", "agc" -> "S", "agt" -> "S", "aga" -> "R", "agg" -> "R", "gta" -> "V", "gtc" -> "V", "gtg" -> "V", "gtt" -> "V", "gca" -> "A", "gcc" -> "A", "gcg" -> "A", "gct" -> "A", "gac" -> "D", "gat" -> "D", "gaa" -> "E", "gag" -> "E", "gga" -> "G", "ggc" -> "G", "ggg" -> "G", "ggt" -> "G"};\)\)

The way we do this is first convert the DNA string into characters using the Character function, partition the subsequent list into sublist of 3 bases, rejoin the individual characters in the sublists using StringJoin,  and then apply the codon rules to each sublist. Thus the predicted amino acid sequence generated by bases at positions 687-3158  is

predAminoAcidSeq1 = 
      StringJoin[Partition[Characters[baseSeq1], 
              3] /. {x_, y_, z_} :> StringJoin[x, y, z] /. 
          CodonRules]

If we drop the "stop" codon, we can test to see whether our predicted amino acid sequence matches the one given in the XML file

AminoAcidSeq1 == StringDrop[predAminoAcidSeq1, -1]

Let us examine the output for the second gene REV7

{AminoAcidSeq2, loc2, comp2} = geneData["REV7"]

In this case the start and stop bases used for the amino acid sequence are defined in terms of the complement DNA structure. Thus to extract out the appropriate bases we first take the complement of our DNA sequence, select the bases at positions 3300 through 4037, and then reverse the segment. Our base sequence is then

baseSeq2 = 
        StringReverse[
          StringTake[DNAComplement[DNAseq], {3300, 4037}]];

Applying the codon rules to this sequence of bases gives

predAminoAcidSeq2 = 
      StringJoin[Partition[Characters[baseSeq2], 
              3] /. {x_, y_, z_} :> StringJoin[x, y, z] /. 
          CodonRules]

Again if we drop the "stop" codon, we can test to see whether our predicted amino acid sequence matches the one given in the XML file for the REV7 gene

AminoAcidSeq2 == StringDrop[predAminoAcidSeq2, -1]

Working with Restriction Enzymes

If foreign DNA  is introduced into a bacterium, a group of enzymes renders the DNA ineffective by cleaving it, thus restricting its activity. These enzymes are called restriction enzymes. There are some 300 such enzymes.

In this example we examine ways to use and manipulate restriction enzymes. We will make use of a data file that can be found on the Rebase web site at the following URL:
		http://rebase.neb.com/rebase/rebase.files.html 
The file we will be working with is called IGSuite on the above site. Since the site does not support XML formats, it is necessary to download the file and then import that file into Mathematica. To keep this notebook self-contained, we will define the following function that list the data for the first 40 enzymes

\(\(myResData = {"AaaI (XmaIII)                     C^GGCCG", "AacI (BamHI)                      GGATCC", "AaeI (BamHI)                      GGATCC", "AagI (ClaI)                       AT^CGAT", "AaqI (ApaLI)                      GTGCAC", "AarI                              CACCTGCNNNN^", "AarI                              ^NNNNNNNNGCAGGTG", "AasI (DrdI)                       GACNNNN^NNGTC", "AatI (StuI)                       AGG^CCT", "AatII                             GACGT^C", "AauI (Bsp1407I)                   T^GTACA", "AbaI (BclI)                       T^GATCA", "AbeI (BbvCI)                      CC^TCAGC", "AbeI (BbvCI)                      GC^TGAGG", "AbrI (XhoI)                       C^TCGAG", "AcaI (AsuII)                      TTCGAA", "AcaII (BamHI)                     GGATCC", "AcaIII (MstI)                     TGCGCA", "AcaIV (HaeIII)                    GGCC", "AccI                              GT^MKAC", "AccII (FnuDII)                    CG^CG", "AccIII (BspMII)                   T^CCGGA", "Acc16I (MstI)                     TGC^GCA", "Acc36I (BspMI)                    ACCTGCNNNN^", "Acc36I (BspMI)                    ^NNNNNNNNGCAGGT", "Acc38I (EcoRII)                   CCWGG", "Acc65I (KpnI)                     G^GTACC", "Acc113I (ScaI)                    AGT^ACT", "AccB1I (HgiCI)                    G^GYRCC", "AccB2I (HaeII)                    RGCGC^Y", "AccB7I (PflMI)                    CCANNNN^NTGG", "AccBSI (BsrBI)                    CCG^CTC", "AccBSI (BsrBI)                    GAG^CGG", "AccEBI (BamHI)                    G^GATCC", "AceI (TseI)                       G^CWGC", "AceII (NheI)                      GCTAG^C", "AceIII                            CAGCTCNNNNNNN^", "AceIII                            ^NNNNNNNNNNNGAGCTG", "AciI                              C^CGC", "AciI                              G^CGG"};\)\)

Manipulating Restriction Enzyme data

The file called myResData contains the restriction enzyme information. We need to edit each one of these lines. For our purposes we will remove the names in parenthesis and then parse the actual recognition sequence information so that it is useful for doing motif searches on a target sequence.

To keep the discussion simple we will first show the steps we take to modify the data in a single row, say row 30

myResData[[30]]

Note that the recognition sequence contains letters that are not the usual DNA bases "A","G","C","T" The  claret symbol "^" indicates where the cut is made. In the above example the AccB2I enzyme cuts the DNA strand at the RGCGCY site between the C and Y bases. The extended alphabet  used in the recognition sequence is defined by the following rules

\(\(RestrictionSiteCodeRules = {"R" -> "[GA]", "Y" -> "[CT]", "M" -> "[AC]", "K" -> "[GT]", "S" -> "[GC]", "W" -> "[AT]", "B" -> "[CGT]", "D" -> "[AGT]", "H" -> "[ACT]", "V" -> "[ACG]", "N" -> "[ACGT]"};\)\)

Thus the letter Y in the recognition sequence RGCGCY stands for either C or T, and the letter R stands for either G or A. So the AccB2I enzyme will cut at all of the following sites in a  DNA strand:
		GGCGCC,GGCGCT, AGCGCC, AGCGCT

We will find it convenient for later processing to convert each string of the data in the rows of myResData into separate elements. The easiest way to do this is to convert the string to a stream and then use ReadList as shown below for row 30. We also remove the name in parenthesis with a simple rule.

s = StringToStream[myResData[[30]]];
enzymeData = 
      ReadList[s, Word, 
          RecordSeparators -> " "] /. {x_String, y_String, 
            z_String} -> {x, z}

It is a simple matter to use the above ideas to write a short function that works on  all elements of our list:

enzymeListData[file_] := 
      Module[{}, 
        Table[ReadList[StringToStream[file[[i]]], Word, 
              RecordSeparators -> " "] /. {x_String, y_String, 
                z_String} -> {x, z}, {i, 1, Length[file]}]]

Here is a list of the first 40 enzymes after the file has been parsed.

myEnzymeList = enzymeListData[myResData];

We will find it convenient to have a few utility functions that operate on this list. First, it will be convenient to have  look-up function : For a given enzyme name, we can read from the file the recognition sequence. Here is a simple function that does this task

enzymeLookup[enzymeName_String, enzymeList_List] := 
      First[Cases[
          enzymeList, {x_String /; StringMatchQ[x, enzymeName], 
              y_String} -> y, ∞]]

In the following example we specify the restriction enzyme "AccB2I"

enzymeLookup["AccB2I", myEnzymeList]

We also need a function that translates the recognition sequence into the 4 nucleotide bases using the restrictionSiteCodeRues rules.

resSiteString[str_String] := 
      StringReplace[str, 
        Flatten[{"^" -> "", 
            RestrictionSiteCodeRules}]]

We can combine this function with the previous one so that by specifying the name of the enzyme we get the code for the recognition sequence

resSiteString[enzymeLookup["AccB2I", myEnzymeList]]

Another  function we will need is one that generates a list of all the possible recognition sequence for a given enzyme.What we need is a  function that can parse the recognition sequence string that is obtained from our resSiteString function. The following function takes the string say GCGC[CG] and returns  a list of the various characters in the string. i.e. GCGCC and GCGCG

\(resSiteList[str_String] := 
      Module[{s, s1, s2, s3, stringCombine, str1, 
          str2}, 
stringCombine[str1_List, str2_List] := 
          Flatten[Outer[StringJoin, str1, str2]]; 
        s = StringToStream[
            StringReplace[
              str, {"[" -> "#[", "]" -> "]#"}]]; 
s1 = 
          ReadList[s, Word, 
            WordSeparators -> {"#"}]; 
s2 = 
          s1 //. {y___, x_String /; StringMatchQ[x, "[*"], z___} :> {y, 
                Characters[
                  StringReplace[x, {"[" -> "", "]" -> ""}]], 
                z}; 
s3 = 
          Map[Flatten[#] &, 
            s2 /. x_String -> {x}]; 
Fold[
          stringCombine, First[s3], Rest[s3]]]\)

Let us try it out on the string "GCGC[CG][AG]T"

resSiteList["GCGC[CG][AG]T"]

Again we can concatenate these functions as follows and from the name of the enzyme we get all possible restriction sites for our enzyme "AccB2I"

resSiteList[
      resSiteString[enzymeLookup["AccB2I", myEnzymeList]]]

For some of the enzymes the number of restriction sites can be quite large. Consider the case for "AceIII"

resSiteString[enzymeLookup["AceIII", myEnzymeList]]

This  enzyme would generate (4^7)=16385 possible sites.

The last utility function that we need is one the determines the location in the target DNA where the cut is made.Thus enzyme "AccB2I" has the recognition sequence RGCGC^Y which means it cuts between the 5th and 6th base in the site. The following function determines the base location to the left of the cut:

cutLocation[str_String] := 
      First[Flatten[StringPosition[str, "^"]]] - 1

Here is a test with the AccB2I enzyme

cutLocation[enzymeLookup["AccB2I", myEnzymeList]]

In the next section we show use these functions to determine the location a particular enzyme cuts a DNA sequence

Motif Searches

Suppose we are given a DNA string and we want to determine the restriction sites on that string for a given enzyme. We will use the following DNA strand for our  example, which can be down loaded from the GenBank site:
			http://www.ncbi.nlm.nih.gov/Genbank/index.html
 For this purpose we are going to work with the file for the Homo Sapiens mRNA for AIRE Protein.. The nucleotide Accession Number  for  file is Z97990, while its unique identifier (UID) is 2665370. We can use Mathematica's XML import capabilities to bring this file into our notebook  in the form of Symbolic XML

xmldna = 
        XML`Parser`XMLGet["http://127.0.0.1/xmldata/Z97990.xml"] ;

We can extract the DNA sequence by using the following pattern in conjunction with Cases

myDNA = 
      First[Flatten[
          Cases[xmldna, 
            XMLElement["IUPACna", _, x_] -> 
              x, ∞]]]

To determine the position of a given motif in our DNA strand we use the StringPosition  function. Thus consider the motif "GGAACC". Its location is

StringPosition[myDNA, {"GGAACC"}]

We can improve  this  search by using the following function which gives the starting position of the DNA segments that match the enzyme's recognition sequence.

positionOfMotifs[str_String, motif__String] := 
      Map[{First[#][[1]], Last[#]} &, 
        Select[Map[{Flatten[StringPosition[str, #], 1], #} &, {motif}], 
          Length[First[#]] > 0 &]]

Thus for the enzyme "AccB2I"  we has the following list of recognition sequences

RecogSeq = 
      resSiteList[
        resSiteString[enzymeLookup["AccB2I", myEnzymeList]]]

We find that there are 3 locations in our DNA string that match  the enzymes recognition sequences

restrictionSites[data_] := 
      positionOfMotifs[data, Apply[Sequence, RecogSeq]]

restrictionSites[myDNA]

Thus the restriction enzyme AccB2I has 3 recognition sequences that have matches in our DNA. To determine the cut locations we use the information that the enzyme cuts between  C and Y bases : RGCGC^Y. Here is a function that does the task

findCutLocations[enzyme_, sites_] := 
      Module[{s, cutLocation, str}, 
        cutLocation[str_String] := 
          First[Flatten[StringPosition[str, "^"]]] - 
            2; 
s = enzymeLookup[enzyme, myEnzymeList]; 
        Map[{First[#] + cutLocation[s]} &, sites]]

The cut positions on this strand of DNA are

strand1Pos = 
      findCutLocations["AccB2I", restrictionSites[myDNA]]

Of course the enzyme will also cut the complement strand of our DNA sequence which can be found using the following function

\(reverseComplement[x_String] := 
      StringReplace[
        x, {"A" -> "T", "T" -> "A", "G" -> "C", "C" -> "G"}]\)

Thus for the complement DNA the cut locations are

myDNAc = reverseComplement[myDNA];

strand2Pos = 
      findCutLocations["AccB2I", restrictionSites[myDNAc]]

In the next section we illustrate how we can display our DNA sequence with the cut locations.

Displaying DNA sequences

Finally we would like to display graphically where the cuts will be made on our DNA sequence. The following method is based on a suggested by David Park of the Mathgroup. Our display function takes several arguments: The first argument defines the DNA sequence; the second argument the complement sequence, the third and fourth arguments are the positions where the cuts occur, the fifth argument is the location of the base pair where we want to start displaying; the sixth argument is the length of the row (i.e the number of base pairs displayed in a row); the last argument is the number of rows.

\(SequenceCartoon[strand1_, strand2_, cut1_, cut2_, \ start_, rowlength_, 
        rows_] := Module[{pos, seq11, seq21, subdata1, subdata2, subdata, \ 
          rowgrid1}, 
seq11 = 
          MapAt[ToLowerCase, Characters[strand1], cut1] /. 
            x_String?LowerCaseQ :> 
              StyleBox[ToUpperCase[x], 
                FontColor -> 
                  RGBColor[0, 0, 1]]; 
seq21 = 
          MapAt[ToLowerCase, Characters[strand2], cut2] /. 
            x_String?LowerCaseQ :> 
              StyleBox[ToUpperCase[x], 
                FontColor -> 
                  RGBColor[0, 0, 1]]; 
rowgrid1[
            subdata_]\  := 
\ \ \ \ \ \ GridBox[{subdata[\([1]\)], \ 
              MapAt[StyleBox[#, FontColor -> RGBColor[0, 0, 1]] &, \ 
                Table["|", \ {rowlength}], 
                Position[subdata[\([1]\)], 
                  StyleBox[x__]]], 
\ \ \ \ \ \ \ \ \ subdata[\([2]\)]}, \ 
            GridBaseline\  -> \ 
              Bottom]; 
\ \ subdata1\  = 
\ \ \ \ \ \ Take[
            seq11, \ {start, \ 
              Min[start\  + \ rows\ rowlength\  - \ 1, \ 
                Length[seq11]]}]; 
\ \ \ \ subdata1\  = \ 
          Partition[subdata1, \ rowlength, \ 
            rowlength, \ {1, \ 
              1}, \ " "]; 
subdata2 = 
\ \ \ \ \ \ Take[seq21, \ {start, \ 
              Min[start\  + \ rows\ rowlength\  - \ 1, \ 
                Length[seq21]]}]; 
\ \ \ \ subdata2\  = \ 
          Partition[subdata2, \ rowlength, \ 
            rowlength, \ {1, \ 1}, \ " "]; 
\ \ \ \ subdata\  = 
          Map[\ {rowgrid1[#]}\  &\ , 
            Transpose[{subdata1, subdata2}]]; 
\ \ \ \ DisplayForm[
          FrameBox[GridBox[subdata, \ RowLines\  -> \ True]]]]\)

The results are displayed below for the first 200 base pairs. The enzyme  "AccB2I" cuts to the right of the bases marked in blue:

SequenceCartoon[myDNA, myDNAc, strand1Pos, strand2Pos, 1, 20, 
      10]

Remarks

It is worth mentioning that the algorithm we used in this example for finding the motifs (the enzymes recognition sequences) is suitable if the number recognition sequences is moderate. If it is large, say in the thousands, then it makes sense to use a divide and conquer algorithm.

Protein 3-D Structure

In this example we are going to investigate how we can use Mathematica to manipulate protein data bank files. The data bank can be found at the following URL:
				http://www.rcsb.org/pdb/

We will a work with the file for the Sugar Binding Protein called "LIGAND-FREE CONGERIN I". For the purpose of this example I will be refer to the file as "1C1F.pdb". See website for this file.

We are going to make use of the ReadList function to import the data into Mathematica

myPDBfile = 
        ReadList[Experimental`FileBrowse[False], Word,  
          WordSeparators -> {"\r"}];

Here is what the first 10 lines of the file look like

Table[myPDBfile[[i]], {i, 1, 10}] // TableForm

The text at the beginning of each line defines the record type. Thus in the above output we have a HEADER, TITLE, COMPND and SOURCE record types. There may be multiple records of a given record type as the above output shows:  there are  3 records of the SOURCE record type, and 4 records of  the COMPD record type.

3-D Structure

The PDB file also contains the coordinates of the atoms in our protein. The record type that contains coordinate data is called ATOM. Let us extract all the atom coordinates from our data file

atomdata = 
        Select[myPDBfile, StringMatchQ[#, "ATOM*"] &];

This is a large file as we can see by calculating the length of the list

Length[atomdata]

Let us examine  a single entry

atomdata[[35]] // TableForm

Let us extract out  a subset of the atom data: the coordinates. In PDB files this data is in columns 31 through 54. Thus using StringTake gives

StringTake[atomdata[[35]], {31, 54}]

We use the following function to extract the data

s = 
      StringToStream[
        StringJoin[
          Table[StringTake[atomdata[[i]], {31, 54}], {i, 1, 
              Length[atomdata]}]]]

We use ReadList to read the data. Since we are dealing with numbers we specify that the each record is a sublist of 3 numbers separated by a Null  string.

atomCoord = 
        ReadList[s, {Number, Number, Number}, 
          RecordSeparators -> {" "}];

Here is an example of the coordinates for atom #2

atomCoord[[2]]

One of the things we can do with our data is plot the data as a 3D line plot using the function ScatterPlot. To use ScatterPlot we need to first load the package Graphics`Graphics3D`

<<Graphics`Graphics3D`

Here is a line plot of our protein

ScatterPlot3D[atomCoord, PlotJoined -> True, 
        Boxed -> False, Axes -> False, 
        PlotStyle -> RGBColor[0, 0, 1]];

We can use the RealTime3D function to interactively view the 3-D structure from different perspectives. Note: we lose the color capabilities

<< RealTime3D`

Here is the plot in 3D

ScatterPlot3D[atomCoord, PlotJoined -> True, 
        Boxed -> False, Axes -> False, 
        PlotStyle -> RGBColor[0, 0, 1]];

<< Default3D`

Concluding Comments/Road Ahead

Mathematica is an attractive programming language for teaching bioinformatic skills

-  Strong string manipulation/processing capabilities
	-  Suitable for web-based scripting (webMathematica,webservices)
	- Modules for common bioinformatic operations easily constructed

- Can handle complex projects involving large programs
	- Can readily handle complex calculations/graphics 
	- Works on all platforms, and has a user friendly interface

At UC Davis we plan to offer a bioinformatics class for our biochemical engineering majors and then extend the class to other majors.

Acknowledgments

Thanks to the Mathgroup for answering numerous queries I had about string manipulation in Mathematica. In particular, I would like to thank David Parks, Allan Hayes, Harmut Wolf and Daniel Lichtblau for material assistance.