How to Annotate BLAST Hits with Gene Ontology Terms in Python

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-11-11

By the end of this guide you will have a working Python script that takes a BLAST results table and produces a spreadsheet where every hit carries its Gene Ontology (GO) terms, term names, ontology branch (biological process, molecular function, or cellular component), and evidence codes. You will run BLAST yourself, parse the tabular output in Python, pull GO annotations from UniProt through its REST API, resolve GO identifiers to human-readable names from the GO Consortium's OBO file, and write a clean CSV you can drop into a notebook or a manuscript supplement.

This is a core pattern in python for bioinformatics. Python dominates this niche because the BLAST command-line toolkit writes machine-readable output, UniProt exposes a clean REST interface, and libraries like Biopython make the glue code short. Published BLAST automation tools follow the same shape. PARSID is a command-line Python script that automates blasting a sequence file against a reference library and filtering matches in five defined steps, writing the results to a spreadsheet [1]. BLAST-QC exists for the same reason, because the authors found that analyzing large BLAST result sets by hand is tedious and that few good parsers were available for download, so researchers were left writing their own [2]. This article gives you that parser and the annotation layer on top of it, in one afternoon of work.

What You Need Before You Start

Have the following installed and on hand:

• NCBI BLAST+ command-line tools (blastp, blastx, makeblastdb). On macOS with Homebrew: brew install blast. On Ubuntu: sudo apt install ncbi-blast+.
• Python 3.9 or newer.
• The Python packages biopython, pandas, requests, and obonet. Biopython provides sequence parsing and BLAST XML parsing, obonet parses OBO format ontology files, pandas handles the tabular work, and requests talks to UniProt. If you want a broader Biopython tutorial before starting, the official Biopython documentation covers SeqIO and the BLAST parsers in its first chapters.
• A protein query file in FASTA format (for example, proteins.faa).
• Either a local reference database to blast against, or access to swissprot for remote-style local searches.

Everything after Step 2 works from the BLAST output file alone, so you can also apply this workflow to BLAST tables you already have or tables exported from a web BLAST run with the "Tabular" output option.

Background: What a BLAST Hit Gives You and What a GO Annotation Adds

A BLAST hit is a similarity statement. It says that your query sequence aligns to a database sequence with a given percent identity, alignment length, and expectation value (E-value, the expected number of chance alignments with a score at least this good in a database of this size). Similarity is not function. A hit tells you where to look for functional evidence, and a GO annotation is that evidence, attached to the database entry by curators or inferred by automated pipelines.

The Gene Ontology is a controlled vocabulary of terms organized in a directed acyclic graph across three branches: biological process (for example, GO:0006915, apoptotic process), molecular function (for example, GO:0003700, DNA-binding transcription factor activity), and cellular component (for example, GO:0005634, nucleus). Each annotation in UniProt carries an evidence code. Experimental codes such as IDA (inferred from direct assay) are the strongest, while IEA (inferred from electronic annotation) means the term was assigned by automated methods and no curator has reviewed it. When you annotate your hits, you should always carry the evidence code into your output table, because a claim resting entirely on IEA terms is a weaker claim than one resting on IDA terms.

The retrieval logic of this workflow matters as much as the code. PARIGA, a server for filtering BLAST results, was built on the observation that as data volume grows, manual inspection becomes tedious and tools for parsing and filtering results by parameters such as reciprocal best hits are required [3]. You will build a small version of that idea: parse once, then filter by percent identity, coverage, and evidence code before you ever commit an annotation to your query.

A Step-by-Step Guide to Annotate BLAST Hits with Gene Ontology Terms in Python

This section is the complete workflow to annotate BLAST hits with Gene Ontology terms in Python, step by step. Each step ends with the expected output so you can confirm you are on track.

Step 1: Set Up the Environment

Create an isolated environment so the bioinformatics packages do not collide with anything else on your machine.

mamba create -n goannot -c conda-forge -c bioconda \
  blast biopython pandas obonet requests python=3.11
mamba activate goannot

