This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Braker3 ETP annotations Filtering

Hi,

I have some fungal genome assemblies (hifi-reads based) which i annotated for structural gene predictions using Braker3 v3.0.8.

For annotations:

  • I have RNASeq data from 4 close strains (4 replicates each so 4x4=16 datasets)
  • I also Fungi proteins from orthodb12
  • Downloaded Protein sequences of 4 close samples (same as in RNASeq data) from NCBI and concatinated them to Fungi_ortho?db12.faa
  • Using Braker3 (which ran in ETP mode) i completed the annotations and i get the following results

From braker.gff3 in sample1

  • genes: 16829
  • mRNA: 31810
  • CDS: 31810
  • proteins: 31810 (from braker.aa file)

I ran BUSCO.v6 in proteins mode using ascomycetes_odb12 (2826 ortho genes) and got following stats

  • C:99.3%[S:4.0%,D:95.4%],F:0.2%,M:0.4%,n:2826
  • 2807 Complete BUSCOs (C)
  • 112 Complete and single-copy BUSCOs (S)
  • 2695 Complete and duplicated BUSCOs (D)
  • 7 Fragmented BUSCOs (F)
  • 12 Missing BUSCOs (M)
  • 2826 Total BUSCO groups searched

I busco results i get too much of duplicated genes and when i checked my .gff3 file i see that each gene has multiple transcripts (mRNA) but from different sources, like AUGUSTUS/GeneMark/gmst. For example the first 2 genes in gff3 look like this.

PB01_chr01  AUGUSTUS    gene    17425   18672   .   -   .   ID=g1;
PB01_chr01  AUGUSTUS    mRNA    17425   18672   1   -   .   ID=g1.t1;Parent=g1;
PB01_chr01  AUGUSTUS    stop_codon  17425   17427   .   -   0   ID=g1.t1.stop1;Parent=g1.t1;
PB01_chr01  AUGUSTUS    CDS 17428   18672   1   -   0   ID=g1.t1.CDS1;Parent=g1.t1;
PB01_chr01  AUGUSTUS    exon    17428   18672   .   -   .   ID=g1.t1.exon1;Parent=g1.t1;
PB01_chr01  AUGUSTUS    start_codon 18670   18672   .   -   0   ID=g1.t1.start1;Parent=g1.t1;
PB01_chr01  GeneMark.hmm3   mRNA    17425   18672   .   -   .   ID=g1.t2;Parent=g1;
PB01_chr01  GeneMark.hmm3   stop_codon  17425   17427   .   -   0   ID=g1.t2.stop1;Parent=g1.t2;
PB01_chr01  GeneMark.hmm3   CDS 17425   18672   .   -   0   ID=g1.t2.CDS1;Parent=g1.t2;
PB01_chr01  GeneMark.hmm3   exon    17425   18672   .   -   0   ID=g1.t2.exon1;Parent=g1.t2;
PB01_chr01  GeneMark.hmm3   start_codon 18670   18672   .   -   0   ID=g1.t2.start1;Parent=g1.t2;
PB01_chr01  AUGUSTUS    gene    21361   22664   .   +   .   ID=g2;
PB01_chr01  AUGUSTUS    mRNA    21361   22664   0.78    +   .   ID=g2.t1;Parent=g2;
PB01_chr01  AUGUSTUS    start_codon 21361   21363   .   +   0   ID=g2.t1.start1;Parent=g2.t1;
PB01_chr01  AUGUSTUS    CDS 21361   22256   1   +   0   ID=g2.t1.CDS1;Parent=g2.t1;
PB01_chr01  AUGUSTUS    exon    21361   22256   .   +   .   ID=g2.t1.exon1;Parent=g2.t1;
PB01_chr01  AUGUSTUS    intron  22257   22306   1   +   .   ID=g2.t1.intron1;Parent=g2.t1;
PB01_chr01  AUGUSTUS    CDS 22307   22559   0.78    +   1   ID=g2.t1.CDS2;Parent=g2.t1;
PB01_chr01  AUGUSTUS    exon    22307   22559   .   +   .   ID=g2.t1.exon2;Parent=g2.t1;
PB01_chr01  AUGUSTUS    intron  22560   22628   0.78    +   .   ID=g2.t1.intron2;Parent=g2.t1;
PB01_chr01  AUGUSTUS    CDS 22629   22661   0.78    +   0   ID=g2.t1.CDS3;Parent=g2.t1;
PB01_chr01  AUGUSTUS    exon    22629   22661   .   +   .   ID=g2.t1.exon3;Parent=g2.t1;
PB01_chr01  AUGUSTUS    stop_codon  22662   22664   .   +   0   ID=g2.t1.stop1;Parent=g2.t1;
PB01_chr01  gmst    mRNA    21361   22664   .   +   .   ID=g2.t2;Parent=g2;
PB01_chr01  gmst    start_codon 21361   21363   57.616532   +   0   ID=g2.t2.start1;Parent=g2.t2;
PB01_chr01  gmst    CDS 21361   22256   57.616532   +   0   ID=g2.t2.CDS1;Parent=g2.t2;
PB01_chr01  gmst    exon    21361   22256   57.616532   +   0   ID=g2.t2.exon1;Parent=g2.t2;
PB01_chr01  gmst    intron  22257   22306   57.616532   +   0   ID=g2.t2.intron1;Parent=g2.t2;
PB01_chr01  gmst    CDS 22307   22664   57.616532   +   1   ID=g2.t2.CDS2;Parent=g2.t2;
PB01_chr01  gmst    exon    22307   22664   57.616532   +   1   ID=g2.t2.exon2;Parent=g2.t2;
PB01_chr01  gmst    stop_codon  22662   22664   57.616532   +   0   ID=g2.t2.stop1;Parent=g2.t2;

