6. Tutorial
This tutorial guides you through sample identification using GPID with a test dataset, step by step.
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
Alternatively, you can download it manually by clicking this link.
Extract example files using:
tar -zxvf example_data.tar.gz
If you downloaded the file to your local machine, you can also manually extract 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 the same 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
Contains a folder for the test sample CQL_2, which in turn contains .FNA files with gene sequences for this sample.
Example contents of two genes, ANGIO353g4527 and ANGIO353g4989:
File ANGIO353g4527.FNA:
>CQL_2
TTGGTTATTGGTGGAGGGGGAAGGGAACATGCATGCTATGCTTTGAAGCGATCTCCCTCA
TGTGATGCTGTATTTTGTGCTCCCGGCAATGCGGGGATATCCAGCTCAGGGGATGCAACT
TGTGTCACAGACTTGGACATTTTAGATGGGGAAGCTGTGATCTCCTTCTGCCGCAAGTGG
File ANGIO353g4989.FNA:
>CQL_2
GGGAAAGCTCGTAGTGACATGACTGTGGAGGAAAGAATTGGTGCAGCAGACCGCAGAAAG
ATTGACGGAAATGCCTTCTTTAAGGAGGAGAAACTGGAAGAGGCCATGCAGCAGTATGAA
This folder contains the three required calibration files.
All calibration files were produced using the calibration script and a set of test samples of known identity, see Method calibration.
File calibration_gene_performance.csv:
This file contains the gene performance for each gene (id_correct_pct), i.e. the percentage of test samples that were correctly identified to species.
| gene | id_correct_pct |
|---|---|
| ANGIO353g4471 | 56.67 |
| ANGIO353g4527 | 65.22 |
| ANGIO353g4691 | 65.11 |
| ANGIO353g4724 | 85 |
| ANGIO353g4744 | 47.56 |
File calibration_filtering_thresholds.csv:
This file contains the thresholds that were identified as optimal for the given lineage.
| simimilarity | length | gap | mismatch | evalue | bitscore | gene_performance | parliament_size |
|---|---|---|---|---|---|---|---|
| 98 | 100 | 1 | 5 | 1e-60 | 200 | 30 | 10 |
File 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.
For example, an identification supported by between 40 and 60% of genes has a probability of 82.14% of being correct (correct to species level), 17.86% of being close (correct to a user-defined group of closely related species), and 0% of being wrong (neither correct nor close).
| support | correct | close | 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 |
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 |
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 dataset, you need to perform method calibration. 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 |
This only needs to be done once.
To run the GeneParliamentID pipeline using the example data, execute the following command from the example_data directory:
gpid -i samples/CQL_2/ -r reference/ -g calibration/calibration_gene_performance.csv -t calibration/calibration_filtering_thresholds.csv -c calibration/calibration_confidence_support.csv -s 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 (in different formats) with top 10 identifications:
CQL_2_gpid.jpgCQL_2_gpid.pdfCQL_2_gpid.svg
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.