blast-align-tree

A pipeline to identify BLAST hits and perform phylogenetic analysis across multiple queries and local genome databases.

📖 Introduction

A common task in bioinformatics is to find similar genes across a set of genomes and compare them using phylogenetic methods. As an alternative to using online tools for such analyses, researchers may wish to download genomes of interest for local BLAST and downstream analyses. Homolog curation, header parsing, tree construction, and annotation of motifs, domains, and other datasets (e.g. gene expression) can give insights into a gene family of interest.

⚙️ Installation

Step 1. Clone the repo

The recommended install path starts with a clone — it gives you the conda environment YAMLs, helper scripts, and hmm_files/ locally so the commands below run as-is.

git clone https://github.com/steinbrennerlab/blast-align-tree.git
cd blast-align-tree

Step 2. Create the conda environment (non-Python tools)

blast-align-tree calls several external CLIs that pip can't install:

BLAST+ (makeblastdb, blastdbcmd, tblastn, blastp, psiblast)
Clustal Omega (clustalo) — default aligner
MAFFT (≥ 7) — alternative aligner
trimAl — alignment cleanup
FastTree and/or RAxML-NG — tree inference
R with ggtree, ape, phytools, ggplot2, optparse, treeio, tidytree, broom
HMMER (hmmscan, hmmpress) — optional, only needed for --hmm

Environment YAML files for each platform live under environments/. Create and activate the env with conda, mamba, or micromamba (mamba / micromamba are noticeably faster):

Linux:

mamba env create -f environments/bat-environment-linux.yml
mamba activate bat

macOS (Apple Silicon or Intel):

mamba env create -f environments/bat-environment-ARMorIntel-mac.yml
mamba activate bat

Windows:

mamba env create -f environments/bat-environment-windows.yml
mamba activate bat

Platform-specific gaps

Only the Linux YAML is a complete one-shot install. The macOS and Windows YAMLs each cover a different half of the stack; fill in the missing pieces after creating the env.

macOS — R / plotting stack is not in the YAML. The mac YAML ships the CLI tools (BLAST+, MAFFT, Clustal Omega, trimAl, FastTree, RAxML-NG, HMMER) plus Python, but not R or any of the tree-plotting packages — bioconda's R / Bioconductor coverage is patchy on Apple Silicon (osx-arm64), so R is installed via CRAN instead. Install R separately (e.g. brew install r or from CRAN) and then run the bundled installer to pull ggtree, ggplot2, ape, phytools, tidytree, treeio, optparse, broom, and Biostrings:

Rscript environments/install_r_deps.R

Windows — external CLI tools and one R package are not in the YAML. Bioconda doesn't build BLAST+, MAFFT, Clustal Omega, trimAl, FastTree, RAxML-NG, or HMMER for Windows, so the Windows YAML only provisions the majority of the R stack plus Python. Install the R package 'Biostrings' separately or run the bundled installer. Then install the CLI tools from their vendors (add each to PATH after install):

NCBI BLAST+ (ncbi-blast-*+-x64-win64.exe)
MAFFT (all-in-one Windows build)
Clustal Omega (Windows build) — or compile from source
FastTree (FastTree.exe)
trimAl (Windows build) — or compile from source
HMMER — use win-hmmer

Unfortunately Windows users cannot use RAxML-NG for tree construction. Use WSL2 instead.

If you'd rather avoid hand-installing these, run the Linux YAML under WSL2 (Ubuntu) and run the pipeline from there — everything is installable through conda/mamba environment via bioconda on WSL2.

3. Install the Python package

From inside the activated env from Step 2:

pip install blast-align-tree

This installs three console commands:

Command	Purpose
`blast-align-tree`	Run the pipeline (BLAST → align → tree → visualize)
`bat-genome-selector`	Tkinter GUI for building `blast-align-tree` commands
`blast-align-tree-fetch`	Download bundled genome FASTAs into `./genomes/`

If you're developing against the repo, use an editable install instead:

pip install -e .

Verify everything is wired up:

blast-align-tree --check-env

Skipping the clone

If you don't want to clone the repo, you can grab just the environment YAML for your platform directly from GitHub and run mamba env create -f <file> against it:

Linux: bat-environment-linux.yml
macOS (Apple Silicon or Intel): bat-environment-ARMorIntel-mac.yml
Windows: bat-environment-windows.yml

