2. Method calibration
When running the GeneParliamentID pipeline for a lineage for the first time, method calibration using samples of known identity is highly recommended to identify the optimal pipeline parameters for this lineage, which will increase the accuracy of identification. The base command is:
gpid calibrate
Method calibration is structured into five subcommands:
gpid calibrate prepare: Prepare input files for calibration by matching each calibration sample against the reference dataset
gpid calibrate alignments: Identify optimal alignment filtering thresholds
gpid calibrate genes: Estimate gene performance and identify optimal gene threshold
gpid calibrate parliament: Identify optimal minimum parliament size threshold
gpid calibrate combine: Combine manually selected thresholds in a calibration file for subsequent use
For method calibration, you need to provide a dataset with samples of known identity, the Calibration dataset. This needs to be provided as a directory containing multiple unaligned genes, each comprising all retrieved sequences for the calibration samples.
- Calibration samples should be expert-identified, e.g. by taxonomists.
- Voucher specimens of the calibration samples should ideally be deposited in a museum or herbarium to allow validation of the identification.
- Reflect the taxonomic diversity of species across the lineage of interest. This is to ensure that the method is calibrated for the entire lineage, including taxonomically challenging species.
- Reflect different data qualities. This helps to identify minimum thresholds for data quality.
- Sufficient number of calibration samples. The more samples, the more reliable the method calibration will be.
- Gene sequences need to be fasta files ending with a
.FNAsuffix. - Gene names need to match the gene names in the reference dataset.
- Each calibration sample name needs to start with genus and species name
<Genus>_<species>, followed by a unique identifier for the sample. - Use underscores
_as separators in the sample names. - Sample names in each file need to be identical.
Gene1.FNA:
>Entandrophragma_angolense_Calibration_CQL_38
CTCCCATATAGGGGTGCTTGGTTGTGGGTGGGTTCTGAGATGATTCATTTAATGGAACTT
CCAAATCCAGATCCCTTAACTGGACGACCGGCACATGGTGGTCGAGATCGTCATACTTGT
ATTGCGATTCGAGATGTGTCTAAACTCAAAATAATCCTTGATAAAGCAGGTATTCCTTAC
>Entandrophragma_bussei_Calibration_CQL_EB63
CTCCCATATAGGGGTGCTTGGTTGTGGGTGGGTTCTGAGATGATTCATTTAATGGAACTT
CCAAATCCAGACCCCTTAACTGGACGACCGGCACATGGTGGTCGAGATCGTCATACTTGT
ATTGCGATTCGAGATGTGTCTAAACTCAAAATAATCCTTGATAAAGCAGGTATTCCTTAC
>Entandrophragma_candollei_Calibration_CQL_72
CTTCCGTATAGGGGTGCTTGGTTGTGGGTGGGTTCTGAGATGATTCATTTAATGGAACTT
CCAAATCCAGACCCCTTAACGGGACGACCGGCACATGGTGGTCGAGATCGTCATACTTGT
ATTGCGATTCGAGATGTGTCTAAACTGAAAATAATCCTTGACAAAGCAGGTGTTCCTTAC
Gene2.FNA:
>Entandrophragma_angolense_Calibration_CQL_38
GTTGGGGATCTACACGGAGACCTTGATCAAGCGAGATGTGCACTTGAGATTGCTGGTGTC
>Entandrophragma_bussei_Calibration_CQL_EB63
GTTGGGGATCTACACGGAGACCTTGATCAAGCGAGATGTGCACTTGAGATTGCTGGTGTC
>Entandrophragma_candollei_Calibration_CQL_72
GCTGGGGATCTACATGGAGACCTTGATCAAGCAAGATGTGCACTTGAGATTGCTGGTGTC
Gene3.FNA:
>Entandrophragma_angolense_Calibration_CQL_38
GAATCTCAAGGATTAACATCTCGACGATTGTGTCTTACATGTATATGTTCAACTCTAGCT
CTGATTAACAGTTCCGGCACGTTGGTTTCTGTACAAAAGGCAATTGCTTTGGAAGGAAAA
>Entandrophragma_bussei_Calibration_CQL_EB63
GATATGTGTGGTGGTACAGGAAAATGGAAAGCTCTCAACAGAAAACGTGCTAAAGATGTT
TACGAGTTTACAGAATGTCCAAATTGTTATGGTCGTGGGAAACTTGTGTGTCCGGTTTGC
>Entandrophragma_candollei_Calibration_CQL_72
GAATCTCAGATATCAACATCTCGCCGTTGGTGCCTTACGTGTATACTTACATGTATATGT
TCAACTCTAGCTCTGATTAACAGTTCCGGCACATTGGTTTCTGTACAAAAGGCAATTGCT
To prepare the calibration dataset for method calibration, run:
gpid calibrate prepare -r <reference dataset directory> -i <calibration dataset directory>
Required options:
-r Reference directory containing one FASTA file per gene and the corresponding BLAST databases. This is usually prepared with gpid reference, see 1. Reference construction.
-i Calibration dataset directory containing one FASTA file per gene
This matches each sample in the calibration dataset against all samples in the reference dataset using BLAST, gene by gene. The results from this are used as input for the subsequent method calibration analyses.
Two output files are created:
calibration/preparations/calibration_blast.tsvcalibration/preparations/calibration_prepared.rds
Method calibration is performed in three steps:
- Alignment filtering thresholds
- Gene performance threshold
- Parliament size threshold
For each calibration step, one or multiple figures and corresponding tables are produced to enable an informed decision on the parameters optimal for the given dataset. Using the samples of known identity, a large number of different settings is explored to assess the effect of different pipeline parameters and filtering thresholds on overall accuracy (the percentage of correctly identified samples) and retrievability (the percentage of samples passing the data filtering thresholds).
After each calibration step, the optimal thresholds need to be manually specified. This optimal threshold is then fixed for the remaining calibration steps. For example, once optimal alignment filtering thresholds have been specified, these thresholds are then applied to all subsequent calibration steps.
To conduct method calibration for the alignment filtering thresholds, run:
gpid calibrate alignments [-i <prepared calibration RDS>]
The optional flag -i specifies the RDS file generated by gpid calibrate prepare. By default, the input file created by gpid calibrate prepare is used.
This command explores the impacts of different filtering thresholds for the following BLAST alignment variables, assessing each variable independently:
- Minimum alignment similarity
- Minimum alignment length
- Maximum alignment gap openings
- Maximum alignment mismatches
- Maximum E-value
- Minimum Bit-score
For each variable, a pdf figure and csv file are written to calibration/tests/alignments/.
Based on the figures produced, the selected optimal thresholds for the dataset need to be manually specified for each variable by replacing NA with your selected value in the file calibration/manual_input_needed/calibration_alignments.tsv that is automatically produced.
For similarity, use values between 0 and 100. For E-value, use values between 0 and 1 using scientific notation, e.g. 1e-60. For all other values, use integer numbers >= 0.
| parameter | value |
|---|---|
| min_similarity | NA |
| min_length | NA |
| max_gapopens | NA |
| max_mismatches | NA |
| max_evalue | NA |
| min_bitscore | NA |
Guidelines:
- To decide which alignment filtering thresholds are optimal, we recommend balancing accuracy with retrievability.
- Keep in mind that stricter filtering thresholds often lead to higher accuracy but lower retrievability.
Example:
The below figure shows the impact of increasing minimum thresholds for alignment length on accuracy and retrievability.
Accuracy (solid line) initially slightly increases at a threshold of 100 bp, but decreases again at higher thresholds.
Retrievability (dashed line) remains initially at maximum until 800 bp, and decreases at higher thresholds.
Overall, a threshold of 100 bp is selected as optimal because it optimises accuracy whilst maintaining maximum retrievability.
To run method calibration for gene performance thresholds, run:
gpid calibrate genes [-i <prepared calibration RDS>] -a <alignment thresholds TSV>
Required:
-a Alignment thresholds TSV produced and manually edited after gpid calibrate alignments.
The default output path from gpid calibrate alignments is: calibration/manual_input_needed/calibration_alignments.tsv
Optional:
-i Intermediate RDS produced by gpid calibrate prepare
Default: calibration/preparations/calibration_prepared.rds
This command conducts two actions:
- Calculate gene performance for each gene, i.e. the percentage of samples correctly identified to species.
The gene performances are saved in the file:
calibration/calibration_gene_performance.csv
Example:
| gene | performance |
|---|---|
| Gene1 | 56.67 |
| Gene2 | 85.22 |
| Gene3 | 47.56 |
| ... | ... |
- Explore the impact of different minimum thresholds of gene performance on overall accuracy of identification and retrievability.
The results of this analysis are written as a figure and table to the directory:
calibration/tests/genes/
Based on the outputs produced, the selected optimal gene performance threshold needs to be manually entered in the automatically generated file calibration/manual_input_needed/calibration_genes.tsv by replacing NA with the selected threshold.
Use a value between 0 and 100 for the gene performance threshold.
| parameter | value |
|---|---|
| min_gene_performance | NA |
To perform method calibration for parliament size, i.e. the number of genes in the Gene Parliament, run:
gpid calibrate parliament [-i <prepared calibration RDS>] -a <alignment thresholds TSV> -g <gene threshold TSV>
Required:
-a Alignment thresholds file produced and manually edited after gpid calibrate alignments.
-g Gene threshold file produced and manually edited after gpid calibrate genes.
Optional:
-i Intermediate RDS file produced by gpid calibrate prepare
This step assesses the impact of different minimum thresholds of parliament size on overall accuracy of identification and retrievability.
The results of this analysis are written as a figure and table to the directory:
calibration/tests/parliament/
Based on the outputs produced, the selected optimal parliament size threshold needs to be manually entered in the automatically generated file calibration/manual_input_needed/calibration_parliament.tsv by replacing NA with the selected threshold.
Use an integer number >= 0 for the parliament size threshold.
| parameter | value |
|---|---|
| min_parliament_size | NA |
To combine the optimal thresholds selected during method calibration into a single calibration file, run:
gpid calibrate combine -a <alignment thresholds TSV> -g <gene threshold TSV> -p <parliament threshold TSV>
Required:
-a Alignment thresholds TSV produced and manually edited after gpid calibrate alignments
Default output path from gpid calibrate alignments:
calibration/manual_input_needed/calibration_alignments.tsv
-g Gene threshold TSV produced and manually edited after gpid calibrate genes
Default output path from gpid calibrate genes:
calibration/manual_input_needed/calibration_genes.tsv
-p Parliament threshold TSV produced and manually edited after gpid calibrate parliament
Default output path from gpid calibrate parliament:
calibration/manual_input_needed/calibration_parliament.tsv
This command produces the following calibration file that can be used as input for Method validation and Sample identification:
calibration/calibration_filtering_thresholds.csv
Example calibration file:
| min_similarity | min_length | max_gapopens | max_mismatches | max_evalue | min_bitscore | min_gene_performance | min_parliament_size |
|---|---|---|---|---|---|---|---|
| 98 | 100 | 1 | 5 | 1e-60 | 200 | 30 | 10 |
This file contains the thresholds that were identified as optimal for the given lineage.
Each variable is preceded by min or max, depending on whether the threshold is a minimum or maximum.
It is possible to run GeneParliamentID without conducting method calibration by using dummy calibration files with maximally lenient thresholds that bypass all filtering steps. This is generally not recommended as it will most likely result in considerably reduced accuracy of identification compared to running the pipeline with optimal parameters. However, it may be justified e.g. for a first explorative analysis or if a calibration dataset is not available.
See Bypass method calibration and validation.
Upon completion of Method calibration, continue to 3. Method validation.