Even if i look at the protein sequences predicted, they actually look the same, for example for the same sample above, protein sequences look like this

>g1.t1
MGYPYTRQSVARLAKLNDWDAFVEGAYERAFQQECLGTAIQADSRVWLELRPGVADLFAA
CKTCYMDFLANEIFANEYIPSVPPPGPNHQWTCALSQLSVRWALEAAVSQQDHAVFIEAV
RTISGLSPCTSAGITDGRWFTLIRDCPKFSICEGCYAGAIKSRGLDRFFVEIPLQPEELQ
TPMLCSFCPSAPRFERLFQKLYEAIDTDNFDAFSDFAVKFCRVPLCPRIGAWQNAKWWGY
PGLLFCEECYHDFVADTTLGNKLPMNGVVYKEYQMCQIWSPRMRRIWKEVCEAGPPGSAE
SDSALEEFKTFAAQRMSVYDQTIRQIEFLKQMQQIKNREAAFQGVMSVQYQGISSIASWG
SRDPYKYGNTSMGWWDNRYGAEASRRLDVMSSGFRNANNMSREIIRLRGIWETVE
>g1.t2
MGYPYTRQSVARLAKLNDWDAFVEGAYERAFQQECLGTAIQADSRVWLELRPGVADLFAA
CKTCYMDFLANEIFANEYIPSVPPPGPNHQWTCALSQLSVRWALEAAVSQQDHAVFIEAV
RTISGLSPCTSAGITDGRWFTLIRDCPKFSICEGCYAGAIKSRGLDRFFVEIPLQPEELQ
TPMLCSFCPSAPRFERLFQKLYEAIDTDNFDAFSDFAVKFCRVPLCPRIGAWQNAKWWGY
PGLLFCEECYHDFVADTTLGNKLPMNGVVYKEYQMCQIWSPRMRRIWKEVCEAGPPGSAE
SDSALEEFKTFAAQRMSVYDQTIRQIEFLKQMQQIKNREAAFQGVMSVQYQGISSIASWG
SRDPYKYGNTSMGWWDNRYGAEASRRLDVMSSGFRNANNMSREIIRLRGIWETVE*
>g2.t1
MDPNLELYRSILDLGPKERRQRMQHLPKEELIRVKSIVEREKWIQMLETAVAGRDLVELA
FTDPREIQENPPFQKALLGRACYPDDENNMVKRITKGLRKNGESLIHTVASFDGPTYPAI
TKDAWILVYCDLFYIDGNNMTLHEVYTSRLQEEELQTRTEQAREVARHDDLKKARRNAKW
MIPALGRLSDEELSQSEYDFSNTLHEIWKQVSHAPSTWIQHILDAQQPWGFTYYKTKQVE
EKYGRTWKDTWIMIIDMPQQSWSSIHCQGKVHEFMELKTEDWAPPPTYEGLTEDDAFRKH
FREHRKSLSSPGILQNTFIVIPIELIPDDPDDDELDLLWVWAYDADWDSSSEEIICNGEK
YQGRIKVPLYALEAWFYAARWEGELEDWDHEPYV
>g2.t2
MDPNLELYRSILDLGPKERRQRMQHLPKEELIRVKSIVEREKWIQMLETAVAGRDLVELA
FTDPREIQENPPFQKALLGRACYPDDENNMVKRITKGLRKNGESLIHTVASFDGPTYPAI
TKDAWILVYCDLFYIDGNNMTLHEVYTSRLQEEELQTRTEQAREVARHDDLKKARRNAKW
MIPALGRLSDEELSQSEYDFSNTLHEIWKQVSHAPSTWIQHILDAQQPWGFTYYKTKQVE
EKYGRTWKDTWIMIIDMPQQSWSSIHCQGKVHEFMELKTEDWAPPPTYEGLTEDDAFRKH
FREHRKSLSSPGILQNTFIVIPIELIPDDPDDDELDLLWVWAYDADWDSSSEEIICNGEK
YQGRIKVPLYALEAWFYAARWEGVSLRDMWLKAQTHEDNLWICHSKELEDWDHEPYV*

