Files
semantica/cookbook/use_cases/biomedical/02_Genomic_Variant_Analysis.ipynb
T

27 KiB

Open In Colab

Genomic Variant Analysis - Graph Analytics & Pathway Analysis

Overview

This notebook demonstrates genomic variant analysis using Semantica's modular architecture with focus on graph analytics, pathway analysis, and temporal knowledge graphs. The pipeline analyzes genomic data to extract variant entities, build temporal genomic knowledge graphs, and analyze disease associations through reasoning.

Key Features

  • Graph Analytics Focus: Emphasizes graph reasoning, centrality measures, and pathway analysis
  • Temporal Analysis: Builds temporal genomic knowledge graphs to track variant evolution
  • Disease Association: Analyzes relationships between variants, genes, and diseases
  • Pathway Analysis: Uses graph traversal to identify biological pathways
  • Impact Prediction: Predicts variant impact using graph-based reasoning

What You'll Learn

  • How to use Semantica modules directly for genomic analysis
  • How to ingest genomic data from multiple sources
  • How to extract variant, gene, and disease entities
  • How to build temporal knowledge graphs
  • How to perform graph analytics (centrality, communities)
  • How to use temporal queries for variant evolution
  • How to analyze pathways using reasoning
  • How to visualize and export genomic knowledge graphs

Pipeline Flow

graph LR
    A[Data Ingestion] --> B[Text Processing]
    B --> C[Entity Extraction]
    C --> D[Relationship Extraction]
    D --> E[Deduplication]
    E --> F[Temporal KG]
    F --> G[Graph Analytics]
    F --> H[Temporal Queries]
    G --> I[Pathway Analysis]
    H --> I
    I --> J[Disease Associations]
    J --> K[Visualization]

Installation

Install Semantica and required dependencies:

In [ ]:
%pip install -qU semantica networkx matplotlib plotly pandas groq sentence-transformers

Configuration & Setup

Set up environment variables and configuration constants.

In [ ]:
import os

os.environ["GROQ_API_KEY"] = os.getenv("GROQ_API_KEY", "")
In [ ]:
CHUNK_SIZE = 1000
CHUNK_OVERLAP = 200
TEMPORAL_GRANULARITY = "day"

Ingesting Genomic Data from Multiple Sources

Ingest data from comprehensive genomic sources including PubMed RSS feeds, preprint servers, and journal feeds.

In [ ]:
from semantica.ingest import FeedIngestor, FileIngestor
import os
from contextlib import redirect_stderr
from io import StringIO

os.makedirs("data", exist_ok=True)

feed_sources = [
    # PubMed RSS Feeds (simplified, working format)
    ("PubMed - Genetics", "https://pubmed.ncbi.nlm.nih.gov/rss/search/1?term=genetics&limit=10"),
    ("PubMed - Genomics", "https://pubmed.ncbi.nlm.nih.gov/rss/search/1?term=genomics&limit=10"),
    ("PubMed - Variant Analysis", "https://pubmed.ncbi.nlm.nih.gov/rss/search/1?term=variant+analysis&limit=10"),
    ("PubMed - GWAS", "https://pubmed.ncbi.nlm.nih.gov/rss/search/1?term=GWAS&limit=10"),
    ("PubMed - Genomic Medicine", "https://pubmed.ncbi.nlm.nih.gov/rss/search/1?term=genomic+medicine&limit=10"),
    ("PubMed - Precision Medicine", "https://pubmed.ncbi.nlm.nih.gov/rss/search/1?term=precision+medicine&limit=10"),
    ("PubMed - Pharmacogenomics", "https://pubmed.ncbi.nlm.nih.gov/rss/search/1?term=pharmacogenomics&limit=10"),
    
    # Nature Feeds (working format)
    ("Nature Genetics", "https://www.nature.com/subjects/genetics.rss"),
    ("Nature - Genomics", "https://www.nature.com/subjects/genomics.rss"),
    
    # PLOS Journals (working Atom feeds)
    ("PLOS Genetics", "https://journals.plos.org/plosgenetics/feed/atom"),
    ("PLOS ONE - Genetics", "https://journals.plos.org/plosone/feed/atom"),
    
    # Other working feeds
    ("Genome Research", "https://genome.cshlp.org/rss/current.xml"),
]

