examples/pooling-imputing-ex.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "outputs": [],
   "source": [
    "import sys, os\n",
    "\n",
    "# force PYTHONPATH to look into the project directory for modules\n",
    "rootdir = os.path.dirname(os.path.dirname(os.getcwd()))\n",
    "sys.path.insert(0, rootdir)\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Pooling-imputation performance from real bead chip data\n",
    "\n",
    "Performs pooling simulation and imputation on data from the chromosome 20 of samples in the 1000 Genomes Project.\n",
    "The 1000 first markers have been selected and all unrelated individuals (2504 samples). \n",
    "These samples are randomly subdivided into a reference panel REF and a study population IMP to be pooled and imputed.\n",
    "Imputation is done with Beagle 4.1 and the default parameters.\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Check number of samples\n"
     ]
    },
    {
     "data": {
      "text/plain": "['2504']"
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('Check number of samples')\n",
    "%sx bcftools query -l ALL.chr20.snps.gt.vcf.gz | wc -l"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Check number of intersected markers\n"
     ]
    },
    {
     "data": {
      "text/plain": "['1000']"
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('Check number of intersected markers')\n",
    "%sx bcftools view -H ALL.chr20.snps.gt.vcf.gz | wc -l"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shuffle split into file for REF and IMP populations\n",
      "/home/camille/genotypooler/examples/REF.chr20.snps.gt.vcf.gz:\r\n",
      " File created? -> True\n",
      "/home/camille/genotypooler/examples/REF.chr20.snps.gt.vcf.gz:\r\n",
      " File indexed? -> True\n",
      "/home/camille/genotypooler/examples/IMP.chr20.snps.gt.vcf.gz:\r\n",
      " File created? -> True\n",
      "/home/camille/genotypooler/examples/IMP.chr20.snps.gt.vcf.gz:\r\n",
      " File indexed? -> True\n"
     ]
    }
   ],
   "source": [
    "print('Shuffle split into file for REF and IMP populations')\n",
    "from genotypooler.poolSNPs import poolvcf\n",
    "from genotypooler.poolSNPs.pooler import Design\n",
    "\n",
    "ds = Design()\n",
    "dm = ds.matrix\n",
    "\n",
    "sfsp = poolvcf.ShuffleSplitVCF(dm, 'ALL.chr20.snps.gt.vcf.gz', stu_size=0.1, wd=os.getcwd())\n",
    "sfsp.split_file('chr20.snps.gt.vcf.gz')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Apply pooling on the study population\n"
     ]
    },
    {
     "data": {
      "text/plain": "['model name\\t: Intel(R) Core(TM) i7-10850H CPU @ 2.70GHz']"
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('Apply pooling on the study population')\n",
    "# get processor characteristics on Linux-based OS\n",
    "%sx cat /proc/cpuinfo  | grep 'name'| uniq\n",
    "# model name\t: Intel(R) Core(TM) i7-7600U CPU @ 2.80GHz"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [
    {
     "data": {
      "text/plain": "['',\n '*******************************************************************************',\n 'Input file = IMP.chr20.snps.gt.vcf.gz',\n 'Output file = IMP.chr20.pooled.snps.gl.vcf.gz',\n '*******************************************************************************',\n '',\n 'Pooling data in /home/camille/genotypooler/examples/IMP.chr20.snps.gt.vcf.gz....',\n '1 variants processed in 000.01 sec..............................................',\n '',\n 'Writing metadata in /home/camille/genotypooler/examples/IMP.chr20.pooled.snps.gl.vcf',\n '',\n 'Writing data in /home/camille/genotypooler/examples/IMP.chr20.pooled.snps.gl.vcf',\n 'Writing data in /home/camille/genotypooler/examples/IMP.chr20.pooled.snps.gl.vcf: Done',\n '/home/camille/genotypooler/examples/IMP.chr20.pooled.snps.gl.vcf.gz:',\n ' File created? -> True',\n '/home/camille/genotypooler/examples/IMP.chr20.pooled.snps.gl.vcf.gz:',\n ' File indexed? -> True',\n 'Time for pooling 1000 variants = 8.324971047000872 sec',\n '',\n 'Time elapsed -->  8.341530646001047']"
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#pooling simulation\n",
    "%sx python3 -u pooling-ex.py IMP.chr20.snps.gt.vcf.gz IMP.chr20.pooled.snps.gl.vcf.gz GP\n",
    "# Time for pooling 1000 variants = 6.537596869000481 sec"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Impute missing genotypes in the pooled file\n"
     ]
    },
    {
     "data": {
      "text/plain": "['Contigs in the reference file',\n '.................................................................................',\n 'Chromosome  20    Startpos = 60343    Endpos = 99420',\n '',\n '',\n 'Check FORMAT field in files for imputation',\n '.................................................................................',\n 'FORMAT in reference panel:  GT',\n 'FORMAT in target:  GL',\n '',\n '',\n 'Check number of samples and number of markers in files for imputation',\n '.................................................................................',\n 'reference:',\n '2264',\n '',\n 'target:',\n '240',\n '',\n '',\n 'Phase reference and target with BEAGLE',\n '.................................................................................',\n 'Beagle .jar file used at: ../bin/beagle.11Mar19.69c.jar',\n '',\n 'FORMAT in the phased ref file: GT',\n 'beagle.11Mar19.69c.jar (version 4.1)',\n 'Copyright (C) 2014-2015 Brian L. Browning',\n 'Enter \"java -jar beagle.11Mar19.69c.jar\" for a summary of command line arguments.',\n 'Start time: 03:53 PM CET on 15 Nov 2022',\n '',\n 'Command line: java -Xmx3904m -jar beagle.jar',\n '  impute=false',\n '  gtgl=IMP.chr20.pooled.snps.gl.vcf.gz',\n '  out=IMP.chr20.pooled.unphased',\n '',\n 'No genetic map is specified: using 1 cM = 1 Mb',\n '',\n 'reference samples:       0',\n 'target samples:        240',\n '',\n 'Window 1 [ 20:60343-99420 ]',\n 'target markers:       1000',\n '',\n 'Starting burn-in iterations',\n '',\n 'Window=1 Iteration=1',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 2      max edges/level: 10',\n 'mean edges/node:  1.269  mean count/edge: 240',\n '',\n 'Window=1 Iteration=2',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 3      max edges/level: 10',\n 'mean edges/node:  1.253  mean count/edge: 160',\n '',\n 'Window=1 Iteration=3',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 3      max edges/level: 9',\n 'mean edges/node:  1.201  mean count/edge: 160',\n '',\n 'Window=1 Iteration=4',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 3      max edges/level: 10',\n 'mean edges/node:  1.196  mean count/edge: 160',\n '',\n 'Window=1 Iteration=5',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 3      max edges/level: 10',\n 'mean edges/node:  1.178  mean count/edge: 160',\n '',\n 'Window=1 Iteration=6',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 3      max edges/level: 10',\n 'mean edges/node:  1.170  mean count/edge: 160',\n '',\n 'Window=1 Iteration=7',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 3      max edges/level: 10',\n 'mean edges/node:  1.150  mean count/edge: 160',\n '',\n 'Window=1 Iteration=8',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 3      max edges/level: 11',\n 'mean edges/node:  1.152  mean count/edge: 160',\n '',\n 'Window=1 Iteration=9',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 3      max edges/level: 11',\n 'mean edges/node:  1.139  mean count/edge: 160',\n '',\n 'Window=1 Iteration=10',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 3      max edges/level: 11',\n 'mean edges/node:  1.131  mean count/edge: 160',\n '',\n 'Starting phasing iterations',\n '',\n 'Window=1 Iteration=11',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     1 second',\n 'DAG statistics',\n 'mean edges/level: 9      max edges/level: 17',\n 'mean edges/node:  1.057  mean count/edge: 53',\n '',\n 'Window=1 Iteration=12',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     2 seconds',\n 'DAG statistics',\n 'mean edges/level: 11     max edges/level: 24',\n 'mean edges/node:  1.042  mean count/edge: 44',\n '',\n 'Window=1 Iteration=13',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     1 second',\n 'DAG statistics',\n 'mean edges/level: 12     max edges/level: 22',\n 'mean edges/node:  1.039  mean count/edge: 40',\n '',\n 'Window=1 Iteration=14',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     1 second',\n 'DAG statistics',\n 'mean edges/level: 13     max edges/level: 26',\n 'mean edges/node:  1.034  mean count/edge: 37',\n '',\n 'Window=1 Iteration=15',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     1 second',\n 'DAG statistics',\n 'mean edges/level: 13     max edges/level: 25',\n 'mean edges/node:  1.034  mean count/edge: 37',\n '',\n 'Number of markers:                1000',\n 'Total time for building model: 4 seconds',\n 'Total time for sampling:       7 seconds',\n 'Total run time:                12 seconds',\n '',\n 'End time: 03:53 PM CET on 15 Nov 2022',\n 'beagle.11Mar19.69c.jar (version 4.1) finished',\n '',\n 'beagle.11Mar19.69c.jar (version 4.1)',\n 'Copyright (C) 2014-2015 Brian L. Browning',\n 'Enter \"java -jar beagle.11Mar19.69c.jar\" for a summary of command line arguments.',\n 'Start time: 03:53 PM CET on 15 Nov 2022',\n '',\n 'Command line: java -Xmx3904m -jar beagle.jar',\n '  impute=false',\n '  gt=IMP.chr20.pooled.unphased.vcf.gz',\n '  out=IMP.chr20.pooled.phased',\n '',\n 'No genetic map is specified: using 1 cM = 1 Mb',\n '',\n 'reference samples:       0',\n 'target samples:        240',\n '',\n 'Window 1 [ 20:60343-99420 ]',\n 'target markers:       1000',\n '',\n 'Starting burn-in iterations',\n '',\n 'Window=1 Iteration=1',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 5      max edges/level: 15',\n 'mean edges/node:  1.119  mean count/edge: 96',\n '',\n 'Window=1 Iteration=2',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 6      max edges/level: 18',\n 'mean edges/node:  1.085  mean count/edge: 80',\n '',\n 'Window=1 Iteration=3',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 6      max edges/level: 15',\n 'mean edges/node:  1.069  mean count/edge: 80',\n '',\n 'Window=1 Iteration=4',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 7      max edges/level: 15',\n 'mean edges/node:  1.062  mean count/edge: 69',\n '',\n 'Window=1 Iteration=5',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 7      max edges/level: 16',\n 'mean edges/node:  1.061  mean count/edge: 69',\n '',\n 'Window=1 Iteration=6',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 7      max edges/level: 16',\n 'mean edges/node:  1.060  mean count/edge: 69',\n '',\n 'Window=1 Iteration=7',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 7      max edges/level: 15',\n 'mean edges/node:  1.054  mean count/edge: 69',\n '',\n 'Window=1 Iteration=8',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 9      max edges/level: 16',\n 'mean edges/node:  1.050  mean count/edge: 53',\n '',\n 'Window=1 Iteration=9',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 8      max edges/level: 16',\n 'mean edges/node:  1.051  mean count/edge: 60',\n '',\n 'Window=1 Iteration=10',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     0 seconds',\n 'DAG statistics',\n 'mean edges/level: 8      max edges/level: 17',\n 'mean edges/node:  1.047  mean count/edge: 60',\n '',\n 'Starting phasing iterations',\n '',\n 'Window=1 Iteration=11',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     2 seconds',\n 'DAG statistics',\n 'mean edges/level: 11     max edges/level: 19',\n 'mean edges/node:  1.034  mean count/edge: 44',\n '',\n 'Window=1 Iteration=12',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     1 second',\n 'DAG statistics',\n 'mean edges/level: 14     max edges/level: 24',\n 'mean edges/node:  1.029  mean count/edge: 34',\n '',\n 'Window=1 Iteration=13',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     1 second',\n 'DAG statistics',\n 'mean edges/level: 13     max edges/level: 24',\n 'mean edges/node:  1.032  mean count/edge: 37',\n '',\n 'Window=1 Iteration=14',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     1 second',\n 'DAG statistics',\n 'mean edges/level: 16     max edges/level: 27',\n 'mean edges/node:  1.026  mean count/edge: 30',\n '',\n 'Window=1 Iteration=15',\n 'Time for building model:         0 seconds',\n 'Time for sampling (singles):     1 second',\n 'DAG statistics',\n 'mean edges/level: 14     max edges/level: 25',\n 'mean edges/node:  1.029  mean count/edge: 34',\n '',\n 'Number of markers:                1000',\n 'Total time for building model: 4 seconds',\n 'Total time for sampling:       6 seconds',\n 'Total run time:                10 seconds',\n '',\n 'End time: 03:53 PM CET on 15 Nov 2022',\n 'beagle.11Mar19.69c.jar (version 4.1) finished',\n 'FORMAT in the phased target file: GT',\n '',\n '',\n 'Deduplicate possibly duplicated markers',\n '.................................................................................',\n 'Lines   total/split/realigned/skipped:\\t1000/0/0/0',\n 'Lines   total/split/realigned/skipped:\\t1000/0/0/0',\n '',\n '',\n 'Unify reference and target markers with CONFORM-GT',\n '.................................................................................',\n 'conform-gt .jar file used at: ../bin/conform-gt.jar',\n '',\n '',\n '',\n 'Impute target from reference with BEAGLE',\n '.................................................................................',\n 'Beagle .jar file used at: ../bin/beagle.11Mar19.69c.jar',\n '',\n 'beagle.11Mar19.69c.jar (version 4.1)',\n 'Copyright (C) 2014-2015 Brian L. Browning',\n 'Enter \"java -jar beagle.11Mar19.69c.jar\" for a summary of command line arguments.',\n 'Start time: 03:53 PM CET on 15 Nov 2022',\n '',\n 'Command line: java -Xmx3904m -jar beagle.jar',\n '  gt=IMP.chr20.pooled.cfgt.vcf.gz',\n '  ref=REF.chr20.phased.vcf.gz',\n '  impute=true',\n '  gprobs=true',\n '  out=IMP.chr20.pooled.imputed',\n '',\n 'No genetic map is specified: using 1 cM = 1 Mb',\n '',\n 'reference samples:    2264',\n 'target samples:        240',\n '',\n 'Window 1 [ 20:60343-99420 ]',\n 'reference markers:    1000',\n 'target markers:        999',\n '',\n 'Imputation time (this window): 0 seconds',\n '',\n 'Number of reference markers:      1000',\n 'Number of target markers:          999',\n 'Total time for imputation:     0 seconds',\n 'Total run time:                1 second',\n '',\n 'End time: 03:53 PM CET on 15 Nov 2022',\n 'beagle.11Mar19.69c.jar (version 4.1) finished',\n 'FORMAT in the imputed target file: GT:DS:GP',\n '',\n '',\n 'Cleaning directory from log files',\n '.................................................................................',\n 'Done.',\n '']"
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('Impute missing genotypes in the pooled file')\n",
    "%sx cp ../bin/beagle_pipeline.sh ./\n",
    "%sx bash beagle_pipeline.sh"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "outputs": [],
   "source": [
    "# Verify files created at the different phasing and imputation steps\n",
    "assert os.path.exists('IMP.chr20.pooled.unphased.vcf.gz')\n",
    "assert os.path.exists('IMP.chr20.pooled.phased.vcf.gz')\n",
    "assert os.path.exists('IMP.chr20.pooled.phased.dedup.vcf.gz')\n",
    "assert os.path.exists('IMP.chr20.pooled.imputed.vcf.gz')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Compute results with customized metrics\n",
    "Show classification task-based metrics and classical genotype imputation metrics"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [
    {
     "data": {
      "text/plain": "['/home/camille',\n '/home/camille/genotypooler/poolSNPs/metrics/quality.py:328: RuntimeWarning: divide by zero encountered in myprodlog_numba',\n '  result = myprodlog_numba(dftrue.values, dfpred.values)',\n '/home/camille/genotypooler/poolSNPs/metrics/quality.py:328: RuntimeWarning: divide by zero encountered in myprodlog_numba',\n '  result = myprodlog_numba(dftrue.values, dfpred.values)',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/scipy/stats/stats.py:3845: PearsonRConstantInputWarning: An input array is constant; the correlation coefficent is not defined.',\n '  warnings.warn(PearsonRConstantInputWarning())',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n '/home/camille/1000Genomes/src/VCFPooling/venv3.6/lib/python3.6/site-packages/sklearn/metrics/_classification.py:1248: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.',\n '  _warn_prf(average, modifier, msg_start, len(result))',\n 'Figure(2000x800)']"
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%sx python3 -u imputation_quality.py ./ IMP.chr20.snps.gt.vcf.gz IMP.chr20.pooled.imputed.vcf.gz ../bin/gt_to_log_gl.sh 'id'"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [],
   "source": [
    "# Verify files created at the different phasing and imputation steps\n",
    "assert os.path.exists('imputation_quality_gtgl.pdf')\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "outputs": [],
   "source": [
    "# # if bcftools is configured for python 2.7 usage \n",
    "# print('Plotting results with bcftools stats')\n",
    "# %sx deactivate\n",
    "# # bcftools stats needs python 2.7\n",
    "# %sx bcftools stats --af-bins 0.01,0.02,0.04,0.08,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.98 --collapse snps -S study.population IMP.chr20.pooled.imputed.vcf.gz IMP.chr20.snps.gt.vcf.gz > filestats.vchk\n",
    "# %sx plot-vcfstats -p bcftoolstats -s filestats.vchk\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}