Representing 
Biological Sequence Data 
in the Wolfram Language

John Cassel (jcassel), Wolfram

ABSTRACT

This talk will introduce BioSequence, a new structure for representing biological sequence data in the Wolfram Language.  BioSequence can represent DNA, RNA, and amino acid sequences.  BioSequence is integrated with the Wolfram Knowledge Base to allow obtaining the BioSequences associated with gene and protein entities, as well as find the genes, proteins, and SNPs associated with a given biological sequence.  This representation supports the cases where sequence terms are not known through rich handling and pattern matching for degenerate terms.  BioSequence can undertake biological conversion processes such as transcription and translation.  Overall, BioSequence provides the basic representation needed as the Wolfram Language extends its support for computational biology.

Outline

Introduction

Defining Biological Sequence Processing

Purpose of this Talk

Very Brief Biology Primer

Near-Term Functionality Under Development

Representing Biological Sequences with BioSequence

Obtaining Biological Sequences for Entities

Finding Features of BioSequences with FindBioSequenceFeatures

Constructing Molecules from Biological Sequences

Extending Existing Functionality to Support BioSequence

Biological Sequence Complements and Conversions

Matching and Validating Biological Sequences

Resolving Degenerate Sequence Terms

A Broader Functional Picture

Many Functionalities to Pursue

Extensions to Represented Types

Indexing Data and Providing Search and Feature Annotation Capabilities

Alignments/Structural Comparisons

Platform and Deployment Considerations

Conclusion

Introduction

Defining Biological Sequence Processing

The analysis for uncovering relevant comparisons between biological systems

differing by taxonomic category, individual, cell type, developmental stage, cellular process, or fidelity to function 

enabled by representing biologically significant sequentially-connected similar molecules as linear textual sequences

Introduction

Purpose of this Talk

To share proposed Wolfram Language features for biological sequence processing
 
while they are under development 
(undocumented/experimental and subject to change)

to elicit collaboration and feedback

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 that indicate the boundaries of genes 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 dictate that the structure folds into different structures appropriate to different functions, yielding the working protein.

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

Near-Term Functionality Under Development

Representing Biological Sequences with BioSequence

Exact Sequences

BioSeqeuence provides a representation DNA, RNA, and amino acid sequences:

BioSequence["DNA","ATCGGGAGTAGCAAAAGCCCTCTA"]

BioSequence["RNA","AUCGGGAGUAGCAAAAGCCCUCUA"]

BioSequence["AminoAcid","ITGTSTVGVGRGVLGDQKNINTTYSTYYYLQDNTRGNGIFTYDAKYRTTLPGSLWADADNQFFASYDAPAVDAHYYAGVTYDYYK"]

Arbitrary Sequences

BioSeqeuence allows for the representation of arbitrary sequences of a given type:

BioSequence["DNA",63]

BioSequence["AminoAcid",21]

Degenerate Sequences

BioSeqeuence allows for mixing exact and degenerate terms as defined by the relevant IUPAC standards:

BioSequence["DNA","ABCDGHKMNRSTVWXY"]

BioSequence["RNA","ABCDGHKMNRSUVWXY"]

BioSequence["AminoAcid","ABCDEFGHIJKLMNOPQRSTUVWXYZ."]

Unspecified Sequences

BioSequence allows users to provide sequences of unspecified type:

BioSequence["Unspecified","ABC"]

If the user provides no type for the sequence, the sequence is of unspecified type:

BioSequence["ABC"]

Unspecified BioSequences can be converted into a supported type if appropriate:

BioSequenceConvert[BioSequence["ATT"],"DNA"]

Compound Sequences

Combinations of multiple types can be incorporated into a given BioSequence:

BioSequence[{BioSequence["AminoAcid","MGLSDG"],BioSequence["Unspecified","ABC"],BioSequence["DNA","ATG"]}]

Near-Term Functionality Under Development

Obtaining Biological Sequences for Entities

We can obtain BioSequences corresponding to genes and proteins from EntityValue:

EntityValue[(Entity["zinc finger protein 501", Entity["Gene",{"ZNF501",{"Species"->"HomoSapiens"}}], "Entity[\"Gene\", {\"ZNF501\", {\"Species\" -> \"HomoSapiens\"}}]", "gene"]),"ReferenceBioSequence"]

BioSequence will construct the most appropriate sequence for the given entity:

BioSequence[(Entity["BRCA1, DNA repair associated", Entity["Gene",{"BRCA1",{"Species"->"HomoSapiens"}}], "Entity[\"Gene\", {\"BRCA1\", {\"Species\" -> \"HomoSapiens\"}}]", "gene"])]

BioSequence[(Entity["myoglobin", Entity["Protein","MB"], "Entity[\"Protein\", \"MB\"]", "protein"])]

Near-Term Functionality Under Development

Finding Features of BioSequences with FindBioSequenceFeatures

Finding genomic position, genes, and SNPs in biological sequences

Find a position on the reference genome of the given species that reliably corresponds to positions on the given sequence:

FindBioSequenceFeatures[BioSequence[],(Entity["human", Entity["Species","Species:HomoSapiens"], "Entity[\"Species\", \"Species:HomoSapiens\"]", "species specification"])]

Find the positions of a particular gene in the given sequence:

FindBioSequenceFeatures[BioSequence[],(Entity["hemoglobin subunit alpha 2", Entity["Gene",{"HBA2",{"Species"->"HomoSapiens"}}], "Entity[\"Gene\", {\"HBA2\", {\"Species\" -> \"HomoSapiens\"}}]", "gene"])]

Find the position of all human genes in the given sequence:

FindBioSequenceFeatures[BioSequence[],EntityClass["Gene",{"Species"->"HomoSapiens"}]]

Find a SNP in the given sequence:

FindBioSequenceFeatures[BioSequence[],(Entity["rs7480526", Entity["SNP","RS7480526"], "Entity[\"SNP\", \"RS7480526\"]", "SNP"])]

Find all standard SNPs on a given sequence:

FindBioSequenceFeatures[BioSequence[],(EntityClass["SNPs", EntityClass["SNP",All], "EntityClass[\"SNP\", All]", "SNPs"])]

Specifying Desired Results

The default results are the features found and the position of those features on the provided sequence:

FindBioSequenceFeatures[BioSequence[],(Entity["hemoglobin subunit alpha 2", Entity["Gene",{"HBA2",{"Species"->"HomoSapiens"}}], "Entity[\"Gene\", {\"HBA2\", {\"Species\" -> \"HomoSapiens\"}}]", "gene"])->{"Feature","Position"}]

“Sequence” can also be specified to extract the subsequence corresponding to the discovered feature from the provided sequence:

FindBioSequenceFeatures[BioSequence[],(Entity["hemoglobin subunit alpha 2", Entity["Gene",{"HBA2",{"Species"->"HomoSapiens"}}], "Entity[\"Gene\", {\"HBA2\", {\"Species\" -> \"HomoSapiens\"}}]", "gene"])->{"Feature","Position","Sequence"}]

Near-Term Functionality Under Development

Constructing Molecules from Biological Sequences

We can obtain molecules for DNA, RNA, and amino acids, including capping specifications as desired:

SequenceMolecule[BioSequence["DNA","ATGCC"]]

SequenceMolecule[BioSequence["DNA","ATGCC"],"BothCapsDNA"]

SequenceMolecule[BioSequence["RNA","AUCG"]]

SequenceMolecule[BioSequence["AminoAcid","LI"]]

Near-Term Functionality Under Development

Extending Existing Functionality to Support BioSequence

RandomInstance

Generate three random DNA sequences of length 50:

RandomInstance[BioSequence["DNA",50],3]

RandomInstance also allows sampling from more precisely specified sequences:

Dataset[Tally[RandomInstance[BioSequence["DNA","ST"],50]]]

RandomInstance also allows sampling from weighted choices:

RandomInstance[BioSequence["DNA",{0.2,0.8}->{"A","C"},9],6]

Random instances of BioSequence always resolves degenerate terms:

RandomInstance[BioSequence["DNA",{0.2,0.8}->{"A","B"},9],6]

SequenceAlignment

SequenceAlignment works on BioSequences:

SequenceAlignment[BioSequence[],BioSequence[]]

With scoring rules that match degenerate bases as similar, “GG” matches “NN”, and “N” matches “T”, allowing a potentially perfect alignment.

SequenceAlignment[BioSequence[],BioSequence[],SimilarityRules->"SimilarDegenerateBases"]

Similarity Measures

SmithWatermanSimilarity also now works over BioSequences and shows the improvement of considering similar degenerate bases as matching:

