Question about genome annotation QC

Hi,
Thanks for putting this tutorial online it is very helpful.
I am trying to annotate a genome and I'd love to have your input and suggestions.

My PI sequenced a marine mammal using PacBio, but the reads came short ~ 5 kb, know, shocking, and low input ~ 40 GB for a 2.4 GB genome. Thus 2 runs were done, total ~ 90 GB. The assembly done with HiFiasm, after haplotype purging and decontamination, not so bad:
N of seqs: 397,
length: min : 6,208, average: 6,032,705.3, max: 59,387,461.
The BUSCO Score was pretty good:
C:98.9%[S:97.6%,D:1.3%],F:0.3%,M:0.8%,n:13727,E:2.8%
AFter this, I scaffolded it using an individual sequenced and scaffolded with HiC by DNAZoo/same species.

Here is my issue: finding a "good & reliable & complete" annotation.
My PI sequenced 4 tissues, regular RNA-seq. I used EviAnn, with protein data: Uniref90- carnivore, BUSCO-db mammals, 10 closely related species, 10 non-carnivore ref species, and transcript data from closely related species, as suggested.
The BUSCO proteome is lower:
C:90.8%[S:25.9%,D:65.0%],F:2.1%,M:7.1%,n:13727
This is with full proteome, I used agat to get the longest isoform per gene, and the BUSCO was lower:
C:89.1%[S:88.4%,D:0.8%],F:2.7%,M:8.2%,n:13727

I asked Aleksey Zimin, author of EviAnn, about the lower proteome score, and he replied this:
"BUSCO score for proteome will never be as high as for the genome. For a protein to be annotated it must have an intact ORF. BUSCO counts an ortholog to be a match even if only a small part of the protein aligns to the genome translated in all six frames, allowing for in-frame stop codons, not requiring ORF completeness. The organism may not even have a valid functional protein that matches the ortholog in BUSCO, but it may have a small remnant of that protein in the DNA as pseudogene that lost function."

Do you know another way to get the longest isoform per gene different than with Agat?
I can't believe that with the longest isoform per gene the BUSCO score is lower.

Then, I used the new released AI-annotator/ Tiberius, which in terms of BUSCO got higher scores
C:95.5%[S:94.7%,D:0.8%],F:1.1%,M:3.4%,n:13727
same with longest isoform per gene:
C:95.5%[S:94.7%,D:0.8%],F:1.1%,M:3.4%,n:13727

One thing is the BUSCO Score, the other is the number of genes identified.
EviAnn found: 21283, while Tiberius found: 19564. The difference in number of genes seems not unsignificant, but the BUSCO score confused me about which one is the best~/~more complete.

Any other thing that you suggest me to try?
According to the authors of Tiberius, it is better than Braker3.

Thanks;

Hi - glad you have been finding our tutorial helpful! Very cool that you tried Tiberius - that one has been on our to-do list but I was having trouble with installation. That's not a bad BUSCO score, especially for the lower number of genes.

The BUSCO score for your assembly sounds great. The score for the annotation has to be lower be definition, but 90% seems quite low. I'd like to confirm the evidence you used for EviAnn: protein databases from related species and transcript data from closely related species? Did you not say you also had RNA-seq from the species you sequenced? If you used protein and transcript data from closely related species and not the species you sequenced, that could partly explain the lower score. In my experience, protein data doesn't seem to create the most thorough annotations and RNA-seq that's not similar enough to your species might get ditched.

Have you tried Liftoff with a closely related species? We have great success getting high BUSCO scores from that and there are some good whale RefSeq genomes available. Liftoff is super easy to run, it rarely takes over an hour so highly recommend.

As for finding the longest isoform, I personally don't think it's super suprising that your score is lower. The longest isoform doesn't necessarily mean it's the "best" gene model to keep, and sometimes these annotators can squish multiple things together that don't belong therefore creating a very long gene model that won't match anything. We have been using Mikado to pick the "best" gene model, but another option that I haven't yet explored is OMAnnotator. Both of these tools essentially take multiple GFFs as input and create a consensus annotation. These will probably help improve your BUSCO score rather than relying on a single tool.

As for the total number of genes detected, it sounds like Tiberius is more precise than EviAnn. It's hard to say which annotation is "better" as it's hard to measure genome annotation quality, but the EviAnn annotation seems to be missing some key gene models, whereas Tiberius may have some false negatives in there, not sure.

If I were you, I would keep the EviAnn and Tiberius annotations, give Liftoff a go, and then try combining the annotations with Mikado or OMAnnotator and see what happens with your BUSCO scores and other quality metrics.

Hi
Thank you very much for the quick reply. I have been trying to figure out why the annotation part is so confusing, genome assembly was pretty straightforward. I did a bat assembly before, and I think that I got the opposite of what is expected - my genome assembly had a BUSCO C of 95% but the proteome BUSCO C score was 98. I used GeMoMa at the time for annotation, and I was super happy that BUSCO was higher from the proteome.

For the annotation, we have RNA-seq from the same species seven tissues (NCBI). On top of that, I used transcript, and protein data from related and "distant" species.
For this sister species, I will try liftoff, any parameter you recommend to try? I remember reading one from humans T2T, and they used liftoff with the -copy parameter + -sc 0.95 _ extra genes should be at least 95% similar. The sister species was annotated by a company in germany (Goenomics) - you might have heard off, using an old-fashioned way - augustus plus Stringie basically.
The species I am working, Weddel seal, had its genome sequenced - short-reads, very fragmented, and assembled using RNA-seq from 4 tissues, while 7 are now available on NCBI.
https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000349705.1/
The genome annotation on NCBI states that it has 22214 protein-coding genes, which makes it a little-tiny-bit less confusing, and makes the EviAnn annotation - look better, though the lower BUSCO C score.

Which other quality metrics should I use? I was basically using the gene number and the BUSCO, and somehow using the Omamer assessment.

We plan to do ortholog detection - using Orthofinder, and it needs a single isoform per gene - for all species, otherwise multiple isoforms will inflate the output not only during the clustering step in Orthofinder but also during later (PAML) when using dm/ds to detect selection with the multiple alignments output from Orthofinder, which we plan to do next.

If you want I could share my conda environment clone for you to run there. Installation of Tiberius is very complicated and ridiculously painful - though it is supposed to be simple. you need an old version of python - but for some reason, conda keeps on updating it... to analyze seqs > 200k , otherwise you have to chop your assembly, which sounds very silly and ridiculous. What I found incredible, is that they recommed using GPU, which might seem reasonable, but I don't have access to one; my first attempt using a 100 GB node took 35 days to run, but when I used a node with 500 GB of RAM, it took less than a day. Can't explain the reduction in runtime.

@14zac2

Hey!

I did a bat assembly before, and I think that I got the opposite of what is expected - my genome assembly had a BUSCO C of 95% but the proteome BUSCO C score was 98. -

