Genome Annotation - Pgenerosa_v074 Using GenSAS

In our various attempts to get the Panopea generosa genome annotated in such a manner that we’re comfortable with (the previous annotation attempts we’re lacking any annotations in almost all of the largest scaffolds, which didn’t seem right), Steven stumbled across GenSAS, a web/GUI-based genome annotation program, so we gave it a shot.

This version of the genome annotation will be referred to as:

Panopea-generosa-vv0.74.a3

I uploaded the following to the GenSAS website to potentially use as “evidence files”:

Transcriptome FastA files (links to notebook entries):

TransDecoder protein FastA files (links to notebook entries)

Repeats Files

RepeatModeler library

RESULTS

This took way longer than I was expecting! This took nearly an entire month (the majority of that time was running Augustus ab initio gene prediction, which took ~3 weeks):

GenSAS project summary processes and runtimes screencap

Output folder:

20190710_Pgenerosa_v074_gensas_annotation/

Merged GFF (SwissProt IDs in Column 9 - from BLASTp and DIAMOND):

Panopea-generosa-vv0.74.a3-merged-2019-09-24-9-20-04.gff3

Feature counts:

awk 'NR>3 { print $3 }' Panopea-generosa-vv0.74.a3-merged-2019-09-03-6-14-33.gff3 | sort | uniq -c

 192022 CDS
 192022 exon
  45748 gene
  45748 mRNA

BUSCO assessment:

68.4% complete BUSCOs present in predicted genes

Individual feature GFFs were made with the following shell commands:

features_array=(CDS exon gene mRNA)

for feature in ${features_array[@]}
do
output="Panopea-generosa-vv0.74.a3.${feature}.gff3"
input="Panopea-generosa-vv0.74.a3-merged-2019-09-24-9-20-04.gff3"
head -n 3 Panopea-generosa-vv0.74.a3-merged-2019-09-24-9-20-04.gff3 \
>> ${output}
awk -v feature="$feature" '$3 == feature {print}' ${input} \
>> ${output}
done

SwissProt functional annotations (tab-delimited text):

Pfam annotations (tab-delimited text):

Panopea-generosa-vv0.74.a3.5d65aaa449919-pfam.tab

Grabbed the top 10 most abundant Pfam Accessions to see how things looked:

awk '{print $2}' Panopea-generosa-vv0.74.a3.5d65aaa449919-pfam.tab | sort | uniq -c | sort -nr | head

Feature Count	Pfam Accession	Pfam
1062	PF00078.22	Reverse transcriptase (RNA-dependent DNA polymerase)
370	PF00665.21	Integrase
353	PF13358.1	DDE superfamily endonuclease
264	PF03372.18	Endonuclease/Exonuclease/phosphatase family
213	PF14529.1	Endonuclease-reverse transcriptase
208	PF00643.19	B-box zinc finger
199	PF07690.11	Major Facilitator Superfamily
158	PF00001.16	7 transmembrane receptor (rhodopsin family)
148	PF12796.2	Ankyrin repeat
136	PF00096.21	Zinc finger

A couple of interesting things that I notice from this table:

The four of the top five most abundant are involved in DNA transposition.
A rhodopsin protein family appears in the Top 10 most abundant Pfams?! Proteins in this family are involved in light detection…

InterProScan annotations (tab-delimited text):

Panopea-generosa-vv0.74.a3.5d65aaa22961e-interproscan.tab
- Contains gene ontology (GO) terms

Project Summary file (TEXT):

t_5d25021b0d20b-5d6843e21278a-publish-project-summary.txt

=================================
 Project Summary
---------------------------------
# Project Information
  Project Name         : Pgenerosa_v074
  Owner                : kubu4
  Create Date          : 2019-07-09 14:07:39

# Project Properties
  Genus                : Panopea
  Species              : generosa
  Project Type         : invertebrate
  Prefix               : PGEN_
  Common Name          : Pacific geoduck
  Genetic Code         : Standard Code

