# Python Bio SQL For Bioinformatics

> Query genomic data with SQLite (sqlite3, pandas.read_sql_query): SELECT/WHERE/GROUP BY/HAVING, INNER/LEFT JOIN, subqueries, parameterized inserts. Use when writing SQL over gene/variant tables or building a SQLite database.

- Skill: `pavel-kravchenko/python-bio-sql-for-bioinformatics` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/python-bio-sql-for-bioinformatics`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/python-bio-sql-for-bioinformatics/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Data & Analytics
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/python-bio-sql-for-bioinformatics

---


# SQL for Bioinformatics

## When to Use

- Storing gene annotations, variant calls, or expression results in a relational schema instead of flat files.
- Writing SQL joins/aggregations to answer questions like "genes with pathogenic variants AND high tumor expression."
- Building a local queryable cache of data pulled from Ensembl, UCSC, NCBI, or dbSNP (all backed by relational DBs).
- Loading a pandas DataFrame from a SQL query (`pd.read_sql_query`) for downstream plotting/stats.
- Avoiding SQL-injection bugs when a query needs to include a user- or pipeline-supplied gene name or position.

## Version Compatibility

- Python ≥3.10, `sqlite3` (stdlib, no install needed), pandas ≥2.0.
- Same SQL patterns apply to PostgreSQL/MySQL via `sqlalchemy` + `pd.read_sql_query`, with minor dialect differences (e.g. `SERIAL` vs `INTEGER PRIMARY KEY AUTOINCREMENT`).

## Prerequisites

- `pip install pandas numpy` (sqlite3 ships with Python).
- Basic familiarity with relational schema design (primary/foreign keys) and pandas DataFrames.

## Key Concepts

**JOIN type determines which rows survive:** `INNER JOIN` — only matched rows; `LEFT JOIN` — all left rows, NULL for unmatched right side.

**`HAVING` vs `WHERE`:** `WHERE` filters rows before grouping; `HAVING` filters groups after aggregation. `WHERE COUNT(*) > 3` is invalid SQL — use `HAVING`.

**Subqueries vs JOINs:** A JOIN is usually more readable and faster. A subquery with `IN (SELECT ...)` is clearer when you only need membership in a small inner set, not its columns.

**Always use parameterized queries:** Never build SQL with f-strings/`.format()` containing external input — that's SQL injection. Use `cursor.execute("... WHERE gene = ?", (gene_name,))`.

**Goal:** build a small relational schema (genes, variants, expression, pathways) and populate it so it can be queried like a mini Ensembl/dbSNP mirror.

**Approach:** create tables with `executescript`, then bulk-load rows with `executemany` and parameterized `?` placeholders.

```python
import sqlite3
import pandas as pd
import numpy as np


def build_demo_db() -> sqlite3.Connection:
    """Create an in-memory SQLite DB with genes/variants/expression/pathway tables."""
    conn = sqlite3.connect(":memory:")
    conn.executescript("""
        CREATE TABLE genes (
            gene_id    INTEGER PRIMARY KEY,
            symbol     TEXT NOT NULL,
            chromosome TEXT,
            start_pos  INTEGER,
            end_pos    INTEGER,
            biotype    TEXT
        );
        CREATE TABLE variants (
            variant_id  INTEGER PRIMARY KEY,
            gene_id     INTEGER REFERENCES genes(gene_id),
            position    INTEGER,
            ref_allele  TEXT,
            alt_allele  TEXT,
            clinical_significance TEXT
        );
        CREATE TABLE expression (
            expr_id   INTEGER PRIMARY KEY,
            gene_id   INTEGER REFERENCES genes(gene_id),
            tissue    TEXT,
            condition TEXT,
            tpm       REAL
        );
        CREATE TABLE pathways (
            pathway_id   INTEGER PRIMARY KEY,
            pathway_name TEXT
        );
        CREATE TABLE gene_pathway (
            gene_id    INTEGER REFERENCES genes(gene_id),
            pathway_id INTEGER REFERENCES pathways(pathway_id)
        );
    """)

    genes = [
        (1, 'BRCA1', 'chr17', 43044295, 43125483, 'protein_coding'),
        (2, 'TP53',  'chr17',  7661779,  7687538, 'protein_coding'),
        (3, 'EGFR',  'chr7',  55019017, 55207337, 'protein_coding'),
        (4, 'MYC',   'chr8', 127735434, 127742951, 'protein_coding'),
        (5, 'KRAS',  'chr12', 25204789, 25250936, 'protein_coding'),
        (6, 'PTEN',  'chr10', 89692905, 89728532, 'protein_coding'),
        (7, 'RB1',   'chr13', 47775885, 47954065, 'protein_coding'),
    ]
    conn.executemany("INSERT INTO genes VALUES (?,?,?,?,?,?)", genes)

    variants = [
        (1, 1, 43045629, 'A', 'T', 'pathogenic'),
        (2, 2, 7674220,  'C', 'T', 'pathogenic'),
        (3, 3, 55181320, 'G', 'A', 'likely_pathogenic'),
        (4, 5, 25245347, 'G', 'T', 'pathogenic'),
        (5, 6, 89711933, 'T', 'A', 'benign'),
    ]
    conn.executemany("INSERT INTO variants VALUES (?,?,?,?,?,?)", variants)

    rng = np.random.default_rng(42)
    expr_rows, eid = [], 1
    for gid in range(1, 8):
        for tissue in ['liver', 'kidney', 'brain']:
            for condition in ['normal', 'tumor']:
                base = rng.uniform(10, 200)
                tpm = round(base * (1.8 if condition == 'tumor' else 1.0) + rng.normal(0, 5), 2)
                expr_rows.append((eid, gid, tissue, condition, max(tpm, 0.1)))
                eid += 1
    conn.executemany("INSERT INTO expression VALUES (?,?,?,?,?)", expr_rows)

    conn.executemany("INSERT INTO pathways VALUES (?,?)",
                      [(1, 'DNA repair'), (2, 'Cell cycle'), (3, 'Apoptosis')])
    conn.executemany("INSERT INTO gene_pathway VALUES (?,?)",
                      [(1, 1), (2, 1), (2, 2), (4, 2), (3, 3), (2, 3)])
    conn.commit()
    return conn


