Bioinformatics in the Wolfram Language

John B. Cassel, 
Wolfram|Alpha Scientific Content Group

Join the conversation! #WolframTechConf

ABSTRACT

Bioinformatics computes inferences about the integrated structure and functioning of biological components using obtainable data, such as sequences and micro-array concentrations.  The Wolfram Language’s existing computable data provides a rich starting point for bioinformatics, with both baseline data, such as genome and protein sequences, and also further information, such as functional annotation, protein domains, and secondary structure.  However, the Wolfram Language is not limited to these resources, and we’ll see how to search for and import further bioinformatics data.  From here, we’ll see how the Wolfram Language can process this data via sequence alignment, clustering,  motif discovery, and other strategies to understand component function and relationships.  Finally, we will see how the Wolfram Language allows us to demonstrate models of the discovered biological system dynamics.

Outline

Introduction

Accessing Built-in Computable Data

Gathering and Importing External Data

Analyzing Sequences

Analyzing and Modeling Interactions

Conclusion

Introduction

Defining Bioinformatics

The process of computing inferences 
about the integrated structure and functioning of biological components 
using obtainable data

Introduction

Purpose of this Talk

To provide a necessarily partial demonstration of
how existing Wolfram Language features facilitate bioinformatics.

Introduction

Limitations

This talk focuses on 
briefly covering basic applications that illustrate functionality 
rather than 
providing a comprehensive data-to-dynamics workflow

Introduction

Very Brief Biology Primer

Deoxyribonucleic acid (DNA) molecules are formed of two strands that wrap around each other, only binding together only if residues along their length match into complementary base pairs.  This correspondence allows the maintenance and duplication of heritable biological instruction.

ListAnimate@{,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,}

