AnnoAudit - Annotation Auditor

AnnoAudit is a robust Nextflow pipeline designed to evaluate the quality of genomic annotations through a multifaceted approach.

Overview of the workflow

The workflow assess the annotation quality based on different criteria:

• Protein evidence support
• RNASeq evidence support
• Statistics of the predictions (i.e., gene length, exon number, etc.)
• Ortholog analysis (BUSCO, OMArk)

Input data

• Reference genome genome.[.fna, .fa, .fasta]
• Annotation output annotation.gff
• RNAseq data listed in a metadata csv file. Input type can be mixed between long and short reads, with the option of single-end read. The input file should follow the format below:

sample_id,R1_path,R2_path,read_type
SAM1,/path/to/R1,,long             # For long reads
SAM2,/path/to/R1,/path/to/R2,short # For PE reads
SAM3,/path/to/R1,,short            # For SE reads

• Protein reference data in fasta format for evaluation, if not given, then the Uniprot-SwissProt will be downloaded and used.

Pipeline steps

The main pipeline is divided into five different subworkflows.

• General statistics: Calculate the statistics obtained from the GFF file.
• RNASeq analysis: Map the RNASeq data to the genome (or with provided mapping bam file) to generate exon, intron, transcript coverage.
• Ortholog analysis: Compare the predicted proteome to known database using BUSCO and OMArk (OMA database).
• Protein analysis: Blast the predicted proteome to a known database (could be of relative species) to obtain best reciprocal hits (BRH), then generate statistics based on the BRH results.

Output data

• Output text file contain the statistic calculated from the input GFF file:
  • General statistics
  • BUSCO
  • OMArk
  • PSAURON
  • Best reciprocal hits
  • RNASeq analysis

Prerequisites

The following programs are required to run the workflow and the listed version were tested.

nextflow v23.04.0 or higher

singularity

docker (have not been tested but in theory should work fine)

Installation

Simply get the code from github or workflowhub and directly use it for the analysis with nextflow.

git clone https://github.com/ERGA-consortium/pipelines

Running AnnoAudit

Before running the pipeline (IMPORTANT)

One thing with Nextflow is that it is running off a Java Virtual Machine (JVM), and it will try to use all available memory for Nextflow even though it is unnecessary (for workflow management and job control). This will cause much trouble if you run a job on an HPC cluster. Thus, to minimize the effect of it, we need to limit the maximum memory the JVM can use.

export NFX_OPTS="-Xms=512m -Xmx=3g"

-Xms is the lower limit, which is set as 512 MB. -Xmx is the upper limit, which in this case is set as 3 GB. Please modify this according to your situation.

How to run the code

nextflow run main.nf --genome genome.fasta \
      --gff annotation.gff3 \
      --rnaseq metadata.csv [--genome_bam path/to/the/mapped/bam]\
      --outdir OUTDIR_NAME \
      --taxon_id 9606 [Optional] \
      --ncbi_query_email xxxx \
      --rm -resume

Other parameters for running the analysis

Input parameter:
--genome                  Draft genome fasta file contain the assembled contigs/scaffolds.
--gff                     Annotation file that need to be evaluated.
--genome_bam              BAM file contain the mapped information from the RNASeq to the genome FASTA.
--rnaseq                  A metadata CSV file following the pattern: sample_id,R1_path,R2_path,read_type. Required if `genome_bam` is not provided.
--taxon_id                Taxon ID for identifying BUSCO lineage and download protein data from NCBI if needed.
--ncbi_query_email        Email for querying protein from NCBI database.

Optional input:
--protein                  Fasta file containing translated protein sequences from the GFF for running evaluation. If not specified, the workflow will automatically extract it from the
 `genome` and `gff`.
--ref_protein              Fasta file containing the reference protein sequences to be used for evaluation. Ideally this should come from the same species and/or closely related specie
s. If not provided, the workflow will download the proteome from NCBI or using Uniprot SwissProt database.
--lineage                  Lineage information providing for BUSCO, if not provided, the workflow will automatically search for the closest lineage. Example: eudicots_odb10.
--genetic_code             Genetic code for translation of protein.
--stranding                Strandness of the RNASeq reads used for extraction of junction position using `regtools`.

Database input:
--odb_version              odb version to choose to run BUSCO, option: odb12, odb10. [default: odb12]
--busco_database           Pathway to the BUSCO databse store locally. [default: null]
--oma_database             Pathway to the OMA database, if not specified, the workflow will download it automatically. [default: null]
--ref_protein              Pathway to the reference proteome for comparison. [default: null]
--ncbi_query_count         Number of protein to extract from the NCBI database. [default: 100000]
--ncbi_query_batch         Number of protein to query for each batch. [default: 1000]

Output option:
--pdf                      Output PDF name. [default: AnnoAudit_Report.pdf]
--outdir                   Output directory. [default: /env/export/bigtmp2/pdoan/evaluate_pipeline]
--tracedir                 Pipeline information. [default: /env/export/bigtmp2/pdoan/evaluate_pipeline/pipeline_info]
--publish_dir_mode         Option for nextflow to move data to the output directory. [default: copy]
--tmpdir                   Database directory. [default: /env/export/bigtmp2/pdoan/evaluate_pipeline/tmpdir]