That is indeed interesting. Did you happen to investigate the TE content and/or intron lengths of the species? Do you also happen to know what the assembly QV was for the species (i.e., predicted base-pair accuracy?) If there is weird splicing, TE interruption in genes and exons, super super long introns, or inaccurate base-pairs around splice sites BUSCO-assembly could assign the gene as fragmented while the busco-protein would see all of the exons, not care about anything else, and say it's "complete". Short read assemblies have huge issues deriving LINE elements which like to live in introns, and bats have accumulated LINE elements in recent evolution (Morales et al., 2025) (and LINEs are super common in mammal generally). If your assembly didn't have long reads, I'd guess that the bat annotation BUSCO was greater than the assembly BUSCO because of the inherent challenges in correctly assembling TEs with short read data.

Then, I used the new released AI-annotator/ Tiberius, which in terms of BUSCO got higher scores

My (basic) understanding is that Tiberius will do a good job with Busco genes. Tiberius is a deep learning classifier that is trained on 37 mammalian genome annotations. These genomes were selected based on their phylogenetic diversity and if they have an annotation BUSCO of >90%. They also picked the longest isoform for training (which is why it focuses on transcripts and intron-exon chains) and will not provide alternative isoforms. In a sense, TIBERIUS learns the features needed to find the single best BUSCO gene in a mammalian genome at the cost of some sensitivity for species specific stuff and having multiple transcripts per-gene. As Zoe said, we’re a big fan of integrating each GFF and asking an algorithm (we use Mikado in GenAnT) do decide which transcript is best at a per-transcript level.

As for Tiberius env. I’m okay, thank you! We’re trying to make GenAnT more container driven with the hopes of making a comprehensive annotation environment that anyone can install so I was going to try https://hub.docker.com/r/larsgabriel23/tiberius @14zac2 ?

liftoff with the -copy parameter + -sc 0.95 _ extra genes should be at least 95% similar.

That’s what we use. We’re starting to play with liftoff parameters and the relationship between species relatedness, genome assembly quality, and when to use liftoff vs. TOGA. Any of our findings at this stage are totally anecdotal though (but basically liftoff default does a remarkably good job for the same genus and can even do a good job for the same family). Once we try this more systematically, I’ll have a better idea here. There is also LiftON https://pubmed.ncbi.nlm.nih.gov/38798552/ which integrates liftoff and miniprot and was designed to help for slightly more distant species.

Which other quality metrics should I use? I was basically using the gene number and the BUSCO, and somehow using the Omamer assessment.

Good question, the written tutorial associated with this github is currently in press at Nat Protocols and should be out any day now and we try our best to answer this there. I’ll try and summarize below:
• BUSCO is good albeit a bit limited
• Total gene number can be useful when relating to a sister species. Like if you have >20% fewer or more genes than a sister species than it could we worth looking into. More doesn’t strictly mean better though as it can be a sign of gene fragmentation.
• GFFcompare can also be useful, especially because you have the NCBI annotation. For example, you can use liftoff to carry the NCBI annotation over with liftoff (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000349705.1/) and directly compare assembly precision and recall. Of course your annotation may be better than theirs as you have more RNA-seq but if either metric is terrible it’s a sign something went wrong. You can also use their annotation as a layer of evidence for your annotations with Mikado, which I really like to do when improving assemblies and annotations.
• I’m personally a huge fan of looking on the genome browser. Throw on some bam files and pick some genomic regions that are syntenic across mammals to make sure the gene order is good (TP53 is a good one). You can also pick some big genes like titin/timeless and see if they’re fragmented into smaller chunks compared to a high-quality sister species. Looking at assembly/annotation together it can be useful to check for some hard to annotate genes (Sox2, XIST, Sry, opposite-strand RNA’s). But spotting those (or not) can tell you a bit about if you
• Specific to if you use Mikado but I found the Mikado metrics file to be quite valuable. Every transcript has the values in this list associated with it: https://mikado.readthedocs.io/en/stable/Scoring_files/#metrics. https://github.com/BaderLab/GenAnT/blob/main/6-AfterAnnotation/PostprocessReformatFinalGFF.R in the tutorial talks about how to add it to a table of gene symbols + a GFF

We plan to do ortholog detection - using Orthofinder, and it needs a single isoform per gene - for all species, otherwise multiple isoforms will inflate the output not only during the clustering step in Orthofinder but also during later (PAML) when using dm/ds to detect selection with the multiple alignments output from Orthofinder, which we plan to do next.

I haven’t added it to GenAnT but this is the script I wrote to take the gene symbol table from GenAnT and filter it to help with Orthofinder/HyPy type analyses. Basically the idea is that it keeps genes with a unique symbol. Genes that have duplicates (ABC-copy1, ABC-copy2) are then checked if. If “copy2” is low quality or single exon, ABC-copy1 is kept and ABC-copy2 is removed. If both ABC-copy1 and ABC-copy2 are multi exonic and high quality, then they are both removed as they’re paralogs falsely given the same gene symbol or a gene copy event. Either way it wouldn’t be helpful for Orthofinder.

library(pbapply)
tab <- read.table("mHetGla1.1.pri_gene_symbols_metrics.tsv",header=TRUE,as.is=TRUE,sep="\t")
tab <- tab[!is.na(tab$hierarchical_gene),]

tab_fix_copy <- c()

uni <- unique(tab$hierarchical_gene)

single_copy_table <- function(i) {
  print(i)
  tab1 <- tab[tab$hierarchical_gene == i,]
  if(nrow(tab1) == 1) {
    return(tab1)
  } else {
    tab1 <- tab1[tab1$hierarchical_gene == i & tab1$is_complete == "True" & tab1$score > 16 & tab1$highest_cds_exon_number > 1,]
    if(nrow(tab1) == 1) {
      return(tab1)
    }
  }
}	
single_copy_genes <- as.data.frame(do.call("rbind",pblapply(uni,single_copy_table)))

write.table(single_copy_genes, file = "mHetGlaPri.single.copy.tsv",
            quote=FALSE,row.names = FALSE,col.names = TRUE,sep="\t")

Hi,

I wanted to ask you another question.
I was reading and I noticed that for non-coding RNA annotation, you mentioned only using Infernal. I have seen people using tRNAscan-SE for tRNA. Is there any reason why you did not use it?

Thanks

Hello! We've tried tRNAscan before, but largely recommend Infernal because it detects more kinds of ncRNA. In the end, it reduces the number of tools you have to use to try to find all kinds of different ncRNA.

Hello! We've tried tRNAscan before, but largely recommend Infernal because it detects more kinds of ncRNA. In the end, it reduces the number of tools you have to use to try to find all kinds of different ncRNA.

like Mirmachine (which we use), tRNAscan-SE uses Infernal1.1 as it's engine, the core difference is in the pre-trained tRNA models for tRNAscan-SE. When going through the paper my interpretation was that these models are super helpful for mt-tRNA in animals, but they didn't do much to discuss chromosomal tRNA. With enough interest I could test incorporating it into our workflows.

If you end up using tRNAscan-SE for tRNA specifically and Infernal for other ncRNA familes, I'd think about removing the tRNAs called by Infernal to avoid redundancy.

