predictions with RNA evidence

[xvazquezc]

Hi there,

I'm annotating a genome and I have RNAseq data. I already have assembled transcripts (via Trinity) and I was wondering if there was any preference and/or differences if using the --transcripts or --bam with funannotate2 predict.

In the same line, are the --weights applicable to the evidences or just to the gene predictors? I've only seen how to specify them in #43 , e.g. here. So wondering about the appropriate usage.

Cheers,
Xabi

[xvazquezc]

Hi again,

After overcoming the issue described in #93, I ran the prediction to completion like this (btw, conda install with funannotate2 v26.6.21 in Rocky Linux):

funannotate2 predict --cpus $NCPUS --memory-limit 120 \
--fasta "${MSKGENOME}" \
--species "${SPECIES}" \
--strain $ID \
--locus-tag $ID \
--input-dir $TRAINOUT \
--external $HLXROUT \
--proteins ${FUNANNOTATE2_DB}/uniprot_sprot.fasta \
--transcripts $TRANSCRIPTS \
--out $PREDICTOUT

This was after running funannotate2 train (note the genome is basically T2T):

funannotate2 train \
--fasta "${MSKGENOME}" \
--species "${SPECIES}" \
--min-contig-length 50000 \
--genemark-mode guided \
--augustus-species aspergillus_fumigatus \
--busco-lineage dothideomycetes \
--strain $ID \
--out $TRAINOUT \
--cpus $NCPUS

Looking at the resuts I was wondering how the model weighting is managed. For example, these were the models scoring (tbh, surprised by the low values for gapmm2-gene:

{
  "augustus": {
    "busco": 0.9813542688910697,
    "train": {
      "tool": "augustus",
      "model": "908c72a4-3a9a-419b-ba56-9b3f1f5c37ab",
      "n_test_genes": 200,
      "ref_genes_found": 200,
      "ref_genes_missed": 0,
      "extra_query_genes": 154,
      "average_aed": 0.11906437352332498,
      "nucleotide_sensitivity": 0.9071259860882048,
      "nucleotide_precision": 0.9068441978367979,
      "exon_sensitivity": 0.6522842639593909,
      "exon_precision": 0.6594085271750247,
      "gene_sensitivity": 1.0,
      "gene_precision": 0.3418803418803419
    }
  },
  "gapmm2-gene": {
    "busco": 0.48380765456329733
  },
  "genemark": {
    "busco": 0.9960745829244357,
    "train": {
      "tool": "genemark",
      "model": "gmhmm.mod",
      "n_test_genes": 200,
      "ref_genes_found": 200,
      "ref_genes_missed": 0,
      "extra_query_genes": 217,
      "average_aed": 0.10446726703519064,
      "nucleotide_sensitivity": 0.8814030905419991,
      "nucleotide_precision": 0.9487698187664156,
      "exon_sensitivity": 0.6363636363636364,
      "exon_precision": 0.6772126022126022,
      "gene_sensitivity": 1.0,
      "gene_precision": 0.2668918918918919
    }
  },
  "glimmerhmm": {
    "busco": 0.9283611383709519,
    "train": {
      "tool": "glimmerhmm",
      "model": "train",
      "n_test_genes": 200,
      "ref_genes_found": 200,
      "ref_genes_missed": 0,
      "extra_query_genes": 147,
      "average_aed": 0.12894121641414472,
      "nucleotide_sensitivity": 0.8979005764266919,
      "nucleotide_precision": 0.8888833915488055,
      "exon_sensitivity": 0.6096938775510204,
      "exon_precision": 0.6076166180758017,
      "gene_sensitivity": 1.0,
      "gene_precision": 0.30660377358490565
    }
  },
  "helixer": {
    "busco": 0.9990186457311089
  },
  "miniprot-gene": {
    "busco": 0.13052011776251227
  "snap": {
    "busco": 0.9087340529931305,
    "train": {
      "tool": "snap",
      "model": "snap-trained.hmm",
      "n_test_genes": 200,
      "ref_genes_found": 200,
      "ref_genes_missed": 0,
      "extra_query_genes": 241,
      "average_aed": 0.12325902516548182,
      "nucleotide_sensitivity": 0.8673618102316233,
      "nucleotide_precision": 0.9277418342864765,
      "exon_sensitivity": 0.6005154639175257,
      "exon_precision": 0.6240549828178694,
      "gene_sensitivity": 1.0,
      "gene_precision": 0.22258064516129034
    }
  }
}

Led to the following weights:

{"augustus-hiq": 4, "genemark": 3, "helixer": 2, "augustus": 2, "gapmm2-gene": 1, "snap": 1, "glimmerhmm": 1, "miniprot-gene": 1}

And consensus gene models:

[["genemark", 5340], ["augustus-hiq", 2053], ["helixer", 370], ["augustus", 199], ["gapmm2-gene", 179], ["glimmerhmm", 178], ["snap", 116], ["miniprot-gene", 1]]

I don't know if the low weighting of the transcript evidence is bc of the way the transcripts were provided (Trinity-assembled transcripts) but I found it a bit surprising to see the best evidence source so lowly scored. Also, HelixerLite predictions had the highest busco score (0.999...). Wondering if there are fewer parameters to evaluate the sources/models when they don't come from the training models generated with funannotate2 train.

I leave the funannotate2 predict log in case it helps

[nextgenusfs]

The scores are generated from an evaluation of BUSCO models. I wouldn't necessarily expect the Trinity transcripts to capture all of the BUSCO gene models, ie they would all need to be decently expressed to be properly assembled/mapped. The evidence is also used at consensus prediction step to guide the choice of the model at each locus.

Since you passed in helixerlite results via the external models, it doesn't have the same training data so it can't use the exact same information to automatically determine the weighting, I'd probably set it to 4, you can do that with the CLI --weights helixer:4.

[xvazquezc]

By changing the HelixerLite weights, other than most of the gene predictions being selected from that model, there is a big drop on the UTR regions predicted compared to not setting any weights.

EDIT: I also see a drop in the predicted UTRs when providing the bam file with the mapped transcripts. Not sure how much of a big deal it might be but I thought worth pointing out.