Then pip install blast-align-tree inside the env. Note that some tutorial commands (hmmpress hmm_files/kinase.hmm, the scripts/populate_manifest.py helper, etc.) assume a repo clone.

4. Build HMM profiles with `hmmpress`

The repo ships .hmm input files under hmm_files/ (e.g. hmm_files/kinase.hmm) but not the .h3* binary indices — those are build artifacts and are deleted from the repo. Before using --hmm, build the index locally:

hmmpress hmm_files/kinase.hmm

This produces kinase.hmm.h3f, .h3i, .h3m, .h3p alongside the input.

Working directory

blast-align-tree is filesystem-driven: it reads genome FASTAs from ./genomes/ and writes run outputs into the current working directory. If you cloned the repo (the recommended install), run everything from the repo root — ./genomes/ and ./datasets/ are already in place.

If you installed via pip without cloning, pick any project folder instead:

mkdir ~/bat-project && cd ~/bat-project

You can keep several projects side-by-side; each has its own ./genomes/ and its own run outputs. Pipeline commands invoked from any other directory won't find your genomes.

Fetching bundled genome databases

A clone ships the plant default set already in ./genomes/ (🌱🌿🫘🫛🍅🍃) plus 📊 Klepikova in ./datasets/. Animal and fungal genomes are hosted on GitHub Releases instead of being checked in — they're too large to bundle in the clone and exceed PyPI's per-file limit, so they're fetched on demand:

blast-align-tree-fetch                       # default plant set 🌱🌿🫘🫛🍅🍃📊
blast-align-tree-fetch --all                 # everything in the manifest 🌍
blast-align-tree-fetch --list                # show available genomes + sizes
blast-align-tree-fetch human_cds mouse_cds   # fetch specific genomes by name

Pip-only installs (no clone) use the same command to pull the plant set into the project directory.

blast-align-tree ships a manifest (blast_align_tree/data/genomes_manifest.json inside the installed package) listing each hosted genome with its URL and sha256 checksum. The manifest is read-only from a user's perspective — you don't edit it directly.

Files land in ./genomes/. Animal and fungal genomes sort into ./genomes/animals/. Downloads are checksum-verified against the manifest and decompressed automatically, so re-running is safe — already present files are skipped if their hash matches. After each successful download blast-align-tree-fetch runs makeblastdb on the FASTA (the nucleotide / protein mode is auto-detected), so the files are ready for the pipeline with no extra step. Pass --no-index to skip that.

The default set 🌱🌿🫘🫛🍅🍃📊 is:

TAIR10 CDS — Arabidopsis thaliana coding sequences
TAIR10 proteins — Arabidopsis thaliana proteome
Pvul218 CDS — Phaseolus vulgaris (common bean) coding sequences, Phytozome genome ID 218, v1.0
Vung469 CDS — Vigna unguiculata (cowpea) coding sequences, Phytozome genome ID 469, v1.1
NbLab360 CDS — Nicotiana benthamiana coding sequences (LAB360 v103)
NbLab360 proteins — Nicotiana benthamiana proteome (LAB360 v103)
📊 Klepikova atlas subset — Arabidopsis expression overlay dataset (lands in ./datasets/, not ./genomes/)

--all 🌍 additionally pulls the rest of the hosted lineup (fetched into ./genomes/animals/): human CDS, mouse CDS, rat CDS, chimp CDS, zebrafish CDS, fruit fly CDS, C. elegans CDS, yeast ORFs (S. cerevisiae S288C).

Run blast-align-tree-fetch --list for the full lineup and sizes.

Using your own genome files instead of (or alongside) the hosted set? Drop any FASTA into ./genomes/ — no manifest edit, no reinstall needed. See Adding a new genome below.

🖥️ GUI: `bat-genome-selector`

The easiest way to build a valid blast-align-tree invocation is the Tkinter GUI. Launch it from a directory that contains a genomes/ folder:

bat-genome-selector

Key features