feed_ingestor = FeedIngestor()
all_documents = []

print(f"Ingesting from {len(feed_sources)} feed sources...")
for i, (feed_name, feed_url) in enumerate(feed_sources, 1):
    try:
        with redirect_stderr(StringIO()):
            feed_data = feed_ingestor.ingest_feed(feed_url, validate=False)
        
        feed_count = 0
        for item in feed_data.items:
            if not item.content:
                item.content = item.description or item.title or ""
            if item.content:
                if not hasattr(item, 'metadata'):
                    item.metadata = {}
                item.metadata['source'] = feed_name
                all_documents.append(item)
                feed_count += 1
        
        if feed_count > 0:
            print(f"  [{i}/{len(feed_sources)}] {feed_name}: {feed_count} documents")
    except Exception as e:
        print(f"  [{i}/{len(feed_sources)}] {feed_name}: Failed")
        continue

# Always include fallback variant data for demonstration
variant_data = """
Variant rs699 is located in the AGT gene and associated with hypertension.
Variant rs7412 in APOE gene is linked to Alzheimer's disease risk.
BRCA1 variant c.5266dupC increases breast cancer susceptibility.
CFTR variant F508del causes cystic fibrosis.
Variant rs1800566 in NAT2 gene affects drug metabolism.
Variant rs1042713 in ADRB2 gene is associated with asthma response.
TP53 variant R273H is linked to multiple cancer types.
Variant rs1799853 in CYP2C9 gene affects warfarin metabolism.
Variant rs1057910 in CYP2C9 affects phenytoin metabolism.
Variant rs9923231 in VKORC1 gene influences warfarin dosing.
"""

with open("data/variants.txt", "w") as f:
    f.write(variant_data)

file_ingestor = FileIngestor()
fallback_docs = file_ingestor.ingest("data/variants.txt")
all_documents.extend(fallback_docs)

documents = all_documents
print(f"\nTotal ingested: {len(documents)} documents")

Normalizing and Chunking Genomic Documents

Clean and normalize text, then split into chunks using entity-aware chunking to preserve variant/gene entity boundaries.

In [ ]:
from semantica.normalize import TextNormalizer
from semantica.split import TextSplitter

normalizer = TextNormalizer()
splitter = TextSplitter(
    method="entity_aware",
    ner_method="spacy",
    chunk_size=CHUNK_SIZE,
    chunk_overlap=CHUNK_OVERLAP
)

print(f"Normalizing {len(documents)} documents...")
normalized_documents = []
for i, doc in enumerate(documents, 1):
    normalized_text = normalizer.normalize(
        doc.content if hasattr(doc, 'content') else str(doc),
        clean_html=True,
        normalize_entities=True,
        remove_extra_whitespace=True,
        lowercase=False
    )
    normalized_documents.append(normalized_text)
    if i % 50 == 0 or i == len(documents):
        print(f"  Normalized {i}/{len(documents)} documents...")

print(f"Chunking {len(normalized_documents)} documents...")
chunked_documents = []
for i, doc_text in enumerate(normalized_documents, 1):
    try:
        with redirect_stderr(StringIO()):
            chunks = splitter.split(doc_text)
        chunked_documents.extend(chunks)
    except Exception:
        simple_splitter = TextSplitter(method="recursive", chunk_size=CHUNK_SIZE, chunk_overlap=CHUNK_OVERLAP)
        chunks = simple_splitter.split(doc_text)
        chunked_documents.extend(chunks)
    if i % 50 == 0 or i == len(normalized_documents):
        print(f"  Chunked {i}/{len(normalized_documents)} documents ({len(chunked_documents)} chunks so far)")

print(f"Created {len(chunked_documents)} chunks from {len(normalized_documents)} documents")
In [ ]:
from semantica.semantic_extract import NERExtractor

# Using spaCy ML method (similar to Drug Discovery Pipeline)
entity_extractor = NERExtractor(method="ml", model="en_core_web_sm")

all_entities = []
print(f"Extracting entities from {len(chunked_documents)} chunks...")

