Lesson 2 of 13 · 11 min

Collecting and Curating Your Sequences

A phylogeny is only as good as the sequences that go into it. Use blastp when both your query and the database are protein, and use tblastn to hunt for homologs in nucleotide or genome databases, where the protein query is compared against subjects translated in all six reading frames. Start from one well-annotated query protein and search outward from there.

Searching With BLAST

bash

# query = one reviewed rbcL protein (RuBisCO large subunit)
# -remote runs on NCBI servers, so no local database is needed
blastp \
  -query rbcL_query.faa \
  -db nr \
  -remote \
  -evalue 1e-20 \
  -max_target_seqs 500 \
  -outfmt "6 saccver pident length qcovs evalue bitscore staxids sscinames" \
  -out hits.tsv

NP_051067.1	99.790	476	100	0.0	972	3702	Arabidopsis thaliana
YP_009179418.1	92.017	476	100	0.0	889	4113	Solanum tuberosum
ABK20720.1	88.235	475	99	0.0	845	161934	Beta vulgaris
AGT21077.1	85.014	474	98	0.0	812	4577	Zea mays
CAA25948.1	74.262	474	96	0.0	708	3055	Chlamydomonas reinhardtii

Not every hit belongs in your dataset. Keep matches with high query coverage and a small e-value, and discard short, partial alignments that will only inject gaps into the tree. Deduplicate the accessions before downloading so you never fetch the same protein twice.

bash

# keep hits with >=90% query coverage (column 4 = qcovs), unique accessions only
awk -F'\t' '$4 >= 90 {print $1}' hits.tsv | sort -u > keep.acc

# download those proteins from Entrez (NCBI EDirect: epost then efetch)
cat keep.acc | epost -db protein | efetch -format fasta > blast_hits.faa

Curated Sets From UniProt

bash

# reviewed (Swiss-Prot) rbcL proteins across green plants (Viridiplantae, taxid 33090)
# --data-urlencode lets curl encode the query for you
curl -s -G "https://rest.uniprot.org/uniprotkb/stream" \
  --data-urlencode "query=gene:rbcL AND reviewed:true AND taxonomy_id:33090" \
  --data "format=fasta" \
  -o uniprot_rbcL.faa

The ingroup is the set of taxa whose relationships you want to resolve, so sample it broadly and evenly rather than piling up sequences from one over-studied genus. The outgroup is one or a few sequences that sit just outside the ingroup, and because it roots the tree you need at least one. Aim for even taxon sampling across the ingroup paired with a single, well-placed outgroup.

Tip: Choose an outgroup that is close enough to align cleanly but clearly outside the ingroup. Too close and it may nest inside your tree; too distant and long-branch attraction plus alignment noise will distort the root.

bash

# combine both sources, then curate into a clean multi-FASTA
cat blast_hits.faa uniprot_rbcL.faa > all.faa

# remove exact-duplicate sequences (e.g. the same protein from both sources)
seqkit rmdup -s all.faa > nr.faa

# keep only near-full-length sequences (rbcL is ~475 aa); drop fragments
seqkit seq -m 440 -M 500 nr.faa > sized.faa

# rewrite headers into tree-friendly tips: Genus_species_Accession
# expr 1 = NCBI headers (... [Genus species]); expr 2 = UniProt headers (... OS=Genus species OX=...)
sed -E \
  -e 's/^>([^ .]+)[^ ]* .*\[([A-Za-z]+) ([A-Za-z]+).*/>\2_\3_\1/' \
  -e 's/^>[a-z]+\|([^|]+)\|[^ ]* .*OS=([A-Za-z]+) ([A-Za-z]+).*/>\2_\3_\1/' \
  sized.faa > clean.faa

seqkit stats clean.faa

Try it yourself: Pick a protein you know well (for example rbcL, cytochrome b, or elongation factor Tu). Run a remote blastp for it, add three to five reviewed UniProt sequences from your target clade, and choose one outgroup from a sister group. Then dedup and length-filter into a single multi-FASTA and confirm every header is a unique Genus_species_Accession with no spaces or punctuation before you move on to alignment.

Homology, Orthology, and Paralogy Multiple Sequence Alignment and Why It Matters