qc_plots.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Goal: Prototype some QC plots for use with MultiQC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "# Collect per-exon read counts from runs of the Dragen pipeline:\n",
    "mkdir -p data\n",
    "\n",
    "# From a recent run of our 218-gene GOAL-probe heme panel on 8 FFPE/Blood samples using KAPA LP+Cap kits:\n",
    "(echo -ne \"Chrom\\tStart\\tEnd\\tGene\\t\" && ls /mnt/pns/vcfs/220426_A01718_0009_AHWVGHDRXY/Heme/*{E,F}/*.target_bed_read_cov_report.bed | cut -f7 -d/ | perl -pe 's/.*(\\d)E$/S$1_50ng/; s/.*(\\d)F$/S$1_100ng/;' | paste -s) > data/goal_kapa_heme_exon_covg.txt\n",
    "ls /mnt/pns/vcfs/220426_A01718_0009_AHWVGHDRXY/Heme/*{E,F}/*.target_bed_read_cov_report.bed | xargs paste | cut -f1-4,$(seq -s, 6 8 62) | tail -n+2 >> data/goal_kapa_heme_exon_covg.txt\n",
    "\n",
    "# From a recent run of our 218-gene GOAL-probe heme panel on 8 FFPE/Blood samples using Twist LP+Cap kits:\n",
    "(echo -ne \"Chrom\\tStart\\tEnd\\tGene\\t\" && ls /mnt/pns/vcfs/220426_A01718_0009_AHWVGHDRXY/Heme/*{A,C}/*.target_bed_read_cov_report.bed | cut -f7 -d/ | perl -pe 's/.*(\\d)A$/S$1_50ng/; s/.*(\\d)C$/S$1_25ng/;' | paste -s) > data/goal_heme_exon_covg.txt\n",
    "ls /mnt/pns/vcfs/220426_A01718_0009_AHWVGHDRXY/Heme/*{A,C}/*.target_bed_read_cov_report.bed | xargs paste | cut -f1-4,$(seq -s, 6 8 62) | tail -n+2 >> data/goal_heme_exon_covg.txt\n",
    "\n",
    "# From a recent run of our 220-gene Twist-probe heme panel on 8 FFPE/Blood samples using Twist LP+Cap kits:\n",
    "(echo -ne \"Chrom\\tStart\\tEnd\\tGene\\t\" && ls /mnt/pns/vcfs/220915_A01718_0013_AHWVHFDRXY/Heme/*{A,C}/*.target_bed_read_cov_report.bed | cut -f7 -d/ | perl -pe 's/.*(\\d)A$/S$1_50ng/; s/.*(\\d)C$/S$1_25ng/;' | paste -s) > data/ucla_heme_exon_covg.txt\n",
    "ls /mnt/pns/vcfs/220915_A01718_0013_AHWVHFDRXY/Heme/*{A,C}/*.target_bed_read_cov_report.bed | xargs paste | cut -f1-4,$(seq -s, 6 8 62) | tail -n+2 >> data/ucla_heme_exon_covg.txt\n",
    "\n",
    "# From a recent run of the Twist Exome v2 panel on 8 Clinical blood samples:\n",
    "(echo -ne \"Chrom\\tStart\\tEnd\\tGene\\t\" && ls /mnt/pns/vcfs/220426_A01718_0009_AHWVGHDRXY/Exome/Exome*A/*.target_bed_read_cov_report.bed | cut -f7 -d/ | perl -pe 's/^Exome00(\\d)A$/S$1_50ng/' | paste -s) > data/twist_v2_exon_covg.txt\n",
    "ls /mnt/pns/vcfs/220426_A01718_0009_AHWVGHDRXY/Exome/Exome*A/*.target_bed_read_cov_report.bed | xargs paste | cut -f1-4,$(seq -s, 6 8 62) | tail -n+2 >> data/twist_v2_exon_covg.txt\n",
    "\n",
    "# Make a subset for the genes in the 220-gene heme panel:\n",
    "cut -f1 ~/src/heme-panel-design/data/exon_targets_gene_list.txt > data/ucla_heme_220.txt\n",
    "head -n1 data/twist_v2_exon_covg.txt > data/twist_v2_exon_covg_heme_genes.txt\n",
    "grep -wf data/ucla_heme_220.txt data/twist_v2_exon_covg.txt >> data/twist_v2_exon_covg_heme_genes.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load necessary packages and data\n",
    "import pandas as pd\n",
    "import plotly.express as px\n",
    "dfk = pd.read_csv(\"data/goal_kapa_heme_exon_covg.txt\", sep='\\t')\n",
    "dfg = pd.read_csv(\"data/goal_heme_exon_covg.txt\", sep='\\t')\n",
    "dfh = pd.read_csv(\"data/ucla_heme_exon_covg.txt\", sep='\\t')\n",
    "dfe = pd.read_csv(\"data/twist_v2_exon_covg_heme_genes.txt\", sep='\\t')\n",
    "\n",
    "# Convert read counts to read depth per exon, and then divide that by overall mean to get fold difference\n",
    "dfk.iloc[:, list(range(4,12))] = dfk.iloc[:, list(range(4,12))].apply(lambda x: 150*x/(dfk.End - dfk.Start))\n",
    "dfk.iloc[:, list(range(4,12))] = dfk.iloc[:, list(range(4,12))].apply(lambda x: x/x.mean())\n",
    "dfg.iloc[:, list(range(4,12))] = dfg.iloc[:, list(range(4,12))].apply(lambda x: 150*x/(dfg.End - dfg.Start))\n",
    "dfg.iloc[:, list(range(4,12))] = dfg.iloc[:, list(range(4,12))].apply(lambda x: x/x.mean())\n",
    "dfh.iloc[:, list(range(4,12))] = dfh.iloc[:, list(range(4,12))].apply(lambda x: 150*x/(dfh.End - dfh.Start))\n",
    "dfh.iloc[:, list(range(4,12))] = dfh.iloc[:, list(range(4,12))].apply(lambda x: x/x.mean())\n",
    "dfe.iloc[:, list(range(4,12))] = dfe.iloc[:, list(range(4,12))].apply(lambda x: 150*x/(dfe.End - dfe.Start))\n",
    "dfe.iloc[:, list(range(4,12))] = dfe.iloc[:, list(range(4,12))].apply(lambda x: x/x.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a violin plot:\n",
    "fig = px.violin(dfk, x=dfk.columns[list(range(4,12))], range_x=[-0.5,17], title=\"Exon depth uniformity of 218 heme genes using KAPA LP/Cap kits and GOAL probes\", labels={\"variable\": \"\", \"value\": \"Fold difference from mean of exon depths\"}, width=1200, height=1200, template=\"seaborn\", box=True, points=False)\n",
    "fig.add_vrect(x0=0.5, x1=2, opacity=0.3)\n",
    "fig.add_annotation(x=0.5, y=7.4, text=\"We want most exons within<br>0.5x-2x the mean exon depth\", showarrow=True, arrowhead=4)\n",
    "fig.update_layout(violingap=0, xaxis=dict(tickmode='linear', dtick=1), yaxis=dict(categoryorder=\"category descending\"))\n",
    "fig.show(renderer=\"png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a violin plot:\n",
    "fig = px.violin(dfg, x=dfg.columns[list(range(4,12))], range_x=[-0.5,17], title=\"Exon depth uniformity of 218 heme genes using Twist LP/Cap kits and GOAL probes\", labels={\"variable\": \"\", \"value\": \"Fold difference from mean of exon depths\"}, width=1200, height=1200, template=\"seaborn\", box=True, points=False)\n",
    "fig.add_vrect(x0=0.5, x1=2, opacity=0.3)\n",
    "fig.add_annotation(x=0.5, y=7.4, text=\"We want most exons within<br>0.5x-2x the mean exon depth\", showarrow=True, arrowhead=4)\n",
    "fig.update_layout(violingap=0, xaxis=dict(tickmode='linear', dtick=1), yaxis=dict(categoryorder=\"category descending\"))\n",
    "fig.show(renderer=\"png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a violin plot:\n",
    "fig = px.violin(dfh, x=dfh.columns[list(range(4,12))], range_x=[-0.5,17], title=\"Exon depth uniformity of 220 heme genes using Twist LP/Cap kits and Twist probes\", labels={\"variable\": \"\", \"value\": \"Fold difference from mean of exon depths\"}, width=1200, height=1200, template=\"seaborn\", box=True, points=False)\n",
    "fig.add_vrect(x0=0.5, x1=2, opacity=0.3)\n",
    "fig.add_annotation(x=0.5, y=7.4, text=\"We want most exons within<br>0.5x-2x the mean exon depth\", showarrow=True, arrowhead=4)\n",
    "fig.update_layout(violingap=0, xaxis=dict(tickmode='linear', dtick=1), yaxis=dict(categoryorder=\"category descending\"))\n",
    "fig.show(renderer=\"png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a violin plot:\n",
    "fig = px.violin(dfe, x=dfe.columns[list(range(4,12))], range_x=[-0.5,17], title=\"Exon depth uniformity of 220 heme genes using Twist Exome v2 on 8 CES samples\", labels={\"variable\": \"\", \"value\": \"Fold difference from mean of exon depths\"}, width=1200, height=1200, template=\"seaborn\", box=True, points=False)\n",
    "fig.add_vrect(x0=0.5, x1=2, opacity=0.3)\n",
    "fig.add_annotation(x=0.5, y=7.4, text=\"We want most exons within<br>0.5x-2x the mean exon depth\", showarrow=True, arrowhead=4)\n",
    "fig.update_layout(violingap=0, xaxis=dict(tickmode='linear', dtick=1), yaxis=dict(categoryorder=\"category descending\"))\n",
    "fig.show(renderer=\"png\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.5 ('bio')",
   "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.5"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "8e9fddb4e4315a775da86d9ed8228b3875b23df377a010a2ede7c3212891566c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}