I tried to look for any post which might indicate if i need to do any curation of braker3 results before moving down to functional annotation, but couldn't find any.

Braker3 has an optiton where it takes a BUSCO_db as input and runs compleasm. But a quick ChatGPT to understand this option showed that it is only to find missing BUSCO genes, and tends to make results BUSCO-ish. Tho i tried it but it didnot give any huge difference ( ithought it might help with multiple transcripts and duplications)

From braker.gff3 in sample1 after compleasm

  • genes: 16838
  • mRNA: 31997
  • CDS: 31997
  • proteins: 31997 (from braker.aa file)

and busco scores after brraker3 with compleasm

  • C:99.5%[S:3.0%,D:96.4%],F:0.2%,M:0.3%,n:2826
  • 2811 Complete BUSCOs (C)
  • 86 Complete and single-copy BUSCOs (S)
  • 2725 Complete and duplicated BUSCOs (D)
  • 7 Fragmented BUSCOs (F)
  • 8 Missing BUSCOs (M)
  • 2826 Total BUSCO groups searched

Any help in this ith what should i do next is really appreciated. I have tried AGAT for longest isoforms and tho results look good but i am not sure is this is something like Best-Practice in fungal genomes. After AGAT results look like this

From original (no compleasm) braker.gff3 after running AGAT longest isoform

  • genes: 16829
  • mRNA: 16829
  • CDS: 46261
  • proteins: 16829 (from braker.aa file)

and BUSCO results

  • C:99.3%[S:97.6%,D:1.8%],F:0.2%,M:0.5%,n:2826
  • 2807 Complete BUSCOs (C)
  • 2757 Complete and single-copy BUSCOs (S)
  • 50 Complete and duplicated BUSCOs (D)
  • 6 Fragmented BUSCOs (F)
  • 13 Missing BUSCOs (M)
  • 2826 Total BUSCO groups searched

for agat i used these commands

# 1. AGAT: get statistics of the Braker3 gene predictions
agat_sp_statistics.pl -i "${currentBRAKERGFF}" -g "${currentFASTA}" -o "${brakerQCout}/${currentSAMPLE}.braker.ETP.gff3.stats"

