easy_config_prep.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import glob\n",
    "import os\n",
    "import yaml\n",
    "import re\n",
    "import copy\n",
    "from io import StringIO\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def open_sample_sheet(sample_sheet_fp, lanes=False, skiprows=18):\n",
    "    \"\"\"Read in an IGM sample sheet and return a pandas DF with primary data table\"\"\"\n",
    "    sample_sheet = pd.read_excel(sample_sheet_fp, skiprows = skiprows, header=1)\n",
    "    \n",
    "    # Remove trailing whitespace from column names and homogenize case\n",
    "    sample_sheet.columns = [x.strip().lower() for x in sample_sheet.columns]\n",
    "    \n",
    "    if lanes:\n",
    "        sample_sheet = sample_sheet.loc[sample_sheet['lane'].isin(lanes)]\n",
    "    if 'sample name' not in sample_sheet.columns and 'sample_id' in sample_sheet.columns:\n",
    "        sample_sheet['sample name'] = sample_sheet['sample_id']\n",
    "\n",
    "    return(sample_sheet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def open_sequencing_manifest(manifest_fp, lanes='False', sheetname='Sample Information', skiprows=18):\n",
    "    \"\"\"Read in an IGM sequencing manifest and return a pandas DF with primary data table\"\"\"\n",
    "    sample_sheet = pd.read_excel(manifest_fp, sheetname=sheetname, skiprows=skiprows, header=1)\n",
    "    \n",
    "    # Remove trailing whitespace from column names and homogenize case\n",
    "    sample_sheet.columns = [x.strip().lower() for x in sample_sheet.columns]\n",
    "    \n",
    "    if lanes:\n",
    "        sample_sheet = sample_sheet.loc[sample_sheet['lane'].isin(lanes)]\n",
    "        \n",
    "    return(sample_sheet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_read(sample, seq_dir, read):\n",
    "    \"\"\"Function to pull a given read based on sample name from the reads directory\"\"\"\n",
    "    reads = glob.glob(os.path.join(seq_dir, \"{0}_*_{1}_*.fastq.gz\".format(sample, read)))\n",
    "    if len(reads) == 1:\n",
    "        return(reads[0])\n",
    "    elif len(reads) > 1:\n",
    "        raise ValueError('Too many reads found for {0} in {1}:\\n'\n",
    "                         'read_str: {2}'.format(sample, seq_dir, reads))\n",
    "    elif len(reads) < 1:\n",
    "        raise ValueError('Too few reads found for {0} in {1}:\\n'\n",
    "                         'read_str: {1}\\n'.format(sample, seq_dir, reads))   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_sample_dict(sample_sheet, seq_dir,\n",
    "                        forward = 'R1',\n",
    "                        reverse = 'R2',\n",
    "                        adaptor = '$CONDA_ENV_PATH/share/trimmomatic-*/adapters/TruSeq3-PE-2.fa',\n",
    "                        phred = 'phred33',\n",
    "                        sample_header = 'Sample_Prefix',\n",
    "                        sample_name = 'Sample',\n",
    "                        tolerate_losses = True):\n",
    "    \n",
    "    samples_pe = {'samples_pe': {}}\n",
    "    for x in sample_sheet.index:\n",
    "        try:\n",
    "            samples_pe['samples_pe'][sample_sheet.loc[x, sample_name]] = {\n",
    "               'forward': get_read(sample_sheet.loc[x, sample_header], seq_dir, forward),\n",
    "               'reverse': get_read(sample_sheet.loc[x, sample_header], seq_dir, reverse),\n",
    "               'adaptor': adaptor,\n",
    "               'phred': phred\n",
    "            }\n",
    "        except ValueError:\n",
    "            if tolerate_losses:\n",
    "                print('No sequences found for %s' % sample_sheet.loc[x, sample_name])\n",
    "            else:\n",
    "                raise ValueError\n",
    "    \n",
    "    return(samples_pe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def format_yaml(RUN, samples_pe,\n",
    "                   TMP_DIR_ROOT = '/localscratch',\n",
    "                   host_db = '/home/jgsanders/ref_data/genomes/Homo_sapiens_Bowtie2_v0.1/Homo_sapiens',\n",
    "                   envs = {'KNEAD_ENV': 'source activate kneaddata',\n",
    "                           'BOWTIE_ENV': 'source activate kneaddata',\n",
    "                           'QC_ENV': 'source activate kneaddata',\n",
    "                           'TRIM_ENV': 'source activate kneaddata',\n",
    "                           'HUMANN2_ENV': 'source activate kneaddata',\n",
    "                           'METAPHLAN_ENV': 'source activate kneaddata'},\n",
    "                   software = {'gzip': 'gzip',\n",
    "                               'trimmomatic': 'trimmomatic',\n",
    "                               'mash': '/home/jgsanders/git_sw/git_bin/mash-Linux64-v1.1.1/mash',\n",
    "                               'seqtk': '/home/jgsanders/git_sw/seqtk/seqtk'},\n",
    "                   params = {'TRIMMOMATIC': {'QUAL': 'LEADING:20 TRAILING:20 AVGQUAL:30 MINLEN:32 TOPHRED33',\n",
    "                                              'ILLUMINACLIP': '2:30:10'},\n",
    "                             'MASH': {'REFSEQ_DB': '/home/jgsanders/ref_data/mash/RefSeqSketchesDefaults.msh',\n",
    "                                      'OTHER': '-r -m 2 -k 21 -s 1000'},\n",
    "                             'HUMANN2': {'METAPHLAN_DIR': '/home/jgsanders/share/metaphlan2',\n",
    "                                         'HUMANN2_NT_DB': '/home/jgsanders/ref_data/humann2/chocophlan',\n",
    "                                         'HUMANN2_AA_DB': '/home/jgsanders/ref_data/humann2/uniref',\n",
    "                                         'NORMS': ['cpm','relab'],\n",
    "                                         'OTHER': ''}},\n",
    "                   default_flow_style = False):\n",
    "    \n",
    "    config_str = ''\n",
    "\n",
    "    config_str += yaml.dump({'TMP_DIR_ROOT': TMP_DIR_ROOT}, default_flow_style = default_flow_style)\n",
    "\n",
    "    config_str += yaml.dump({'RUN': RUN}, default_flow_style = default_flow_style)\n",
    "    \n",
    "    config_str += yaml.dump({'HOST_DB': host_db}, default_flow_style = default_flow_style)\n",
    "    \n",
    "    config_str += yaml.dump({'ENVS': envs}, default_flow_style = default_flow_style)            \n",
    "    \n",
    "    config_str += yaml.dump({'SOFTWARE': software}, default_flow_style = default_flow_style)\n",
    "\n",
    "    config_str += yaml.dump({'PARAMS': params}, default_flow_style = default_flow_style)\n",
    "\n",
    "    config_str += yaml.dump(samples_pe, default_flow_style = default_flow_style)\n",
    "    \n",
    "    return(config_str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fps = ['Bs_S106_L001_R1_001.fastq.gz',\n",
    "       'Bs_S106_L001_R2_001.fastq.gz',\n",
    "       'Bs_S106_L002_R1_001.fastq.gz',\n",
    "       'Bs_S106_L002_R2_001.fastq.gz',\n",
    "       'Vf_S104_L001_R1_001.fastq.gz',\n",
    "       'Vf_S104_L001_R2_001.fastq.gz',\n",
    "       'Vf_S104_L002_R1_001.fastq.gz',\n",
    "       'Vf_S104_L002_R2_001.fastq.gz']\n",
    "\n",
    "exp_df = pd.DataFrame.from_dict({'Extension': {0: 'fastq.gz',\n",
    "  1: 'fastq.gz',\n",
    "  2: 'fastq.gz',\n",
    "  3: 'fastq.gz',\n",
    "  4: 'fastq.gz',\n",
    "  5: 'fastq.gz',\n",
    "  6: 'fastq.gz',\n",
    "  7: 'fastq.gz'},\n",
    " 'File': {0: 'Bs_S106_L001_R1_001.fastq.gz',\n",
    "  1: 'Bs_S106_L001_R2_001.fastq.gz',\n",
    "  2: 'Bs_S106_L002_R1_001.fastq.gz',\n",
    "  3: 'Bs_S106_L002_R2_001.fastq.gz',\n",
    "  4: 'Vf_S104_L001_R1_001.fastq.gz',\n",
    "  5: 'Vf_S104_L001_R2_001.fastq.gz',\n",
    "  6: 'Vf_S104_L002_R1_001.fastq.gz',\n",
    "  7: 'Vf_S104_L002_R2_001.fastq.gz'},\n",
    " 'Index': {0: 'S106',\n",
    "  1: 'S106',\n",
    "  2: 'S106',\n",
    "  3: 'S106',\n",
    "  4: 'S104',\n",
    "  5: 'S104',\n",
    "  6: 'S104',\n",
    "  7: 'S104'},\n",
    " 'Lane': {0: 'L001',\n",
    "  1: 'L001',\n",
    "  2: 'L002',\n",
    "  3: 'L002',\n",
    "  4: 'L001',\n",
    "  5: 'L001',\n",
    "  6: 'L002',\n",
    "  7: 'L002'},\n",
    " 'Read': {0: 'R1',\n",
    "  1: 'R2',\n",
    "  2: 'R1',\n",
    "  3: 'R2',\n",
    "  4: 'R1',\n",
    "  5: 'R2',\n",
    "  6: 'R1',\n",
    "  7: 'R2'},\n",
    " 'Run': {0: '001',\n",
    "  1: '001',\n",
    "  2: '001',\n",
    "  3: '001',\n",
    "  4: '001',\n",
    "  5: '001',\n",
    "  6: '001',\n",
    "  7: '001'},\n",
    " 'Sample': {0: 'Bs',\n",
    "  1: 'Bs',\n",
    "  2: 'Bs',\n",
    "  3: 'Bs',\n",
    "  4: 'Vf',\n",
    "  5: 'Vf',\n",
    "  6: 'Vf',\n",
    "  7: 'Vf'}})\n",
    "\n",
    "def parse_ilm_fps_to_df(fps,\n",
    "                        pattern = '^((.+?)_(S\\d+)_(L\\d+)_(R[12])_(\\d+)\\.(.+))$',\n",
    "                        pattern_names = ['File','Sample','Index','Lane','Read','Run','Extension']):\n",
    "    \n",
    "    p = re.compile(pattern)\n",
    "\n",
    "    df = pd.DataFrame(columns = pattern_names)\n",
    "\n",
    "    for f in fps:\n",
    "        m = p.match(f)\n",
    "\n",
    "        if m:\n",
    "            df = df.append(dict(zip(pattern_names, m.groups())), ignore_index = True)\n",
    "\n",
    "    return(df)\n",
    "\n",
    "obs_df = parse_ilm_fps_to_df(fps)\n",
    "pd.util.testing.assert_frame_equal(obs_df.sort_index(axis=1), exp_df.sort_index(axis=1))\n",
    "\n",
    "df = exp_df\n",
    "seq_dir = './example/reads'\n",
    "exp_dict = {'Bs': {'forward': ['./example/reads/Bs_S106_L001_R1_001.fastq.gz',\n",
    "                               './example/reads/Bs_S106_L002_R1_001.fastq.gz'],\n",
    "                   'reverse': ['./example/reads/Bs_S106_L001_R2_001.fastq.gz',\n",
    "                               './example/reads/Bs_S106_L002_R2_001.fastq.gz']},\n",
    "            'Vf': {'forward': ['./example/reads/Vf_S104_L001_R1_001.fastq.gz',\n",
    "                               './example/reads/Vf_S104_L002_R1_001.fastq.gz'],\n",
    "                   'reverse': ['./example/reads/Vf_S104_L001_R2_001.fastq.gz',\n",
    "                               './example/reads/Vf_S104_L002_R2_001.fastq.gz']}}\n",
    "\n",
    "def get_sample_reads_df(df, seq_dir):\n",
    "    sample_reads_dict = {}\n",
    "    \n",
    "    samples = list(set(df['Sample']))\n",
    "    \n",
    "    for s in samples:\n",
    "        f_fps = sorted(list(df.loc[(df['Sample'] == s) & (df['Read'] == 'R1'),'File'].values))\n",
    "        r_fps = sorted(list(df.loc[(df['Sample'] == s) & (df['Read'] == 'R2'),'File'].values))\n",
    "\n",
    "        sample_reads_dict[s] = {'forward': [os.path.join(seq_dir, x) for x in f_fps],\n",
    "                                'reverse': [os.path.join(seq_dir, x) for x in r_fps]}\n",
    "    \n",
    "    return(sample_reads_dict)\n",
    "\n",
    "obs_dict = get_sample_reads_df(df, seq_dir)\n",
    "assert(exp_dict == obs_dict)\n",
    "\n",
    "def get_sample_paths(seq_dir):\n",
    "    fps = os.listdir(seq_dir)\n",
    "    \n",
    "    files_df = parse_ilm_fps_to_df(fps)\n",
    "    sample_reads_dict = get_sample_reads_df(files_df, seq_dir)\n",
    "    \n",
    "    return(sample_reads_dict)\n",
    "\n",
    "exp_db_dict = {'Bs': {'filter_db': '/home/jgsanders/ref_data/genomes/Homo_sapiens_Bowtie2_v0.1/Homo_sapiens',\n",
    "                      'forward': ['./example/reads/Bs_S106_L001_R1_001.fastq.gz',\n",
    "                                  './example/reads/Bs_S106_L002_R1_001.fastq.gz'],\n",
    "                      'reverse': ['./example/reads/Bs_S106_L001_R2_001.fastq.gz',\n",
    "                                  './example/reads/Bs_S106_L002_R2_001.fastq.gz']},\n",
    "               'Vf': {'filter_db': '/home/jgsanders/ref_data/genomes/Homo_sapiens_Bowtie2_v0.1/Homo_sapiens',\n",
    "                      'forward': ['./example/reads/Vf_S104_L001_R1_001.fastq.gz',\n",
    "                                  './example/reads/Vf_S104_L002_R1_001.fastq.gz'],\n",
    "                      'reverse': ['./example/reads/Vf_S104_L001_R2_001.fastq.gz',\n",
    "                                  './example/reads/Vf_S104_L002_R2_001.fastq.gz']}}\n",
    "\n",
    "exp_db_dict_update = {'Bs': {'filter_db': '/home/jgsanders/ref_data/genomes/Homo_sapiens_Bowtie2_v0.1/Homo_sapiens',\n",
    "                      'forward': ['./example/reads/Bs_S106_L001_R1_001.fastq.gz',\n",
    "                                  './example/reads/Bs_S106_L002_R1_001.fastq.gz'],\n",
    "                      'reverse': ['./example/reads/Bs_S106_L001_R2_001.fastq.gz',\n",
    "                                  './example/reads/Bs_S106_L002_R2_001.fastq.gz']},\n",
    "               'Vf': {'filter_db': '/home/jgsanders/ref_data/genomes/mouse/mouse',\n",
    "                      'forward': ['./example/reads/Vf_S104_L001_R1_001.fastq.gz',\n",
    "                                  './example/reads/Vf_S104_L002_R1_001.fastq.gz'],\n",
    "                      'reverse': ['./example/reads/Vf_S104_L001_R2_001.fastq.gz',\n",
    "                                  './example/reads/Vf_S104_L002_R2_001.fastq.gz']}}\n",
    "\n",
    "exp_db_dict_partial = {'Bs': {'filter_db': None,\n",
    "                      'forward': ['./example/reads/Bs_S106_L001_R1_001.fastq.gz',\n",
    "                                  './example/reads/Bs_S106_L002_R1_001.fastq.gz'],\n",
    "                      'reverse': ['./example/reads/Bs_S106_L001_R2_001.fastq.gz',\n",
    "                                  './example/reads/Bs_S106_L002_R2_001.fastq.gz']},\n",
    "               'Vf': {'filter_db': '/home/jgsanders/ref_data/genomes/mouse/mouse',\n",
    "                      'forward': ['./example/reads/Vf_S104_L001_R1_001.fastq.gz',\n",
    "                                  './example/reads/Vf_S104_L002_R1_001.fastq.gz'],\n",
    "                      'reverse': ['./example/reads/Vf_S104_L001_R2_001.fastq.gz',\n",
    "                                  './example/reads/Vf_S104_L002_R2_001.fastq.gz']}}\n",
    "\n",
    "def add_filter_db(sample_fp_dict, db_fp, samples = None):\n",
    "    \n",
    "    if samples is None:\n",
    "        samples = sample_fp_dict.keys()\n",
    "    \n",
    "    samples_dict = copy.deepcopy(sample_fp_dict)\n",
    "    \n",
    "    for s in samples_dict:\n",
    "        if s in samples:\n",
    "            samples_dict[s]['filter_db'] = db_fp\n",
    "        elif 'filter_db' in samples_dict[s]:\n",
    "            continue\n",
    "        else:\n",
    "            samples_dict[s]['filter_db'] = None\n",
    "        \n",
    "\n",
    "        \n",
    "    return(samples_dict)\n",
    "\n",
    "obs_dict = add_filter_db(exp_dict, '/home/jgsanders/ref_data/genomes/Homo_sapiens_Bowtie2_v0.1/Homo_sapiens')\n",
    "\n",
    "assert(obs_dict == exp_db_dict)\n",
    "\n",
    "obs_dict = add_filter_db(exp_db_dict, '/home/jgsanders/ref_data/genomes/mouse/mouse', samples = ['Vf'])\n",
    "\n",
    "assert(obs_dict == exp_db_dict_update)\n",
    "\n",
    "obs_dict = add_filter_db(exp_dict, '/home/jgsanders/ref_data/genomes/mouse/mouse', samples = ['Vf'])\n",
    "\n",
    "assert(obs_dict == exp_db_dict_partial)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Steps: \n",
    "\n",
    "1. Get samples\n",
    "    - read in sample sheet\n",
    "    - guess samples from sequence folder\n",
    "2. Get params\n",
    "    - filter db per sample\n",
    "    - samples for binning\n",
    "    - samples for abundance profiling\n",
    "    - tmp dir\n",
    "    - envs\n",
    "    - software\n",
    "    - assembler\n",
    "    - trimmer\n",
    "    - params\n",
    "        * atropos\n",
    "        * maxbin \n",
    "        * humann2\n",
    "        * metaphlan\n",
    "        * kraken\n",
    "        * shogun?\n",
    "3. format yaml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: get samples\n",
    "\n",
    "First, we need to read in a set of samples for analysis. \n",
    "\n",
    "There are two options for this: guess the sample names, or read them in from a sample sheet or manifest. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Option 1: guess sample names from sequence output folder. \n",
    "\n",
    "Enter the correct sequencing directory below, and a list of some example filenames should pop up. Make sure these look correct. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 180864\r\n",
      "-rw-r--r--+ 1 jonsanders  staff  14312564 Aug  7 05:23 Bs_S106_L001_R1_001.fastq.gz\r\n",
      "-rw-r--r--+ 1 jonsanders  staff   9523184 Aug  7 05:23 Bs_S106_L001_R2_001.fastq.gz\r\n",
      "-rw-r--r--+ 1 jonsanders  staff  10451788 Aug  7 05:23 Bs_S106_L002_R1_001.fastq.gz\r\n",
      "-rw-r--r--+ 1 jonsanders  staff   9305476 Aug  7 05:23 Bs_S106_L002_R2_001.fastq.gz\r\n",
      "-rw-r--r--+ 1 jonsanders  staff  13733195 Aug  7 05:23 Vf_S104_L001_R1_001.fastq.gz\r\n",
      "-rw-r--r--+ 1 jonsanders  staff   8845407 Aug  7 05:23 Vf_S104_L001_R2_001.fastq.gz\r\n",
      "-rw-r--r--+ 1 jonsanders  staff  16096558 Aug  7 05:23 Vf_S104_L002_R1_001.fastq.gz\r\n",
      "-rw-r--r--+ 1 jonsanders  staff  10323953 Aug  7 05:23 Vf_S104_L002_R2_001.fastq.gz\r\n"
     ]
    }
   ],
   "source": [
    "seq_dir = './example/reads'\n",
    "!ls -l {seq_dir} | head"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If that looks correct, execute the following cell. It should result in a dictionary associating samples to filepaths.\n",
    "\n",
    "To validate, a portion of the dictionary is converted to yaml and printed below: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bs:\n",
      "  forward:\n",
      "  - ./example/reads/Bs_S106_L001_R1_001.fastq.gz\n",
      "  - ./example/reads/Bs_S106_L002_R1_001.fastq.gz\n",
      "  reverse:\n",
      "  - ./example/reads/Bs_S106_L001_R2_001.fastq.gz\n",
      "  - ./example/reads/Bs_S106_L002_R2_001.fastq.gz\n",
      "Vf:\n",
      "  forward:\n",
      "  - ./example/reads/Vf_S104_L001_R1_001.fastq.gz\n",
      "  - ./example/reads/Vf_S104_L002_R1_001.fastq.gz\n",
      "  reverse:\n",
      "  - ./example/reads/Vf_S104_L001_R2_001.fastq.gz\n",
      "  - ./example/reads/Vf_S104_L002_R2_001.fastq.gz\n",
      "\n"
     ]
    }
   ],
   "source": [
    "sample_paths = get_sample_paths(seq_dir)\n",
    "\n",
    "print(yaml.dump(sample_paths, default_flow_style = False))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Option 2: read sample names and prefixes from sample sheet\n",
    "\n",
    "This option is in development "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Lane</th>\n",
       "      <th>Sample_ID</th>\n",
       "      <th>Sample_Name</th>\n",
       "      <th>Sample_Plate</th>\n",
       "      <th>Sample_Well</th>\n",
       "      <th>I7_Index_ID</th>\n",
       "      <th>index</th>\n",
       "      <th>I5_Index_ID</th>\n",
       "      <th>index2</th>\n",
       "      <th>Sample_Project</th>\n",
       "      <th>Description</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>Bs</td>\n",
       "      <td>Bs</td>\n",
       "      <td>Example Plate 1</td>\n",
       "      <td>A1</td>\n",
       "      <td>iTru7_101_01</td>\n",
       "      <td>ACGTTACC</td>\n",
       "      <td>iTru5_01_A</td>\n",
       "      <td>ACCGACAA</td>\n",
       "      <td>Example Project</td>\n",
       "      <td>Bs.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>Vf</td>\n",
       "      <td>Vf</td>\n",
       "      <td>Example Plate 1</td>\n",
       "      <td>B1</td>\n",
       "      <td>iTru7_101_02</td>\n",
       "      <td>CTGTGTTG</td>\n",
       "      <td>iTru5_01_B</td>\n",
       "      <td>AGTGGCAA</td>\n",
       "      <td>Example Project</td>\n",
       "      <td>Vf.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>Bs</td>\n",
       "      <td>Bs</td>\n",
       "      <td>Example Plate 1</td>\n",
       "      <td>A1</td>\n",
       "      <td>iTru7_101_01</td>\n",
       "      <td>ACGTTACC</td>\n",
       "      <td>iTru5_01_A</td>\n",
       "      <td>ACCGACAA</td>\n",
       "      <td>Example Project</td>\n",
       "      <td>Bs.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2</td>\n",
       "      <td>Vf</td>\n",
       "      <td>Vf</td>\n",
       "      <td>Example Plate 1</td>\n",
       "      <td>B1</td>\n",
       "      <td>iTru7_101_02</td>\n",
       "      <td>CTGTGTTG</td>\n",
       "      <td>iTru5_01_B</td>\n",
       "      <td>AGTGGCAA</td>\n",
       "      <td>Example Project</td>\n",
       "      <td>Vf.1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Lane Sample_ID Sample_Name     Sample_Plate Sample_Well   I7_Index_ID  \\\n",
       "0     1        Bs          Bs  Example Plate 1          A1  iTru7_101_01   \n",
       "1     1        Vf          Vf  Example Plate 1          B1  iTru7_101_02   \n",
       "2     2        Bs          Bs  Example Plate 1          A1  iTru7_101_01   \n",
       "3     2        Vf          Vf  Example Plate 1          B1  iTru7_101_02   \n",
       "\n",
       "      index I5_Index_ID    index2   Sample_Project Description  \n",
       "0  ACGTTACC  iTru5_01_A  ACCGACAA  Example Project        Bs.1  \n",
       "1  CTGTGTTG  iTru5_01_B  AGTGGCAA  Example Project        Vf.1  \n",
       "2  ACGTTACC  iTru5_01_A  ACCGACAA  Example Project        Bs.1  \n",
       "3  CTGTGTTG  iTru5_01_B  AGTGGCAA  Example Project        Vf.1  "
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# read in sample sheet\n",
    "\n",
    "def read_sample_sheet(f, sep='\\t'):\n",
    "    data = False\n",
    "    data_lines = ''\n",
    "    for line in f:\n",
    "        if data:\n",
    "            if line.startswith('[') or line.strip() == '':\n",
    "                data = False\n",
    "                continue\n",
    "            data_lines += line\n",
    "        elif line.startswith('[Data]'):\n",
    "            data = True\n",
    "            continue\n",
    "        else:\n",
    "            continue\n",
    "            \n",
    "    data_df = pd.read_csv(StringIO(data_lines), sep = '\\t', comment='#')\n",
    "    return(data_df)\n",
    "\n",
    "with open('./example/example_sample_sheet.txt', 'r') as f:\n",
    "    ss_df = read_sample_sheet(f)\n",
    "\n",
    "# column for read prefixes\n",
    "# column for sample names\n",
    "# aggregate into dict\n",
    "\n",
    "ss_df.head()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>File</th>\n",
       "      <th>Sample</th>\n",
       "      <th>Index</th>\n",
       "      <th>Lane</th>\n",
       "      <th>Read</th>\n",
       "      <th>Run</th>\n",
       "      <th>Extension</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Bs_S106_L001_R1_001.fastq.gz</td>\n",
       "      <td>Bs</td>\n",
       "      <td>S106</td>\n",
       "      <td>L001</td>\n",
       "      <td>R1</td>\n",
       "      <td>001</td>\n",
       "      <td>fastq.gz</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Bs_S106_L001_R2_001.fastq.gz</td>\n",
       "      <td>Bs</td>\n",
       "      <td>S106</td>\n",
       "      <td>L001</td>\n",
       "      <td>R2</td>\n",
       "      <td>001</td>\n",
       "      <td>fastq.gz</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Bs_S106_L002_R1_001.fastq.gz</td>\n",
       "      <td>Bs</td>\n",
       "      <td>S106</td>\n",
       "      <td>L002</td>\n",
       "      <td>R1</td>\n",
       "      <td>001</td>\n",
       "      <td>fastq.gz</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Bs_S106_L002_R2_001.fastq.gz</td>\n",
       "      <td>Bs</td>\n",
       "      <td>S106</td>\n",
       "      <td>L002</td>\n",
       "      <td>R2</td>\n",
       "      <td>001</td>\n",
       "      <td>fastq.gz</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Vf_S104_L001_R1_001.fastq.gz</td>\n",
       "      <td>Vf</td>\n",
       "      <td>S104</td>\n",
       "      <td>L001</td>\n",
       "      <td>R1</td>\n",
       "      <td>001</td>\n",
       "      <td>fastq.gz</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Vf_S104_L001_R2_001.fastq.gz</td>\n",
       "      <td>Vf</td>\n",
       "      <td>S104</td>\n",
       "      <td>L001</td>\n",
       "      <td>R2</td>\n",
       "      <td>001</td>\n",
       "      <td>fastq.gz</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Vf_S104_L002_R1_001.fastq.gz</td>\n",
       "      <td>Vf</td>\n",
       "      <td>S104</td>\n",
       "      <td>L002</td>\n",
       "      <td>R1</td>\n",
       "      <td>001</td>\n",
       "      <td>fastq.gz</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Vf_S104_L002_R2_001.fastq.gz</td>\n",
       "      <td>Vf</td>\n",
       "      <td>S104</td>\n",
       "      <td>L002</td>\n",
       "      <td>R2</td>\n",
       "      <td>001</td>\n",
       "      <td>fastq.gz</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           File Sample Index  Lane Read  Run Extension\n",
       "0  Bs_S106_L001_R1_001.fastq.gz     Bs  S106  L001   R1  001  fastq.gz\n",
       "1  Bs_S106_L001_R2_001.fastq.gz     Bs  S106  L001   R2  001  fastq.gz\n",
       "2  Bs_S106_L002_R1_001.fastq.gz     Bs  S106  L002   R1  001  fastq.gz\n",
       "3  Bs_S106_L002_R2_001.fastq.gz     Bs  S106  L002   R2  001  fastq.gz\n",
       "4  Vf_S104_L001_R1_001.fastq.gz     Vf  S104  L001   R1  001  fastq.gz\n",
       "5  Vf_S104_L001_R2_001.fastq.gz     Vf  S104  L001   R2  001  fastq.gz\n",
       "6  Vf_S104_L002_R1_001.fastq.gz     Vf  S104  L002   R1  001  fastq.gz\n",
       "7  Vf_S104_L002_R2_001.fastq.gz     Vf  S104  L002   R2  001  fastq.gz"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seq_dir = './example/reads/'\n",
    "\n",
    "fps = os.listdir(seq_dir)\n",
    "\n",
    "files_df = parse_ilm_fps_to_df(fps)\n",
    "\n",
    "files_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bs.1:\n",
      "  forward:\n",
      "  - ./example/reads/Bs_S106_L001_R1_001.fastq.gz\n",
      "  - ./example/reads/Bs_S106_L002_R1_001.fastq.gz\n",
      "  reverse:\n",
      "  - ./example/reads/Bs_S106_L001_R2_001.fastq.gz\n",
      "  - ./example/reads/Bs_S106_L002_R2_001.fastq.gz\n",
      "Vf.1:\n",
      "  forward:\n",
      "  - ./example/reads/Vf_S104_L001_R1_001.fastq.gz\n",
      "  - ./example/reads/Vf_S104_L002_R1_001.fastq.gz\n",
      "  reverse:\n",
      "  - ./example/reads/Vf_S104_L001_R2_001.fastq.gz\n",
      "  - ./example/reads/Vf_S104_L002_R2_001.fastq.gz\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def get_sample_reads_df_from_sample_sheet(ss_df, seq_dir, sample_col='Description', prefix_col='Sample_ID'):\n",
    "    sample_reads_dict = {}\n",
    "    \n",
    "    samples = list(set(ss_df[sample_col]))\n",
    "    \n",
    "    fps = os.listdir(seq_dir)\n",
    "\n",
    "    files_df = parse_ilm_fps_to_df(fps)\n",
    "    \n",
    "    for s in samples:\n",
    "        # get the subset of sample sheet rows for that sample\n",
    "        sample_rows = ss_df.loc[ss_df[sample_col] == s,]\n",
    "        \n",
    "        f_fps = []\n",
    "        r_fps = []\n",
    "        \n",
    "        # iter across subset table and find file corresponding to each row\n",
    "        for idx, row in sample_rows.iterrows():\n",
    "            lane = 'L{0:03d}'.format(row['Lane']) \n",
    "            \n",
    "            f_fps.extend(files_df.loc[(files_df['Sample'] == row['Sample_ID']) & \n",
    "                                      (files_df['Lane'] == lane) &\n",
    "                                      (files_df['Read'] == 'R1'), 'File'].values)\n",
    "            \n",
    "            r_fps.extend(files_df.loc[(files_df['Sample'] == row['Sample_ID']) & \n",
    "                                      (files_df['Lane'] == lane) &\n",
    "                                      (files_df['Read'] == 'R2'), 'File'].values)\n",
    "        \n",
    "        f_fps = sorted(f_fps)\n",
    "        r_fps = sorted(r_fps)\n",
    "        \n",
    "        sample_reads_dict[s] = {'forward': [os.path.join(seq_dir, x) for x in f_fps],\n",
    "                                'reverse': [os.path.join(seq_dir, x) for x in r_fps]}\n",
    "    \n",
    "    return(sample_reads_dict)\n",
    "\n",
    "sample_reads_dict = get_sample_reads_df_from_sample_sheet(ss_df, seq_dir,  sample_col='Description', prefix_col='Sample_ID')\n",
    "print(yaml.dump(sample_reads_dict, default_flow_style = False))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: define parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    - filter db per sample\n",
    "    - samples for binning\n",
    "    - samples for abundance profiling\n",
    "    - tmp dir\n",
    "    - envs\n",
    "    - software\n",
    "    - assembler\n",
    "    - trimmer\n",
    "    - params\n",
    "        * atropos\n",
    "        * maxbin \n",
    "        * humann2\n",
    "        * metaphlan\n",
    "        * kraken\n",
    "        * shogun?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "config_dict = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add filter databases per sample\n",
    "\n",
    "These are the filepaths to the databases used for host filtering. \n",
    "\n",
    "The method 'add host filter db' adds the provided filepath to the sample dictionary you created above.\n",
    "\n",
    "By default, all samples get same database.\n",
    "\n",
    "You can also provide a list of specific sample names to update if you want to set different filter dbs for some samples. To do this, just re-execute `add_filter_db` providing a list of samples to update; the existing `filter_db` paths will remain unchanged if they are not in this list. \n",
    "\n",
    "Example:\n",
    "\n",
    "```\n",
    "filter_db1 = '/home/jgsanders/ref_data/genomes/Homo_sapiens_Bowtie2_v0.1/Homo_sapiens'\n",
    "filter_db2 = '/home/jgsanders/ref_data/genomes/mouse/mouse'\n",
    "\n",
    "samples_dict = add_filter_db(sample_paths, filter_db1)\n",
    "samples_dict = add_filter_db(samples_dict, filter_db2, samples = ['Vf'])\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bs:\n",
      "  filter_db: /home/jgsanders/ref_data/genomes/Homo_sapiens_Bowtie2_v0.1/Homo_sapiens\n",
      "  forward:\n",
      "  - ./example/reads/Bs_S106_L001_R1_001.fastq.gz\n",
      "  - ./example/reads/Bs_S106_L002_R1_001.fastq.gz\n",
      "  reverse:\n",
      "  - ./example/reads/Bs_S106_L001_R2_001.fastq.gz\n",
      "  - ./example/reads/Bs_S106_L002_R2_001.fastq.gz\n",
      "Vf:\n",
      "  filter_db: /home/jgsanders/ref_data/genomes/Homo_sapiens_Bowtie2_v0.1/Homo_sapiens\n",
      "  forward:\n",
      "  - ./example/reads/Vf_S104_L001_R1_001.fastq.gz\n",
      "  - ./example/reads/Vf_S104_L002_R1_001.fastq.gz\n",
      "  reverse:\n",
      "  - ./example/reads/Vf_S104_L001_R2_001.fastq.gz\n",
      "  - ./example/reads/Vf_S104_L002_R2_001.fastq.gz\n",
      "\n"
     ]
    }
   ],
   "source": [
    "filter_db = '/home/jgsanders/ref_data/genomes/Homo_sapiens_Bowtie2_v0.1/Homo_sapiens'\n",
    "\n",
    "samples_dict = add_filter_db(sample_paths, filter_db)\n",
    "\n",
    "print(yaml.dump(samples_dict, default_flow_style = False))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add samples for binning\n",
    "\n",
    "These are the samples from which assembled contigs will be binned into putative draft genomes.\n",
    "\n",
    "By default, bin all samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# binning_samples = ['Vf','Bs']\n",
    "binning_samples = list(samples_dict.keys())\n",
    "\n",
    "config_dict['binning_samples'] = binning_samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add samples for abundance profiling\n",
    "\n",
    "These are the samples from which abundance information will be used for binning of `binning_samples`\n",
    "\n",
    "By default, use all samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#abundance_samples = ['Vf','Bs']\n",
    "\n",
    "# This line filters out BLANK* samples from abundance mapping\n",
    "abundance_samples = [x for x in samples_dict if not x.startswith('BLANK')]\n",
    "\n",
    "# This line uses all samples for abundance mapping\n",
    "abundance_samples = list(samples_dict.keys())\n",
    "\n",
    "config_dict['abundance_samples'] = abundance_samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add trimmer to use\n",
    "\n",
    "These are the read trimmers to use for qc.\n",
    "\n",
    "Currently, the available options are:\n",
    "\n",
    "```\n",
    "- atropos\n",
    "- skewer\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "config_dict['trimmer'] = 'atropos'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add assemblers to use\n",
    "\n",
    "These are the assemblers to use for assembly.\n",
    "\n",
    "Currently, the available options are:\n",
    "\n",
    "```\n",
    "- spades\n",
    "- metaspades\n",
    "- megahit\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "assemblers = ['megahit','metaspades']\n",
    "\n",
    "config_dict['assemblers'] = assemblers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add assembly to use for binning\n",
    "\n",
    "This is the assembler to use for binning.\n",
    "\n",
    "Currently, the available options are:\n",
    "\n",
    "```\n",
    "- spades\n",
    "- metaspades\n",
    "- megahit\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "config_dict['mapping_assembler'] = 'metaspades'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add temporary directory\n",
    "\n",
    "This is the directory used during execution of rules. Should ideally be local scratch storage. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tmp_dir_root = '/localscratch'\n",
    "\n",
    "config_dict['tmp_dir_root'] = tmp_dir_root"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add environments\n",
    "\n",
    "These are the conda environments to use for various rules. \n",
    "\n",
    "`env` is the primary `snakemake_assemble` environment. Change others as necessary for other rules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "config_dict['env'] = 'source activate snakemake_assemble'\n",
    "config_dict['anvi_env'] = 'source activate sn_anvio'\n",
    "config_dict['metaphlan_env'] = 'source activate humann2'\n",
    "config_dict['humann2_env'] = 'source activate humann2'\n",
    "config_dict['kraken_env'] = 'source activate kraken'\n",
    "config_dict['shogun_env'] = 'source activate shogun'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add software paths\n",
    "\n",
    "Some software binaries are specified directly. \n",
    "\n",
    "Also, the installation folder for the snakemake_assemble repository is provided as `snakemake_folder`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "software = {'snakemake_folder': '/home/jgsanders/git_sw/snakemake_assemble',\n",
    "            'seqtk': 'seqtk'}\n",
    "\n",
    "config_dict['software'] = software"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### anvio resources\n",
    "\n",
    "Sets paths to resources for anvio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "resources =  {'centrifuge_base': '/home/jgsanders/miniconda/envs/anvio2/centrifuge',\n",
    "              'centrifuge_models': '/home/jgsanders/miniconda/envs/anvio2/centrifuge/b+h+v/b+h+v'}\n",
    "config_dict['resources'] = resources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add parameters definitions\n",
    "\n",
    "These are parameters passed to the various tools used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "params = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Atropos params\n",
    "\n",
    "Sets quality and adapter trimming parameters. \n",
    "\n",
    "Some suggested defaults: \n",
    "Nextera: `-a CTGTCTCTTATACACATCT -A CTGTCTCTTATACACATCT -q 15 --minimum-length 100 --pair-filter any`\n",
    "KapaHyperPlus: `-a GATCGGAAGAGCACACGTCTGAACTCCAGTCAC -A GATCGGAAGAGCGTCGTGTAGGGAAAGGAGTGT -q 15 --minimum-length 100 --pair-filter any`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "params['atropos'] = ' -a GATCGGAAGAGCACACGTCTGAACTCCAGTCAC -A GATCGGAAGAGCGTCGTGTAGGGAAAGGAGTGT -q 15 --minimum-length 100 --pair-filter any'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### maxbin params\n",
    "\n",
    "Sets binning parameters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "params['maxbin'] = '-plotmarker'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Kraken params\n",
    "\n",
    "Sets parameters for Kraken taxonomic profiling. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "params['kraken'] = {'kraken_db': '/home/qiz173/Databases/Kraken/stdb'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### SHOGUN params\n",
    "\n",
    "Sets parameters for Shogun taxonomic profiling. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "params['shogun'] = \"--utree_indx /home/qiz173/Databases/SHOGUN/annotated/utree/stdb.ctr\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### MetaPhlAn2 params\n",
    "\n",
    "Sets parameters for MetaPhlAn2 taxonomic profilling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "params['metaphlan'] = {'metaphlan_dir': '/home/jgsanders/git_sw/metaphlan2'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### humann2\n",
    "\n",
    "Sets various parameters for HUMAnN2 funcitonal profiling. \n",
    "\n",
    "`humann2_aa_db`: Amino acid database for translated amino acid search\n",
    "`humann2_nt_db`: ChocoPhlAn database for nucleotide search\n",
    "`metaphan_dir`: path to metaphlan2 directory install, should also have the default metaphlan2 db\n",
    "`norms`: normalizations for which to output tables\n",
    "`other`: any other HUMAnN2 parameters to use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "humann2 = {'humann2_aa_db': 'humann2_aa_db: /home/jgsanders/ref_data/humann2/uniref',\n",
    "           'humann2_nt_db': '/home/jgsanders/ref_data/humann2/chocophlan',\n",
    "           'metaphlan_dir': '/home/jgsanders/share/metaphlan2',\n",
    "           'norms': ['cpm','relab'],\n",
    "           'other': ''}\n",
    "\n",
    "params['humann2'] = humann2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### other\n",
    "\n",
    "Various other legacy parameters that I need to fix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "params['cram'] = ''\n",
    "params['skewer'] = '-x ./adapters/HyperPlus.fa -n -l 100 -m any -q 15'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: format yaml\n",
    "\n",
    "Now the entire configuration dictionary can be formatted as a yaml file and saved for use in the pipeline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "config_dict['params'] = params\n",
    "config_dict['samples'] = samples_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Format as yaml and write to filepath:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "config_yaml = yaml.dump(config_dict, default_flow_style = False)\n",
    "\n",
    "with open('config.yaml', 'w') as f:\n",
    "    f.write(config_yaml)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the complete config yaml:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "abundance_samples:\n",
      "- Bs\n",
      "- Vf\n",
      "anvi_env: source activate sn_anvio\n",
      "assemblers:\n",
      "- megahit\n",
      "- metaspades\n",
      "binning_samples:\n",
      "- Bs\n",
      "- Vf\n",
      "env: source activate snakemake_assemble\n",
      "humann2_env: source activate humann2\n",
      "kraken_env: source activate kraken\n",
      "mapping_assembler: metaspades\n",
      "metaphlan_env: source activate humann2\n",
      "params:\n",
      "  atropos: ' -a GATCGGAAGAGCACACGTCTGAACTCCAGTCAC -A GATCGGAAGAGCGTCGTGTAGGGAAAGGAGTGT\n",
      "    -q 15 --minimum-length 100 --pair-filter any'\n",
      "  cram: ''\n",
      "  humann2:\n",
      "    humann2_aa_db: 'humann2_aa_db: /home/jgsanders/ref_data/humann2/uniref'\n",
      "    humann2_nt_db: /home/jgsanders/ref_data/humann2/chocophlan\n",
      "    metaphlan_dir: /home/jgsanders/share/metaphlan2\n",
      "    norms:\n",
      "    - cpm\n",
      "    - relab\n",
      "    other: ''\n",
      "  kraken:\n",
      "    kraken_db: /home/qiz173/Databases/Kraken/stdb\n",
      "  maxbin: -plotmarker\n",
      "  metaphlan:\n",
      "    metaphlan_dir: /home/jgsanders/git_sw/metaphlan2\n",
      "  shogun: --utree_indx /home/qiz173/Databases/SHOGUN/annotated/utree/stdb.ctr\n",
      "  skewer: -x ./adapters/HyperPlus.fa -n -l 100 -m any -q 15\n",
      "resources:\n",
      "  centrifuge_base: /home/jgsanders/miniconda/envs/anvio2/centrifuge\n",
      "  centrifuge_models: /home/jgsanders/miniconda/envs/anvio2/centrifuge/b+h+v/b+h+v\n",
      "samples:\n",
      "  Bs:\n",
      "    filter_db: /home/jgsanders/ref_data/genomes/Homo_sapiens_Bowtie2_v0.1/Homo_sapiens\n",
      "    forward:\n",
      "    - ./example/reads/Bs_S106_L001_R1_001.fastq.gz\n",
      "    - ./example/reads/Bs_S106_L002_R1_001.fastq.gz\n",
      "    reverse:\n",
      "    - ./example/reads/Bs_S106_L001_R2_001.fastq.gz\n",
      "    - ./example/reads/Bs_S106_L002_R2_001.fastq.gz\n",
      "  Vf:\n",
      "    filter_db: /home/jgsanders/ref_data/genomes/Homo_sapiens_Bowtie2_v0.1/Homo_sapiens\n",
      "    forward:\n",
      "    - ./example/reads/Vf_S104_L001_R1_001.fastq.gz\n",
      "    - ./example/reads/Vf_S104_L002_R1_001.fastq.gz\n",
      "    reverse:\n",
      "    - ./example/reads/Vf_S104_L001_R2_001.fastq.gz\n",
      "    - ./example/reads/Vf_S104_L002_R2_001.fastq.gz\n",
      "shogun_env: source activate shogun\n",
      "software:\n",
      "  seqtk: seqtk\n",
      "  snakemake_folder: /home/jgsanders/git_sw/snakemake_assemble\n",
      "tmp_dir_root: /localscratch\n",
      "trimmer: atropos\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(config_yaml)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}