# How to Use NCBI CD-Search with Custom Databases: A Guide By the end of this guide you will be able to do three things. First, you will understand how the NCBI Conserved Domain Search (CD-Search) works and what it can and cannot do on the web. Second, you will build your own conserved domain database from protein alignments you trust, using PSSMs (position-specific score matrices) and the `makeprofiledb` utility from the NCBI BLAST+ suite. Third, you will search that custom database with `rpsblast+` for protein queries and `rpstblastn` for nucleotide queries, and you will know how to verify that the results are correct. To follow along you need a Unix-like environment (Linux or macOS, or Windows with WSL), the BLAST+ command-line tools version 2.26 or newer, a multiple-sequence aligner such as MAFFT or MUSCLE, a few protein sequences in FASTA format, and roughly 1 GB of free disk space. Everything in this guide uses free, official NCBI tools and standard open-source software. ## What CD-Search and the CDD Actually Do A conserved domain is a stretch of [protein sequence](/blog/guides/protein-sequence) whose three-dimensional structure and function have been preserved across evolution. A kinase domain, a helix-turn-helix DNA-binding motif, and an aminoacyl-tRNA synthetase catalytic core are all conserved domains. Detecting them is often the fastest way to go from an unknown sequence to a defensible functional hypothesis. CD-Search is the NCBI service that compares a query protein against the Conserved Domain Database (CDD). CDD is a curated NCBI databases resource that contains two kinds of domain models: models written and curated by NCBI staff, and models imported from partner resources such as Pfam, TIGRFAMs, COG, and SMART. Each model is stored as a [multiple sequence alignment](/blog/guides/multiple-sequence-alignment) and a precomputed PSSM. A PSSM, or position-specific score matrix, is a table of log-odds scores with one column per position in the alignment and one row per amino acid. Conserved catalytic positions get large positive scores, tolerant positions get scores near zero, and positions that reject certain residues get negative scores. A PSSM captures far more family-specific information than a single consensus sequence, which is why profile methods find remote homologs that simple pairwise search misses. CD-Search runs the query against the database of PSSMs using RPS-BLAST, which stands for reverse position-specific BLAST. RPS-BLAST is the profile-aware counterpart of BLASTP. The subject is a profile rather than a single sequence. The web version of CD-Search at the NCBI structure site performs exactly this comparison and returns domain footprints, E-values, bit scores, and links to the parent domain models on the CDD website. The core fact that shapes this whole guide: the public CD-Search web service only searches the official CDD. If your family of interest is a novel lineage-specific protein, a viral gene of your own annotation, or a model that CDD does not yet cover, you need a custom database. Building one is a local command-line job, and the rest of this article walks through it step by step. ## Where CD-Search Fits Among NCBI Databases If you are still learning how to use NCBI resources, the CDD occupies a specific slot. The [Protein database](/blog/guides/protein-database) stores sequences. BLAST compares sequences to sequences. CDD stores alignments of homologous domains and their profiles. CD-Search is the query interface to CDD. When CD-Search hits a domain, the domain model page links onward to structure records in MMDB, publications, Gene records, and superfamily trees. In practice, a CD-Search hit is often the highest-value single annotation you can attach to a new protein, because the domain model carries curated notes on active-site residues and functional sites. For large query sets, NCBI provides the Batch CD-Search service, an HTTP API that accepts a file of proteins and returns the same results in a tabular form. This service, like the web form, searches the official CDD. It is still worth knowing because it is the natural benchmark for validating a custom database, as described later in this guide. ## Web Service or Local Custom Database: Which One You Need ![Decision path asking whether the family is missing from CDD, whether paralog models or offline runs are needed, ending at local build or web service.](/article-images/illustrations/cd-search-custom-database-web-or-local.webp) Use the web service when your goal is annotating sequences against well-characterized domains. It requires no installation and the models are maintained by NCBI. Build and search a custom database when any of the following applies: 1. Your protein family is absent from or poorly represented in CDD. 2. You need domain models tuned to your own data, for example a paralog-specific model that separates two nearly identical domains. 3. You need reproducible, versioned, offline annotation for a pipeline, and you cannot depend on a remote service. 4. You want to control the model composition, score scale, and reporting format. You cannot upload your own PSSMs to the public CD-Search service. This is a deliberate design boundary, and the official NCBI path for custom models is the local toolchain covered below. ## Prerequisites and Installation Install BLAST+ with conda, which also avoids permission problems with system package managers: ```bash conda create -n rps -c conda-forge -c bioconda blast mafft muscle conda activate rps ``` Confirm the tools exist: ```bash rpsblast+ -version makeprofiledb -version rpstblastn -version mafft --version ``` You should see version numbers for the BLAST+ programs. Next, fetch the CDD source utility package from the NCBI FTP site, which includes helper scripts such as `msa2pssm.pl` for building PSSMs directly from alignments: ```bash wget https://ftp.ncbi.nlm.nih.gov/pub/mmdb/cdd/cdd_src.tar.gz tar -xzf cdd_src.tar.gz ls cdd_src/ ``` If you prefer not to use helper scripts, you can generate PSSMs with PSI-BLAST alone, and this guide shows that route as the main path because it depends only on BLAST+. ## Building a Custom CD-Search Database Step by Step ![Five-step flow from choosing seed sequences through alignment, PSSM conversion, makeprofiledb compilation, and seed round-trip testing.](/article-images/illustrations/cd-search-custom-database-build-pipeline.webp) The workflow has five stages: pick seed sequences, align them, convert the alignment to a PSSM, compile PSSMs into an RPS database, and test the database. The worked example builds a tiny two-model database, but the commands scale to hundreds of models. ### Step 1: Choose Seed Sequences and Download Them Select representative proteins for your family. Aim for sequences that span the diversity you care about and exclude near-identical duplicates, because redundancy flattens the profile. For a small curated family, five to twenty well-chosen seeds are enough to start. Pull them from [NCBI Protein](/blog/guides/ncbi-protein) with EFetch, which is part of the E-utilities API: ```bash curl "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=protein&id=WP_003042283.1,WP_010927252.1&rettype=fasta&retmode=text" > seeds.fasta ``` Name each FASTA header deliberately, for example: ``` >HTH_arp1|model_v1|arabinose_responsive_HTH MKKVTIGK... ``` The identifier of the seed sequence becomes the identifier of the domain model in search output. A clean header like `HTH_arp1` will make downstream parsing trivial. Headers with spaces or unusual punctuation cause parsing pain later, so keep them to letters, digits, underscores, and pipes. ### Step 2: Build a Multiple Sequence Alignment Align the seeds with MAFFT: ```bash mafft --auto seeds.fasta > hth.aln ``` Check the alignment in a viewer such as AliView or Jalview before proceeding. The region you intend to model should be aligned without obvious frame shifts. Trim columns outside the domain if the seeds include long unrelated tails: ```bash # Example: keep columns 15 to 104 of the alignment using seqkit seqkit subseq --range 15:104 hth.aln > hth_core.aln ``` Profile quality depends directly on alignment quality in the modeled region. Misaligned columns produce distorted PSSM scores and false confidence around positions that should be variable. ### Step 3: Convert the Alignment to a PSSM The most reliable route is PSI-BLAST with the alignment supplied as an input MSA. PSI-BLAST reads the alignment, builds the initial profile, optionally refines it against a reference database, and writes the PSSM: ```bash psiblast \ -in_msa hth_core.aln \ -db swissprot \ -num_iterations 2 \ -out_pssm hth.pssm \ -out hth_psiblast.log ``` Parameters that matter here: - `-in_msa` supplies your alignment instead of a single query. The first sequence in the file acts as the master, so put your designated representative first. - `-db` is the database used for profile refinement. A curated protein set such as Swiss-Prot gives controlled refinement. If you want a profile derived only from your alignment, set `-num_iterations 1`. - `-num_iterations 2` performs one refinement pass after building the initial profile. - `-out_pssm` writes the PSSM in ASN.1 ASCII format, which is the input format `makeprofiledb` expects. If you have many alignments, loop over them: ```bash for aln in models/*.aln do base=$(basename "$aln" .aln) psiblast -in_msa "$aln" -db swissprot -num_iterations 2 \ -out_pssm "models/${base}.pssm" -out "logs/${base}.log" done cat models/*.pssm > all_models.pssm ``` Concatenating the individual PSSM files into one input file is how `makeprofiledb` receives multiple models in a single call. ### Step 4: Create the Database with makeprofiledb Compile the PSSMs into an RPS database: ```bash makeprofiledb \ -in all_models.pssm \ -out my_domains \ -title "My domain models v1.0" \ -dbtype rps \ -threshold 9.82 \ -scale 100.0 ``` You should see progress output ending without error, and a set of database files will appear: ``` my_domains.rps my_domains.pog my_domains.pdb my_domains.phr my_domains.pin my_domains.psq ``` Two parameters deserve explanation. `-dbtype rps` tells the builder to treat the input as reverse position-specific profiles, which is mandatory for RPS-BLAST searches. `-scale 100.0` sets the bit-score scaling to the convention used by CDD itself, so your bit scores are comparable in magnitude to those reported by the public CD-Search service. Omitting `-scale` is one of the most common mistakes with custom databases and it makes your scores inconsistent with everything else. The threshold value shown follows the CDD convention. If the command fails with a message in the BLAST+ bracketed style, such as an error naming the input file, the usual cause is an empty or malformed PSSM file. Check that each PSSM file is non-empty and that the PSI-BLAST step in Step 3 completed without reporting zero iterations. ### Step 5: Test the Database Against Its Own Seeds A database that cannot recover its own training sequences is broken. Search the seeds back: ```bash rpsblast+ \ -query seeds.fasta \ -db my_domains \ -evalue 1e-2 \ -outfmt "6 qseqid qstart qend sseqid evalue bitscore" \ -out roundtrip.tsv ``` Expected output in `roundtrip.tsv` looks like this: ``` HTH_arp1 3 96 HTH_arp1 2.1e-28 126 seed2 3 96 HTH_arp1 4.7e-25 119 seed3 5 98 HTH_arp1 9.3e-22 111 cupin_A 22 138 cupin_A 6.8e-31 138 ``` Each seed should hit its own model with a strong score, typically many orders of magnitude below the E-value cutoff. If a seed fails to recover its own model, the alignment or PSSM for that model is defective, and you should fix that model before annotating anything else. ## Searching Your Custom Database ### Protein Queries with rpsblast+ With the database verified, search real query proteins: ```bash rpsblast+ \ -query predicted_proteins.faa \ -db my_domains \ -evalue 1e-3 \ -outfmt "6 qseqid qstart qend sseqid evalue bitscore score length" \ -out domain_hits.tsv ``` Output rows mirror the fields you requested: ``` protein_0421 41 212 cupin_A 3.9e-19 94.5 94.5 172 protein_0421 1 38 HTH_arp1 2.6e-04 31.0 31.0 38 protein_1107 88 340 HTH_arp1 1.1e-27 121.0 121.0 253 ``` The `qstart` and `qend` columns give the domain footprint on the query, which is what you need to annotate multi-domain proteins. A protein with hits from two different models at different coordinates, like `protein_0421` above, has a two-domain architecture. Useful parameter notes: - `-evalue` controls the hit cutoff. Start at 1e-2 for exploration, then tighten to 1e-3 or 1e-5 for reporting. Unlike BLASTP, RPS-BLAST profiles carry a lot of information per position, so genuine family members often score far better than 1e-5. - `-outfmt 5` produces XML if you need structured output for parsing libraries. - `-seg no` disables low-complexity filtering if you are specifically hunting repetitive domains, but leave the default on for general annotation. - Add `-num_threads 4` (or however many cores you can spare) on large query sets. ### Nucleotide Queries with rpstblastn To find domains in raw contigs without predicting genes first, use `rpstblastn`, which translates the nucleotide query in all six frames: ```bash rpstblastn \ -query assembly_contigs.fna \ -db my_domains \ -evalue 1e-3 \ -seg yes \ -outfmt "6 qseqid qstart qend sframe sseqid evalue bitscore" \ -out contig_domain_hits.tsv ``` The `-seg yes` flag turns on compositional adjustment for low-complexity translated sequence, which reduces spurious hits in repeat-rich regions. The `sframe` column tells you which reading frame carried the hit, useful when you reconstruct gene models manually. ### Reading and Interpreting the Output The most important fields are the pair of E-value and bit score, the domain footprint, and the model identifier. Two interpretation rules apply. First, E-values in RPS-BLAST are computed against the size of your specific database. Adding fifty models to your database changes E-values for every search against it. When you report E-values from a custom database, state the database version, or your numbers will not be reproducible by anyone else. Second, related models can produce overlapping hits on the same query region, the classic behavior of any profile library. When two models of the same superfamily both hit a region, keep the hit with the better bit score and record the alternative as a secondary annotation rather than discarding it silently. ## Using the Batch Web CD-Search API Even though the public service cannot search custom models, the Batch CD-Search API is part of the toolset because it serves two roles in a custom-database project. It annotates large query sets against the official CDD quickly, and it provides an independent benchmark against which you can sanity-check your own models. Submit a query set with `curl`: ```bash curl -F queries=@batch_proteins.faa \ -F dbs=cdd \ -F evalue=0.01 \ "https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi" > job_response.txt ``` The immediate response reports the job state and an identifier: ``` #Status:Pending SearchId=QID202611110045 ``` Poll the same endpoint with the identifier until the status reads Success: ```bash curl "https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi?cdsid=QID202611110045" ``` A successful response returns tabular hits with columns for the query, the query start and end, E-value, bit score, the domain accession, the short model name, an incomplete flag, and the superfamily. Respect the service by submitting jobs sequentially and batching your queries into single files rather than thousands of tiny requests. As with all NCBI services, heavy automated use should follow the guidance posted on the NCBI website. ## Key Options Reference | Command | Parameter | What it does | Typical value | |--|--|--|--| | `psiblast` | `-in_msa` | Uses your alignment as the profile source | aligned FASTA file | | `psiblast` | `-num_iterations` | Number of profile refinement rounds | 1 or 2 | | `psiblast` | `-out_pssm` | Writes PSSM in ASN.1 ASCII for `makeprofiledb` | one file per model | | `makeprofiledb` | `-dbtype` | Declares database type | `rps` | | `makeprofiledb` | `-scale` | Bit-score scale | `100.0` (CDD convention) | | `makeprofiledb` | `-threshold` | Model score threshold at build time | `9.82` (CDD convention) | | `rpsblast+` | `-evalue` | Hit significance cutoff | `1e-2` to `1e-5` | | `rpsblast+` | `-outfmt` | Output format | `6` with chosen fields | | `rpstblastn` | `-seg` | Low-complexity filtering for translated queries | `yes` | ## Verifying That Your Results Are Right A custom database is only as trustworthy as its validation. Run three checks before you use it in production annotation. **Round-trip recovery.** As shown in Step 5, every seed must hit its own model strongly. This catches corrupted PSSMs, truncated alignments, and database build failures. **Cross-check against public CD-Search.** Take ten query proteins that your custom database annotates, run them through the public CD-Search web interface or the Batch API, and compare domain footprints where both services have coverage. Your custom model should roughly agree with the closest official model on domain boundaries. Systematic disagreement in footprint boundaries usually means your alignment contained flanking sequence you did not intend to model. **Independent-family search.** Search a curated set of proteins known to be unrelated to your family, such as a random sample of Swiss-Prot sequences. A healthy custom database returns no hits, or only weak borderline hits, against unrelated proteins. Clusters of false positives indicate a PSSM with distorted scores, often caused by a small or highly gapped alignment. Record the model construction date, the reference database used for PSI-BLAST refinement, and the `makeprofiledb` parameters in a README that ships with the database. Anyone who reruns your annotation later, including you, will need those details. ## Common Mistakes and Limitations These are the errors that actually appear in practice, with their causes and fixes. **Expecting the web page to accept custom models.** The CD-Search web form and Batch API search only the official CDD. If you need your own models, the local toolchain in this guide is the supported route. **Building a PSSM from a single sequence.** A profile built from one sequence carries no position-specific information beyond trivial derivations, and it will not outperform a simple BLAST search. Models need an alignment of genuinely related sequences. If only one seed exists, gather homologs first with PSI-BLAST and align them. **Feeding a sloppy alignment into the PSSM.** Frameshifted columns, untrimmed fusion proteins, and fragments aligned end to end all corrupt the profile. Inspect the alignment visually before Step 3. **Skipping `-scale 100.0`.** Without the CDD scaling convention, your bit scores occupy a different range than every CDD-derived score you will ever compare them to. Set it at build time. **Using `rpsblast+` on nucleotide queries.** The database contains protein profiles, so a nucleotide query produces a molecule-type mismatch error and no results. Use `rpstblastn`, which translates the query in six frames. **Leaving generic FASTA headers on seeds.** If every seed is named `seq1`, `seq2`, and so on, every hit in your output will carry an uninformative model identifier. Rename seeds before building. **Treating E-values as transferable between databases.** E-value depends on database size and composition. An E-value of 1e-5 against a 40-model custom database is not equivalent to 1e-5 against the full CDD. Report the database version alongside any E-value. **Overlapping hits reported as independent domains.** Related models hit the same region. Establish a rule, for example one model per query region based on the best bit score, and apply it consistently in post-processing. **Never validating.** The single most damaging mistake is putting an unvalidated custom database into an annotation pipeline. The three checks in the previous section take under an hour for a small database and catch most construction errors. **A structural limitation to keep in mind.** A custom PSSM is only as good as the alignment and sequence sampling behind it. A model built from one narrow clade will miss divergent relatives, and a model built from too-diverse sequences will lose family specificity. When your family grows or your sampling changes, rebuild the models and version the database rather than editing in place. ## How a PSSM Actually Encodes a Family, and How to Read One Most guides treat the PSSM as a black box that `psiblast` fills and `makeprofiledb` consumes. Opening the box changes how you build models, because nearly every construction mistake in practice traces back to a misunderstanding of what the matrix records. A PSSM produced by PSI-BLAST has one column per position of your alignment and, in the ASCII form that feeds `makeprofiledb`, a set of scores per position covering the twenty standard amino acids. Each score is a log-odds value: it compares the probability of observing a given residue at that position in your family against the probability of observing it by chance in background protein sequence. Three qualitative regimes matter when you inspect a model: - Strongly conserved catalytic or structural positions produce large positive scores for one or two residues and large negative scores for everything else. These columns do most of the discrimination work. - Positions under weak constraint produce small scores spread across several residues. These columns contribute little to the total score and are where profile drift, if it happened during refinement, shows up first. - Positions that tolerate two distinct residue classes, for example a small hydrophobic versus a small polar residue, produce a bimodal pattern. This is normal and is one of the genuine advantages of profiles over consensus sequences, which must pick a single letter per column. A schematic example makes the idea concrete. The following table is illustrative, not a real output dump, but it shows the shape you should expect in a healthy model: | Position | A | C | D | E | F | ... | Interpretation | |-|-|-|-|-|-|-|-| | 1 | 3 | -4 | -3 | -2 | -5 | ... | Strongly prefers Ala | | 2 | -1 | 4 | -2 | -3 | -2 | ... | Strongly prefers Cys | | 3 | 0 | -1 | 1 | 0 | -1 | ... | Weakly constrained | | 4 | -6 | -7 | 5 | 3 | -8 | ... | Acidic residue required | When you inspect a real PSSM from your own build, look for three warning signs. First, a model in which nearly every column is weakly constrained usually means the alignment was too diverse or contained non-members, and the profile will retrieve loosely related proteins along with your family. Second, a model in which nearly every column is strongly positive for one residue suggests redundancy in the seed set: near-identical sequences reinforce each other and the profile becomes an overfit fingerprint of a narrow clade. Third, columns that appear to reject almost everything, with uniformly negative scores across all residues, often indicate alignment artifacts such as a frameshifted region, where no real residue pattern exists to model. Two structural details of the format are worth knowing because they affect decisions elsewhere in the workflow. The first sequence in the input alignment acts as the master: gaps are inserted relative to it, and the coordinate system of the final model follows it. If your master sequence carries a large insertion relative to the rest of the family, you get a model padded with gap columns that contribute nothing. Choose a representative, complete seed as the first sequence in every alignment file. Second, the profile records what the alignment plus any refinement iterations produced, not what you intended. If a refinement round pulled in a misleading relative, that relative is now baked into the matrix, and no post-hoc adjustment fixes it. Rebuild instead. ## Building the Seed Set: The Selection Strategy That Determines Model Quality The existing walkthrough tells you what to do once you have seeds. This section covers the harder upstream question: which sequences to put in the alignment, and which to deliberately leave out. Start from a discovery search, not from memory. Take one or two confident representatives of your family and search them against a broad database with PSI-BLAST or plain BLASTP, collecting hits at a significance threshold you would be comfortable defending in a methods section. The goal at this stage is a candidate pool, not a final seed set. Then curate the pool against four criteria: 1. **Coverage of diversity.** The seeds should span the range of organisms, or the range of sequence divergence, within which you want the model to detect members. If your family has clear clades and you want one model per clade, plan separate alignments per clade rather than one giant alignment. 2. **Completeness.** Drop fragments, partial genes at contig edges, and sequences with obvious sequencing artifacts. A fragment aligned against full-length seeds introduces long gap columns that distort neighboring positions. 3. **Domain boundaries.** Where possible, extract only the domain region for the alignment rather than aligning full-length proteins whose flanking domains will misalign against each other. Profiles model what you give them; unmodeled flanking sequence is noise. 4. **Redundancy control.** Near-identical sequences add nothing and skew residue frequencies. Collapse the pool at an identity cutoff you choose and document. For a family-level model, many practitioners cluster in the high-80s to high-90s percent identity range and keep one representative per cluster; for a subfamily-specific model, keep more of the within-group diversity. The exact number matters less than applying it consistently and recording it. A practical curation pass that uses only BLAST+ itself: ```bash # Make a searchable database from the candidate pool makeblastdb -in candidates.faa -dbtype prot -out candb # Find near-duplicates of each candidate blastp -query candidates.faa -db canddb \ -outfmt "6 qseqid sseqid pident length" \ -out self_hits.tsv ``` Sort the output by identity and inspect pairs with very high identity over most of the domain. Remove one member of each near-identical pair, preferring the entry with the cleaner annotation, a complete-length sequence, and a stable accession. Two selection errors deserve special warnings because they are common and subtle. Including a single distant outgroup in an otherwise coherent seed set flattens the entire profile: the shared signal between the close members gets diluted across columns that now have to accommodate the outgroup, and family specificity drops. Conversely, sampling only from one well-studied lineage makes a model that is sharp within that lineage and effectively blind elsewhere. The honest fix for both problems is the same: decide explicitly whether you are building a broad family model or a narrow subfamily model, sample accordingly, and, when in doubt, build both variants and let the validation step in a later section decide. ## Scaling Up: Organizing a Multi-Model Database Build The worked example in the main walkthrough builds two models by hand. Real projects grow to dozens or hundreds of models, and ad-hoc command entry stops working around model number five. The fix is a fixed directory layout and a build script, which together turn the pipeline into something a colleague, or you in six months, can rerun exactly. A layout that works well: ``` domaindb/ v1.0/ seeds/ one FASTA per model: hth_arp1.seeds.faa, cupin_A.seeds.faa aln/ one MAFFT output per model, trimmed pssm/ one PSSM per model logs/ psiblast and mafft logs, named per model db/ makeprofiledb output files scripts/ README.md ``` The one-file-per-model convention matters more than any particular layout. When each model has its own seed file, alignment, PSSM, and log, a failing model is isolated: you can fix and rebuild it without touching anything else, and the build script can skip models that already have valid outputs. A build loop with failure checks looks like this: ```bash #!/usr/bin/env bash set -euo pipefail for seeds in seeds/*.seeds.faa; do model=$(basename "$seeds" .seeds.faa) mafft -auto "$seeds" > "aln/${model}.aln" 2> "logs/${model}.mafft.log" psiblast -in_msa "aln/${model}.aln" -db swissprot \ -num_iterations 2 \ -out_pssm "pssm/${model}.pssm" \ -out "logs/${model}.psiblast.log" # Fail fast if the PSSM came out empty or truncated if [ ! -s "pssm/${model}.pssm" ]; then echo "EMPTY PSSM: ${model}" >&2 exit 1 fi done cat pssm/*.pssm > db/all_models.pssm makeprofiledb -in db/all_models.pssm \ -out db/my_domains_v1 \ -title "My domain models v1.0" \ -dbtype rps \ -threshold 9.82 \ -scale 100.0 ``` Three properties of this script earn their keep in practice. The `set -euo pipefail` line makes the whole run stop on the first failure, so you never compile a database that silently dropped models. The size check on each PSSM catches the most common silent failure, a PSI-BLAST run that produced no usable output, at the moment it happens rather than at search time. And every artifact carries its model name, so the trail from a suspicious hit back to the alignment that produced it is one directory listing long. When the model count passes a few dozen, consider generating the build script from a manifest table instead of editing the script directly. A simple CSV with columns for model name, seed accessions, iteration count, and any trimming coordinates becomes the single source of truth, and the script becomes a dumb executor. This also gives you, almost for free, the metadata table that belongs in your README. ## Alignments for Profiles: Refinement, Trimming, and Quality Control Profile construction is unusually unforgiving of alignment errors because every column becomes a scoring position. A misaligned column does not merely look wrong in a viewer; it becomes a position in the model that rewards a residue pattern no real protein satisfies. Inspect every alignment before generating its PSSM. In a viewer such as AliView or Jalview, check for these specific defects: - **Staggered blocks**, where one sequence appears shifted by one or two residues relative to the rest across a stretch. This is the signature of a frameshifted sequence or of a divergent member the aligner could not place. Remove or realign that sequence. - **Columns that are gaps in all but one or two sequences.** These usually indicate insertions in outlier sequences. Ask whether the insertion is biologically part of the family or an idiosyncrasy of one seed; if the latter, remove the seed or trim the insertion. - **Overhanging ends**, where some seeds extend well beyond the shared core. Decide explicitly whether the model should cover those regions. Usually the right move is to trim everything to the shared core, because unaligned overhangs produce noisy terminal columns. - **A member that aligns sensibly to one half of the core and badly to the other half.** This often means the "family" is really two related families, or the seed is a fusion or a domain swap. Split the model or drop the seed. For trimming, aligner-independent coordinates are more reproducible than hand-drawn boundaries. Anchor the model region to positions in the master sequence, record those coordinates in your manifest, and apply them programmatically, as the existing guide shows with `seqkit subseq`. If the master sequence changes in a later version of the model, re-derive the coordinates rather than copying them forward. One decision worth making deliberately is the choice of aligner mode. For a family of clearly homologous, moderately divergent proteins, a single-pass global mode such as the MAFFT default auto setting is appropriate and fast. For a family with long insertions in some members or with only distant mutual similarity, MAFFT's iterative refinement modes, which are documented in the MAFFT documentation, spend more compute to place divergent sequences more sensibly. Whatever you choose, record the exact MAFFT version and mode in the model's log directory; alignment output is not reproducible across versions, and your PSSM inherits that non-reproducibility. Finally, resist the temptation to hand-edit the alignment and skip re-verification. Manual adjustments are sometimes correct and often valuable, but they must be saved, versioned, and noted, because an untracked hand edit is indistinguishable from an accident six months later. ## PSI-BLAST Refinement Choices: When to Iterate and When Not To The walkthrough uses `-num_iterations 2` as the default. This section explains the decision behind that number, because it is the single highest-leverage choice in profile construction. With `-in_msa`, PSI-BLAST does two conceptually distinct things. On the first iteration, it converts your alignment into a profile and, if a database is supplied, searches that database. On each subsequent iteration, it incorporates the sequences it found above its inclusion threshold into the profile and searches again. The `-out_pssm` file always reflects the profile as of the last completed iteration. The tradeoff is straightforward. Refinement rounds broaden the profile with real homologs you did not have in your seed set, which typically improves sensitivity. But each round also risks profile drift: if the search database contains a sequence that is similar to your family without being a member, for example a promiscuous domain with shared catalytic chemistry, it gets incorporated and degrades specificity. Drift is worse on large, poorly curated databases and milder on curated sets, which is why the walkthrough recommends Swiss-Prot as the refinement database rather than `nr`. The parameter that governs what gets incorporated is the inclusion E-value threshold, exposed by PSI-BLAST as `-inclusion_ethresh`. A stricter value lets fewer database sequences into the profile per round; a looser value admits more, raising both the sensitivity gain and the drift risk. If you refine, keep the value strict and inspect the log to see exactly which sequences were incorporated at each round. The log lists them, and reading that list is the difference between controlled refinement and uncontrolled drift. Three build recipes cover almost all real needs: - **Pure alignment profile.** `-num_iterations 1`, database still required syntactically but effectively unused. The PSSM derives only from your seeds. Most specific, least sensitive. Use when you trust your seed sampling and prioritize precision, for example when discriminating between two paralogs. - **One refinement round on a curated database.** `-num_iterations 2` against Swiss-Prot. The walkthrough default. A good balance for most families. - **Multiple refinement rounds.** Reserve for families where you know your seed sampling is thin and the surrounding database is clean. Check the incorporation list at every round and stop the moment a questionable sequence appears. The recommendation that follows from this: build the pure-alignment version and the refined version of each important model, run both through your validation checks, and keep the variant that wins on your own hold-out data. It costs little and replaces a matter of taste with a measurement. ## Hold-Out Validation: Estimating Sensitivity and Specificity Before You Trust the Database Round-trip recovery, the check in the main walkthrough, verifies construction. It does not estimate how the model will behave on sequences it has never seen, and it is systematically optimistic: a seed that helped build the profile will almost always score against it. Hold-out validation closes that gap. The procedure is leave-one-out testing at the model level: 1. For each model, remove one seed from the seed set, ideally one that is moderately divergent from the others. 2. Rebuild that model's alignment and PSSM from the reduced seed set. 3. Search the held-out seed, and a small panel of other held-out family members if you have them, against the rebuilt model. 4. Record whether the held-out sequence is recovered, at what E-value, and whether its footprint matches the known domain boundaries. A shell sketch of the rebuild step: ```bash # Given held-out sequence ID held_out_1, filter it from the seed set grep -v -A1 "held_out_1" seeds/cupin_A.seeds.faa \ | grep -v '^-$' > aln/cupin_A_loo.faa mafft -auto aln/cupin_A_loo.faa > aln/cupin_A_loo.aln psiblast -in_msa aln/cupin_A_loo.aln -db swissprot \ -num_iterations 2 -out_pssm pssm/cupin_A_loo.pssm \ -out logs/cupin_A_loo.log ``` Then search the held-out sequence against a small database containing only the leave-one-out model. Repeat across models. A model that recovers all of its own seeds but fails its own hold-out is overfit: it has memorized the seeds rather than learned the family. The remedy is usually more diverse seeds or, if that is impossible, a stricter refinement configuration. Complement sensitivity testing with a specificity panel. Assemble fifty to a hundred proteins you are confident are unrelated to any model in your database, search them, and count hits. A handful of borderline hits at the edge of your working cutoff is normal; clusters of comfortable hits on unrelated proteins mean at least one model is diffuse, and you can identify the culprit model directly from the `sseqid` column of the output. Record the results as a small table per model version: number of seeds, hold-out recovery yes or no, lowest hold-out E-value, and count of unrelated-set hits. This table, plus the README, is what turns a pile of PSSM files into a defensible, versioned annotation resource. ## Troubleshooting the Toolchain: Symptoms, Causes, and Fixes The errors that come up with this toolchain cluster into a small number of patterns. The table below covers the ones seen most often. Exact error wording varies between BLAST+ versions, so match on the symptom and the stage rather than the string. | Stage | Symptom | Likely cause | Fix | |-|-|-|-| | `psiblast -in_msa` | Fails immediately, message about the input alignment | Alignment is not truly aligned: rows of unequal length, or the file is plain FASTA with no gaps | Verify all rows have equal length; check the file actually came from the aligner, not a raw FASTA | | `psiblast` | Log reports zero or very few sequences incorporated | Refinement database too small or misnamed, or the family is genuinely rare in it | Confirm the database exists and is versioned; accept low incorporation if it is real | | `psiblast` | PSSM file empty or missing | Run errored before the write, or `-out_pssm` path is wrong | Read the full log; check directory permissions | | `makeprofiledb` | Fails naming the input file | Empty or malformed concatenated PSSM file; a truncated PSSM from a failed PSI-BLAST run | Check each PSSM is non-empty; rerun failed models before concatenating | | `makeprofiledb` | Fails with a parameter complaint | Missing or misspelled flag, or `-dbtype` not set to `rps` | Compare the command against the reference table in the main guide | | `rpsblast+` | Error about molecule type | Nucleotide query against a protein-profile database | Use `rpstblastn` for nucleotide queries | | `rpsblast+` | Database open failure | Wrong path, or the database was built with an incompatible BLAST+ version | Confirm all database files are present in the directory; rebuild the database with the same BLAST+ version used for searching | | `rpsblast+` | Valid run, but no hits on proteins you know are family members | E-value cutoff too strict, or the model misses the divergent clade | Loosen to 1e-2 temporarily; if the member still does not appear, the model needs broader seeds | | `rpstblastn` | Many short spurious hits in repeat-rich contigs | Low-complexity translated sequence, short chance ORFs | Keep `-seg yes`; post-filter by alignment length and E-value together | | Any search | Hits show absurdly large E-values despite good alignments | Expected with a very small custom database; E-values scale with database size | Report bit scores and the database version; do not compare these E-values to CDD values | Two of these deserve elaboration. The BLAST+ version mismatch is the most confusing for newcomers: a database built under one BLAST+ release may refuse to open under another, with an error that looks like file corruption when it is really a format incompatibility. The reliable policy is to pin your BLAST+ version in your environment specification and rebuild the database whenever you upgrade the toolchain, then rerun the round-trip check. The no-hits-on-known-members case is the one worth taking seriously rather than working around. Loosening the E-value until the member appears is a diagnostic step, not a fix. If a true family member is invisible even at 1e-2, the profile genuinely lacks the signal, and the correct response is to revisit the seed set and the alignment for that model. ## Post-Processing Hits: Collapsing Overlaps and Producing Annotation Files Raw RPS-BLAST output is a hit list, not an annotation. Two transformations turn it into something usable, and both are simple enough to script reliably. The first transformation is overlap collapsing. Related models hit the same query region, and a naive listing reports several "domains" where there is one. A consistent rule, best hit per query region, applied per query protein: ```python import csv from collections import defaultdict rows = list(csv.DictReader(open("domain_hits.tsv"), delimiter="\t")) # Group hits per query, sort by start coordinate per_query = defaultdict(list) for r in rows: per_query[r["qseqid"]].append(r) kept = [] for q, hits in per_query.items(): hits.sort(key=lambda r: int(r["qstart"])) for h in hits: h_start, h_end = int(h["qstart"]), int(h["qend"]) # Does this hit overlap a already-kept, stronger hit? conflict = any(h_start <= int(k["qend"]) and h_end >= int(k["qstart"]) for k in kept if k["qseqid"] == q) if not conflict: kept.append(h) # Hits overlapping a stronger hit are dropped here; extend this # to record them as secondary annotations if you need them. ``` In production, sort candidate keepers by bit score first so that the strongest hit in a region claims it, and log every discarded overlapping hit to a side file rather than discarding silently. The side file costs nothing and rescues you when a discarded "overlap" turns out to be a genuinely interesting second domain that your interval logic misjudged. The second transformation is emission into a standard annotation format. A GFF3 record for each accepted hit needs the query as the sequence identifier, the footprint as start and end, the model as the feature type or attribute, and the E-value and bit score as attributes: ``` protein_0421 customCDD domain 41 212 . + . model=cupin_A;evalue=3.9e-19;bitscore=94.5 ``` Note that the strand column is only meaningful for `rpstblastn` output, where you derive it from the sign of the frame field. For `rpsblast+` protein queries, fill it with a placeholder and record that choice in your pipeline documentation. When you also run the official Batch CD-Search against the same query set, as recommended for benchmarking, you end up with two annotation layers. Merge them with an explicit precedence rule rather than ad-hoc judgment: for overlapping footprints, prefer the official curated model when one exists, because it carries curated active-site and functional-site annotations; prefer your custom model where official coverage is absent or where your model was built specifically to discriminate subfamilies that official models lump together. Whichever rule you choose, write it down and apply it in code. ## PSSM Profiles and Profile HMMs: How This Workflow Compares to HMMER The most common question from people setting up custom domain detection is whether to use this PSSM-based workflow or to build profile hidden Markov models with HMMER. The honest answer is that they overlap in purpose and differ in mechanics, and the right choice depends on your family and your constraints. The structural difference is how each model treats insertions and deletions. A PSSM scores aligned columns independently and has no explicit model of insert and delete events; gap handling in the underlying search is position-independent in character. A profile HMM carries explicit match, insertion, and deletion states per position, which makes it inherently better at modeling families where lengths vary substantially between members. As a rule of thumb drawn from broad community experience: for short, compact domains with moderate divergence, a well-built PSSM performs competently; for long domains with variable loop regions or for very diverse superfamilies, profile HMMs generally have the edge. The practical differences favoring the PSSM route are real nonetheless: - **CDD compatibility.** If you want scores and conventions comparable to the NCBI CD-Search ecosystem, including the bit-score scale set by `-scale 100.0`, the PSSM route is the native one. - **Translated search convenience.** `rpstblastn` gives you six-frame translated search against your profiles with one command, which is operationally convenient for contig-level screening. The HMMER equivalent workflow exists but requires more plumbing. - **Toolchain simplicity.** Everything in this guide comes from one BLAST+ installation plus an aligner. There is value in a dependency tree that small. The pragmatic production pattern used by many groups is to build both representations from the same curated alignment and require agreement for high-confidence calls. Where the PSSM and the HMM both hit a region, annotate with confidence; where only one hits, flag for review. This costs one extra model build per family and buys you a genuine cross-check, which is more valuable than any single method's tuning. ## Where This Workflow Shines: Practical Applications Custom domain databases earn their keep in several recurring situations, and knowing them helps you decide how much effort to invest in model building. **Annotating understudied organisms and viruses.** When you sequence an organism from a poorly sampled lineage, official domain libraries cover its proteins thinly. A custom database built from the lineage's own protein families, even a modest one covering your twenty most important families, turns a sea of hypothetical proteins into structured annotation. Because the models come from the lineage's own diversity, they detect members that broad, generalist models miss. **Subfamily discrimination.** Official superfamily-level models deliberately lump related subfamilies, which is correct for general-purpose annotation and useless when the question is which paralog you are looking at. A custom database with one carefully built model per paralog, validated with cross-hits checked between the paralog models, answers questions official resources are not designed to answer. **Reproducible pipeline annotation.** Large annotation pipelines need annotation steps that produce byte-identical results given identical inputs. A remote web service with models that update underneath you cannot guarantee that. A versioned local database can, which is often the deciding factor in production genomics. **Screening and surveillance.** Public health, food safety, and biodefense applications often need to track a specific defined family across thousands of submitted genomes. A purpose-built database with a handful of validated models is faster, more targeted, and easier to audit than a general search against everything. **Quality control of automated annotation.** Running a custom database over a genome's predicted proteins, with models for families where misannotation is known to be rampant, catches the cases where a gene caller or a transfer-based annotator assigned a function the domain evidence does not support. **Teaching and method development.** The workflow is small enough to run entirely on a laptop, which makes it an excellent vehicle for teaching profile methods. Students can see the full chain from alignment to matrix to search result, which is invisible when the search happens behind a web form. ## Common Misconceptions About Custom Domain Databases Separate from the procedural mistakes listed in the main guide, several conceptual misunderstandings recur, and they lead to bad decisions even when every command is typed correctly. **"A profile is just a weighted consensus."** No. A consensus sequence collapses each column to one residue and then scores queries against a single sequence. A PSSM retains the full residue distribution per column and scores each query residue against that distribution. The difference is exactly why profiles detect remote homologs that consensus-based methods miss, and it is also why removing near-identical seeds changes a PSSM far less than it changes a consensus: the distribution is already an aggregate. **"The E-value is the probability that the hit is random."** It is not. The E-value is the expected number of hits with scores at least this good that you would see by chance in a database of this size and composition. It is a count, not a probability, and it depends on the search space. This is the mathematical reason, not merely a reporting convention, behind the rule that E-values from a custom database cannot be compared to E-values from the full CDD. **"Bit scores and E-values behave the same across databases."** They do not. Bit scores are normalized scores, and with consistent scaling, which is what `-scale 100.0` provides, they are far more transferable between databases than E-values. If you need one number to compare across your custom database and the official service, prefer the bit score; report the E-value only alongside the database it was computed against. **"A custom database is automatically more sensitive for my family."** Only if it is built and validated well. A carelessly built custom model can be less sensitive than the closest official model, and a moderately built one is frequently less specific. Custom means controllable, not better. The validation steps exist precisely because the default outcome is not guaranteed. **"More seeds always make a better model."** Seeds help when they add diversity and hurt when they add redundancy or non-members. Twenty carefully chosen seeds outperform a hundred scraped hits, and the alignment inspection effort per seed is the real cost. **"`rpstblastn` replaces gene prediction."** It does not. Six-frame translated search finds domain footprints on contigs, which are evidence, not gene models. Footprints that straddle a sequencing error, hits in overlapping frames, and footprints without flanking coding sequence all occur. Use translated domain hits to guide and verify gene calling, not to bypass it. **"Validation is a one-time step."** Models degrade in relevance as sequence databases grow and your taxonomic scope shifts. Any change to seeds, aligner, refinement database, or BLAST+ version invalidates prior validation. Re-run the checks with every versioned rebuild. ## Reproducibility: Versioning, Documentation, and Containers A domain database that exists only on one machine, in one directory, with parameters recalled from shell history, is not a scientific resource. Three cheap practices fix this. **Version everything, in the artifact.** The version belongs in the `makeprofiledb` title, the database directory name, the README, and ideally the filename your pipeline references. When a paper or a pipeline report cites "annotation with domain database v1.2", every parameter needed to reproduce v1.2 should be in the repository alongside the build script. **Write the README before the crisis.** A useful README for a custom database contains: the model manifest, meaning model name, seed accessions, and trimming coordinates; the exact build commands or the script itself; the aligner name and version; the refinement database and iteration count per model; the `makeprofiledb` parameters including the threshold and scale values; the BLAST+ version used to build the database; the validation summary table per model; and the date. Record checksums for the database files with `md5sum db/* > checksums.md5` so that anyone can verify they have the same artifact you validated. **Containerize the toolchain.** A conda environment file or a container image pinning BLAST+, MAFFT, and any helper scripts removes the largest source of irreproducibility, tool version drift. Combined with the build script and manifest, a container turns the entire database into a buildable artifact: one command regenerates it from seeds on any machine. Note that the database binaries themselves should still be rebuilt inside the target environment rather than copied between environments, because of the BLAST+ version compatibility issue described in the troubleshooting section. A final practice that pays off disproportionately: put the round-trip check into version control as an automated test. When a pull request adds a new model, the test builds the database and verifies that every seed recovers its model. This converts validation from a chore someone must remember into a gate that cannot be bypassed. ## Managing Database Growth: Adding Models, Retiring Models, and Knowing When to Split Databases are living resources, and the decisions you make as yours grows determine whether it stays trustworthy. **Adding a model.** Always add models through a full versioned rebuild: new seeds, alignment, PSSM, concatenate with all existing PSSMs, run `makeprofiledb` with unchanged parameters, rerun validation. Never incrementally patch a built database. The rebuild is cheap, the version bump preserves traceability of old results, and the fresh validation run rechecks the entire database rather than only the newcomer. **Retiring a model.** Retire, do not delete. A model that fails hold-out validation, or that repeatedly produces false positives on your specificity panel, should be removed from the active PSSM set but kept in the repository with a note explaining why. Silent deletion breaks result traceability; a documented retirement preserves the audit trail. **Splitting an overloaded model.** The clearest signal that one model should become two is cross-hitting between known distinct subgroups: members of subgroup A hit the subgroup B model and vice versa, usually with close scores. When that happens, split the alignment, rebuild one model per subgroup with subgroup-specific seeds, and revalidate. Expect each child model to be sharper within its subgroup and blinder outside it, which is exactly the trade you are choosing. **Merging near-duplicate models.** The opposite failure also occurs, especially after splitting zealously: two models that overlap so heavily that every query hits both. If two models agree on nearly every hit, one of them is redundant. Merge their seed sets, rebuild, and validate. **Rebuilding versus editing.** There is no editing. The inputs are alignments and parameters; the database is always derived. This discipline feels heavyweight the first two times and automatic by the fifth, and it is what makes the entire resource auditable. Every question of the form "which version of the model produced this annotation" should be answerable from the version string and the README alone. ## A Closing Checklist for Production Use Before a custom CD-Search database goes into real annotation work, confirm each item on this list: - Every model has a documented seed set, alignment provenance, and PSSM build log. - `-scale 100.0` and the CDD-convention threshold were set at build time. - Every seed recovers its own model in the round-trip check. - Hold-out validation was run, and the results are recorded per model. - A panel of unrelated proteins produces no convincing hits. - The database version appears in the title, the directory name, and the README. - The BLAST+ version used for building matches the version used for searching. - A written rule governs overlap collapsing and the precedence of custom versus official CDD models when results are merged. - The toolchain is pinned in an environment specification or container. - The round-trip validation runs automatically whenever the database is rebuilt. A database that passes this checklist is a small but genuine scientific resource: versioned, validated, and reproducible. Everything after that is routine annotation, which is exactly what you want the hard part to have been. ## Frequently Asked Questions ### Can the NCBI CD-Search web page search a custom database? No. The web service and the Batch CD-Search API search only the official Conserved Domain Database. Custom models must be built locally with `makeprofiledb` and searched locally with `rpsblast+` or `rpstblastn`. ### Which tool do I need, rpsblast or rpsblast+? Use `rpsblast+` from the BLAST+ suite, which is actively maintained and installed by conda with the `blast` package. The older legacy `rpsblast` binary is obsolete. ### How many seed sequences does a custom domain model need? There is no fixed minimum, but profiles improve substantially with an alignment of several distinct homologs that span the diversity you want to detect. Avoid padding the alignment with near-identical sequences, which adds nothing and can distort pseudocount behavior. ### Can I combine my custom models with official CDD models in one search? The cleanest approach is to build one PSSM set that contains both your models and the alignments underlying the official models you need, then compile it into a single database. Alternatively, search your custom database and the official CDD database separately and merge the results in post-processing. ### Why are my E-values different from the values on the NCBI website? E-values depend on the size and composition of the database being searched. A small custom database produces different E-values than the full CDD for the same query and model. Always report the database version with your numbers. ### Can I search nucleotide sequences against my custom domain database? Yes. Use `rpstblastn`, which translates the nucleotide query in all six frames before comparing it to the protein profiles in your database. ### What E-value cutoff should I use? Use 1e-2 while exploring results, then tighten to 1e-3 or 1e-5 for final annotation. Because profile hits are information-rich, true family members typically score far better than the cutoff, and most of your decision-making will concern borderline hits near the threshold. ### How do I update my database when I add new models? Rebuild from scratch. Generate the new PSSMs, concatenate all PSSMs including the old ones, rerun `makeprofiledb` with the same parameters, and rerun the validation checks. Version the database in its title and file name so old results remain traceable to the exact database that produced them.