SmithWatermanSimilarity[BioSequence[],BioSequence[]]

SmithWatermanSimilarity[BioSequence[],BioSequence[],SimilarityRules->"SimilarDegenerateBases"]

NeedlemanWunschSimilarity also has this upgraded functionality:

NeedlemanWunschSimilarity[BioSequence[],BioSequence[]]

NeedlemanWunschSimilarity[BioSequence["DNA", "AAATTCCAAANNTNCCAAAA"],BioSequence[],SimilarityRules->"SimilarDegenerateBases"]

Longest Common Sequences and Subsequences

Longest common subsequences are formed by finding longest contiguous sequence of bases:

LongestCommonSubsequencePositions[BioSequence[],BioSequence[]]

LongestCommonSubsequence[BioSequence[],BioSequence[]]

Longest common sequences are formed by finding longest sequence of contiguous or disjoint bases:

LongestCommonSequencePositions[BioSequence[],BioSequence[]]

LongestCommonSequence[BioSequence[],BioSequence[]]

Near-Term Functionality Under Development

Biological Sequence Complements and Conversions

Complements

BioSequenceComplement will produce a sequence of complement bases:

BioSequenceComplement[BioSequence["DNA","AGG"]]

BioSequenceReverseComplement will produce a sequence of complement bases in reverse order, appropriately reflecting a complimentary sequence:

BioSequenceReverseComplement[BioSequence["DNA","AGG"]]

Transcription

BioSequenceTranscribe transcribes DNA into RNA:

BioSequenceTranscribe[BioSequence["DNA","GATATC"]]

BioSequenceInverseTranscribe@%

Translation

BioSequenceTranslate performs translation into the corresponding AminoAcid sequence for either DNA or RNA sequences:

translatedSequence=BioSequenceTranslate[BioSequence["DNA","CTGATATAG"]]

BioSequenceInverseTranslate@translatedSequence

Take[BioSequenceInverseTranslate[translatedSequence,All],UpTo[5]]

With BioSequenceInverseTranslationCount, you’ll know how many results will be returned by BioSequenceInverseTranslate with All:

BioSequenceInverseTranslationCount@BioSequence["AminoAcid","LI."]

You can use any genetic translation table entity to specify the desired codon translation:

BioSequenceTranslate[BioSequence["DNA","CTGATATAG"],BioSequenceTranslationTable->(Entity["NCBI::VertebrateMitochondrial::4p4", Entity["GeneticTranslationTable","NCBI::VertebrateMitochondrial::4p4"], "Entity[\"GeneticTranslationTable\", \"NCBI::VertebrateMitochondrial::4p4\"]", "GeneticTranslationTable"])]

Offsets and Frame Adjustments in Translation

If the sequence is not framed for a correct codon translation, translation may generate a sequence with drastically different biological meaning:

BioSequenceTranslate[BioSequence["DNA","ACTGATATAG"]]

One easy way to handle this is merely to skip up to a start codon, using the sequence itself to shift into the appropriate frame.

BioSequenceTranslate[BioSequence["DNA","ACTGATATAC"],Automatic,"ToStartCodon"]

A negative compensating offset can be provided to move the sequence backward on the frame, dropping the initial terms:

BioSequenceTranslate[BioSequence["DNA","ACTGATATAG"],Automatic,-1]

Alternatively, a positive offset moves the sequence forward along the frame, introducing degenerate terms to the start of the sequence:

BioSequenceTranslate[BioSequence["DNA","ACTGATATAG"],Automatic,2]

Conversion in General

BioSequenceConvert provides a simple interface for conversions between types:

BioSequenceConvert[BioSequence["DNA","TCC"],"AminoAcid"]

BioSequenceConvert[BioSequence["AminoAcid","S"],"DNA",All]

Use BioSequenceConversionCount before calling BioSequenceConvert with All to be sure the size of the output will be manageable:

BioSequenceConversionCount[BioSequence["AminoAcid","HNDKRALL"],"RNA"]

Near-Term Functionality Under Development

Matching and Validating Biological Sequences

Validating Sequences with BioSequenceQ

BioSequenceQ determines if a given argument is a BioSequence of the given specification:

BioSequenceQ[BioSequence[]]

BioSequenceQ[BioSequence[], "FullySpecified"]

BioSequenceQ[BioSequence[], "DNA"]

