# Bio Applied Capstone Project

> BLAST-identify unknown DNA/CDS with Biopython, QC/translate sequences, build NJ/UPGMA trees, and scan protein motifs. Use for sequence-to-discovery capstones, unknown-sequence ID, or FASTA-BLAST-tree pipelines.

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

---


# From Sequence to Discovery: Integrative Bioinformatics Capstone

## When to Use

- Identifying a batch of unknown DNA/CDS sequences and building an evidence chain from BLAST hit to phylogenetic placement
- Building NJ/UPGMA trees from a small set of candidate/orthologous sequences and comparing topologies
- Scanning translated proteins for a known short motif (e.g. CXXCH heme attachment) or measuring per-position conservation
- Running a full pipeline: QC -> BLAST ID -> translate -> tree -> motif scan -> GO/pathway context -> publication figure
- Estimating a synonymous/nonsynonymous (dN/dS-style) signal between two close CDS to argue purifying selection

## Version Compatibility

Python >= 3.10, biopython >= 1.81 (`Bio.Align.PairwiseAligner`, `Bio.Phylo.TreeConstruction`), pandas >= 2.0, matplotlib >= 3.7. NCBI BLAST web API (`Bio.Blast.NCBIWWW.qblast`) as of 2024-2025.

## Prerequisites

- `pip install biopython pandas matplotlib numpy`
- Internet access for `Bio.Entrez` / `Bio.Blast.NCBIWWW` calls; set `Entrez.email` before any NCBI query
- Familiarity with FASTA/CDS basics, BLAST output, and distance-based tree construction (see `bio-sequence-io-read-sequences`, `bio-database-access-blast-searches`, `bio-phylogenetics-tree-io`)

## Pipeline Overview

```text
Unknown DNA -> QC/clean -> BLAST ID -> translate -> MSA/tree -> motif/structure -> GO/pathway -> figure
```

## Step 1-2: QC, Translate, and BLAST-Identify

