module-8-variantcalling/Variant_calling_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ce25a930",
   "metadata": {},
   "source": [
    "### Variant calling module\n",
    "\n",
    "**CMM262, Winter 2022**\n",
    "\n",
    "Kyle Gaulton, kgaulton@gmail.com\n",
    "<br>\n",
    "<br>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a46ff4ea",
   "metadata": {},
   "source": [
    "In this walkthrough we will be calling, filtering and annotating genetic variants from a sequence alignment file\n",
    "\n",
    "**Required Files in resources/:**<br>\n",
    "*Human hg38 chr20 reference*<br>\n",
    "chr20.fa.gz \n",
    "chr20.dict,chr20.fa.fai,chr20.fa.gzi \n",
    "<br><br>\n",
    "*Variant call sets*<br>\n",
    "resources_broad_hg38_v0_1000G_phase1.snps.high_confidence.hg38.chr20.vcf.gz \n",
    "resources_broad_hg38_v0_1000G_phase1.snps.high_confidence.hg38.chr20.vcf.gz.tbi \n",
    "resources_broad_hg38_v0_1000G_omni2.5.hg38.vcf.gz \n",
    "resources_broad_hg38_v0_1000G_omni2.5.hg38.vcf.gz.tbi \n",
    "resources_broad_hg38_v0_hapmap_3.3.hg38.vcf.gz \n",
    "resources_broad_hg38_v0_hapmap_3.3.hg38.vcf.gz.tbi \n",
    "<br><br>\n",
    "*Annotation scripts*<br>\n",
    "annovar/\n",
    "   table_annovar.pl\n",
    "   annotate_variation.pl\n",
    "   humandb/*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36a95ac9",
   "metadata": {},
   "source": [
    "**Download and prepare .bam file for genotyping**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da6d56d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "/opt/conda/envs/r-bio/bin/samtools view -h -b ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000_genomes_project/data/GBR/HG00249/alignment/HG00249.alt_bwamem_GRCh38DH.20150718.GBR.low_coverage.cram chr20:30000000-30500000 > HG00249.bam\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a1b5f98",
   "metadata": {},
   "outputs": [],
   "source": [
    "# index bam file\n",
    "/opt/conda/envs/r-bio/bin/samtools index HG00249.bam\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77b11f44",
   "metadata": {},
   "outputs": [],
   "source": [
    "# list directory contents\n",
    "ls -la\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bee9747",
   "metadata": {},
   "outputs": [],
   "source": [
    "# view just the reads in the first 1000 bases of the file\n",
    "/opt/conda/envs/r-bio/bin/samtools view -h HG00249.bam chr20:30000000-30001000\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5da8dfa9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# summarize alignments\n",
    "/opt/conda/envs/r-bio/bin/samtools flagstat HG00249.bam\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9e040fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "## commands to fix and sort the alignment for duplicate marking/removal\n",
    "/opt/conda/envs/r-bio/bin/samtools sort -n -o HG00249.sort.bam HG00249.bam  ## commands to fix and sort the alignment for duplicate marking\n",
    "/opt/conda/envs/r-bio/bin/samtools fixmate -m HG00249.sort.bam HG00249.sort.fixed.bam\n",
    "/opt/conda/envs/r-bio/bin/samtools sort -o HG00249.resort.bam HG00249.sort.fixed.bam\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91edc95b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter alignments based on quality score of 30\n",
    "/opt/conda/envs/r-bio/bin/samtools view -b -q 30 -o HG00249.filter.bam HG00249.resort.bam\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "018512a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# index bam\n",
    "/opt/conda/envs/r-bio/bin/samtools index HG00249.filter.bam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3be10390",
   "metadata": {},
   "outputs": [],
   "source": [
    "# view stats for bam with quality filter\n",
    "/opt/conda/envs/r-bio/bin/samtools flagstat HG00249.filter.bam\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d306374c",
   "metadata": {},
   "outputs": [],
   "source": [
    "## remove duplicates (probably better to just mark them)\n",
    "/opt/conda/envs/r-bio/bin/samtools markdup -r HG00249.filter.bam HG00249.rmdup.bam\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "022b97b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "## index bam\n",
    "/opt/conda/envs/r-bio/bin/samtools index HG00249.rmdup.bam\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cc9f518",
   "metadata": {},
   "outputs": [],
   "source": [
    "## view stats for bam with duplicates removed\n",
    "/opt/conda/envs/r-bio/bin/samtools flagstat HG00249.rmdup.bam\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f02ea05b",
   "metadata": {},
   "outputs": [],
   "source": [
    "## view pileup of read counts per base\n",
    "/opt/conda/envs/r-bio/bin/samtools mpileup -f chr20.fa.gz HG00249.rmdup.bam\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f8141d7",
   "metadata": {},
   "source": [
    "**Call genetic variants with bcftools**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "939088ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "## call variants using bcftools and output to VCF\n",
    "/opt/conda/envs/py-bio/bin/bcftools mpileup -Ou -f chr20.fa.gz HG00249.rmdup.bam | /opt/conda/envs/py-bio/bin/bcftools call -mv -Ov -o HG00249.bcftools.vcf\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0016ac14",
   "metadata": {},
   "outputs": [],
   "source": [
    "## filter variant calls by quality score\n",
    "/opt/conda/envs/py-bio/bin/bcftools view -i '%QUAL>=20' HG00249.bcftools.vcf > HG00249.bcftools.filter.vcf\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d58f797b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# view top of file\n",
    "head -n 5000 HG00249.bcftools.filter.vcf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3df7a1d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# view summary of variant calls\n",
    "/opt/conda/envs/py-bio/bin/bcftools stats HG00249.bcftools.filter.vcf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bc4951f",
   "metadata": {},
   "source": [
    "**Call genetic variants using GATK**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90048563",
   "metadata": {},
   "outputs": [],
   "source": [
    "# see all of the tools available in GATK\n",
    "/opt/conda/envs/r-bio/bin/gatk --list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09e701fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# base recalibration\n",
    "/opt/conda/envs/r-bio/bin/gatk BaseRecalibrator -I HG00249.rmdup.bam -R chr20.fa.gz --known-sites resources_broad_hg38_v0_1000G_phase1.snps.high_confidence.hg38.chr20.vcf.gz -O recal.table\n",
    "/opt/conda/envs/r-bio/bin/gatk ApplyBQSR -R chr20.fa.gz -I HG00249.rmdup.bam --bqsr-recal-file recal.table -O HG00249.rmdup.recal.bam\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5736de1",
   "metadata": {},
   "outputs": [],
   "source": [
    "head -n 142 recal.table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7297998b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create initial call set using GATK HaplotypeCaller\n",
    "/opt/conda/envs/r-bio/bin/gatk HaplotypeCaller -I HG00249.rmdup.recal.bam -O HG00249.gatk.vcf -R chr20.fa.gz\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39bdef18",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check stats of initial VCF file\n",
    "/opt/conda/envs/py-bio/bin/bcftools stats HG00249.gatk.vcf\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "292134b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# variant quality recalibration and filtering\n",
    "/opt/conda/envs/r-bio/bin/gatk VariantRecalibrator -R chr20.fa.gz -V HG00249.gatk.vcf --resource:hapmap,known=false,training=true,truth=true,prior=15.0 resources_broad_hg38_v0_hapmap_3.3.hg38.vcf.gz --resource:omni,known=false,training=true,truth=false,prior=12.0 resources_broad_hg38_v0_1000G_omni2.5.hg38.vcf.gz --resource:1000G,known=false,training=true,truth=false,prior=10.0 resources_broad_hg38_v0_1000G_phase1.snps.high_confidence.hg38.chr20.vcf.gz -an QD -an FS -mode SNP -O recal.var --tranches-file output.tranches --rscript-file output.plots.R\n",
    "/opt/conda/envs/r-bio/bin/gatk ApplyVQSR -R chr20.fa.gz -V HG00249.gatk.vcf -O HG00249.gatk.filter.vcf --truth-sensitivity-filter-level 90.0 --tranches-file output.tranches --recal-file recal.var -mode SNP\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42c62d77",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check stats of filtered VCF file\n",
    "/opt/conda/envs/py-bio/bin/bcftools stats HG00249.gatk.filtered.vcf\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e19ab8e",
   "metadata": {},
   "source": [
    "**Annotate genetic variants**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7d2ec63",
   "metadata": {},
   "outputs": [],
   "source": [
    "# run variant annotation\n",
    "perl annovar/table_annovar.pl HG00249.bcftools.filter.vcf annovar/humandb/ -buildver hg38 -out HG00249 -remove -protocol refGene -operation g -nastring . -vcfinput\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b50fde8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Should produce a VCF with the annotations included and a text file of annotations\n",
    "ls -la *multianno*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb1c8e62",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pull out all variants in the promoter region of a gene\n",
    "grep 'upstream' HG00249.hg38_multianno.txt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "740af22d",
   "metadata": {},
   "source": [
    "**Convert genotypes to tab-delimited file**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75760030",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compress and tabix VCFs\n",
    "bgzip HG00249.bcftools.filter.vcf\n",
    "tabix -vcf HG00249.bcftools.filter.vcf.gz\n",
    "\n",
    "bgzip HG00249.gatk.filter.vcf\n",
    "tabix -vcf HG00249.gatk.filter.vcf.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "131fca32",
   "metadata": {},
   "outputs": [],
   "source": [
    "# output text file\n",
    "/opt/conda/envs/variant_calling/bin/vcf2tsv -g HG00249.gatk.filter.vcf.gz > HG00249.gatk.filter.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "925b5d44",
   "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
}