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.

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.

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.

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.