# 2. AGAT: keep only the longest isoform per gene from the Braker3 gff3 file
agat_sp_keep_longest_isoform.pl -g "${currentBRAKERGFF}" -o "${brakerQCout}/${currentSAMPLE}.braker.ETP.longestIsoform.gff3"

# 3. extract the corresponding CDS sequences for the longest isoform genes using AGAT
agat_sp_extract_sequences.pl -g "${brakerQCout}/${currentSAMPLE}.braker.ETP.longestIsoform.gff3" \
    -f "${currentFASTA}" -o "${brakerQCout}/${currentSAMPLE}.braker.ETP.longestIsoform.cds" -t cds

# 4. extract the corresponding protein sequences for the longest isoform genes using AGAT
agat_sp_extract_sequences.pl -g "${brakerQCout}/${currentSAMPLE}.braker.ETP.longestIsoform.gff3" \
    -f "${currentFASTA}" -o "${brakerQCout}/${currentSAMPLE}.braker.ETP.longestIsoform.aa" -t cds -p

If you can assist or point me in right direction, it will be great.

braker3 annotations fungi
• 1,280 views • link
written by SomeOne • SomeOne

i annotated for structural gene predictions

Can you clarify where all this analysis leading. Work you describe above is thorough but these are predictions. They remain just that, until they are proven by experimental evidence. No amount of bioinformatics juggling is going to get to the "truth" or "gold standard", if that is what you are looking for.

Your best bet would be to do some RNAseq experiments and get experimental data so you can compare that data to validate these predictions.

• 0 views • link
written by GenoMax GenoMax

From the completed annotations i will move toward the pangenome side. But also look for genes which are unique to each specie/strain. because some of the sample are clinical and some environmental. So having diuplicates maked this confusing.

• 0 views • link
written by SomeOne • SomeOne

From a quick dialog with GPT it looks like braker should not produce redundant annotations. genemark annotations are created for training and final gene models should only be from AUGUSTUS. In the example you show above more than one program models appear to be present. Is that the actual output of braker or have you done something different to create a merged annotation file.

Technically you don't need annotations (and you do have good gene models) to create the pangenome. Consider using a RefSeq genome as an anchor.

• 0 views • link
written by GenoMax GenoMax

From what i learned from Braker3 docs, If you provide RNASeq data and Protein sequences as evidance, it will (coppied from github docs)

  1. This will call GeneMark-ETP, which will use RNA-Seq and protein hints for training GeneMark-ETP
  2. Subsequently, AUGUSTUS is trained on 'high-confindent' genes (genes with very high extrinsic evidence support) from the GeneMark-ETP prediction and a set of genes is predicted by AUGUSTUS
  3. In a last step, the predictions of AUGUSTUS and GeneMark-ETP are combined using TSEBRA.

Thats why it has multiple sources in final annotation files.

My guess is that i can use the Output files from AUGUSTUS skipping the step3 results.

• 0 views • link
written by SomeOne • SomeOne

My guess is that i can use the Output files from AUGUSTUS skipping the step3 results.

Possibly. Looking at your example data it appears that AUGUSTUS and GeneMark results are almost identical except difference in CDS coordinates. gene entry coordinates are identical.

Where did the proteins set come for the BRAKER run. Were they from OrthoDB or from your assembled genome?

Annotation is hard to do as you are finding out. $1K genome requires a $100K annotation effort.

• 0 views • link
written by GenoMax GenoMax

Where did the proteins set come for the BRAKER run. Were they from OrthoDB or from your assembled genome?

I downloaded Fungi proteins from Ortho_db12 and also added some previously annotated proteins from close references, (downloaded sequences from NCBI-RefSeq and concatinated them to orthodb.aa file) at end the protein sequences file was a mix of diverse and also really close sequences.

• 0 views • link
written by SomeOne • SomeOne

Annotation is hard to do as you are finding out. $1K genome requires a $100K annotation effort.

i am finding way too much about the whole Genome assembly to annotations, in the course work i only knew about Velvet for genome assemblies

• 0 views • link
written by SomeOne • SomeOne

0 answers

No answers yet.

Log in to answer this question.

More posts like this