Conditioning options:
--rnaseq_single             If specify, will run `featureCounts` in single read mode, this is necessary if the mapped RNASeq is single-ended. [default: false]
--run_blast                 If specify, will use `blast` for running best reciprocal hits instead of DIAMOND. [default: false]
--query_ncbi_prot           If specify, will download the reference proteome from NCBI, other wise, will use the provided proteom or Uniprot SwissProt. [default: true]
--cds_only                  If specify, only extracting information from the GFF file using the CDS line. [default: "False"]

--help                   Print help message.

Execution/Engine profiles:
The pipeline supports profiles to run via different Executers and Engines e.g.: -profile local,conda

Executer (choose one):
  local
  slurm

Engines (choose one):
  docker
  singularity
  apptainer

Per default: -profile slurm,singularity is executed.

Example output

Below is the sample output of this workflow. The example PDF output is located in assest folder.

|General Statistics                 | Value           |
-------------------------------------------------------
|num_genes                          | 36391           |
|num_genes_without_introns          | 12968 (35.64%)  |
|mean_gene_length                   | 2359.57         |
|median_gene_length                 | 1562.0          |
|num_exons                          | 149725          |
|mean_exons_per_gene                | 4.11            |
|median_exons_per_gene              | 2.0             |
|num_exon_3n                        | 76783 (51.28%)  |
|num_exon_3n1                       | 36932 (24.67%)  |
|num_exon_3n2                       | 36010 (24.05%)  |
|mean_cds_length                    | 1091.4          |
|median_cds_length                  | 873.0           |
|total_cds_length                   | 39717145        |
|percentage_cds_coverage            | 10.64%          |
|num_introns                        | 113334          |
|mean_intron_length                 | 407.2           |
|median_intron_length               | 149.0           |
|short_intron_<120_3n0_without_stop | 4324 (3.82)%    |
|long_intron_>120_3n0_without_stop  | 1185 (1.05)%    |
|short_intron_<120_3n1_without_stop | 4205 (3.71)%    |
|long_intron_>120_3n1_without_stop  | 1291 (1.14)%    |
|short_intron_<120_3n2_without_stop | 4319 (3.81)%    |
|long_intron_>120_3n2_without_stop  | 1249 (1.10)%    |
|short_intron_<120_3n0_with_stop    | 12073 (10.65)%  |
|long_intron_>120_3n0_with_stop     | 20332 (17.94)%  |
|short_intron_<120_3n1_with_stop    | 11652 (10.28)%  |
|long_intron_>120_3n1_with_stop     | 20486 (18.08)%  |
|short_intron_<120_3n2_with_stop    | 11733 (10.35)%  |
|long_intron_>120_3n2_with_stop     | 20485 (18.07)%  |

|BUSCO                              | Value           |
-------------------------------------------------------
|lineage_dataset                    | poales_odb10    |
|complete                           | 97.6%           |
|single_copy                        | 95.8%           |
|multi_copy                         | 1.8%            |
|fragmented                         | 0.2%            |
|missing                            | 2.2%            |
|num_markers                        | 4896            |
|domain                             | eukaryota       |


|OMARK                              | Value           |
-------------------------------------------------------
|OMA_clade                          | Oryza           |
|num_conserved_hogs                 | 15087           |
|single                             | 13316 (88.26%)  |
|duplicated                         | 1353 (8.97%)    |
|duplicated_unexpected              | 1101 (7.30%)    |
|duplicated_expected                | 252 (1.67%)     |
|missing                            | 418 (2.77%)     |
|num_proteins_in_proteome           | 36387           |
|total_consistent                   | 30365 (83.45%)  |
|consistent_partial_hits            | 1803 (4.96%)    |
|consistent_fragmented              | 1625 (4.47%)    |
|total_inconsistent                 | 2283 (6.27%)    |
|inconsistent_partial_hits          | 517 (1.42%)     |
|inconsistent_fragmented            | 1444 (3.97%)    |
|total_contaminants                 | 0 (0.00%)       |
|contaminants_partial_hits          | 0 (0.00%)       |
|contaminants_fragmented            | 0 (0.00%)       |
|total_unknown                      | 3739 (10.28%)   |

|PSAURON                            | Value           |
-------------------------------------------------------
|psauron_score                      | 83.8            |
|true_count                         | 30494           |
|false_count                        | 5893            |
|median_score                       | 0.98278         |
|max_score                          | 1.0             |
|min_score                          | 0.00022         |


|Best Reciprocal Hits               | Value           |
-------------------------------------------------------
|num_best_reciprocal_hits           | 29185           |
|num_splitting_genes_08             | 932 (3.19%)     |
|num_splitting_genes_05             | 0 (0.0%)        |
|num_fusion_genes_12                | 437 (1.5%)      |
|num_fusion_genes_15                | 482 (1.65%)     |
|KL_divergence_normalized           | 0.0105          |
|JS_divergence_normalized           | 0.0023          |
|Wasserstein_distance               | 2.480915        |


|RNASeq                             | Value           |
-------------------------------------------------------
|mapping_rate                       | 96.27%          |
|primary_mapping_rate               | 95.83%          |
|properly_paired                    | 92.47%          |
|num_gene_unsupported               | 9445 (25.95%)   |
|num_exon_unsupported               | 20232 (13.51%)  |
|num_intron_supported               | 107202          |
|num_intron_supported_canonical     | 107131 (99.93%) |
|num_intron_supported_non_canonical | 71 (0.07%)      |

Performance of the workflow on assessing annotation

To be added

Future work

• Adding other plots for easier evaluation
• Perform comparative performance with different genomes