--- name: mims-harvard/tooluniverse-variant-to-mechanism source: https://app.decimal.ai/s/mims-harvard-tooluniverse-variant-to-mechanism@1/SKILL.md source_sha256: d7279d444581 --- # Variant-to-Mechanism Analysis Skill Trace the full causal chain from a genetic variant to its disease mechanism: regulatory context, target gene(s), molecular pathways, and phenotypic consequences. Integrates 7+ databases across 3 evidence layers (regulatory, molecular, disease) to build an evidence-graded mechanistic model. **IMPORTANT**: Always use English terms in tool calls. Respond in the user's language. --- ## LOOK UP, DON'T GUESS When uncertain about any scientific fact, SEARCH databases first (PubMed, UniProt, ChEMBL, ClinVar, etc.) rather than reasoning from memory. A database-verified answer is always more reliable than a guess. --- ## When to Use This Skill Apply when users: - Ask "how does rs7903146 cause type 2 diabetes?" - Want to trace a GWAS variant to its biological mechanism - Need to connect a non-coding variant to downstream pathways - Ask "what gene does this variant affect and through what pathway?" - Want a mechanistic narrative for a regulatory variant - Need to identify the causal gene for a GWAS locus - Ask "what is the functional impact of this intronic SNP?" **NOT for** (use other skills instead): - Coding variant pathogenicity / ACMG classification -> Use `tooluniverse-variant-interpretation` - Pure regulatory element annotation without mechanism -> Use `tooluniverse-regulatory-genomics` - Pure GWAS hit listing without mechanism -> Use `tooluniverse-gwas-snp-interpretation` - Pharmacogenomic variant annotation -> Use `tooluniverse-pharmacogenomics` - Gene-disease association without variant context -> Use `tooluniverse-gene-disease-association` --- ## Input Parameters | Parameter | Required | Description | Example | |-----------|----------|-------------|---------| | **variant** | Yes | rsID or genomic coordinates | `"rs7903146"` or `"10:112998590:C:T"` | | **trait** | No | Disease/trait context (helps prioritize) | `"type 2 diabetes"` | | **tissue** | No | Tissue of interest for eQTL/expression | `"pancreas"` | --- ## Workflow Overview **Phase 1** Variant Characterization (VEP, population frequency, CADD) -> **Phase 2** Regulatory Context (GWAS, RegulomeDB, ENCODE, cCREs) -> **Phase 3** Target Gene Identification (GTEx eQTL, OpenTargets L2G) -> **Phase 4** Gene Function & Pathways (STRING, Reactome, GO) -> **Phase 5** Disease Connection (OpenTargets, GenCC, DisGeNET) -> **Phase 6** Mechanistic Synthesis (causal chain, evidence grading, confidence scoring) --- ## Phase 1: Variant Characterization Resolve variant identity and get basic annotations. ```python from tooluniverse import ToolUniverse tu = ToolUniverse() tu.load_tools() # Step 1: VEP annotation vep = tu.tools.EnsemblVEP_annotate_rsid(variant_id="rs7903146") # NOTE: param is `variant_id`, NOT `rsid` # Response format variable: list, {data, metadata}, or {error} -- handle all three # Extract: consequence_type, gene, chromosome, position, alleles # Step 2: Population frequency and pathogenicity scores myvar = tu.tools.MyVariant_query_variants( query="rs7903146", # NOTE: param is `query`, NOT `variant_id` fields="dbsnp,gnomad_genome,cadd,clinvar" ) # Returns: {hits: [{_id, dbsnp, gnomad_genome: {af: {af, ...ancestry}}, cadd: {phred}, clinvar}]} # gnomAD AF: frequency in general population # CADD phred: >=20 = top 1% deleterious, >=30 = top 0.1% # Step 3: Confirm gene context gwas_snp = tu.tools.gwas_search_snps(rs_id="rs7903146") # Returns: SNP location, mapped genes, functional class # Optional one-call shortcut (when a GRCh38 coordinate is available): favor = tu.tools.FAVOR_annotate_variant(variant="19-44908822-C-T") # Returns freq (BRAVO/gnomAD-by-ancestry/1000G) + CADD/SIFT/PolyPhen/AlphaMissense # + conservation + ClinVar + regulatory annotation in a single call — a fast way to # populate most of Phase 1 (and the regulatory block of Phase 2) before drilling in. ``` --- ## Phase 2: Regulatory Context Assess the variant's regulatory environment. ### 2a: GWAS Associations ```python # All trait associations for this variant gwas = tu.tools.gwas_search_associations(rs_id="rs7903146", size=50) # Returns: {data: [{disease_trait, p_value, or_per_copy_number, beta, study_accession}]} # NOTE: Use gwas_search_associations, NOT gwas_get_associations_for_trait (BROKEN) ``` ### 2b: RegulomeDB Score ```python regulome = tu.tools.RegulomeDB_query_variant(rsid="rs7903146") # Returns: {status, data: {score/ranking, probability, features}} # Score 1a-2a = strong regulatory evidence # Uses GRCh38 (NOT hg19) ``` ### 2c: ENCODE Chromatin Context ```python # Check enhancer (H3K27ac), promoter (H3K4me3), and accessibility in disease-relevant tissue encode_h3k27ac = tu.tools.ENCODE_search_histone_experiments( histone_mark="H3K27ac", biosample_term_name="pancreas", limit=10) encode_h3k4me3 = tu.tools.ENCODE_search_histone_experiments( histone_mark="H3K4me3", biosample_term_name="pancreas", limit=10) encode_atac = tu.tools.ENCODE_search_chromatin_accessibility( biosample_term_name="pancreas", limit=10) ``` ### 2d: cCRE Overlap (if coordinates known) ```python # Check candidate cis-regulatory elements ccres = tu.tools.UCSC_get_encode_cCREs( chrom="chr10", start=112998000, end=112999000 ) # Returns: cCRE type (PLS=promoter, pELS=proximal enhancer, dELS=distal enhancer, CTCF-only) ``` ### 2e: OLS Trait Resolution (if trait provided) ```python # Resolve trait to ontology ID for downstream tools ols = tu.tools.ols_search_terms(query="type 2 diabetes", ontology="efo") # Returns: {status, data: [{iri, label, short_form, ontology_name}]} # For OpenTargets: use MONDO_0005148 for T2D (NOT EFO_0001360) ``` --- ## Phase 3: Target Gene Identification Identify which gene(s) the variant regulates. This is the critical link. ### 3a: GTEx eQTLs (primary evidence) ```python # Find eQTL associations for the nearest gene eqtls = tu.tools.GTEx_query_eqtl(gene_symbol="TCF7L2", size=100) # Returns: {status, data: [{snpId, geneSymbol, tissueSiteDetailId, pValue, nes}]} # Filter for the specific rsID in the results # NES > 0: alt allele increases expression; NES < 0: decreases # Check expression levels to prioritize tissues expression = tu.tools.GTEx_get_median_gene_expression(gene_symbol="TCF7L2") # Returns: median TPM across GTEx tissues # IMPORTANT: GTEx uses v8 data (v10 endpoints may return empty) ``` ### 3b: OpenTargets Locus-to-Gene (L2G) ```python # Get credible sets and L2G predictions credible = tu.tools.OpenTargets_get_variant_credible_sets( variant_id="10_112998590_C_T" # chr_pos_ref_alt format ) # Returns: credible set membership, L2G scores for causal genes # L2G > 0.5: high confidence causal gene # L2G 0.1-0.5: moderate confidence ``` ### 3c: Nearest Gene Mapping ```python # Get gene details for the candidate target gene_info = tu.tools.MyGene_query_genes( query="symbol:TCF7L2", species="human", fields="symbol,ensembl.gene,entrezgene,name,summary,go", size=5 ) # Extract Ensembl ID for OpenTargets queries # Filter by exact symbol match (first hit may not be correct gene) ``` ### How to Reason About Target Gene Identification The nearest gene is often NOT the causal gene for non-coding variants. Regulatory elements can act over hundreds of kilobases, skipping intervening genes entirely. Before looking at tools, reason about which gene is the likely target: - **eQTL in disease-relevant tissue is the strongest evidence.** If a variant is an eQTL for Gene X in pancreas and the trait is diabetes, that is a direct mechanistic link. An eQTL in an unrelated tissue is weaker but still informative. - **L2G scores integrate multiple signals** (distance, chromatin interaction, eQTL) into a single prediction. L2G > 0.5 is high confidence, but it is a computational prediction, not experimental proof. Combine it with eQTL data -- when both agree, confidence is high. - **Proximity alone is unreliable.** Many GWAS loci have the causal gene 2-3 genes away from the lead SNP. Never default to "nearest gene" without checking eQTL and L2G. - **When multiple candidate genes have evidence**, present all of them ranked by strength. The answer may genuinely be uncertain. --- ## Phase 4: Gene Function & Pathways Once target gene(s) identified, characterize molecular function. ```python # STRING PPI network — NOTE: param is `identifiers` (string), NOT `protein_ids` ppi = tu.tools.STRING_get_interaction_partners( identifiers="TCF7L2", species=9606, required_score=700) # STRING functional enrichment (GO, KEGG, Reactome) enrichment = tu.tools.STRING_get_functional_enrichment( identifiers="TCF7L2", species=9606) # Reactome pathway enrichment — space-separated STRING, NOT array reactome = tu.tools.ReactomeAnalysis_pathway_enrichment( identifiers="TCF7L2 CTNNB1 APC") # UniProt function (returns list of strings, NOT dict) uniprot_func = tu.tools.UniProt_get_function_by_accession(accession="Q9NQB0") # PANTHER GO enrichment — GO:0008150=BP, GO:0003674=MF, GO:0005575=CC panther = tu.tools.PANTHER_enrichment( gene_list="TCF7L2,CTNNB1,APC,LEF1", organism=9606, annotation_dataset="GO:0008150") ``` --- ## Phase 5: Disease Connection Link the molecular mechanism to disease phenotype. ```python # Disease associations for target gene (requires Ensembl ID) ot_diseases = tu.tools.OpenTargets_get_diseases_phenotypes_by_target_ensembl( ensemblId="ENSG00000148737") # Specific gene-disease evidence (both IDs must be pre-resolved) ot_evidence = tu.tools.OpenTargets_target_disease_evidence( ensemblId="ENSG00000148737", efoId="MONDO_0005148") # GenCC curated validity (Definitive/Strong/Moderate/Limited/Disputed/Refuted) gencc = tu.tools.GenCC_search_gene(gene_symbol="TCF7L2") # DisGeNET scored associations (requires DISGENET_API_KEY) disgenet = tu.tools.DisGeNET_search_gene(gene="TCF7L2", limit=20) # OpenTargets GWAS studies — use MONDO IDs, NOT EFO ot_gwas = tu.tools.OpenTargets_search_gwas_studies_by_disease( diseaseIds=["MONDO_0005148"], size=20) ``` --- ## Phase 6: Mechanistic Synthesis ### How to Build and Evaluate the Causal Chain The goal is an explicit chain: **Variant -> changed protein or expression -> altered pathway -> disease phenotype.** Each arrow is a link that needs evidence. A chain with all links supported is strong. A chain with a missing link is a hypothesis, not a conclusion. **Start with the variant's location.** Is this coding or non-coding? Coding variants directly change the protein -- a missense variant swaps one amino acid, a nonsense variant truncates it. Non-coding variants change REGULATION -- they affect how much protein is made, not what protein is made. This distinction determines your entire analysis path. **For coding variants**, ask: Is this position conserved across species? Is it in a functional domain (active site, binding interface, transmembrane helix)? A variant in a conserved active site residue is more likely pathogenic than one in a disordered loop. Check ClinVar and CADD for existing assessments, but reason about the structural context yourself. **For non-coding variants**, ask: Is there an eQTL linking this variant to a nearby gene? Is the variant in a TFBS, enhancer, or promoter (from RegulomeDB and ENCODE)? A variant in open chromatin with an eQTL in disease-relevant tissue has a clear regulatory mechanism. A variant in a gene desert with no eQTL is much harder to interpret -- acknowledge this uncertainty. **Evaluate each link independently.** For each arrow in the chain, ask: what is the evidence? Is it replicated? Is it from the right tissue? Then assess the chain as a whole: - All links well-supported with convergent evidence from multiple databases = **Established** mechanism - Most links strong but one relies on a single source = **Strong** mechanism - Evidence is mixed or only moderate across links = **Moderate** mechanism - One or more links have no direct evidence, only inference = **Preliminary** mechanism - The connection rests on proximity or text-mining alone = **Speculative** mechanism **The evidence hierarchy matters.** Clinical observation (ClinVar pathogenic with clinical data) is stronger than functional assays (eQTL, chromatin marks), which are stronger than computational predictions (CADD score, L2G model). A CADD score alone is weak evidence. Multiple convergent computational predictions are moderate. An eQTL replicated across studies in disease-relevant tissue is strong. **Always consider alternative mechanisms.** If eQTL points to Gene A but L2G favors Gene B, present both models. If the variant could act through splicing rather than expression, note that. Honest uncertainty is more useful than false confidence. --- ## Evidence Grading Grade each piece of evidence by how directly it supports a mechanistic link: - **T1 (Strongest)**: Replicated GWAS association combined with functional validation -- an eQTL in the right tissue plus regulatory chromatin marks at the variant position. This is convergent evidence from independent sources. - **T2**: Curated biological knowledge -- the target gene is in a known disease-relevant pathway (Reactome, GenCC Definitive). This does not prove the variant acts through this pathway, but it makes the mechanism biologically plausible. - **T3**: Computational predictions -- L2G scores, PPI network enrichment, GO term enrichment. These are hypothesis-generating, not confirmatory. Useful when T1/T2 evidence is sparse. - **T4**: Single-source or indirect evidence -- nearest gene assignment, literature co-mention, text-mining associations. These should never be the sole basis for a mechanistic claim. --- ## Common Mistakes and Fallbacks **Parameter pitfalls** (correct param names are shown in code blocks above): - `EnsemblVEP_annotate_rsid` uses `variant_id`, NOT `rsid` - `MyVariant_query_variants` uses `query`, NOT `variant_id` - `ReactomeAnalysis_pathway_enrichment` takes space-separated STRING, NOT array - OpenTargets disease search uses `queryString`, NOT `query` - Use MONDO IDs (e.g., MONDO_0005148 for T2D), NOT EFO IDs in OpenTargets - `gwas_get_associations_for_trait` is BROKEN; use `gwas_search_associations` - ENCODE biosample names are tissues/cell lines, not disease names **When data is missing**, adapt the analysis: - No eQTL: check additional tissues, use L2G alone, note reduced confidence - No GWAS hits: variant may be sub-threshold or in LD with the true signal - No RegulomeDB data: check ENCODE directly for chromatin marks at the position - No L2G data: fall back to eQTL + nearest gene with appropriate caveats - Gene not in STRING/Reactome: use UniProt function + GO annotations - Multiple candidate genes: present all ranked by evidence, do not force a single answer --- ## Reasoning Through Common Scenarios ### Scenario 1: Non-coding variant with clear regulatory signal When VEP shows an intronic or intergenic variant, GWAS shows strong trait associations, and RegulomeDB returns a high score (1a-2a), the path is straightforward: the variant disrupts a regulatory element. The key question becomes WHICH GENE it regulates. Check eQTL in disease-relevant tissue first, then L2G. When both converge on the same gene, follow the pathway analysis through to disease. This is the "happy path" and yields Established confidence. ### Scenario 2: Non-coding variant with ambiguous target gene When the variant is intergenic with multiple nearby genes and eQTL/L2G evidence is weak or conflicting, resist the temptation to pick the nearest gene. Instead: check eQTLs for ALL candidate genes in the locus, compare L2G scores, and consider whether one candidate has known disease-relevant biology. Present multiple models ranked by evidence. Report as Preliminary or Moderate. ### Scenario 3: No functional signal at all When there is no eQTL, no L2G support, and minimal regulatory annotation, acknowledge the gap honestly. The variant may act through a mechanism not captured by current databases (cell-type-specific regulation, splicing QTL, 3D genome interactions). State what was checked, what was negative, and what experiments would resolve the question. This is a valid and useful finding. --- ## Programmatic Access (Beyond Tools) When ToolUniverse tools return partial annotations, aggregate from multiple sources directly: ```python import requests, pandas as pd rsid = "rs12345" # Ensembl VEP REST — functional consequence prediction vep = requests.get(f"https://rest.ensembl.org/vep/human/id/{rsid}?content-type=application/json").json() # GTEx eQTL — tissue-specific expression effects eqtl = requests.get(f"https://gtexportal.org/api/v2/association/singleTissueEqtl", params={"snpId": rsid, "datasetId": "gtex_v8"}).json() # ClinVar bulk (variant_summary.txt.gz, ~200MB, all variants) # Download once, filter locally: # df = pd.read_csv("https://ftp.ncbi.nlm.nih.gov/pub/clinvar/tab_delimited/variant_summary.txt.gz", sep="\t") # hits = df[df["RS# (dbSNP)"] == int(rsid.replace("rs",""))] # Combine into unified annotation DataFrame annotations = {"rsid": rsid, "vep_consequence": vep[0]["most_severe_consequence"] if vep else None, "eqtl_tissues": len(eqtl.get("singleTissueEqtl", [])), "sources_checked": ["VEP", "GTEx", "ClinVar"]} ``` See `tooluniverse-data-wrangling` skill for API patterns and error handling. --- ## Limitations - GTEx eQTLs are bulk tissue (v8 data); cell-type-specific and splicing effects may be missed. - L2G and CADD are computational predictions, not experimental proof. - GWAS associations are correlative; LD tagging means the lead SNP may not be the causal variant. - Long-range regulatory effects (>1Mb) and 3D genome interactions are poorly captured by current tools. - OpenTargets variant credible sets require chr_pos_ref_alt format; get alleles from VEP first.