**Goal:** trim sequencing artifacts, verify the reading frame, translate to protein, and confirm identity for at least a few representatives via real BLAST (don't assume the rest by similarity alone).
**Approach:** strip leading/trailing `N`, check `ATG` start and length % 3 == 0 before translating; BLAST only 2-3 representatives (NCBI queries take 1-5 min) and cache the XML.

```python
from Bio.Seq import Seq
from Bio.SeqUtils import gc_fraction
from Bio.Blast import NCBIWWW, NCBIXML
import pandas as pd

def qc_and_translate(raw_sequences: dict) -> tuple[dict, pd.DataFrame]:
    """Trim N's, validate ORF frame, translate to protein.

    raw_sequences: {sample_id: raw_dna_string}
    Returns (protein_seqs, qc_dataframe).
    """
    cleaned, proteins, rows = {}, {}, []
    for name, seq_str in raw_sequences.items():
        raw = seq_str.upper()
        stripped = raw.strip('N')                    # drop sequencing-artifact N runs
        clean = Seq(stripped)
        gc = gc_fraction(clean) * 100
        starts_atg = str(clean)[:3] == 'ATG'
        trim_len = len(clean) - (len(clean) % 3)      # force in-frame length
        proteins[name] = str(clean[:trim_len].translate())
        cleaned[name] = clean
        rows.append({'sample': name, 'raw_len': len(raw), 'clean_len': len(clean),
                     'GC%': round(gc, 1), 'starts_ATG': starts_atg,
                     'in_frame': len(clean) % 3 == 0})
    qc_df = pd.DataFrame(rows).set_index('sample')
    flagged = qc_df[(~qc_df['starts_ATG']) | (~qc_df['in_frame'])]
    if len(flagged):
        print(f"Flag for manual review before translating: {list(flagged.index)}")
    return proteins, qc_df


def blast_identify(seq: str, program: str = 'blastn', database: str = 'nt',
                    top_n: int = 3, cache_path: str = 'blast_cache.xml') -> list[dict]:
    """BLAST one sequence against NCBI and return the top hits.

    Real API is Bio.Blast.NCBIWWW.qblast(program, database, sequence) --
    there is no Bio.Blast.blast() convenience function (a common LLM
    hallucination). Always cache the XML so re-parsing doesn't re-query.
    """
    result_handle = NCBIWWW.qblast(program, database, seq)
    with open(cache_path, 'w') as out:
        out.write(result_handle.read())
    with open(cache_path) as f:
        record = NCBIXML.read(f)
    hits = []
    for alignment in record.alignments[:top_n]:
        hsp = alignment.hsps[0]
        hits.append({'title': alignment.title, 'e_value': hsp.expect,
                     'pct_identity': 100 * hsp.identities / hsp.align_length})
    return hits
```

## Step 3-4: Distance Trees (NJ / UPGMA)

**Goal:** build and compare Neighbor-Joining and UPGMA trees from a pairwise identity distance matrix.
**Approach:** for equal-length, indel-free CDS an index-by-index identity distance is a valid shortcut; for real orthologs with indels, align first (MUSCLE/Clustal Omega via `Bio.Align.Applications`, or `mafft --auto`) and run `DistanceCalculator('identity')` on the resulting `MultipleSeqAlignment` instead.

```python
from Bio.Phylo.TreeConstruction import DistanceMatrix, DistanceTreeConstructor

def build_trees(seqs: dict[str, str]):
    """Build NJ and UPGMA trees from an identity-based distance matrix.

    seqs: {name: sequence}, all sequences assumed equal length / gap-free.
    Returns (nj_tree, upgma_tree).
    """
    def pairwise_distance(s1: str, s2: str) -> float:
        n = min(len(s1), len(s2))
        matches = sum(1 for i in range(n) if s1[i] == s2[i])
        return 1.0 - matches / n

    names = list(seqs.keys())
    matrix = [[0.0 if i == j else pairwise_distance(seqs[names[i]], seqs[names[j]])
               for j in range(i + 1)] for i in range(len(names))]
    dm = DistanceMatrix(names, matrix)
    constructor = DistanceTreeConstructor()
    return constructor.nj(dm), constructor.upgma(dm)
```

UPGMA assumes a molecular clock (constant substitution rate, forces an ultrametric tree); NJ does not. Compare both topologies before trusting either branch length as literal divergence time.

## Step 5-6: Motif Scan and dN/dS Signal

**Goal:** locate a short functional motif (e.g. the CXXCH heme-attachment signature) in translated proteins, and classify DNA differences between the two most divergent sequences as synonymous/nonsynonymous.
**Approach:** use a regex/PROSITE-style pattern for the motif (BLAST is overkill for a short positional signature); classify codon-by-codon by re-translating each codon pair.

```python
import re
from collections import Counter
from Bio.Seq import Seq

def find_motif_and_dnds(protein_seqs: dict[str, str], seq_a: str, seq_b: str) -> dict:
    """Scan a CXXCH-style motif in each protein, then classify DNA differences
    between two same-frame CDS as synonymous/nonsynonymous per codon.

    protein_seqs: {species: protein_string}; seq_a/seq_b: two in-frame CDS strings.
    """
    for species, prot in protein_seqs.items():
        m = re.search(r'C.{2}CH', prot)
        print(f"{species}: {'FOUND at ' + str(m.start() + 1) if m else 'NOT FOUND'}")

    def classify(c1: str, c2: str) -> str:
        if c1 == c2:
            return 'same'
        aa1, aa2 = str(Seq(c1).translate()), str(Seq(c2).translate())
        return 'syn' if aa1 == aa2 else 'nonsyn'

    counts = Counter()
    for k in range(min(len(seq_a), len(seq_b)) // 3):
        counts[classify(seq_a[3 * k:3 * k + 3], seq_b[3 * k:3 * k + 3])] += 1
    ratio = counts['nonsyn'] / counts['syn'] if counts['syn'] else float('inf')
    print(f"syn={counts['syn']} nonsyn={counts['nonsyn']} nonsyn/syn={ratio:.3f}")
    return dict(counts)
```

A nonsyn/syn ratio well below 1, concentrated at codon position 3 (wobble), is the classic signature of purifying selection on a conserved protein. For structure, load a real reference (e.g. PDB `1YCC`, yeast iso-1-cytochrome c) with `Bio.PDB.PDBParser` and confirm residue identities computationally rather than citing a remembered residue number.

## Step 7: GO/Pathway Enrichment and Figures

Run enrichment on the identified gene set with `goatools` or the Enrichr API (see `bio-pathway-analysis-go-enrichment`), applying Benjamini-Hochberg FDR correction since you're testing thousands of terms. Assemble a multi-panel `matplotlib` figure (QC table, trees, variation map, motif hits) as the capstone deliverable.

## Pitfalls

- `NCBIWWW.qblast(program, database, sequence)` is the real signature — not a `Bio.Blast.blast()` convenience function, and not `qblast(sequence, program, database)`; check argument order against actual Biopython docs.
- BLASTing only 2-3 representatives and inferring the rest "by alignment similarity" risks silently propagating a paralog/pseudogene/contaminant misassignment onto sequences that were never independently confirmed.
- The identity-distance shortcut (`1 - matches/len`) only works index-by-index on equal-length, indel-free sequences; real orthologs need a true MSA first, or the naive comparison silently misaligns everything downstream.
- UPGMA's ultrametric assumption (molecular clock) can differ from NJ's topology — report both, don't cherry-pick the one matching a prior expectation.
- Verify any "known" motif/residue number (CXXCH, Met80, etc.) against the actual structure/sequence file, not from memory — residue numbering off-by-one and wrong-chain citations are common.
- Multiple testing: GO/pathway enrichment needs FDR (Benjamini-Hochberg) correction across all tested terms, not raw p-values.

## See Also

- `bio-database-access-blast-searches`
- `bio-phylogenetics-modern-tree-inference`
- `bio-pathway-analysis-go-enrichment`
- `bio-structural-biology-structure-io`