2. Method calibration
Before conducting Method calibration, make sure to have completed 1. Reference construction.
Method calibration is used to identify the optimal parameters for identification using GPID for any lineage and set of genes. Conducting method calibration can substantially improve the accuracy of identification.
Method calibration involves running GPID using samples of known species identity. It explores a large number of different settings 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). For each parameter that can be changed, a plot of accuracy and retrievability depending on the filtering threshold is produced. This enables you to make an informed decision on the optimal filtering thresholds:
- Select optimal alignment filtering thresholds
- Calculate gene performance (percentage of correct identifications) and select optimal gene performance threshold
- Select minimum parliament size threshold
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.
For method calibration, you need to provide a calibration dataset with samples of known identity. 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. While we can't offer a definitive number, we found that approximately 100 samples can allow good insights into the effect of different calibration decisions.
- Gene sequences need to be fasta files ending with a
.FNAsuffix. - Gene names need to match the gene names in the reference dataset.
- Each test 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_Test_CQL_38
CTCCCATATAGGGGTGCTTGGTTGTGGGTGGGTTCTGAGATGATTCATTTAATGGAACTT
CCAAATCCAGATCCCTTAACTGGACGACCGGCACATGGTGGTCGAGATCGTCATACTTGT
ATTGCGATTCGAGATGTGTCTAAACTCAAAATAATCCTTGATAAAGCAGGTATTCCTTAC
>Entandrophragma_bussei_Test_CQL_EB63
CTCCCATATAGGGGTGCTTGGTTGTGGGTGGGTTCTGAGATGATTCATTTAATGGAACTT
CCAAATCCAGACCCCTTAACTGGACGACCGGCACATGGTGGTCGAGATCGTCATACTTGT
ATTGCGATTCGAGATGTGTCTAAACTCAAAATAATCCTTGATAAAGCAGGTATTCCTTAC
>Entandrophragma_candollei_Test_CQL_72
CTTCCGTATAGGGGTGCTTGGTTGTGGGTGGGTTCTGAGATGATTCATTTAATGGAACTT
CCAAATCCAGACCCCTTAACGGGACGACCGGCACATGGTGGTCGAGATCGTCATACTTGT
ATTGCGATTCGAGATGTGTCTAAACTGAAAATAATCCTTGACAAAGCAGGTGTTCCTTAC
Gene2.FNA:
>Entandrophragma_angolense_Test_CQL_38
GTTGGGGATCTACACGGAGACCTTGATCAAGCGAGATGTGCACTTGAGATTGCTGGTGTC
>Entandrophragma_bussei_Test_CQL_EB63
GTTGGGGATCTACACGGAGACCTTGATCAAGCGAGATGTGCACTTGAGATTGCTGGTGTC
>Entandrophragma_candollei_Test_CQL_72
GCTGGGGATCTACATGGAGACCTTGATCAAGCAAGATGTGCACTTGAGATTGCTGGTGTC
Gene3.FNA:
>Entandrophragma_angolense_Test_CQL_38
GAATCTCAAGGATTAACATCTCGACGATTGTGTCTTACATGTATATGTTCAACTCTAGCT
CTGATTAACAGTTCCGGCACGTTGGTTTCTGTACAAAAGGCAATTGCTTTGGAAGGAAAA
>Entandrophragma_bussei_Test_CQL_EB63
GATATGTGTGGTGGTACAGGAAAATGGAAAGCTCTCAACAGAAAACGTGCTAAAGATGTT
TACGAGTTTACAGAATGTCCAAATTGTTATGGTCGTGGGAAACTTGTGTGTCCGGTTTGC
>Entandrophragma_candollei_Test_CQL_72
GAATCTCAGATATCAACATCTCGCCGTTGGTGCCTTACGTGTATACTTACATGTATATGT
TCAACTCTAGCTCTGATTAACAGTTCCGGCACATTGGTTTCTGTACAAAAGGCAATTGCT
As a preparation for method calibration, we need to match all calibration dataset against the reference dataset using BLAST, gene by gene.
To do this, 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 creates two output files that are required for the subsequent method calibration steps:
calibration/preparations/calibration_blast.tsvcalibration/preparations/calibration_prepared.rds
Method calibration is performed in three steps. 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.
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.
Adjusting thresholds
- Thresholds are always defined as a sequence (seq) of values between a minimum value and maximum value, with the specified step size
- E.g.,
seq(50,100,1)creates thresholds from 50 to 100 in steps of 1 - Thresholds can be adjusted for each variable as desired by changing the minimum, maximum, or step size
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 optimal thresholds for the dataset need to be manually specified for each variable in the file:
calibration/manual_input_needed/calibration_alignments.tsv.
Guidelines:
- To decide on 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.
- Thresholds established in the main manuscript across rattans and mahoganies are specified to provide guidance.
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.
This step first calculates gene performance for each gene, i.e. the percentage of samples correctly identified to species.
Second, the script explores the impact of different minimum thresholds of gene performance on overall accuracy of identification and retrievability.
Based on the outputs produced, the optimal gene performance threshold for the dataset needs to be specified.
This step assesses the impact of different minimum thresholds of parliament size, i.e. the number of genes in the Gene Parliament, on overall accuracy of identification and retrievability.
Based on the outputs produced, the optimal parliament size threshold for the dataset needs to be specified.
Based on the optimal thresholds and parameters specified in the script, the following three calibration files are produced as the last step of the script:
Gene performance table calibration_gene_performance.csv:
This file contains in the first column the gene name and in the second column the performance of the gene, i.e. the percentage of test samples that were correctly identified to species.
Example file:
| gene | performance |
|---|---|
| Gene1 | 56.67 |
| Gene2 | 65.22 |
| Gene3 | 65.11 |
| Gene4 | 85 |
| Gene5 | 47.56 |
Filtering thresholds table calibration_filtering_thresholds.csv:
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.
Example 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 |
Confidence estimates table calibration_confidence_support.csv:
This file contains the probability of the top identification being correct, close or wrong, depending on the percentage of genes supporting the identification.
Example file:
| range_support | probability_correct | probability_close | probability_wrong |
|---|---|---|---|
| [0,20] | 0 | 0 | 100 |
| (20,40] | 78.57 | 21.43 | 0 |
| (40,60] | 82.14 | 17.86 | 0 |
| (60,80] | 100 | 0 | 0 |
| (80,100] | 100 | 0 | 0 |
These calibration files can be directly used as input in the GeneParliamentID pipeline.
It is possible to run GeneParliamentID without conducting method calibration by providing "dummy" calibration files with maximally lenient thresholds that pass all filtering steps.
WARNING: This will most likely result in considerably reduced accuracy of identification.
Doing this may be justified e.g. if a test dataset is not available or when conducting a first explorative analysis.
Dummy gene performance table:
To include all genes in the analyses, set performance to 100 for all genes.
| gene | performance |
|---|---|
| Gene1 | 100 |
| Gene2 | 100 |
| Gene3 | 100 |
| ... | 100 |
| Gene353 | 100 |
Dummy filtering thresholds table:
To disable all filtering, set minimum thresholds to 0 and maximum thresholds to 99999:
| min_similarity | min_length | max_gapopens | max_mismatches | max_evalue | min_bitscore | min_gene_performance | min_parliament_size |
|---|---|---|---|---|---|---|---|
| 0 | 0 | 99999 | 99999 | 99999 | 0 | 0 | 0 |
Dummy confidence estimates table:
To include dummy confidence estimates, specify a single bin from 0 to 100, with NA for all categories:
| range_support | probability_correct | probability_close | probability_wrong |
|---|---|---|---|
| [0,100] | NA | NA | NA |
Upon completion of Method calibration, continue to 3. Method validation.