If you prefer pip for the Python side after installing BLAST+ separately:

python -m venv .venv
source .venv/bin/activate
pip install biopython pandas requests obonet

Verify the tools:

blastp -version
python -c "import Bio, pandas, requests, obonet; print('ok')"

Expected output:

blastp: 2.15.0+
  Package: blast 2.15.0, build Jul  6 2023 21:16:05
ok

Any version of BLAST+ from the last several years supports the flags used below, and any Python from 3.9 onward runs the code unchanged.

Step 2: Run BLAST with Machine-Readable Output

The single most useful decision in this workflow is to make BLAST emit tabular output with exactly the columns you need, using -outfmt 6 plus a column list. Do not parse the human-readable default report. Do not parse plain HTML. Structured output at the source is the difference between a ten-line parser and a hundred-line one, and every published Python BLAST tool cited here assumes structured input: PARSID consumes BLAST output with cutoff parameters and writes a spreadsheet [1], BlaSTorage parses multi-file BLAST results into a queryable object database [4], and BLAST-QC reads raw BLAST output directly to produce condensed summaries [2].

Build a database if you are using your own reference sequences:

makeblastdb -in uniprot_sprot.fasta -dbtype prot -out db/uniprot_sprot

Then run the search:

blastp -query proteins.faa -db db/uniprot_sprot -out blast_hits.tsv \
  -outfmt "6 qseqid sseqid pident length qcovs evalue bitscore stitle" \
  -evalue 1e-5 -max_target_seqs 5 -num_threads 8

Each parameter matters:

• -outfmt 6 selects tabular output. The quoted column list controls which fields appear, in order.
• qseqid and sseqid are the query and subject identifiers, which you need to join hits to annotations.
• pident is percent identity across the aligned region.
• length is the alignment length in residues.
• qcovs is query coverage per subject, the fraction of your query covered by that single hit, reported as a percentage.
• evalue and bitscore rank hit quality.
• stitle carries the subject definition line, which often contains the protein name and, when you blast against Swiss-Prot, the full description.
• -evalue 1e-5 is a conventional screening cutoff for "this similarity is unlikely to be chance," used widely in practice and adjustable to your data.
• -max_target_seqs 5 limits reported hits per query.

A caution on -max_target_seqs: this flag affects the heuristic search process itself, not just the number of rows written, and its effects have been a subject of active discussion among BLAST users. Torkian and colleagues made exactly this point when justifying standalone parsing as the reliable path for large BLAST datasets [2]. The practical consequence is that the set of hits BLAST reports at a given -max_target_seqs value may not be the true top N hits. Report the parameter in your methods, and if a hit list drives an important conclusion, run a sensitivity check at a higher value.

If your subject database is Swiss-Prot, sseqid will look like sp|P04637|TP53. If it is NCBI nr or RefSeq, it will look like WP_012345678.1 or ref|NP_000537.3|. Step 4 handles both shapes.

Step 3: Parse the BLAST Table in Python

Load the table with pandas, naming the columns explicitly to match the -outfmt order.

import pandas as pd

COLS = ["qseqid", "sseqid", "pident", "length",
        "qcovs", "evalue", "bitscore", "stitle"]

hits = pd.read_csv("blast_hits.tsv", sep="\t", header=None, names=COLS)
print(hits.head())

Expected output (your values will differ):

    qseqid          sseqid  pident  length  qcovs       evalue  bitscore  \
0  prot_01   sp|P04637|TP53    99.4     393    393  0.000000e+00     782.1
1  prot_01   sp|P56537|MDM4     12.1      87     84  4.200000e-06      33.5
2  prot_02  sp|P24941|CDK2A    100.0     298    298  0.000000e+00     610.0

Now apply the filters that matter for annotation transfer. A common, defensible rule set is at least 30 percent identity, at least 70 percent query coverage, and alignment length at least 100 residues for protein queries. Adjust to your domain. This is also the point where more advanced filtering strategies from the literature become available. PARIGA lets users retrieve reciprocal best hits and cases where a query matches a target in multiple regions, using logical operations over stored results without re-running the search [3]. A reciprocal best hit filter, which keeps a hit only if the subject's own best match is the query, is a cheap way to enrich for orthologs before you transfer annotations.

