4. Manual: Taxonomic classification

4.1 Diamond

Now we will annotate the aggregated contigs by assigning taxonomic classifications to them based on sequence similarity to known proteins in a database using diamond blastx.

Please take note of the following blastx parameters: -f (--outfmt), -b (--block-size) and -c (--index-chunks), documentation can be found here.

The -f 6 parameter will ensure the output is in a tabular format. The 6 may be followed by a space-separated list of various keywords, each specifying a field of the output.

The -b parameter is the main parameter for controlling the program’s memory and disk space usage. Bigger numbers will increase the use of memory and temporary disk space, but also improve performance. The program can be expected to use roughly six times this number of memory (in GB). The default value is -b 2. The parameter can be decreased for reducing memory use, as well as increased for better performance (values of >20 are not recommended).

The -c parameter controls the number of chunks for processing the seed index. This option can be additionally used to tune the performance. The default value is -c 4, while setting this parameter to -c 1 instead will improve the performance at the cost of increased memory use.

diamond blastx \
-q {input} \
-d {db} \
-o {output} \
-f 6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore staxids \
--threads {threads} \
-b 10 -c 1

{input} is the contig file created in either step 3.3 or step 3.1, depending on your amount of samples.
{db} is the protein database to be searched against.
{output} is a .tsv file containing the annotation results.

Multiple different databases can be found in:

/mnt/viro0002-data/workgroups_projects/Bioinformatics/DB

4.2 Split annotation files

We will split the combined annotation file back into individual annotation files for each sample. This step can be seen as optional if you are dealing with a single sample.

Modify and run:

mkdir -p tmp_split

sed 's/_NODE/|NODE/' {input} | awk -F'|' '{
    identifier = $1;  # Construct the identifier using the first two fields

    output_file = "tmp_split/" identifier;  # Construct the output filename

    if (!seen[identifier]++) {
        close(output_file);  # Close the previous file (if any)
        output = output_file;  # Update the current output file
    }

    print $2 > output;  # Append the line to the appropriate output file
}'

for file in tmp_split/*; do
    mkdir -p {output}/$(basename "$file")/;
    mv "$file" {output}/$(basename "$file")/diamond_output.tsv;
done

rmdir tmp_split

{input} is the combined annotation file from step 4.1.
{output} is a directory path. This directory will be automatically filled with subdirectories for each sample. In each subdirectory you will find a diamond_out.tsv file.

4.3 Parsing diamond output

Now we will process the DIAMOND output files with a custom Python script called post_process_diamond_v1.0.py. This script will further enrich taxonomic information for each contig based on the DIAMOND alignment results. If a contig has multiple matches in the database, it will select the best hit based on a combined score of bitscore and length. Lastly, it separates the contigs into two lists: those that were successfully annotated and unannotated.

This python script utilizes the biopython library.

python scripts/post_process_diamond_v1.0.py \
-i {input.annotation} \
-c {input.contigs} \
-o {output.annotated} \
-u {output.unannotated}

{input.annotation} is the annotation file step 4.2.
{input.contigs} are the contigs from the SPAdes step 3.1.
{output.annotated} is a .tsv file with a set of annotated contigs.
{output.unannotated} is a .tsv file with a set of unannotated contig IDs.

4.4 Optional: read classification with Kaiju

Classifying contigs is great because it provides higher accuracy annotations and ignores low noise from random sequencing errors. However, you might still want to taxonomically classify your reads to potentially help detect rare species that are unable to be assembled due to low read depth.

Kaiju is a fast taxonomic classification program which translates reads into amino acid sequences and matches them against reference protein databases.

Kaiju’s paired-end mode assumes perfectly synchronized forward and reverse files, it cannot accept a third singleton file in the same command. The workaround for this is to simply run Kaiju twice: once for the paired-end reads, and once for the singletons. Because Kaiju outputs results line-by-line, you can simply cat the two output files together.

# Run Kaiju for R1 and R2 reads
kaiju -t {kaiju_db}/nodes.dmp \
      -f {kaiju_db}/kaiju_db_*.fmi \
      -i {input.R1} \
      -j {input.R2} \
      -o {kaiju_out}.pe \
      -a greedy \
      -z {threads}

# Run kaiju for singleton reads
kaiju -t {kaiju_db}/nodes.dmp \
      -f {kaiju_db}/kaiju_db_*.fmi \
      -i {input.S} \
      -o {kaiju_out}.se \
      -a greedy \
      -z {threads}

# Merge paired and singleton results
cat {output.kaiju_out}.pe {output.kaiju_out}.se > {output.kaiju_out}

{kaiju_db} is the directory of your Kaiju database. This directory should contain a kaiju_db_*.fmi, nodes.dmp, and names.dmp file.
{input.R1}, {input.R2} and {input.S} are the host-filtered FASTQ files (R1, R2, and S) from step 2.3.
{kaiju_out}.pe and {kaiju_out}.se are the output tables.
{kaiju_db}/kaiju_db_*.fmi: Replace this with your actual .fmi file depending on your selected database.

The base Kaiju output is a simplistic table containing only 3 columns:

Either C or U, indicating whether the read is classified or unclassified.
Name of the read
NCBI taxon identifier of the assigned taxon

To make it more informative, we can create Krona files with the kaiju2krona command.

kaiju2krona -t {kaiju_db}/nodes.dmp \
            -n {kaiju_db}/names.dmp \
            -i {output.kaiju_out} \
            -o {output.kaiju_krona}

{output.kaiju_out} is the merged output from the base Kaiju runs.
{output.kaiju_krona} is a tab separated krona file that can be used for standard Krona visualization.

Lastly, we can use the kaiju-addTaxonNames command to add taxonomic names of various levels to our Kaiju results, and then bin reads based on family, genus and species level. This way, you can easily access reads for your virus of interest. We can use the scripts/bin_reads.py script to bin the reads on different taxonomic levels.

# Add superkingdom, phylum, class, order, family, genus, species
kaiju-addTaxonNames -t {kaiju_db}/nodes.dmp \
                    -n {kaiju_db}/names.dmp \
                    -i {output.kaiju_out} \
                    -o {output.kaiju_names} \
                    -r superkingdom,phylum,class,order,family,genus,species
                    
# Bin reads on family, genus and species level
python scripts/bin_reads.py \
  {output.kaiju_names} \
  {input.R1} \
  {input.R2} \
  {output_bin.dir} \
  {input.S}

{output.kaiju_out} is the merged output from the base Kaiju runs.
{output.kaiju_names} is the Kaiju output with extra columns for different taxonomic levels.
{output_bin}.dir is the directory where family, genus, species subdirs are created with their corresponding FASTA files.

Note

We can now move on to the final steps where we will create various files needed for downstream analysis.