newdocs/API-analysis/cookbook-digest_genomes.ipynb

{
 "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 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\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. "
   ]
  },
  {
   "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",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "fio = open(genome)\n",
    "scaffolds = fio.read().split(\">\")[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "ordered = sorted(scaffolds, key=lambda x: len(x), reverse=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "38760089"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(ordered[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "digest.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 23984\r\n",
      "-rw-r--r-- 1 deren 12270437 Dec 27 14:36 \u001b[0m\u001b[01;31mamaranthus-digest_R1_.fastq.gz\u001b[0m\r\n",
      "-rw-r--r-- 1 deren 12285598 Dec 27 14:36 \u001b[01;31mamaranthus-digest_R2_.fastq.gz\u001b[0m\r\n"
     ]
    }
   ],
   "source": [
    "ll digested_genomes/"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}