Skip to content

EVidenceModeler (EVM)

Jump to bottom
swarbred edited this page Apr 29, 2024 · 2 revisions

EVidenceModeler (EVM) commands

1) Partition the EVM inputs

cd /home/train/Annotation_workshop/EVM
partition_EVM_inputs.pl --genome Inputs/Reference/Athaliana_447_TAIR10_Chr3_clean.fa --gene_predictions Inputs/Models/gene_models_list_A.gff3 --protein_alignments Inputs/Evidence/protein_alignments.gff3 --transcript_alignments Inputs/Evidence/transcript_alignments.gff3 --segmentSize 500000 --overlapSize 50000 --partition_listing partitions_list.out

2) Extract just the partitions that contain genes in our example region

grep "Chr3_900001-1400000\|Chr3_1350001-1850000" partitions_list.out > my_partitions_list.out

3) Create the EVM commands file

write_EVM_commands.pl --genome Inputs/Reference/Athaliana_447_TAIR10_Chr3_clean.fa --gene_predictions Inputs/Models/gene_models_list_A.gff3 --protein_alignments Inputs/Evidence/protein_alignments.gff3 --transcript_alignments Inputs/Evidence/transcript_alignments.gff3 --weights $PWD/Inputs/Configs/weights.txt --search_long_introns 5000 --output_file_name evm.out  --partitions my_partitions_list.out > commands.list

4) We can run the EVM commands for the two partitions

Scripts/run.sh commands.list

5) Gather together the EVM outputs

recombine_EVM_partial_outputs.pl --partitions my_partitions_list.out --output_file_name evm.out

6) Convert the EVM output to gff3

convert_EVM_outputs_to_GFF3.pl  --partitions my_partitions_list.out --output evm.out --genome Inputs/Reference/Athaliana_447_TAIR10_Chr3_clean.fa

7) Generate the annotation stats for the EVM models

mikado util stats Chr3/evm.out.gff3 evm.out.gff3.stats

8) Summarise the annotation stats

parse_mikado_stats evm.out.gff3.stats > evm.out.gff3.stats.summary

9) Compare annotation stats between the reference annotation and the EVM models

Scripts/paste_mikado_summary_stats.sh Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf.stats.summary evm.out.gff3.stats.summary
	Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf.stats.summary	evm.out.gff3.stats.summary
Number of genes                   122.00   106.00
Number of Transcripts             221.00   106.00
Transcripts per gene              1.78     1.00
Number of monoexonic genes        23.00    23.00
Monoexonic transcripts            29.00    23.00
Transcript mean size cDNA (bp)    1806.26  1332.25
Transcript median size cDNA (bp)  1714.00  1174.00
Min cDNA                          162.00   195.00
Max cDNA                          5588.00  5112.00
Total exons                       1579.00  645.00
Exons per transcript              7.14     6.08
Exon mean size (bp)               252.81   218.94
CDS mean size (bp)                190.14   218.94
Transcript mean size CDS (bp)     1288.83  1332.25
Transcript median size CDS (bp)   1209.00  1174.00
Min CDS                           108.00   195.00
Max CDS                           5112.00  5112.00
Intron mean size (bp)             155.22   149.45
5'UTR mean size (bp)              250.33   0.00
3'UTR mean size (bp)              267.10   0.00

10) Run mikado compare to to see how well the EVM models match the reference annotation

mikado compare -r Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf -p Chr3/evm.out.gff3 -o mikado_compare.evm

11) Repeat but limiting the comparison to just CDS features

mikado compare -r Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf -p Chr3/evm.out.gff3 -eu -o mikado_compare_eu.evm

Add UTRs and splice variaints using Mikado

12) Run mikado configure

mikado configure --full --list Inputs/Configs/list.txt --reference Inputs/Reference/Athaliana_447_TAIR10_Chr3_clean.fa --scoring plant.yaml --copy-scoring plant.yaml --only-reference-update configuration.toml

13) Run mikado prepare

mikado prepare -p 1 --json-conf configuration.toml

14) Identify the models with no CDS features and generate a cDNA fasta file of these transcripts

gffread --nc -T -o mikado_prepared.nc.gtf -w mikado_prepared.nc.fasta -g Inputs/Reference/Athaliana_447_TAIR10_Chr3_clean.fa  mikado_prepared.gtf

15) Run prodigal to call ORFs

prodigal -g 1 -f gff -i mikado_prepared.nc.fasta -o mikado_prepared.cds.gff

16) Load ORFs and junctions

mikado serialise -p 1 --json-conf configuration.toml --junctions Inputs/Mikado/portcullis.pass.junctions.regionA.bed --orfs mikado_prepared.cds.gff

17) Run mikado pick to select gene models

mikado pick -p 4 --output-dir Mikado --json-conf configuration.toml

18) Run mikado compare to to see how well the EVM-mikado models match the reference annotation

mikado compare -r Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf -p Mikado/mikado.loci.gff3 -o mikado_compare.evm-mikado

19) Compare mikado compare results between EVM and EVM-mikado models