for i, chunk in enumerate(chunked_documents, 1):
    chunk_text = chunk.text if hasattr(chunk, 'text') else str(chunk)
    try:
        entities = entity_extractor.extract_entities(chunk_text)
        all_entities.extend(entities)
    except Exception:
        continue
    
    if i % 20 == 0 or i == len(chunked_documents):
        remaining = len(chunked_documents) - i
        print(f"  Processed {i}/{len(chunked_documents)} chunks ({len(all_entities)} entities found, {remaining} remaining)")

# Filter entities - spaCy returns standard types, map to genomic categories
# Look for variant patterns (rs numbers, c. notation, etc.)
variants = [
    e for e in all_entities 
    if (e.text.startswith("rs") or 
        "c." in e.text.lower() or 
        "variant" in e.text.lower() or
        e.label == "PRODUCT" and any(kw in e.text.lower() for kw in ["rs", "variant", "mutation"]))
]

# Look for gene patterns (gene names, protein names)
genes = [
    e for e in all_entities 
    if (e.label == "ORG" or 
        e.label == "PRODUCT" or
        any(kw in e.text.lower() for kw in ["gene", "protein", "enzyme", "receptor", "kinase"]))
]

# Look for disease patterns
diseases = [
    e for e in all_entities 
    if (e.label == "ORG" or
        any(kw in e.text.lower() for kw in ["disease", "syndrome", "disorder", "cancer", "hypertension", "alzheimer"]))
]

print(f"Extracted {len(variants)} variants, {len(genes)} genes, {len(diseases)} diseases")

Extracting Genomic Relationships

Extract relationships between variants, genes, and diseases to understand genomic associations.

In [ ]:
from semantica.semantic_extract import RelationExtractor

# Using spaCy dependency parsing (similar to Drug Discovery Pipeline)
relation_extractor = RelationExtractor(method="dependency", model="en_core_web_sm")

all_relationships = []
print(f"Extracting relationships from {len(chunked_documents)} chunks...")

for i, chunk in enumerate(chunked_documents, 1):
    chunk_text = chunk.text if hasattr(chunk, 'text') else str(chunk)
    try:
        relationships = relation_extractor.extract_relations(
            chunk_text,
            entities=all_entities,
            relation_types=["associated_with", "located_in", "causes", "increases_risk", "affects", "linked_to"]
        )
        all_relationships.extend(relationships)
    except Exception:
        continue
    
    if i % 20 == 0 or i == len(chunked_documents):
        print(f"  Processed {i}/{len(chunked_documents)} chunks ({len(all_relationships)} relationships found)")

print(f"Extracted {len(all_relationships)} relationships")

Building Temporal Genomic Knowledge Graph

Construct a temporal knowledge graph from extracted entities and relationships to enable time-aware analysis and variant evolution tracking.

Conflict Detection and Resolution

Detect and resolve conflicts in genomic variant data from multiple research sources.

  • Detection Method: Entity and relationship conflict detection identifies discrepancies in variant-gene-disease associations across sources
  • Resolution Strategy: Credibility-weighted resolution prioritizes higher-credibility sources (e.g., Nature Genetics over preprints)
  • Use Case: Handles conflicting information when multiple sources report different variant associations, disease risks, or gene locations
In [ ]:
from semantica.conflicts import ConflictDetector, ConflictResolver

# Initialize with best strategies for genomic analysis
detector = ConflictDetector()
resolver = ConflictResolver(default_strategy="credibility_weighted")

# Convert entities to format expected by detector
entities = [
    {
        "id": ent.text if hasattr(ent, 'text') else str(ent),
        "name": ent.text if hasattr(ent, 'text') else str(ent),
        "type": ent.label if hasattr(ent, 'label') else "ENTITY",
        "confidence": getattr(ent, 'confidence', 1.0),
        "source": ent.metadata.get("source", "unknown") if hasattr(ent, 'metadata') and ent.metadata else "unknown"
    }
    for ent in all_entities if hasattr(ent, 'text') or hasattr(ent, 'label')
]

