Protein domains on gene models

JBrowse v5 beta

The tutorials target the JBrowse v5 beta, and the v4.3.0 release on the download page lacks some of what they show. To install the beta of JBrowse Web, run npm install -g @jbrowse/cli@next, then jbrowse create jbrowse2 --branch v5.0.0-beta.17. Desktop beta builds are coming soon.

A gene model shows where the exons are, and Pfam names the protein domains they encode, such as a kinase or a DNA-binding fold. Domain callers report each domain in residues along the protein, so the domains need placing back on the exons before a gene track can draw them. Done for every transcript, the placement shows which exons encode each domain and which isoforms leave one out. We work through the honey bee RefSeq annotation and the vav gene (LOC413488), whose isoforms differ at their first domain. We:

translate the gene models and call Pfam domains on the proteins with hmmscan
place each domain on the exons that encode it with jbrowse make-protein-features
draw the domains on the gene track's isoforms, colored by Pfam family

Prerequisites

a JBrowse to paste the tracks into (Web or Desktop), or the honey bee genome on genomes.jbrowse.org
gffread, which translates each transcript's CDS
HMMER (hmmscan, hmmpress), which calls the domains
htslib (bgzip, tabix), which indexes the domain file
node, for the JBrowse CLI
samtools, for the FASTA index the build script writes

# Debian/Ubuntu
sudo apt install gffread hmmer tabix samtools nodejs
# or conda, which has every tool above
conda install -c bioconda -c conda-forge gffread hmmer htslib samtools \
  nodejs
# any platform
npm install -g @jbrowse/cli

The data

The figures read the honey bee Amel_HAv3.1 assembly and its NCBI RefSeq annotation from genomes.jbrowse.org, and a Pfam domain file built from that annotation with the commands below. Nothing to download to follow the page: the track configs read the domain file and the gene models by URL. The commands take your own genome FASTA and GFF3 the same way.

What it is	File
The honey bee genome, which gffread reads the CDS from	GCF_003254395.2_Amel_HAv3.1_genomic.fna.gz
Its RefSeq gene models, which the domains are placed on	GCF_003254395.2_Amel_HAv3.1_genomic.gff.gz
The Pfam HMM library, about 400 MB, which hmmscan searches	Pfam-A.hmm.gz

Calling Pfam domains on the proteins

We'll start from the two files any annotation project has, the genome FASTA and the gene models, and translate every transcript's CDS into a protein FASTA:

# -y writes each protein named by its transcript ID; -S spells a stop codon `*`,
# which HMMER reads and gffread's default `.` is not. gffread wants the FASTA
# uncompressed.
gffread -S -y proteins.fa -g genome.fa genes.gff3

hmmscan searches each protein against every Pfam family and writes one row per domain hit, in residue coordinates along the protein. The library needs pressing once:

hmmpress Pfam-A.hmm

# --cut_ga: each Pfam family's curated gathering threshold, the cutoff InterPro
# uses, in place of one E-value for every family. hmmscan keeps about two cores
# busy whatever --cpu asks for, so the build script splits proteins.fa into
# chunks and scans them at once.
hmmscan --cut_ga --cpu 2 --noali --domtblout pfam.domtblout Pfam-A.hmm proteins.fa >/dev/null

An annotation pipeline that already ran InterProScan has the same calls in its TSV output (-f tsv), which the next step reads in place of pfam.domtblout.

Placing the domains on the exons

jbrowse make-protein-features walks each transcript's CDS in the direction of translation and places each domain on the bases that encode it, on either strand and across introns. It writes a bgzipped, indexed BED with one row per domain per transcript:

# matches each protein to its transcript by transcript ID, CDS protein_id or
# CDS ID, so gffread, NCBI and Ensembl protein names all work. Where HMMER hits
# on one protein overlap, the lowest E-value is kept.
jbrowse make-protein-features genes.gff3 pfam.domtblout --out pfam.bed.gz

A domain crossing an intron is one row with a block per exon, as BED12 draws a spliced transcript.
The row's name is the Pfam family; accession, description, transcript, proteinStart, proteinEnd and evalue follow as named columns, and JBrowse shows them on click.
--source Pfam keeps only the Pfam rows of an InterProScan TSV.

Drawing the domains on the genes

JBrowse loads the domain file as a track, with disableGeneHeuristic to stop it reading a stranded BED12 row as a coding transcript:

Config file

Goes in the tracks array of config.json. See Tracks.

{
  "type": "FeatureTrack",
  "trackId": "pfam_domains",
  "name": "Pfam domains",
  "assemblyNames": ["GCF_003254395.2"],
  "adapter": {
    "type": "BedTabixAdapter",
    "uri": "https://jbrowse.org/demos/protein_domains/pfam.bed.gz",
    "disableGeneHeuristic": true
  }
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "FeatureTrack",
  "trackId": "pfam_domains",
  "name": "Pfam domains",
  "assemblyNames": ["GCF_003254395.2"],
  "adapter": {
    "type": "BedTabixAdapter",
    "uri": "https://jbrowse.org/demos/protein_domains/pfam.bed.gz",
    "disableGeneHeuristic": true
  }
}'

In the app

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:

