4. Taxonomic profiling

1. Taxonomic classification using a reference database of protein sequences using Kaiju

makeDB.sh -r

I] In ~/practical copy the data for this exercise from the shared disk:

cp -R /net/software/practical/4 .

II] Change the permissions of this directory:

chmod -R a+rwx 4

III] In ~/practical/4/ create a directory named kaiju and enter this directory:

mkdir kaiju
cd kaiju

IV] Make symbolic links to the preprocessed FASTQ files from the 3. Trimming & filtering (the paired FASTQ from repair). You might have named the after you own choice - which is ok. In the following examples, we call these input files sample1_R1_pair.fastq.gz and sample1_R2_pair.fastq.gz

ln -s ~/practical/3/repair/sample1_R* .

V] View the Kaiju run options (or arguments), type:

kaiju -h

VI] Run Kaiju against the ProGenomes database:

kaiju -t /net/software/databases/proGenomes/nodes.dmp -f /net/software/databases/proGenomes/kaiju_db.fmi -i <(gunzip -c sample1_R1_pair.fastq.gz) -j <(gunzip -c sample1_R2_pair.fastq.gz) -z 4 -x -s 75 -a greedy -e 5 -o sample1_proGenomes_kaiju.out -v

VII] View the top of the Kaiju output:

head sample1_proGenomes_kaiju.out

❓Which NCBI taxon does the first sequence classified (C) read represent?

91803	root	cellular organisms	Bacteria	FCB group	Bacteroidetes/Chlorobi group	Bacteroidetes	Bacteroidia	Bacteroidales	Bacteroidaceae	Bacteroides	Bacteroides ovatus
87964	root	cellular organisms	Bacteria	FCB group	Bacteroidetes/Chlorobi group	Bacteroidetes	Bacteroidia	Bacteroidales	Bacteroidaceae	Bacteroides	Bacteroides thetaiotaomicron

VIII] Convert the Kaiju output file into the proper Krona Tools format using kaiju2krona:

kaiju2krona -t /net/software/databases/proGenomes/nodes.dmp -n /net/software/databases/proGenomes/names.dmp -i sample1_proGenomes_kaiju.out -o sample1_proGenomes_kaiju.krona

IX] The file sample1_proGenomes_kaiju.krona can then be imported into Krona and converted into an HTML file using Krona’s ktImportText program. Activate the Conda Krona environment first!

ktImportText -o sample1_proGenomes_kaiju.html sample1_proGenomes_kaiju.krona

X] Open the HTML file you just generated in a web browser.

❓How many reads did not get assign to bacteria?

❓What seems to be the most abundant genera?

❓How many percent of the reads have mapped to this genera?

❓Locate the genus Escherichia coli. How many of the reads in your sample originate from this specie and how many percent does this species represent from the total reads (Hint: you can type the name in the search field)?

The program kaiju2table can convert Kaiju’s tab-separated output file into a summary report file for a given taxonomic rank, e.g., species.

XI] Print a report at species level with species that comprise at least 1 percent of the total reads. Name the output file containing the report sample1_proGenomes_kaiju_species.summary:

kaiju2table -t /net/software/databases/proGenomes/nodes.dmp -n /net/software/databases/proGenomes/names.dmp -r species -m 0.1 -o sample1_proGenomes_species.summary sample1_proGenomes_kaiju.out

❓ What does the -r and -m flag do?

❓How many percent of the reads belongs to the most abundant species in the sample?

❓ How many of the reads in your sample originate from Lactobacillus rhamnosus?

❓ Look at the last line in the report, how many percent of the reads could not be classified (does not match anything in the database)?

❓ How many percent of the classified reads could not be assign at species level?

XII] Create a report on genus level. Name the output file containing the report sample1_proGenomes_genus.summary

❓How many percent of the classified reads could not be assign at genus level?

❓ What is the most abundant genera, and how many percent of the reads have mapped to this genera?

2. Taxonomic classification using a K-mer based method

I] Create a directory named kraken in ~/practical/4/ and make symbolic links to the preprocessed FASTQ files from the 3. Trimming & filtering - the paired FASTQ from repair (sample1_R1_pair.fastq.gz and sample1_R2_pair.fastq.gz).

II] Activate the correct Conda environment and view the Kraken run options (or arguments), type:

kraken -h

In order to classify reads with Kraken, you need to tell the program where the database is locate. We have downloaded the reduced database, formatted it and placed it here: /net/software/databases/minikraken/.

II] Run Kraken on the two FASTQ files and name the output sample1_kraken.out. Remember that it is a paired end sample:

kraken --db /net/software/databases/minikraken/ --paired sample1_R1_pair.fastq.gz sample1_R2_pair.fastq.gz --output sample1_kraken.out --threads 4

