NCBI Blast

Taking a set of unknown sequences and working out what genes they are

How to read this module

This is a walkthrough, not a script — the code is here to read and adapt, and the page does not execute anything. To run it you need NCBI BLAST+ and somewhere with a bit of memory; see Computing Resources.

The idea

You have a few thousand sequences and no idea what they are — off a de novo assembly of a non-model organism, say, with no annotation to inherit.

The strategy is comparison: search each unknown sequence against a database of proteins whose function is known, and take the annotation from the best match.

The database here is UniProt/Swiss-Prot, the manually reviewed part of UniProt. Smaller than the full database and far better curated — the right trade for annotation.

The three steps

flowchart LR
  A["Swiss-Prot<br/>FASTA"] -->|makeblastdb| B["BLAST<br/>database"]
  C["Your<br/>sequences"] --> D
  B --> D["blastx"]
  D --> E["Hit table<br/>outfmt 6"]
  E -->|join in R| F["Gene names,<br/>GO terms"]

There is a script that does all of this

sr320/workflow-annotation wraps the entire workflow in one command:

bash blast2slim.sh -i transcripts.fasta

It downloads and indexes Swiss-Prot, picks blastx or blastp for you, pulls UniProt annotations over the REST API, maps GO terms to GO-Slim categories with GOATOOLS, and writes TSVs ready for downstream use. --diamond swaps in DIAMOND when BLAST+ is too slow for the dataset.

Work through the steps below first so you know what it is doing — then use the script for real work rather than retyping this by hand.

Everything below is those three steps. Set your BLAST+ location once if it is not already on your PATH:

export BLAST=/home/shared/ncbi-blast-2.11.0+/bin   # adjust to your machine
$BLAST/blastx -version                             # check what you actually have

Step 1 — Build the database

Download Swiss-Prot and index it. Record the release — annotation results shift between releases, so it belongs in your methods.

cd ../data
RELEASE=r2025_03    # check https://www.uniprot.org/help/release-statistics

curl -O https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
mv uniprot_sprot.fasta.gz uniprot_sprot_${RELEASE}.fasta.gz
gunzip -k uniprot_sprot_${RELEASE}.fasta.gz

mkdir -p ../blastdb
makeblastdb \
-in ../data/uniprot_sprot_${RELEASE}.fasta \
-dbtype prot \
-out ../blastdb/uniprot_sprot_${RELEASE}

Using a variable for the release is deliberate: type it by hand in both places and they will eventually disagree, leaving you with a database built from a different file than you think — and nothing will error.

-out gives a prefix, not a filename. Several index files are written from it, and that same prefix is what you pass to -db.

Step 2 — Run BLAST

Pick the right program

Program	Query	Database
`blastn`	nucleotide	nucleotide
`blastp`	protein	protein
`blastx`	nucleotide, translated	protein
`tblastn`	protein	nucleotide, translated

Transcripts against a protein database means blastx.

curl -k https://eagle.fish.washington.edu/cnidarian/Ab_4denovo_CLC6_a.fa \
> ../data/Ab_4denovo_CLC6_a.fa

grep -c ">" ../data/Ab_4denovo_CLC6_a.fa    # how many sequences?

blastx \
-query ../data/Ab_4denovo_CLC6_a.fa \
-db ../blastdb/uniprot_sprot_${RELEASE} \
-out ../output/Ab_4-uniprot_blastx.tab \
-evalue 1E-20 \
-num_threads 20 \
-max_target_seqs 1 \
-outfmt 6

-evalue 1E-20 — how many hits this good you would expect by chance. Smaller is stricter.
-num_threads 20 — set to what the machine has (nproc).
-max_target_seqs 1 — best hit only.
-outfmt 6 — tabular; everything downstream expects this.

Tip

Run head -100 of your input through the whole thing first. You will find path typos in 30 seconds rather than 3 hours.

Read the output

head -2 ../output/Ab_4-uniprot_blastx.tab
wc -l ../output/Ab_4-uniprot_blastx.tab

There is no header row, so you have to know the 12 columns:

#	1	2	3	4	5	6	7–8	9–10	11	12
	`qseqid`	`sseqid`	`pident`	`length`	`mismatch`	`gapopen`	`qstart` `qend`	`sstart` `send`	`evalue`	`bitscore`

Expect fewer hits than sequences. Anything without a match above threshold produces no line — normal for a non-model assembly, not a failure.