filtered = hits[
    (hits["pident"] >= 30) &
    (hits["qcovs"] >= 70)
].copy()
print(len(hits), "raw hits,", len(filtered), "after filtering")

For very large result sets, tens of millions of lines and up, a plain pandas load can strain memory. BlaSTorage was designed for exactly this case, parsing multi-file BLAST results into a purpose-built object database with a complete API and benchmarked at speeds comparable to a C++ parser [4]. BLAST-QC similarly targets fast, workflow-friendly summarization of large BLAST outputs [2]. For most annotation projects under a few million hits, pandas with dtype={"evalue": "float64"} and engine="c" is enough.

Step 4: Extract UniProt Accessions from Subject IDs

GO annotations are keyed to database accessions. UniProt is the practical hub because it cross-references GO directly and exposes the mappings through a documented REST API. Extract the accession part of each sseqid:

import re

def to_uniprot_acc(sseqid: str) -> str | None:
    # UniProt FASTA headers look like sp|P04637|TP53 or tr|A0A123|A0A123_HUMAN
    m = re.match(r"^(?:sp|tr)\|([^|]+)\|", sseqid)
    if m:
        return m.group(1)
    # Already a bare accession
    if re.match(r"^[A-Z]\d[A-Z0-9]{3}\d$", sseqid):
        return sseqid
    return None

filtered["uniprot_acc"] = filtered["sseqid"].map(to_uniprot_acc)
print(filtered[["qseqid", "sseqid", "uniprot_acc"]].head())

Expected output:

    qseqid          sseqid uniprot_acc
0  prot_01   sp|P04637|TP53      P04637
1  prot_01   sp|P56537|MDM4      P56537
2  prot_02  sp|P24941|CDK2A      P24941

If your hits come from NCBI databases instead, the identifiers are NCBI accessions, not UniProt accessions. You have two clean options. First, run BLAST against Swiss-Prot or UniProtKB as the subject database, which is what the examples above assume. Second, keep the NCBI accessions and query NCBI E-utilities for the associated Gene records and their GO products, which is more code and more fragile. For GO annotation, blasting against UniProtKB is the simpler and more reproducible route. One more reason to prefer curated subject databases: sequences deposited in GenBank can carry outdated taxonomic identifications, which is the motivation behind curated reference libraries in automated identification workflows [1]. An outdated species label on your subject does not corrupt GO terms directly, but it corrupts any downstream interpretation that depends on which organism the function came from.

Step 5: Fetch GO Annotations from UniProt

UniProt serves each entry as JSON at a predictable URL. In the JSON, GO terms appear in the uniProtKBCrossReferences list with database == "GO", and each cross-reference carries Evidence and GoAspect properties. Fetch with a cache and a polite delay so you do not hammer the service.

import json
import time
import requests
from pathlib import Path

CACHE_FILE = Path("go_cache.json")
cache = json.loads(CACHE_FILE.read_text()) if CACHE_FILE.exists() else {}

def fetch_go(acc: str) -> list[dict]:
    """Return a list of GO terms with evidence codes for one accession."""
    if acc in cache:
        return cache[acc]
    url = f"https://rest.uniprot.org/uniprotkb/{acc}.json"
    r = requests.get(url, timeout=30)
    r.raise_for_status()
    data = r.json()
    terms = []
    for xref in data.get("uniProtKBCrossReferences", []):
        if xref.get("database") == "GO":
            props = {p["key"]: p["value"] for p in xref.get("properties", [])}
            terms.append({
                "go_id": xref["id"],
                "aspect": props.get("GoAspect", ""),   # P, F, or C
                "evidence": props.get("Evidence", ""), # e.g. IEA or IDA:xxx
            })
    cache[acc] = terms
    time.sleep(0.3)
    return terms

accs = sorted({a for a in filtered["uniprot_acc"] if a})
go_by_acc = {a: fetch_go(a) for a in accs}
CACHE_FILE.write_text(json.dumps(cache))

Expected content of go_cache.json for one accession, abbreviated:

{
  "P04637": [
    {"go_id": "GO:0005634", "aspect": "C", "evidence": "IDA"},
    {"go_id": "GO:0003700", "aspect": "F", "evidence": "IDA"},
    {"go_id": "GO:0006915", "aspect": "P", "evidence": "IDA"}
  ]
}

Cache to disk rather than to memory alone. UniProt occasionally merges or deletes accessions during its release cycle, and a cached response plus a log of failed accessions makes your run reproducible across weeks of work. If you prefer tabular responses or need to fetch tens of thousands of accessions in one request, the UniProt search endpoint accepts a query of OR-joined accessions with a fields=accession,go_id parameter and returns TSV, which you can regex-parse for GO: tokens. JSON is the more structured choice here, though parse cost varies by language and implementation, and a cross-language benchmark of native JSON parsers found measurable differences in parsing speed and memory use between implementations, including Python's [5]. For annotation-scale jobs the difference is irrelevant. For gigabyte-scale JSON, prefer UniProt's TSV or compressed streaming formats.

Step 6: Map GO IDs to Names with an OBO File

GO identifiers are stable, but readers want names. Download the GO Consortium's OBO file once and resolve every identifier locally. The go-basic.obo file is filtered to exclude terms you should not use directly, which makes it the right file for annotation work.

curl -O http://purl.obolibrary.org/obo/go/go-basic.obo

import obonet

graph = obonet.read_obo("go-basic.obo")

def term_name(go_id: str) -> str:
    node = graph.nodes.get(go_id)
    return node["name"] if node else go_id  # fall back to the raw ID

print(term_name("GO:0006915"))

Expected output:

apoptotic process

The OBO file is a few hundred megabytes in memory and takes a minute or two to parse, so parse it once and reuse the graph. Keep the file dated in your methods, because term names and parentage evolve with each GO release even though identifiers do not. This is the same principle that motivates ontology-aware data pipelines in other fields. The O-RAW radiomics workflow was built in Python specifically to attach a standardized lexicon to computed features and emit a W3C-compliant semantic graph, in order to make multi-center results comparable [6], and its technical report frames the same goal in terms of FAIR data principles, that is, data that are findable, accessible, interoperable, and reusable [7]. Your output CSV is a lightweight version of that idea: identifiers plus names plus evidence in one interoperable table.

Step 7: Build the Final Annotation Table

Explode the per-accession term lists into one row per query-hit-term, and write the deliverable.

ASPECT = {"P": "biological_process",
          "F": "molecular_function",
          "C": "cellular_component"}

rows = []
for _, hit in filtered.iterrows():
    acc = hit["uniprot_acc"]
    if not acc:
        continue
    for term in go_by_acc.get(acc, []):
        rows.append({
            "query": hit["qseqid"],
            "subject": hit["sseqid"],
            "subject_accession": acc,
            "pident": hit["pident"],
            "qcovs": hit["qcovs"],
            "evalue": hit["evalue"],
            "go_id": term["go_id"],
            "go_name": term_name(term["go_id"]),
            "aspect": ASPECT.get(term["aspect"], term["aspect"]),
            "evidence": term["evidence"],
        })

annot = pd.DataFrame(rows)
annot.to_csv("blast_go_annotations.csv", index=False)
print(annot.head())

Expected output:

    query       subject subject_accession  pident  qcovs        evalue  \
0  prot_01  sp|P04637|TP53            P04637    99.4    393  0.000000e+00
1  prot_01  sp|P04637|TP53            P04637    99.4    393  0.000000e+00
2  prot_01  sp|P04637|TP53            P04637    99.4    393  0.000000e+00
3  prot_02  sp|P24941|CDK2A            P24941   100.0    298  0.000000e+00

         go_id             go_name               aspect evidence