❓ What does the --paired flag do?

III] Look at the first lines in output file (head sample1_kraken.out).

IV] In order to make a Krona chart of the Kraken output you need to remove some of the columns in the sample1_kraken.out file:

cut -f2,3 sample1_kraken.out > sample1_kraken.krona.in

V] Make a Krona chart from the file sample1_kraken.krona.in using Krona’s ktImportText program (activate the Krona environment):

ktImportTaxonomy sample1_kraken.krona.in -o sample1_kraken.html

VI] Open the Krona chart from the Kraken analysis.

❓How many percent of the reads belongs to the most abundant species in the sample?

❓What is the most abundant genera, and how of the reads have mapped to this genera?

❓Does any of the reads map to the specie Escherichia coli? How many of the reads in your sample originate from this specie?

VII] Generate a taxonomic report of the Kraken output. NB! It is important that the database used must be the same as the one used to generate the output file (activate the Kraken environment):

kraken-report --db /net/software/databases/minikraken/ sample1_kraken.out > sample1_kraken.rpt

VIII] Open the Kraken report in a text editor (e.g. Gedit)

❓How many percent of the reads were unclassified?

❓How many percent of the reads belongs to the most abundant species in the sample?

XI] Compare this report with the report from the run against the smaller database (minikraken).

❓Against which database does the analysis produce least unclassified reads and how many reads are unclassified?

❓Are there any major differences in the taxa that are found using the two different databases?

X] Make a Krona chart from the sample1_stor_kraken.out

❓ Can you see any taxonomic differences between the analysis against Minikraken and the complete Kraken database in the Krona plot?

3. Taxonomic classification using clade specific markes

I] Create a directory named metaphlan in ~/practical/4/ and make symbolic links to the preprocessed FASTQ files from the sample1_R1_pair.fastq.gz and sample1_R2_pair.fastq.gz from the exercise 3. Trimming & filtering - the paired FASTQ from repair

II] View the MetaPhlAn3 run options (or arguments), type:

metaphlan2.py -h

III] Run MetaPhlAn3 on the two FASTQ files and name the output sample1_metaphlan.txt. Remember that it is a paired end sample:

metaphlan sample1_R1_pair.fastq.gz,sample1_R2_pair.fastq.gz --bowtie2db /opt/databases/metaphlan3_db --bowtie2out sample1.bz2 --nproc 16 --input_type fastq -o sample1_metaphlan.txt

# --bowtie2db, location of the database

❓ What does the --nproc flag do?

IV] Look at the start of this file:

bzcat sample1.bz2 | head

V] Look at the first lines in output file:

head sample1_metaphlan.txt

After three lines with information on how the output was generated, the abundance profile should look similar to this - a text file with four columns:

The first column lists clades, ranging from taxonomic kingdoms (Bacteria, Archaea, etc.) through species. The taxonomic level of each clade is prefixed to indicate its level:

Kingdom: k__, Phylum: p__, Class: c__, Order: o__, Family: f__, Genus: g__, Species: s__

The second column list the NCBI taxonomy of the taxa

The third column list the relative abundances of the identified taxa. Since sequence-based profiling is relative and does not provide absolute cellular abundance measures, clades are hierarchically summed. Each level will sum to 100%; that is, the sum of all kingdom-level clades is 100%, the sum of all genus-level clades (including unclassified) is also 100%, and so forth.

In order to make a Krona chart of the MetaPhlAn3 output you need to convert it using the script metaphlan2krona.py. However, this script was written for MetaPhlAn2 which generate an abundance profile output with only two columns:

In order for metaphlan2krona.py to work properly, you first need to do some manipulations to the MetaPhlAn3 output to simulate a MetaPhlAn2 output

V] Remove column 2 and 4 from the sample1_metaphlan.txt:

cut -f1,3 sample1_metaphlan.txt > sample1_metaphlan_ok.txt

VI] Convert the output to Krona format:

metaphlan2krona.py -p sample1_metaphlan_ok.txt -k sample1_metaphlan.krona

VII] Then make the HTML plot in the correct conda environment:

ktImportText -o sample1_metaphlan.html sample1_metaphlan.krona

❓ Can you see any major taxonomic differences between this taxonomic profile, and the profiles you got from the analysis using Kaiju and Kraken?

❓How many of the reads in your sample originate from Escherichia coli?

4. Visualisation of taxonomic profiles with MEGAN

I] Create a directory named megan in ~/practical/4/ and make symbolic links to the Kaiju output against RefSeq (sample1_refseq_kaiju.out):

-ln -s ~/practical/4/optional_refseq_kaiju/sample1_refseq_kaiju.out

II] Look at the first part of the file sample1_refseq_kaiju.out:

head sample1_refseq_kaiju.out