{
  "type": "FeatureTrack",
  "trackId": "pfam_domains",
  "name": "Pfam domains",
  "assemblyNames": ["GCF_003254395.2"],
  "adapter": {
    "type": "BedTabixAdapter",
    "uri": "https://jbrowse.org/demos/protein_domains/pfam.bed.gz",
    "disableGeneHeuristic": true
  }
}

A gene track names that track in proteinFeatures and draws each domain over the CDS that encodes it, colored by Pfam family, with a key naming the families on screen. This one reads the same RefSeq GFF as the honey bee hub's gene track:

Config file

Goes in the tracks array of config.json. See Tracks.

{
  "type": "FeatureTrack",
  "trackId": "genes_with_pfam",
  "name": "RefSeq genes with Pfam domains",
  "assemblyNames": ["GCF_003254395.2"],
  "adapter": {
    "type": "Gff3TabixAdapter",
    "uri": "https://jbrowse.org/hubs/genark/GCF/003/254/395/GCF_003254395.2/GCF_003254395.2_Amel_HAv3.1_genomic.gff.gz",
    "index": {
      "location": {
        "uri": "https://jbrowse.org/hubs/genark/GCF/003/254/395/GCF_003254395.2/GCF_003254395.2_Amel_HAv3.1_genomic.gff.gz.csi"
      },
      "indexType": "CSI"
    }
  },
  "displayDefaults": {
    "proteinFeatures": {
      "track": "pfam_domains",
      "color": { "field": "name", "legend": { "title": "Pfam family" } }
    }
  }
}

CLI (add-track-json)

jbrowse add-track-json '{
  "type": "FeatureTrack",
  "trackId": "genes_with_pfam",
  "name": "RefSeq genes with Pfam domains",
  "assemblyNames": ["GCF_003254395.2"],
  "adapter": {
    "type": "Gff3TabixAdapter",
    "uri": "https://jbrowse.org/hubs/genark/GCF/003/254/395/GCF_003254395.2/GCF_003254395.2_Amel_HAv3.1_genomic.gff.gz",
    "index": {
      "location": {
        "uri": "https://jbrowse.org/hubs/genark/GCF/003/254/395/GCF_003254395.2/GCF_003254395.2_Amel_HAv3.1_genomic.gff.gz.csi"
      },
      "indexType": "CSI"
    }
  },
  "displayDefaults": {
    "proteinFeatures": {
      "track": "pfam_domains",
      "color": { "field": "name", "legend": { "title": "Pfam family" } }
    }
  }
}'

In the app

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:

{
  "type": "FeatureTrack",
  "trackId": "genes_with_pfam",
  "name": "RefSeq genes with Pfam domains",
  "assemblyNames": ["GCF_003254395.2"],
  "adapter": {
    "type": "Gff3TabixAdapter",
    "uri": "https://jbrowse.org/hubs/genark/GCF/003/254/395/GCF_003254395.2/GCF_003254395.2_Amel_HAv3.1_genomic.gff.gz",
    "index": {
      "location": {
        "uri": "https://jbrowse.org/hubs/genark/GCF/003/254/395/GCF_003254395.2/GCF_003254395.2_Amel_HAv3.1_genomic.gff.gz.csi"
      },
      "indexType": "CSI"
    }
  },
  "displayDefaults": {
    "proteinFeatures": {
      "track": "pfam_domains",
      "color": { "field": "name", "legend": { "title": "Pfam family" } }
    }
  }
}

For your own data, point the first uri at your pfam.bed.gz, the gene adapter at your sorted, indexed GFF3, and assemblyNames at your assembly.

To follow along on genomes.jbrowse.org, open the honey bee genome, add both tracks, and go to NC_037640.1:1,971,400-1,976,800. Clicking a domain shows its transcript, accession and residue range.

Honey bee vav: its RefSeq isoforms with their Pfam domains drawn over the CDS, colored by family. The top isoform skips the exon encoding CH.

Reproduce it end to end

build_protein_domains.sh runs every step above on your own genome FASTA and GFF3. The script:

Translates each transcript's CDS with gffread.
Downloads and presses the Pfam library once, then calls domains with hmmscan, on chunks of the proteins at once.
Places the domains on the exons with jbrowse make-protein-features.
Builds a JBrowse whose gene track draws the domains.

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_protein_domains.sh
bash build_protein_domains.sh genome.fa genes.gff3   # builds ./protein_domains_build/jbrowse2
npx --yes serve protein_domains_build/jbrowse2       # then open the printed URL

The script reads genome.fa and genes.gff3 plain or gzipped, and THREADS=32 bash build_protein_domains.sh ... scans more chunks at once.

Citations

Mistry et al. (2021). Pfam: The protein families database in 2021
Eddy (2011). Accelerated profile HMM searches
Pertea and Pertea (2020). GFF utilities: GffRead and GffCompare
Wallberg et al. (2019). A hybrid de novo genome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds

Feedback on this tutorial is welcome: contact us.

Protein domains on gene models

Prerequisites#

The data#

Calling Pfam domains on the proteins#

Placing the domains on the exons#

Drawing the domains on the genes#

Reproduce it end to end#

See also#

Citations#