work/4_ref_extraction/1_ref_extraction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract amplicon sequence data from a set of genomes using seekdeep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# panel data\n",
    "PANEL = '../../data/panel_extended_info.csv'\n",
    "# reference data\n",
    "GENOME_DIR = '../../../data/genome_vobs/'\n",
    "# primers file\n",
    "PRIMERS = os.path.join(GENOME_DIR, 'primers.txt')\n",
    "# output for extraction\n",
    "OUT_DIR = '../../../data/phylo_ampl/vobs_genomes'\n",
    "# extraction parameters\n",
    "READLEN = '150'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare files for extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "target\tforward\treverse\r\n",
      "0\tTGTSTACGGTCTGAAGAACATc\tTTATCCGGCTCCAAGTTAAGG\r\n",
      "1\tGAGCGtGCGGCcAAGATG\tACAgACCGACGTTAATGGC\r\n",
      "2\tCAGTCAAATTTCCAGACAATCT\tCGGAAGTGCATTTGAAGG-AAaA\r\n",
      "3\tGaTATAAATTGTCGATCACACAAACT\tTGCATTTATCGTAGTACAATCTCA\r\n",
      "4\tATGcTBGTCATgATGATGATCT\tCCGATCCACGATAAGGAGTAC\r\n",
      "5\tGCTGGCGCATAATTATCaCAAA\ttTTCCACTTCATCGCTCGC\r\n",
      "6\tGCAAAATTTCCGTCCCATTA\tTGTAATTAGCTGTGTCTTGTG\r\n",
      "7\tGTcTCgGAGCACATYGTG\tTCGTACTTCATTATTCTTTGGACTG\r\n",
      "8\tAGTGRCTCCAGACGGTgTT\tCCAAGGATTTGCTACTACCAcT\r\n"
     ]
    }
   ],
   "source": [
    "# primers file from panel\n",
    "# ambiguous chars and gaps are acceptable\n",
    "p = pd.read_csv(PANEL, dtype=str)[['Primary_ID', 'F', 'R']]\n",
    "p.columns = ['target','forward','reverse']\n",
    "p.to_csv(PRIMERS, sep='\\t', index=False)\n",
    "! head {PRIMERS}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "anopheles-albimanus-steclachromosomesaalbs2.fasta\r\n",
      "anopheles-aquasalis-A_aquasalis_v1.0-scaffolds.fasta\r\n",
      "anopheles-arabiensis-dongolascaffoldsaarad1.fasta\r\n",
      "anopheles-arabiensis-sharakhov.fasta\r\n",
      "anopheles-atroparvus-ebrochromosomesaatre3.fasta\r\n",
      "anopheles-christyi-achkn1017scaffoldsachra1.fasta\r\n",
      "anopheles-coluzzii-mali-nihscaffoldsacolm1.fasta\r\n",
      "anopheles-cracens-ASM209184v1-scaffolds.fasta\r\n",
      "anopheles-culicifacies-37scaffoldsacula1.fasta\r\n",
      "anopheles-darlingi-coariscaffoldsadarc3.fasta\r\n",
      "anopheles-dirus-wrair2scaffoldsadirw1.fasta\r\n",
      "anopheles-epiroticus-epiroticus2scaffoldsaepie1.fasta\r\n",
      "anopheles-farauti-far1scaffoldsafarf2.fasta\r\n",
      "anopheles-funestus-fumozchromosomesafunf3.fasta\r\n",
      "anopheles-funestus-fumozscaffoldsafunf1.fasta\r\n",
      "anopheles-gambiae-pestchromosomesagamp4.fasta\r\n",
      "anopheles-gambiae-pimperenascaffoldsagams1.fasta\r\n",
      "anopheles-koliensis-akwgs3contigsjxxb01.fasta\r\n",
      "anopheles-maculatus-maculatus3scaffoldsamacm1.fasta\r\n",
      "anopheles-melas-cm1001059ascaffoldsamelc2.fasta\r\n",
      "anopheles-merus-mafscaffoldsamerm2.fasta\r\n",
      "anopheles-minimus-minimus1scaffoldsaminm1.fasta\r\n",
      "anopheles-punctulatus-apwgs2contigsjxxa01.fasta\r\n",
      "anopheles-quadriannulatus-sangwescaffoldsaquas1.fasta\r\n",
      "anopheles-sinensis-chinascaffoldsasinc2.fasta\r\n",
      "anopheles-sinensis-sinensisscaffoldsasins2.fasta\r\n",
      "anopheles-stephensi-indianscaffoldsastei2.fasta\r\n",
      "anopheles-stephensi-sda-500scaffoldsastes1.fasta\r\n"
     ]
    }
   ],
   "source": [
    "# reference genomes have to be with fasta extension\n",
    "# naming of reference genome is then used for species inference\n",
    "! ls {GENOME_DIR} | grep '.fasta'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Version 2.6.4\r\n",
      "Programs\r\n",
      "Use SeekDeep [PROGRAM] --help to see more details about each program\r\n",
      "Commands are not case sensitive\r\n",
      "SeekDeep\r\n",
      "1) extractor\r\n",
      "2) extractorPairedEnd\r\n",
      "3) makeSampleDirectories\r\n",
      "4) processClusters\r\n",
      "5) qluster\r\n",
      "SeekDeepServer\r\n",
      "1) genProjectConfig\r\n",
      "2) popClusteringViewer\r\n",
      "SeekDeepUtils\r\n",
      "1) dryRunQualityFiltering\r\n",
      "2) genTargetInfoFromGenomes\r\n",
      "3) rBind\r\n",
      "4) replaceUnderscores\r\n",
      "5) runMultipleCommands\r\n",
      "6) setupTarAmpAnalysis\r\n"
     ]
    }
   ],
   "source": [
    "! SeekDeep --version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning, no sequences extracted for P1\r\n",
      "Warning, no sequences extracted for P2\r\n"
     ]
    }
   ],
   "source": [
    "# this takes quite some time\n",
    "! rm -rf {OUT_DIR}\n",
    "! SeekDeep genTargetInfoFromGenomes \\\n",
    "        --genomeDir {GENOME_DIR} \\\n",
    "        --primers {PRIMERS} \\\n",
    "        --numThreads 2 \\\n",
    "        --pairedEndLength {READLEN} \\\n",
    "        --dout {OUT_DIR} "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "no sequences for Plasmodium primers - that's a good sign\n",
    "\n",
    "## Output folder structure\n",
    "- `{target}` dir:\n",
    "    - `{target}`.fasta - per-genome sequences, same seqeunces for multiple genomes - reflected in header\n",
    "    - `{target}`_primersRemoved.fasta - same, without primers\n",
    "    - `extractionCounts.tab.txt` - correct counts for primer hits and extraction counts\n",
    "    - `genomeLocations/{genome}.bed` - per-genome amplicon coordinates\n",
    "    - `refAlignments/{genome}_[forward|reverse]Primer.sorted.bam` - bowtie2 alignments of primers\n",
    "    - `[forward|reverse]Primer.fasta` - all primers sequences resulting from ambiguous bases\n",
    "- `forSeekDeep` dir - additional input files for SeekDeep:\n",
    "    - `lenCutOffs.txt`: target minlen maxlen - from alignments to reference genomes\n",
    "    - `overlapStatuses.txt` target overlap status\n",
    "    - `refSeqs/{target}.fasta` same as `{target}/{target}.fasta`\n",
    "- `locationsByGenome/{genome}.bed` - all amplicon locations for each genome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "target\tstatus\r\n",
      "16\tNoOverLap\r\n",
      "28\tNoOverLap\r\n"
     ]
    }
   ],
   "source": [
    "# targets with unexpected overlap\n",
    "# NoOverLap is driven by one or more sequences being longer than 2 * READ_LEN\n",
    "! grep -i -v R1EndsInR2 {OUT_DIR}/forSeekDeep/overlapStatuses.txt"
   ]
  }
 ],
 "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
}