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

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.

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.

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.

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