ORI (Oxford nanopore Reads Identification) is a software using long nanopore reads to identify bacteria present in a sample at the strain level.

There are two sub-parts in the ORI program: (1) the creation of the index containing the reference genomes of the interest species and (2) the query of this index with long reads from Nanopore sequencing in order to identify the strain(s).

The index is based on the structure implemented in HowDeSBT [1] modified in order to use qgrams (word from spaced seeds) instead of kmers.

As previously said, we replaced kmers by qgrams. To get qgrams we use spaced seeds which introduce don’t care positions in kmers that won’t be disturbed by sequencing errors. To select the best seed pattern for classification of long reads we used the iedera software [2][3]. The best seed for our classification tools seems to be the following: size 15 and weight 13, 111111001111111. This seed is then applied to words of size 15 resulting in qgrams used instead of kmers.

The preconstructed indexes for Streptococcus thermophilus strains are available in the ORI github in the directory: preconstructed_indexes directory of ORI.

_{1. Robert S Harris and Paul Medvedev, Improved representation of sequence bloom trees, Bioinformatics, btz662 <10.1093/bioinformatics/btz662>}

_{2. Noe L., Best hits of 11110110111: model-free selection and parameter-free sensitivity calculation of spaced seeds, Algorithms for Molecular Biology, 12(1). 2017 http://doi.org/10.1186/s13015-017-0092-1}

_{3. Kucherov G., Noe L., Roytberg, M., A unifying framework for seed sensitivity and its application to subset seeds, Journal of Bioinformatics and Computational Biology, 4(2):553-569, 2006 http://doi.org/10.1142/S0219720006001977}

The easiest way to install ORI is through conda.

conda create -p ori_env
conda activate ori_env
conda install -c gsiekaniec -c conda-forge ori

The seedfile.txt can be found in the seed directory of ORI.

You can also create the seedfile.txt in your chosen repertory:

cd path/to/repertory
touch seedfile.txt
echo "111111001111111" > seedfile.txt

⚠️ Caution: reference genomes must be in .fasta or .fna and in the same directory.

If you have the fasta files distributed in several subfolders, you have to redirect them to a single directory.

Repertory containing fasta files (must be the current directory) to run the following scripts:

cd path/to/the/fasta/files/repertory

As the last quantification step requires the size of the genomes, it is preferable to calculate it when we have our genomes:

ORI.py length -g path/to/the/genomes -o length.txt

The ORI.py length step also allows to calculate an effective size for the bloom filters. Use the -seed_size and -false_positive_rate options . The bloom filter size is given in the bf_min_size.txt file.

Then we create the bloom filters for all genomes:

howdesbt makebfQ --k=15 --qgram=path/to/seedfile.txt --bits=0.25G *.fasta

Then we get the names of the bf (bloom filter) files used to create the tree:

ls *.bf > leafname

Now that the bloom filters are created it is no longer necessary to keep the fastas files. If it is not necessary to keep them, they can be deleted to save space. In addition, if the fasta files cannot be completely downloaded on the machine due to lack of space, it is possible to download them little by little and create the filters as you go by deleting the fasta files once in the form of a filter (.bf).

howdesbt distance --list=leafname

⚠️ Caution: the threshold parameter t depends on the genetic proximity of your strains (based on Hamming distances matrix). To adapt this parameter to your bacterial species, just launch the following command once (ORI.py threshold). It gives a figure threshold.png as output containing the distribution of the distances between the strains of the index (distances are multiplied by 1e05 in the figure). More generally, if your cluster of strains is too large and gives you to many possibilities of identification, try a lower t value (e.g. i know that the strains number 205, 51 and 55 are really closed on a phylogenic tree, but a bit farther away to strains 54 and 78; if a threshold of 0.0002 (default value) gives you a cluster containing the 5 strains, you can lower to -t 0.0001 to obtained two separated clusters).

ORI.py threshold -m hamming_matrix.tsv -t 0.0002

The default 0.0002 value is the value used to merge Streptococcus thermophilus strains using filters of size 0.5G. This value must be modified in the case of using another species and/or another filter size.

Once you have defined your own t value, merge your strains in adapted clusters:

howdesbt distance --list=leafname --threshold=0.0002 --matrix=hamming_matrix.bin --merge

ORI.py clean_merge -n leafname -r path/to/repository/with/bf/files -o list_number_file.txt

ls *.bf > leafname_merge

Since the genomes of some strains have been merged, the size of these clusters must also be recalculated:

ORI.py merge_length -b leafname_merge -l length.txt -c list_number_file.txt -o merge_length.txt

To run these commands you must be in the directory containing the .bf files.

howdesbt cluster --list=leafname/or/leafname_merge --tree=union.sbt --nodename=node{number} --cull

howdesbt build --howde --tree=union.sbt --outtree=howde.sbt

Once the compressed bloom filters have been created, we can delete those that are not compressed:

ls | grep -Pv 'detbrief.rrr.' | grep '.bf' | xargs rm --

In order to facilitate identification it may be wise to remove reads of too poor quality. For this it is possible to use:

ORI.py suppr_bad_reads -fq fastq -q min_quality_value -l min_length_value

As we show in the publication (coming soon), ORI's identification is better with 4000 reads than with 16000 due to noise related to sequencing errors. It is therefore advisable to reduce the number of reads with for example:

head -n 16000 fastq_file_better_than_number.fastq > fastq_file_4000_reads.fq

Then we can start the identification:

howdesbt queryQ --sort --qgram=path/to/seedfile.txt --tree=path/to/howde.sbt --threshold=0.5  fastq_file_4000_reads > path/to/results_howde.txt

ORI.py matrix -f path/to/results/from/HowDeSBT -l path/to/leafname/or/leafname_merge -o path/to/results/matrix.tsv

ORI.py identification -m path/to/matrix.tsv -f path/to/results/from/HowDeSBT -le path/to/length.txt/or/merge_length.txt -l path/to/leafname/or/leafname_merge -c path/to/clingo/or/$(which clingo)(with the conda installation)

If the genomes of close strains were not merged during the creation of the index:

num=0; for i in `cat path/to/leafname`;do echo -e "${num} \t ${i}" >> lnf.txt; let num++; done

The results are not very readable (especially in case of merge of close strains), it is possible to have cleaner results:

ORI.py beautiful_results -f path/to/results/from/ORI -n path/to/lnf.txt/or/list_number_file.txt --pie_chart