{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d015209b-0d23-4c85-994c-79b3535c75a3",
   "metadata": {},
   "source": [
    "# How to get some statistics about the SNP calling done ?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18a5bfca-2b96-49b1-bdc0-4d174f9f7c3e",
   "metadata": {},
   "source": [
    "## I. Get some basic stats from vcf files\n",
    "\n",
    "Be careful, the kernel of yhis jupyterbook is __python__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "304d9ab0-294c-4307-8cb7-a14154fb10cb",
   "metadata": {},
   "source": [
    "#### First go into the directory that contains the vcf file to analyze  </span> \n",
    "\n",
    "* List the content of this directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "701c59bb-2116-4af6-87fd-a63c73a96c60",
   "metadata": {},
   "outputs": [],
   "source": [
    "cd  /home/jovyan/SV_DATA/RESULTS/POP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0c0a6de-6fc7-48bf-ab3a-278590402d1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%ls /home/jovyan/SV_DATA/RESULTS/POP"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa5e1b2b-82f7-4710-b84f-06545badbfa9",
   "metadata": {},
   "source": [
    "#### Count the number of variants with `bcftools stat`\n",
    "\n",
    "* run the bcftools stats on the vcf file and save the result into the file `rice-CHR6.200000.vcf.stat`\n",
    "* check that the file have been correctly created and display the 35 first lines of this file\n",
    "* How many samples were used for this SNP analysis ?\n",
    "* How many SNPs were detected ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e005df0c-3ac1-4c29-930e-8c4cda3036a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "bcftools stats all_30samples.vcf >  all_30samples.bcfstat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d79f4a9-afaf-4b0a-9a4a-b3d46fcca1e0",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Display the first lines of the stat file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6350543e-0c09-42fa-a4ac-9bb93d03f711",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "head -n 35  all_30samples.bcfstat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74b53089-03fd-4d0a-a710-712acb2c8268",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Generating statistics from a VCF to determine how to set filters on it - `vcftools`\n",
    "\n",
    "We will generate more statistics from a VCF using vcftools  , a very useful and fast program for handling vcf files \n",
    "to easily calculate these statistics in order to better define filters we will apply and to get an idea of how to set such filtering thresholds. \n",
    "\n",
    "The main information we will consider are:\n",
    "* Depth: Usually, we filter SNP with a minimum and maximum depth. We use a minimum depth cutoffs to remove false positive calls and to keep higher quality calls too. \n",
    "A maximum cut off allow to remove regions with very, very high read depths such as repetitive regions.\n",
    "* Quality Genotype quality : With this filter, we should not trust any genotype with a Phred score below 20 which suggests a less than 99% accuracy.\n",
    "* Minor allele frequency MAF can cause big problems with SNP calls - and also inflate statistical estimates downstream. Ideally you want an idea of the distribution of your \n",
    "allelic frequencies but 0.05 to 0.10 is a reasonable cut-off. You should keep in mind however that some analyses, particularly demographic inference can be biased by MAF thresholds.\n",
    "* Missing data How much missing data are you willing to tolerate? It will depend on the study but typically any site with >25% missing data should be dropped.\n",
    "* biallelic, heterozygosity...\n",
    "\n",
    "In this training, we will just display quality distribution... but you should do on each value filterd.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9641f84-af45-4137-8f60-6016edd56596",
   "metadata": {},
   "source": [
    "#### Extracting quality per site `--site-quality`\n",
    "\n",
    "* run vcftool with the correct options\n",
    "* check that the files have been created and display the first lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d16f68ef-956a-4385-917b-ad8285965021",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "vcftools --vcf all_30samples.vcf  --site-quality --out qual"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf7e3619-df35-4f8a-a7f4-0c490d69a7b7",
   "metadata": {},
   "source": [
    "#### Display the first lines of the file just created by vcftools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "861692d3-23ed-4848-893a-77981a00fef2",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "ls -lrt\n",
    "head -n20  *qual"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d6e6f7d-37f0-48c7-994c-0f624d957a43",
   "metadata": {},
   "source": [
    "\n",
    "__Calculate allele frequency__\n",
    "\n",
    "* --freq2 : outputs the frequencies without information about the alleles\n",
    "* --freq would return their identity. \n",
    "* --max-alleles 2 to exclude sites that have more than two alleles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8d02f46-2fb0-462a-9077-49f1aa614f66",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "vcftools --bcf all_samples.bcf --freq --out AF --max-alleles 2\n",
    "vcftools --bcf all_samples.bcf --freq2 --out AF2 --max-alleles 2\n",
    "ls -lrt\n",
    "head *.frq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2de89286-b5c4-4c32-964e-45aca7ae12bc",
   "metadata": {},
   "source": [
    "### Generating density plot QUAL\n",
    "\n",
    "#### Plotting quality per site "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "307877db-c85c-47de-a16f-6c38ab9e3e70",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "\n",
    "qual_file=\"qual.lqual\"\n",
    "df_qual = pd.read_csv(qual_file, index_col=False, sep=\"\\t\")\n",
    "print(df_qual)\n",
    "df_qual['QUAL'].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3dc4b15-7c35-4b40-acb7-30aeeece268b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot with seaborn\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "sns.kdeplot(x=\"QUAL\", data=df_qual)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "788df901-4e83-4102-95c5-c097e6d72171",
   "metadata": {},
   "source": [
    "### III - FILTERING VCF "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ab0fde3-23e2-46c9-a54a-cd84c6960fc5",
   "metadata": {},
   "source": [
    "#### Which filters ?\n",
    "\n",
    "Here, we will apply on the vcf the following filters :\n",
    "\n",
    "* QUAL > 300\n",
    "* DP > 30 \n",
    "\n",
    "The threshodl of each filter depends on the SNP analysis (sample number, sequencing depth).\n",
    "\n",
    "In a second step, according the following analysys (eg: population genomics), usually, we could apply other filters such as :\n",
    "* removing missing data\n",
    "* keeping only biallellic\n",
    "* heterozygosity\n",
    "    \n",
    "#### Select only the SNPs (and remove the INDELs) \n",
    "* Filter vcf to keep only SNPs\n",
    "* Check that the new vcf has been created\n",
    "* Get the number of polymorphisms in the new vcf file "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52221eb2-8721-4547-b6a9-0b676bd6db3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "vcftools --vcf all_30samples.vcf --remove-indels --recode --recode-INFO-all --out all_30samples_onlySNPs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0e51c8f-ab8a-4f86-88a7-404d663552d5",
   "metadata": {},
   "source": [
    "#### Get the SNP count of this new vcf file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a91a9c9e-6069-4543-afff-df4fdb07a0a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "ls -lrth\n",
    "bcftools stats all_30samples_onlySNPs.recode.vcf | head -n35"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90e77c6a-3d4c-4d6e-bf30-60bfbe237765",
   "metadata": {},
   "source": [
    "#### Compress the vcf and generate the index of the compressed vcf - `tabix -p vcf vcf_file`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2af6a39-3eac-40cc-bc0d-cff8cfa5f35f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "bgzip all_30samples_onlySNPs.recode.vcf\n",
    "ls -lrth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3adc808-bde3-4baa-8afe-ce61140cedfb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57232bd3-aa17-44d5-839f-380cc1e653b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "tabix all_30samples_onlySNPs.recode.vcf.gz\n",
    "ls -lrt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67d62c00-4c7d-47d2-8334-69c78d78bbb3",
   "metadata": {},
   "source": [
    "#### Applying filters on QUAL\n",
    "\n",
    "Be careful, it miss one filter in the following command"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c360b424-da46-4580-8e06-d5dccf988025",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "#MAF=0.1\n",
    "#MISS=0.9\n",
    "QUAL=200\n",
    "#MIN_DEPTH=10\n",
    "#MAX_DEPTH=50\n",
    "\n",
    "vcftools --gzvcf all_30samples_onlySNPs.recode.vcf.gz \\\n",
    "--remove-indels --minQ $QUAL --recode --stdout | gzip -c > \\\n",
    "all_30samples_onlySNPs.recode.QUAL.vcf.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d3ce9eb-10ee-4f88-a402-4289f137c98c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "bcftools stats all_30samples_onlySNPs.recode.QUAL.vcf.gz | head -n 35"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "850f66c5-a78a-448d-bbc2-12a485c3d1bc",
   "metadata": {},
   "source": [
    "#### _Commands to apply other filters_\n",
    "\n",
    "##### remove missing data if more than 5% of missing data at a site (12 samples over 228 samples) `vcftools`\n",
    "\n",
    "`vcftools --vcf input.vcf --max-missing-count 12 --remove-filtered-all --recode --recode-INFO-all  --out ouput.onlySNP.12na`\n",
    "\n",
    "##### keep only biallelic sites\n",
    "`vcftools --vcf input.vcf --min-alleles 2 --max-alleles 2 --remove-filtered-all --recode --recode-INFO-all --out output.onlySNP.filteredPASS.na.225.Homoz.biallelic`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "358a403f-8c44-4853-9d32-037dda58b8cb",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}