Working directory. Shows the current project path beneath the logo; click it to open the folder in Explorer, Finder, or your Linux file manager.
Auto-discovery. Scans ./genomes/ (recursively) for .fa, .faa, .fas, .fasta, .fna files and ignores BLAST index sidecars.
Header auto-detection. Peeks at the first FASTA record in each database and suggests a plausible -hdr token (e.g. gene:, locus=, polypeptide=), with a live "Parsed name" preview so you can see exactly what will appear on the tree.
Per-row controls. One row per genome: include checkbox, query column, -hdr, -hdr_sfx, -n (hits to keep), nucleotide/protein type, and a "build DB" shortcut that runs makeblastdb when the BLAST indices are missing.
Bulk actions. Select All Hits, Deselect All Hits, Clear Fields, Refresh, plus a Default -n spinbox and Set All -n button.
Options panel. Aligner (Clustal Omega or MAFFT + mode), tree builder (FastTree or RAxML), BLAST type (tblastn/blastp), thread count.
Advanced panel (collapsible): outgroups (-add, -add_db) and the tip to root the tree on (-a), AA slice (-aa, single range applied to all queries, or one range per query — see the tutorial section below), motif patterns (regex or PROSITE, overlap toggle), and HMM profiles (--hmm).
Generate Command / Copy to Clipboard. Produces a ready-to-paste blast-align-tree … command.
Recent Runs tab. Lists past ENTRY/runs/TIMESTAMP/ directories in the current working directory so you can quickly jump back to prior results.

📚 Full Tutorial

The tutorial commands below use the default files included in the repo clone. If you start from a new directory, blast-align-tree-fetch will download the genomes and sample dataset, but not environment and hmm files. Make sure to set up all packages, as described in the Installation section Steps 1-4.

Generate a simple tree

The example below runs the pipeline for a SERK query and redraws the resulting tree with a new subnode/outgroup. It searches three genomes from the default fetched set: 🌱 TAIR10cds, 🫘 Pvul218cds, and 🫛 Vung469cds.

blast-align-tree -q AT4G33430.1 -qdbs TAIR10cds.fa \
                 -n 15 15 15 \
                 -dbs TAIR10cds.fa Pvul218cds.fa Vung469cds.fa \
                 -hdr gene: polypeptide= locus=

The pipeline creates a folder AT4G33430.1/ in your working directory. The timestamped run root keeps the tree PDFs. Newick tree files, gene lists, alignment FASTAs, mappings, features, BLAST hit FASTAs, and per-genome summaries go under AT4G33430.1/runs/<TIMESTAMP>/genes_alignments_trees/.

After the run finishes, the pipeline prints a re-draw hint. For example, to reroot on an outgroup (-a AT5G10290) and zoom in on a subnode (-n 45). (Note that you must replace the path to your R script and timestamp with details from your specific run.)

Rscript "<bundled-visualize_tree.r>" -e AT4G33430.1 -b SERK_tree \
        --subdir "runs/<TIMESTAMP>" -a AT5G10290 -n 45

The bundled path is printed for you at the end of each pipeline run, already wrapped in double quotes — keep the quotes when copy-pasting, especially on Windows, where unquoted paths can cause Rscript to segfault if they contain spaces or backslashes that the shell misparses.

Run blast-align-tree for a different gene, ACC Oxidase

Find 15 homologs of Arabidopsis ACC Oxidase 1 from three plant genomes using tblastn against complete CDS databases. -q specifies the query locus, -qdbs the database it lives in, -dbs the databases to search, and -hdr the regex tokens used to parse gene names out of each database's FASTA headers.

blast-align-tree -q AT2G19590.1 -qdbs TAIR10cds.fa \
                 -n 15 15 15 \
                 -dbs TAIR10cds.fa Pvul218cds.fa Vung469cds.fa \
                 -hdr gene: polypeptide= locus=

This creates AT2G19590.1/ with tree PDFs at the timestamped run root and Newick tree files, alignment files, BLAST hit FASTAs, and per-genome summaries under genes_alignments_trees/.

A powerful feature of ggtree is the ability to plot associated data. Each run produces two complementary tree PDFs:

Text version — gene symbols and dataset values printed as labels next to each tip.
Heatmap version — the same tree with associated data rendered as a coloured heatmap alongside the tips.

By default, both include expression data from the Klepikova Arabidopsis expression atlas (headers are matched to the AtGenExpress / eFP browser tissue naming). The overlay TSV is fetched automatically into ./datasets/ by blast-align-tree-fetch. The screenshot below shows the heatmap version:

A separate PDF with .MSA.pdf appended shows a cartoon alignment — useful for spotting large differences in domain architecture. Open the underlying FASTA files in genes_alignments_trees/ to inspect the alignment in detail.

All three PDF versions include a branch-length scale bar labeled in substitutions/site by default, sized to the displayed tree (including subtrees). Trees without branch lengths or with a zero branch-length span omit the scale bar.

Redraw the ACC Oxidase tree

You can re-run visualize_tree.r at any time to produce new PDFs. The pipeline prints a ready-to-edit Rscript … command at the end of each run; copy it and tweak options such as:

-b <NAME> — filename stem for the new PDFs
-a <ID> — reroot on this outgroup (also available on the pipeline itself as blast-align-tree -a <ID>, see below)
-n <NODE> — draw a subtree at this node (use --help for the full option list)
-k 1 — show bootstraps
-l 0 — hide node number labels
-m 2 — enlarge gene-symbol text

For example, reroot the default ACO tree on JRG21 (AT2G38240) and zoom into the ACO clade at node 58:

Rscript "<bundled-visualize_tree.r>" -e AT2G19590.1 -b ACO_v3 \
        --subdir "runs/<TIMESTAMP>" -a AT2G38240 -n 58 -k 1 -l 0 -m 2

The -n option is especially helpful for extracting a subset of the tree as a FASTA. Sequences are listed in the FASTA in the same order as the tree, and trimAl is used to strip blank-only alignment columns — useful for a quick view of conserved residues (e.g. the ACO active site) in a viewer like AliView.

BLASTP instead of TBLASTN

Use --blast_type blastp against protein databases. The example below pulls 10 NIMIN-1 homologs from the Arabidopsis and Nicotiana benthamiana proteomes:

blast-align-tree --blast_type blastp \
                 -q AT1G02450.1 -qdbs TAIR10protein.fa \
                 -n 10 10 \
                 -dbs TAIR10protein.fa NbLab360.v103.gff3.CDS.fasta.AA.fasta \
                 -hdr gene: id

Multiple queries

You can pass several query sequences with -q; the pipeline extracts each from the database listed at the matching position in -qdbs, de-duplicates, then searches each -dbs entry. The example below uses three queries drawn from two databases and searches two other databases:

blast-align-tree -q AT5G45250.1 Phvul.007G077500.1 AT5G17890.1 \
                 -qdbs TAIR10cds.fa Pvul218cds.fa TAIR10cds.fa \
                 -n 3 4 \
                 -dbs TAIR10cds.fa Vung469cds.fa \
                 -hdr gene: locus=

Identifier collisions and de-duplication

-hdr turns a FASTA description into the short identifier used as a tree tip label, and that mapping is not one-to-one. -hdr gene: against TAIR10cds collapses AT2G27490.1 and AT2G27490.4 onto AT2G27490; 10,690 of the 48,321 records in that database share a locus identifier with another record. Two genomes can also independently use the same gene symbol.

The pipeline never resolves these silently. After the BLAST searches finish, but before results from several queries are merged, every identifier claimed by more than one source record is reported and resolved by an explicit rule:

Situation	What happens
Several isoforms or duplicated loci in one database share an identifier	The longest amino-acid sequence is retained; ties break on the lexicographically smallest source ID. The others are dropped and logged with the ID that replaced them.
Two queries hit different records that share an identifier	Same rule, no separate prompt.
Two databases use the same identifier	Both records are kept, each tagged with its genome (`RPS5_TAIR10cds`, `RPS5_NbLab360`). Nothing is dropped, so a genome can never vanish from the tree.
The same record is hit by several queries	Not a collision. Counted as overlap and reported separately at the end of the run.

Note that the retained isoform is the longest one, which is not necessarily the isoform you passed as -q.

With a terminal attached, each query's collisions are shown for confirmation before anything is merged:

  [ids] query 'AT1G71830.1': 6 identifier(s) claimed by more than one hit
    AT2G13800  <- keep AT2G13800.1 (601 aa)
        drop AT2G13800.3 (601 aa)
        drop AT2G13800.2 (484 aa)
        Deduplicate these 6 identifier(s) for 'AT1G71830.1'? [Y/n] (10s -> Y)

Answering n keeps every record instead, disambiguated as AT2G13800__1, AT2G13800__2, … Those suffixed labels will not join to --datasets tables keyed on the bare identifier, which is the trade-off for retaining all isoforms.

Each prompt waits 10 seconds and then takes its default, so a long run left unattended finishes rather than stalling at the question. Prompts that timed out are named in the end-of-run summary and in the # unanswered prompts (default taken): line of deduplication_log.tsv, so a default that was taken is never mistaken for one that was agreed to. Raise or lower the wait with CONFIRM_TIMEOUT_SECONDS in blast_align_tree/cli.py; use --duplicates auto to skip the prompts altogether.

Use --duplicates to control this:

Value	Behaviour
`ask` (default)	Confirm each query's collisions, waiting 10 s per prompt before taking the default. Falls back to `auto` when there is no terminal, so batch and HPC runs never block.
`auto`	Apply the rules without prompting and warn at the end.
`fail`	Stop the run and list the collisions, so you can pick a more specific `-hdr` / `-hdr_sfx`.

Every run writes deduplication_log.tsv next to the tree PDFs in the timestamped run folder — one row per affected record, with the stage, query, database, identifier, action, source ID, amino-acid length, the original FASTA header, and the reason:

stage         query        database      identifier  action   source_id     aa_len
within_query  AT1G71830.1  TAIR10cds.fa  AT2G13800   kept     AT2G13800.1   601
within_query  AT1G71830.1  TAIR10cds.fa  AT2G13800   dropped  AT2G13800.3   601

The end of the run summarises the collisions and, separately, how much the queries overlapped:

  Identifier collisions
    within a query          6  (isoforms / duplicated loci sharing one identifier)
    between queries         0  (same identifier, different source records)
    between databases       0  (genome tag appended, nothing dropped)
    records dropped         9  (longest amino-acid sequence retained)

  Overlap between queries
    AT1G71830.1 ∩ AT4G33430.1: 12 shared of 15/15 hits in TAIR10cds.fa

Internal stop codons and reading frame

tblastn hits are retrieved as whole nucleotide records and translated in forward frame +1 of the record as stored in the database. The frame and strand are never inferred from the data. This is correct for CDS databases (all bundled genomes) and wrong for transcript databases whose entries carry a 5′ UTR.

--internal-stops controls what happens when a record contains an in-frame stop:

Value	Behaviour
`truncate` (default)	Translation ends at the first in-frame stop. The reported protein is the truncated product the locus encodes. Residue p maps to nucleotides `3p−2..3p` of the record.
`readthrough`	Every codon kept in register; each in-frame stop written as `X`. Use when downstream domain content matters, e.g. comparing domain architecture across a degraded locus.
`excise`	Legacy v1.0 behaviour: stop codons deleted and the flanking sequence joined. Produces a protein the genome does not encode and shifts downstream coordinates by one residue per excised stop. Kept only to reproduce older runs.

A well-formed CDS — begins with ATG, ends with a stop, no internal stops — gives identical output under all three. The nucleotide outputs (all_hits.nt.fa, <entry>.nt.parse.merged.fa) are always the untransformed records as returned by blastdbcmd, whatever the policy.

Every tblastn run writes translation_report.tsv beside the tree PDFs, one row per sequence:

Column	Meaning
`identifier`, `source_id`, `database`, `queries`	which record this is, and which queries found it
`action`	`kept`, or `dropped` if identifier reconciliation discarded it
`nt_len`, `n_codons`	length of the retrieved record
`frame`, `strand`	always `+1` and the stored strand, stated per row so the assumption is never implicit
`start_codon`, `terminal_stop`	whether the record looks like a complete CDS
`n_internal_stops`, `first_stop_aa_pos`	how many stops, and where the first one is (1-based)
`aa_len_reported`	length actually written, under the policy in force
`aa_len_ranking`	length excluding stops; identical under every policy, and what isoform ranking uses — so the stop policy cannot change which isoform wins a collision
`aa_after_first_stop`	codons remaining in frame after the first internal stop
`longest_downstream_orf_aa`	longest stop-free run after it — what matters for judging whether a real ORF survives downstream
`best_frame`	the frame that reads furthest, from a six-frame diagnostic translation
`flags`	see below
`policy`	the `--internal-stops` value in force

Flags:

Flag	Meaning
`internal_stop`	at least one in-frame stop before the end of the record
`no_start_codon` / `no_terminal_stop`	the record does not look like a complete CDS
`len_not_multiple_of_3`	trailing bases could not be translated
`frame_mismatch:<frame>`	another frame reads substantially further
`possible_utr`	that frame is a forward frame — the database may hold transcripts with 5′ UTRs rather than CDS
`possible_wrong_strand`	that frame is a reverse frame
`orf_mostly_downstream`	more coding potential after the first stop than before it; the truncated protein is unlikely to be the real product
`empty_translation`	no complete codon could be translated

