this is for holding javascript data
Ino de Bruijn Add kraken kmer origin counting
almost 10 years ago
Commit id: 150004ccbcfeff5a6d4190ec7de385c7691bce40
deletions | additions
diff --git a/notebooks/kraken.ipynb b/notebooks/kraken.ipynb
new file mode 100644
index 0000000..02d0c09
--- /dev/null
+++ b/notebooks/kraken.ipynb
...
{
"metadata": {
"name": "",
"signature": "sha256:4359cc7d2281c0abc72356e35eb921bdb3bf886f0c079d739c5949d5bdc64bd8"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Metagenomic Assembly - Assessing Chimericity in Contigs with Kraken\n",
"In de novo metagenomic assembly it is possible that different genomes end up getting mixed in one contig, because they are so similar. How likely is this to occur?\n",
"\n",
"[Kraken](https://github.com/DerrickWood/kraken) performs taxonomic classification using kmers. You can build a database of your own references in fasta format as long as they have a GI number. Here we will use Kraken to determine where the kmers in a wrongly assembled contig are coming from. \n",
"\n",
"Starting data:\n",
"- 52 reference genomes\n",
"- Velvet assembly with kmer=31 and no scaffolding using only Illumina paired reads that align to those 52 genomes\n",
"\n",
"## Installing Kraken\n",
"This went pretty smoothly. I just cloned the repository and ran the installation script. Then copied the binaries to my home dir which is in my ``$PATH`` variable.\n",
"\n",
"```\n",
"$ cd ~/src\n",
"$ git clone https://github.com/DerrickWood/kraken\n",
"$ ./install_kraken.sh ../kraken-install\n",
"$ cd ../kraken-install\n",
"$ cp kraken kraken-build ~/bin\n",
"$ kraken --version\n",
"Kraken version 0.10.5-beta\n",
"Copyright 2013-2014, Derrick Wood ([email protected])\n",
"```\n",
"\n",
"For your own database you need to install jellyfish 1, which luckily was already installed on our server:\n",
"```\n",
"$ module load bioinfo-tools\n",
"$ module load jellyfish/1.1.11\n",
"```\n",
"\n",
"## Building the database of the references\n",
"\n",
"### Getting GI numbers for the genomes you are interested in\n",
"The reference fasta file of the 52 genomes that I had unfortunately did not contain the [GI number](http://www.ncbi.nlm.nih.gov/Sitemap/sequenceIDs.html) of NCBI, which is necessary for Kraken. It did contain some fairly descriptive names that could be used to query NCBI. Therefore I used those names to retrieve a ``original_name,gi_number`` csv file:\n",
"\n",
"```\n",
"for g in $(grep '^>' ../Mircea_07102013_selected_refs.fasta | sed 's/>//')\n",
"do \n",
" specific_g=$(echo $g | sed 's/tokodaii/tokodaii_str_7/' | sed 's/SulfuriYO3AOP1/Sulfuri_YO3AOP1/' | sed 's/_/+/g')\n",
" echo $g,$(\n",
" curl -sg 'http://www.ncbi.nlm.nih.gov/nuccore/?term=%28'$specific_g'+%29+AND+%22ddbj+embl+genbank%22[Filter]' | \n",
" grep -o -e '\\(complete genome\\|, whole genome shotgun sequencing project\\|chromosome 1, complete sequence\\|S2, complete sequence\\)\\(.\\)\\{500\\}' | \n",
" grep -o --color=always -e 'nuccore/.*report=genbank' -m1 | \n",
" cut -f2 -d'/' | cut -f1 -d'?'\n",
" )\n",
"done > rename-fasta-ids.csv\n",
"```\n",
"\n",
"Kudos if you can decipher that! We iterate over all the reference names in the fastafile. The first line in the loop changes some reference names so they actually give a result while searching NCBI - no love for corner cases. Then all the names in the fasta file contained underscores instead of spaces so I changed that to a ``+`` which is the character NCBI uses for a space in a query. NCBI is queried using the result. A grep expression is used to search for complete genomes in the retrieved HTML - again a bunch of corner case killed a simple generalization - and the GI id is obtained after cutting out the other text.\n",
"\n",
"After that I used [this](https://github.com/inodb/biorhino-tools/blob/92a5da5719804261fb55aa6c3e22a20774af78a0/scripts/br-fasta-replace-ids.py) script to replace the fasta ids.\n",
"\n",
"```\n",
"br-fasta-replace-ids.py rename-fasta-ids.csv ../Mircea_07102013_selected_refs.fasta -d, > Mircea_07102013_selected_refs_gi.fasta\n",
"```\n",
"\n",
"### Actually performing the build commands\n",
"```\n",
"kraken-build --download-taxonomy --db Mircea_07102013_selected_refs\n",
"kraken-build --add-to-library Mircea_07102013_selected_refs_gi.fasta --db Mircea_07102013_selected_refs\n",
"kraken-build --build --kmer-len 31 --db Mircea_07102013_selected_refs\n",
"```\n",
"\n",
"## Classifying\n",
"After the database is built, classifying the assembly was easy and went really fast:\n",
"```\n",
"kraken --threads 16 ma-contigs-velvetnoscaf31.fa --db Mircea_07102013_selected_refs > ma-contigs-velvetnoscaf31-kraken.tsv\n",
"```\n",
"\n",
"The resulting file has the following format:\n",
"\n",
"```\n",
"OUTPUT FORMAT\n",
"=============\n",
"\n",
"Each sequence classified by Kraken results in a single line of\n",
"output. Output lines contain five tab-delimited fields; from\n",
"left to right, they are:\n",
"\n",
" 1) \"C\"/\"U\": one letter code indicating that the sequence was\n",
" either classified or unclassified.\n",
" 2) The sequence ID, obtained from the FASTA/FASTQ header.\n",
" 3) The taxonomy ID Kraken used to label the sequence; this is\n",
" 0 if the sequence is unclassified.\n",
" 4) The length of the sequence in bp.\n",
" 5) A space-delimited list indicating the LCA mapping of each k-mer\n",
" in the sequence. For example, \"562:13 561:4 A:31 0:1 562:3\"\n",
" would indicate that:\n",
" - the first 13 k-mers mapped to taxonomy ID #562\n",
" - the next 4 k-mers mapped to taxonomy ID #561\n",
" - the next 31 k-mers contained an ambiguous nucleotide\n",
" - the next k-mer was not in the database\n",
" - the last 3 k-mers mapped to taxonomy ID #562\n",
"```\n",
"\n",
"## The fun part\n",
"So now that we know for each contig where each kmers come from we can have some fun plotting the results.\n",
"\n",
"### Import Kraken taxonomy file for tax_id to rank mappings"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import pandas as pd"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def strip(text):\n",
" try:\n",
" return text.strip()\n",
" except AttributeError:\n",
" return text\n",
"\n",
"nodes = pd.read_csv(\"/media/milou/glob/src/kraken-install/databases/Mircea_07102013_selected_refs/taxonomy/nodes.dmp\", sep=\"|\", \n",
" skipinitialspace=True, index_col=False,\n",
" names=[\"tax_id\",\"parent_tax_id\",\"rank\",\"embl_code\",\n",
" \"division_id\",\"inherited_div_flag\",\"genetic_code_id\",\n",
" \"inherited_GC_flag\", \"mitochondrial genetic_code_id\",\n",
" \"inherited_MGC_flag\", \"Genbank_hidden_flag\",\"hidden_subtree_flag\",\n",
" \"comments\"],\n",
" converters={\"parent_tax_id\":strip,\"rank\":strip,\"embl_code\":strip})"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 60
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"nodes.head()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" | tax_id | parent_tax_id | rank | embl_code | division_id | inherited_div_flag | genetic_code_id | inherited_GC_flag | mitochondrial genetic_code_id | inherited_MGC_flag | Genbank_hidden_flag | hidden_subtree_flag | comments |
|---|
| 0 | 1 | 1 | no rank | | 8 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | \\t\\t |
|---|
| 1 | 2 | 131567 | superkingdom | | 0 | 0 | 11 | 0 | 0 | 0 | 0 | 0 | \\t\\t |
|---|
| 2 | 6 | 335928 | genus | | 0 | 1 | 11 | 1 | 0 | 1 | 0 | 0 | \\t\\t |
|---|
| 3 | 7 | 6 | species | AC | 0 | 1 | 11 | 1 | 0 | 1 | 1 | 0 | \\t\\t |
|---|
| 4 | 9 | 32199 | species | BA | 0 | 1 | 11 | 1 | 0 | 1 | 1 | 0 | \\t\\t |
|---|
\n",
"5 rows \u00d7 13 columns
\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" | classified | seq_id | tax_id | length | lca |
|---|
| 0 | True | NODE_1_length_4050_cov_8.365432 | 395495 | 4080 | [(395495, 4050)] |
|---|
| 1 | True | NODE_2_length_84_cov_2.071429 | 436114 | 114 | [(436114, 84)] |
|---|
| 2 | True | NODE_3_length_532_cov_4.377820 | 436114 | 562 | [(436114, 532)] |
|---|
| 3 | True | NODE_4_length_336_cov_9.851191 | 395495 | 366 | [(395495, 336)] |
|---|
| 4 | True | NODE_5_length_217_cov_7.797235 | 395495 | 247 | [(395495, 217)] |
|---|
\n",
"5 rows \u00d7 5 columns
\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" | class | family | forma | genus | infraclass | infraorder | kingdom | no rank | non-existent | order | parvorder | phylum | species | species group | species subgroup | strain | subclass | subfamily | subgenus | subkingdom | |
|---|
| 14316 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 54117 | 0 | 0 | 0 | 0 | ... |
|---|
| 33063 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 62 | 0 | 0 | 0 | 0 | 0 | 0 | 53541 | 0 | 0 | 0 | 0 | ... |
|---|
| 42327 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 90 | 0 | 0 | 0 | 0 | 0 | 0 | 66552 | 0 | 0 | 0 | 0 | ... |
|---|
| 43511 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 56783 | 0 | 0 | 0 | 0 | ... |
|---|
| 202248 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 61 | 0 | 0 | 0 | 0 | 0 | 0 | 60540 | 0 | 0 | 0 | 0 | ... |
|---|
\n",
"5 rows \u00d7 31 columns
\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" | strain | uknowit | ugotit | loveit | 5 |
|---|
| 0 | 0 | 0 | 0 | 0 | 0 |
|---|
\n",
"1 rows \u00d7 5 columns
\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
"\n",
"1 rows \u00d7 1 columns
\n",
"