0  GO:0005634             nucleus  cellular_component      IDA
1  GO:0003700  DNA-binding transcription factor activity  molecular_function  IDA
2  GO:0006915    apoptotic process          biological_process   IDA
3  GO:0005634             nucleus  cellular_component      IDA

If you want one compact row per query instead of one row per term, group by query and join terms with a separator:

compact = (annot
           .groupby(["query", "subject"], as_index=False)
           .agg(go_terms=("go_id", lambda s: "|".join(sorted(set(s)))),
                n_terms=("go_id", "nunique"),
                best_pident=("pident", "max")))
compact.to_csv("blast_go_compact.csv", index=False)

A useful trick for annotation transfer is to restrict terms to the strongest experimental evidence before summarizing, so the headline table reflects curator-reviewed biology and the full table keeps everything:

experimental = annot[annot["evidence"].str.startswith(("IDA", "IMP", "IPI", "IGI", "IEP", "HDA", "HMP", "HGI", "HEP"))]

Parameters That Matter, Summarized

How to Check That Your Annotations Are Correct

Verify the pipeline at three levels before trusting the output.

First, verify the join. Pick one query with a confident hit, open its subject accession on the UniProt website, and confirm that the GO terms in your CSV match the "Gene Ontology" section of the web entry, including the evidence tags. Any mismatch means the accession extraction or the JSON path is wrong.

Second, verify the identifier hygiene. Run two asserts:

import re
assert annot["go_id"].str.fullmatch(r"GO:\d{7}").all(), "malformed GO ID found"
assert annot["subject_accession"].notna().all(), "unparsed accession found"
print(annot["evidence"].value_counts())

The evidence distribution is a quick reality check. If nearly everything is IEA, your annotations are automated inferences and you should say so in any figure caption. If a curated database returns nothing but IEA, suspect that you are hitting the wrong accession namespace.

Third, verify the biology at the extremes. Queries annotated to nucleus, cytoplasm, and a handful of process terms that agree with the subject protein's known name are reassuring. A query whose best hit is a transmembrane receptor but whose top annotation is a purely nuclear process usually signals a low-quality hit that slipped past your filters, and rereading its pident and qcovs will typically show why.

Common Mistakes and Limitations

These are the errors that show up repeatedly when people annotate BLAST hits with GO terms in Python, with the exact failure you will see and the fix.

• Parsing the wrong output format. If you saved BLAST output with the default pairwise formatting and try to parse it as a table, pandas raises pandas.errors.EmptyDataError: No columns to parse from file or produces a single garbage column. Rerun BLAST with -outfmt 6 and a column list. The reverse also bites: Biopython's NCBIXML parser expects -outfmt 5 XML and raises a parse error on tabular files.
• Forgetting that -max_target_seqs shapes the search. The parameter influences the BLAST heuristic, not merely the report length, a point raised explicitly in the BLAST-QC paper's discussion of BLAST parameters [2]. Hits at a low setting may not be the true best N. Report the value and re-run at a higher setting when a specific hit list matters.
• Treating hit identity as query identity. A GO term on a subject sequence is evidence about the subject. Transferring it to your query is an inference whose strength depends on percent identity, coverage, and whether the hit is an ortholog. Filters in Step 3 make the transfer defensible, reciprocal best hit logic makes it stronger, and neither makes it a measurement. PARIGA's reciprocal best hit operations exist precisely to support this kind of claim construction [3].
• Mixing accession namespaces. An NCBI WP_ accession sent to the UniProt REST endpoint returns requests.exceptions.HTTPError: 404 Client Error: Not Found for a reason that looks mysterious if you forgot the namespace mismatch. Either blast against UniProtKB as in this guide, or use a proper cross-reference service between the namespaces. Do not guess.
• Ignoring evidence codes. A table full of IEA terms is a table of computational predictions. Keep the column, summarize experimental codes separately, and let readers see the difference.
• Using a stale OBO file. GO terms become obsolete and merge into replacements. obonet will not contain obsolete identifiers removed from go-basic.obo, and your term_name fallback then silently writes the raw ID. That fallback is the correct behavior, but check the counts: a large number of raw-ID fallbacks means your OBO file predates your annotation source or vice versa.
• Fetching without a cache. A re-run against tens of thousands of accessions with no disk cache wastes an hour and adds avoidable load to UniProt's service. Write the cache file on every increment.
• Misreading GenBank subject metadata. Subject titles in GenBank-derived databases can carry outdated organism identifications, which is a documented motivation for curated reference libraries in automated BLAST workflows [1]. If your downstream claim depends on taxonomy, verify it rather than trusting the stitle string.
• One-size-fits-all cutoffs. The identity and coverage thresholds in Step 3 are starting points. Viral and highly diverged protein families need looser, structurally informed criteria, and shallow coverage can still be meaningful for domain-level function. State your thresholds and why.
• Over-automating the interpretation. PARSID's authors designed their five-step workflow with cutoff parameters and supplementary metadata in a CSV precisely so routine decisions are explicit rather than buried in code [1]. Whatever thresholds you choose, write them into your methods text.