conn = build_demo_db()
```

**Goal:** filter, aggregate, and join across the schema — the operations that make SQL worth using over flat-file grep.

**Approach:** push filtering/grouping into the database with `pd.read_sql_query`, and let `HAVING` filter on the aggregate rather than the raw column.

```python
def query_examples(conn: sqlite3.Connection) -> dict[str, pd.DataFrame]:
    """Run representative SELECT/JOIN/subquery patterns and return each result as a DataFrame."""
    results = {}

    # WHERE + ORDER BY on a computed column
    results["long_genes_chr17"] = pd.read_sql_query("""
        SELECT symbol, chromosome, (end_pos - start_pos) AS length
        FROM genes
        WHERE chromosome = 'chr17' AND (end_pos - start_pos) > 20000
        ORDER BY length DESC
    """, conn)

    # GROUP BY + HAVING: filters on the aggregate, not the raw rows
    results["high_tumor_expr"] = pd.read_sql_query("""
        SELECT g.symbol, ROUND(AVG(e.tpm), 2) AS avg_tumor_tpm
        FROM genes g
        JOIN expression e ON g.gene_id = e.gene_id
        WHERE e.condition = 'tumor'
        GROUP BY g.symbol
        HAVING AVG(e.tpm) > 50
        ORDER BY avg_tumor_tpm DESC
    """, conn)

    # LEFT JOIN: keep genes with zero variants (INNER JOIN would drop them)
    results["variant_counts"] = pd.read_sql_query("""
        SELECT g.symbol, COUNT(v.variant_id) AS n_variants
        FROM genes g
        LEFT JOIN variants v ON g.gene_id = v.gene_id
        GROUP BY g.symbol
        ORDER BY n_variants DESC
    """, conn)

    # Subquery: genes that satisfy two independent conditions
    results["tumor_pathogenic"] = pd.read_sql_query("""
        SELECT symbol FROM genes
        WHERE gene_id IN (
            SELECT gene_id FROM expression
            WHERE condition = 'tumor'
            GROUP BY gene_id
            HAVING AVG(tpm) > 50
        )
        AND gene_id IN (
            SELECT gene_id FROM variants WHERE clinical_significance = 'pathogenic'
        )
    """, conn)

    return results


tables = query_examples(conn)
print(tables["tumor_pathogenic"])
```

**Goal:** write results back to the database and read them out with a safe, parameterized filter.

**Approach:** `CREATE TABLE IF NOT EXISTS`, bulk `executemany`, then `execute` with a `?` placeholder — never an f-string — for the caller-supplied value.

```python
def save_de_results(conn: sqlite3.Connection, de_data: list[tuple]) -> None:
    """Persist differential-expression results (gene_id, log2fc, pvalue, padj, significant)."""
    conn.execute("""
        CREATE TABLE IF NOT EXISTS de_results (
            gene_id INTEGER REFERENCES genes(gene_id),
            log2fc REAL, pvalue REAL, padj REAL, significant INTEGER
        )
    """)
    conn.executemany("INSERT INTO de_results VALUES (?,?,?,?,?)", de_data)
    conn.commit()


def fetch_gene_variants(conn: sqlite3.Connection, gene_symbol: str) -> pd.DataFrame:
    """Look up variants for one gene safely (parameterized -- never f-string the symbol in)."""
    return pd.read_sql_query(
        "SELECT v.* FROM variants v JOIN genes g ON v.gene_id = g.gene_id WHERE g.symbol = ?",
        conn, params=(gene_symbol,)
    )


save_de_results(conn, [(1, 2.3, 0.001, 0.01, 1), (2, 1.8, 0.005, 0.03, 1)])
assert not fetch_gene_variants(conn, "BRCA1").empty
conn.close()
```

## Pitfalls

- **JOIN type**: INNER loses unmatched rows; LEFT preserves them — choose deliberately, especially when counting "genes with 0 variants."
- **HAVING vs WHERE**: putting an aggregate in `WHERE` raises an error; use `HAVING` after `GROUP BY`.
- **SQL injection**: f-string/`.format()` SQL with external input is a real vulnerability — always use `?` placeholders and `params=`.
- **Off-by-one coordinates**: Python ranges are half-open `[start, stop)`; genomic coordinates (GFF/VCF) are often 1-based, BED is 0-based — check before comparing across formats.
- **In-memory DB scope**: `sqlite3.connect(":memory:")` is per-connection; a second `connect(":memory:")` call gets an empty, unrelated database.

## See Also

- `bio-expression-matrix-counts-ingest` — loading count/TPM matrices before pushing them into SQL tables.
- `bio-variant-calling-vcf-basics` — parsing VCF into the tabular form used by the `variants` table here.
- `bio-database-access-entrez-fetch` — pulling gene/variant records from NCBI to populate a local DB.
- `polars` — a faster DataFrame-native alternative to SQL joins for larger local datasets.