Step 3 — Add annotation

BLAST told you which protein each sequence matches, not what it does. Join to a UniProt table for that.

Swiss-Prot IDs arrive as sp|P12345|NAME_SPECIES; split on the pipe to expose the accession:

tr '|' '\t' < ../output/Ab_4-uniprot_blastx.tab \
> ../output/Ab_4-uniprot_blastx_sep.tab

That turns column 2 into three, shifting everything after it right by two — accession lands in column 3, E-value moves from 11 to 13.

library(tidyverse)

bltabl <- read.csv("../output/Ab_4-uniprot_blastx_sep.tab", sep = "\t", header = FALSE)
spgo   <- read.csv("https://gannet.fish.washington.edu/seashell/snaps/uniprot_table_r2023_01.tab",
                   sep = "\t", header = TRUE)

annot_tab <- left_join(bltabl, spgo, by = c("V3" = "Entry")) |>
  select(V1, V3, V13,
         Protein.names, Organism,
         Gene.Ontology..biological.process., Gene.Ontology.IDs)

write.table(annot_tab, "../output/blast_annot_go.tab",
            sep = "\t", row.names = FALSE, quote = FALSE)

left_join keeps every BLAST hit even where UniProt has no matching row, so misses show as NA rather than disappearing — a silently shrinking table is the easiest way to lose data without noticing.

That is the workflow. annot_tab has a gene name, organism, and GO terms for every sequence that matched.

Note

Most hits will be to human, mouse, and other heavily studied organisms. That is a property of Swiss-Prot, not a discovery about your animal.

Variation — starting from a reference proteome

More often you are in a species that already has a reference, and the goal is to attach GO terms to an existing gene set for enrichment testing later. Same three steps, two differences.

Work backwards from your identifiers. Whatever you annotate must produce IDs matching your count matrix, or the final join silently fails. For a Crassostrea virginica matrix built against the RefSeq GFF, genes look like gene-LOC111099033, and the file carrying those in its headers is translated_cds.faa.

Use the current assembly

Older notebooks use GCF_002022765.2 (C_virginica-3.0), which is now suppressed in RefSeq. The current reference is GCF_053477285.1 (ASM5347728v1). See Bivalves.

cd ../data
ASM=GCF_053477285.1_ASM5347728v1
curl -O https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/053/477/285/${ASM}/${ASM}_translated_cds.faa.gz
gunzip -k ${ASM}_translated_cds.faa.gz

These are already proteins, so use blastp:

blastp \
-query ../data/${ASM}_translated_cds.faa \
-db ../blastdb/uniprot_sprot_${RELEASE} \
-out ../output/Cvir_transcds-uniprot_blastp.tab \
-evalue 1E-20 -num_threads 40 -max_target_seqs 1 -outfmt 6

The gene ID lives in the FASTA header ([gene=LOC111126949]), not in the BLAST output, so pull it out and pair it with the Swiss-Prot accession:

g.spid <- left_join(blast, cdsftab, by = "V1") |>
  mutate(gene = str_extract(V2.y, "(?<=\\[gene=)\\w+"),          # from the header
         SPID = str_extract(V2.x, "(?<=\\|)[^\\|]*(?=\\|)")) |>  # from sp|P12345|NAME
  distinct(gene, SPID)

write.table(g.spid["SPID"], "../output/SPID.txt",
            sep = "\t", row.names = FALSE, quote = FALSE)

Then map those accessions to GO terms at https://www.uniprot.org via Retrieve/ID mapping, adding Gene Ontology IDs under Customize columns, and join the result back on SPID. The output — each LOC gene with its GO terms — is what an enrichment test expects.

When it doesn’t work

Empty output — the -db prefix does not match what makeblastdb -out wrote, or the path is wrong for where you are running.
Join gives all NA — keys do not match. head() both key columns; usually an identifier still has a sp| prefix or version suffix attached.
Row count grows after a join — duplicate keys on the right-hand table.
Running for hours — check input size and -num_threads, then go back to the head -100 test.

Next steps

With gene-to-GO mappings you can run functional enrichment on a gene set — the differentially expressed genes from an experiment, for instance — and ask which biological processes are over-represented.

Now that you know what each step does, use workflow-annotation rather than assembling this by hand each time. It also produces GO-Slim terms, which the manual workflow above does not.

For the animal behind the example, see Bivalves.