BioSequenceQ[BioSequence[], "FullySpecifiedDNA"]

Matching Sequences with BioSequenceMatchQ

BioSequenceMatchQ allows patterns to match degenerate bases:

BioSequenceMatchQ[BioSequence["DNA","AGGGCSCGSGCGSGCSGGCSGC"],"A"~~("S"..)~~"Y"]

Near-Term Functionality Under Development

Resolving Degenerate Sequence Terms

Single and Complete Resolution

BioSequenceResolveDegenerateTerms resolves unspecified bases into an arbitrary, repeatable instance with compatible bases.  For DNA:

BioSequenceResolveDegenerateTerms[BioSequence["DNA","WST"]]

All potential resolutions can be obtained:

BioSequenceResolveDegenerateTerms[BioSequence["DNA","WST"],All]

Avoiding Combinatorial Explosion

BioSequenceDegenerateTermResolutionCount counts the number of expanded results, to help avoid calling All on BioSequenceResolveDegenerateTerms if the number of results would be too large:

BioSequenceDegenerateTermResolutionCount[BioSequence["DNA","WSABNBNAVNNVBDHBV"]]

A Broader Functional Picture

Many Functionalities to Pursue

Extensions to Represented Types

Additional Data Sources and Formats

Data Management and Preprocessing Functionality

Indexing Data and Providing Search and Feature Annotation Capabilities

Data Validation/Anomaly Detection

Alignments/Structural Comparisons

Clustering

Dimensionality Reduction

Model Building and Providing Classification Abilities

Pattern Discovery/Motif Finding

Annotation/Enrichment

Visualization

Explicit Models for “Spaces” of Biological Sequences

Platform and Deployment Considerations

Modeling Dynamics

Optimization

A Broader Functional Picture

Extensions to Represented Types

The behavior of BioSequence could be extended by users defining additional BioSequenceTypes:

BioSequenceQ[BioSequence["tRNA","GCIAUΨ"]]

With a means to also specify conversions, users might enjoy integrations with existing functionality, such as BioSequenceConvert making organism-specific wobble complements between our RNA and ribosome-specific tRNA definitions.

Type extension will help support new research in synthetic biology, such as  Hachimoji DNA.

A Broader Functional Picture

Indexing Data and Providing Search and Feature Annotation Capabilities

Potential types of indexing and search to be supported include:

• sequences, prebuilt and included: currently GenomeLookup, SNP finding in development as shown, potentially including ProteinLookup and others

• by feature or category, prebuilt and included:  such as finding SNPs by popular category (features, diseases, heritage, evolutionary markers), or genes by associated structure or functionality

• index creation tools and search over user-created indexes: for user-provided data or data from their preferred sources

• external resources through utilizing APIs: for users to expose existing search services to the Wolfram environment

A Broader Functional Picture

Alignments/Structural Comparisons

Key functionalities in structural comparison include:

• Many-to-many sequence-to-feature alignment:  for example, short-read alignments (such as RNA-seq).  Consider a “TabulateBioSequenceFeatures” function over a corpus of sequences.

• Sequence assembly: both reference-based and denovo assembly

• Multiple sequence alignment:  extending from pairwise to multiple alignment, with continued gap penalty support

• Many-to-many alignment comparisons: for example, differential expression

• Enrichment analysis:  enrichment analysis concerns annotating expression levels with categories, whether associated function or pathways, and using statistical tests to find if particular categories are significantly more or less expressed compared to the baseline

A Broader Functional Picture

Platform and Deployment Considerations

Biological sequence processing could be undertaken in all forms of deployment: desktop, server, virtualized, distributed, parallelized, containerized, and embedded systems.

Containerized deployment is particularly popular, as components of functionality can be added with modularity to existing workflows.

WL notebooks, whether on the desktop or cloud, are ideal as:

•  a prototyping platform from which the developed approach might then be containerized.

• an environment for interactive reports generated as the result of a containerized workflow.

Conclusion

We are starting to fill in how the Wolfram Language can support biological sequence processing.

There are many directions in which this support will develop and your suggestions, both broad and detailed, can yet be incorporated into the design.

Conclusion

Thank you!

Please feel free to email jcassel@wolfram.com to express your interests in biological sequence processing and other biology/bioinformatics/biotechnology topics.

Thanks to the biologists and other scientists I work with for their input and support.  Any remaining errors and omissions are my own.