6. Tutorial
This tutorial guides you through sample identification using GPID with a test dataset, step by step.
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
Conduct Method calibration
Before running GeneParliamentID for the first time, we strongly recommend conducting Method calibration to increase and assess the accuracy of identification.
This only needs to be conducted once for any lineage or dataset.
Run GeneParliamentID pipeline
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.