# Input FASTA
  Filename           : Pgenerosa_v074.fa
  Filesize           : 913.68 MB
  Number of Sequence : 18

=================================
 Job Information
---------------------------------
# Official Gene Set
  >PASA Refinement
  - version : 2.3.3
  - Transcripts FASTA file : Trinity.fasta

  # The source Job of the refinement job
    >Augustus-01
    - version : 3.3.1
    - Species : fly
    - Report genes on : both
    - Allowed gene structure : partial
    - cDNA (transcripts) sequences : Trinity.fasta
    - Protein sequences : 20180827_trinity_geoduck.fasta.transdecoder.fa


# The consensus mask Job
  >Masked Repeat Consensus

  # The source jobs for consensus mask job
    >RepeatMasker
    >RepeatModeler

  # Family copy number summary
    Family	Copy Numbers
    DNA	675
    DNA/Academ	1327
    DNA/Crypton	344
    DNA/Ginger	130
    DNA/Kolobok-T2	141
    DNA/Maverick	942
    DNA/MuLE-MuDR	201
    DNA/MuLE-NOF?	142
    DNA/P	167
    DNA/PIF-Harbinger	227
    DNA/RC	587
    DNA/Sola	508
    DNA/TcMar-Fot1	117
    DNA/TcMar-Mariner	6734
    DNA/TcMar-Tc1	3718
    DNA/hAT-Tip100	516
    DNA/hAT-hAT5	1037
    Type:DNA	17513
    LINE	883
    LINE/CR1	5204
    LINE/CR1-Zenon	14653
    LINE/I	980
    LINE/I-Nimb	1119
    LINE/L1	4031
    LINE/L1-Tx1	6620
    LINE/L2	8879
    LINE/L2-Hydra	113
    LINE/Penelope	1026
    LINE/RTE-X	21214
    Type:LINE	64722
    SINE/tRNA-Core-L2	41152
    Type:SINE	41152
    LTR/Caulimovirus	140
    LTR/DIRS	448
    LTR/Gypsy	1031
    LTR/Ngaro	343
    LTR/Pao	82
    Type:LTR	2044
    Type:EVERYTHING_TE	125431
    Type:Simple_repeat	19235
    Type:Unknown	1465471

# The functional Jobs on the OGS
  >InterProScan
  - version : 5.29-68.0

  >Pfam
  - version : 1.6
  - E-value Sequence : 1
  - E-value Domain : 10

  >SignalP
  - version : 4.1
  - Organism group : euk
  - Method : best
  - D-cutoff for SignalP-noTM networks : 0.45
  - D-cutoff for SignalP-TM networks : 0.50
  - Minimal predicted signal peptide length : 10
  - Truncate to sequence length : 70

  >BLAST protein vs protein (blastp)_SP01
  - version : 2.7.1
  - Protein Data Set : SwissProt
  - Maximum HSP Distinace : 30000
  - Output type : tab
  - Matrix : BLOSUM62
  - Expect : 1e-8
  - Word Size : 3
  - Gap Open : 11
  - Gap Extend : 1

  >DIAMOND Functional_SP01
  - version : 0.9.22
  - Protein Data Set : SwissProt

  >BLAST protein vs protein (blastp)
  - version : 2.7.1
  - Protein Data Set : 20180827_trinity_geoduck.fasta.transdecoder.fa
  - Maximum HSP Distinace : 30000
  - Output type : tab
  - Matrix : BLOSUM62
  - Expect : 1e-8
  - Word Size : 3
  - Gap Open : 11
  - Gap Extend : 1

  >DIAMOND Functional
  - version : 0.9.22
  - Protein Data Set : 20180827_trinity_geoduck.fasta.transdecoder.fa

Overall, this annotation is much more believable than the previous MAKER annotations, due to the fact that GenSAS actually predicts genes to exist on all the scaffolds (unlike MAKER)! As such, this will likely become the canonical P.generosa genome going forward. With that being said, we still should manually curate this when we have the time to see how well the predictions line up with the evidence.