KRAKEN2 with complete database

[rx32940]

• in the previous analysis, only minikraken database was used for computational purpose
• however, for the accuracy of the evaluation, both standard database for kraken2 and a custom database built with bacteria, archaea, viral and rat reference sequence is going to be built for the best kraken2 performance.
  • rat reference genome from NCBI refseq : https://www.ncbi.nlm.nih.gov/genome/?term=rat

[rx32940]

Kraken2 standard library: default kmer length 35, minimizer length 31. NCBI taxonomy, bacterial, archaeal, and viral domains, human and a collection of known vectors (UniVec_Core) Kraken2 Custom library: bacteria, archaea, viral and rat reference.

• in the previous analysis, library was used. code to build the kraken2 both standard and custom database: https://github.com/rx32940/Lepto-Metagenomics/blob/f913e682aff7f038670b821a12fd4154109aae6c/run-kraken2.sh

[rx32940]

run kraken2/bracken in phylum and genus level: dir for kraken2 output with standard db:

/scratch/rx32940/kraken/output/kraken_out

dir for kraken2 output with standard db after bracken estimation:

/scratch/rx32940/kraken/output/bracken_out

dir for kraken2 output with custom db:

/scratch/rx32940/kraken/output/custom_out

dir for kraken2 output with custom db after bracken estimation:

/scratch/rx32940/kraken/output/custom_bracken

code used to run kraken2 and bracken in phylum and genus levels

[rx32940]

• by using standard library, large portion of reads was classified as homo, for example, 96.571% in sample R22.K's were classified into homo on the genus level. (53.15% reads were unclassified)
• in this library, rat, which is the host genome was not included.
• it is not reasonable to have reads classified under human, because these are rats samples
• reads classified into human can be due to similarities between rats and human genome (eukaryotes)
• we will see the classification with custom database before we reach a conclusion

output in Dropbox:

$HOME/Dropbox/5. Rachel's projects/Metagenomic_Analysis/Kraken2-standard/standard/phylum(genus)/results

output on sapelo2:

/scratch/rx32940/kraken/output/bracken_out/genus(phylum)

	Classfied (%)	Unclassfied (%)
R22.K	46.85	53.15
R22.L	24.93	75.07
R22.S	46.44	53.56
R26.K	52.83	47.17
R26.L	33.54	66.46
R26.S	47.44	52.56
R27.K	56.64	43.36
R27.L	35.5	74.5
R27.S	42.27	57.73
R28.K	25.73	74.27
R28.L	28.1	71.9
R28.S	26.79	73.21

[rx32940]

with the custom database, rattus replaced homo became the taxa with the highest reads abundance. the percentage of total reads been classified has also increased significantly.

• from this example, we can tell the deficiency of taxonomic profiling. reads can be easily classified under the taxa that have genome similar to the actual taxa the read belongs to.
• these will cause some of the identified taxa in the metagenome samples not making biological sense.
• false positives

output in Dropbox:

$HOME/Dropbox/5. Rachel's projects/Metagenomic_Analysis/Kraken2-standard/custom/phylum(genus)/results

output on sapelo2:

/scratch/rx32940/kraken/output/custom_bracken/genus(phylum)

	Classfied (%)	Unclassfied (%)
R22.K	70.92	29.08
R22.L	30.43	69.57
R22.S	62.96	37.04
R26.K	70.03	29.97
R26.L	44.85	55.15
R26.S	63.24	36.76
R27.K	69.66	30.34
R27.L	32.45	67.55
R27.S	61.73	38.27
R28.K	86.29	13.71
R28.L	83.42	16.58
R28.S	83.75	16.25

[rx32940]

compare the improvement with using the minikraken2 library:

Sample ID	Classified	Unclassified
R22.K	14.72%	85.28%
R22.L	6.03%	93.97%
R22.S	13.46%	86.54%
R26.K	14.45%	85.55%
R26.L	7.55%	92.45%
R26.S	10.83%	89.17%
R27.K	13.85%	86.15%
R27.L	6.62%	93.38%
R27.S	10.89%	89.11%
R28.K	8.58%	91.42%
R28.L	7.45%	92.55%
R28.S	6.52%	93.48%

[rx32940]

Relative Abundance of each genus identified in the metagenomic samples with custom Kraken2 library

• other: summed fraction of genus other than top 10 most abundant genus found in each sample
• top 10 genus for each sample shown in the bar plot with individual color code
• Rattus fraction identified included

Absolute abundance of each genus identified in the metagenomic samples with custom Kraken2 library

[rx32940]

• we can conclude from relative abundance plot that Lung samples has highest composition of microbes.

• from the absolute abundance plot we are able to conclude that sample with better quality sequences tend to have more reads identified by Kraken2

• [x] figure out why Shigella is missing in relative abundance _solved: ordered top 10 by read counts instead of proportion, because fraction rounds up_

• [x] change relative to absolute in title solved

• [x] check the quality of R28 sequences. why is R28 Lung and spleen have such low microbe composition. while the composition seem normal for kidney, however, mainly responsible by the composition of Leptospira. _solved:The effective reads for R28 samples are three times higher than that of R22, which is the second best with effect reads_

[rx32940]

Relative and absolute abundance

[rx32940]

code for barplots in the last two issues

• to plot genus, replace all "phylum" with "genus", all "Chordata" with "Rattus"
• run relative and absolute separately

[rx32940]

• The classification system for viruses for both KRAKEN2 and CLARK is very tricky.

• Due to the differences in taxonomy level between viruses and bacteria, viruses can't be expected to classified under the target phylum, genus or species like bacteria does.

• however, both software requires the classification levels indicated for analysis (Kraken2 didn't, but required for bracken for abundance estimation and more accurate results)

• Both Softwares chose to have viruses classified with a Genus and under it, species

• However, while Clark under the genus level classification classified a large percentage of reads to Pandoravirus, none of the top 10 abundant taxa classified by Kraken2 under genus or species level were viruses.

• a large number of reads were classified as Pandoravirus with KRAKEN2, but only after many more bacteria taxa.

  • [X] maybe try to do top 20 most abundant taxa (code: First part of the code before plotting)
  • Pandoravirus found, more overlaps between taxa identified by two software, but most abundant taxa are still very diverged.

• [ ] maybe try to find top 20 taxa with clark results too, then do comprehensive comparison between identified taxa

/Users/rx32940/Dropbox/5.Rachel-projects/Metagenomic_Analysis/Kraken2-standard/custom/genus/top_20_genus_for_each_sample.csv

• Roseolovirus, even though was found in Clark analysis as one of the most abundant taxa, was not classified by kraken2 (do exist in Kraken2 database).