One honest limitation to state once: similarity-based annotation describes what your sequence probably does, and any individual biological claim, particularly one heading toward an experiment or a publication, deserves targeted verification against the primary literature for that protein family.

Frequently Asked Questions

Can I annotate nucleotide BLAST hits the same way?

Yes, as long as you translate the identifiers correctly. Use blastx to compare translated nucleotide queries against protein databases, then the workflow is identical from Step 2 onward, because GO terms attach to protein function rather than to nucleotide sequence.

Do I need UniProt accessions, or can I use NCBI protein IDs?

You can use NCBI IDs, but you then need a cross-reference step to reach GO annotations, since UniProt's REST interface is keyed to UniProt accessions. The simplest reproducible path is to blast against UniProtKB directly so the accessions are native to your results.

What if my BLAST hit has no GO terms in UniProt?

Some entries, particularly in TrEMBL, have few or no GO annotations. Fall back to the InterPro or Gene3D cross-references on the same entry, keep the raw ID mapping in a separate column, and flag those rows as unannotated rather than leaving blanks.

Is it valid to transfer the hit's GO terms to my query sequence?

It is a standard and defensible inference when the hit passes identity and coverage thresholds, but it remains an inference. Report the underlying hit statistics alongside the transferred terms so readers can judge the strength of the transfer.

How do I handle thousands of queries without getting rate limited?

Batch your accessions into OR-joined queries against the UniProt search endpoint, request TSV or compressed output, keep the disk cache, and add a short sleep between chunks. For truly large jobs, download a full UniProt GO annotation file once and join locally instead of querying the API per entry.

Should I use blastp, blastx, or tblastn for this workflow?

Use blastp for protein queries against protein databases, blastx when your query is nucleotide, and tblastn when your database is nucleotide but your query is protein. The annotation steps do not change, only the identifier handling in Step 4 does.

Can I run GO enrichment testing on the annotated hits?

Yes, once every query carries GO terms, the annotated table is exactly the input that enrichment tools expect. Tools such as goatools in Python take a background set and a test set of GO identifiers and report overrepresented terms, which is the natural next analysis after this workflow.

Do I need an internet connection, or can I run everything locally?

You need connectivity for the UniProt fetches in Step 5 unless you pre-download UniProt annotation files, and everything else, including BLAST, OBO parsing, and table building, runs fully offline. The disk cache means a connected first run makes all later runs local.

Sources

1. PARSID: a Python script for automatic analysis of local BLAST results for a rapid molecular taxonomic identification
2. BLAST-QC: automated analysis of BLAST results
3. The PARIGA Server for Real Time Filtering and Analysis of Reciprocal BLAST Results
4. BlaSTorage: a fast package to parse, manage and store BLAST results
5. A Performance Analysis of Native JSON Parsers in Java, Python, MS.NET Core, JavaScript, and PHP
6. [[OA071] O-RAW: Ontology-guided radiomics analysis workflow](https://doi.org/10.1016/J.EJMP.2018.06.143)
7. Technical Note: Ontology-guided radiomics analysis workflow (O-RAW)