(Image from https://en.wikipedia.org/wiki/DNA#/media/File:ADN_animation.gif; public domain.)

Along these strands, triples of DNA base pairs form sequences in three letter blocks, called codons.  Through codons that indicate beginnings and endings, reading frames are delimited.  If regulatory regions around the reading frames are bound appropriately, these codons are transcribed from DNA into messenger RNA (mRNA).  Then, this mRNA delivers these sequences to ribosomes which translate them into the right amino acids for a new protein.  The chemical properties of those acids fold into different structures appropriate to different functions.

(Image from https://en.wikipedia.org/wiki/File:Gene_structure _prokaryote _ 2_annotated.svg; see link for usage details)

Built-in Computable Data

Overviews from Wolfram|Alpha

Starting a line with ==, or calling WolframAlpha, allows us to see an overview of available results through full Wolfram|Alpha output.

Genes

mus81 is a gene for a protein that, as part of a broader complex, repairs nicked Holliday junctions, a form of damage to a structure made during DNA recombination:

Partial Sequences

For a given sequence of base pairs, we can see which amino acids it composes and where it lives on the genome:

Proteins

In addition to the information about the mus81 gene, information about the mus81 protein is also available:

Metabolic Pathways

Metabolic pathways are the key cycles by which food and fuel is converted into the structures and energy of biological systems.  Here is the very widespread TCA cycle:

Built-in Computable Data

Exposure through Entities

The entity framework delivers a consistent experience by providing a uniform way to get the latest data in computable form.

Genes

As mentioned above, mus81 is the gene for a protein that, as part of a broader complex, repairs nicked Holliday junctions, a form of damage to a structure made during DNA recombination.

mus81Gene=%

Here we see those properties in computable form:

displayEntityValueDataSet[entity_]:=Module[{entityProperties=EntityProperties[entity]},
Dataset`$DatasetTargetRowCount=Length[entityProperties];
Dataset[Select[<|"Property"->#,"Value"->EntityValue[entity,#]|>&/@entityProperties,Head[#["Value"]]=!=Missing&]]
];
displayEntityValueDataSet[mus81Gene]

We can still retrieve the chromosomal context elided above:

EntityValue[mus81Gene,"Context"]

Proteins

We can also retrieve the associate protein entity for mus81:

mus81Protein=%

We see that we have a fair amount of data available for the mus81: protein, including , associated cellular components, biological processes, and molecular functions, and also domains with positions:

displayEntityValueDataSet[mus81Protein]

aconitase1 is an enzyme that catalyzes a step of the TCA cycle, and is thus quite important and well known.

aconitaseProtein=%

Given the prominence of aconitase1, more data is available to describe it than mus81.

primaryPDBUrl[protein_]:="http://www.rcsb.org/pdb/download/downloadFile.do?fileFormat=pdb&compression=NO&structureId="<>EntityValue[protein,"PrimaryPDBID"];displayEntityValueDataSet[aconitaseProtein]

Let’s have a look at the Robinson (ribbon) diagram:

aconitaseProtein[(EntityProperty["ribbon diagram", EntityProperty["Protein","StructureDiagram"], "EntityProperty[\"Protein\", \"StructureDiagram\"]"])]

Clearly, there’s a lot of secondary structure here, so let’s have a look at a larger version, in which we can see a number of α-helixes and β-sheets:

Import[EntityValue[aconitaseProtein,"PrimaryPDBURL"],"PDB",ImageSize->Large]

We can see where in the sequence these secondary structures occur:

aconitaseProtein[(EntityProperty["secondary structure rules", EntityProperty["Protein","SecondaryStructureRules"], "EntityProperty[\"Protein\", \"SecondaryStructureRules\"]"])]

SNPs

Single nucleotide polymorphism (SNPs) are variations of a particular genetic base between individuals, and may indicate traits, diseases, heritage, and evolutionary history.

rs1029=Interpreter["SNP"]["RS1029"]

displayEntityValueDataSet[rs1029]

Built-in Computable Data

Data Functions

To explore particular domains, built-in functions allow a different path for exploration.

GenomeData

A number of properties are available to describe genomes:

genomeDataProperties=GenomeData["Properties"]

Here is a different way to list the molecular functions of mus81.

GenomeData["MUS81","MolecularFunctions"]

Gathering all of the values into a Dataset is straightforward:

Dataset`$DatasetTargetRowCount=Length[genomeDataProperties];
Dataset[<|DeleteCases[
#->GenomeData["MUS81",#]&/@genomeDataProperties,
_->_Missing]|>]

GenomeLookup

Look up where along the human genome particular base-pair sequences can be found:

GenomeLookup["GCTATGCGATGAT"]

We can fetch a particular range of the sequence to confirm:

GenomeData[{{"Chromosome8",1},{82498080,82498092}}]

GenomeData[{{"Chromosome9",-1},{46880660,46880672}}]

ProteinData

Find a protein sequence by name:

ProteinData["MUS81"]

Other properties are available:

proteinDataProperties=ProteinData["Properties"]

List the properties available for this protein:

Dataset[<|DeleteCases[
#->ProteinData["MUS81",#]&/@proteinDataProperties,
_->_Missing]|>]

Gathering and Importing External Data

There are thousands of different bioinformatics sources online.  Though we’re going to only cover a few here, the sources here are significant, broadly applicable, and demonstrate similarities to how you’ll access data from other sources.

NCBI Entrez E-Utilities

Other Notable External Data Sources

Import File Formats

Gathering and Importing External Data

NCBI Entrez E-Utilities

The National Center for Biotechnology Information (NCBI) is a key source for curated bioinformatics information.  One way to access NCBI data is through the Entrez system.  The NCBI’s Entrez E-Utilities are an organized suite of interfaces for accessing NCBI Entrez resources in a uniform way.

Resource Functions

We’ve provided access to a number of these E-Utilities into resource functions to ease gathering this data. Provided services include:

EInfo:  gathering metadata about NCBI databases and the kinds of fields and links they contain.

eInfo= ResourceFunction["EInfoNCBI"]

EGQuery:  given a term, return the count of records for each NCBI database

egQuery=ResourceFunction["EGQueryNCBI"]

ESearch:  given queries of field content, gather IDs for the data

eSearch=ResourceFunction["ESearchNCBI"]

ESummary:  given IDs, return summaries of the corresponding data

eSummary=ResourceFunction["ESummaryNCBI"]

EFetch:  given IDs, return the full corresponding data

eFetch=ResourceFunction["EFetchNCBI"]

ESpell:  given a word, return a spelling correction to that word if necessary

eSpell=ResourceFunction["ESpellNCBI"]

Let’s use these tools to look at retrieving data for MUS81.

EInfo

EInfo is a service for gathering metadata about NCBI databases and the kinds of fields and links they contain.

See the list of databases:

eInfo["databases"]

Now we can get a dataset of fields in the “gene” database.  Later, we’ll use the GENE field.

eInfo["gene","fields"]

EGQuery

The EGQuery service will return the count of records containing a given term for each NCBI database.  Here, we find that MUS81 has numerous occurrences in the “gene” database:

egQuery["MUS81"]

ESearch

ESearch provides IDs for the data that match the fields of the given queries.  Here we are looking to records related to the A2M gene in the GENE field of the gene database:

mus81Ids=eSearch["gene","GENE","MUS81"]

ESummary

ESummary provides summaries of the data corresponding to the provided ids.  Let’s look at the summary of one of the results provided in the last query:

eSummary["gene",mus81Ids[[2]]]

EFetch

EFetch will return the full contents for the data associated with the given id.  
For the valid return types for each fetch, see Table 1 of the “The E-utilities In-Depth: Parameters, Syntax and More” (https://www.ncbi.nlm.nih.gov/books/NBK25499/) .

efetchMus81GeneResponse=eFetch["gene", mus81Ids[[2]], "text", "gene_table"];

Let’s only see the first 30 lines of the result:

Column[Flatten[efetchMus81GeneResponse[[1;;30]]]]

ESpell

On a different topic, ESpell provides a nice spell-checking service for biological terms.  Here, we see that “mitocondrial” should be spelled “mitochondrial”:

eSpell["mitocondrial"]

Of course, “mitochondrial” itself is fine:

eSpell["mitochondrial"]

For more information

An additional slide has been provided at the end of this presentation to demonstrate accessing these services at a lower-level than resource functions.  For more about the E-Utilities services themselves, see the “The E-utilities In-Depth: Parameters, Syntax and More” (https://www.ncbi.nlm.nih.gov/books/NBK25499/) section of the Entrez Programming Utilities Help.  For further examples of using Entrez resources in the Wolfram Language, see George Mias’s book “Mathematica for Bioinformatics:  A Wolfram Language Approach to Omics”.

Gathering and Importing External Data

Other External Data Sources

There are a great number of other bioinformatics sources out there.  We’re only going to look at two of the biggest and longest-lasting.

The RCSB Protein Data Bank is an import resource for information about proteins:

Import["http://www.rcsb.org/pdb/download/downloadFile.do?fileFormat=pdb&compression=NO&structureId=1tf6","PDB",ImageSize-> Small]

The European Bioinformatics Institute (EBI) is the EU equivalent to NCBI.  It has extensive resources including its Ensembl and UniProt projects which provide nice REST interfaces:

Dataset/@URLExecute["https://www.ebi.ac.uk/proteins/api/proteins/interaction/Q01705",{"RequestMethod"->"GET"},"RawJSON"]

Gathering and Importing External Data

Import File Formats

GenBank

GenBank is the standard plaintext format for storing and exchanging annotated DNA sequences.

Here is an example of GenBank importing where everything goes right:

Import["ExampleData/rCRS2006.gbk", "Elements"]

Dataset[<|#->Import["ExampleData/rCRS2006.gbk", #]&/@Import["ExampleData/rCRS2006.gbk", "Elements"]|>]

Even when a file isn’t formatted entirely to the standard, you’ll still be able to do the analysis you need.  This record might look fine:

mus81NuccoreIds=eSearch["nuccore","GENE","MUS81"];
mus81GenBankText=eFetch["nuccore",mus81NuccoreIds[[1]],"text","gb"];
StringTake[mus81GenBankText,1000]

When we import it, we see that something is missing:

ImportString[mus81GenBankText,{"GenBank","Elements"}]

It turns out the GenBankID is missing:

ImportString[mus81GenBankText,{"GenBank","GenBankID"}]

However, one can access the rest of the data successfully:

Dataset[<|ImportString[mus81GenBankText,{"GenBank","Features"}]|>]

FASTA

FASTA is the standard format for storing and exchanging DNA and protein sequences.

mus81FastaText=eFetch["nuccore",mus81NuccoreIds[[12]],"text","fasta_cds_na"];

ImportString[mus81FastaText,{"FASTA", "Header"}]

ImportString[mus81FastaText,{"FASTA", "Sequence"}]

PDB

PDB is commonly used for storing structure information for large biological molecules such as proteins and nucleic acids:

pdb1tf6Elements=Import["http://www.rcsb.org/pdb/download/downloadFile.do?fileFormat=pdb&compression=NO&structureId=1tf6",{"PDB","Elements"}]

pdb1tf6Results=Import["http://www.rcsb.org/pdb/download/downloadFile.do?fileFormat=pdb&compression=NO&structureId=1tf6",{"PDB",#}]&/@pdb1tf6Elements;

Dataset[<|MapThread[Rule,{pdb1tf6Elements,pdb1tf6Results}]|>]

Other Supported Bioinformatics File Formats

FASTQ:  an extended version of FASTA that includes encoded probability estimates for if each base is correct.

NEXUS:  a commonly used format for storage and exchange of phylogenetic data capable of representing DNA and protein sequences, taxa distances, alignment scores, and phylogenetic trees.

AgilentMicroarray:  contains labeled data, process information, and overall statistics for Agilent microarray results.

Affymetrix:  contains labeled data, process information, and confidence statistics for Affymetrix microarray results.

SFF:  contains DNA sequences, flow values, and quality estimates from pyrosequencing runs, along with associated process information.

Sequence Analysis

Sequence Alignment and Clustering

Let’s first have a toy example of sequence alignment to see how it works:

SequenceAlignment["ATAC","ATGGC"]

EME1 Gene

The EME1 gene has homologs across many species.  Wolfram|Alpha has EME1 sequence data for humans, chimpanzees, mice, and rats.

WolframAlpha["EME1 gene sequence"]

WolframAlpha["EME1 gene mouse sequence"]

Gathering Gene Data

Gather the gene sequences for each of the organisms we’ve selected:

organismLabels={"Mouse","Rat","Chimpanzee","Human"};

eme1OrganismGenes=Interpreter["Gene"]["EME1 gene"<>#]&/@organismLabels

eme1GeneSequences=EntityValue[#,"ReferenceSequence"]&/@eme1OrganismGenes;

StringTake[#,5]&/@eme1GeneSequences

Add an indexing function for later convenience:

organismNameToIndex=Association@MapThread[Rule,{organismLabels,Range[Length[organismLabels]]}]

Defining a Distance Metric

Smith-Waterman similarity is the standard local similarity value score in sequence alignment.  The potential size of the score depends on the size of the sequence, so we will need to normalize it against the potential maximum of perfect alignment.

SmithWatermanSimRatio[seq1_,seq2_]:=SmithWatermanSimilarity[seq1,seq2]/SmithWatermanSimilarity[seq1,seq1]

Clustering expects similarity to be a distance: i.e. smaller the more similar the sequences are:

SmithWatermanRatioDistance[seq1_,seq2_]:=1-SmithWatermanSimRatio[seq1,seq2]

Sequence Alignment

From the sequence alignment, it appears that the EME1 gene sequence for humans and chimpanzees are much more similar than those for rats and mice.

eme1MouseRatAlignment=SequenceAlignment[eme1GeneSequences[[organismNameToIndex["Mouse"]]],eme1GeneSequences[[organismNameToIndex["Rat"]]]];
eme1MouseRatAlignment[[1;;20]]

eme1ChimpHumanAlignment=SequenceAlignment[eme1GeneSequences[[organismNameToIndex["Chimpanzee"]]],eme1GeneSequences[[organismNameToIndex["Human"]]]];
eme1ChimpHumanAlignment[[1;;4]]

Here’s a small tool to explore these alignments more interactively:

exploreSequenceAlignment[seqAlign_,shownSize_:5]:=With[{seqLength=Length[seqAlign]},
Manipulate[Grid[{If[ListQ[#],Column[#,Frame->All,Spacings->1],#]&/@seqAlign[[x;;x+shownSize]]},Frame->All],{x,1,seqLength-shownSize,1}]
]

Using this tool, we can see that the start of the sequence appears well aligned, but by the end, there is less agreement.

exploreSequenceAlignment[eme1ChimpHumanAlignment]

Sequence Clustering

First, we’ll calculate all similarity ratios:

allGeneEME1SimilarityRatios=Apply[SmithWatermanSimRatio,Partition[Tuples[{eme1GeneSequences,eme1GeneSequences}],Length[eme1GeneSequences]],{2}];

Unsurprisingly, we see two clear clusters for similarity: mouse/rat and human/chimpanzee:

axisTicks=MapThread[List,{Range[Length[organismLabels]],organismLabels}];
ArrayPlot[allGeneEME1SimilarityRatios,PlotLabel->"EME1 Sequence Similarity",PlotLegends->Automatic,FrameTicks->{axisTicks,axisTicks}]

styledOrganismLabels=Style[#,Bold,Black]&/@organismLabels;
BarChart3D[allGeneEME1SimilarityRatios,ChartLayout->"Grid",ChartLabels->{styledOrganismLabels,styledOrganismLabels},PlotLabel->"EME1 Sequence Similarity",ViewPoint->{0.7 π,-0.45 π,0.35π}]

By taking the dendrogram, we can see the relative tightness of these clusters.

Dendrogram[MapThread[Rule,{eme1GeneSequences,organismLabels}],DistanceFunction->SmithWatermanRatioDistance,PlotLabel->"EME1 Sequence Similarity",ClusterDissimilarityFunction->"Complete"]

Sequence Analysis

Finding Sequence Motifs

Unlike pairwise sequence alignment, the Wolfram Language does not have primitives for finding motifs in multiple sequences.  However, it does have sequence manipulation functions that allow motifs to be extracted without too much trouble.

Mus81 protein sequences

Here I load some FASTA protein sequences referred to in Enzlin, J. H., & Schärer, O. D. (2002). The active site of the DNA repair endonuclease XPF–ERCC1 forms a highly conserved nuclease motif. The EMBO Journal, 21(8), 2045–2053. http://doi.org/10.1093/emboj/21.8.2045

proteinSequenceByAccession[accession_String]:=Module[{fetchResult,index},
fetchResult=eFetch["protein",First[eSearch["protein","PACC",accession]],"text","fasta"];
If[ListQ[fetchResult],
fetchResult=StringRiffle[Flatten[fetchResult, 1], "\n"];
];
index=FirstPosition[First@*StringSplit/@ImportString[fetchResult,{"FASTA","Header"}],accession][[1]];
ImportString[fetchResult,{"FASTA","Sequence"}][[index]]
];

mus81AccessionRules={"Yeast"->"NP_010674.3","Fruitfly"->"AAF45571.1","Mouse"->"AAL28066.1","Human"->"AAL28065.1"};
aminoAcidMus81Sequences=proteinSequenceByAccession[Last[#]]&/@mus81AccessionRules;
{yeastAminoAcidMus81Seq,fruitflyAminoAcidMus81Seq,mouseAminoAcidMus81Sequence,humanAminoAcidMus81Seq}=aminoAcidMus81Sequences

Calculating Sequence Motifs

First, we find the longest non-contiguous sequence in common between all of these protein sequences:

mus81CommonAminoAcidSeq=Fold[LongestCommonSequence,First@aminoAcidMus81Sequences,Rest@aminoAcidMus81Sequences];

Next, we find local alignments between that common subsequence and each of the protein sequences under study:

mus81CommonAminoAcidAlignments=SequenceAlignment[#,mus81CommonAminoAcidSeq,Method->"Local"]&/@aminoAcidMus81Sequences;

For each of these alignments, we’ll replace single-character disagreements with an underscore and longer disagreements with a plus:

localMotifs=(If[ListQ[#],If[StringLength[#[[1]]]===1,"_","+"],#]&/@#)&/@mus81CommonAminoAcidAlignments;

Grid[localMotifs]

Given the above align, we’ll take the minimum matched sequence for each group, splitting + sequences into separate patterns.

StringSplit[StringJoin[First@*MinimalBy[StringLength]/@Transpose[localMotifs]],"+"]

Given the incredible simplicity of this approach, we still find most of the well-conserved IVERKR_DDL motif from  Enzlin, J. H., & Schärer, O. D. (2002):

Using External Services

A wide variety of external services are available for multiple-sequence alignment at EMBL-EBI:  https://www.ebi.ac.uk/seqdb/confluence/pages/viewpage.action?pageId=54646455

Let’s take a look at an example using Clustal Omega, a standard middle-of-the-road choice.  First, we’ll see the parameters we need to supply.

clustaloBaseURL="https://www.ebi.ac.uk/Tools/services/rest/clustalo";
fastaSeqSubmission=StringRiffle[Flatten[MapThread[List,{StringJoin[">",Last[#]]&/@mus81AccessionRules,aminoAcidMus81Sequences}]],"\n"];

URLExecute[clustaloBaseURL<>"/parameters",{"RequestMethod"->"GET"}]

Now, let’s submit our sequences to run an alignment job.

submissionResult=URLRead[HTTPRequest[clustaloBaseURL<>"/run",<|Method->"POST","Body"->{"sequence"->fastaSeqSubmission,"email"->$MyEmail}|>],{"StatusCode","Headers","Body"}]

jobID="Body"/.submissionResult

Submission status

URLExecute[clustaloBaseURL<>"/status/"<>jobID,{"RequestMethod"->"GET"}]

Available submission result types:

resultTypes=URLExecute[clustaloBaseURL<>"/resulttypes/"<>jobID,{"RequestMethod"->"GET"}];
resultTypesResponseXML=ImportString[resultTypes,"XML","ReadDTD"->False];
Grid[Function[{result},Flatten[xmlTagExtractor[#,result]&/@{"identifier","description"}]]/@Cases[resultTypesResponseXML,XMLElement["type",___],Infinity],Frame->All]

Submission results

clustalOmegaResult=URLExecute[clustaloBaseURL<>"/result/"<>jobID<>"/aln-clustal_num",{"RequestMethod"->"GET"}]

Rather than trying to interpret this wall of data, let’s extract some motifs from the alignment results:

stringMask[givenString_String,maskString_String]:=Module[{rawChars,positions},
rawChars=Table["_",{StringLength[givenString]}];
Scan[(rawChars[[#]]=StringPart[givenString,#])&/@First/@StringPosition[maskString,"*"]];
StringJoin[rawChars]
];
gatherClusterOmegaMotifs[clustalOmegaResult_]:=Module[{splitLines,currentResultLines={},resultLines={}, currentLine,positions,rawMotifs},
splitLines=Append[StringSplit[clustalOmegaResult,"\n"][[2;;-1]],""];
If[StringLength[#]>0,
currentLine=StringReplace[StringDrop[#,17],
"\t"~~DigitCharacter..->""];
currentResultLines=Append[currentResultLines,currentLine],
If[Length[currentResultLines]>0,
If[Length[resultLines]>0,
resultLines=MapThread[StringJoin,{resultLines,currentResultLines}],
resultLines=currentResultLines
];
currentResultLines={}
];
]&/@splitLines;
positions=StringPosition[Last[resultLines],Except[WhitespaceCharacter]..,Overlaps->False];
rawMotifs=Select[Transpose[StringTake[#,positions]&/@resultLines],StringContainsQ["*"]@*Last];
Select[StringTrim[stringMask[First@#,Last@#],"_"]&/@rawMotifs,StringLength[#]>3&]
];
gatherClusterOmegaMotifs[clustalOmegaResult]

Here we see that we also find the IVERKR_DDL motif, this time including the trailing L.

Analyzing and Modeling Interactions

Bioinformatics is not only about understanding the functionality of particular biological structures, but also understanding the interaction between structures.  Here, we discover the insight that comes from discovering frequent structures of interaction and investigating the resulting dynamics.  This work is based on Shen-Orr, S. S., Milo, R., Mangan, S., & Alon, U. (2002). Network motifs in the transcriptional regulation network of Escherichia coli. Nature genetics, 31(1), 64 and related resources.

E. Coli Regulatory Interactions Overview

Shown below is a graph of the gene regulation interactions in E. Coli.  While this appears complicated, we will be able to get some insights about why it has these structures by discovering common substructures, even when these substructures are very small.

ecoliRegulatoryDataset=ResourceData["Transcriptional Regulation Network of Escherichia coli"];

createRegulatoryGraph[regulatoryInteractions_Dataset,options___]:=
createRegulatoryGraph[{#["RegulatingTranscriptionFactor"],#["RegulatedOperon"],#["RegulationType"]}&/@Normal[regulatoryInteractions],options];
createRegulatoryGraph[regulatoryInteractions_List,options___]:=
Graph[Normal[DirectedEdge[#[[1]],#[[2]]]&/@regulatoryInteractions],
EdgeStyle->Normal[(#[[1]]->#[[2]])->#[[3]]/.{"activate"->Red,"inhibit"->Blue,"dual"->Purple}&/@regulatoryInteractions],options];
regulatoryGraph=createRegulatoryGraph[ecoliRegulatoryDataset]

Analyzing and Modeling Interactions

Auto-regulation as an Interaction Motif

Looking at the previous graph, it seems as though a node providing negative feedback to itself occurs relatively commonly:

createRegulatoryGraph[{{"1","1","2"}}]

Self-regulation versus Erdos-Renyi random likelihood

We expect the number of self-edges to be the probability of a given edge  ((E/(N^2)))times the number of opportunities to make that edge (N), or  ((E/N)).

Length@Union[Normal[ecoliRegulatoryDataset[All,"RegulatedOperon"]]]

nodeCount=Length@Union[Normal[ecoliRegulatoryDataset[All,"RegulatedOperon"]]]

edgeCount=Length[ecoliRegulatoryDataset]

We expect the number of self-edges to be the probability of a given edge  ((E/(N^2)))times the number of opportunities to make that edge (N), or  ((E/N)).

(Subscript[n, self,random])=N[edgeCount/nodeCount]

Due to results not shown here, as the number of nodes and edges grow large, we can approximate the standard deviation of an expected edge as the square root of the expected self-edge count.

(Subscript[σ, self,random])=N[(Sqrt[edgeCount/nodeCount])]

The actual number of self-regulators in the graph is 59:

selfRegulators=Select[ecoliRegulatoryInteractions,#[[1]]===#[[2]]&];

(Subscript[n, self])=Length[selfRegulators]

Therefore, the actual number is more than 46 standard deviations away from what we might expect randomly:

(((Subscript[n, self])-(Subscript[n, self,random]))/(Subscript[σ, self,random]))

Modeling Auto-regulation Dynamics

First, let’s look at the case without self-regulation.

Assume that mRNA production and degradation are in quasi-steady state, and that transcription factor degradation is negligible at the considered timescale.

Given α, the maximal production rate;  and K,  the half-saturation concentration, we can represent   p, the concentration of transcription factor P, as (d/(dt))p(t)=α (1/(1+p(t)/K))

DSolve[p'[t]==(α(1/(1+(p[t]/K)))),p[t],t]

Assuming initial concentration (Subscript[c, 1])=0,  and that α=1 and K=1:

Plot[-1+(Sqrt[1+2 t ]),{t,0.0,4.0}]

Now, let’s look at the auto-regulation case.  Here, we add a negative regulation with a binding factor of σ, leading to (d/(dt))p(t)=α (1/(1+p(t)/K))-σ p(t)

DSolve[p'[t]==(α(1/(1+(p[t]/K))))-(σ p[t]),p[t],t]

Let us arbitrarily set K=1:

DSolve[p'[t]==(α(1/(1+(p[t]))))-(σ p[t]),p[t],t]

To get an equivalent reaction by a particular time,  we need to compensate for the negative regulation with a higher initial production rate.  We see that  σ=0.31 allows the maximum production rate to double while still reaching an equivalent production.

Manipulate[Plot[{-1+(Sqrt[1+2  t ]),InverseFunction[((ArcTan[(((Sqrt[σ]) (1+2 #1))/(Sqrt[-4 α-σ]))])/((Sqrt[-4 α-σ]) (Sqrt[σ])))+((Log[α-σ #1 (1+#1)])/(2 σ))&][-t+1]},{t,0,4}],{α,1.75,2.25},{σ,0.25,0.35}]

This faster response time allows for faster adaptation to new conditions.  Though not demonstrated here, the presence of feedback also allows for a greater robustness to noise in concentration of other activating factors.

Analyzing and Modeling Interactions

Three-node Interaction Motifs

The next simplest possible substructure is if three nodes are connected in an interesting way.  This section looks to discover interesting three-node connection motifs, check if they are significant, and if so determine their behavior.

Finding three-node interaction motif candidates

First, we isolate all of the subgraphs of three connected nodes:

(* remove self-edges and activation *)
nonSelfRegulators={#["RegulatingTranscriptionFactor"],#["RegulatedOperon"]}&/@Normal[Select[ecoliRegulatoryDataset[All,{"RegulatingTranscriptionFactor","RegulatedOperon"}],#["RegulatingTranscriptionFactor"]=!=#["RegulatedOperon"]&]];
(*find two edges that share one node *)
threeVertexEdgePairs=Union[Flatten[Function[{firstEdge},
Sort[{firstEdge,#}]&
/@Select[nonSelfRegulators,And[firstEdge=!=#,AnyTrue[SameQ@@@(Tuples[{firstEdge,#}]),Identity]]&]
]/@nonSelfRegulators,1]];
(* get only the three vertexes *)
connectedThreeVertexes=Union@*Flatten/@threeVertexEdgePairs;
(* find the edges connecting each of those three vertexes *)
connectedTripleGroups=Function[{connectedThree},Select[nonSelfRegulators,ContainsAll[connectedThree,#]&]]/@connectedThreeVertexes;
exampleGraphs=Graph@(DirectedEdge@@@#)&/@connectedTripleGroups[[1;;5]]

Find potential motifs by gathering each graph into isomorphic sets:

(* see for each graph, see if it is isomorphic to a representative of an existing group; if so, add it to that group; otherwise, start a new group *)
findIsomorphicGroups[groupsOfEdgePairs_]:=Module[
{candidateGraph,isomorphicGroups=<||>,isomorphicGraph},
Do[
candidateGraph=Graph@(DirectedEdge@@@edgePairGroup);
isomorphicGraph=SelectFirst[Keys[isomorphicGroups],IsomorphicGraphQ[candidateGraph,#]&];
If[Head[isomorphicGraph]===Missing,
isomorphicGroups[candidateGraph]={edgePairGroup},
isomorphicGroups[isomorphicGraph]=Append[isomorphicGroups[isomorphicGraph],edgePairGroup]
],
{edgePairGroup,groupsOfEdgePairs}
];
isomorphicGroups
];
isomorphicGroups=findIsomorphicGroups[connectedTripleGroups];
Length/@isomorphicGroups

The only interaction we find of more than two edges per three nodes is a “feed-forward” regulatory arrangement.  Are these arrangements significant?

Feed-forward motif likelihood versus Erdos-Renyi random likelihood

We’ve found 126 feed-forward subgraphs.

->126

(Subscript[n, feedforward])=126

The probability of a random edge is the number of edges divided among all pairs.

(Subscript[p, edge])=(edgeCount/(nodeCount^2))

There is only arrangement (a) of a feed-forward graph in which that feed-forward arrangement is maintained:

feedForwardIsomorphicPermulations=1

Given N nodes and n nodes to choose, there are ~(N^n)possible selections, and the odds of choosing the right placement for g edges is  ((Subscript[p, edge])^g).  Therefore, the 
(Subscript[n, arrangement,random]) ≈ a (N^n) ((Subscript[p, edge])^g)

(Subscript[n, feedforward,random])=N[feedForwardIsomorphicPermulations ((nodeCount^3) ((Subscript[p, edge])^3))]

As before, we can approximate the standard deviation of an expected motif as the square root of the expected count:

(Subscript[σ, feedforward,random])=(Sqrt[(Subscript[n, feedforward,random])])

Given this, we have over 64 standard deviations more feed-forward subgraphs than would be expected randomly:

(((Subscript[n, feedforward])-(Subscript[n, feedforward,random]))/(Subscript[σ, feedforward,random]))

Dynamic motifs among feed-forward loops

First, we’ll great a way to map all of the edges of a given subgraph together with common names:

mapTripleToCommonPoints[ffEdgeGroup_]:=Module[{tallyFirst,tallyLast},
tallyFirst=Tally[First/@ffEdgeGroup];
tallyLast=Tally[Last/@ffEdgeGroup];
<|SelectFirst[tallyFirst,Last[#]===2&][[1]]->"z",SelectFirst[tallyFirst,Last[#]===1&][[1]]->"y",SelectFirst[tallyLast,Last[#]===2&][[1]]->"x"|>
];
feedForwardEdgeGroups=isomorphicGroups[SelectFirst[Keys[isomorphicGroups],IsomorphicGraphQ[,#]&]];
mapTripleToCommonPoints[feedForwardEdgeGroups[[1]]]

Then, we’ll apply that transformation to all of the feedforward subgraphs, and tally up each subgraph, which is now only distinguished by the types of edge interaction:

unlabeledToLabeledInteraction=<|Function[{rowAssoc},Rule[rowAssoc/@{"RegulatingTranscriptionFactor","RegulatedOperon"},rowAssoc/@{"RegulatingTranscriptionFactor","RegulatedOperon","RegulationType"}]]/@Normal[ecoliRegulatoryDataset]|>;Grid[Transpose[{createRegulatoryGraph[#[[1]](*,VertexLabels->"Name"*)],#[[2]]}&/@Tally[Module[{groupMap},
groupMap=mapTripleToCommonPoints[#];
Sort[(Lookup[groupMap,#,#]&/@unlabeledToLabeledInteraction[#])&/@#]
]&/@feedForwardEdgeGroups]],Frame->All]

We see that the most frequent feed-forward loop motif is coherent (all paths effect child nodes in the same way) and consists of all activations.

Feed-forward coherent loop dynamics

As before, assume that mRNA production and degradation are in quasi-steady state, and that transcription factor degradation is negligible at the considered timescale.

Given (Subscript[α, x]), the maximal production rate;  and (Subscript[K, x]),  the half-saturation concentration, we can represent x, the concentration of transcription factor X of regulatory node x, as (d/(dt))x(t)=α (1/(1+x(t)/K))

(Subscript[K, xy])=(1/2);
(Subscript[K, yz])=((1;)/2)

DSolve[x'[t]==α/(1+(x[t]/k)),x[t],t,Assumptions->And[t≥0,α>0,k>0]]

Let us assume that the initial concentration is zero,  and that (Subscript[α, x])=1 and (Subscript[K, x])=0.5:

xSolved[t_,α_,k_]:=-k+(Sqrt[(k^2)+2 k t α])

xSolved[t_]:=xSolved[t,1,(1/2)]

Plot[xSolved[t],{t,0,2}]

Let’s make similar assumptions as above, and observe that production of Y is dependent upon X for activation.

Given (Subscript[α, y]), the maximal production rate;  and (Subscript[K, y]),  the half-saturation concentration, we can represent y, the concentration of transcription factor Y of regulatory node y, as (d/(dt))y(t)=(Subscript[α, y]) ((x(t)/(Subscript[K, y]))/(1+x(t)/(Subscript[K, y])))

DSolve[y'[t]==α ((x[t]/ky)/(1+(x[t]/ky))),y[t],t]

Let us assume that the initial concentration is zero, and that (Subscript[α, y])=1 and (Subscript[K, y])=(1/2):

DSolve[y'[t]==((xSolved[t]/((1/2)))/(1+(xSolved[t]/((1/2))))),y[t],t]

ySolved[t_]:=1.` t-(0.5) (Sqrt[1+4` t])+0.5

Plot[{ySolved[t]},{t,0,2},PlotLabels->{"y"}]

Let’s make similar assumptions as above, and observe that the production of Z is dependent upon X and Y for activation.

Given (Subscript[α, z]), the maximal production rate;  and (Subscript[K, yz]),  the half-saturation concentration, we can represent y, the concentration of transcription factor Y of regulatory node y, as (d/(dt))z(t)=(Subscript[α, z]) ((x(t)/(Subscript[K, xz]) y(t)/(Subscript[K, yz]))/(1+x(t)/(Subscript[K, xz])+y(t)/(Subscript[K, yz])))

DSolve[z'[t]==α ((((x[t]/kxz) (y[t]/kyz)))/(1+(x[t]/kxz)+(y[t]/kyz)+((x[t]/kxz) (y[t]/kyz)))),z[t],t]

Let’s make similar assumptions as above, and observe that the production of Z is dependent upon X and Y for activation.

Let us assume that the initial concentration is zero, and that (Subscript[α, z])=1, (Subscript[K, xz])=(1/2), and (Subscript[K, yz])=(1/2):

kxz=0.5;
kyz=0.5;
αz=1.0;

DSolve[z'[t]==αz ((((xSolved[t]/kxz) (ySolved[t]/kyz)))/(1+(xSolved[t]/kxz)+(ySolved[t]/kyz)+((xSolved[t]/kxz) (ySolved[t]/kyz)))),z[t],t]

DSolve[{z'[t]== ((((xSolved[t]/((1/2))) (ySolved[t]/((1/2)))))/(1+(xSolved[t]/((1/2)))+(ySolved[t]/((1/2)))+((xSolved[t]/((1/2))) (ySolved[t]/((1/2)))))),z[0]==0},z[t],t,Assumptions->t≥0]

zSolved[t_]:=-1.` (-0.32671320486001376`-1.` t+1.` (Sqrt[0.25` +1.` t])+(0.` +0.5 I) ArcTan[(0.9428090415820632` +0.6666666666666666` I) (Sqrt[0.25` +1.` t])]-0.4999999999999999` ArcTanh[(0.6666666666666666` +0.9428090415820632` I) (Sqrt[0.25` +1.` t])]+0.25` Log[1.5` +1.9999999999999996` t+1.9999999999999996` (t^2)])

With all of these results in place, we can see that growth in the concentration of Z is quite delayed.  This has the effect that noise in the production of X leads to very little Z being produced, and that Z will be robust to temporary bursts of X production.

Plot[{xSolved[t],ySolved[t],zSolved[t]},{t,0,1},PlotLabels->{"x","y","z"}]

Instead of Z requiring both X and Y for production, either X or Y might have been sufficient for activation.  In this circumstance, Z will react immediately to concentrations of X, but will be robust to temporary gaps in the production of X.

What wasn’t covered

SystemModeler integration and related support for systems biology

Biochemistry data

RandomChoice[ChemicalData["Biomolecules"],20]

Hidden Markov Models

Any number of other topics...

Conclusion

Wolfram Language support for bioinformatics is broad,

The Wolfram Language has suitable constructs for all phases of bioinformatics projects:
data gathering, analysis, modeling, visualization, etc.

but not comprehensive

There are many ways the Wolfram Language can be extended to support bioinformatics:
further data sources, more import formats, other analysis approaches, 
deeper modeling support, more visualizations, further linkage, and further support materials

and instead open to new directions

We are always interested in feedback; 
there are many directions for Wolfram Language bioinformatics development, 
and which to pursue depends on which our collaborators can put to use

Thank you!

Please feel free to email jcassel@wolfram.com to express your interests in bioinformatics and other biology topics

Acknowledgements

Special thanks to fellow members of the W|A scientific content team

Appendix:  Additional NCBI Service Access Details

This appendix provides the internals to how we accomplished accessing NCBI services, which might be useful for accessing other services not demonstrated in this talk, such as ELink and EPost.

EInfo Service Details

EInfo supports JSON results, which is convenient as those results are easily converted into nested associations.

See the list of databases:

eInfoDatabaseList=URLExecute["https://eutils.ncbi.nlm.nih.gov/entrez/eutils/einfo.fcgi",{"RequestMethod"->"POST","retmode"->"json"},"RawJSON"];

eInfoDatabaseList["einforesult"]["dblist"]

Now we can get the list of fields in the “gene” database.  Later, we’ll use the GENE field.

geneDBInfo=URLExecute["https://eutils.ncbi.nlm.nih.gov/entrez/eutils/einfo.fcgi",{"RequestMethod"->"POST","retmode"->"json","db"->"gene"},"RawJSON"]["einforesult"];

Dataset[geneDBInfo["dbinfo"]["fieldlist"]]

EGQuery Service Details

The EGQuery service returns the count of records containing a given term for each NCBI database.  Unfortunately, we see JSON isn’t supported here.

URLExecute["https://eutils.ncbi.nlm.nih.gov/entrez/eutils/egquery.fcgi",{"RequestMethod"->"POST","retmode"->"json","term"->"MUS81"},"RawJSON"]

We can use XML instead; it’ll require just a bit more extraction.  First we can define a nice little extraction utility for extracting the content within a particular tag:

(RefLink[Needs, paclet:ref/Needs])["XML`"];
xmlTagExtractor[tag_,text_]:=Cases[text,(XMLElement[tag,{},{content_}]:>content),Infinity]

Then we’ll map that utility over each result item to get the fields we want, and find that MUS81 has numerous occurrences in the “gene” database:

mus81EGQuery=URLExecute["https://eutils.ncbi.nlm.nih.gov/entrez/eutils/egquery.fcgi",{"RequestMethod"->"POST","retmode"->"json","term"->"MUS81"}];
mus81EGQueryResponseXML=ImportString[mus81EGQuery,"XML","ReadDTD"->False];
counts=Grid[(Sort[Function[{result},Flatten[xmlTagExtractor[#,result]&/@{"DbName","MenuName","Count"}]]/@Cases[mus81EGQueryResponseXML,XMLElement["ResultItem",___],Infinity]]),
Frame->All]

ESearch Service Details

ESearch provides IDs for the data that match the fields of the given queries.  Here we are looking to records related to the A2M gene in the GENE field of the gene database:

eSearchResult=URLExecute["https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi",{"RequestMethod"->"POST","retmode"->"json","db"->"gene","term"->"MUS81[GENE]"},"RawJSON"]["esearchresult"];

Dataset[eSearchResult]

mus81Ids=eSearchResult["idlist"]

ESummary Service Details

ESummary provides summaries of the data corresponding to the provided ids.  Let’s look at the summary of one of the results provided in the last query:

eSummaryMus81=URLExecute["https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi",{"RequestMethod"->"POST","retmode"->"json","db"->"gene","id"->mus81Ids[[1]]},"RawJSON"];

Dataset[eSummaryMus81["result"][mus81Ids[[1]]]]

EFetch Service Details

EFetch will return the full contents for the data associated with the given id.  
For the valid return types for each fetch, see Table 1 of the “The E-utilities In-Depth: Parameters, Syntax and More” (https://www.ncbi.nlm.nih.gov/books/NBK25499/) .

efetchMus81GeneResponseFlat=URLExecute["https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi",{"RequestMethod"->"POST","retmode"->"text","rettype"->"gene_table","db"->"gene","id"->mus81Ids[[1]],"complexity"->"0"}];

Column[Flatten[efetchMus81GeneResponseFlat[[1;;30]]]]

ESpell Service Details

On a different topic, ESpell provides a nice spell-checking service for biological terms.  Here, we see that “mitocondrial” should be spelled “mitochondrial”:

spellCheckResult=URLExecute["https://eutils.ncbi.nlm.nih.gov/entrez/eutils/espell.fcgi",{"RequestMethod"->"POST","db"->"pubmed","term"->"mitocondrial"}];
responseXML=ImportString[spellCheckResult,"XML","ReadDTD"->False];
Flatten[xmlTagExtractor[#,responseXML]&/@{"Query","CorrectedQuery"},1]

For more details

For more about the E-Utilities service, see the “The E-utilities In-Depth: Parameters, Syntax and More” (https://www.ncbi.nlm.nih.gov/books/NBK25499/) section of the Entrez Programming Utilities Help.
For further examples of using Entrez resources in the Wolfram Language, see George Mias’s book “Mathematica for Bioinformatics:  A Wolfram Language Approach to Omics”.