# Convert relationships to format expected by detector
relationships = [
    {
        "id": f"{rel.subject.text}_{rel.object.text}_{rel.predicate}" if hasattr(rel, 'subject') else f"{i}",
        "source_id": rel.subject.text if hasattr(rel, 'subject') else str(rel.get("source", "")),
        "target_id": rel.object.text if hasattr(rel, 'object') else str(rel.get("target", "")),
        "type": rel.predicate if hasattr(rel, 'predicate') else rel.get("type", "related_to"),
        "confidence": getattr(rel, 'confidence', 1.0),
        "properties": rel.metadata if hasattr(rel, 'metadata') else {},
        "source": rel.metadata.get("source", "unknown") if hasattr(rel, 'metadata') and rel.metadata else "unknown"
    }
    for i, rel in enumerate(all_relationships) if hasattr(rel, 'subject') or isinstance(rel, dict)
]

# Detect both entity and relationship conflicts
print(f"Detecting conflicts in {len(entities)} entities, {len(relationships)} relationships...")

# Detect entity conflicts
entity_conflicts = detector.detect_conflicts(entities)
print(f"Detected {len(entity_conflicts)} entity conflicts")

# Detect relationship conflicts
relationship_conflicts = detector.detect_relationship_conflicts(relationships)
print(f"Detected {len(relationship_conflicts)} relationship conflicts")

# Resolve entity conflicts
if entity_conflicts:
    resolver.resolve_conflicts(entity_conflicts, strategy="credibility_weighted")
    print(f"Resolved {len(entity_conflicts)} entity conflicts")

# Resolve relationship conflicts
if relationship_conflicts:
    resolver.resolve_conflicts(relationship_conflicts, strategy="credibility_weighted")
    print(f"Resolved {len(relationship_conflicts)} relationship conflicts")

# GraphBuilder will use resolve_conflicts=True to apply resolutions automatically
if entity_conflicts or relationship_conflicts:
    print("Conflicts resolved. GraphBuilder will use cleaned data.")
else:
    print("No conflicts detected. Data is clean.")
In [ ]:
from semantica.kg import GraphBuilder

# Conflicts already detected and resolved in previous cell
# Enable temporal features for genomic variant tracking
graph_builder = GraphBuilder(
    resolve_conflicts=False,  # Conflicts already handled
    enable_temporal=True,
    temporal_granularity=TEMPORAL_GRANULARITY
)

print(f"Building temporal knowledge graph from {len(all_entities)} entities, {len(all_relationships)} relationships...")
kg = graph_builder.build({
    "entities": all_entities,
    "relationships": all_relationships
})

entities_count = len(kg.get('entities', []))
relationships_count = len(kg.get('relationships', []))
print(f"Graph: {entities_count} entities, {relationships_count} relationships")

Analyzing Graph Structure

Perform comprehensive graph analytics including centrality measures, community detection, and connectivity analysis.

In [ ]:
from semantica.kg import GraphAnalyzer, CentralityCalculator, CommunityDetector

graph_analyzer = GraphAnalyzer()
centrality_calc = CentralityCalculator()
community_detector = CommunityDetector()

analysis = graph_analyzer.analyze_graph(kg)

degree_centrality = centrality_calc.calculate_degree_centrality(kg)
betweenness_centrality = centrality_calc.calculate_betweenness_centrality(kg)
closeness_centrality = centrality_calc.calculate_closeness_centrality(kg)

communities = community_detector.detect_communities(kg, method="louvain")
connectivity = graph_analyzer.analyze_connectivity(kg)

print(f"Graph analytics:")
print(f"  - Communities: {len(communities)}")
print(f"  - Connected components: {len(connectivity.get('components', []))}")
print(f"  - Graph density: {analysis.get('density', 0):.3f}")

Temporal Graph Queries

Query the temporal knowledge graph at specific time points, analyze temporal evolution, and detect temporal patterns.

In [ ]:
from semantica.kg import TemporalGraphQuery

temporal_query = TemporalGraphQuery(
    enable_temporal_reasoning=True,
    temporal_granularity=TEMPORAL_GRANULARITY
)

# Query variants at specific time point
query_results = temporal_query.query_at_time(
    kg,
    query="Variant",
    at_time="2024-01-01"
)

