2. Method calibration
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 the GeneParliamentID pipeline with 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 thresholds and parameters in four steps (see below), one at a time:
- 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.
requires a calibration dataset of samples with known species identity.
Guidelines for assembling a good test dataset:
- Test samples should be expert-identified, e.g. by taxonomists.
- Voucher specimens of the test 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 test 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.
Directory containing multiple unaligned genes, each comprising all retrieved sequences for the test samples.
Requirements:
- 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.
Example contents of test sample directory
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 test samples against the reference dataset using BLAST, gene by gene.
Enable conda environment
conda activate gpid
Specify path to test sample directory (CHANGE TO YOUR PATH)
test_dir=<PATH/TO/TEST_DIR>
Specify path to reference directory (CHANGE TO YOUR PATH)
ref_dir=<PATH/TO/REF_DIR>
Provide list with all genes
genelist=genelist.txt
- The file
genelist.txtneeds to contain the names of all genes that were retrieved, one per line. For example:Gene1 Gene2 Gene3
The following script matches the test samples against the reference samples for each gene in the gene list. It then retrieves the best match per gene and sample, and combines them in a single file blast.tsv that serves as input for method calibration.
for gene in `cat $genelist`; do
# Make BLAST database
makeblastdb -in "$ref_dir"/"$gene".FNA -parse_seqids -dbtype nucl
# Match query samples against BLAST database
blastn -query "$test_dir"/"$gene".FNA -db "$ref_dir"/"$gene".FNA -task megablast -outfmt "6 qseqid sseqid pident length mismatch gapopen evalue bitscore" -max_target_seqs 1000000 |
# Sort by query (k1,1) and then by bit score (k8,8), descending and numerical (rn), shuffling identical reference queries (k2,2R)
sort -t $'\t' -k1,1 -k8,8rn -k2,2R |
# Get first match per query (which is the one with the highest bit score)
awk '!seen[$1]++' |
# Add gene name as the first column
awk -v genename="$gene" '{print genename "\t" $0}' ; done |
# Add header row and combine individual searches in a single file
sed '1i gene\tquery\ttarget\tpident\tlength\tmismatch\tgapopen\tevalue\tbitscore' > blast.tsv
Method calibration is performed using the interactive script calibration.R, which contains detailed step-by-step instructions. 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.
The script requires the file blast.tsv as an input. Optionally, a file specifying for each species a species group can be provided.
Places in the script where user input is required are highlighted using INPUT REQUIRED, the remaining script should not be changed.
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
In step 1, the script 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
Based on the figures produced, the optimal thresholds for the dataset need to be specified for each variable.
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 |