Flagged records are reported, never re-framed — the frame actually used is always +1. A database that trips possible_utr on many hits probably holds transcripts rather than CDS.

  [translation] 10 sequences translated in frame +1 (strand: + (record as stored in the database; not inferred))
    Internal stop policy: truncate - translation ends at the first in-frame stop; the truncated protein is reported
    7 of 10 sequences contain an internal stop -> truncated at the first stop
    5 flagged orf_mostly_downstream: more coding potential after the first stop than before it

    !! 4 sequences read further in another frame.
       Frame +1 is used regardless - these were NOT re-framed.
       2 flagged possible_utr: the database may hold transcripts with 5' UTRs rather than CDS.
       2 flagged possible_wrong_strand.
           4  SynthTranscripts.fa
       Check these databases before interpreting the affected proteins.

Rooting the tree on an outgroup with `-a`

Outgroup handling is two steps: pull the sequence into the run with -add/-add_db, then root the tree on it with -a/--reroot. -a takes a single tip and is passed straight through to visualize_tree.r, so the first-pass PDFs come out already rooted — no redraw needed:

blast-align-tree -q AT2G19590.1 -qdbs TAIR10cds.fa                  -n 15 15 -dbs TAIR10cds.fa Vung469cds.fa                  -hdr gene: locus=                  -add AT2G38240 -add_db TAIR10cds.fa                  -a AT2G38240

The ID must match the tip label as it appears in the tree, which is whatever -hdr parsing leaves behind — gene: strips the isoform suffix, so the tip is AT2G38240, not AT2G38240.1. If the ID is not a tip, the run says so, suggests the near match, and draws the tree unrooted rather than failing at the last step; you can then reroot with the Rscript … redraw command printed at the end of the run.

Slicing query amino-acid ranges with `-aa`

-aa trims each query to a sub-range (0-based, Python-style: start-inclusive, end-exclusive) before BLAST. The hits are not trimmed, only the query sequences.

Two forms are supported:

Single range, applied to every query — two bare integers:
```
-aa 10 200
```
One range per query — START:END tokens, one per entry in -q (matching order). Use - to skip slicing for a particular query:
```
blast-align-tree \
    -q   Phvul.007G077500.1 Phvul.002G196200.1 Phvul.004G100000.1 Phvul.010G073300.1 \
    -qdbs Pvul218cds.fa     Pvul218cds.fa      Pvul218cds.fa      Pvul218cds.fa \
    -n   15 15 \
    -dbs TAIR10cds.fa Vung469cds.fa \
    -hdr gene: locus= \
    -aa  705:885 701:1164 903:1104 861:1086
```
This is handy when you want the same homologous sub-region from each query (for example a kinase domain whose ungapped coordinates differ between sequences). Mix and match skips with -, e.g. -aa 705:885 - - 861:1086.

Adding a new genome

The manifest used by blast-align-tree-fetch covers only the hosted set shipped with the package. Using your own FASTA does not require editing the manifest — just drop the file into ./genomes/ and build a BLAST database. Both the pipeline and bat-genome-selector auto-discover any .fa, .faa, .fas, .fasta, or .fna file in ./genomes/ and its subfolders (so ./genomes/mygroup/foo.fa works the same as ./genomes/foo.fa).

For each new genome you need a local BLAST database.

For CDS files:

makeblastdb -in GenomeCDS.fa -parse_seqids -dbtype nucl

For protein files:

makeblastdb -in GenomeProteins.fa -parse_seqids -dbtype prot

Example: add Nicotiana tabacum and rebuild an earlier tree

Download an annotated tobacco proteome into ./genomes/:

N. tabacum v4.5 from Sol Genomics → Nitab-v4.5_proteins_Edwards2017.fasta

Build the BLAST database:

cd genomes
makeblastdb -in Nitab-v4.5_proteins_Edwards2017.fasta -parse_seqids -dbtype prot
cd ..

Inspecting the first record of Nitab-v4.5_proteins_Edwards2017.fasta shows that headers use a gene id followed by a description — -hdr id keeps just the first token.

Now build a SOBIR1 homolog tree across Arabidopsis, N. benthamiana, and the freshly added tobacco proteome:

blast-align-tree --blast_type blastp \
                 -q AT2G31880.1 -qdbs TAIR10protein.fa \
                 -n 10 10 10 \
                 -dbs TAIR10protein.fa \
                       NbLab360.v103.gff3.CDS.fasta.AA.fasta \
                       Nitab-v4.5_proteins_Edwards2017.fasta \
                 -hdr gene: id id \
				 --hmm kinase.hmm

