4snp_calling_GBSworkflow.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/usr/bin/sh: line 0: source: activate: file not found\r\n"
     ]
    }
   ],
   "source": [
    "!source activate py36"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "\n",
    "#sys.path.append('/home/faske/g/anaconda3/envs/py34/lib/python3.4/site-packages')\n",
    "sys.path.append('/data/gpfs/assoc/parchmanlab/tfaske/anaconda3/envs/py36/lib/python3.6/site-packages')\n",
    "sys.path.append(\"/data/gpfs/assoc/parchmanlab/tfaske/ipynb/include_utils\")\n",
    "\n",
    "import ipyparallel as ipp\n",
    "import os, time\n",
    "import include_utils as u\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import numbers\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as mpatches\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.colors as mcolors\n",
    "#import vcf\n",
    "from sklearn import preprocessing\n",
    "from subprocess import Popen, PIPE, call, check_output\n",
    "import seaborn as sns\n",
    "from IPython.display import FileLink\n",
    "import urllib.request as urllib2\n",
    "import dill\n",
    "import traceback\n",
    "from pandas import Series, DataFrame\n",
    "import gzip\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore',category=pd.io.pytables.PerformanceWarning)\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "from Bio import SeqIO\n",
    "#import pysam\n",
    "from collections import OrderedDict, namedtuple, Counter\n",
    "import operator\n",
    "import multiprocessing as mp\n",
    "import shutil\n",
    "import tempfile\n",
    "#from ipyparallel import Client\n",
    "import scandir\n",
    "import glob\n",
    "from Bio.SeqIO.QualityIO import FastqGeneralIterator\n",
    "import pickle\n",
    "import re\n",
    "from itertools import chain\n",
    "#import Levenshtein as lv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "root = \"/data/gpfs/home/tfaske/d/rabbit/full\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "snp_dir = os.path.join(root,'SNPcall')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {
    "ExecuteTime": {}
   },
   "outputs": [],
   "source": [
    "bwa_dir = os.path.join(snp_dir,'bwa')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/data/gpfs/assoc/denovo/tfaske/rabbit/full/SNPcall/bwa\n"
     ]
    }
   ],
   "source": [
    "cd $bwa_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "ExecuteTime": {}
   },
   "outputs": [],
   "source": [
    "bam_files = []\n",
    "files = !find . -type f -name '*sorted.bam'\n",
    "files = [os.path.abspath(x) for x in files if 'bam' in x]\n",
    "for x in files:\n",
    "    bam_files.append(x)\n",
    "bam_files = sorted(bam_files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "ExecuteTime": {}
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(608,\n",
       " '/data/gpfs/assoc/denovo/tfaske/rabbit/full/SNPcall/bwa/EN_AH_10_sorted.bam')"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(bam_files), bam_files[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {}
   },
   "outputs": [],
   "source": [
    "samtools = \"samtools\"\n",
    "bcftools = \"bcftools\"\n",
    "java = \"java\"\n",
    "perl = \"perl\"\n",
    "#perl = \"/home/cfriedline/gpfs/opt/ActivePerl-5.18/bin/perl\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {}
   },
   "outputs": [],
   "source": [
    "assembly = \"/data/gpfs/assoc/denovo/tfaske/rabbit/full/assembly/reference.fasta\"\n",
    "\n",
    "#indexed done by dDocent "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {}
   },
   "outputs": [],
   "source": [
    "!$samtools faidx {assembly}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculates coverage from bam_files\n",
    "also outputs a file called \"bam_coverage.csv\" in the bam_coverage folder  \n",
    "\n",
    "this will make a file for each sorted.bam and then consolidate in order to run efficiently\n",
    "\n",
    "first run to get size and remove empties: ls --sort=size -lh *sorted_bam "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/data/gpfs/assoc/denovo/tfaske/rabbit/full/SNPcall\n"
     ]
    }
   ],
   "source": [
    "cd $snp_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mkdir: cannot create directory ‘bam_coverage’: File exists\r\n"
     ]
    }
   ],
   "source": [
    "!mkdir bam_coverage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir bam_coverage/shdir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "### select options\n",
    "account = 'cpu-s5-denovo-0'\n",
    "partition = 'cpu-core-0'\n",
    "jobname = 'bam_cov'\n",
    "time = '1-00:00:00' #time limit 1 days\n",
    "cpus = 1\n",
    "mem_cpu = 2500\n",
    "email = 'tfaske@nevada.unr.edu'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "def write_bamcov_ind_sh(account,partition,time,jobname,cpus,mem_cpu,email,bam_files):\n",
    "    for i in range(0,len(bam_files)):\n",
    "        b = bam_files[i].split('/')[10].split('_sorted.bam')[0] #set this \n",
    "        with open(\"bam_coverage/shdir/run_bamcov_%s.sh\" % (b), \"w\") as o:\n",
    "            o.write(\"\"\"#!/usr/bin/env bash\n",
    "#SBATCH --account=%s\n",
    "#SBATCH --partition=%s\n",
    "#SBATCH --time=%s\n",
    "#SBATCH --ntasks 1\n",
    "#SBATCH --cpus-per-task %d\n",
    "#SBATCH --mem-per-cpu=%d\n",
    "#SBATCH --job-name %s_bamcov\n",
    "#SBATCH --output bam_coverage/%s_bamcov.txt\n",
    "#SBATCH --mail-user=%s \\n\\n\n",
    "    \n",
    "samtools depth -a %s\"\"\" % (account,partition,time,int(cpus),int(mem_cpu),b,b,email,bam_files[i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "write_bamcov_ind_sh(account,partition,time,jobname,cpus,mem_cpu,email,bam_files)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### finds all bamcov slurm scripts and writes bash script to sbatch them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "bamcov_files = []\n",
    "os.chdir('{}/{}'.format(root, 'SNPcall/bam_coverage/shdir'))\n",
    "files = !find . -name '*.sh'\n",
    "files = [os.path.abspath(x) for x in files]\n",
    "for x in files:\n",
    "        bamcov_files.append(x)\n",
    "bamcov_files = sorted(bamcov_files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(608,\n",
       " '/data/gpfs/assoc/denovo/tfaske/rabbit/full/SNPcall/bam_coverage/shdir/run_bamcov_EN_AH_1.sh')"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(bamcov_files),bamcov_files[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def write_bash_bamcov_sh(sh_files):\n",
    "    with open(\"SNPcall/run_bash_bamcov.sh\", \"w\") as o:\n",
    "        o.write(\"\"\"#!/usr/bin/env bash \\n\\n\"\"\")\n",
    "        for f in sh_files:\n",
    "            o.write(\"sbatch %s \\n\" % (f))    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/data/gpfs/assoc/denovo/tfaske/rabbit/full\n"
     ]
    }
   ],
   "source": [
    "cd $root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "write_bash_bamcov_sh(bamcov_files)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### run_bash_bamcov.sh\n",
    "    cd /data/gpfs/home/tfaske/d/rabbit/full/SNPcall \n",
    "    source activate py36\n",
    "    sbatch run_bash_bamcov.sh\n",
    "    \n",
    "### Read in all bamcov files and summarize "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "cov_files = []\n",
    "os.chdir('{}/{}'.format(root, 'SNPcall/bam_coverage'))\n",
    "files = !find . -name '*cov.txt'\n",
    "files = [os.path.abspath(x) for x in files]\n",
    "for x in files:\n",
    "        cov_files.append(x)\n",
    "cov_files = sorted(cov_files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(608,\n",
       " '/data/gpfs/assoc/denovo/tfaske/rabbit/full/SNPcall/bam_coverage/EN_AH_10_bamcov.txt')"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(cov_files),cov_files[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "__init__() got an unexpected keyword argument 'bam'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                           Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-77-40395681b620>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      9\u001b[0m     \u001b[0mcov_list\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcoverage\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m     \u001b[0;31m#print(float(coverage[0]),b)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m \u001b[0mcov_df\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbam\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mbam_names\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcoverage\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcov_list\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     12\u001b[0m \u001b[0mcov_df\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhead\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: __init__() got an unexpected keyword argument 'bam'"
     ]
    }
   ],
   "source": [
    "bam_names = []\n",
    "cov_list = []\n",
    "for i in range(0,len(cov_files)):\n",
    "    cov_file = cov_files[i]\n",
    "    c = cov_file.split('/')[10].split('_bamcov.txt')[0] #set this\n",
    "    b = c + '_sorted.bam'\n",
    "    coverage = !awk '{{ total += $$3 }} END {{ print total/NR }}' \"$cov_file\"\n",
    "    bam_names.append(b)\n",
    "    cov_list.append(float(coverage[0]))\n",
    "    #print(float(coverage[0]),b)\n",
    "cov_df = pd.DataFrame(bam=bam_names,coverage=cov_list)\n",
    "cov_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>bam</th>\n",
       "      <th>coverage</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>EN_AH_10_sorted.bam</td>\n",
       "      <td>32.7632</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>EN_AH_11_sorted.bam</td>\n",
       "      <td>30.5986</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>EN_AH_12_sorted.bam</td>\n",
       "      <td>35.1714</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>EN_AH_13_sorted.bam</td>\n",
       "      <td>33.0546</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>EN_AH_14_sorted.bam</td>\n",
       "      <td>22.6255</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   bam  coverage\n",
       "0  EN_AH_10_sorted.bam   32.7632\n",
       "1  EN_AH_11_sorted.bam   30.5986\n",
       "2  EN_AH_12_sorted.bam   35.1714\n",
       "3  EN_AH_13_sorted.bam   33.0546\n",
       "4  EN_AH_14_sorted.bam   22.6255"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cov_dict = {\"bam\":bam_names,'coverage':cov_list}\n",
    "cov_df = pd.DataFrame(cov_dict)\n",
    "cov_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "cov_out = os.path.join(snp_dir,'bam_coverage/bam_coverage.csv')\n",
    "cov_df.to_csv(path_or_buf=cov_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "count    608.000000\n",
       "mean      29.271217\n",
       "std        6.462445\n",
       "min        2.030350\n",
       "25%       26.070175\n",
       "50%       30.365850\n",
       "75%       33.688850\n",
       "max       46.151000\n",
       "Name: coverage, dtype: float64"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cov_df.coverage.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select good coverage bam files and move them to new dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9, 20, 10)"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(cov_df[cov_df.coverage < 10]),len(cov_df[cov_df.coverage < 15]),len(cov_df[cov_df.coverage < 11])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(599, 'EN_AH_10_sorted.bam')"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "good_bam = cov_df.bam[cov_df.coverage >= 10]\n",
    "good_bam = good_bam.tolist()\n",
    "len(good_bam),good_bam[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(599, '/data/gpfs/home/tfaske/d/rabbit/full/SNPcall/bwa/EN_AH_11_sorted.bam')"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "good_bam_files = []\n",
    "for i in range(0,len(good_bam)):\n",
    "    bam = bwa_dir + '/' + good_bam[i]\n",
    "    good_bam_files.append(bam)\n",
    "len(good_bam_files), good_bam_files[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/data/gpfs/assoc/denovo/tfaske/rabbit/full/SNPcall\n"
     ]
    }
   ],
   "source": [
    "cd $snp_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir 'good_bams'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "good_bam_dir = os.path.join(snp_dir,'good_bams')\n",
    "assert good_bam_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(0,len(good_bam_files)):\n",
    "    good_bam = str(good_bam_files[i])\n",
    "    !cp $good_bam $good_bam_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/data/gpfs/assoc/denovo/tfaske/rabbit/full/SNPcall/good_bams\n"
     ]
    }
   ],
   "source": [
    "cd $good_bam_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "good_bam_cp_files = []\n",
    "files = !find . -type f -name '*sorted.bam'\n",
    "files = [os.path.abspath(x) for x in files if 'bam' in x]\n",
    "for x in files:\n",
    "    good_bam_cp_files.append(x)\n",
    "good_bam_cp_files = sorted(good_bam_files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(599, '/data/gpfs/home/tfaske/d/rabbit/full/SNPcall/bwa/EN_AH_10_sorted.bam')"
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(good_bam_cp_files),good_bam_cp_files[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SNP CALL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/data/gpfs/assoc/denovo/tfaske/rabbit/full/SNPcall\n"
     ]
    }
   ],
   "source": [
    "cd $snp_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/data/gpfs/assoc/denovo/tfaske/rabbit/full/SNPcall'"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "### select options\n",
    "account = 'cpu-s5-denovo-0'\n",
    "partition = 'cpu-core-0'\n",
    "jobname = 'samtools_ERNA'\n",
    "time = '4-00:00:00' #time limit 3 day2\n",
    "cpus = 64\n",
    "mem_cpu = 2500\n",
    "email = 'tfaske@nevada.unr.edu'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "def write_callsnps_sh(account,partition,jobname,time,cpus,mem_cpu,email,assembly,good_bam_dir,snp_dir):\n",
    "    #print(account)\n",
    "    #print(partition)\n",
    "    #print(cpus)\n",
    "    #print(email)\n",
    "    #print(cmds[0])\n",
    "    with open(\"./run_callsnps.sh\", \"w\") as o:\n",
    "        o.write(\"\"\"#!/usr/bin/env bash\n",
    "#SBATCH --account=%s\n",
    "#SBATCH --partition=%s\n",
    "#SBATCH --time=%s\n",
    "#SBATCH --ntasks 1\n",
    "#SBATCH --cpus-per-task %d\n",
    "#SBATCH --mem-per-cpu=%d\n",
    "#SBATCH --job-name %s_callsnps\n",
    "#SBATCH --output output_callsnps.txt\n",
    "#SBATCH --mail-type=BEGIN,END\n",
    "#SBATCH --mail-user=%s \\n\\n\"\"\" % (account,partition,time,int(cpus),int(mem_cpu),jobname,email))\n",
    "\n",
    "        o.write(\"samtools mpileup -t DP,AD,ADF,ADR,SP,INFO/AD,INFO/ADF,INFO/ADR -Iugf %s %s/*sorted.bam | bcftools call -f GP,GQ -vmO z -o %s/samtools1-9.%s.vcf.gz\" %\n",
    "            (assembly,good_bam_dir,snp_dir,jobname))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [],
   "source": [
    "write_callsnps_sh(account,partition,jobname,time,cpus,mem_cpu,email,assembly,good_bam_dir,snp_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# run_call_snps.sh \n",
    "    cd /data/gpfs/home/tfaske/d/rabbit/full/SNPcall\n",
    "    source activate py36\n",
    "    sbatch run_callsnps.sh"
   ]
  },
  {
   "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}