3.light_contigs_quality.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "67552fbb",
   "metadata": {},
   "source": [
    "\n",
    "# <span style=\"color:#006E7F\">__Introduction to Oxford Nanopore Data Analysis__ <a class=\"anchor\"></span>  \n",
    "\n",
    "\n",
    "Created by J. Orjuela (DIADE-IRD), F. Sabot (DIADE-IRD) and G. Sarah (AGAP-INRAE) - Septembre 2021 Formation SouthGreen\n",
    "\n",
    "Adapted by J. Orjuela (DIADE-IRD), F. Sabot (DIADE-IRD) - Novembre 2022\n",
    "    \n",
    "Adapted by J. Orjuela (DIADE-IRD) - mai 2023\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "689c1b2d",
   "metadata": {
    "tags": []
   },
   "source": [
    "# <span style=\"color:#006E7F\">__TP3 - CONTIGS QUALITY__ <a class=\"anchor\" id=\"data\"></span>  \n",
    "    \n",
    "# <span style=\"color: #4CACBC;\"> 1. Perform a comparison of assemblies using QUAST</span>  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a778bfc9",
   "metadata": {},
   "source": [
    "The metrics of assemblies can be evaluated using a quality assessment tool such as [QUAST](http://quast.bioinf.spbau.ru/manual.html).\n",
    "\n",
    "In part 2, we have assembled raw data using various assemblers and also polished and corrected it. \n",
    "\n",
    "Aggregate all assemblies you have produced for your sample in a new folder named \"AGGREGATED\".\n",
    "\n",
    "We're going to compare it with QUAST."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e448eb4f-2e13-4ada-969a-08772c33dcdc",
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir -p ~/work/RESULTS/AGGREGATED\n",
    "cd ~/work/RESULTS/AGGREGATED"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe6ac271-8998-45f7-b9b5-9ad9d6316a9b",
   "metadata": {},
   "source": [
    "### Order assemblies obtained in the samples "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb3f3d31-7cb4-4968-bb13-46d685040d34",
   "metadata": {},
   "source": [
    "#### sample 4222"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ab80ab6-6f8d-4d8e-849c-22952c64521c",
   "metadata": {},
   "outputs": [],
   "source": [
    "ln -s ~/work/RESULTS/4222_FLYE/assembly.fasta 4222_FLYE.fasta\n",
    "ln -s ~/work/RESULTS/4222_FLYE_MEDAKA/MEDAKA_CONSENSUS/consensus.fasta 4222_FLYE_MEDAKA.fasta\n",
    "ln -s ~/work/RESULTS/4222_RAVEN/assembly.fasta 4222_RAVEN.fasta\n",
    "ln -s ~/work/RESULTS/4222_RAVEN_MEDAKA/MEDAKA_CONSENSUS/consensus.fasta 4222_RAVEN_MEDAKA.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c3a0ae8-4607-4d47-997a-110964d21109",
   "metadata": {},
   "source": [
    "#### sample B8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76d945f6-7aac-4523-93f6-0587f409c7f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "ln -s ~/work/RESULTS/B8_FLYE/assembly.fasta B8_FLYE.fasta\n",
    "ln -s ~/work/RESULTS/B8_FLYE_MEDAKA/MEDAKA_CONSENSUS/consensus.fasta B8_FLYE_MEDAKA.fasta\n",
    "ln -s ~/work/RESULTS/B8_RAVEN/assembly.fasta B8_RAVEN.fasta\n",
    "ln -s ~/work/RESULTS/B8_RAVEN_MEDAKA/MEDAKA_CONSENSUS/consensus.fasta B8_RAVEN_MEDAKA.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68d9b023-7941-4196-a4ca-598685cf3702",
   "metadata": {},
   "source": [
    "#### sample G11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f253b4f7-6d62-49e5-afaf-b23329bb5592",
   "metadata": {},
   "outputs": [],
   "source": [
    "ln -s ~/work/RESULTS/G11_FLYE/assembly.fasta G11_FLYE.fasta\n",
    "ln -s ~/work/RESULTS/G11_FLYE_MEDAKA/MEDAKA_CONSENSUS/consensus.fasta G11_FLYE_MEDAKA.fasta\n",
    "ln -s ~/work/RESULTS/G11_RAVEN/assembly.fasta G11_RAVEN.fasta\n",
    "ln -s ~/work/RESULTS/G11_RAVEN_MEDAKA/MEDAKA_CONSENSUS/consensus.fasta G11_RAVEN_MEDAKA.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23b3e1ad-5a00-4ef4-a1ec-af4b64e5cb6f",
   "metadata": {},
   "source": [
    "#### check symbolic links"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14784ede-a468-49c9-9d07-4c56104401fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "ls -l ~/work/RESULTS/AGGREGATED"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a5ccbc0",
   "metadata": {},
   "source": [
    "### <span style=\"color: #4CACBD;\"> 1.1  Run QUAST on CLONE assemblies and compare them </span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "210bf799",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cd ~/work/RESULTS/AGGREGATED/\n",
    "time quast.py *.fasta -o QUAST"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f251f7c9",
   "metadata": {},
   "source": [
    "#### Go to the QUAST directory and check the file content.\n",
    "\n",
    "#### Looking at the output statistics, what are the main differences between the assemblies?\n",
    "\n",
    "#### Look total sizes, number of contigs, N50 statistics.\n",
    "\n",
    "#### For each sample what is the best one ?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e54cf5cb",
   "metadata": {},
   "source": [
    "### <span style=\"color: #4CACBD;\"> 1.2  Comparison of the assemblies against the reference genome of closely organism </span>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f0309a7",
   "metadata": {},
   "source": [
    "It's possible to compare assemblies against the reference genome of closely related organism. \n",
    "\n",
    "Use corrected assemblies and compare it with the reference genome"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d86ff5c7",
   "metadata": {},
   "source": [
    "### Compare best algae assemblies against the reference sequence using QUAST\n",
    "\n",
    "create for that a QUAST_REF directory into the AGGREGATED one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1ca734f-d663-4f2f-856a-ca6ea02af627",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "c482f4ba",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Which of the assemblies are closest in size relative to the reference ?\n",
    "\n",
    "#### Does this have the largest N50 and fewest number of contigs?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07ff361a",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### What assembly does obtain the best contiguity metrics?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30db29bd",
   "metadata": {
    "tags": []
   },
   "source": [
    "## <span style=\"color: #4CACBD;\"> 2. Assessing gene space using BUSCO </span>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae690225",
   "metadata": {},
   "source": [
    "Benchmarking Universal Single-Copy Orthologs ([BUSCO](https://busco.ezlab.org/busco_userguide.html)) helps to check if you have a good assembly, by searching the expected single-copy lineage-conserved orthologs in any newly-sequenced genome from an appropriate phylogenetic clade."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84bbbda5",
   "metadata": {},
   "source": [
    "Calculate the gene space for each assembly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed5fd510",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mkdir -p ~/work/RESULTS/AGGREGATED/BUSCO\n",
    "cd ~/work/RESULTS/AGGREGATED/BUSCO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40200c0c",
   "metadata": {},
   "source": [
    "Lineage can be choose in busco database using --list-datasets parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4851997d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# busco env\n",
    "conda activate busco"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d69a0fd-5301-4237-95d1-bc670dbdbfc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "busco --list-datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f870bc44",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ASSEMBLY=\"/home/jovyan/work/RESULTS/AGGREGATED/4222_FLYE_MEDAKA.fasta\"\n",
    "LINEAGE=chlorophyta_odb10\n",
    "busco -i \"$ASSEMBLY\" -l \"$LINEAGE\" -c \"${CPUS:-8}\" -m genome -o 4222_FLYE_MEDAKA_BUSCO \n",
    "# deactivate busco env\n",
    "conda deactivate"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b7b858a",
   "metadata": {
    "tags": []
   },
   "source": [
    "### What does the gene space look like for this assembly?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c9cbb00-d2d4-4d53-8523-20b9770e9c95",
   "metadata": {},
   "source": [
    "### Run BUSCO for the best assemblies ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38866ce2-d5de-415d-9b0c-c41c51082750",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "57067289-72d8-4d8a-8a1e-98b7ad023012",
   "metadata": {
    "tags": []
   },
   "source": [
    "## <span style=\"color: #4CACBD;\"> 3. Read alignment statistics - remapping</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e937007-8a83-473e-b589-56ac7cb303b1",
   "metadata": {},
   "source": [
    "Read congruency is an important measure in determining assembly accuracy. Clusters of read pairs that align incorrectly are strong indicators of mis-assembly.\n",
    "\n",
    "How well do the reads align back to the draft assemblies? Use minimap2 and samtools to assess the basic alignment statistics.\n",
    "\n",
    "Make a folder for your results.\n",
    "\n",
    "We will use 4222_FLYE_MEDAKA.fasta assembly and ONT READS as well ILLUMINA READS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "918a39b0-17cb-459b-8e32-375f46bf9187",
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir -p ~/work/RESULTS/REMAPPING\n",
    "cd ~/work/RESULTS/REMAPPING/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31f2f97f-8fa7-4fbe-8e37-2f195210c506",
   "metadata": {},
   "outputs": [],
   "source": [
    "ASSEMBLY=\"/home/jovyan/work/RESULTS/AGGREGATED/4222_FLYE_MEDAKA.fasta\"\n",
    "ONT=\"/home/jovyan/work/DATA/4222_RB2.fastq.gz\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a987f87d-f1ad-4caf-8188-7f99e13e732a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#symbolic link from last assembly in the current repertory REMAPPING\n",
    "ln -s ${ASSEMBLY} 4222_FLYE_MEDAKA.fasta\n",
    "# overwrite ASSEMBLY variable\n",
    "ASSEMBLY=4222_FLYE_MEDAKA.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3e56be3-f4ac-43a8-ada9-49d806ce3833",
   "metadata": {},
   "source": [
    "## Mapping assemblies vs ONT reads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eda47b7a-cf84-4e4a-83bc-9af06d621721",
   "metadata": {},
   "outputs": [],
   "source": [
    "minimap2 -ax map-ont -t 4 ${ASSEMBLY} ${ONT} | samtools sort -@ 1 -T \"${ASSEMBLY/.fasta/}\" -O BAM -o \"${ASSEMBLY/.fasta/_ONT_minimap2.bam}\" -\n",
    "samtools index \"${ASSEMBLY/.fasta/_ONT_minimap2.bam}\"\n",
    "samtools view -F 0x904 -c \"${ASSEMBLY/.fasta/_ONT_minimap2.bam}\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff9bb810-a6af-4915-9dc9-233b090ab837",
   "metadata": {},
   "source": [
    "### What is the percentage of aligned ONT reads in your clone assembly?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "211c4705-669d-463c-b084-14c9006ec4ec",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "2b1ca487-8991-440b-8c71-daea695edff9",
   "metadata": {},
   "source": [
    "## <span style=\"color: #4CACBD;\"> [OPTIONAL] Blobtools </span>\n",
    "\n",
    "During the sequence quality assessment stage we tried to discern whether contamination was present. Sometimes this is not feasible at the read level. By plotting Contig GC content vs Contig Read Coverage we can look for clusters of contigs that share similar coverage. The appearance of multiple clusters can indicate multiple organisms. Occasionally, contigs can also be taxonomically classified, providing further evidence for contaminants.\n",
    "\n",
    "First we need to download some files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05ee6194-ff4b-4ab6-8461-501155974f51",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cd ~/work/DATA\n",
    "wget --no-check-certificat -rm -nH --cut-dirs=1 --reject=\"index.html*\" https://itrop.ird.fr/ont-training/blobtools.tar.gz\n",
    "tar zxvf blobtools.tar.gz\n",
    "wget --no-check-certificat -rm -nH --cut-dirs=1 --reject=\"index.html*\" https://itrop.ird.fr/ont-training/testBacteria.dmnd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23b91930-8a7f-41b7-9c52-0c300ee87be7",
   "metadata": {},
   "source": [
    "Run Blobtools on each assembly. Blobtools requires both a BAM file as input and blast output for the classification step.2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afe2d3cd-f298-4add-b62d-9e32fbd52a52",
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir -p ~/work/RESULTS/BLOBTOOLS\n",
    "cd ~/work/RESULTS/BLOBTOOLS"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f723a846-7ea6-4793-9f53-1d7370d92500",
   "metadata": {},
   "source": [
    "### Blastx using diamond "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ed88289-4d5c-4e66-a334-6626718ff803",
   "metadata": {},
   "source": [
    "Run diamond in blastx mode using assembled contigs vs a pre-formatted diamond bacteria database (protein)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84a9cde6-448f-4c2e-a8a4-fd82ff8f952a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#symbolic link from last assembly in the current repertory BLOBTOOLs\n",
    "ln -s ~/work/RESULTS/4222_FLYE_MEDAKA/MEDAKA_CONSENSUS/consensus.fasta 4222_FLYE_MEDAKA.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "798fd16e-2d48-4fb7-ae16-a1650b1850cd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "time diamond blastx --query 4222_FLYE_MEDAKA.fasta --db ~/work/DATA/testBacteria.dmnd --outfmt 6 qseqid staxids bitscore qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore --sensitive  --max-target-seqs 1  --evalue 1e-25  --threads 4  --out diamond.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99478f50-b7c9-408d-94e6-7cb3d96e203e",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "head diamond.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84e9d3ba-9d8b-4785-b888-a46b15895a4e",
   "metadata": {},
   "source": [
    "### Run blobtools "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ec54503-62b8-4d38-be45-01a98d6065d7",
   "metadata": {},
   "source": [
    "Specifically to Blobtools: * nodes and names from ncbi taxdump database can be download from [here](https://github.com/DRL/blobtools#download-ncbi-taxdump-and-create-nodesdb).\n",
    "\n",
    "In this training nodes and names are available on the DATA directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8200d28-5b8d-4839-b5cb-97d9cf9a6fb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ASSEMBLY=4222_FLYE_MEDAKA.fasta\n",
    "BAM=~/work/RESULTS/REMAPPING/${ASSEMBLY/.fasta/_ONT_minimap2.bam}\n",
    "DIAMONDX=~/work/RESULTS/BLOBTOOLS/diamond.csv\n",
    "BLOB_NODES=~/work/DATA/blobtools/nodes.dmp\n",
    "BLOB_NAMES=~/work/DATA/blobtools/names.dmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7c7b83b-5ece-4234-9444-cd281a6543c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "blobtools create -i ${ASSEMBLY} -b ${BAM} -t ${DIAMONDX} -o quality --names ${BLOB_NAMES} --nodes ${BLOB_NODES} --db nodesDB.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89243aaf-04c6-475c-8a7a-620e9213b914",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "blobtools view -i quality.blobDB.json --cov -o output;\n",
    "blobtools plot -i quality.blobDB.json;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de79c4e2-50dd-4055-aaa7-8265ac6d5db5",
   "metadata": {},
   "source": [
    "#### Is there contamination in the assembly?\n",
    "\n",
    "#### Do any assemblies show strange clustering?\n",
    "\n",
    "#### Why might coverage vary across contigs within an assembly?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4254dab5-6412-4a6d-a02a-0849097463f7",
   "metadata": {},
   "source": [
    "The Blobplots all indicate a single cluster. Some contigs show fairly high coverage in comparison to the rest of the genome which could be repetitive elements in the genome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd77a472-d13d-4291-b40f-53a2511362c8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cee151e-64a5-4e95-8bab-457d5cd7f573",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f12a4805-1f16-4d43-82e7-2d6e17bae2af",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Bash",
   "language": "bash",
   "name": "bash"
  },
  "language_info": {
   "codemirror_mode": "shell",
   "file_extension": ".sh",
   "mimetype": "text/x-sh",
   "name": "bash"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}