run_EIGENSTRAT.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EIGENSTRAT / SMARTPCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Goal: To examine patterns of population structure in duplicated loci."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bring:\n",
    "1. *.haplotypes.tsv file (produced by [Stacks](http://creskolab.uoregon.edu/stacks/))\n",
    "    * used for dominant coding of each allele\n",
    "\n",
    "2. genotype file in *.ped format\n",
    "    * normal codominant genotypes of biallelic SNPs, used to compare \n",
    "    \n",
    "\n",
    "### Take away:\n",
    "* PCAs of genetic data:\n",
    "    - full data coded as codominant genotypes (normal bi-allelic SNPs)\n",
    "    - full data coded as dominant alleles\n",
    "    - paralogs coded as dominant alleles\n",
    "\n",
    "### Programs used:\n",
    "* [EIGENSOFT](http://www.hsph.harvard.edu/alkes-price/software/) (specifically SMARTPCA)\n",
    " - Patterson N, Price AL, Reich D (2006) Population Structure and Eigenanalysis. PLoS Genet 2(12): e190. doi: 10.1371/journal.pgen.0020190\n",
    "    \n",
    "### Steps\n",
    "1. a\n",
    "2. b\n",
    "3. c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Python imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os.path\n",
    "import collections\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import random\n",
    "from IPython.core.pylabtools import figsize\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Dominance coding\n",
    "Haplotypes is a class to store the haplotypes assigned to a group of indidividuals at a particular catalog ID.  Each catalog ID is split into one 'locus' per observed haplotyped.  These loci are named as \"catID_haplotype\"\n",
    "\n",
    "Genotypes of these dominant loci are scored as: \n",
    "    - 1 = allele is present\n",
    "    - 0 = allele is absent\n",
    "    - 9 = no call"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class Haplotypes(object):\n",
    "    def __init__(self, catID, ind_names, haplotypes, missing_codes = [\"\", \"-\", \"No Call\", 'Invalid', 'consensus']):\n",
    "        self.catID = catID\n",
    "        self.missing_codes = missing_codes\n",
    "        self.haplotype_of_ind = collections.OrderedDict(zip(ind_names, haplotypes))\n",
    "        self.set_alleles()\n",
    "        self.split_catID_haplotypes()\n",
    "        \n",
    "    def set_alleles(self):    \n",
    "        seen = set()\n",
    "        for ind, haplo in self.haplotype_of_ind.items():\n",
    "            if haplo in self.missing_codes:\n",
    "                haplo = \"-\"\n",
    "            else:\n",
    "                haplo_set = frozenset(haplo.split(\"/\"))\n",
    "                for allele in haplo_set:\n",
    "                    if allele not in self.missing_codes:\n",
    "                        seen.add(allele)\n",
    "        self.alleles = seen\n",
    "        \n",
    "    def split_catID_haplotypes(self):\n",
    "        self.dom_coding_of_allele = dict()\n",
    "        for allele in self.alleles:\n",
    "            dom_coding = list()\n",
    "            for xx in self.haplotype_of_ind.values():\n",
    "                if xx is '-':\n",
    "                    dom_coding.append(9)\n",
    "                elif allele in xx.split('/'):\n",
    "                    dom_coding.append(1)\n",
    "                else:\n",
    "                    dom_coding.append(0)\n",
    "            self.dom_coding_of_allele['{}_{}'.format(self.catID, allele)] = dom_coding      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trying out the Haplotypes class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CMHAMM10_0002 GG/TA\n",
      "CMHAMM10_0005 GG/TA\n",
      "CMHAMM10_0008 GG/TT\n",
      "CMHAMM10_0011 GG/TA\n",
      "CMHAMM10_0012 GG/TA\n",
      "CMHAMM10_0014 GG/TA\n",
      "CMHAMM10_0015 GG/TA\n",
      "CMHAMM10_0016 GG/TA\n",
      "CMHAMM10_0017 GG/TA\n",
      "CMHAMM10_0018 GG/TA\n",
      "CMHAMM10_0022 GG/TA\n",
      "CMHAMM10_0024 -\n",
      "\n",
      "\n",
      "01_TT [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 9]\n",
      "01_GG [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 9]\n",
      "01_TA [1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 9]\n"
     ]
    }
   ],
   "source": [
    "INDNAMES = ['CMHAMM10_0002','CMHAMM10_0005','CMHAMM10_0008','CMHAMM10_0011','CMHAMM10_0012',\n",
    "            'CMHAMM10_0014','CMHAMM10_0015','CMHAMM10_0016', 'CMHAMM10_0017','CMHAMM10_0018','CMHAMM10_0022',\n",
    "            'CMHAMM10_0024']\n",
    "HAPLOS = ['GG/TA','GG/TA','GG/TT','GG/TA','GG/TA','GG/TA','GG/TA','GG/TA','GG/TA','GG/TA','GG/TA','-']\n",
    "\n",
    "AA = Haplotypes(catID='01', ind_names=INDNAMES, haplotypes=HAPLOS)\n",
    "\n",
    "for ind, haplotypes in AA.haplotype_of_ind.items():\n",
    "    print ind, haplotypes\n",
    "    \n",
    "print '\\n'\n",
    "\n",
    "for locus, genotypes in AA.dom_coding_of_allele.items():\n",
    "    print locus, genotypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### creating the input files for SMARTPCA - dominance coding\n",
    "This function will encode genotypes as a set of dominant alleles.\n",
    "\n",
    "- haplotypes_in = haplotypes.tsv (from Stacks)\n",
    "\n",
    "- indiv_in = ped file of the individuals to retain\n",
    "\n",
    "- catID_in = map or snplist file of the catIDs to retain\n",
    "\n",
    "- genotype_out = genotype input for SMARTPCA\n",
    "\n",
    "- genotype_out = map input for SMARTPCA (list of loci/alleles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def write_dom_EIGENSTRAT_files(haplotypes_in, ind_in, catID_in, genotype_out, map_out):\n",
    "    with open(ind_in) as ind_INFILE:\n",
    "        ind_to_keep = [line.split()[1] for line in ind_INFILE]\n",
    "        print('Inds to keep: {}'.format(len(ind_to_keep)))\n",
    "\n",
    "    with open(catID_in) as catID_INFILE:\n",
    "        try: # this works on map files - catId is found by spliting the snpID\n",
    "            catId_to_keep = [line.split()[1].split('_')[0] for line in catID_INFILE]\n",
    "        except IndexError: # this works on snplist files - each line is a catID to keep\n",
    "            catId_to_keep = [line.strip() for line in catID_INFILE]\n",
    "        print('catIDs to keep: {}'.format(len(catId_to_keep)))\n",
    "\n",
    "    with open(haplotypes_in) as INFILE:\n",
    "        HEADER = next(INFILE)\n",
    "        ind_names = HEADER.strip().split('\\t')[2:]\n",
    "        ind_keeps = [xx in ind_to_keep for xx in ind_names]\n",
    "        kept_names = [xx for xx in ind_names if xx in ind_to_keep]\n",
    "        with open(genotype_out, 'w') as geno_OUTFILE:\n",
    "            with open(map_out, 'w') as map_OUTFILE:\n",
    "                for line in INFILE:\n",
    "                    catID = line.split('\\t')[0]\n",
    "                    if catID in catId_to_keep:\n",
    "                        haplos = line.strip().split('\\t')[2:]\n",
    "                        kept_haplos = [xx for (cnt,xx) in enumerate(haplos) if ind_keeps[cnt]]\n",
    "                        haps = Haplotypes(catID, kept_names, kept_haplos)\n",
    "                        for allele, calls in haps.dom_coding_of_allele.items():\n",
    "                            geno_OUTFILE.write(''.join([str(xx) for xx in calls]))\n",
    "                            geno_OUTFILE.write('\\n')\n",
    "                            map_OUTFILE.write('1\\t{}\\t1\\t1\\n'.format(allele))\n",
    "    return(len(catId_to_keep))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice, these are based on **haplotypes** - rather than individual snps the codominant data is based on.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data sets\n",
    "\n",
    "listed as: 'LOCUS_SET' . 'ENCODING' . 'FULL_SUBSAMPLE' . 'EXT'\n",
    "\n",
    "LOCUS_SETS:\n",
    "    1. complete (codominant)\n",
    "    2. complete (dominant alleles)\n",
    "    3. mapped non-paralogs (codominant)\n",
    "    4. mapped non-paralogs (dominant alleles)\n",
    "    5. mapped paralogs (dominant alleles)\n",
    "    \n",
    "    and subsamples of each down to match the number of loci in smallest set (mapped paralogs)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copy data into local dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "base_haplotypes_file = os.path.join('data', 'batch_4', 'all.haplotypes.tsv')\n",
    "base_individual_file = os.path.join('data', 'batch_4', 'pop_genotypes', 'complete.ped')\n",
    "complete_loci_file = os.path.join('data', 'batch_4', 'pop_genotypes', 'complete.map')\n",
    "mapped_loci_file = os.path.join('data', 'batch_4', 'pop_genotypes', 'on_map.map')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import shutil\n",
    "\n",
    "shutil.copy(base_haplotypes_file, os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.haplotypes.tsv'))\n",
    "shutil.copy(base_individual_file, os.path.join('results', 'batch_4', 'EIGENSOFT', 'individuals.ped'))\n",
    "shutil.copy(complete_loci_file,   os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.map'))\n",
    "\n",
    "shutil.copy(mapped_loci_file , os.path.join('results', 'batch_4', 'EIGENSOFT', 'on_map.map'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### edit the complete.map file, setting chrom=0 to chrom 40\n",
    "SMARTPCA removes loci with chrom=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with open(os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.map')) as INFILE:\n",
    "    with open(os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.map.tmp'), 'w') as OUTFILE:\n",
    "        for line in INFILE:\n",
    "            if int(line.strip().split()[0]) == 0:\n",
    "                OUTFILE.write('40\\t{}\\t{}\\t{}\\n'.format(line.strip().split()[1], line.strip().split()[2], line.strip().split()[3]))\n",
    "            else:\n",
    "                OUTFILE.write(line)\n",
    "\n",
    "shutil.copy(os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.map.tmp'), \n",
    "            os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.map'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inds to keep: 174\n",
      "catIDs to keep: 13407\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "13407"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "write_dom_EIGENSTRAT_files(\n",
    "    haplotypes_in = os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.haplotypes.tsv'),\n",
    "    ind_in        = os.path.join('results', 'batch_4', 'EIGENSOFT', 'individuals.ped'),\n",
    "    catID_in      = os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.map'),\n",
    "    genotype_out  = os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.dom.txt'),\n",
    "    map_out       = os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.dom.map')\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set of all mapped loci mapped loci"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inds to keep: 174\n",
      "catIDs to keep: 7259\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "7259"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "write_dom_EIGENSTRAT_files(\n",
    "    haplotypes_in = os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.haplotypes.tsv'),\n",
    "    ind_in        = os.path.join('results', 'batch_4', 'EIGENSOFT', 'individuals.ped'),\n",
    "    catID_in      = os.path.join('results', 'batch_4', 'EIGENSOFT', 'on_map.map'),\n",
    "    genotype_out  = os.path.join('results', 'batch_4', 'EIGENSOFT', 'on_map.dom.txt'),\n",
    "    map_out       = os.path.join('results', 'batch_4', 'EIGENSOFT', 'on_map.dom.map')\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pick out paralogs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>contig</th>\n",
       "      <th>resolved_locus</th>\n",
       "      <th>stacks_CatID</th>\n",
       "      <th>stacks_SNP</th>\n",
       "      <th>LEPname</th>\n",
       "      <th>LEP_LG</th>\n",
       "      <th>cM_OLD</th>\n",
       "      <th>paper1_LG</th>\n",
       "      <th>cM</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>c4311</td>\n",
       "      <td>4311</td>\n",
       "      <td>28282</td>\n",
       "      <td>28282_88</td>\n",
       "      <td>4584</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>c4311</td>\n",
       "      <td>4311</td>\n",
       "      <td>28282</td>\n",
       "      <td>28282_91</td>\n",
       "      <td>4584</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>c56875</td>\n",
       "      <td>56875</td>\n",
       "      <td>39970</td>\n",
       "      <td>39970_17</td>\n",
       "      <td>7377</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   contig resolved_locus  stacks_CatID stacks_SNP  LEPname  LEP_LG  cM_OLD  \\\n",
       "0   c4311           4311         28282   28282_88     4584       1       0   \n",
       "1   c4311           4311         28282   28282_91     4584       1       0   \n",
       "2  c56875          56875         39970   39970_17     7377       1       0   \n",
       "\n",
       "   paper1_LG  cM  \n",
       "0          1   0  \n",
       "1          1   0  \n",
       "2          1   0  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "linkage_map = pd.read_csv(os.path.join('linkage_map', 'LEPmap', 'with_paralogs', 'final', 'PS_chum_map_2015.txt'), sep = '\\t')\n",
    "linkage_map.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "paralogs = pd.read_csv(os.path.join('linkage_map','chum_paralogs.txt'), header = None)\n",
    "paralogs.columns = ['old_CatID']\n",
    "paralogs['contig'] = ['c' + str(xx) for xx in paralogs['old_CatID']]\n",
    "paralogs['paralog'] = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "catID_paralogs = set(pd.merge(linkage_map, paralogs)['stacks_CatID'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with open(os.path.join('linkage_map','mapped_chum_paralogs.catIDs'), 'w') as OUTFILE:\n",
    "    for xx in catID_paralogs:\n",
    "        OUTFILE.write('{}\\n'.format(xx))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Just the mapped paralogs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inds to keep: 174\n",
      "catIDs to keep: 1214\n"
     ]
    }
   ],
   "source": [
    "num_subsample = write_dom_EIGENSTRAT_files(\n",
    "    haplotypes_in = os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.haplotypes.tsv'),\n",
    "    ind_in        = os.path.join('results', 'batch_4', 'EIGENSOFT', 'individuals.ped'),\n",
    "    catID_in      = os.path.join('linkage_map', 'mapped_chum_paralogs.catIDs'),\n",
    "    genotype_out  = os.path.join('results', 'batch_4', 'EIGENSOFT', 'paralogs.dom.txt'),\n",
    "    map_out       = os.path.join('results', 'batch_4', 'EIGENSOFT', 'paralogs.dom.map')\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Subsampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def subsample_map_file(infile, outfile, num_keep = 1000):\n",
    "    with open(infile) as INFILE:\n",
    "        num_loci = sum(1 for line in INFILE)\n",
    "    lines_to_keep = random.sample(range(num_loci), num_keep)\n",
    "    with open(outfile, 'w') as OUTFILE:\n",
    "        with open(infile) as INFILE:\n",
    "            idx = 0\n",
    "            for line in INFILE:\n",
    "                if idx in lines_to_keep:\n",
    "                    OUTFILE.write(line)\n",
    "                idx += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "subsample_map_file(infile  = os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.map'), \n",
    "                   outfile = os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.subsample.map'), \n",
    "                   num_keep = num_subsample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "subsample_map_file(infile  = os.path.join('results', 'batch_4', 'EIGENSOFT', 'on_map.map'), \n",
    "                   outfile = os.path.join('results', 'batch_4', 'EIGENSOFT', 'on_map.subsample.map'), \n",
    "                   num_keep = num_subsample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def subsample_DOMINANT_map_file(infile, outfile, num_keep = 1000):\n",
    "    with open(infile) as INFILE:\n",
    "        loci = set([int(line.split()[1].split('_')[0]) for line in INFILE])\n",
    "    loci_to_keep = random.sample(loci, num_keep)\n",
    "    #print(loci_to_keep)\n",
    "    with open(outfile, 'w') as OUTFILE:\n",
    "        with open(infile) as INFILE:\n",
    "            for line in INFILE:\n",
    "                current_locus = int(line.split()[1].split('_')[0])\n",
    "                if current_locus in loci_to_keep:\n",
    "                    OUTFILE.write(line)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "subsample_DOMINANT_map_file(\n",
    "        infile = os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.dom.map'), \n",
    "        outfile =  os.path.join('results', 'batch_4', 'EIGENSOFT', 'complete.dom.subsample.map'), \n",
    "        num_keep = num_subsample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "subsample_DOMINANT_map_file(\n",
    "        infile = os.path.join('results', 'batch_4', 'EIGENSOFT', 'on_map.dom.map'), \n",
    "        outfile =  os.path.join('results', 'batch_4', 'EIGENSOFT', 'on_map.dom.subsample.map'), \n",
    "        num_keep = num_subsample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def subsample_DOMINANT_allele_file(infile, outfile, old_map, new_map):\n",
    "    with open(old_map) as OLD_MAP:\n",
    "        with open(new_map) as NEW_MAP:\n",
    "            OLD_LINES = OLD_MAP.readlines()\n",
    "            NEW_LINES = NEW_MAP.readlines()\n",
    "            lines_to_keep = ([OLD_LINES.index(xx) for xx in NEW_LINES])\n",
    "    with open(infile) as INFILE:\n",
    "        with open(outfile, 'w') as OUTFILE:\n",
    "            line_index = 0\n",
    "            for line in INFILE:\n",
    "                if line_index in lines_to_keep:\n",
    "                    OUTFILE.write(line)\n",
    "                line_index += 1\n",
    "        \n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "subsample_DOMINANT_allele_file('/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.txt',\n",
    "                               '/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.subsample.txt', \n",
    "                               '/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.map',\n",
    "                              '/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.subsample.map')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "subsample_DOMINANT_allele_file('/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.txt',\n",
    "                               '/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.subsample.txt', \n",
    "                               '/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.map',\n",
    "                              '/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.subsample.map')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### subsample the genotype file using PLINK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "univseral_plink_commands = \"--allow-extra-chr --allow-no-sex --write-snplist --autosome-num 50\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### subsample complete"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PLINK v1.90b3q 64-bit (29 May 2015)        https://www.cog-genomics.org/plink2\n",
      "(C) 2005-2015 Shaun Purcell, Christopher Chang   GNU General Public License v3\n",
      "Logging to results/batch_4/EIGENSOFT/complete.codom.subsample.log.\n",
      "Options in effect:\n",
      "  --allow-extra-chr\n",
      "  --allow-no-sex\n",
      "  --autosome-num 50\n",
      "  --file results/batch_4/EIGENSOFT/complete\n",
      "  --out results/batch_4/EIGENSOFT/complete.codom.subsample\n",
      "  --recode\n",
      "  --thin-count 1214\n",
      "  --write-snplist\n",
      "\n",
      "32127 MB RAM detected; reserving 16063 MB for main workspace.\n",
      ".ped scan complete (for binary autoconversion).\n",
      "Performing single-pass .bed write (13407 variants, 174 samples).\n",
      "--file: results/batch_4/EIGENSOFT/complete.codom.subsample-temporary.bed +\n",
      "results/batch_4/EIGENSOFT/complete.codom.subsample-temporary.bim +\n",
      "results/batch_4/EIGENSOFT/complete.codom.subsample-temporary.fam written.\n",
      "13407 variants loaded from .bim file.\n",
      "174 samples (0 males, 0 females, 174 ambiguous) loaded from .fam.\n",
      "Ambiguous sex IDs written to\n",
      "results/batch_4/EIGENSOFT/complete.codom.subsample.nosex .\n",
      "--thin-count: 12193 variants removed (1214 remaining).\n",
      "Using 1 thread (no multithreaded calculations invoked).\n",
      "Before main variant filters, 174 founders and 0 nonfounders present.\n",
      "Calculating allele frequencies... 0%\b\b1%\b\b2%\b\b3%\b\b4%\b\b5%\b\b6%\b\b7%\b\b8%\b\b9%\b\b10%\b\b\b11%\b\b\b12%\b\b\b13%\b\b\b14%\b\b\b15%\b\b\b16%\b\b\b17%\b\b\b18%\b\b\b19%\b\b\b20%\b\b\b21%\b\b\b22%\b\b\b23%\b\b\b24%\b\b\b25%\b\b\b26%\b\b\b27%\b\b\b28%\b\b\b29%\b\b\b30%\b\b\b31%\b\b\b32%\b\b\b33%\b\b\b34%\b\b\b35%\b\b\b36%\b\b\b37%\b\b\b38%\b\b\b39%\b\b\b40%\b\b\b41%\b\b\b42%\b\b\b43%\b\b\b44%\b\b\b45%\b\b\b46%\b\b\b47%\b\b\b48%\b\b\b49%\b\b\b50%\b\b\b51%\b\b\b52%\b\b\b53%\b\b\b54%\b\b\b55%\b\b\b56%\b\b\b57%\b\b\b58%\b\b\b59%\b\b\b60%\b\b\b61%\b\b\b62%\b\b\b63%\b\b\b64%\b\b\b65%\b\b\b66%\b\b\b67%\b\b\b68%\b\b\b69%\b\b\b70%\b\b\b71%\b\b\b72%\b\b\b73%\b\b\b74%\b\b\b75%\b\b\b76%\b\b\b77%\b\b\b78%\b\b\b79%\b\b\b80%\b\b\b81%\b\b\b82%\b\b\b83%\b\b\b84%\b\b\b85%\b\b\b86%\b\b\b87%\b\b\b88%\b\b\b89%\b\b\b90%\b\b\b91%\b\b\b92%\b\b\b93%\b\b\b94%\b\b\b95%\b\b\b96%\b\b\b97%\b\b\b98%\b\b\b99%\b\b\b\b done.\n",
      "Total genotyping rate is 0.928028.\n",
      "1214 variants and 174 samples pass filters and QC.\n",
      "Note: No phenotypes present.\n",
      "List of variant IDs written to\n",
      "results/batch_4/EIGENSOFT/complete.codom.subsample.snplist .\n",
      "--recode to results/batch_4/EIGENSOFT/complete.codom.subsample.ped +\n",
      "results/batch_4/EIGENSOFT/complete.codom.subsample.map ... 0%\b\b1%\b\b2%\b\b3%\b\b4%\b\b5%\b\b6%\b\b7%\b\b8%\b\b9%\b\b10%\b\b\b11%\b\b\b12%\b\b\b13%\b\b\b14%\b\b\b15%\b\b\b16%\b\b\b17%\b\b\b18%\b\b\b19%\b\b\b20%\b\b\b21%\b\b\b22%\b\b\b23%\b\b\b24%\b\b\b25%\b\b\b26%\b\b\b27%\b\b\b28%\b\b\b29%\b\b\b30%\b\b\b31%\b\b\b32%\b\b\b33%\b\b\b34%\b\b\b35%\b\b\b36%\b\b\b37%\b\b\b38%\b\b\b39%\b\b\b40%\b\b\b41%\b\b\b42%\b\b\b43%\b\b\b44%\b\b\b45%\b\b\b46%\b\b\b47%\b\b\b48%\b\b\b49%\b\b\b50%\b\b\b51%\b\b\b52%\b\b\b53%\b\b\b54%\b\b\b55%\b\b\b56%\b\b\b57%\b\b\b58%\b\b\b59%\b\b\b60%\b\b\b61%\b\b\b62%\b\b\b63%\b\b\b64%\b\b\b65%\b\b\b66%\b\b\b67%\b\b\b68%\b\b\b69%\b\b\b70%\b\b\b71%\b\b\b72%\b\b\b73%\b\b\b74%\b\b\b75%\b\b\b76%\b\b\b77%\b\b\b78%\b\b\b79%\b\b\b80%\b\b\b81%\b\b\b82%\b\b\b83%\b\b\b84%\b\b\b85%\b\b\b86%\b\b\b87%\b\b\b88%\b\b\b89%\b\b\b90%\b\b\b91%\b\b\b92%\b\b\b93%\b\b\b94%\b\b\b95%\b\b\b96%\b\b\b97%\b\b\b98%\b\b\b99%\b\b\bdone.\n"
     ]
    }
   ],
   "source": [
    "!plink --file results/batch_4/EIGENSOFT/complete {univseral_plink_commands} \\\n",
    "    --thin-count {num_subsample} --recode --out results/batch_4/EIGENSOFT/complete.codom.subsample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### subsample on_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PLINK v1.90b3q 64-bit (29 May 2015)        https://www.cog-genomics.org/plink2\n",
      "(C) 2005-2015 Shaun Purcell, Christopher Chang   GNU General Public License v3\n",
      "Logging to results/batch_4/EIGENSOFT/on_map.codom.subsample.log.\n",
      "Options in effect:\n",
      "  --allow-extra-chr\n",
      "  --allow-no-sex\n",
      "  --autosome-num 50\n",
      "  --exclude results/batch_4/EIGENSOFT/paralog_SNPs.txt\n",
      "  --file results/batch_4/EIGENSOFT/complete\n",
      "  --not-chr 40\n",
      "  --out results/batch_4/EIGENSOFT/on_map.codom.subsample\n",
      "  --recode\n",
      "  --thin-count 1214\n",
      "  --write-snplist\n",
      "\n",
      "32127 MB RAM detected; reserving 16063 MB for main workspace.\n",
      ".ped scan complete (for binary autoconversion).\n",
      "Performing single-pass .bed write (7259 variants, 174 samples).\n",
      "--file: results/batch_4/EIGENSOFT/on_map.codom.subsample-temporary.bed +\n",
      "results/batch_4/EIGENSOFT/on_map.codom.subsample-temporary.bim +\n",
      "results/batch_4/EIGENSOFT/on_map.codom.subsample-temporary.fam written.\n",
      "7259 variants loaded from .bim file.\n",
      "174 samples (0 males, 0 females, 174 ambiguous) loaded from .fam.\n",
      "Ambiguous sex IDs written to\n",
      "results/batch_4/EIGENSOFT/on_map.codom.subsample.nosex .\n",
      "--exclude: 6251 variants remaining.\n",
      "Warning: At least 263 duplicate IDs in --exclude file.\n",
      "--thin-count: 5037 variants removed (1214 remaining).\n",
      "Using 1 thread (no multithreaded calculations invoked).\n",
      "Before main variant filters, 174 founders and 0 nonfounders present.\n",
      "Calculating allele frequencies... 0%\b\b1%\b\b2%\b\b3%\b\b4%\b\b5%\b\b6%\b\b7%\b\b8%\b\b9%\b\b10%\b\b\b11%\b\b\b12%\b\b\b13%\b\b\b14%\b\b\b15%\b\b\b16%\b\b\b17%\b\b\b18%\b\b\b19%\b\b\b20%\b\b\b21%\b\b\b22%\b\b\b23%\b\b\b24%\b\b\b25%\b\b\b26%\b\b\b27%\b\b\b28%\b\b\b29%\b\b\b30%\b\b\b31%\b\b\b32%\b\b\b33%\b\b\b34%\b\b\b35%\b\b\b36%\b\b\b37%\b\b\b38%\b\b\b39%\b\b\b40%\b\b\b41%\b\b\b42%\b\b\b43%\b\b\b44%\b\b\b45%\b\b\b46%\b\b\b47%\b\b\b48%\b\b\b49%\b\b\b50%\b\b\b51%\b\b\b52%\b\b\b53%\b\b\b54%\b\b\b55%\b\b\b56%\b\b\b57%\b\b\b58%\b\b\b59%\b\b\b60%\b\b\b61%\b\b\b62%\b\b\b63%\b\b\b64%\b\b\b65%\b\b\b66%\b\b\b67%\b\b\b68%\b\b\b69%\b\b\b70%\b\b\b71%\b\b\b72%\b\b\b73%\b\b\b74%\b\b\b75%\b\b\b76%\b\b\b77%\b\b\b78%\b\b\b79%\b\b\b80%\b\b\b81%\b\b\b82%\b\b\b83%\b\b\b84%\b\b\b85%\b\b\b86%\b\b\b87%\b\b\b88%\b\b\b89%\b\b\b90%\b\b\b91%\b\b\b92%\b\b\b93%\b\b\b94%\b\b\b95%\b\b\b96%\b\b\b97%\b\b\b98%\b\b\b99%\b\b\b\b done.\n",
      "Total genotyping rate is 0.92925.\n",
      "1214 variants and 174 samples pass filters and QC.\n",
      "Note: No phenotypes present.\n",
      "List of variant IDs written to\n",
      "results/batch_4/EIGENSOFT/on_map.codom.subsample.snplist .\n",
      "--recode to results/batch_4/EIGENSOFT/on_map.codom.subsample.ped +\n",
      "results/batch_4/EIGENSOFT/on_map.codom.subsample.map ... 0%\b\b1%\b\b2%\b\b3%\b\b4%\b\b5%\b\b6%\b\b7%\b\b8%\b\b9%\b\b10%\b\b\b11%\b\b\b12%\b\b\b13%\b\b\b14%\b\b\b15%\b\b\b16%\b\b\b17%\b\b\b18%\b\b\b19%\b\b\b20%\b\b\b21%\b\b\b22%\b\b\b23%\b\b\b24%\b\b\b25%\b\b\b26%\b\b\b27%\b\b\b28%\b\b\b29%\b\b\b30%\b\b\b31%\b\b\b32%\b\b\b33%\b\b\b34%\b\b\b35%\b\b\b36%\b\b\b37%\b\b\b38%\b\b\b39%\b\b\b40%\b\b\b41%\b\b\b42%\b\b\b43%\b\b\b44%\b\b\b45%\b\b\b46%\b\b\b47%\b\b\b48%\b\b\b49%\b\b\b50%\b\b\b51%\b\b\b52%\b\b\b53%\b\b\b54%\b\b\b55%\b\b\b56%\b\b\b57%\b\b\b58%\b\b\b59%\b\b\b60%\b\b\b61%\b\b\b62%\b\b\b63%\b\b\b64%\b\b\b65%\b\b\b66%\b\b\b67%\b\b\b68%\b\b\b69%\b\b\b70%\b\b\b71%\b\b\b72%\b\b\b73%\b\b\b74%\b\b\b75%\b\b\b76%\b\b\b77%\b\b\b78%\b\b\b79%\b\b\b80%\b\b\b81%\b\b\b82%\b\b\b83%\b\b\b84%\b\b\b85%\b\b\b86%\b\b\b87%\b\b\b88%\b\b\b89%\b\b\b90%\b\b\b91%\b\b\b92%\b\b\b93%\b\b\b94%\b\b\b95%\b\b\b96%\b\b\b97%\b\b\b98%\b\b\b99%\b\b\bdone.\n"
     ]
    }
   ],
   "source": [
    "!plink --file results/batch_4/EIGENSOFT/complete --not-chr 40 \\\n",
    "    --exclude results/batch_4/EIGENSOFT/paralog_SNPs.txt {univseral_plink_commands} \\\n",
    "    --thin-count {num_subsample} --recode --out results/batch_4/EIGENSOFT/on_map.codom.subsample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### TODO: I want to set up a way to better manage SmartPCA runs \n",
    "I want to be able to:\n",
    "   * downsample\n",
    "   * subset loci by paralog/non"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Steps to run SmartPCA\n",
    "\n",
    "   * create / subsample files\n",
    "\n",
    "   * create input files\n",
    "       * genotype file\n",
    "       * snp file\n",
    "       * individual file\n",
    "   * create parfiles for each analysis\n",
    "   * call SmartPca"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Subsample the map files, retianing XX loci at random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#cd {os.path.join('data','batch_4','pop_genotypes')}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Smartpca\n",
    "Smartpca, within EIGENSOFT, is used for PCA analysis. Smartpca uses a parameter file to specify program options."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Notes\n",
    "\n",
    "   * You may need to check your in input \\*.map files, EIGENSOFT is excluding loci on chr=0.  If you want to retain unplaced loci, set chr (column 1) to 1.\n",
    "\n",
    "   * smartpca also complains about the  6th column of the .ped file (phenotype) when the value is '-9' (ie. missing).  We can set this to 1 or 0 to remove this error, specifing cases or controls. see [FAQ](http://www.hsph.harvard.edu/alkes-price/eigensoft-frequently-asked-questions/)\n",
    "   * Patterson (2006) notes that the normalization ('usenorm') should **not** be applied to dominance-coded data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "SMARTPCA_path = '/home/ipseg/Programs/EIGENSOFT/EIG6.0.1/bin/smartpca'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### write parfile(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def function_to_write_parfiles(base_dir):\n",
    "    ## TODO, this is a placeholder\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "complete_codom_PARFILE_path = '/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.parfile'\n",
    "\n",
    "with open(complete_codom_PARFILE_path, 'w') as OUTFILE:\n",
    "    OUTFILE.write('''\\\n",
    "genotypename:  /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.ped\n",
    "snpname:       /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.map\n",
    "indivname:     /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/individuals.ped\n",
    "\n",
    "evecoutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.evec\n",
    "evaloutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.eval\n",
    "grmoutname:      /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.rel\n",
    "snpweightoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.snpweights\n",
    "deletesnpoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.snpremoved\n",
    "outlieroutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.ind_outliers\n",
    "\n",
    "fastmode:      NO\n",
    "usenorm:       YES\n",
    "missingmode:   NO\n",
    "outliermode:   2 \n",
    "numchrom:      50\n",
    "numoutlieriter: 5 \n",
    "numoutlierevec: 10 \n",
    "lsqproject: NO \\\n",
    "    ''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "complete_dom_PARFILE_path = '/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.parfile'\n",
    "\n",
    "with open(complete_dom_PARFILE_path, 'w') as OUTFILE:\n",
    "    OUTFILE.write('''\\\n",
    "genotypename:  /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.txt\n",
    "snpname:       /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.map\n",
    "indivname:     /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/individuals.ped\n",
    "\n",
    "evecoutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.evec\n",
    "evaloutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.eval\n",
    "grmoutname:      /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.rel\n",
    "snpweightoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.snpweights\n",
    "deletesnpoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.snpremoved\n",
    "outlieroutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.ind_outliers\n",
    "\n",
    "fastmode:      NO\n",
    "usenorm:       NO\n",
    "missingmode:   NO\n",
    "outliermode:   2 \n",
    "numchrom:      50\n",
    "numoutlieriter: 5 \n",
    "numoutlierevec: 10 \n",
    "lsqproject: NO \\\n",
    "    ''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "on_map_codom_PARFILE_path = '/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.parfile'\n",
    "\n",
    "with open(on_map_codom_PARFILE_path, 'w') as OUTFILE:\n",
    "    OUTFILE.write('''\\\n",
    "genotypename:  /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.ped\n",
    "snpname:       /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.map\n",
    "indivname:     /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.ped\n",
    "\n",
    "evecoutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.evec\n",
    "evaloutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.eval\n",
    "grmoutname:      /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.rel\n",
    "snpweightoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.snpweights\n",
    "deletesnpoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.snpremoved\n",
    "outlieroutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.ind_outliers\n",
    "\n",
    "fastmode:      NO\n",
    "usenorm:       YES\n",
    "missingmode:   NO\n",
    "outliermode:   2 \n",
    "numchrom:      50\n",
    "numoutlieriter: 0 \n",
    "numoutlierevec: 0 \n",
    "lsqproject: NO \\\n",
    "    ''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "on_map_dom_PARFILE_path = '/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.parfile'\n",
    "\n",
    "with open(on_map_dom_PARFILE_path, 'w') as OUTFILE:\n",
    "    OUTFILE.write('''\\\n",
    "genotypename:  /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.txt\n",
    "snpname:       /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.map\n",
    "indivname:     /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.ped\n",
    "\n",
    "evecoutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.evec\n",
    "evaloutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.eval\n",
    "grmoutname:      /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.rel\n",
    "snpweightoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.snpweights\n",
    "deletesnpoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.snpremoved\n",
    "outlieroutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.ind_outliers\n",
    "\n",
    "fastmode:      NO\n",
    "usenorm:       NO\n",
    "missingmode:   NO\n",
    "outliermode:   2 \n",
    "numchrom:      50\n",
    "numoutlieriter: 5 \n",
    "numoutlierevec: 10 \n",
    "lsqproject: NO \\\n",
    "    ''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "paralogs_dom_PARFILE_path = '/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/paralogs.dom.parfile'\n",
    "\n",
    "with open(paralogs_dom_PARFILE_path, 'w') as OUTFILE:\n",
    "    OUTFILE.write('''\\\n",
    "genotypename:  /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/paralogs.dom.txt\n",
    "snpname:       /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/paralogs.dom.map\n",
    "indivname:      /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.ped\n",
    "\n",
    "evecoutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/paralogs.dom.evec\n",
    "evaloutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/paralogs.dom.eval\n",
    "grmoutname:      /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/paralogs.dom.rel\n",
    "snpweightoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/paralogs.dom.snpweights\n",
    "deletesnpoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/paralogs.dom.snpremoved\n",
    "outlieroutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/paralogs.dom.ind_outliers\n",
    "\n",
    "fastmode:      NO\n",
    "usenorm:       NO\n",
    "missingmode:   NO\n",
    "outliermode:   2 \n",
    "numchrom:      50\n",
    "numoutlieriter: 5 \n",
    "numoutlierevec: 10 \n",
    "lsqproject: NO \\\n",
    "    ''')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## subsampled data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "complete_codom_subsample_PARFILE_path = '/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.subsample.parfile'\n",
    "\n",
    "with open(complete_codom_subsample_PARFILE_path, 'w') as OUTFILE:\n",
    "    OUTFILE.write('''\\\n",
    "genotypename:  /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.subsample.ped\n",
    "snpname:       /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.subsample.map\n",
    "indivname:     /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.ped\n",
    "\n",
    "evecoutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.subsample.evec\n",
    "evaloutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.subsample.eval\n",
    "grmoutname:      /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.subsample.rel\n",
    "snpweightoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.subsample.snpweights\n",
    "deletesnpoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.subsample.snpremoved\n",
    "outlieroutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.subsample.ind_outliers\n",
    "\n",
    "fastmode:      NO\n",
    "usenorm:       YES\n",
    "missingmode:   NO\n",
    "outliermode:   2 \n",
    "numchrom:      50\n",
    "numoutlieriter: 5 \n",
    "numoutlierevec: 10 \n",
    "lsqproject: NO \\\n",
    "    ''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "complete_dom_subsample_PARFILE_path = '/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.subsample.parfile'\n",
    "\n",
    "with open(complete_dom_subsample_PARFILE_path, 'w') as OUTFILE:\n",
    "    OUTFILE.write('''\\\n",
    "genotypename:  /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.subsample.txt\n",
    "snpname:       /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.subsample.map\n",
    "indivname:     /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.ped\n",
    "\n",
    "evecoutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.subsample.evec\n",
    "evaloutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.subsample.eval\n",
    "grmoutname:      /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.subsample.rel\n",
    "snpweightoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.subsample.snpweights\n",
    "deletesnpoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.subsample.snpremoved\n",
    "outlieroutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.subsample.ind_outliers\n",
    "\n",
    "fastmode:      NO\n",
    "usenorm:       NO\n",
    "missingmode:   NO\n",
    "outliermode:   2 \n",
    "numchrom:      50\n",
    "numoutlieriter: 5 \n",
    "numoutlierevec: 10 \n",
    "lsqproject: NO \\\n",
    "    ''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "on_map_codom_subsample_PARFILE_path = '/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.subsample.parfile'\n",
    "\n",
    "with open(on_map_codom_subsample_PARFILE_path, 'w') as OUTFILE:\n",
    "    OUTFILE.write('''\\\n",
    "genotypename:  /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.subsample.ped\n",
    "snpname:       /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.subsample.map\n",
    "indivname:     /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.ped\n",
    "\n",
    "evecoutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.subsample.evec\n",
    "evaloutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.subsample.eval\n",
    "grmoutname:      /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.subsample.rel\n",
    "snpweightoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.subsample.snpweights\n",
    "deletesnpoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.subsample.snpremoved\n",
    "outlieroutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.subsample.ind_outliers\n",
    "\n",
    "fastmode:      NO\n",
    "usenorm:       YES\n",
    "missingmode:   NO\n",
    "outliermode:   2 \n",
    "numchrom:      50\n",
    "numoutlieriter: 5 \n",
    "numoutlierevec: 10 \n",
    "lsqproject: NO \\\n",
    "    ''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "on_map_dom_subsample_PARFILE_path = '/home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.subsample.parfile'\n",
    "\n",
    "with open(on_map_dom_subsample_PARFILE_path, 'w') as OUTFILE:\n",
    "    OUTFILE.write('''\\\n",
    "genotypename:  /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.subsample.txt\n",
    "snpname:       /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.subsample.map\n",
    "indivname:     /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.ped\n",
    "\n",
    "evecoutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.subsample.evec\n",
    "evaloutname:    /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.subsample.eval\n",
    "grmoutname:      /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.subsample.rel\n",
    "snpweightoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.subsample.snpweights\n",
    "deletesnpoutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.subsample.snpremoved\n",
    "outlieroutname: /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.subsample.ind_outliers\n",
    "\n",
    "fastmode:      NO\n",
    "usenorm:       NO\n",
    "missingmode:   NO\n",
    "outliermode:   2 \n",
    "numchrom:      50\n",
    "numoutlieriter: 5 \n",
    "numoutlierevec: 10 \n",
    "lsqproject: NO \\\n",
    "    ''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### call SMARTPCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "!{SMARTPCA_path} -p {complete_codom_PARFILE_path} > \\\n",
    " /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.logfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "!{SMARTPCA_path} -p {complete_dom_PARFILE_path} > \\\n",
    " /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.logfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "!{SMARTPCA_path} -p {on_map_codom_PARFILE_path}> \\\n",
    " /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.logfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "!{SMARTPCA_path} -p {on_map_dom_PARFILE_path} > \\\n",
    " /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.logfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "!{SMARTPCA_path} -p {paralogs_dom_PARFILE_path}> \\\n",
    " /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/paralogs.dom.logfile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### on subsample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "!{SMARTPCA_path} -p {complete_codom_subsample_PARFILE_path}> \\\n",
    " /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.codom.subsample.logfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "!{SMARTPCA_path} -p {on_map_codom_subsample_PARFILE_path}> \\\n",
    " /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.codom.subsample.logfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "!{SMARTPCA_path} -p {complete_dom_subsample_PARFILE_path}> \\\n",
    " /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/complete.dom.subsample.logfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "!{SMARTPCA_path} -p {on_map_dom_subsample_PARFILE_path}> \\\n",
    " /home/ipseg/Desktop/waples/chum_populations/results/batch_4/EIGENSOFT/on_map.dom.subsample.logfile"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}