III] Add taxonomy labels (names not just numbers) to each assigned read using kaiju-addTaxonNames, but do not specify the taxonomic rank, since some reads have been assigned to species level, while other have been assign to lower taxonomic levels.

kaiju-addTaxonNames -t /net/software/databases/refseq_kaiju/nodes.dmp -n /net/software/databases/refseq_kaiju/names.dmp -i sample1_refseq_kaiju.out -o sample1_refseq_kaiju.out.taxa

OUTPUT: sample1_refseq_kaiju.out.taxa

IV] Look at the first part of the file sample1_refseq_kaiju.out.taxa.

❓ Do you spot any changes in this file relative to sample1_refseq_kaiju.out?

V] Summarise all taxa in
sample1_refseq_kaiju.out.taxa using this command:

awk -F '\t' 'BEGIN{FS=OFS="\t"}{print $8}' sample1_refseq_kaiju.out.taxa | sort | uniq -c | sed -E 's/^ *//; s/ /\t/' | awk 'BEGIN {FS="\t"; OFS="\t"} {print $2, $1}' > sample1_refseq_kaiju_megan.txt

# awk remove all columns except column 8 containing the taxonomic label
# sort will sort the taxa
# uniq will count the number of occurrences of each taxa, collapse identical taxa into one and write a second row with the number of occurrences for each taxa
# sed will change the delimiter to tab
# awk will change order of column 1 (counts) and column 2 (taxa)

OUTPUT: sample1_refseq_kaiju_megan.txt

❓ How many unique taxa are there in the sample?

VI] Open the file sample1_refseq_kaiju_megan.txt in a text editor. The first lines should look something like this:

	513462
Acaryochloris marina MBIC11017	1
Acetoanaerobium sticklandii	6
Acetobacteraceae	1

VII] Add “NA” to all the unassigned reads in the first row, second column in a text editor:

NA	513462
Acaryochloris marina MBIC11017	1
Acetoanaerobium sticklandii	6
Acetobacteraceae	1

VIII] Sort sample1_refseq_kaiju_megan.txt so that the most abundant taxa are at the top:

# -t specify tab as delimiter
# -k which column to sort by
# -nr numerical values with the most abundant at the top

IX] Open sample1_refseq_kaiju_megan_sorted.txt in a text editor and remove all taxa with less than 100 counts (sequence reads). Save the new file as sample1_refseq_kaiju_megan_100.txt

❓How many taxa are left?

X] Open MEGAN6:

MEGAN

XII] Press ‘File’, select ‘Import’ and ‘CSV Format’:

XIII] Select the sample1_refseq_kaiju_megan_100.txt to load the data into MEGAN:

XIV] Select ‘Class Count’, ‘Tab’ and ‘Taxonomy’. Press ‘Apply’:

A new window with the taxonomic assignments on a hierarchical tree structure will appear.

In addition there is some information about the file you just imported in the “Message window”:

XV] Click on the circle representing Bacteria. Two numbers will appear:

1. Assigned means how many reads that were classified to the Bacteria taxa - that is on domain level

2. Summed means how many reads in total in the sample that are assigned to this AND deeper hierarchical levels of this domain

❓How many reads are assigned to the Bacteria level and not placed further down in the taxonomic hierarchy (eg. at family or genus level)?

XV] Open sample1_refseq_kaiju_megan_100.txt.

❓How many reads are assigned to the Bacteria level by Kaiju? Is it the same number as in MEGAN?

XVI] Press the “Rank” button, and expand the tree to the Family level.

❓What does the different sizes of the circles represent?

❓How many reads are in total assigned to the family Bifidobacteriaceae?

XVII] Click on the “Rank” button and select Species. MEGAN will try to expand the taxonomic tree all the way to species level.

The tree will expand to species level, and taxa assigned to Species level will be highlighted. Notice that some taxa are not highlighted, eg. Campylobacter.

❓ Why do you think there are no leaves on species level for the genus Campylobacter?

XVIII] Click on the circle representing Bifidobacterium longum.

❓ How many reads are in total were assigned at subspecies level?

You can double check this by looking into your sample1_refseq_kaiju_megan_100.txt file

XIX] Either way you choose, select all leaves at Phylum level and click on the “Draw data as heatmaps”.

❓ What does the darker green colour mean?

XX] Click on the symbol “Show charts” and select “Show bar charts”. In the new window showing the distribution of reads at the level of the sample that appears, click on the “Classes” tab in the left part of the window (next to “Series”).

XXI] Try to change the appearance of the taxonomic profile by selecting some of the other chart views on the top, for example “Bar chart”

❓ Approximately how many percent of the reads belong to Actinobacteria? Hint: Press the % sign at the top.

XXII] Make a word cloud of the taxonomic profile before you close this MEGAN.