2. Method calibration
This page walks you through the process of selecting the optimal parameters for the GeneParliamentID pipeline for any lineage, step by step.
For method calibration, you require a test dataset (see below) and a reference dataset.
Method calibration requires a test dataset of samples with known species identity.
Gene sequences of the test samples can be retrieved from targeted sequence capture data using, for example, HybPiper or CAPTUS.
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.
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
- Estimate confidence of identification depending on the percentage of genes supporting the top identification
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 allow full flexibility in selecting parameters, method calibration is conducted as an interactive R script calibration.R.
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
Directory containing multiple unaligned genes, each comprising all retrieved sequences for the reference samples.
The structure and requirements are the same as for the test dataset.
Requirements:
- Gene sequences need to be unaligned and in
.FNAformat. - Gene names need to match the gene names in the test dataset.
- Each reference 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 reference directory
Gene1.FNA:
>Entandrophragma_angolense_CQL_22
CTCCCATATAGGGGTGCTTGGTTGTGGGTGGGTTCTGAGATGATTCATTTAATGGAACTT
CCAAATCCAGATCCCTTAACTGGACGACCGGCACATGGTGGTCGAGATCGTCATACTTGT
ATTGCGATTCGAGATGTGTCTAAACTCAAAATAATCCTTGATAAAGCAGGTATTCCTTAC
>Entandrophragma_bussei_CQL_EB42
CTCCCATATAGGGGTGCTTGGTTGTGGGTGGGTTCTGAGATGATTCATTTAATGGAACTT
CCAAATCCAGACCCCTTAACTGGACGACCGGCACATGGTGGTCGAGATCGTCATACTTGT
ATTGCGATTCGAGATGTGTCTAAACTCAAAATAATCCTTGATAAAGCAGGTATTCCTTAC
>Entandrophragma_candollei_CQL_13
CTTCCGTATAGGGGTGCTTGGTTGTGGGTGGGTTCTGAGATGATTCATTTAATGGAACTT
CCAAATCCAGACCCCTTAACGGGACGACCGGCACATGGTGGTCGAGATCGTCATACTTGT
ATTGCGATTCGAGATGTGTCTAAACTGAAAATAATCCTTGACAAAGCAGGTGTTCCTTAC
Gene2.FNA:
>Entandrophragma_angolense_CQL_22
GTTGGGGATCTACACGGAGACCTTGATCAAGCGAGATGTGCACTTGAGATTGCTGGTGTC
>Entandrophragma_bussei_CQL_EB42
GTTGGGGATCTACACGGAGACCTTGATCAAGCGAGATGTGCACTTGAGATTGCTGGTGTC
>Entandrophragma_candollei_CQL_13
GCTGGGGATCTACATGGAGACCTTGATCAAGCAAGATGTGCACTTGAGATTGCTGGTGTC
Gene3.FNA:
>Entandrophragma_angolense_CQL_22
GAATCTCAAGGATTAACATCTCGACGATTGTGTCTTACATGTATATGTTCAACTCTAGCT
CTGATTAACAGTTCCGGCACGTTGGTTTCTGTACAAAAGGCAATTGCTTTGGAAGGAAAA
>Entandrophragma_bussei_CQL_EB42
GATATGTGTGGTGGTACAGGAAAATGGAAAGCTCTCAACAGAAAACGTGCTAAAGATGTT
TACGAGTTTACAGAATGTCCAAATTGTTATGGTCGTGGGAAACTTGTGTGTCCGGTTTGC
>Entandrophragma_candollei_CQL_13
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.
This step calculates the percentage of test samples with correct, close or wrong identification depending on the percentage of genes supporting the identification.
First, the top identification for each sample is retrieved.
Second, the top identifications are assigned to "bins" depending on the percentage of genes supporting the top identification, and for each bin the percentage of samples with correct, close and wrong identifications is calculated.
- The range of support is first divided into equal-sized bins
- The top identifications are assigned to bins depending on the percentage of genes supporting the identification
- For each bin, the percentage of correct, close and wrong identifications is estimated
By default, confidence estimates are produced for between 1 and 10 bins (the number of bins can be adjusted freely). For example, if the number of bins is set to 5, five categories of0-20%,>20-40%,>40-60%,>60-80%and>80-100%are created. A top identification supported by 37% of genes would be assigned the category>20-40%.
For each number of bins, a plot showing support for each bin is produced.
Based on these results, the optimal number of bins for the dataset needs to be specified.
Guidelines:
- Each bin should contain multiple samples to allow an estimate of the accuracy of identification for the bin
- There is a trade-off between too few bins and too many bins
- Too few bins: loss of resolution
- Too many bins: some bins may be unrepresentative due to the low number of samples, or even contain no samples at all
- Increasing the number of test samples may allow to also increase the number of bins and therefore the granularity of confidence assessments
Example:
The below figures shows confidence estimates of the same test data depending on whether support is divided into 3, 5 or 7 bins. Note that there are very few samples falling into the first (leftmost) bin.
The figure with 5 bins includes samples in all bins and contains more detailed information than the figure with 3 bins. While the figure with 7 bins contains even more detailed information, there are no samples in the first bin, and confidence can therefore not be estimated for this bin.
In this case, estimating confidence using 5 bins might be the preferred choice for a maximum granularity in the confidence assessments whilst being able to assess confidence for all bins.
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 |