paste mikado_compare.evm.stats mikado_compare.evm-mikado.stats |expand -t 60
Command line:                                               Command line:
/usr/local/bin/mikado compare -r Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf -p Chr3/evm.out.gff3 -o mikado_compare.evm                                    /usr/local/bin/mikado compare -r Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf -p Mikado/mikado.loci.gff3 -o mikado_compare.evm-mikado
221 reference RNAs in 122 genes                             221 reference RNAs in 122 genes
106 predicted RNAs in  106 genes                            147 predicted RNAs in  106 genes
--------------------------------- |   Sn |   Pr |   F1 |    --------------------------------- |   Sn |   Pr |   F1 |
                        Base level: 61.58  97.72  75.55                             Base level: 78.72  93.05  85.29
            Exon level (stringent): 51.39  71.78  59.90                 Exon level (stringent): 53.94  67.50  59.96
              Exon level (lenient): 76.06  97.76  85.55                   Exon level (lenient): 81.72  93.32  87.14
                 Splice site level: 86.56  99.17  92.43                      Splice site level: 90.85  96.14  93.42
                      Intron level: 90.65  99.07  94.67                           Intron level: 93.74  95.46  94.59
                 Intron level (NR): 84.23  99.07  91.05                      Intron level (NR): 88.80  94.62  91.62
                Intron chain level: 38.54  86.75  53.37                     Intron chain level: 49.48  72.09  58.68
           Intron chain level (NR): 37.89  86.75  52.75                Intron chain level (NR): 48.95  72.09  58.31
      Transcript level (stringent): 0.90  1.89  1.22              Transcript level (stringent): 0.90  1.36  1.09
  Transcript level (>=95% base F1): 3.17  6.60  4.28          Transcript level (>=95% base F1): 29.41  43.54  35.11
  Transcript level (>=80% base F1): 26.24  53.77  35.27       Transcript level (>=80% base F1): 46.61  68.71  55.54
         Gene level (100% base F1): 1.64  1.89  1.75                 Gene level (100% base F1): 1.64  1.89  1.75
        Gene level (>=95% base F1): 5.74  6.60  6.14                Gene level (>=95% base F1): 47.54  54.72  50.88
        Gene level (>=80% base F1): 46.72  53.77  50.00             Gene level (>=80% base F1): 74.59  85.85  79.82

#   Matching: in prediction; matched: in reference.         #   Matching: in prediction; matched: in reference.

            Matching intron chains: 72                                  Matching intron chains: 93
             Matched intron chains: 74                                   Matched intron chains: 95
   Matching monoexonic transcripts: 14                         Matching monoexonic transcripts: 15
    Matched monoexonic transcripts: 14                          Matched monoexonic transcripts: 16
        Total matching transcripts: 86                              Total matching transcripts: 108
         Total matched transcripts: 88                               Total matched transcripts: 111

          Missed exons (stringent): 438/901  (48.61%)                 Missed exons (stringent): 415/901  (46.06%)
           Novel exons (stringent): 182/645  (28.22%)                  Novel exons (stringent): 234/720  (32.50%)
            Missed exons (lenient): 192/802  (23.94%)                   Missed exons (lenient): 147/804  (18.28%)
             Novel exons (lenient): 14/624  (2.24%)                      Novel exons (lenient): 47/704  (6.68%)
                    Missed introns: 100/634  (15.77%)                           Missed introns: 71/634  (11.20%)
                     Novel introns: 5/539  (0.93%)                               Novel introns: 32/595  (5.38%)

       Missed transcripts (0% nF1): 25/221  (11.31%)               Missed transcripts (0% nF1): 24/221  (10.86%)
        Novel transcripts (0% nF1): 3/106  (2.83%)                  Novel transcripts (0% nF1): 5/147  (3.40%)
             Missed genes (0% nF1): 19/122  (15.57%)                     Missed genes (0% nF1): 18/122  (14.75%)
              Novel genes (0% nF1): 3/106  (2.83%)                        Novel genes (0% nF1): 2/106  (1.89%)

20) Generate the annotation stats for the EVM-mikado models

mikado util stats Mikado/mikado.loci.gff3 mikado.loci.gff3.stats

21) Summarise the annotation stats

parse_mikado_stats mikado.loci.gff3.stats > mikado.loci.gff3.stats.summary

22) Compare annotation stats between the reference, EVM and EVM-mikado annotations

Scripts/paste_mikado_summary_stats.sh Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf.stats.summary evm.out.gff3.stats.summary mikado.loci.gff3.stats.summary
	Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf.stats.summary	evm.out.gff3.stats.summary	mikado.loci.gff3.stats.summary
Number of genes                   122.00   106.00   106.00
Number of Transcripts             221.00   106.00   147.00
Transcripts per gene              1.78     1.00     1.39
Number of monoexonic genes        23.00    23.00    17.00
Monoexonic transcripts            29.00    23.00    18.00
Transcript mean size cDNA (bp)    1806.26  1332.25  1824.03
Transcript median size cDNA (bp)  1714.00  1174.00  1710.00
Min cDNA                          162.00   195.00   243.00
Max cDNA                          5588.00  5112.00  5609.00
Total exons                       1579.00  645.00   1006.00
Exons per transcript              7.14     6.08     6.84
Exon mean size (bp)               252.81   218.94   266.53
CDS mean size (bp)                190.14   218.94   208.31
Transcript mean size CDS (bp)     1288.83  1332.25  1343.40
Transcript median size CDS (bp)   1209.00  1174.00  1209.00
Min CDS                           108.00   195.00   195.00
Max CDS                           5112.00  5112.00  5112.00
Intron mean size (bp)             155.22   149.45   153.11
5'UTR mean size (bp)              250.33   0.00     212.05
3'UTR mean size (bp)              267.10   0.00     268.58

After the above commands to cleanup i.e. return to only the original files.

IMPORTANT: Before deleting, please make sure that you are in the EVM folder under /home/train/Annotation_workshop/EVM. Otherwise, other data files will be deleted

rm -v !("Inputs"|"Scripts"|"commands.txt"|"Example_output"); rm -rf Chr3 Mikado

Clone this wiki locally