Genomic Features, Biological Pathways and Phylogenetic Trees in Wolfram Language

Keiko Hirayama, Wolfram|Alpha

#WolframTechConf

Genomic analysis functions overview

Wolfram Function Repository

A public resource that hosts an expanding collection of contributed standalone functions

Biology functions suitable for immediate use in any Wolfram Language computation

Genome feature reference

Retrieve features that overlap a given genomic region

Key genomic features, such as nucleotide sequences, genes, and variants, are effectively studied using public resource function

The EnsemblGenomeRegion function, based on the Ensembl database allows to retrieve features that overlap a specified genomic region

Query regions are specified by the chromosome name and integer coordinates

ResourceFunction["EnsemblGenomeRegion"]

Find a nucleotide sequence for a specified region of human chromosome 6:

["6",{33570001,33600000},"BioSequence"]

Use the Species option to specify the organism of the genomics feature:

["6",{33570001,33600000},"BioSequence","Species"->"Bos taurus"]

Find motifs in the specified genome region:

motif=["6",{33570001,33600000},"Feature"->"motif"]

Explore motif sequences and associated transcription factors:

motifseq=motif[GroupBy["TranscriptionFactorComplex"],All,ResourceFunction["EnsemblGenomeRegion"][#SeqRegionName,{#Start,#End},"BioSequence","Strand"->#Strand]&][]

[#[[2]],PlotLabel->#[[1]]]&/@Normal@Normal@motifseq//GraphicsColumn

Genome feature reference

Find protein coding genes that overlap a specified region of human chromosome 17:

genereg=["17",{43047001,43048000},"Feature"->"gene","BioType"->"protein_coding"]

BioDBnetGeneData accesses information on genes available in major biological databases.

Retrieve a dataset of all available properties for a gene entity:

brca1=[genereg[1,"ExternalName"]]

To visualize the protein structure, retrieve the BioMolecule associated with the specified PDB identifier:

bioMolecule=BioMolecule@brca1["PDBID",1]

Visualize the protein structure:

BioMoleculePlot3D[bioMolecule]

Genome feature reference

In addition to genes and proteins, we can also identify variant features that overlap the same genomic region:

variant=["17",{43047001,43048000},"Feature"->"variation"]

From the set of all variants, select those classified as pathogenic:

pathogenic=variant[Select[MemberQ[#"ClinicalSignificance","pathogenic"]&]]

The NCBIGenomicSNPData function provides detailed information on the identified single nucleotide polymorphism (SNP).

Find the variant with clinical significance information:

[pathogenic[[3,1]],"ClinicalSignificance"]

Use the EnsemblHomology function to get homology information for a selected human gene. Homologous genes share evolutionary conserved genomic sequences that often serve similar functions:

homolog=[(Entity["human", Entity["TaxonomicSpecies","HomoSapiens::4pydj"], "Entity[\"TaxonomicSpecies\", \"HomoSapiens::4pydj\"]", "species specification"]),genereg[1,"ExternalName"]]

Compare the homologous peptide sequences of lowland gorilla and Sumatran orangutan:

Diff[Sequence@@Normal@homolog[Select[!FreeQ[#,"gorilla gorilla"|"pongo abelii"]&]][All,"AlignedSequence"],"CharacterAlignment",]

Gene tree analysis

Homological relations can be visualized using a tree structure

The gene tree illustrates the evolutionary relationships among a group of homologous genes that share a common evolutionary origin

Retrieve the dataset of the gene tree associated with CLOCK genes, which are involved in circadian rhythms:

["CLOCK","Dataset"]

Retrieve the Tree object for the CLOCK genes to analyze the relations of associated species.

["CLOCK","Tree"]

Gene tree analysis

Tree grafter

In addition to a known gene, the provided protein sequence can also be classified

PANTHERTreeGrafter is based on the PANTHER tree grafter tool, which classifies the identified protein sequence by determining its most likely position in the protein family tree.

Retrieve the annotation result for a protein sequence:

[BioSequence[],"Dataset"]

Highlight the location of the grafted sequence on the protein family tree. Each node is labeled with the corresponding accession from the dataset:

[BioSequence[],"TreeGraphic"]

Genomic sequence analysis

Another common way to analyze biological sequence is to use the Basic Local Alignment Search Tool (BLAST)

BLAST is a tool provided by the NCBI for aligning query sequences against those present in a selected target database

Perform a sequence similarity search against the human genome using a provided nucleotide sequence:

querySequence=BioSequence[];
alignedSequenceSet=[querySequence["SequenceString"],{"Species"->(Entity["human", Entity["TaxonomicSpecies","HomoSapiens::4pydj"], "Entity[\"TaxonomicSpecies\", \"HomoSapiens::4pydj\"]", "species specification"]),MaxItems->3}]

The result includes the accession to the identified nucleotide sequences and position information as well as aligned sequences.

Compare the original with the aligned sequence and identify the variant:

Diff[alignedSequenceSet[1,"Sequence"],querySequence["SequenceString"],]

Genomic sequence analysis

In order to find the known SNPs (single nucleotide polymorphisms) associated with the identified variant, we need to compute the variant position on the coding sequence.

The coding sequence can be retrieved by accessing the NCBI nucleotide database, using the NCBIEntrezData function:

codingseq=[alignedSequenceSet[1,"RefSeqAccession"],"EFetch",{"Database"->"nucleotide","RetType"->"fasta_cds_na","RetMode"->"BioSequence"}]

Diff[codingseq["SequenceString"],querySequence["SequenceString"],]

Here is the position of the variant within the full nucleotide sequence, including the non-coding region:

variantpos=Diff[alignedSequenceSet[1,"Sequence"],querySequence["SequenceString"]]["Edits"]

Compute the length of the non-coding sequence at the beginning:

alignment=SequenceAlignment[codingseq["SequenceString"],querySequence["SequenceString"]];
alignment[[1]]
noncodinglen=alignment[[1,2]]//StringLength

Find the variant position on the coding sequence:

allelepos=variantpos[[1,1,1]]-noncodinglen

Using this position information, find the change in the amino acid using the BioSequenceTranslate:

Diff[BioSequenceTranslate[BioSequence[codingseq["SequenceString"]]],BioSequenceTranslate[BioSequence[StringReplacePart[codingseq["SequenceString"],variantpos[[1,-1,-1]],{allelepos,allelepos}]]],]

E (glutamate) is replaced by K (lysine).

Genomic sequence analysis

Using the identified nucleotide sequence accession and variant position formatted in HGVS nomenclature, search for potential known clinical variant in the NCBI ClinVar database:

accession=alignedSequenceSet[1,"RefSeqAccession"]["ExternalID"]
hgvs="c."<>ToString[allelepos]<>variantpos[[1,1,2]]<>">"<>variantpos[[1,2,2]]

clinvar=[{accession,hgvs},"ESearch",{"Database"->"clinvar"}]

Get the detail on the identified clinical variant:

clinvarset=[Normal@clinvar["IDList"],"ESummary",{"Database"->"clinvar"}]

Further information on the identified variant, such as the population frequency and clinical significance, is retrieved using the NCBIGenomicSNPData function:

varID=clinvarset[1,"variation_set"][1,2,-1,-1];

[varID,"FrequencyData"]

Retrieve clinical significance information:

[varID,"ClinicalSignificance"]

The variant may influence alcohol sensitivity and susceptibility to alcohol-related conditions.

Phenotype analysis

Phenotypes are traits that are associated with the genomic variants

The EnsemblPhenotype function helps to identify genomic features associated with specific traits

Find single nucleotide polymorphisms (SNPs) associated with coffee consumption:

snps=["coffee consumption","GenomicFeatures"]

Analyze SNPs affecting the coffee consumption rate:

coffeConsumptionRateAffectingSNP=snps[Select[!MissingQ[#BetaCoefficient]&&!StringFreeQ[#BetaCoefficient,"unit"|"cups"]&]][All,][GroupBy[First->Last],MeanAround][ReverseSortBy[#/.a_Around:>a[[1]]&]]

Plot SNPs associated with the change in coffee consumption rate:

ListPlot[coffeConsumptionRateAffectingSNP,]

Use the NCBIGenomicSNPData function to find population frequency information for one of the notable SNPs.

Compare allele frequencies across 12 diverse populations:

populationfrequency=[Keys[coffeConsumptionRateAffectingSNP][1],"ALFAFrequencyData"]

BarChart[populationfrequency[All,#"AltAlleleFrequency"&],]

Phenotype analysis

Find genes associated with the coffee consumption rate:

coffeConsumptionRateAffectingGene=Dataset[GroupBy[Flatten@Normal@snps[Select[!MissingQ[#AssociatedGene]&&!MissingQ[#BetaCoefficient]&&!StringFreeQ[#BetaCoefficient,"unit"|"cups"]&]][All,][ReverseSortBy[Mean[#]&]],Last->First,MeanAround]]

Plot genes associated with the change in coffee consumption rate:

ListPlot[coffeConsumptionRateAffectingGene,]

These studies are effective in highlighting SNPs and genes associated with known traits.

Enrichment analysis and biological pathways

The DAVID (Database for Annotation, Visualization, and Integrated Discovery) provides annotation tools to help interpret the biological functions of genes

The DAVIDGeneEnrichmentAnalysis identifies enriched (overrepresented) biological themes, such as Gene Ontology (GO) terms and biological pathways, for a given set of genes (e.g. those found to be differentially expressed) and classifies them into enriched functional groups

For a list of selected NCBI Entrez genes, identify enriched Gene Ontology terms:

goGenes=[{(ExternalIdentifierURL[ExternalIdentifier["EntrezGeneID","581"], "Entrez Gene", "581", https://www.ncbi.nlm.nih.gov/gene/581]),(ExternalIdentifierURL[ExternalIdentifier["EntrezGeneID","637"], "Entrez Gene", "637", https://www.ncbi.nlm.nih.gov/gene/637]),(ExternalIdentifierURL[ExternalIdentifier["EntrezGeneID","2810"], "Entrez Gene", "2810", https://www.ncbi.nlm.nih.gov/gene/2810]),(ExternalIdentifierURL[ExternalIdentifier["EntrezGeneID","7157"], "Entrez Gene", "7157", https://www.ncbi.nlm.nih.gov/gene/7157])},"Annotation"->"GOTERM_BP_3"]

Compare count of genes and p-values across associated GO terms:

With[{goPValue=},
Legended[BarChart[MapThread[Style[#1,ColorData["LightTemperatureMap"][Rescale[-Log10[#2],MinMax[-Log10/@goPValue[[All,3]]]]]]&,{goPValue[[All,2]],goPValue[[All,3]]}],],]]

The orange to blue gradient reflects p-values, where orange corresponds to more statistically significant results.

Enrichment analysis and biological pathways

For the same set of genes, identify enriched KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways instead:

enrichedKEGGPath=[{(ExternalIdentifierURL[ExternalIdentifier["EntrezGeneID","581"], "Entrez Gene", "581", https://www.ncbi.nlm.nih.gov/gene/581]),(ExternalIdentifierURL[ExternalIdentifier["EntrezGeneID","637"], "Entrez Gene", "637", https://www.ncbi.nlm.nih.gov/gene/637]),(ExternalIdentifierURL[ExternalIdentifier["EntrezGeneID","2810"], "Entrez Gene", "2810", https://www.ncbi.nlm.nih.gov/gene/2810]),(ExternalIdentifierURL[ExternalIdentifier["EntrezGeneID","7157"], "Entrez Gene", "7157", https://www.ncbi.nlm.nih.gov/gene/7157])},"Annotation"->{"KEGG_PATHWAY"}]

We can visualize how these provided genes comprises the biological pathways.

First, find the gene symbols for the provided Entrez genes:

geneSymbols=[{(ExternalIdentifierURL[ExternalIdentifier["EntrezGeneID","581"], "Entrez Gene", "581", https://www.ncbi.nlm.nih.gov/gene/581]),(ExternalIdentifierURL[ExternalIdentifier["EntrezGeneID","637"], "Entrez Gene", "637", https://www.ncbi.nlm.nih.gov/gene/637]),(ExternalIdentifierURL[ExternalIdentifier["EntrezGeneID","2810"], "Entrez Gene", "2810", https://www.ncbi.nlm.nih.gov/gene/2810]),(ExternalIdentifierURL[ExternalIdentifier["EntrezGeneID","7157"], "Entrez Gene", "7157", https://www.ncbi.nlm.nih.gov/gene/7157])},"ESummary","Database"->"gene"][All,"name"]//Normal

Then visualize the relevant genes on the KEGGPathway graph:

keggPathGraph=["hsa","04115","Graph",VertexStyle->Prepend[Map[#->Yellow&,geneSymbols],_->LightBlue],]

The graph highlights the genes involved in the p53 signaling pathway.

Enrichment analysis and biological pathways

Enriched pathways are analyzed based on expression levels of genes

Expression analysis measures the overall activity of identified biological pathways based on the expression levels of the genes that comprise those pathways

The ReactomeEnrichmentAnalysis queries the Reactome database for expression analysis

Query Reactome for pathway analysis using a sample dataset of Ensembl genes and their expression levels:

geneexpressiondata=

Get a summary result including a token of the analysis and statistics:

expressionsummary=[geneexpressiondata,"Summary"]

Retrieve a dataset of enriched Reactome pathways:

expressionpath=[geneexpressiondata]

Use a token from the analysis and plug it into the ReactomePathways function to retrieve the result overlaid on a selected pathway:

[expressionpath[1,1]["ExternalID"],"Image","Token"->expressionsummary["token"],]

The expressed gene are colored on the pathway graphic. The biological pathway analysis gives insight on how differentially expressed genomic features may affect the biological activities.

Summary

Wolfram Language offers programmatic access to a wide range of biological datasets, including genes, homology, and variants

The functions available in the Wolfram Function Repository enable effective analysis of gene expression, phylogenetic relationships, biological pathways, and other fundamental genomic features