kraken2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6a99628a-adb7-444a-bf7a-ebd5808dd814",
   "metadata": {},
   "source": [
    "# Taxonomic Profiling using Kraken2\n",
    "\n",
    "There are a number of pre-computed databases here: https://benlangmead.github.io/aws-indexes/k2 . For demonstration purposes, let us use one of the smallest ones : Standard-8. This has already been downloaded in data/kraken2_db/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "283c7186-c3dc-430c-8103-7ae6786d8ac3",
   "metadata": {},
   "outputs": [],
   "source": [
    "! mkdir -p results/kraken2_output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31d6ca59-6cb8-4c8e-b4c9-4514e8e4d527",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "To run on one sample, you would do :\n",
    "\n",
    "`kraken2 --paired --db data/kraken2_db/standard-8/ data/metagenome_samples/FH1_1.fastq.gz data/metagenome_samples/FH1_2.fastq.gz --output results/kraken2_output/FH1`\n",
    "\n",
    "To run on all at once, use the script below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "13345d1f-66d0-4a8e-b50d-7a81f81d79c9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing sample: BH1\n",
      "Loading database information... done.\n",
      "1000000 sequences (302.00 Mbp) processed in 26.460s (2267.6 Kseq/m, 684.81 Mbp/m).\n",
      "  705810 sequences classified (70.58%)\n",
      "  294190 sequences unclassified (29.42%)\n",
      "Processing sample: BH2\n",
      "Loading database information... done.\n",
      "1000000 sequences (302.00 Mbp) processed in 33.129s (1811.1 Kseq/m, 546.94 Mbp/m).\n",
      "  336398 sequences classified (33.64%)\n",
      "  663602 sequences unclassified (66.36%)\n",
      "Processing sample: BH3\n",
      "Loading database information... done.\n",
      "1000000 sequences (302.00 Mbp) processed in 43.397s (1382.6 Kseq/m, 417.54 Mbp/m).\n",
      "  208838 sequences classified (20.88%)\n",
      "  791162 sequences unclassified (79.12%)\n",
      "Processing sample: BH4\n",
      "Loading database information... done.\n",
      "1000000 sequences (302.00 Mbp) processed in 41.553s (1443.9 Kseq/m, 436.07 Mbp/m).\n",
      "  261480 sequences classified (26.15%)\n",
      "  738520 sequences unclassified (73.85%)\n",
      "Processing sample: FH1\n",
      "Loading database information... done.\n",
      "1000000 sequences (302.00 Mbp) processed in 47.106s (1273.7 Kseq/m, 384.67 Mbp/m).\n",
      "  570952 sequences classified (57.10%)\n",
      "  429048 sequences unclassified (42.90%)\n",
      "Processing sample: FH2\n",
      "Loading database information... done.\n",
      "1000000 sequences (302.00 Mbp) processed in 41.208s (1456.0 Kseq/m, 439.72 Mbp/m).\n",
      "  452272 sequences classified (45.23%)\n",
      "  547728 sequences unclassified (54.77%)\n",
      "Processing sample: FH3\n",
      "Loading database information... done.\n",
      "1000000 sequences (302.00 Mbp) processed in 42.042s (1427.1 Kseq/m, 430.99 Mbp/m).\n",
      "  517900 sequences classified (51.79%)\n",
      "  482100 sequences unclassified (48.21%)\n",
      "Processing sample: FH4\n",
      "Loading database information... done.\n",
      "1000000 sequences (302.00 Mbp) processed in 34.320s (1748.2 Kseq/m, 527.97 Mbp/m).\n",
      "  731162 sequences classified (73.12%)\n",
      "  268838 sequences unclassified (26.88%)\n"
     ]
    }
   ],
   "source": [
    "!for sample in $(ls data/metagenome_samples/*_1.fastq.gz | sed 's/.*\\///' | sed 's/_1.*//'); do \\\n",
    "    echo \"Processing sample: ${sample}\"; \\\n",
    "    kraken2 --paired --db data/kraken2_db/standard-8/  --output results/kraken2_output/${sample}.kraken.output --report results/kraken2_output/${sample}.kraken.report \\\n",
    "    data/metagenome_samples/${sample}_1.fastq.gz data/metagenome_samples/${sample}_2.fastq.gz;\\\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "01c944fb-a4c6-476a-8891-50cbe22f3eb8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 29.42\t294190\t294190\tU\t0\tunclassified\n",
      " 70.58\t705810\t485\tR\t1\troot\n",
      " 70.50\t704961\t2859\tR1\t131567\t  cellular organisms\n",
      " 68.09\t680862\t29099\tD\t2\t    Bacteria\n",
      " 63.41\t634069\t38696\tP\t1224\t      Pseudomonadota\n",
      " 55.35\t553477\t2461\tC\t28216\t        Betaproteobacteria\n",
      " 54.62\t546158\t9330\tO\t80840\t          Burkholderiales\n",
      " 51.18\t511791\t12223\tF\t506\t            Alcaligenaceae\n",
      " 49.69\t496927\t34939\tG\t222\t              Achromobacter\n",
      " 44.69\t446934\t1893\tG1\t2626865\t                unclassified Achromobacter\n"
     ]
    }
   ],
   "source": [
    "!head  results/kraken2_output/BH1.report \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21efe1b1-19d0-405a-8b0e-dc974b51e074",
   "metadata": {},
   "source": [
    "bracken -d data/kraken2_db/standard-8/ -l G -i results/kraken2_output/BH1.report -o results/kraken2_output/BH1.bracken.output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8f9115c6-85df-4ba9-b421-196a986c8902",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " >> Checking for Valid Options...\n",
      " >> Running Bracken \n",
      "      >> python src/est_abundance.py -i results/kraken2_output/BH1.kraken.report -o results/kraken2_output/BH1.bracken.output -k data/kraken2_db/standard-8/database100mers.kmer_distrib -l G -t 0\n",
      "PROGRAM START TIME: 01-17-2025 02:26:48\n",
      ">> Checking report file: results/kraken2_output/BH1.kraken.report\n",
      "BRACKEN SUMMARY (Kraken report: results/kraken2_output/BH1.kraken.report)\n",
      "    >>> Threshold: 0 \n",
      "    >>> Number of genuses in sample: 1118 \n",
      "\t  >> Number of genuses with reads > threshold: 1118 \n",
      "\t  >> Number of genuses with reads < threshold: 0 \n",
      "    >>> Total reads in sample: 1000000\n",
      "\t  >> Total reads kept at genuses level (reads > threshold): 598537\n",
      "\t  >> Total reads discarded (genuses reads < threshold): 0\n",
      "\t  >> Reads distributed: 107257\n",
      "\t  >> Reads not distributed (eg. no genuses above threshold): 16\n",
      "\t  >> Unclassified reads: 294190\n",
      "BRACKEN OUTPUT PRODUCED: results/kraken2_output/BH1.bracken.output\n",
      "PROGRAM END TIME: 01-17-2025 02:26:49\n",
      "  Bracken complete.\n",
      " >> Checking for Valid Options...\n",
      " >> Running Bracken \n",
      "      >> python src/est_abundance.py -i results/kraken2_output/BH2.kraken.report -o results/kraken2_output/BH2.bracken.output -k data/kraken2_db/standard-8/database100mers.kmer_distrib -l G -t 0\n",
      "PROGRAM START TIME: 01-17-2025 02:26:49\n",
      ">> Checking report file: results/kraken2_output/BH2.kraken.report\n",
      "BRACKEN SUMMARY (Kraken report: results/kraken2_output/BH2.kraken.report)\n",
      "    >>> Threshold: 0 \n",
      "    >>> Number of genuses in sample: 1456 \n",
      "\t  >> Number of genuses with reads > threshold: 1456 \n",
      "\t  >> Number of genuses with reads < threshold: 0 \n",
      "    >>> Total reads in sample: 1000000\n",
      "\t  >> Total reads kept at genuses level (reads > threshold): 218922\n",
      "\t  >> Total reads discarded (genuses reads < threshold): 0\n",
      "\t  >> Reads distributed: 117473\n",
      "\t  >> Reads not distributed (eg. no genuses above threshold): 3\n",
      "\t  >> Unclassified reads: 663602\n",
      "BRACKEN OUTPUT PRODUCED: results/kraken2_output/BH2.bracken.output\n",
      "PROGRAM END TIME: 01-17-2025 02:26:49\n",
      "  Bracken complete.\n",
      " >> Checking for Valid Options...\n",
      " >> Running Bracken \n",
      "      >> python src/est_abundance.py -i results/kraken2_output/BH3.kraken.report -o results/kraken2_output/BH3.bracken.output -k data/kraken2_db/standard-8/database100mers.kmer_distrib -l G -t 0\n",
      "PROGRAM START TIME: 01-17-2025 02:26:49\n",
      ">> Checking report file: results/kraken2_output/BH3.kraken.report\n",
      "BRACKEN SUMMARY (Kraken report: results/kraken2_output/BH3.kraken.report)\n",
      "    >>> Threshold: 0 \n",
      "    >>> Number of genuses in sample: 1624 \n",
      "\t  >> Number of genuses with reads > threshold: 1624 \n",
      "\t  >> Number of genuses with reads < threshold: 0 \n",
      "    >>> Total reads in sample: 1000000\n",
      "\t  >> Total reads kept at genuses level (reads > threshold): 102386\n",
      "\t  >> Total reads discarded (genuses reads < threshold): 0\n",
      "\t  >> Reads distributed: 106329\n",
      "\t  >> Reads not distributed (eg. no genuses above threshold): 123\n",
      "\t  >> Unclassified reads: 791162\n",
      "BRACKEN OUTPUT PRODUCED: results/kraken2_output/BH3.bracken.output\n",
      "PROGRAM END TIME: 01-17-2025 02:26:49\n",
      "  Bracken complete.\n",
      " >> Checking for Valid Options...\n",
      " >> Running Bracken \n",
      "      >> python src/est_abundance.py -i results/kraken2_output/BH4.kraken.report -o results/kraken2_output/BH4.bracken.output -k data/kraken2_db/standard-8/database100mers.kmer_distrib -l G -t 0\n",
      "PROGRAM START TIME: 01-17-2025 02:26:49\n",
      ">> Checking report file: results/kraken2_output/BH4.kraken.report\n",
      "BRACKEN SUMMARY (Kraken report: results/kraken2_output/BH4.kraken.report)\n",
      "    >>> Threshold: 0 \n",
      "    >>> Number of genuses in sample: 1656 \n",
      "\t  >> Number of genuses with reads > threshold: 1656 \n",
      "\t  >> Number of genuses with reads < threshold: 0 \n",
      "    >>> Total reads in sample: 1000000\n",
      "\t  >> Total reads kept at genuses level (reads > threshold): 129064\n",
      "\t  >> Total reads discarded (genuses reads < threshold): 0\n",
      "\t  >> Reads distributed: 132413\n",
      "\t  >> Reads not distributed (eg. no genuses above threshold): 3\n",
      "\t  >> Unclassified reads: 738520\n",
      "BRACKEN OUTPUT PRODUCED: results/kraken2_output/BH4.bracken.output\n",
      "PROGRAM END TIME: 01-17-2025 02:26:50\n",
      "  Bracken complete.\n",
      " >> Checking for Valid Options...\n",
      " >> Running Bracken \n",
      "      >> python src/est_abundance.py -i results/kraken2_output/FH1.kraken.report -o results/kraken2_output/FH1.bracken.output -k data/kraken2_db/standard-8/database100mers.kmer_distrib -l G -t 0\n",
      "PROGRAM START TIME: 01-17-2025 02:26:50\n",
      ">> Checking report file: results/kraken2_output/FH1.kraken.report\n",
      "BRACKEN SUMMARY (Kraken report: results/kraken2_output/FH1.kraken.report)\n",
      "    >>> Threshold: 0 \n",
      "    >>> Number of genuses in sample: 1643 \n",
      "\t  >> Number of genuses with reads > threshold: 1643 \n",
      "\t  >> Number of genuses with reads < threshold: 0 \n",
      "    >>> Total reads in sample: 1000000\n",
      "\t  >> Total reads kept at genuses level (reads > threshold): 469557\n",
      "\t  >> Total reads discarded (genuses reads < threshold): 0\n",
      "\t  >> Reads distributed: 101394\n",
      "\t  >> Reads not distributed (eg. no genuses above threshold): 1\n",
      "\t  >> Unclassified reads: 429048\n",
      "BRACKEN OUTPUT PRODUCED: results/kraken2_output/FH1.bracken.output\n",
      "PROGRAM END TIME: 01-17-2025 02:26:50\n",
      "  Bracken complete.\n",
      " >> Checking for Valid Options...\n",
      " >> Running Bracken \n",
      "      >> python src/est_abundance.py -i results/kraken2_output/FH2.kraken.report -o results/kraken2_output/FH2.bracken.output -k data/kraken2_db/standard-8/database100mers.kmer_distrib -l G -t 0\n",
      "PROGRAM START TIME: 01-17-2025 02:26:50\n",
      ">> Checking report file: results/kraken2_output/FH2.kraken.report\n",
      "BRACKEN SUMMARY (Kraken report: results/kraken2_output/FH2.kraken.report)\n",
      "    >>> Threshold: 0 \n",
      "    >>> Number of genuses in sample: 1649 \n",
      "\t  >> Number of genuses with reads > threshold: 1649 \n",
      "\t  >> Number of genuses with reads < threshold: 0 \n",
      "    >>> Total reads in sample: 1000000\n",
      "\t  >> Total reads kept at genuses level (reads > threshold): 382610\n",
      "\t  >> Total reads discarded (genuses reads < threshold): 0\n",
      "\t  >> Reads distributed: 69660\n",
      "\t  >> Reads not distributed (eg. no genuses above threshold): 2\n",
      "\t  >> Unclassified reads: 547728\n",
      "BRACKEN OUTPUT PRODUCED: results/kraken2_output/FH2.bracken.output\n",
      "PROGRAM END TIME: 01-17-2025 02:26:50\n",
      "  Bracken complete.\n",
      " >> Checking for Valid Options...\n",
      " >> Running Bracken \n",
      "      >> python src/est_abundance.py -i results/kraken2_output/FH3.kraken.report -o results/kraken2_output/FH3.bracken.output -k data/kraken2_db/standard-8/database100mers.kmer_distrib -l G -t 0\n",
      "PROGRAM START TIME: 01-17-2025 02:26:51\n",
      ">> Checking report file: results/kraken2_output/FH3.kraken.report\n",
      "BRACKEN SUMMARY (Kraken report: results/kraken2_output/FH3.kraken.report)\n",
      "    >>> Threshold: 0 \n",
      "    >>> Number of genuses in sample: 1633 \n",
      "\t  >> Number of genuses with reads > threshold: 1633 \n",
      "\t  >> Number of genuses with reads < threshold: 0 \n",
      "    >>> Total reads in sample: 1000000\n",
      "\t  >> Total reads kept at genuses level (reads > threshold): 431225\n",
      "\t  >> Total reads discarded (genuses reads < threshold): 0\n",
      "\t  >> Reads distributed: 86644\n",
      "\t  >> Reads not distributed (eg. no genuses above threshold): 31\n",
      "\t  >> Unclassified reads: 482100\n",
      "BRACKEN OUTPUT PRODUCED: results/kraken2_output/FH3.bracken.output\n",
      "PROGRAM END TIME: 01-17-2025 02:26:51\n",
      "  Bracken complete.\n",
      " >> Checking for Valid Options...\n",
      " >> Running Bracken \n",
      "      >> python src/est_abundance.py -i results/kraken2_output/FH4.kraken.report -o results/kraken2_output/FH4.bracken.output -k data/kraken2_db/standard-8/database100mers.kmer_distrib -l G -t 0\n",
      "PROGRAM START TIME: 01-17-2025 02:26:51\n",
      ">> Checking report file: results/kraken2_output/FH4.kraken.report\n",
      "BRACKEN SUMMARY (Kraken report: results/kraken2_output/FH4.kraken.report)\n",
      "    >>> Threshold: 0 \n",
      "    >>> Number of genuses in sample: 1394 \n",
      "\t  >> Number of genuses with reads > threshold: 1394 \n",
      "\t  >> Number of genuses with reads < threshold: 0 \n",
      "    >>> Total reads in sample: 1000000\n",
      "\t  >> Total reads kept at genuses level (reads > threshold): 629323\n",
      "\t  >> Total reads discarded (genuses reads < threshold): 0\n",
      "\t  >> Reads distributed: 101838\n",
      "\t  >> Reads not distributed (eg. no genuses above threshold): 1\n",
      "\t  >> Unclassified reads: 268838\n",
      "BRACKEN OUTPUT PRODUCED: results/kraken2_output/FH4.bracken.output\n",
      "PROGRAM END TIME: 01-17-2025 02:26:51\n",
      "  Bracken complete.\n"
     ]
    }
   ],
   "source": [
    "!for sample in $(ls data/metagenome_samples/*_1.fastq.gz | sed 's/.*\\///' | sed 's/_1.*//'); do \\\n",
    "    bracken -d data/kraken2_db/standard-8/ -l G -i results/kraken2_output/${sample}.kraken.report -r 100 -o results/kraken2_output/${sample}.bracken.output;\\\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "340360bb-5f33-481d-8b11-4c56ce9e5133",
   "metadata": {},
   "source": [
    "We could just proceed with the output file, but so that we can reuse the statistical analysis script used for metaphlan, we will convert bracken outout to match the metaphlan format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4e61ed19-e1c7-4ab9-9777-35ee2f05865c",
   "metadata": {},
   "outputs": [],
   "source": [
    "!for sample in $(ls data/metagenome_samples/*_1.fastq.gz | sed 's/.*\\///' | sed 's/_1.*//'); do \\\n",
    "    kreport2mpa.py -r results/kraken2_output/${sample}.kraken_bracken_genuses.report -o results/kraken2_output/${sample}_bracken.mpa; \\\n",
    "done "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "dc4dda83-233c-4eb9-b1f4-ac6809faeb50",
   "metadata": {},
   "outputs": [],
   "source": [
    "!for sample in $(ls data/metagenome_samples/*_1.fastq.gz | sed 's/.*\\///' | sed 's/_1.*//'); do \\\n",
    "    rm -f results/kraken2_output/${sample}_abundance.mpa;\\\n",
    "    sum=$(grep -vP \"\\\\|\" results/kraken2_output/${sample}_bracken.mpa | cut -f 2 | awk '{sum += $1} END {printf \"%.2f\", sum/100}'); \\\n",
    "    awk -v sum=\"$sum\" 'BEGIN {FS=\"\\t\"; OFS=\"\\t\"} {print $1, $2/sum}' results/kraken2_output/${sample}_bracken.mpa >> results/kraken2_output/${sample}_abundance.mpa;\\\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1005af33-bc89-44fb-8ae2-6be5d805e503",
   "metadata": {},
   "source": [
    "Merge "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "0cccf20f-2bcb-49f5-808f-fe1132a7022d",
   "metadata": {},
   "outputs": [],
   "source": [
    "! rm -f results/kraken2_output/merged_kraken_abundance_table.txt\n",
    "! merge_metaphlan_tables.py results/kraken2_output/*abundance.mpa >> results/kraken2_output/merged_kraken_abundance_table.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "17faf466-8ac9-4804-9d3a-ba17a49a369d",
   "metadata": {},
   "outputs": [],
   "source": [
    "! sed -i '1s/_abundance//g' results/kraken2_output/merged_kraken_abundance_table.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3432410f-9099-463e-ae94-432246a5cc22",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}