Hi,
Thanks for the suggestions and comments. Highly appreciate.
I hope you don't mind me asking a new question. I tried using gene model filtering with Mikado but the last step failed several times, and I feel that I should quit. The HPC where I work doesn't have much assistance, and I wanted to move on for non-coding annotation. I would appreciate your guidance on how to proceed.
I used EarlGrey to identify repeats, and I tried retrieving RNA from the gff, and found none. I tried Liftoff with the previous annotation, but for some reason the gff failed every time.
I am now creating the blast.bed after blasting the genome to the RFam db.
Thanks!

Hey, sorry you're having so much trouble with Mikado and the other tools! That sounds super frustrating. We're happy to help with these issues.

Are you able to provide more details about why these different steps might be failing?

For Mikado: Did the first three steps run correctly? And is there an error message that comes up in the "Mikado pick" log file?
For EarlGrey: Did it complete a successful run and spit out a GFF file at all? If so, it would be super helpful to see a snippit of the GFF.
For Liftoff: Liftoff is usually fairly reliable so I'm sure we can make it work. What GFF are you using as the old annotation? Is it from RefSeq or Ensembl? Is Liftoff giving you any specific error messages?

Thanks!
Zoe

@14zac2

Hi @14zac2
For Mikado, there is a problem with the bed file and Mikado serialize failed because of some model errors and some old Linux libraries.

I used this code for EarlGrey:

earlGrey -g LW3.0.all.scaff.fasta \
        -s Lepwed  -o LepWed-3_earlgrey -t 100 -d yes -e yes -m yes

I didn't have any issues with it. I was just surprised that it didn't have any RNA.
Here are a few lines;

Lepwed.3_1      0       2330    Simple_repeat   2135    +
Lepwed.3_1      2656    3000    LINE/L1 543     -
Lepwed.3_1      3006    4891    LINE/L1 1143    +
Lepwed.3_1      7666    7879    LINE/L1 1394    +
Lepwed.3_1      8606    8797    LTR/ERVL-MaLR   432     -
Lepwed.3_1      10720   12191   Simple_repeat   507     +
Lepwed.3_1      14166   14276   Simple_repeat   30      +
Lepwed.3_1      15031   15373   LINE/L1 611     +
Lepwed.3_1      15373   15990   LINE/L1 821     -
Lepwed.3_1      16325   16480   LINE/L1 438     -
Lepwed.3_1      16507   16838   LINE/L1 1592    -
Lepwed.3_1      16938   17150   LINE/L1 741     -
Lepwed.3_1      19630   19881   LINE/L1 391     +
Lepwed.3_1      20929   21169   Simple_repeat   92      +
Lepwed.3_1      23659   23909   LINE/L1 322     +

For Liftoff I basically downloaded the NCBI GFF of my species, GCA_000349705.1, Leptonychotes weddellii (Weddell seal). I was able to transfer the annotation that I did from the sister species, which doesn't have non-coding annotations, and I did understand that Liftoff only transfer protein-coding genes.
The error log is this:

2025-07-15 22:27:32,340 - INFO - Populating features
Populating features table and first-order relations: 0 features
Populating features table and first-order relations: 1000 features
Populating features table and first-order relations: 2000 features
Populating features table and first-order relations: 3000 features
Populating features table and first-order relations: 4000 features
Populating features table and first-order relations: 5000 features....
Populating features table and first-order relations: 1236000 features
Populating features table and first-order relations: 1237000 features
Populating features table and first-order relations: 1238000 features
2025-07-15 22:39:20,475 - INFO - Populating features table and first-order relations: 1238453 features
2025-07-15 22:39:20,476 - INFO - Updating relations
2025-07-15 22:41:16,511 - INFO - Creating relations(parent) index
2025-07-15 22:41:21,526 - INFO - Creating relations(child) index
2025-07-15 22:41:28,042 - INFO - Creating features(featuretype) index
2025-07-15 22:41:31,105 - INFO - Creating features (seqid, start, end) index
2025-07-15 22:41:35,767 - INFO - Creating features (seqid, start, end, strand) index
2025-07-15 22:41:40,908 - INFO - Running ANALYZE features
[M::main::23.3211.00] loaded/built the index for 203 target sequence(s)
[M::mm_mapopt_update::27.7541.00] mid_occ = 471
[M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 203
[M::mm_idx_stat::30.2671.00] distinct minimizers: 99587303 (41.38% are singletons); average occurrences: 4.626; average spacing: 5.328; total length: 2454374544
[M::worker_pipeline::122.76022.08] mapped 10336 sequences
[M::worker_pipeline::198.93425.36] mapped 12698 sequences
[M::worker_pipeline::208.82025.22] mapped 2402 sequences
[M::main] Version: 2.30-r1287
[M::main] CMD: minimap2 -o intermediate_files/reference_all_to_target_all.sam -a --end-bonus 5 --eqx -N 50 -p 0.5 -t 32 /data/common/juanpablo.aguilar/data/compleasm/HydLep.f1x.fasta.mmi intermediate_files/reference_all_genes.fa
[M::main] Real time: 209.811 sec; CPU: 5267.905 sec; Peak RSS: 16.628 GB
Traceback (most recent call last):
File "/data/common/juanpablo.aguilar/apps/miniconda3/envs/NGS-basic/bin/liftoff", line 10, in
sys.exit(main())
^^^^^^
File "/data/common/juanpablo.aguilar/apps/miniconda3/envs/NGS-basic/lib/python3.12/site-packages/liftoff/run_liftoff.py", line 12, in main
run_all_liftoff_steps(args)
File "/data/common/juanpablo.aguilar/apps/miniconda3/envs/NGS-basic/lib/python3.12/site-packages/liftoff/run_liftoff.py", line 24, in run_all_liftoff_steps
feature_db, feature_hierarchy, ref_parent_order = liftover_types.lift_original_annotation(ref_chroms, target_chroms,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/data/common/juanpablo.aguilar/apps/miniconda3/envs/NGS-basic/lib/python3.12/site-packages/liftoff/liftover_types.py", line 15, in lift_original_annotation
align_and_lift_features(ref_chroms, target_chroms, args, feature_hierarchy, liftover_type, unmapped_features,
File "/data/common/juanpablo.aguilar/apps/miniconda3/envs/NGS-basic/lib/python3.12/site-packages/liftoff/liftover_types.py", line 23, in align_and_lift_features
aligned_segments= align_features.align_features_to_target(ref_chroms, target_chroms, args,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/data/common/juanpablo.aguilar/apps/miniconda3/envs/NGS-basic/lib/python3.12/site-packages/liftoff/align_features.py", line 28, in align_features_to_target
return parse_all_sam_files(feature_hierarchy, unmapped_features, liftover_type, sam_files)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/data/common/juanpablo.aguilar/apps/miniconda3/envs/NGS-basic/lib/python3.12/site-packages/liftoff/align_features.py", line 118, in parse_all_sam_files
aligned_segments = parse_alignment(file, feature_hierarchy, unmapped_features, liftover_type)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/data/common/juanpablo.aguilar/apps/miniconda3/envs/NGS-basic/lib/python3.12/site-packages/liftoff/align_features.py", line 132, in parse_alignment
aln_id = add_alignment(ref_seq, align_count_dict, search_type, name_dict,aln_id, feature_hierarchy,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/data/common/juanpablo.aguilar/apps/miniconda3/envs/NGS-basic/lib/python3.12/site-packages/liftoff/align_features.py", line 149, in add_alignment
aligned_blocks = get_aligned_blocks(ref_seq, aln_id, feature_hierarchy, search_type)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/data/common/juanpablo.aguilar/apps/miniconda3/envs/NGS-basic/lib/python3.12/site-packages/liftoff/align_features.py", line 170, in get_aligned_blocks
parent = feature_hierarchy.parents[liftoff_utils.convert_id_to_original(alignment.query_name)]
~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
KeyError: 'gene-12S'

I tried filtering it, several ways, with AGAT, with gffread, nothing worked.

I would like just to get the non-coding genes. I used Tiberius for ab initio prediction, followed by EviAnn with RNA-seq, protein, and transcript evidence, and found 22k genes. What confused me is that using full proteome the BUSCO complete score (from carnivora db 12) was high - 95.8%, but since it has all the proteins from isoforms, the duplicates were ridiculously high: 69.9%. So, I used AGAT - the only way I know, to get the longest isoform per gene, and the duplicate scores reduced nicely, 1.0%, but so did the complete 93.9%. I feel that I am trying too hard, since I am working with a marine mammal, and BUSCO is based on "limited species coverage".

My next goal is to be able to get full length mRNA sequencing, but it might take a while to do, since PacBio is very expensive.

Thanks!

@risserlin

Hey!

Sorry to answer your questions with more questions but:

Did you install mikado on your own or following the setup scripts here? @risserlin built a new mikado container for GenAnT which I strongly recommend.

singularity build mikado_gat.sif docker://risserlin/mikado:ubuntu22_mikado2.3.2

I don't think I'd be too worried about it. the "RNAs" in EarlGrey are TE derived tRNA snRNA snoRNA etc. I'm not sure what repeat library your command would default to but it may have redundancies with other repeat elements. For example the tRNA containing element may be called another SINEs and LTR retrotransposons, or maybe your species doesn't have many evolutionarily "young" TE derived ncRNAs.
Same question as (1). Did you use our helper script for processing your reference genome? It tries to smooth out a lot of the weirdness that can come from GFF formats not being super strict (even if they're all in NCBI).

https://github.com/BaderLab/GenAnT/blob/main/setup/reference_directory_refseq.sh

3a) As for liftoff transferring non-coding models, the answer is slightly annoying. liftoff, by default, transfers features you'd expect to see in PCGs (gene, mRNA / transcript, exon, CDS, five_prime_UTR, three_prime_UTR, start_codon, stop_codon) and I believe it reads the "type" column of the gff to filter what to transfer. You can expand this to include other types with the -f parameter but personally I've found it pretty buggy. The easiest workaround is to apply liftoff to the reference gtf file instead of the gff from NCBI since it uses a much stricter type formatting (i.e., all unique genes are given the "gene" type regardless of if they're PcGs). You can also just use an ensembl annotation but I don't know if that's available for your species. If you go the gtf route, liftoff will return a gtf so don't forget to convert it into a gff with gffread or AGAT if you want to put it into mikado.

Hi!

Hey!

Sorry to answer your questions with more questions but:

Did you install mikado on your own or following the setup scripts here? @risserlin built a new mikado container for GenAnT which I strongly recommend.
singularity build mikado_gat.sif docker://risserlin/mikado:ubuntu22_mikado2.3.2

I wouldn't call it on my own haha. I followed the instructions on the mikado website, since the authors might have changed something, so to at least us the most up-to-date version - despite it using Transdecoder output and not TD2.

I don't think I'd be too worried about it. the "RNAs" in EarlGrey are TE derived tRNA snRNA snoRNA etc. I'm not sure what repeat library your command would default to but it may have redundancies with other repeat elements. For example the tRNA containing element may be called another SINEs and LTR retrotransposons, or maybe your species doesn't have many evolutionarily "young" TE derived ncRNAs.

That is how I feel. I would not based noncoding annotation from an app that is for transposable elements.

Same question as (1). Did you use our helper script for processing your reference genome? It tries to smooth out a lot of the weirdness that can come from GFF formats not being super strict (even if they're all in NCBI).
https://github.com/BaderLab/GenAnT/blob/main/setup/reference_directory_refseq.sh

Is the script to process the reference genome or the annotation GFF file?

3a) As for liftoff transferring non-coding models, the answer is slightly annoying. liftoff, by default, transfers features you'd expect to see in PCGs (gene, mRNA / transcript, exon, CDS, five_prime_UTR, three_prime_UTR, start_codon, stop_codon) and I believe it reads the "type" column of the gff to filter what to transfer. You can expand this to include other types with the -f parameter but personally I've found it pretty buggy. The easiest workaround is to apply liftoff to the reference gtf file instead of the gff from NCBI since it uses a much stricter type formatting (i.e., all unique genes are given the "gene" type regardless of if they're PcGs). You can also just use an ensembl annotation but I don't know if that's available for your species. If you go the gtf route, liftoff will return a gtf so don't forget to convert it into a gff with gffread or AGAT if you want to put it into mikado.

Any program do you trust or recommend for making a gff a gtf?
Indeed, Ensembl doesn't have my species.
I was mostly curious to see if Mikado does some magic and the BUSCO proteome score increases - the BUSCO genome complete was 98.1%, single 96.9%.

Hey!

RE: Mikado

Okay great. Using the linked container made a massive difference for us. I'd definitely install with that and give it another try if you're interested in mikado.

Is the script to process the reference genome or the annotation GFF file?

It's for GFFs but has a couple of basic cleaning steps for a reference fa that we've added to deal with edge cases (e.g., cleaning fasta headers).

Any program do you trust or recommend for making a gff a gtf?

Like you I think that agat_convert_sp_gxf2gxf.pl is the way to go. gffread can drop ID features.

I was mostly curious to see if Mikado does some magic and the BUSCO proteome score increases - the BUSCO genome complete was 98.1%, single 96.9%.

That would be cool! Would that mainly integrating liftoff and EviAnn? I guess once liftoff works you can guess if it does by running gffcompare between the annotations to get an idea of what may be missing from the other. If there is a reasonable amount then mikado may be able to fill in the gap.

Best!
Dustin

Hey!

RE: Mikado

Okay great. Using the linked container made a massive difference for us. I'd definitely install with that and give it another try if you're interested in mikado.

The biggest problem that I see is that the bed file had some issues that I am trying to figure out with the authors of mikado; 1, issue is running it - I did an installation that worked the configure and prepare steps, but I rebuilt the environment and now I can't install it -following the original instructions...
Kinda exhausted with all the problems. I have the singularity container, let me know any parameter that I should use to run mikado.

Is the script to process the reference genome or the annotation GFF file?

It's for GFFs but has a couple of basic cleaning steps for a reference fa that we've added to deal with edge cases (e.g., cleaning fasta headers).

Any program do you trust or recommend for making a gff a gtf?

Like you I think that agat_convert_sp_gxf2gxf.pl is the way to go. gffread can drop ID features.

I was mostly curious to see if Mikado does some magic and the BUSCO proteome score increases - the BUSCO genome complete was 98.1%, single 96.9%.

That would be cool! Would that mainly integrating liftoff and EviAnn? I guess once liftoff works you can guess if it does by running gffcompare between the annotations to get an idea of what may be missing from the other. If there is a reasonable amount then mikado may be able to fill in the gap.

Best! Dustin

I would use the annotation from Tiberius and with EviAnn. I am not that worried with the Liftoff because the BUSCO score of the publicly available genome is 89%, while ours is 98%, so I don't think it will be very significant impact.
For example, my biggest question is what is affected after doing AGAT longest isoform extraction that reduces the BUSCO score.
Any suggestions appreciated.

Kinda exhausted with all the problems. I have the singularity container, let me know any parameter that I should use to run mikado.

The best place to see how we deal with mikado would be here: https://github.com/BaderLab/GenAnT/tree/main/3-CombineAndFilter
All of the scripts we reference are in https://github.com/BaderLab/GenAnT/tree/main/scripts

I would use the annotation from Tiberius and with EviAnn. I am not that worried with the Liftoff because the BUSCO score of the publicly available genome is 89%, while ours is 98%, so I don't think it will be very significant impact.

"This is with full proteome, I used agat to get the longest isoform per gene, and the BUSCO was lower"

I assume you are referencing this?
The BUSCO proteome is lower:
C:90.8%[S:25.9%,D:65.0%],F:2.1%,M:7.1%,n:13727
This is with full proteome, I used agat to get the longest isoform per gene, and the BUSCO was lower:
C:89.1%[S:88.4%,D:0.8%],F:2.7%,M:8.2%,n:13727
This is something I think mikado actually can help with. The answer is that longest isoform does not necessarily mean top isoform due to technical issues stuff like transcript chimeras, utr extensions or real biological events like rare exons from RNA-seq data etc. You can ask this with the "metrics file" mikado returned.

For example, for the first 10K genes in my naked mole-rat assembly, I asked if the longest isoform for the gene has the top aggregate score from their huge metrics list. The longest isoform was the top transcript 90.21% of the time. I'd guess your drop in busco score is that there are a handful of long isoforms that have some issue, making them not pass BUSCO, however you have a different isoform in the same gene that is correct, letting the aggregate busco proteome have a higher total score.

Best!
Dustin

@DustinSokolowski

@DustinSokolowski is the master of answering questions BUT I just wanted to chip in to add that if you're using the RefSeq Weddell Seal as a reference, you can get both the GFF and GTF from ftp: https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/349/705/GCF_000349705.1_LepWed1.0/GCF_000349705.1_LepWed1.0_genomic.gff.gz

Hello hello!

I'm glad that you're finding a good combination of tools that is working for you. It's possible that you have a higher BUSCO score and more contamination genes because that particular gene set contains more false positive gene models. Is there also a larger number of genes/features in the one with the more contamination genes?

It's important to remember that BUSCO is a rather limited quality assessment tool. Although it is great for assessing highly-conserved regions, it is unclear how good the rest of the genome/annotation is. Any tool that doesn't use RNA-seq directly is probably going to perform best on conserved regions whereas RNA-seq wouldn't have that bias.

Good question about the gene symbols! I think it's a matter of personal preference. I prefer OrthoFinder because if you use human gene symbols, you will get more gene symbols overall that are easier to read. If you use LiftOff, you are limited to the gene symbols that came from your closely-related species (and some of these gene symbols can have weird names).

Cheers,
Zoe

Hi!
Yes, the Mikado gene set have more features/genes, the contamination level was ~20% which was surprisingly high.
Unfortunately, I work with marine mammals and when the project was planned it only included standard RNA-seq. I am now convincing my boss to do long RNA-seq to have better gene model prediction. I used EviAnn and it uses Stringtie and recommends doing tissue-by-tissue gene prediction to minimize mixing transcripts from different tissues and creating chimeric transcripts. Furthermore, for my species, the project only collected skin tissue, which limits the accuracy and power - though a previous project sequenced RNA from 7 other tissues - including placenta!. So, there is some "hard evidence" for finding genes.
Hopefully, we can find funds to do PacBio Kinnex to confirm and find new transcripts.

Funnily, the human T2T assembly reported 20k genes and my species has 24k haha. Science is cool!

Thanks!

Hey!

I also like tissue-by-tissue stringTie and is what GenAnT uses in that step by recommending you merge each bam file by tissue before running the tool. Braker also internally runs stringtie on each individual inputted bam file and then merges afterwards. As you know, using publicly available RNA-seq is definitely encouraged when performing gene annotations so that's great!

Hi.
I appreciate all the help during the annotation process, and I hope you don't mind me your names in the acknowledgment section of the paper. Please let me know how should I address you, your preferred name.
Thank you very much!

Hey!

That's very generous, thank you!! My preferred name is Dustin J Sokolowski. I think Zoe would want Zoe A. Clarke.

Thanks again!
Dustin

It's been our pleasure! Thanks for the acknowledgement! If you found our general advice helpful, citing our paper would be awesome, too, if you can :)

Sure! It is! Thanks for the advice, help and hints, Dr. Dustin J Sokolowski and Dr. Zoe A. Clarke! Juan Pablo Aguilar

…

________________________________ From: Zoe Clarke ***@***.***> Sent: Wednesday, April 22, 2026 9:29:57 PM To: BaderLab/GenAnT ***@***.***> Cc: Aguilar Cabezas, Juan Pablo ***@***.***>; Mention ***@***.***> Subject: [External] Re: [BaderLab/GenAnT] Question about genome annotation QC (Issue #1) Use caution with links and attachments. [https://github.com/proxy/avatars.githubusercontent.com/u/41270914?s=20&v=4]14zac2 left a comment (BaderLab/GenAnT#1)<#1 (comment)> It's been our pleasure! Thanks for the acknowledgement! If you found our general advice helpful, citing our paper would be awesome, too, if you can :) — Reply to this email directly, view it on GitHub<#1?email_source=notifications&email_token=AJWD2VNG57XKSRYYKIALEQ34XFWZLA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTIMZQGEYDSMZVHE3KM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#issuecomment-4301093596>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AJWD2VLWP3V2N2SC25HG6QL4XFWZLAVCNFSM6AAAAACRDT7GY6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHM2DGMBRGA4TGNJZGY>. You are receiving this because you were mentioned.Message ID: ***@***.***>

Hi,
Sorry to bother you one more time. I annotated the genome using all the PCG but later I found that some of them might be erroneous models (ID'd by OMark) so I felt that just deleted them after using them during annotation might be inappropriate.
So, I deleted those, but I also forgot that EviAnn does annotate lncRNA, and it found 1702 . I did a quick look at the EviAnn annotation file, and I noticed something interesting and complicated/confusing at the same time; EviAnn outputted for lncRNAs a "gene" and also like "mRNAs" and "Exons" - based on RNA-seq data:

lncRNA with 2 "mRNA"s

Lepwed.3_1      EviAnn  gene    53764485        53778317        .       -       .       ID=XLOC_001061U_lncRNA;geneID=XLOC_001061U_lncRNA;gene_biotype=lncRNA
Lepwed.3_1      EviAnn  lnc_RNA 53764485        53773740        .       -       .       ID=XLOC_001061U_lncRNA-mRNA-1;Parent=XLOC_001061U_lncRNA;EvidenceTranscriptID=MSTRG_00021992:2:5.518105;geneID=XLOC_001061U_lncRNA
Lepwed.3_1      EviAnn  exon    53764485        53764827        .       -       .       Parent=XLOC_001061U_lncRNA-mRNA-1;ID=XLOC_001061U_lncRNA-mRNA-1:exon:1
Lepwed.3_1      EviAnn  exon    53773338        53773740        .       -       .       Parent=XLOC_001061U_lncRNA-mRNA-1;ID=XLOC_001061U_lncRNA-mRNA-1:exon:2
Lepwed.3_1      EviAnn  lnc_RNA 53764537        53778317        .       -       .       ID=XLOC_001061U_lncRNA-mRNA-2;Parent=XLOC_001061U_lncRNA;EvidenceTranscriptID=MSTRG_00003899:2:6.192656;geneID=XLOC_001061U_lncRNA
Lepwed.3_1      EviAnn  exon    53764537        53764827        .       -       .       Parent=XLOC_001061U_lncRNA-mRNA-2;ID=XLOC_001061U_lncRNA-mRNA-2:exon:1
Lepwed.3_1      EviAnn  exon    53773338        53773559        .       -       .       Parent=XLOC_001061U_lncRNA-mRNA-2;ID=XLOC_001061U_lncRNA-mRNA-2:exon:2
Lepwed.3_1      EviAnn  exon    53773656        53773754        .       -       .       Parent=XLOC_001061U_lncRNA-mRNA-2;ID=XLOC_001061U_lncRNA-mRNA-2:exon:3
Lepwed.3_1      EviAnn  exon    53777273        53778317        .       -       .       Parent=XLOC_001061U_lncRNA-mRNA-2;ID=XLOC_001061U_lncRNA-mRNA-2:exon:4

Lepwed.3_1      EviAnn  gene    172194737       172211863       .       +       .       ID=XLOC_000691U_lncRNA;geneID=XLOC_000691U_lncRNA;gene_biotype=lncRNA
Lepwed.3_1      EviAnn  lnc_RNA 172194737       172211863       .       +       .       ID=XLOC_000691U_lncRNA-mRNA-1;Parent=XLOC_000691U_lncRNA;EvidenceTranscriptID=MSTRG_00021368:3:21.453440;geneID=XLOC_000691U_lncRNA
Lepwed.3_1      EviAnn  exon    172194737       172194833       .       +       .       Parent=XLOC_000691U_lncRNA-mRNA-1;ID=XLOC_000691U_lncRNA-mRNA-1:exon:1
Lepwed.3_1      EviAnn  exon    172194918       172195020       .       +       .       Parent=XLOC_000691U_lncRNA-mRNA-1;ID=XLOC_000691U_lncRNA-mRNA-1:exon:2
Lepwed.3_1      EviAnn  exon    172198370       172198461       .       +       .       Parent=XLOC_000691U_lncRNA-mRNA-1;ID=XLOC_000691U_lncRNA-mRNA-1:exon:3
Lepwed.3_1      EviAnn  exon    172211001       172211863       .       +       .       Parent=XLOC_000691U_lncRNA-mRNA-1;ID=XLOC_000691U_lncRNA-mRNA-1:exon:4
Lepwed.3_1      EviAnn  lnc_RNA 172194751       172202244       .       +       .       ID=XLOC_000691U_lncRNA-mRNA-2;Parent=XLOC_000691U_lncRNA;EvidenceTranscriptID=MSTRG_00009089:6:9.305950;geneID=XLOC_000691U_lncRNA
Lepwed.3_1      EviAnn  exon    172194751       172194833       .       +       .       Parent=XLOC_000691U_lncRNA-mRNA-2;ID=XLOC_000691U_lncRNA-mRNA-2:exon:1
Lepwed.3_1      EviAnn  exon    172194918       172195020       .       +       .       Parent=XLOC_000691U_lncRNA-mRNA-2;ID=XLOC_000691U_lncRNA-mRNA-2:exon:2
Lepwed.3_1      EviAnn  exon    172198370       172198461       .       +       .       Parent=XLOC_000691U_lncRNA-mRNA-2;ID=XLOC_000691U_lncRNA-mRNA-2:exon:3
Lepwed.3_1      EviAnn  exon    172201572       172202244       .       +       .       Parent=XLOC_000691U_lncRNA-mRNA-2;ID=XLOC_000691U_lncRNA-mRNA-2:exon:4

lncRNA with 3 "mRNA"s

grep "XLOC_018268U" LW2_lncRNA.gff
Lepwed.3_8      EviAnn  gene    74515540        74517499        .       -       .       ID=XLOC_018268U_lncRNA;geneID=XLOC_018268U_lncRNA;gene_biotype=lncRNA
Lepwed.3_8      EviAnn  lnc_RNA 74515540        74517497        .       -       .       ID=XLOC_018268U_lncRNA-mRNA-1;Parent=XLOC_018268U_lncRNA;EvidenceTranscriptID=MSTRG_00323542:8:209.581955;geneID=XLOC_018268U_lncRNA
Lepwed.3_8      EviAnn  exon    74515540        74515938        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-1;ID=XLOC_018268U_lncRNA-mRNA-1:exon:1
Lepwed.3_8      EviAnn  exon    74516170        74516219        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-1;ID=XLOC_018268U_lncRNA-mRNA-1:exon:2
Lepwed.3_8      EviAnn  exon    74516529        74516590        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-1;ID=XLOC_018268U_lncRNA-mRNA-1:exon:3
Lepwed.3_8      EviAnn  exon    74516674        74516737        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-1;ID=XLOC_018268U_lncRNA-mRNA-1:exon:4
Lepwed.3_8      EviAnn  exon    74516822        74517046        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-1;ID=XLOC_018268U_lncRNA-mRNA-1:exon:5
Lepwed.3_8      EviAnn  exon    74517430        74517497        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-1;ID=XLOC_018268U_lncRNA-mRNA-1:exon:6
Lepwed.3_8      EviAnn  lnc_RNA 74515540        74517497        .       -       .       ID=XLOC_018268U_lncRNA-mRNA-2;Parent=XLOC_018268U_lncRNA;EvidenceTranscriptID=MSTRG_00323544:7:102.003525;geneID=XLOC_018268U_lncRNA
Lepwed.3_8      EviAnn  exon    74515540        74515938        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-2;ID=XLOC_018268U_lncRNA-mRNA-2:exon:1
Lepwed.3_8      EviAnn  exon    74516155        74516219        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-2;ID=XLOC_018268U_lncRNA-mRNA-2:exon:2
Lepwed.3_8      EviAnn  exon    74516529        74516590        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-2;ID=XLOC_018268U_lncRNA-mRNA-2:exon:3
Lepwed.3_8      EviAnn  exon    74516674        74516737        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-2;ID=XLOC_018268U_lncRNA-mRNA-2:exon:4
Lepwed.3_8      EviAnn  exon    74516822        74517046        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-2;ID=XLOC_018268U_lncRNA-mRNA-2:exon:5
Lepwed.3_8      EviAnn  exon    74517430        74517497        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-2;ID=XLOC_018268U_lncRNA-mRNA-2:exon:6
Lepwed.3_8      EviAnn  lnc_RNA 74515543        74517499        .       -       .       ID=XLOC_018268U_lncRNA-mRNA-3;Parent=XLOC_018268U_lncRNA;EvidenceTranscriptID=MSTRG_00323543:7:63.662434;geneID=XLOC_018268U_lncRNA
Lepwed.3_8      EviAnn  exon    74515543        74515938        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-3;ID=XLOC_018268U_lncRNA-mRNA-3:exon:1
Lepwed.3_8      EviAnn  exon    74516529        74516590        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-3;ID=XLOC_018268U_lncRNA-mRNA-3:exon:2
Lepwed.3_8      EviAnn  exon    74516674        74516737        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-3;ID=XLOC_018268U_lncRNA-mRNA-3:exon:3
Lepwed.3_8      EviAnn  exon    74516822        74517046        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-3;ID=XLOC_018268U_lncRNA-mRNA-3:exon:4
Lepwed.3_8      EviAnn  exon    74517430        74517499        .       -       .       Parent=XLOC_018268U_lncRNA-mRNA-3;ID=XLOC_018268U_lncRNA-mRNA-3:exon:5

This information is interesting but now I am unsure what and how to "include" it during the annotation, like, should I only use the "gene" START and STOP positions, or should I use the "mRNA"s START/STOP positions? The "exon" regions are interesting, but I don't know if it makes sense to add it or not.

I would be extremely grateful for your input on what should be relevant and "valid"/ useful to use during the ncRNA annotation, and how to proceed.

I feel that I am so close but the end point seems like too far away haha.

THANK YOU VERY MUCH!

Hello, hello!

Genome annotation will never end if you don't force it to ;)

It looks like you have some alternatively spliced lncRNAs. I'm not super familiar with EviAnn - is the issue just that mRNA is appearing in the ID name?

I know that feeling. We used regular 2x150 bp RNA-seq and I am REALLY interested in getting long-read transcript data to improve the annotation!

The issue is that I am not sure what to add to the noncoding RNA-seed haha. I just GOoGle'd quickly and I found a paper actually describing lncRNA alternative splicing. Learning a new thing every day! Great!

I am in doubt whether to use the START-STOP positions of the parent lncRNA "gene" to make the combined bed seed file, or to use the START-STOP positions of the "mRNA" forms of the lncRNA "genes".
I am more leaning towards using the "gene" positions since the "mRNA" positions will be overlapping, but I would like to hear from the experts on the best way to approach this dilemma!

Thanks

Ah, I see! I was confused because you are referring to the lncRNAs as mRNAs. Yes, I would just use the gene coordinates! It should give you the loweset low and the highest high of each position.

Hi.
Sorry to bother you again.
I am trying to assign Gene Symbols, and for some reason the gene names from Orthofinder don't match the ones from the reference.table. I downloaded the Human T2T assembly and gff. Since Orthofinder needs one protein per gene, I used AGAT to make a gff with the longest isoform per gene, and then I used gffread to get the aa sequences:
gffread -g ../T2T.fasta -y Homsap.pep ../longert/T2T_longest_iso.gff

Maybe I am wrong, but I wanted to get as many genes as possible, so since Orthofinder complains when it finds an in-frame stop codons "" , I changed the "" on the aa seqs to X. I did the same with my species aa sequences, and Orthofinder ran fine.

I then created the reference.table as stated here using the total T2T annotation,

and when it tried to join the reference and the orthogroup_pairs tables, the gene Id column is empty:

> head(reference)
    gene orthofinder_gene
1 SCNN1D  rna-NR_037668.3
2  PUSL1  rna-NR_144369.2
3  CCNL2  rna-NR_146722.2
4  CCNL2  rna-NR_146723.2
5  CCNL2  rna-NR_135154.2
6   MIB2  rna-NR_033183.2

> orthofinder_pairs <- dplyr::left_join(orthofinder_pairs, reference, by = "orthofinder_gene")
> head(orthofinder_pairs)
            mikado_id orthofinder_gene gene
1 LOC_00007244-mRNA-1 rna-XM_054371163 <NA>
2 LOC_00022770-mRNA-1 rna-XM_054371163 <NA>
3 LOC_00005837-mRNA-1 rna-XM_054371163 <NA>
4 LOC_00021813-mRNA-1 rna-XM_054371163 <NA>
5 LOC_00015061-mRNA-1 rna-XM_054371163 <NA>
6 LOC_00014576-mRNA-1 rna-XM_054371163 <NA>

I am shocked and confused because I haven't changed anything in the gff. and for some reason the gene names don't match.
I even trimmed the reference gene names, to eliminate everything after the dot, and nothing.

head(reference)
    gene orthofinder_gene
1 SCNN1D  rna-NR_037668.3
2  PUSL1  rna-NR_144369.2
3  CCNL2  rna-NR_146722.2
4  CCNL2  rna-NR_146723.2
5  CCNL2  rna-NR_135154.2
6   MIB2  rna-NR_033183.2
reference$orthofinder_gene <- sub("\\..*", "", reference$orthofinder_gene)
head(reference)
    gene orthofinder_gene
1 SCNN1D    rna-NR_037668
2  PUSL1    rna-NR_144369
3  CCNL2    rna-NR_146722
4  CCNL2    rna-NR_146723
5  CCNL2    rna-NR_135154
6   MIB2    rna-NR_033183

sum(unique(orthofinder_pairs$orthofinder_gene) %in% reference$orthofinder_gene)

Perhaps could you bring some light on what might be going on?
I followed the code word-by-word - to make the table.txt,
as required in the R script:

## Gene name conversion

Read in the text file of gene symbols and transcript IDs to use as a conversion key between the IDs used for OrthoFinder and the gene symbols we would like to eventually be using.

reference <- read.table("reference.table.txt", header = TRUE, sep = "\t")

# Read in the GFF file using rtracklayer
gff <- readGFF("T2T.gff")

gff_transcript <- gff_transcript <- gff[gff$type == "transcript",]

# Get gene ID and transcript ID keys
# if this is a character already it doesn't change anything, but Rtracklayer sometimes loads "Parent" in as a list.

gff_transcript$gene_id <- unlist(gff_transcript$Parent)
gff_transcript$trans_id <- unlist(gff_transcript$ID)

key <- gff_transcript[,c("gene_id","trans_id","gene")]
colnames(key) <- c("geneID","transcriptID","geneSymbol")

write.table(key[,c("geneSymbol","transcriptID")],file="T2T.table.txt",
            quote=FALSE,row.names=FALSE,col.names=TRUE,sep="\t")

Hey! Sorry for the slow reply.

So while I don't think this will solve your issue, I'd be quite concerned with replacing stop codons with "unknown" codons. I would trust an internal stop codon as a real stop and that the rest of the transcript could be at risk of being misannotated. For example, if your reference has two transcripts incorrectly fused together, orthofinder would build false orthogenes where multiple genes in your target species have a 50% reference (human) 100% target match to your reference species, which could be enough to have two different genes get the same symbol. Since human is your reference you can skip any
Of these issues by using MANE: https://www.ncbi.nlm.nih.gov/refseq/MANE/

For your primary issue, are you making the gene symbol table from the main reference or from the agat-transformed reference gff? I'm wondering if agat may be changing what lives in ID/Parent etc as part of in cleaning process, and as a result you're now trying to line up protein ids with gene ids or something like that.

Hi, I did the Orthofinder analysis twice, with only the proteins with START/STOP CODONS and without in-frame STOP CODONS, and all of them.
I did make the gene symbol table from the main reference:

# Read in the GFF file using rtracklayer
gff <- readGFF("T2T.gff")

I just renamed the NCBI file to T2T.gff.
I checked the T2T.longest.gff and it seems fine to me. For instance I checked the protein:

 head Homsap.pep

rna-NM_001005277.1 ```
and the longest.gff it has::

NC_060925.1	BestRefSeq	gene	111940	112877	.	-	.	**ID=gene-OR4F16**;Dbxref=GeneID:81399,HGNC:HGNC:15079;Name=OR4F16;description=olfactory receptor family 4 subfamily F member 16;gbkey=Gene;gene=OR4F16;gene_biotype=protein_coding;gene_synonym=OR1-1,OR7-21
NC_060925.1	BestRefSeq	mRNA	111940	112877	.	-	.	ID=**rna-NM_001005277.1**;**Parent=gene-OR4F16;**Dbxref=GeneID:81399,GenBank:NM_001005277.1,HGNC:HGNC:15079;Name=NM_001005277.1;Note=The RefSeq transcript has 9 substitutions%2C 1 frameshift compared to this genomic sequence;exception=annotated by transcript or proteomic data;gbkey=mRNA;gene=OR4F16;inference=similar to RNA sequence%2C mRNA (same species):RefSeq:NM_001005277.1;product=olfactory receptor family 4 subfamily F member 16;tag=RefSeq Select;transcript_id=NM_001005277.1

I am surprised that the reference table and the Orthofinder results don't match:

head(reference)
gene orthofinder_gene
1 SCNN1D rna-NR_037668
2 PUSL1 rna-NR_144369
3 CCNL2 rna-NR_146722
4 CCNL2 rna-NR_146723
5 CCNL2 rna-NR_135154
6 MIB2 rna-NR_033183
head(orthofinder_pairs)
mikado_id orthofinder_gene
1 LOC_00007244-mRNA-1 rna-XM_054371163
2 LOC_00015061-mRNA-1 rna-XM_054371163
3 LOC_00014576-mRNA-1 rna-XM_054371163
4 LOC_00022770-mRNA-1 rna-XM_054371163
5 LOC_00005837-mRNA-1 rna-XM_054371163
6 LOC_00021813-mRNA-1 rna-XM_054371163

I checked the gff and I noticed something weird. For instance, if I check for instance SCNN1D, the "gene" name, rna-NR_037668, corresponds to the transcript "type",

NC_060925.1     BestRefSeq      transcript      709443  723551  .       +       .       ID=rna-NR_037668.3;Parent=gene-SCNN1D;Dbxref=GeneID:6339,GenBank:NR_037668.3,HGNC:HGNC:10601,MIM:601328;Name=NR_037668.3;gbkey=misc_RNA;gene=SCNN1D;product=sodium channel epithelial 1 subunit delta%2C transcript variant 2;transcript_id=NR_037668.3

while gffread used the mRNA "type" ID.

grep "SCNN1D" ../longert/T2T_longest_iso.gff
NC_060925.1     BestRefSeq      gene    709443  723551  .       +       .       ID=gene-SCNN1D;Dbxref=GeneID:6339,HGNC:HGNC:10601,MIM:601328;Name=SCNN1D;description=sodium channel epithelial 1 subunit delta;gbkey=Gene;gene=SCNN1D;gene_biotype=protein_coding;gene_synonym=dNaCh,ENaCd,ENaCdelta,SCNED
NC_060925.1     BestRefSeq      mRNA    709443  723551  .       +       .       ID=rna-NM_001130413.4;Parent=gene-SCNN1D;Dbxref=GeneID:6339,GenBank:NM_001130413.4,HGNC:HGNC:10601,MIM:601328;Name=NM_001130413.4;gbkey=mRNA;gene=SCNN1D;product=sodium channel epithelial 1 subunit delta%2C transcript variant 1;tag=RefSeq Select;transcript_id=NM_001130413.4

Perhaps I can modify something in the R code to make the reference table?

Lastly, I was reading again the GeneSymbol part, and it says to use the -S function of gffread to retrieve the proteins,

Let's start by using GFFRead to create predicted protein sequences from your target species. -y indicates that protein sequences are being generated and points to the new protein sequence file, -g is your softmasked genome sequence, and -S is the GFF file that you created in step 4.

gffread -y target_proteins.faa -g target_genome.softmasked.fasta -S full_annotation.gff

but it just changes the stop codon symbol from '.' to '*'.
Will Orthofinder run fine with proteins with internal stop codons?
I thought that those proteins should not be included for Orthofinder.

Thank you very much!

Hey!

You can definitely modify the R code to make the reference table.

Modifying the reference table would involve modifying this script https://github.com/BaderLab/GenAnT/blob/main/setup/ncbi_gene_key.R

Change line 28 to whatever column headers you like:
key <- gff_transcript[,c("gene_id","trans_id","gene")] # current default.

The more robust (and probably recommended) way to do it would be edit the gene symbols as a post-processing step. We have created a markdown for this and please see this issue for my answer to a very similar question: #4 (comment).

Regards,
Dustin

Navigation Menu

Repository navigation

Question about genome annotation QC #1

26 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Question about genome annotation QC #1

Description

Activity

14zac2 commented on Jan 8, 2026

desmodus1984 commented on Jan 8, 2026

DustinSokolowski commented on Jan 9, 2026

desmodus1984 commented on Jan 13, 2026

14zac2 commented on Jan 13, 2026

DustinSokolowski commented on Jan 14, 2026

desmodus1984 commented on Feb 11, 2026

14zac2 commented on Feb 11, 2026

desmodus1984 commented on Feb 11, 2026

DustinSokolowski commented on Feb 11, 2026

desmodus1984 commented on Feb 11, 2026

DustinSokolowski commented on Feb 11, 2026

desmodus1984 commented on Feb 11, 2026

DustinSokolowski commented on Feb 11, 2026

14zac2 commented on Feb 11, 2026

26 remaining items

14zac2 commented on Mar 23, 2026

desmodus1984 commented on Mar 25, 2026

DustinSokolowski commented on Mar 26, 2026

desmodus1984 commented on Apr 21, 2026

DustinSokolowski commented on Apr 22, 2026

14zac2 commented on Apr 23, 2026

desmodus1984 commented on Apr 23, 2026

desmodus1984 commented on May 29, 2026

14zac2 commented on May 29, 2026

desmodus1984 commented on May 29, 2026

14zac2 commented on May 30, 2026

desmodus1984 commented on Jun 2, 2026

DustinSokolowski commented on Jun 6, 2026

desmodus1984 commented on Jun 16, 2026

DustinSokolowski commented on Jun 22, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions