4.5 - Classification of sequences with k-mer composition profiles¶
time
- Teaching: 15 minutes
- Teaching: 15 minutes
Objectives and Key points
Objectives¶
- Understand when we may need to use something other than
BLASTfor getting quick classification results. - Know how to run
kraken2using a standard database, and interpret the output files.
Keypoints¶
- k-mer based classification tools provide a rapid means of assigning taxonomy to sequence data.
- The quality of the classification is tool and database dependent, results are also harder to interpret than a BLAST output
- For example, there are no identity/coverage statistics from which judgement calls can be made.
Why not BLAST?¶
As we saw in the level 1 training, running BLASTn to classify sequences is a good way to get annotation information about a sequence, and from this we can infer taxonomic origin of the sequence. However, obtaining classification information is not the primary purpose of BLASTn so we need to perform some manual inspection of the output to make a call on the classification of each sequence.
When we do not care about functional information, and only want classification of sequences, there are a suite of tools designed specifically with this purpose in mind. Rather than perform pairwise sequence alignment between query and target sequences, these tools process their reference database into a compressed form containing sequence fingerprints.
Query sequences are fingerprinted and the results compared to the database providing the most likely, lowest common ancestor for each query sequence.
Options for k-mer classification
How this fingerprinting is done is tool specific, and goes outside the realm of biology. There are a number of tools developed for this kind of work, including;
We will be using kraken2 today, as it is a tool which is used often within our operations.
Getting started with kraken2¶
Navigate to level2/annotation_kraken2/ and run kraken2 with the help paramter (-h).
code
Output
Usage: kraken2 [options] <filename(s)>
Options:
--db NAME Name for Kraken 2 DB
(default: "/opt/nesi/db/Kraken2/standard-2018-09")
--threads NUM Number of threads (default: 1)
--quick Quick operation (use first hit or hits)
--unclassified-out FILENAME
Print unclassified sequences to filename
--classified-out FILENAME
Print classified sequences to filename
--output FILENAME Print output to filename (default: stdout); "-" will
suppress normal output
--confidence FLOAT Confidence score threshold (default: 0.0); must be
in [0, 1].
--minimum-base-quality NUM
Minimum base quality used in classification (def: 0,
only effective with FASTQ input).
--report FILENAME Print a report with aggregrate counts/clade to file
--use-mpa-style With --report, format report output like Kraken 1's
kraken-mpa-report
--report-zero-counts With --report, report counts for ALL taxa, even if
counts are zero
--report-minimizer-data With --report, report minimizer and distinct minimizer
count information in addition to normal Kraken report
--memory-mapping Avoids loading database into RAM
--paired The filenames provided have paired-end reads
--use-names Print scientific names instead of just taxids
--gzip-compressed Input files are compressed with gzip
--bzip2-compressed Input files are compressed with bzip2
--minimum-hit-groups NUM
Minimum number of hit groups (overlapping k-mers
sharing the same minimizer) needed to make a call
(default: 2)
--help Print this message
--version Print version information
If none of the *-compressed flags are specified, and the filename provided
is a regular file, automatic format detection is attempted.
One of the nice perks of working with kraken2 is that the maintainers of the tool provide pre-computed, up-to-date databases for classification. These are released in many different variants - some for use of HPCs and some for local machines, and databases which cover a range of target organisms.
We use PlusPFP database, which contains the following lineages:
- Archaea
- Bacteria
- Fungi
- Human
- Plant
- Protozoa
- Virus
- Plasmid sequences
- Sequencing and transformation vectors (UniVec)
For training, we are using a reduced version of this database (PlusPFP-16) so that our jobs do not need as many resources and do not take too long to run.
This makes it a comprehensive database for most purposes, but if you are tring to annotate insect sequences this is not the right database for you.
Make sure you have the right database for your work
It is critical to know what is in your classification database, as a failure to include the correct reference material in the database will result in sequences being either misclassified, or not classified at all.
As part of the accredited HTS workflow used by the Virology team we keep an up-to-date copy of this database on NeSI, but this will not be appropriate for all teams to use.
Performing basic classification with kraken2¶
To begin, we are going to run kraken2 with a fairly minimal set in input parameters just to see the tool in action.
Complete the following script, then submit your slurm job.
code
The kraken2 output format¶
Once the job is complete, use less or head to take a look at the output file. The output file is a tab-delimited table with the following columns:
input_seqs.out
| Column | Content |
|---|---|
| 1 | A single letter (C or U) designating whether the sequence was classified or not |
| 2 | Sequence name |
| 3 | The NCBI taxonomy ID applied to the sequence (0 if unclassified) Because we used the --use-names flag, this takes the form TAXONOMY NAME (TAXONOMY ID)By default this is just the numeric value |
| 4 | Sequence length For paired-end data, this will be in the form FORWARD|REVERSE |
| 5 | A list of how many k-mers mapped to each taxonomy ID Results are reported as colon-delimited pairings of TAXID:NUMBER |
It is also possible to generate a report file that summarises the number of sequences classified to each taxonomy/rank, using the modified command:
code
Depending on your situation, this may be more useful.
Classification output versus report
When would you need the per-sequence output
If you are screening a sequence set for the specific reads obtained from a particular organism, or are trying to remove contamination from a sample, having the per-sequence results are invaluable.
The content of this file can be easily parsed with the cut and grep tools to extract lists of sequences match or failing to match different search criteria.
When would you need the high-level report
If you trust that your data is uncontaminated, you might be trying ot achive one of the following:
- Detecting the presence of a particular organism in a mixed population (i.e. eDNA surveillance).
- Confirming the identity of a pure culture.
- In this case you might be comfortable taking a majority-rules approach to your data - if most of the sequences are classified as a single species, this is likely the identity of your culture.
Refining the kraken2 classification¶
By default kraken2 does not apply any filtering to the classification data. Reads are classified based on the most commonly reported taxonomy for a sequence, with no regard for how good the evidence for this classification can be. There is a built-in parameter to increase the stringency of classification. The filter is applied during classification and cannot be run after the fact.
This means that if you want to change your filtering criteria you must perform classification again. This differs from tools like BLAST, where a more stringent filter can be applied after running classification.
For example, if we run a
BLASTsearch looking to report all hits with greater than 30% identity, we can later filter the table to keep results with greater than 50% identity by filtering the appropriate column.
Limitations with the kraken2 filtering
It is also difficult to interpret the numeric values of filters.
According to the kraken2 manual, the confidence scoring does not have a simple interpretation.
For example, with a BLAST identitiy of 90%, we can infer that 90% of positions are identical between query and target match, or a bootstrap support of 90% would mean that 90% of the resamplings agree with the full result.
However, there is no such interpretation of what a 0.9 confidence value in kraken2 means, other than it is more strict that a value of 0.5.
We generally don't need to apply much of a cutoff to kraken2 - the tool is already good at identifying sequences which can't be reliably classified.
code
Comparing the outputs¶
Without going into too much detail, we can briefly compare the contents of both output files. For starters, we can use some of the basic bash command line tools to work out how many sequences were in our input and output files.
code
The number of results in each output file is the same, and also the same as the number of sequences in the input file. This is different to BLAST, where an unclassified or unannotated sequence is not reported in the output.
A quick way to look for differences in the classification rate of the unfiltered and filtered kraken2 output is to apply the following commands:
If we point this towards the correct column, we can count how many sequences were successfully classified or not.
Exercise
Look through these tutorial notes, or the kraken2 documentation, to find the appropriate column index for the command above. Once you have determined the correct value for ${n}, complete the command above and run it for both kraken2 outputs.
How many sequences were classified for each run?
As you can see, even a small degree of filtering has a drastic impact on the classification rate. If you are applying a filter value, using a low threshold is recommended to avoid discarding too much usable data.
How reliable are these results?¶
The input sequences in this exercise are simulated Oxford Nanopore reads from the following organisms;
- Tomato (Solanum lycopersicum) (5,447 sequences)
- Xyllela fastidiosa (166 sequences)
- Potato virus X (1 sequence)
If you skim through the results of both output files, you should note the following:
outputs/input_seqs.out
| Lineage | Observations |
|---|---|
| Solanum lycopersicum | 4,827 |
| Homo sapiens | 36 |
| Xylella fastidiosa | 129 |
| Potato virus X | 1 |
outputs/input_seqs.conf_0.1.out
| Lineage | Observations |
|---|---|
| Solanum lycopersicum | 991 |
| Homo sapiens | 14 |
| Xylella fastidiosa | 78 |
| Potato virus X | 1 |
What should be apparent here is that while filtering did remove a lot of the noise from our sample (we lost species we know are not truly present) we also lost a signal that was in the sample.
Recommendation for kraken2
Don't let this put you off using kraken2 - it is an excellent tool for rapidly sifting through HTS data for quick summaries. However, results should only be considered indicative in nature, not definitive.
Use kraken2 to make a quick check for names of interest (regulated organisms, host, contamination) but follow up with more thorough investigation.