{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <span style=\"color:gray\">ipyrad-analysis toolkit: </span> digest genomes\n",
    "\n",
    "The purpose of this tool is to digest a genome file *in silico* using the same restriction enzymes that were used for an empirical data set to attempt to extract homologous data from the genome file. This can be a useful procedure for adding additional outgroup samples to a data set. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Required software"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# conda install ipyrad -c conda-forge -c bioconda"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ipyrad.analysis as ipa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A genome file\n",
    "You will need a genome file in fasta format (optionally it can be gzip compressed). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "genome = \"/home/deren/Downloads/Ahypochondriacus_459_v2.0.fa\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialize the tool (e.g., ddRAD)\n",
    "\n",
    "You can generate single or paired-end data, and you will likely want to restrict the size of selected fragments to be within an expected size selection window, as is typically done in empirical data sets. Here I select all fragments occuring between two restriction enzymes where the intervening fragment is 300-500bp in length. I then ask that the analysis returns the digested fragments as 150bp fastq reads, and to provide 10 copies of each one. I also restrict it to only the first (largest) 12 scaffolds using the 'nscaffolds' arg."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "digest = ipa.digest_genome(\n",
    "    fasta=genome,\n",
    "    name=\"amaranthus-digest\",\n",
    "    workdir=\"digested_genomes\",\n",
    "    re1=\"CTGCAG\",\n",
    "    re2=\"AATTC\",\n",
    "    ncopies=10,\n",
    "    readlen=150,\n",
    "    min_size=300,\n",
    "    max_size=500,  \n",
    "    nscaffolds=12,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "extracted reads from 9058 positions\n"
     ]
    }
   ],
   "source": [
    "digest.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 4520\r\n",
      "-rw-r--r-- 1 deren deren  910924 Aug 30 12:14 amaranthus-digest_R1_.fastq.gz\r\n",
      "-rw-r--r-- 1 deren deren  909659 Aug 30 12:14 amaranthus-digest_R2_.fastq.gz\r\n",
      "-rw-r--r-- 1 deren deren 2793929 Aug 30 12:14 amaranthus-digest-RAD_R1_.fastq.gz\r\n",
      "-rw-r--r-- 1 deren deren      52 Aug 30 12:11 amaranthus-digest-RAD_R2_.fastq.gz\r\n"
     ]
    }
   ],
   "source": [
    "! ls -l digested_genomes/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2 (original RAD data)\n",
    "\n",
    "The original RAD method uses sonication rather than a second restriction digestion to cut all of the fragments down to an appropriate size for sequencing. Thus you only need to provide a single cut site and a selection window. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "digest = ipa.digest_genome(\n",
    "    fasta=genome,\n",
    "    name=\"amaranthus-digest-RAD\",\n",
    "    workdir=\"digested_genomes\",\n",
    "    re1=\"CTGCAG\",\n",
    "    re2=None,\n",
    "    paired=False,\n",
    "    ncopies=10,\n",
    "    readlen=100,\n",
    "    min_size=300,\n",
    "    max_size=500,  \n",
    "    nscaffolds=12,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "extracted reads from 27844 positions\n"
     ]
    }
   ],
   "source": [
    "digest.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}