ICPPB 2024 Workshop:

ICPPB 2024 Workshop: === ###### (Meta)Genomics for pathogen identification Metagenome classification using custom taxonomy databases created with LINtax === ###### Parul Sharma ###### July 2024 In this tutorial we will cover: [TOC] ## LINtax LINtax creates LIN (Life identification number) based taxonomy files for any given list of input genomes. Instead of traditional ranks, LIN-taxonomy is based on 20 thresholds/ranks of average nucleotide identity (ANI). These thresholds of similarity ranges from 70%-99.999%. The current default implementation of LIN thresholds are outlined below: ![image](https://hackmd.io/_uploads/r14M0HUwR.png) Traditional taxonomy considers 'species' as the lowest rank which corresponds to approximately 95% ANI. However, with LIN-taxonomy you get much more resolution and 14 more ranks at the sub-species level which allows for accurate outbreak level detection of pathogens. #### Creating custom taxonomy with LINtax Metagenome classifiers use taxonomy to improve classifications by reducing errors due to genome similarity. All taxonomic classifiers (metagenome classifiers that use taxonomy) make use of the NCBI Taxonomy database which retains the information in 2 files: `names.dmp` : stores all tax ranks and their names `nodes.dmp` : stores all the heirarchial relations between different ranks Using LINtax, users can create this taxonomy file structure for their list of input. Example input-file looks omething like this: | LIN | Species | Strain | FileName | Accession | |-----------------------------------|----------------------|-----------|--------------------------------------------------|-----------------| | 14,1,0,0,0,0,0,0,0,0,6,0,1,0,1,0,0,0,0,0 | Ralstoniasolanacearum | OE1_1 | GCF_001879565.1_ASM187956v1_genomic.fna | GCF_001879565.1 | | 14,1,0,0,0,0,0,0,0,0,6,0,1,0,0,0,0,0,0,0 | Ralstoniasolanacearum | PSS1308 | GCF_001870805.1_ASM187080v1_genomic.fna | GCF_001870805.1 | | 14,1,0,0,0,0,0,0,0,0,2,1,0,0,0,0,0,0,0,0 | Ralstoniasolanacearum | FJAT_1458 | GCF_001887535.1_ASM188753v1_genomic.fna | GCF_001887535.1 | | 14,1,0,0,0,0,0,0,0,0,2,0,0,4,4,0,0,0,0,0 | Ralstoniasolanacearum | Pe_13 | GCF_012062595.1_ASM1206259v1_genomic.fna | GCF_012062595.1 | The current implementation of LINtax works with [Kraken2](https://github.com/DerrickWood/kraken2/wiki/Manual) and [krakenuniq](https://github.com/fbreitwieser/krakenuniq). #### Workflow for analysis ```mermaid graph LR create-taxonomy --> build-db --> classify --> report ``` [Click here](https://github.com/parul-sharma/LINtax) for more information. ## Ralstonia Case study :::info Can we detect the pathogen in a metagenome? Can we identify the correct strain? ::: We will use a pre-computed database to analyze 3 plant metagenome samples. The pre-computed database was build as a [kraken2](https://github.com/DerrickWood/kraken2/wiki/Manual) database using LIN-taxonomy structure with [lintax](https://github.com/parul-sharma/LINtax). To set-up everything for this tutorial: 1. Download the database using: - Kraken2 database ``` wget -O kraken2-db.tar.gz https://osf.io/u4tq8/download ``` - Extract the contents of the file : ``` tar xvf kraken2-db.tar.gz ``` 2. Download the rest of set-up : ``` git clone https://github.com/parul-sharma/LINtax.git ``` - Download lingroups file ``` wget -O LINgroups.txt https://osf.io/download/nqms2/ ``` - Download the samples: download Sample-0 wget -O Sample-0.fasta https://osf.io/r72cv/download - download Sample-II ``` wget -O Sample-II.fasta https://osf.io/xvytc/download ``` - download Sample-IV ``` wget -O Sample-IV.fasta https://osf.io/g4krs/download ``` With everything set-up. We can now analyze our samples. ### Lets start with sample-0 since we have the pre-computed database, we can directly use the 'classify' function to classify the sample (We can also directly use kraken2 commands to classify and then use lintax report to get the final report): :::info Sample-0: ::: Let's see what the input requirements are: ``` lintax classify -h ``` ``` ---------------------------------------------------------------------------------------- LINtax: Create custom taxonomy for metagenome classification ---------------------------------------------------------------------------------------- Usage: lintax classify [options] <input_files> Options: --db LINtax Database to use for classifications (make sure to run build before this step) [REQUIRED] --lins Text file containing lingroups (custom taxa) and their LIN prefixes [REQUIRED] -o or --output LINreport containing the results [REQUIRED] -h or --help Show usage Other options: -c or --confidence Specify the confidence threshold for kraken classification [optional] --paired For paired-end reads Example: for single reads - lintax classify --db /database/ --lins linfile.txt -o sample-results.tsv <sample.fastq> for paired reads - lintax classify --db /database/ --lins linfile.txt -o sample-results.tsv --paired <sample_1.fastq> <sample_2.fastq> ``` This tells us that we need the 'db', 'lins' and sepcify the output file. Let's analyze our sample: ``` lintax classify --db kraken2-db --lins LINgroups.txt -o sample0.LINreport -c 0.1 Sample-0.fasta ``` ``` Loading database information... done. 9204 sequences (37.44 Mbp) processed in 2.539s (217.5 Kseq/m, 884.81 Mbp/m). 0 sequences classified (0.00%) 9204 sequences unclassified (100.00%) ``` The above command shows that 0 sequences in our sample were classified. This is because Sample-0 is our control sample without any pathogenic reads! :::success Sample 0 : Control : No pathogenic reads detected! ::: This is what the final LINreport for this sample looks like (with no classification): ``` LINgroup_Name LINgroup_prefix Assigned_reads Percentage_assigned_reads Unique_Assigned_reads Percentage_unique_assigned_reads Total_reads_length A_Total_reads;B_PhylI 14,1,0,0,0,0,0,0,0,0 A_Total_reads;B_PhylI;C_seq14 14,1,0,0,0,0,0,0,0,0,3 A_Total_reads;B_PhylI;C_seq15 14,1,0,0,0,0,0,0,0,0,2 A_Total_reads;B_PhylI;C_seq34 14,1,0,0,0,0,0,0,0,0,6 A_Total_reads;B_PhylII 14,1,0,0,0,3,0 A_Total_reads;B_PhylII;C_IIA 14,1,0,0,0,3,0,1 A_Total_reads;B_PhylII;C_IIC 14,1,0,0,0,3,0,2 A_Total_reads;B_PhylII;C_IIB 14,1,0,0,0,3,0,0 A_Total_reads;B_PhylII;C_IIB;D_seq1&seq2 14,1,0,0,0,3,0,0,0,0,1,0,0,0,0 ``` :::info Adding -c (or --confidence) option is important! A '0.1' score equates to 10% confidence threshold and tells kraken2 to only accept those taxonomic ranks where this confidence threshold is met. [Read more about confidence thresholds] (https://github.com/DerrickWood/kraken2/wiki/Manual#confidence-scoring) For LINtax, in general, recommended confidence threshold is 0.1 for long reads and 0.3-0.4 for short reads. ::: ### Now let's look at `sampleII` We can run the same classify command to run sampleII with -c 0.1. ``` lintax classify --db kraken2-db --lins LINgroups.txt -o sampleII.LINreport -c 0.1 Sample-II.fasta ``` ``` Loading database information... done. 9204 sequences (37.73 Mbp) processed in 3.188s (173.2 Kseq/m, 710.08 Mbp/m). 812 sequences classified (8.82%) 8392 sequences unclassified (91.18%) ``` The output indicates that 812 reads (8.2%) were classified. Let's look at the LINreport to see what taxonomcic lineages were identified. ###### using `awk '$3 != 0'` to only select for those taxa where assigned reads were more than 0. ``` cat sampleII.LINreport | awk '$3 != 0' ``` ``` LINgroup_Name LINgroup_prefix Assigned_reads Percentage_assigned_reads Unique_Assigned_reads Percentage_unique_assigned_reads Total_reads_length A_Total_reads;B_PhylII 14,1,0,0,0,3,0 732 90.14778325123153 598 73.64532019704434 5557909 A_Total_reads;B_PhylII;C_IIB 14,1,0,0,0,3,0,0 134 16.502463054187192 79 9.729064039408868 708286 A_Total_reads;B_PhylII;C_IIB;D_seq1&seq2 14,1,0,0,0,3,0,0,0,0,1,0,0,0,0 16 1.9704433497536946 14 1.7241379310344827 115923 A_Total_reads;B_PhylII;C_IIB;D_seq1&seq2;E_seq1 14,1,0,0,0,3,0,0,0,0,1,0,0,0,0,0,0 2 0.24630541871921183 2 0.24630541871921183 13897 Total_reads 812 ``` We see that all assigned reads follow the same taxonomic lineage-- leading to Phylotype II sequevar1 !!! This is our infected sample with the Select Agent pathogen! :::success Sample-II: Pathogenic sample; Select Agent Ientified! ::: ### Now let's look at sampleIV ``` lintax classify --db kraken2-db --lins LINgroups.txt -o sampleIV.LINreport -c 0.1 Sample-IV.fasta ``` ``` Loading database information... done. 9204 sequences (37.91 Mbp) processed in 3.206s (172.3 Kseq/m, 709.57 Mbp/m). 789 sequences classified (8.57%) 8415 sequences unclassified (91.43%) ``` We see that 789 reads (8.57% ot total) were classified for this sample. Let's look at the LINreport to see what taxonomcic lineages were identified. ``` cat sampleIV.LINreport | awk '$3!=0' ``` ``` LINgroup_Name LINgroup_prefix Assigned_reads Percentage_assigned_reads Unique_Assigned_reads Percentage_unique_assigned_reads Total_reads_length A_Total_reads;B_PhylIV 14,1,0,0,0,2,0,0,0 585 74.1444866920152 483 61.21673003802282 4654046 A_Total_reads;B_PhylIV;C_seq10 14,1,0,0,0,2,0,0,0,0,0,0 102 12.927756653992395 74 9.37896070975919 680687 Total_reads 789 ``` We find in this sample, all assigned reads follow the taxonomic lineage of Phylotype IV, sequevar 10 group --> This is not a Select agent pathogen. :::success Sample-IV: Pathogenic sample, but not Select Agent. :::

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.