{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Doomsayer Tutorial\n",
    "\n",
    "## Overview\n",
    "\n",
    "This tutorial notebook contains several demos of how to use Doomsayer.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Display doomsayer's help prompt\n",
    "\n",
    "To list all available options for Doomsayer, use the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python doomsayer.py -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download example VCFs+reference genome\n",
    "\n",
    "For the examples provided in this tutorial, you will need to download the following files from the Google Cloud Genomics storage repository\n",
    "\n",
    "- gzipped VCF files for Chr21 and Chr22 from the 1000 Genomes Phase 3 project\n",
    "- GRCh37-lite reference genome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://storage.googleapis.com/genomics-public-data/ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20130502/ALL.chr21.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz\n",
    "!wget https://storage.googleapis.com/genomics-public-data/ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20130502/ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz\n",
    "!wget https://storage.googleapis.com/genomics-public-data/references/GRCh37lite/GRCh37-lite.fa.gz\n",
    "!gunzip GRCh37-lite.fa.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Doomsayer on single VCF\n",
    "\n",
    "The default input for Doomsayer is a single vcf file. The following examples run Doomsayer separately on the 1000 Genomes VCFs for Chromosome 21 and Chromosome 22. To keep the output from both chromosomes, we use the `--projectdir \"1kg_chr[N]\"` option."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python doomsayer.py \\\n",
    "    --mode vcf \\\n",
    "    --input ALL.chr21.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz \\\n",
    "    --fastafile GRCh37-lite.fa \\\n",
    "    --projectdir \"1kg_chr21\" \\\n",
    "    --verbose"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python doomsayer.py \\\n",
    "    --mode vcf \\\n",
    "    --input ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz \\\n",
    "    --fastafile GRCh37-lite.fa \\\n",
    "    --projectdir \"1kg_chr22\" \\\n",
    "    --verbose"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Doomsayer in a pipeline\n",
    "\n",
    "In VCF mode, Doomsayer will accept input piped from other programs by specifying `--input -`, allowing the user to perform any number of pre-processing steps to the VCF file.\n",
    "\n",
    "In the example below, we are pre-filtering the Chr22 VCF with bcftools to include only biallelic SNVs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!bcftools view -m2 -M2 -v snps ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz | \\\n",
    "python doomsayer.py \\\n",
    "  --mode vcf \\\n",
    "  --input - \\\n",
    "  --fastafile GRCh37-lite.fa \\\n",
    "  --verbose"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Doomsayer on multiple VCFs\n",
    "\n",
    "In VCF mode, Doomsayer will also accept a text file containing the file paths of the VCFs we wish to process (one per line).\n",
    "\n",
    "This example runs with the `--cpus` option to enable parallel processing of the VCFs.\n",
    "\n",
    "(this is run without the `--verbose` flag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create input text file\n",
    "!ls *.vcf.gz > vcfs.txt\n",
    "\n",
    "# confirm text file has the VCFs we want\n",
    "!cat vcfs.txt\n",
    "\n",
    "# run doomsayer with VCF list as input\n",
    "!python doomsayer.py \\\n",
    "    --mode vcf \\\n",
    "    --input vcfs.txt \\\n",
    "    --fastafile GRCh37-lite.fa \\\n",
    "    --projectdir \"1kg_combined\" \\\n",
    "    --cpus 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Automatically create filtered output\n",
    "\n",
    "In VCF or text mode, when you specify the `--filterout` option, once the list of outliers has been generated, Doomsayer will create a new VCF with all of the outliers removed. This output is written to `STDOUT` to enable piping to other processing steps; to write to a file, simply add `> output.vcf` at the end of the command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python doomsayer.py \\\n",
    "  --mode vcf \\\n",
    "  --input ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz \\\n",
    "  --fastafile GRCh37-lite.fa \\\n",
    "  --verbose \\\n",
    "  --rank 4 \\\n",
    "  --filterout  > chr22.autofiltered.vcf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Doomsayer on an existing count matrix\n",
    "\n",
    "In some cases, you may wish to re-run the outlier detection or mutation signature analysis using different options (e.g., with a different rank). To avoid having to re-generate the input matrix, you can use Doomsayer in aggregation mode (`--mode agg`), and supply an existing matrix as the input. Here we are using the count matrix we previously generated from the Chr21 VCF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python doomsayer.py \\\n",
    "  --mode agg \\\n",
    "  --input \"1kg_chr21/NMF_M_spectra.txt\" \\\n",
    "  --verbose \\\n",
    "  --rank 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aggregation mode with multiple input matrices\n",
    "\n",
    "You may also encounter situations where you have generated count matrices on partitions of your data and need to recombine them. Doomsayer can aggregate data in two ways:\n",
    "\n",
    "- combining from runs on data split by genomic regions (but in identical samples)\n",
    "- combining from runs on data split by different samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aggregating output from different regions\n",
    "\n",
    "If you need to combine output from multiple previous Doomsayer runs on non-overlapping genomic regions, create a text file named `m_regions.txt` containing the paths to the count matrices (one per line), and use as the input file in aggregation mode (`--mode agg`).\n",
    "\n",
    "This mode is particularly useful if you are working with very large datasets (>10,000 samples), as you can split the input into small regions and spawn many simultaneous Doomsayer runs using Slurm or other workload managers. \n",
    "\n",
    "(though it is often easier to run in multiple VCF mode rather than writing output )\n",
    "\n",
    "Note that the input matrices **must** have identical dimensions. **Samples must be ordered identically in the rows of each matrix.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create text file containing paths of count matrices from the chr21 and chr22 output\n",
    "!find . -type d -name \"1kg_chr*\" -exec find {} -type f -name \"*spectra.txt\" \\; > m_regions.txt\n",
    "\n",
    "# confirm text file has the count matrices we want\n",
    "!cat m_regions.txt\n",
    "\n",
    "# if input file is named \"m_regions.txt,\" doomsayer will add matrices element-wise\n",
    "!python doomsayer.py \\\n",
    "  --mode agg \\\n",
    "  --input m_regions.txt \\\n",
    "  --verbose \\\n",
    "  --rank 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aggregating output from different subsamples\n",
    "\n",
    "In some instances, your data may be split into subsamples--Doomsayer can be run on these separately to generate the subsample count matrices, and then run again in aggregation mode to concatenate these matrices. To aggregate in this way, create a text file named `m_samples.txt` containing the paths to the subsample count matrices (one per line), and use as the input file in aggregation mode (`--mode agg`).\n",
    "\n",
    "We recommend caution when using this mode; data should usually be split into subsamples only after variant calling has been performed over the entire sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create text file containing paths of count matrices from the chr21 and chr22 output\n",
    "!find . -type d -name \"1kg_chr*\" -exec find {} -type f -name \"*spectra.txt\" \\; > m_samples.txt\n",
    "\n",
    "# confirm text file has the count matrices we want\n",
    "!cat m_samples.txt\n",
    "\n",
    "# if input file is named \"m_samples.txt,\" doomsayer will concatenate matrices rowwise\n",
    "# i.e., if matrix 1 is N1xK and matrix 2 is N2xK, the new matrix will be (N1+N2)xK\n",
    "!python doomsayer.py \\\n",
    "  --mode agg \\\n",
    "  --input m_samples.txt \\\n",
    "  --verbose \\\n",
    "  --rank 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate diagnostic report\n",
    "\n",
    "In this example, we run Doomsayer on an existing count matrix, previously generated from the entire 1000 Genomes data, and include the `--report` option to generate a diagnostic report in the output directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python doomsayer.py \\\n",
    "  --mode agg \\\n",
    "  --input tutorial_data/NMF_M_spectra.txt \\\n",
    "  --verbose \\\n",
    "  --rank 4 \\\n",
    "  --report"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[View Report](doomsayer_output/report.html)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}