The run produces a tree PDF with tobacco SOBIR1 homologs slotted in alongside the N. benthamiana and Arabidopsis sequences.

Rebuild the SOBIR1 tree with MAFFT and RAxML

By default the pipeline aligns with Clustal Omega and infers the tree with FastTree, but both are swappable. The command below rebuilds the same SOBIR1 tree (Arabidopsis + N. benthamiana + tobacco) using MAFFT in linsi mode and RAxML-NG for the tree inference:

blast-align-tree --blast_type blastp \
                 --aligner mafft --mafft_mode linsi \
                 --tree_builder RAxML \
                 -q AT2G31880.1 -qdbs TAIR10protein.fa \
                 -n 10 10 10 \
                 -dbs TAIR10protein.fa \
                       NbLab360.v103.gff3.CDS.fasta.AA.fasta \
                       Nitab-v4.5_proteins_Edwards2017.fasta \
                 -hdr gene: id id

RAxML-NG is noticeably slower than FastTree but provides maximum- likelihood branch support via bootstrapping. Comparing the two trees is a quick sanity check that any clades you care about are stable across inference methods.

Use repeated rounds of querying to refine your trees, search different genome versions, and compare aligners / tree builders in order to draw strong conclusions about your gene family of interest.

Name		Name	Last commit message	Last commit date
blast_align_tree		blast_align_tree
datasets		datasets
environments		environments
genomes		genomes
hmm_files		hmm_files
images		images
scripts		scripts
tests		tests
.gitattributes		.gitattributes
.gitignore		.gitignore
LICENSE		LICENSE
README.md		README.md
RELEASE_NOTES_v1.0.0.md		RELEASE_NOTES_v1.0.0.md
RELEASE_NOTES_v1.0.1.md		RELEASE_NOTES_v1.0.1.md
RELEASE_NOTES_v1.0.4.md		RELEASE_NOTES_v1.0.4.md
bat-flowchart.pptx		bat-flowchart.pptx
pyproject.toml		pyproject.toml

Navigation Menu

Repository files navigation

blast-align-tree

📖 Introduction

⚙️ Installation

Step 1. Clone the repo

Step 2. Create the conda environment (non-Python tools)

Platform-specific gaps

3. Install the Python package

Skipping the clone

4. Build HMM profiles with `hmmpress`

Working directory

Fetching bundled genome databases

🖥️ GUI: `bat-genome-selector`

Key features

📚 Full Tutorial

Generate a simple tree

Run blast-align-tree for a different gene, ACC Oxidase

Redraw the ACC Oxidase tree

BLASTP instead of TBLASTN

Multiple queries

Identifier collisions and de-duplication

Internal stop codons and reading frame

Rooting the tree on an outgroup with `-a`

Slicing query amino-acid ranges with `-aa`

Adding a new genome

Example: add Nicotiana tabacum and rebuild an earlier tree

Rebuild the SOBIR1 tree with MAFFT and RAxML

About

Releases

Packages

Contributors

Languages

Latest commit

History

Folders and files

Repository files navigation

blast-align-tree

📖 Introduction

⚙️ Installation

Step 1. Clone the repo

Step 2. Create the conda environment (non-Python tools)

Platform-specific gaps

3. Install the Python package

Skipping the clone

4. Build HMM profiles with hmmpress

Working directory

Fetching bundled genome databases

🖥️ GUI: bat-genome-selector

Key features

📚 Full Tutorial

Generate a simple tree

Run blast-align-tree for a different gene, ACC Oxidase

Redraw the ACC Oxidase tree

BLASTP instead of TBLASTN

Multiple queries

Identifier collisions and de-duplication

Internal stop codons and reading frame

Rooting the tree on an outgroup with -a

Slicing query amino-acid ranges with -aa

Adding a new genome

Example: add Nicotiana tabacum and rebuild an earlier tree

Rebuild the SOBIR1 tree with MAFFT and RAxML

About

Resources

Stars

Watchers

Forks

Releases

Packages

Contributors

Languages

4. Build HMM profiles with `hmmpress`

🖥️ GUI: `bat-genome-selector`

Rooting the tree on an outgroup with `-a`

Slicing query amino-acid ranges with `-aa`