---
author: Sam White
toc-title: Contents
toc-depth: 5
toc-location: left
title: Transcriptome Annotation - M.magister De Novo Transcriptome Assembly Using Trinotate on Mox
date: '2023-11-09'
draft: false
engine: knitr
code-fold: true
execute:
  freeze: auto
categories: 
  - mox
  - 2023
  - Metacarcinus magister
  - Dungenss crab
  - Trinotate
  - transcriptome
  - annotation
---
# Intro

As part of [transcriptome annotation](https://github.com/laurahspencer/DuMOAR/issues/44) (GitHub issue) for [DuMOAR](https://github.com/laurahspencer/DuMOAR) (GitHub repo), I've already run [`TransDecoder`](https://https://github.com/TransDecoder/TransDecoder/wiki.com/Trinotate/Trinotate/wiki) ( run on [20231107](../2023-11-07-Transcriptome-Annotation---M.magister-De-Novo-Transcriptome-Assembly-for-DuMOAR-Project-Using-TransDecoder-on-Mox/index.qmd)), as well as [`DIAMOND`](https://github.com/bbuchfink/diamond) BLASTx (run  [earlier today](../2023-11-09-Transcriptome-Annotation---M.magister-De-Novo-Transcriptome-Assembly-for-DuMOAR-Project-Using-DIAMOND-BLASTx-on-Mox/index.qmd)). The final part of getting the transcriptome annotated was to run [`Trinotate`](https://github.com/Trinotate/Trinotate/wiki). [`Trinotate`](https://github.com/Trinotate/Trinotate/wiki) uses a combination of tools to generate a pretty thorough annotation of transcriptomes. This was run on Mox.

## SLURM script

- [20231109-mmag-trinotate-transcriptome.sh](https://github.com/RobertsLab/sams-notebook/blob/master/sbatch_scripts/20231109-mmag-trinotate-transcriptome.sh) (GitHub)

```bash
#!/bin/bash
## Job Name
#SBATCH --job-name=20231109-mmag-trinotate-transcriptome
## Allocation Definition
#SBATCH --account=srlab
#SBATCH --partition=srlab
## Resources
## Nodes
#SBATCH --nodes=1
## Walltime (days-hours:minutes:seconds format)
#SBATCH --time=4-00:00:00
## Memory per node
#SBATCH --mem=120G
##turn on e-mail notification
#SBATCH --mail-type=ALL
#SBATCH --mail-user=samwhite@uw.edu
## Specify the working directory for this job
#SBATCH --chdir=/gscratch/scrubbed/samwhite/outputs/20231109-mmag-trinotate-transcriptome


# Script to run Trinotate on M.magister transcriptome.


###################################################################################
# These variables need to be set by user

# Input files
## BLASTx
blastx_out="/gscratch/scrubbed/samwhite/outputs/20231109-mmag-transcriptome-diamond-blastx/trinity_denovo.Trinity.blastx.outfmt6"

## TransDecoder
transdecoder_dir="/gscratch/scrubbed/samwhite/outputs/20231107-mmag-transdecoder-transcriptome"
blastp_out="${transdecoder_dir}/blastp_out/blastp.outfmt6"
pfam_out="${transdecoder_dir}/pfam_out/pfam.domtblout"
lORFs_pep="${transdecoder_dir}/trinity_denovo.Trinity.fasta.transdecoder_dir/longest_orfs.pep"

## Transcriptomics
transcriptomes_dir="/gscratch/srlab/sam/data/M_magister/transcriptomes"
trinity_fasta="${transcriptomes_dir}/trinity_denovo.Trinity.fasta"
trinity_gene_map="${transcriptomes_dir}/trinity_denovo.Trinity.fasta.gene_trans_map"

###################################################################################

# Exit script if any command fails
set -e

# Load Python Mox module for Python module availability

module load intel-python3_2017

# SegFault fix?
export THREADS_DAEMON_MODEL=1


wd="$(pwd)"
timestamp=$(date +%Y%m%d)


## Paths to input/output files

## New folders for working directory
rnammer_out_dir="${wd}/RNAmmer_out"
signalp_out_dir="${wd}/signalp_out"
tmhmm_out_dir="${wd}/tmhmm_out"


rnammer_prefix=${trinity_fasta##*/}
prefix="${timestamp}.${rnammer_prefix}.trinotate"

# Output files
rnammer_out="${rnammer_out_dir}/${rnammer_prefix}.rnammer.gff"
signalp_out="${signalp_out_dir}/${prefix}.signalp.out"
tmhmm_out="${tmhmm_out_dir}/${prefix}.tmhmm.out"
trinotate_report="${wd}/${prefix}_annotation_report.txt"

# Paths to programs
rnammer_dir="/gscratch/srlab/programs/RNAMMER-1.2"
rnammer="${rnammer_dir}/rnammer"
signalp_dir="/gscratch/srlab/programs/signalp-4.1"
signalp="${signalp_dir}/signalp"
tmhmm_dir="/gscratch/srlab/programs/tmhmm-2.0c/bin"
tmhmm="${tmhmm_dir}/tmhmm"
trinotate_dir="/gscratch/srlab/programs/Trinotate-v3.1.1"
trinotate="${trinotate_dir}/Trinotate"
trinotate_rnammer="${trinotate_dir}/util/rnammer_support/RnammerTranscriptome.pl"
trinotate_GO="${trinotate_dir}/util/extract_GO_assignments_from_Trinotate_xls.pl"
trinotate_features="${trinotate_dir}/util/Trinotate_get_feature_name_encoding_attributes.pl"
trinotate_sqlite_db="Trinotate.sqlite"

# Generate FastA checksum, for reference if needed.
md5sum ${trinity_fasta} > fasta.checksum.md5

# Make output directories
mkdir "${rnammer_out_dir}" "${signalp_out_dir}" "${tmhmm_out_dir}"

# Copy sqlite database template

cp ${trinotate_dir}/admin/Trinotate.sqlite .

# Run signalp
${signalp} \
-f short \
-n "${signalp_out}" \
${lORFs_pep}

# Run tmHMM
${tmhmm} \
--short \
< ${lORFs_pep} \
> "${tmhmm_out}"

# Run RNAmmer
# Has "||" operator due to previous lack of matches
# Need "||" to continue with annotation.
cd "${rnammer_out_dir}" || exit
${trinotate_rnammer} \
--transcriptome ${trinity_fasta} \
--path_to_rnammer ${rnammer}

cd "${wd}" || exit

# Run Trinotate
## Load transcripts and coding regions into database
${trinotate} \
${trinotate_sqlite_db} \
init \
--gene_trans_map "${trinity_gene_map}" \
--transcript_fasta "${trinity_fasta}" \
--transdecoder_pep "${lORFs_pep}"

## Load BLAST homologies
"${trinotate}" \
"${trinotate_sqlite_db}" \
LOAD_swissprot_blastp \
"${blastp_out}"

"${trinotate}" \
"${trinotate_sqlite_db}" \
LOAD_swissprot_blastx \
"${blastx_out}"

## Load Pfam
"${trinotate}" \
"${trinotate_sqlite_db}" \
LOAD_pfam \
"${pfam_out}"

## Load transmembrane domains
"${trinotate}" \
"${trinotate_sqlite_db}" \
LOAD_tmhmm \
"${tmhmm_out}"

## Load signal peptides
"${trinotate}" \
"${trinotate_sqlite_db}" \
LOAD_signalp \
"${signalp_out}"

## Load RNAmmer
"${trinotate}" \
"${trinotate_sqlite_db}" \
LOAD_rnammer \
"${rnammer_out}"

## Creat annotation report
"${trinotate}" \
"${trinotate_sqlite_db}" \
report \
> "${trinotate_report}"

# Extract GO terms from annotation report
"${trinotate_GO}" \
--Trinotate_xls "${trinotate_report}" \
-G \
--include_ancestral_terms \
> "${prefix}".go_annotations.txt

# Make transcript features annotation map
"${trinotate_features}" \
"${trinotate_report}" \
> "${prefix}".annotation_feature_map.txt

####################################################################

# Capture program options
if [[ "${#programs_array[@]}" -gt 0 ]]; then
  echo "Logging program options..."
  for program in "${!programs_array[@]}"
  do
    {
    echo "Program options for ${program}: "
    echo ""
    # Handle samtools help menus
    if [[ "${program}" == "samtools_index" ]] \
    || [[ "${program}" == "samtools_sort" ]] \
    || [[ "${program}" == "samtools_view" ]]
    then
      ${programs_array[$program]}

    # Handle DIAMOND BLAST menu
    elif [[ "${program}" == "diamond" ]]; then
      ${programs_array[$program]} help

    # Handle NCBI BLASTx menu
    elif [[ "${program}" == "blastx" ]]; then
      ${programs_array[$program]} -help

    # Handle StringTie prepDE script
    elif [[ "${program}" == "prepDE" ]]; then
      python3 ${programs_array[$program]} -h
    fi
    ${programs_array[$program]} -h
    echo ""
    echo ""
    echo "----------------------------------------------"
    echo ""
    echo ""
  } &>> program_options.log || true

    # If MultiQC is in programs_array, copy the config file to this directory.
    if [[ "${program}" == "multiqc" ]]; then
      cp --preserve ~/.multiqc_config.yaml multiqc_config.yaml
    fi
  done
  echo "Finished logging programs options."
  echo ""
fi


# Document programs in PATH (primarily for program version ID)
echo "Logging system $PATH..."
{
date
echo ""
echo "System PATH for $SLURM_JOB_ID"
echo ""
printf "%0.s-" {1..10}
echo "${PATH}" | tr : \\n
} >> system_path.log
echo "Finished logging system $PATH."
```

/----------

# RESULTS

## Runtime

Runtime was just under 3.5hrs:

![Screencap showing Trinotate runtime on Mox of 3hrs, 21mins, 10secs](20231109-mmag-trinotate-transcriptome-runtime.png)

## Files

Output directory:

- [20231109-mmag-trinotate-transcriptome/](https://gannet.fish.washington.edu/Atumefaciens/20231109-mmag-trinotate-transcriptome/)

### Annotation Feature Map
Annotation feature map (10MB; TXT):

- [20231109.trinity_denovo.Trinity.fasta.trinotate.annotation_feature_map.txt](https://gannet.fish.washington.edu/Atumefaciens/20231109-mmag-trinotate-transcriptome/20231109.trinity_denovo.Trinity.fasta.trinotate.annotation_feature_map.txt)

  - [This can be used to update Trinity-based gene expression matrices like so](https://github.com/trinityrnaseq/trinityrnaseq/wiki/Functional-Annotation-of-Transcripts):

    - ```${TRINITY_HOME}/Analysis/DifferentialExpression/rename_matrix_feature_identifiers.pl Trinity_trans.counts.matrix annot_feature_map.txt > Trinity_trans.counts.wAnnot.matrix```

### Gene Ontology (GO) Annotations
Gene ontology (GO) annotations (26MB; TXT):

- [20231109.trinity_denovo.Trinity.fasta.trinotate.go_annotations.txt](https://gannet.fish.washington.edu/Atumefaciens/20231109-mmag-trinotate-transcriptome/20231109.trinity_denovo.Trinity.fasta.trinotate.go_annotations.txt)

```{r GO-annotations, eval=TRUE}
#| label: tbl-GO-annotations
#| tbl-cap: "Gene Ontology annotations"
#| tbl-colwidths: [60,40]
#| warning: false
library(kableExtra)
library(dplyr)

# Read in data frame
go.annotations.df <- read.csv("https://gannet.fish.washington.edu/Atumefaciens/20231109-mmag-trinotate-transcriptome/20231109.trinity_denovo.Trinity.fasta.trinotate.go_annotations.txt",
header = FALSE,
sep = "\t")

# Provide column names
names(go.annotations.df)[1] <- "Gene.ID"
names(go.annotations.df)[2] <- "GO.IDs"

# Create nice HTML table
knitr::kable(head(go.annotations.df))

```
### Full Annotations
Annotation report (107MB; CSV):

- [20231109.trinity_denovo.Trinity.fasta.trinotate_annotation_report.txt](https://gannet.fish.washington.edu/Atumefaciens/20231109-mmag-trinotate-transcriptome/20231109.trinity_denovo.Trinity.fasta.trinotate_annotation_report.txt)

This table is tab-delimited, but uses `^` for intra-field delimitation:


```{r annotation-report, cache=TRUE}
#| label: tbl-annotation-report
#| tbl-cap: "Full annotations"
#| tbl-colwidths: [60,40]
library(kableExtra)

annotations.df <- read.csv("https://gannet.fish.washington.edu/Atumefaciens/20231109-mmag-trinotate-transcriptome/20231109.trinity_denovo.Trinity.fasta.trinotate_annotation_report.txt",
header = TRUE,
sep = "\t")

knitr::kable(head(annotations.df))

```



```{r annotation-counts}
#| cache: true
#| dependson: "annotation-report"
number_of_annotations <- n_distinct(annotations.df$X.gene_id)
```

Number of Trinity genes with annotations: `r number_of_annotations`

### SQlite Database
SQlite database (797MB; SQLITE):

- [Trinotate.sqlite](https://gannet.fish.washington.edu/Atumefaciens/20231109-mmag-trinotate-transcriptome/Trinotate.sqlite)