# Analyze graph evolution
evolution = temporal_query.analyze_evolution(kg)

# Detect temporal patterns
pattern_results = temporal_query.query_temporal_pattern(
    kg,
    pattern="sequence"
)

print(f"Temporal query: {query_results.get('num_relationships', 0)} relationships valid at query time")
print(f"Evolution analysis: {evolution.get('num_relationships', 0)} relationships tracked")
print(f"Temporal patterns detected: {pattern_results.get('num_patterns', 0)}")

Pathway Analysis & Reasoning

Use graph reasoning to find pathways between variants and diseases, and infer biological pathways through logical reasoning.

In [ ]:
from semantica.reasoning import Reasoner
from semantica.kg import GraphAnalyzer

reasoner = Reasoner()
graph_analyzer = GraphAnalyzer()

# Find entities by type
variants = [e for e in kg.get('entities', []) if e.get('type') == 'Variant']
diseases = [e for e in kg.get('entities', []) if e.get('type') == 'Disease']

print(f"Found {len(variants)} variants and {len(diseases)} diseases")

# Find pathways
pathways = []
for variant in variants[:5]:
    variant_id = variant.get('id') or variant.get('name')
    for disease in diseases[:3]:
        disease_id = disease.get('id') or disease.get('name')
        path = graph_analyzer.connectivity_analyzer.calculate_shortest_paths(
            kg, source=variant_id, target=disease_id
        )
        if path.get('exists'):
            pathways.append({
                'variant': variant_id,
                'disease': disease_id,
                'distance': path.get('distance', -1)
            })

# Add rule and infer facts
reasoner.add_rule("IF Variant associated_with Gene AND Gene causes Disease THEN Variant increases_risk Disease")
inferred_facts = reasoner.infer_facts(kg)

print(f"Pathway analysis: {len(pathways)} variant-disease pathways found")
print(f"Inferred facts: {len(inferred_facts)}")

Analyzing Disease Associations

Use graph traversal to find variant-disease associations and calculate association scores.

In [ ]:
from semantica.kg import GraphAnalyzer

graph_analyzer = GraphAnalyzer()

# Find entities by type
variants = [e for e in kg.get('entities', []) if e.get('type') == 'Variant']
diseases = [e for e in kg.get('entities', []) if e.get('type') == 'Disease']

# Find disease associations
disease_associations = []
for variant in variants[:10]:
    variant_id = variant.get('name') or variant.get('id')
    if not variant_id:
        continue
    for disease in diseases[:5]:
        disease_id = disease.get('name') or disease.get('id')
        if not disease_id:
            continue
        path = graph_analyzer.connectivity_analyzer.calculate_shortest_paths(
            kg, source=variant_id, target=disease_id
        )
        if path.get('exists') and path.get('distance', -1) <= 2:
            disease_associations.append({
                'variant': variant_id,
                'disease': disease_id,
                'path_length': path.get('distance', -1),
                'confidence': variant.get('confidence', 1.0)
            })

disease_associations.sort(key=lambda x: x['confidence'], reverse=True)

print(f"Top disease associations:")
for i, assoc in enumerate(disease_associations[:5], 1):
    print(f"{i}. {assoc['variant']} -> {assoc['disease']} (path length: {assoc['path_length']}, confidence: {assoc['confidence']:.3f})")

Visualizing the Temporal Knowledge Graph

Generate an interactive visualization of the temporal genomic knowledge graph.

In [ ]:
from semantica.visualization import TemporalVisualizer

# Visualize temporal dashboard
temporal_viz = TemporalVisualizer()
fig = temporal_viz.visualize_temporal_dashboard(
    kg,
    output="interactive"
)

# Display the figure
fig.show() if fig else None

Exporting Results

Export the temporal knowledge graph to various formats for further analysis or integration with other tools.

In [ ]:
from semantica.export import GraphExporter

exporter = GraphExporter()
exporter.export(kg, output_path="genomic_variant_kg.json", format="json")
exporter.export(kg, output_path="genomic_variant_kg.graphml", format="graphml")

print("Exported knowledge graph to JSON and GraphML formats")