6. Tutorial
This tutorial guides you through sample identification using GPID with a test dataset, step by step:
- Install GPID
- Get example data
- Build reference
- Perform method calibration
- Conduct method validation
- Identify sample
- Interpret results
For details on Reference construction, Method calibration, Method validation, Sample identification and Interpretation of the results, please see the respective other Wiki pages.
We recommend installation of GPID including all dependencies using conda with a new environment:
conda create --name gpid gpid
To activate the GPID conda environment, use:
conda activate gpid
To confirm that the installation has worked and show a help message on how to use GPID, simply run:
gpid
To download the example data folder, run:
wget https://github.com/BenKuhnhaeuser/GPID/blob/main/example_data.tar.gz
Then, extract the files in the folder using:
tar -zxvf example_data.tar.gz
Alternatively, you can download the folder to your machine by clicking this link and extract its contents by double-clicking on the file.
The extracted directory contains the following four folders:
Contains .FNA files for the same genes, each with reference sequences for multiple species.
The header line starts with <Genus>_<species>.
Example content of two genes, ANGIO353g4527 and ANGIO353g4989:
File ANGIO353g4527.FNA:
>Entandrophragma_angolense_CQL_22
AGGGTGGTTGTGTTGGTTATTGGTGGAGGGGGGAGGGAACATGCACTTTGTTATGCTTTG
AAGCGATCTTCCTCATGTGATGCTGTATTTTGTGCTCCTGGAAATGCGGGGATATCCAGC
TCAGGGGATGCAACTTGTATCACGGACCTAGACATTTTAGATGGGGAAGCTGTGATCTCC
>Entandrophragma_bussei_CQL_EB42
GAGAGGGTGGTTGTGTTGGTTATTGGTGGAGGGGGGAGGGAACATGCACTTTGCTATGCT
TTGAAGCGATCTCCCTCATGTGATGCTGTATTTTGTGCTCCCGGAAATGCGGGGATATCC
AGTTCAGGGGATGCAACTTGTATCATGGACCTAGACATTTTAGACGGGGAAGCTGTGATC
>Entandrophragma_candollei_CQL_13
AGGGTGGTTGTGTTGGTTATTGGTGGAGGGGGGAGGGAACATGCACTTTGCTATGCTTTG
AAGCGATCTCCCTCATGTGATGCTGTATTTTGTGCTCCTGGAAATGCGGGGATATCCAGC
TCAGGGGATGCAACTTGTATCACGGACCTAGACACTTTAGACGGGGAAGCTGTGATCTCC
File ANGIO353g4989.FNA:
>Entandrophragma_angolense_CQL_22
AAAGAAATGACTGGCTTGGCTATTGGTGTCTCAAGCATGAAGTCTGGTGAACATGCCCTG
TTACATGTGGGCTGGGAATTGGGTTATGGGAAAGAAGGAAGCTTCTCTTTCCCAAATGTG
>Entandrophragma_bussei_CQL_EB42
TTCCCAAATGTGCCTCCTATGGCAGACTTATTATACGAGGTTGTGCTTATTGGTTTTGAT
GAAACCAAAGAAGGGAAAGCTCGTAGCGACATGACTGTAGAGGAAAGGATTGGTGCAGCA
>Entandrophragma_candollei_CQL_13
AAAGAAATGACTGGCTTGGCTATTGGTGTCTCAAGCATGAAGTCTGGTGAACATGCCCTG
TTACATGTGGGCTGGGAATTGGGTTATGGGAAAGAAGGAAGCTTTTCTTTCCCAAATGTG
Folder with samples used for method calibration.
Contains fasta files for the same genes, each with reference sequences for multiple species. The header line starts with <Genus>_<species>.
Formatting is the same as the reference folder.
Folder with samples used for method validation.
Contains fasta files for the same genes, each with reference sequences for multiple species. The header line starts with <Genus>_<species>.
Formatting is the same as the reference folder.
This optional file specifies for each species a user-defined group of closely related species. This could be at any taxonomic level, for example, family, genus, or based on an infrageneric classification. In this case, we specified the genera as our species groups.
File species_groups.csv:
| genus_species | species_group |
|---|---|
| Entandrophragma_angolense | Entandrophragma |
| Entandrophragma_congoense | Entandrophragma |
| Entandrophragma_bussei | Entandrophragma |
| Khaya_agboensis | Khaya |
| Khaya_nyasica | Khaya |
Contains one folder for each sample, which in turn contains fasta files with gene sequences for this sample.
Example contents of sample folder CQL_2, containing two genes, ANGIO353g4527 and ANGIO353g4989:
File ANGIO353g4527.FNA:
>CQL_2
TTGGTTATTGGTGGAGGGGGAAGGGAACATGCATGCTATGCTTTGAAGCGATCTCCCTCA
TGTGATGCTGTATTTTGTGCTCCCGGCAATGCGGGGATATCCAGCTCAGGGGATGCAACT
TGTGTCACAGACTTGGACATTTTAGATGGGGAAGCTGTGATCTCCTTCTGCCGCAAGTGG
File ANGIO353g4989.FNA:
>CQL_2
GGGAAAGCTCGTAGTGACATGACTGTGGAGGAAAGAATTGGTGCAGCAGACCGCAGAAAG
ATTGACGGAAATGCCTTCTTTAAGGAGGAGAAACTGGAAGAGGCCATGCAGCAGTATGAA
To construct a BLAST database for each gene in the reference directory, run:
gpid reference -r example_data/reference
This command performs the following steps:
- locating FASTA files (.FNA, .fasta, .fa)
- validating FASTA header format
- checking whether BLAST databases already exist
- building missing BLAST databases with
makeblastdb
This only needs to be done once.
To optimise the pipeline parameters for your particular lineage and genes, you need to perform method calibration using a dataset with samples of known identity, the Calibration dataset. This only needs to be done once. Method calibration is done using the following five commands:
To prepare the calibration dataset for method calibration, run:
gpid calibrate prepare -r example_data/reference -i example_data/calibration
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.
To conduct method calibration for the alignment filtering thresholds, run:
gpid calibrate alignments
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/.
Download the files to your local machine and look at them
The produced figures should look as follows:
Based on the figures produced, you need to identify the optimal values for the dataset. Optimal values should increase accuracy as much as possible whilst decreasing retrievability as little as possible. In other words, you are trying to find the value that maximises both accuracy and retrievability.
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.
To open and edit the file, use a word editor such as vim or nano:
vim calibration/manual_input_needed/calibration_alignments.tsv
In this case, we select the following values as our optimal thresholds:
| parameter | value |
|---|---|
| min_similarity | 98 |
| min_length | 100 |
| max_gapopens | 1 |
| max_mismatches | 5 |
| max_evalue | 1d-60 |
| min_bitscore | 200 |
Save and close calibration file.
To run method calibration for the gene performance thresholds, run:
gpid calibrate genes -a calibration/manual_input_needed/calibration_alignments.tsv
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:
vim calibration/manual_input_needed/calibration_genes.tsv
In this case, we select the a gene performance of 30 as our optimal threshold:
| parameter | value |
|---|---|
| min_gene_performance | 30 |
To perform method calibration for parliament size, i.e. the number of genes in the Gene Parliament, run:
gpid calibrate parliament -a calibration/manual_input_needed/calibration_alignments.tsv -g calibration/manual_input_needed/calibration_genes.tsv
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:
vim calibration/manual_input_needed/calibration_parliament.tsv
In this case, we select the a minimum parliament size of 10 genes as our optimal threshold:
| parameter | value |
|---|---|
| min_parliament_size | 10 |
To combine the optimal thresholds selected during method calibration into a single calibration file, run:
gpid calibrate combine -a calibration/manual_input_needed/calibration_alignments.tsv -g calibration/manual_input_needed/calibration_genes.tsv -p 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
The resulting file should look like this:
| 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 |
To assess the accuracy of identification depending on the percentage of genes supporting the identification, you need to provide a validation dataset with samples of known identity, the Validation dataset. The validation dataset should be independent of the calibration dataset to avoid over-fitting. This only needs to be done once.
Method validation is performed in the following steps.
To prepare the validation dataset for method validation, run:
gpid validate prepare -r example_data/reference -i example_data/validation
This matches each sample in the validation 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 validation analyses.
To assess the accuracy of identification depending on support for the identification, run:
gpid validate confidence -g calibration/calibration_gene_performance.csv -t calibration/calibration_filtering_thresholds.csv
This calculates the percentage of test samples with correct, close or wrong identification depending on the percentage of genes supporting the identification, using the gene performance and filtering thresholds files produced during method calibration as inputs. Method validation is conducted in the following steps:
- The top identification for each sample and the percentage of genes supporting this identification are retrieved.
- The percentage of support is divided into equal-sized "bins", and each identification is assigned to one of these bins depending on the percentage of genes supporting the top identification. For example, if the number of bins is set to 5, five categories of
0-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 bin, the percentage of samples with correct, close and wrong identifications is estimated. If a bin contains no samples, the percentages are given as
NA.
Confidence estimates are produced for between 1 and 10 bins. For each number of bins, a plot showing support for each bin is produced. These plots are summarised in a single file:
validation/tests/validate_confidence.pdf
To inspect the plots, download the file validation/tests/validate_confidence.pdf to your local machine. Each plot should look like this:
Based on inspecting the results in validation/tests/validate_confidence.pdf, the optimal number of bins for the dataset needs to be specified.
Guidelines:
- There should be no bins that contain no samples, as the accuracy of identification cannot be assessed in this case.
- 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
In this case, we select five bins as our optimal value. To do this, run:
gpid validate bins -b 5
This command produces a table validation/validation_confidence_support.csv. It details the probability of the top identification being correct, close or wrong depending on the percentage of genes supporting the identification. This file is used as input for sample identification.
To identify a sample (in this case CQL_2) using GPID, run:
gpid identify -i example_data/samples/CQL_2/ -r example_data/reference/ -g calibration/calibration_gene_performance.csv -t calibration/calibration_filtering_thresholds.csv -c calibration/calibration_confidence_support.csv -s example_data/species_groups/species_groups.csv
This will generate the following files in the current directory:
- Gene Parliament table with all identifications:
CQL_2_gpid.csv - Gene Parliament figure with top 10 identifications:
CQL_2_gpid.pdf
To identify a different sample, you would only need to change the sample name: -i example_data/samples/CQL_2/.
The Gene Parliament figure gives a quick overview of the top 10 identifications that were retrieved and their relative support.
In this case, a clear majority of genes (45.1%) support the identification as Entandrophragma angolense, whilst Entandrophragma excelsum (18.85%) and Entandrophragma congoense (16.39%) also get sizeable support. Other species have almost negligible support, but all belong to the same genus Entandrophragma:
To see all identifications that were retrieved, we can have a look at the Gene Parliament table. Importantly, the table contains for the top identification information on the probability of the identification being correct (correct to species level), close (correct to species group) or wrong (neither correct to species nor to species group). This probability is based on the percentage of genes supporting the top identification (Support_pct) and was obtained from the file calibration_confidence_support.csv, which was generated during method calibration using a test dataset.
In this case, the table indicates a probability of 82.14% that the identification is correct to species level, and a further 17.86% (thus totaling 100%) that the identification is close, i.e. correct to species group (in this case genus). The probability that the identification is wrong, i.e. neither correct nor close, is estimated to be 0. Overall, we can be fully confident that the sample belongs to the genus Entandrophragma, and have high confidence that it was taken from the species Entandrophragma_angolense.
| Sample | Rank | Identification | Species_group | Support_pct | Support_count | Parliament_size | Data_checks | ID_correct_pct | ID_close_pct | ID_wrong_pct |
|---|---|---|---|---|---|---|---|---|---|---|
| CQL_2 | 1 | Entandrophragma_angolense | Entandrophragma | 45.08 | 55 | 122 | PASSED | 82.14 | 17.86 | 0 |
| CQL_2 | 2 | Entandrophragma_excelsum | Entandrophragma | 18.85 | 23 | |||||
| CQL_2 | 3 | Entandrophragma_congoense | Entandrophragma | 16.39 | 20 | |||||
| ... |
For more details on how to interpret the Gene Parliament and different scenarios you might encounter, see Interpretation.