SIG_Genome_Cache_Workflow.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Systems Immunogenetics Project\n",
    "\n",
    "## QTL Mapping: Creating the Genome Cache\n",
    "\n",
    "### McWeeney Lab, Oregon Health & Science University\n",
    "\n",
    "#### Author: Michael Mooney (mooneymi@ohsu.edu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "This document will walk through the steps for creating a genome cache for the RIX lines.\n",
    "\n",
    "Required Files:\n",
    "\n",
    "1. This notebook** (`SIG_Genome_Cache_Workflow.ipynb`): [[Download here]](https://raw.githubusercontent.com/biodev/SIG/master/SIG_Genome_Cache_Workflow.ipynb)\n",
    "2. The R script `rix_qtl_mapping_functions.r`: [[Download here]](https://raw.githubusercontent.com/biodev/SIG/master/scripts/rix_qtl_mapping_functions.r)\n",
    "\n",
    "** Note: this notebook can also be downloaded as an R script (only the code blocks seen below will be included): [[Download R script here]](https://raw.githubusercontent.com/biodev/SIG/master/SIG_Genome_Cache_Workflow.r)\n",
    "\n",
    "Required R Libraries:\n",
    "\n",
    "1. `gdata`: [https://cran.r-project.org/web/packages/gdata/index.html](https://cran.r-project.org/web/packages/gdata/index.html)\n",
    "2. `DOQTL`: [https://www.bioconductor.org/packages/release/bioc/html/DOQTL.html](https://www.bioconductor.org/packages/release/bioc/html/DOQTL.html)\n",
    "3. `GenomicRanges`: [https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html](https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html)\n",
    "4. `VariantAnnotation`: [https://bioconductor.org/packages/release/bioc/html/VariantAnnotation.html](https://bioconductor.org/packages/release/bioc/html/VariantAnnotation.html)\n",
    "5. `foreach`: [https://cran.r-project.org/web/packages/foreach/index.html](https://cran.r-project.org/web/packages/foreach/index.html)\n",
    "6. `doParallel`: [https://cran.r-project.org/web/packages/doParallel/index.html](https://cran.r-project.org/web/packages/doParallel/index.html)\n",
    "7. `RColorBrewer`: [https://cran.r-project.org/web/packages/RColorBrewer/index.html](https://cran.r-project.org/web/packages/RColorBrewer/index.html)\n",
    "\n",
    "**All code is available on GitHub: [https://github.com/biodev/SIG](https://github.com/biodev/SIG)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1. Load Necessary R Functions and Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "gdata: read.xls support for 'XLS' (Excel 97-2004) files ENABLED.\n",
      "\n",
      "gdata: read.xls support for 'XLSX' (Excel 2007+) files ENABLED.\n",
      "\n",
      "Attaching package: ‘gdata’\n",
      "\n",
      "The following object is masked from ‘package:stats’:\n",
      "\n",
      "    nobs\n",
      "\n",
      "The following object is masked from ‘package:utils’:\n",
      "\n",
      "    object.size\n",
      "\n",
      "The following object is masked from ‘package:base’:\n",
      "\n",
      "    startsWith\n",
      "\n",
      "Loading required package: BSgenome.Mmusculus.UCSC.mm10\n",
      "Loading required package: BSgenome\n",
      "Loading required package: BiocGenerics\n",
      "Loading required package: parallel\n",
      "\n",
      "Attaching package: ‘BiocGenerics’\n",
      "\n",
      "The following objects are masked from ‘package:parallel’:\n",
      "\n",
      "    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,\n",
      "    clusterExport, clusterMap, parApply, parCapply, parLapply,\n",
      "    parLapplyLB, parRapply, parSapply, parSapplyLB\n",
      "\n",
      "The following object is masked from ‘package:gdata’:\n",
      "\n",
      "    combine\n",
      "\n",
      "The following objects are masked from ‘package:stats’:\n",
      "\n",
      "    IQR, mad, xtabs\n",
      "\n",
      "The following objects are masked from ‘package:base’:\n",
      "\n",
      "    anyDuplicated, append, as.data.frame, as.vector, cbind, colnames,\n",
      "    do.call, duplicated, eval, evalq, Filter, Find, get, grep, grepl,\n",
      "    intersect, is.unsorted, lapply, lengths, Map, mapply, match, mget,\n",
      "    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,\n",
      "    rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,\n",
      "    union, unique, unlist, unsplit\n",
      "\n",
      "Loading required package: S4Vectors\n",
      "Loading required package: stats4\n",
      "Loading required package: IRanges\n",
      "\n",
      "Attaching package: ‘IRanges’\n",
      "\n",
      "The following object is masked from ‘package:gdata’:\n",
      "\n",
      "    trim\n",
      "\n",
      "Loading required package: GenomeInfoDb\n",
      "Loading required package: GenomicRanges\n",
      "Loading required package: Biostrings\n",
      "Loading required package: XVector\n",
      "Loading required package: rtracklayer\n",
      "Loading required package: VariantAnnotation\n",
      "Loading required package: SummarizedExperiment\n",
      "Loading required package: Biobase\n",
      "Welcome to Bioconductor\n",
      "\n",
      "    Vignettes contain introductory material; view with\n",
      "    'browseVignettes()'. To cite Bioconductor, see\n",
      "    'citation(\"Biobase\")', and for packages 'citation(\"pkgname\")'.\n",
      "\n",
      "Loading required package: Rsamtools\n",
      "\n",
      "Attaching package: ‘VariantAnnotation’\n",
      "\n",
      "The following object is masked from ‘package:base’:\n",
      "\n",
      "    tabulate\n",
      "\n",
      "Loading required package: iterators\n"
     ]
    }
   ],
   "source": [
    "## Load R functions and libraries\n",
    "source('./scripts/rix_qtl_mapping_functions.r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2. Create Genome Cache\n",
    "\n",
    "If the RIX list below is complete, you should only have to create the cache once (one for males and one for females). You can then subset the cache based on the particular mapping population being studied (as shown in the QTL Mapping Workflow)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>Mating</th><th scope=col>CC_Mating</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>1</th><td>477x16912  </td><td>CC036xCC051</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>1566x8002  </td><td>CC021xCC032</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>1566x8043  </td><td>CC021xCC023</td></tr>\n",
       "\t<tr><th scope=row>4</th><td>3015x5306  </td><td>CC074xCC062</td></tr>\n",
       "\t<tr><th scope=row>5</th><td>3015x5358  </td><td>CC074xCC058</td></tr>\n",
       "\t<tr><th scope=row>6</th><td>3015x15156 </td><td>CC074xCC002</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|ll}\n",
       "  & Mating & CC\\_Mating\\\\\n",
       "\\hline\n",
       "\t1 & 477x16912   & CC036xCC051\\\\\n",
       "\t2 & 1566x8002   & CC021xCC032\\\\\n",
       "\t3 & 1566x8043   & CC021xCC023\\\\\n",
       "\t4 & 3015x5306   & CC074xCC062\\\\\n",
       "\t5 & 3015x5358   & CC074xCC058\\\\\n",
       "\t6 & 3015x15156  & CC074xCC002\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "      Mating   CC_Mating\n",
       "1  477x16912 CC036xCC051\n",
       "2  1566x8002 CC021xCC032\n",
       "3  1566x8043 CC021xCC023\n",
       "4  3015x5306 CC074xCC062\n",
       "5  3015x5358 CC074xCC058\n",
       "6 3015x15156 CC074xCC002"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## First get all RIX matings\n",
    "rix_list = read.xls('./data/BOBs_Corrected_16_06_22.xlsx', header=T, as.is=T)\n",
    "rix_list = rix_list[,1:2]\n",
    "colnames(rix_list) = c('Mating', 'CC_Mating')\n",
    "head(rix_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "72"
      ],
      "text/latex": [
       "72"
      ],
      "text/markdown": [
       "72"
      ],
      "text/plain": [
       "[1] 72"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>'CC001-Uncb38V01.csv'</li>\n",
       "\t<li>'CC002-Uncb38V01.csv'</li>\n",
       "\t<li>'CC003-Uncb38V01.csv'</li>\n",
       "\t<li>'CC004-TauUncb38V01.csv'</li>\n",
       "\t<li>'CC005-TauUncb38V01.csv'</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 'CC001-Uncb38V01.csv'\n",
       "\\item 'CC002-Uncb38V01.csv'\n",
       "\\item 'CC003-Uncb38V01.csv'\n",
       "\\item 'CC004-TauUncb38V01.csv'\n",
       "\\item 'CC005-TauUncb38V01.csv'\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 'CC001-Uncb38V01.csv'\n",
       "2. 'CC002-Uncb38V01.csv'\n",
       "3. 'CC003-Uncb38V01.csv'\n",
       "4. 'CC004-TauUncb38V01.csv'\n",
       "5. 'CC005-TauUncb38V01.csv'\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] \"CC001-Uncb38V01.csv\"    \"CC002-Uncb38V01.csv\"    \"CC003-Uncb38V01.csv\"   \n",
       "[4] \"CC004-TauUncb38V01.csv\" \"CC005-TauUncb38V01.csv\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Set directory containing founder 36-state probability files, and get filenames\n",
    "cc_dir = '/Users/mooneymi/Documents/SIG/Mapping/CC'\n",
    "cc_prob_files = list.files(cc_dir, pattern=\"CC...-.*b38.*\\\\.csv\")\n",
    "length(cc_prob_files)\n",
    "cc_prob_files[1:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "69"
      ],
      "text/latex": [
       "69"
      ],
      "text/markdown": [
       "69"
      ],
      "text/plain": [
       "[1] 69"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<dl class=dl-horizontal>\n",
       "\t<dt>477</dt>\n",
       "\t\t<dd>'CC036'</dd>\n",
       "\t<dt>16912</dt>\n",
       "\t\t<dd>'CC051'</dd>\n",
       "\t<dt>1566</dt>\n",
       "\t\t<dd>'CC021'</dd>\n",
       "\t<dt>8002</dt>\n",
       "\t\t<dd>'CC032'</dd>\n",
       "\t<dt>8043</dt>\n",
       "\t\t<dd>'CC023'</dd>\n",
       "</dl>\n"
      ],
      "text/latex": [
       "\\begin{description*}\n",
       "\\item[477] 'CC036'\n",
       "\\item[16912] 'CC051'\n",
       "\\item[1566] 'CC021'\n",
       "\\item[8002] 'CC032'\n",
       "\\item[8043] 'CC023'\n",
       "\\end{description*}\n"
      ],
      "text/markdown": [
       "477\n",
       ":   'CC036'16912\n",
       ":   'CC051'1566\n",
       ":   'CC021'8002\n",
       ":   'CC032'8043\n",
       ":   'CC023'\n",
       "\n"
      ],
      "text/plain": [
       "    477   16912    1566    8002    8043 \n",
       "\"CC036\" \"CC051\" \"CC021\" \"CC032\" \"CC023\" "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Get vectors of parental strains\n",
    "rix_strains = unlist(strsplit(rix_list$Mating, 'x'))\n",
    "cc_strains = unlist(strsplit(rix_list$CC_Mating, 'x'))\n",
    "\n",
    "## Check strains (some cc_strains will be duplicated -- re-derived lines)\n",
    "dup_idx = duplicated(rix_strains)\n",
    "cc_strains = cc_strains[!dup_idx]\n",
    "names(cc_strains) = rix_strains[!dup_idx]\n",
    "length(cc_strains)\n",
    "cc_strains[1:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "0"
      ],
      "text/latex": [
       "0"
      ],
      "text/markdown": [
       "0"
      ],
      "text/plain": [
       "[1] 0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Check that all RIX strains mapped to a CC strains\n",
    "sum(is.na(cc_strains))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Set directory containing the 8-state probability files\n",
    "## These files can be created using the collapse_probs() function (see next code block)\n",
    "rix_dir = '/Users/mooneymi/Documents/SIG/Mapping/RIX'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## If you haven't already done so, the following code can be used to \n",
    "## create the 8-state probability files (you should only have to do this once)\n",
    "## File names will be the RIX IDs (not CC IDs). \n",
    "## Note: you may have to make copies of some files for re-derived lines \n",
    "## (e.g. 8042 and 18042 are the same), if this wasn't handled above\n",
    "# for (i in 1:length(cc_strains)) {\n",
    "#     ## Get RIX ID\n",
    "#     rix = names(cc_strains[i])\n",
    "#     \n",
    "#     ## Create and save the 8-state probabilities for each parental strain\n",
    "#     cc_file = file.path(cc_dir, cc_prob_files[grep(paste0(cc_strains[i], '-.*\\\\.csv'), cc_prob_files)])\n",
    "#     cc_8state = collapse_probs(cc_file)\n",
    "#     \n",
    "#     ## Save Y and M chromosomes separately\n",
    "#     autosome_x_markers = cc_8state$marker[!cc_8state$chromosome %in% c('Y', 'M')]\n",
    "#     y_m_markers = cc_8state$marker[cc_8state$chromosome %in% c('Y', 'M')]\n",
    "#     file_name = file.path(rix_dir, paste0(rix, '.csv'))\n",
    "#     write.table(cc_8state[autosome_x_markers, ], file=file_name, row.names=F, col.names=T, sep=',', quote=F)\n",
    "#     file_name = file.path(rix_dir, paste0(rix, '_Y_M.csv'))\n",
    "#     write.table(cc_8state[y_m_markers, ], file=file_name, row.names=F, col.names=T, sep=',', quote=F)\n",
    "# }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 1\n",
      "[1] 2\n",
      "[1] 3\n",
      "[1] 4\n",
      "[1] 5\n",
      "[1] 6\n",
      "[1] 7\n",
      "[1] 8\n",
      "[1] 9\n",
      "[1] 10\n",
      "[1] 11\n",
      "[1] 12\n",
      "[1] 13\n",
      "[1] 14\n",
      "[1] 15\n",
      "[1] 16\n",
      "[1] 17\n",
      "[1] 18\n",
      "[1] 19\n",
      "[1] 20\n",
      "[1] 21\n",
      "[1] 22\n",
      "[1] 23\n",
      "[1] 24\n",
      "[1] 25\n",
      "[1] 26\n",
      "[1] 27\n",
      "[1] 28\n",
      "[1] 29\n",
      "[1] 30\n",
      "[1] 31\n",
      "[1] 32\n",
      "[1] 33\n",
      "[1] 34\n",
      "[1] 35\n",
      "[1] 36\n",
      "[1] 37\n",
      "[1] 38\n",
      "[1] 39\n",
      "[1] 40\n",
      "[1] 41\n",
      "[1] 42\n",
      "[1] 43\n",
      "[1] 44\n",
      "[1] 45\n",
      "[1] 46\n",
      "[1] 47\n",
      "[1] 48\n",
      "[1] 49\n",
      "[1] 50\n",
      "[1] 51\n",
      "[1] 52\n",
      "[1] 53\n",
      "[1] 54\n",
      "[1] 55\n",
      "[1] 56\n",
      "[1] 57\n",
      "[1] 58\n",
      "[1] 59\n",
      "[1] 60\n",
      "[1] 61\n",
      "[1] 62\n",
      "[1] 63\n",
      "[1] 64\n",
      "[1] 65\n",
      "[1] 66\n",
      "[1] 67\n",
      "[1] 68\n",
      "[1] 69\n",
      "[1] 70\n",
      "[1] 71\n",
      "[1] 72\n",
      "[1] 73\n",
      "[1] 74\n",
      "[1] 75\n",
      "[1] 76\n",
      "[1] 77\n",
      "[1] 78\n",
      "[1] 79\n",
      "[1] 80\n",
      "[1] 81\n",
      "[1] 82\n",
      "[1] 83\n",
      "[1] 84\n",
      "[1] 85\n",
      "[1] 86\n",
      "[1] 87\n",
      "[1] 88\n",
      "[1] 89\n",
      "[1] 90\n",
      "[1] 91\n",
      "[1] 92\n",
      "[1] 93\n",
      "[1] 94\n",
      "[1] 95\n",
      "[1] 96\n",
      "[1] 97\n",
      "[1] 98\n",
      "[1] 99\n",
      "[1] 100\n",
      "[1] 101\n",
      "[1] 102\n",
      "[1] 103\n",
      "[1] 104\n",
      "[1] 105\n",
      "[1] 106\n",
      "[1] 107\n",
      "[1] 108\n",
      "[1] 109\n",
      "[1] 110\n",
      "[1] 111\n",
      "[1] 112\n",
      "[1] 113\n",
      "[1] 114\n",
      "[1] 115\n",
      "[1] 116\n",
      "[1] 117\n",
      "[1] 118\n",
      "[1] 119\n",
      "[1] 120\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>120</li>\n",
       "\t<li>8</li>\n",
       "\t<li>76606</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 120\n",
       "\\item 8\n",
       "\\item 76606\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 120\n",
       "2. 8\n",
       "3. 76606\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1]   120     8 76606"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>'samples'</li>\n",
       "\t<li>'founders'</li>\n",
       "\t<li>'markers'</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 'samples'\n",
       "\\item 'founders'\n",
       "\\item 'markers'\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 'samples'\n",
       "2. 'founders'\n",
       "3. 'markers'\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] \"samples\"  \"founders\" \"markers\" "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Create 3D probability array - Males\n",
    "## Note: the number of matings included will be printed\n",
    "model.probs = make_rix_model_probs(rix_list$Mating, rix_dir, 'M')\n",
    "dim(model.probs)\n",
    "names(dimnames(model.probs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>UNC6</th><th scope=col>JAX00000010</th><th scope=col>JAX00240603</th><th scope=col>JAX00240610</th><th scope=col>JAX00240613</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>A</th><td>3.00075e-06</td><td>2.62500e-06</td><td>2.25000e-06</td><td>1.87500e-06</td><td>1.50000e-06</td></tr>\n",
       "\t<tr><th scope=row>B</th><td>3.751e-06</td><td>3.375e-06</td><td>3.000e-06</td><td>2.625e-06</td><td>2.250e-06</td></tr>\n",
       "\t<tr><th scope=row>C</th><td>2.292e-06</td><td>1.875e-06</td><td>1.500e-06</td><td>1.125e-06</td><td>7.500e-07</td></tr>\n",
       "\t<tr><th scope=row>D</th><td>2.292e-06</td><td>1.875e-06</td><td>1.500e-06</td><td>1.125e-06</td><td>7.500e-07</td></tr>\n",
       "\t<tr><th scope=row>E</th><td>0.4999923</td><td>0.4999936</td><td>0.4999947</td><td>0.4999955</td><td>0.4999962</td></tr>\n",
       "\t<tr><th scope=row>F</th><td>0.4999948</td><td>0.4999959</td><td>0.4999970</td><td>0.4999978</td><td>0.4999985</td></tr>\n",
       "\t<tr><th scope=row>G</th><td>7.92e-07</td><td>3.75e-07</td><td>0.00e+00</td><td>0.00e+00</td><td>0.00e+00</td></tr>\n",
       "\t<tr><th scope=row>H</th><td>7.5075e-07</td><td>3.7500e-07</td><td>0.0000e+00</td><td>0.0000e+00</td><td>0.0000e+00</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lllll}\n",
       "  & UNC6 & JAX00000010 & JAX00240603 & JAX00240610 & JAX00240613\\\\\n",
       "\\hline\n",
       "\tA & 3.00075e-06 & 2.62500e-06 & 2.25000e-06 & 1.87500e-06 & 1.50000e-06\\\\\n",
       "\tB & 3.751e-06 & 3.375e-06 & 3.000e-06 & 2.625e-06 & 2.250e-06\\\\\n",
       "\tC & 2.292e-06 & 1.875e-06 & 1.500e-06 & 1.125e-06 & 7.500e-07\\\\\n",
       "\tD & 2.292e-06 & 1.875e-06 & 1.500e-06 & 1.125e-06 & 7.500e-07\\\\\n",
       "\tE & 0.4999923 & 0.4999936 & 0.4999947 & 0.4999955 & 0.4999962\\\\\n",
       "\tF & 0.4999948 & 0.4999959 & 0.4999970 & 0.4999978 & 0.4999985\\\\\n",
       "\tG & 7.92e-07 & 3.75e-07 & 0.00e+00 & 0.00e+00 & 0.00e+00\\\\\n",
       "\tH & 7.5075e-07 & 3.7500e-07 & 0.0000e+00 & 0.0000e+00 & 0.0000e+00\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "1. 3.00075e-06\n",
       "2. 3.751e-06\n",
       "3. 2.292e-06\n",
       "4. 2.292e-06\n",
       "5. 0.4999923325\n",
       "6. 0.499994789\n",
       "7. 7.92e-07\n",
       "8. 7.5075e-07\n",
       "9. 2.625e-06\n",
       "10. 3.375e-06\n",
       "11. 1.875e-06\n",
       "12. 1.875e-06\n",
       "13. 0.499993625\n",
       "14. 0.499995875\n",
       "15. 3.75e-07\n",
       "16. 3.75e-07\n",
       "17. 2.25e-06\n",
       "18. 3e-06\n",
       "19. 1.5e-06\n",
       "20. 1.5e-06\n",
       "21. 0.49999475\n",
       "22. 0.499997\n",
       "23. 0\n",
       "24. 0\n",
       "25. 1.875e-06\n",
       "26. 2.625e-06\n",
       "27. 1.125e-06\n",
       "28. 1.125e-06\n",
       "29. 0.4999955\n",
       "30. 0.49999775\n",
       "31. 0\n",
       "32. 0\n",
       "33. 1.5e-06\n",
       "34. 2.25e-06\n",
       "35. 7.5e-07\n",
       "36. 7.5e-07\n",
       "37. 0.49999625\n",
       "38. 0.4999985\n",
       "39. 0\n",
       "40. 0\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "        markers\n",
       "founders         UNC6 JAX00000010 JAX00240603 JAX00240610 JAX00240613\n",
       "       A 3.000750e-06 0.000002625  0.00000225 0.000001875  0.00000150\n",
       "       B 3.751000e-06 0.000003375  0.00000300 0.000002625  0.00000225\n",
       "       C 2.292000e-06 0.000001875  0.00000150 0.000001125  0.00000075\n",
       "       D 2.292000e-06 0.000001875  0.00000150 0.000001125  0.00000075\n",
       "       E 4.999923e-01 0.499993625  0.49999475 0.499995500  0.49999625\n",
       "       F 4.999948e-01 0.499995875  0.49999700 0.499997750  0.49999850\n",
       "       G 7.920000e-07 0.000000375  0.00000000 0.000000000  0.00000000\n",
       "       H 7.507500e-07 0.000000375  0.00000000 0.000000000  0.00000000"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Check model.probs object\n",
    "model.probs[1,,1:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Get vector of all markers\n",
    "all_markers = dimnames(model.probs)[[3]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "6"
      ],
      "text/latex": [
       "6"
      ],
      "text/markdown": [
       "6"
      ],
      "text/plain": [
       "[1] 6"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Check if any markers sum to 0\n",
    "markers_zero_idx = which(apply(model.probs, 1, colSums) == 0)\n",
    "markers_zero = rep(all_markers, length(samples))[markers_zero_idx]\n",
    "markers_zero = markers_zero[!is.na(markers_zero)]\n",
    "markers_zero = unique(markers_zero)\n",
    "length(markers_zero)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>120</li>\n",
       "\t<li>8</li>\n",
       "\t<li>76600</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 120\n",
       "\\item 8\n",
       "\\item 76600\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 120\n",
       "2. 8\n",
       "3. 76600\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1]   120     8 76600"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Remove markers with all zeros\n",
    "model.probs = model.probs[, , setdiff(all_markers, markers_zero)]\n",
    "dim(model.probs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>UNC6</th><th scope=col>JAX00000010</th><th scope=col>JAX00240603</th><th scope=col>JAX00240610</th><th scope=col>JAX00240613</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>A</th><td>3.00075e-06</td><td>2.62500e-06</td><td>2.25000e-06</td><td>1.87500e-06</td><td>1.50000e-06</td></tr>\n",
       "\t<tr><th scope=row>B</th><td>3.751e-06</td><td>3.375e-06</td><td>3.000e-06</td><td>2.625e-06</td><td>2.250e-06</td></tr>\n",
       "\t<tr><th scope=row>C</th><td>2.292e-06</td><td>1.875e-06</td><td>1.500e-06</td><td>1.125e-06</td><td>7.500e-07</td></tr>\n",
       "\t<tr><th scope=row>D</th><td>2.292e-06</td><td>1.875e-06</td><td>1.500e-06</td><td>1.125e-06</td><td>7.500e-07</td></tr>\n",
       "\t<tr><th scope=row>E</th><td>0.4999923</td><td>0.4999936</td><td>0.4999947</td><td>0.4999955</td><td>0.4999962</td></tr>\n",
       "\t<tr><th scope=row>F</th><td>0.4999948</td><td>0.4999959</td><td>0.4999970</td><td>0.4999978</td><td>0.4999985</td></tr>\n",
       "\t<tr><th scope=row>G</th><td>7.92e-07</td><td>3.75e-07</td><td>1.00e-20</td><td>1.00e-20</td><td>1.00e-20</td></tr>\n",
       "\t<tr><th scope=row>H</th><td>7.5075e-07</td><td>3.7500e-07</td><td>1.0000e-20</td><td>1.0000e-20</td><td>1.0000e-20</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lllll}\n",
       "  & UNC6 & JAX00000010 & JAX00240603 & JAX00240610 & JAX00240613\\\\\n",
       "\\hline\n",
       "\tA & 3.00075e-06 & 2.62500e-06 & 2.25000e-06 & 1.87500e-06 & 1.50000e-06\\\\\n",
       "\tB & 3.751e-06 & 3.375e-06 & 3.000e-06 & 2.625e-06 & 2.250e-06\\\\\n",
       "\tC & 2.292e-06 & 1.875e-06 & 1.500e-06 & 1.125e-06 & 7.500e-07\\\\\n",
       "\tD & 2.292e-06 & 1.875e-06 & 1.500e-06 & 1.125e-06 & 7.500e-07\\\\\n",
       "\tE & 0.4999923 & 0.4999936 & 0.4999947 & 0.4999955 & 0.4999962\\\\\n",
       "\tF & 0.4999948 & 0.4999959 & 0.4999970 & 0.4999978 & 0.4999985\\\\\n",
       "\tG & 7.92e-07 & 3.75e-07 & 1.00e-20 & 1.00e-20 & 1.00e-20\\\\\n",
       "\tH & 7.5075e-07 & 3.7500e-07 & 1.0000e-20 & 1.0000e-20 & 1.0000e-20\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "1. 3.00075e-06\n",
       "2. 3.751e-06\n",
       "3. 2.292e-06\n",
       "4. 2.292e-06\n",
       "5. 0.4999923325\n",
       "6. 0.499994789\n",
       "7. 7.92e-07\n",
       "8. 7.5075e-07\n",
       "9. 2.625e-06\n",
       "10. 3.375e-06\n",
       "11. 1.875e-06\n",
       "12. 1.875e-06\n",
       "13. 0.499993625\n",
       "14. 0.499995875\n",
       "15. 3.75e-07\n",
       "16. 3.75e-07\n",
       "17. 2.25e-06\n",
       "18. 3e-06\n",
       "19. 1.5e-06\n",
       "20. 1.5e-06\n",
       "21. 0.49999475\n",
       "22. 0.499997\n",
       "23. 1e-20\n",
       "24. 1e-20\n",
       "25. 1.875e-06\n",
       "26. 2.625e-06\n",
       "27. 1.125e-06\n",
       "28. 1.125e-06\n",
       "29. 0.4999955\n",
       "30. 0.49999775\n",
       "31. 1e-20\n",
       "32. 1e-20\n",
       "33. 1.5e-06\n",
       "34. 2.25e-06\n",
       "35. 7.5e-07\n",
       "36. 7.5e-07\n",
       "37. 0.49999625\n",
       "38. 0.4999985\n",
       "39. 1e-20\n",
       "40. 1e-20\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "        markers\n",
       "founders         UNC6 JAX00000010  JAX00240603  JAX00240610  JAX00240613\n",
       "       A 3.000750e-06 0.000002625 2.250000e-06 1.875000e-06 1.500000e-06\n",
       "       B 3.751000e-06 0.000003375 3.000000e-06 2.625000e-06 2.250000e-06\n",
       "       C 2.292000e-06 0.000001875 1.500000e-06 1.125000e-06 7.500000e-07\n",
       "       D 2.292000e-06 0.000001875 1.500000e-06 1.125000e-06 7.500000e-07\n",
       "       E 4.999923e-01 0.499993625 4.999947e-01 4.999955e-01 4.999962e-01\n",
       "       F 4.999948e-01 0.499995875 4.999970e-01 4.999978e-01 4.999985e-01\n",
       "       G 7.920000e-07 0.000000375 1.000000e-20 1.000000e-20 1.000000e-20\n",
       "       H 7.507500e-07 0.000000375 1.000000e-20 1.000000e-20 1.000000e-20"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Make all probabilities non-zero\n",
    "model.probs[model.probs == 0] = 1e-20\n",
    "model.probs[1,,1:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Save genome cache - Males\n",
    "mapping_dir = '/Users/mooneymi/Documents/SIG/Mapping'\n",
    "save(model.probs, file=file.path(mapping_dir, 'rix_universal_model_prob_males_27-Jun-2016.rda'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 1\n",
      "[1] 2\n",
      "[1] 3\n",
      "[1] 4\n",
      "[1] 5\n",
      "[1] 6\n",
      "[1] 7\n",
      "[1] 8\n",
      "[1] 9\n",
      "[1] 10\n",
      "[1] 11\n",
      "[1] 12\n",
      "[1] 13\n",
      "[1] 14\n",
      "[1] 15\n",
      "[1] 16\n",
      "[1] 17\n",
      "[1] 18\n",
      "[1] 19\n",
      "[1] 20\n",
      "[1] 21\n",
      "[1] 22\n",
      "[1] 23\n",
      "[1] 24\n",
      "[1] 25\n",
      "[1] 26\n",
      "[1] 27\n",
      "[1] 28\n",
      "[1] 29\n",
      "[1] 30\n",
      "[1] 31\n",
      "[1] 32\n",
      "[1] 33\n",
      "[1] 34\n",
      "[1] 35\n",
      "[1] 36\n",
      "[1] 37\n",
      "[1] 38\n",
      "[1] 39\n",
      "[1] 40\n",
      "[1] 41\n",
      "[1] 42\n",
      "[1] 43\n",
      "[1] 44\n",
      "[1] 45\n",
      "[1] 46\n",
      "[1] 47\n",
      "[1] 48\n",
      "[1] 49\n",
      "[1] 50\n",
      "[1] 51\n",
      "[1] 52\n",
      "[1] 53\n",
      "[1] 54\n",
      "[1] 55\n",
      "[1] 56\n",
      "[1] 57\n",
      "[1] 58\n",
      "[1] 59\n",
      "[1] 60\n",
      "[1] 61\n",
      "[1] 62\n",
      "[1] 63\n",
      "[1] 64\n",
      "[1] 65\n",
      "[1] 66\n",
      "[1] 67\n",
      "[1] 68\n",
      "[1] 69\n",
      "[1] 70\n",
      "[1] 71\n",
      "[1] 72\n",
      "[1] 73\n",
      "[1] 74\n",
      "[1] 75\n",
      "[1] 76\n",
      "[1] 77\n",
      "[1] 78\n",
      "[1] 79\n",
      "[1] 80\n",
      "[1] 81\n",
      "[1] 82\n",
      "[1] 83\n",
      "[1] 84\n",
      "[1] 85\n",
      "[1] 86\n",
      "[1] 87\n",
      "[1] 88\n",
      "[1] 89\n",
      "[1] 90\n",
      "[1] 91\n",
      "[1] 92\n",
      "[1] 93\n",
      "[1] 94\n",
      "[1] 95\n",
      "[1] 96\n",
      "[1] 97\n",
      "[1] 98\n",
      "[1] 99\n",
      "[1] 100\n",
      "[1] 101\n",
      "[1] 102\n",
      "[1] 103\n",
      "[1] 104\n",
      "[1] 105\n",
      "[1] 106\n",
      "[1] 107\n",
      "[1] 108\n",
      "[1] 109\n",
      "[1] 110\n",
      "[1] 111\n",
      "[1] 112\n",
      "[1] 113\n",
      "[1] 114\n",
      "[1] 115\n",
      "[1] 116\n",
      "[1] 117\n",
      "[1] 118\n",
      "[1] 119\n",
      "[1] 120\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>120</li>\n",
       "\t<li>8</li>\n",
       "\t<li>76606</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 120\n",
       "\\item 8\n",
       "\\item 76606\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 120\n",
       "2. 8\n",
       "3. 76606\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1]   120     8 76606"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>'samples'</li>\n",
       "\t<li>'founders'</li>\n",
       "\t<li>'markers'</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 'samples'\n",
       "\\item 'founders'\n",
       "\\item 'markers'\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 'samples'\n",
       "2. 'founders'\n",
       "3. 'markers'\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] \"samples\"  \"founders\" \"markers\" "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Create 3D probability array - Females\n",
    "model.probs = make_rix_model_probs(rix_list$Mating, rix_dir, 'F')\n",
    "dim(model.probs)\n",
    "names(dimnames(model.probs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>UNC6</th><th scope=col>JAX00000010</th><th scope=col>JAX00240603</th><th scope=col>JAX00240610</th><th scope=col>JAX00240613</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>A</th><td>3.00075e-06</td><td>2.62500e-06</td><td>2.25000e-06</td><td>1.87500e-06</td><td>1.50000e-06</td></tr>\n",
       "\t<tr><th scope=row>B</th><td>3.751e-06</td><td>3.375e-06</td><td>3.000e-06</td><td>2.625e-06</td><td>2.250e-06</td></tr>\n",
       "\t<tr><th scope=row>C</th><td>2.292e-06</td><td>1.875e-06</td><td>1.500e-06</td><td>1.125e-06</td><td>7.500e-07</td></tr>\n",
       "\t<tr><th scope=row>D</th><td>2.292e-06</td><td>1.875e-06</td><td>1.500e-06</td><td>1.125e-06</td><td>7.500e-07</td></tr>\n",
       "\t<tr><th scope=row>E</th><td>0.4999923</td><td>0.4999936</td><td>0.4999947</td><td>0.4999955</td><td>0.4999962</td></tr>\n",
       "\t<tr><th scope=row>F</th><td>0.4999948</td><td>0.4999959</td><td>0.4999970</td><td>0.4999978</td><td>0.4999985</td></tr>\n",
       "\t<tr><th scope=row>G</th><td>7.92e-07</td><td>3.75e-07</td><td>0.00e+00</td><td>0.00e+00</td><td>0.00e+00</td></tr>\n",
       "\t<tr><th scope=row>H</th><td>7.5075e-07</td><td>3.7500e-07</td><td>0.0000e+00</td><td>0.0000e+00</td><td>0.0000e+00</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lllll}\n",
       "  & UNC6 & JAX00000010 & JAX00240603 & JAX00240610 & JAX00240613\\\\\n",
       "\\hline\n",
       "\tA & 3.00075e-06 & 2.62500e-06 & 2.25000e-06 & 1.87500e-06 & 1.50000e-06\\\\\n",
       "\tB & 3.751e-06 & 3.375e-06 & 3.000e-06 & 2.625e-06 & 2.250e-06\\\\\n",
       "\tC & 2.292e-06 & 1.875e-06 & 1.500e-06 & 1.125e-06 & 7.500e-07\\\\\n",
       "\tD & 2.292e-06 & 1.875e-06 & 1.500e-06 & 1.125e-06 & 7.500e-07\\\\\n",
       "\tE & 0.4999923 & 0.4999936 & 0.4999947 & 0.4999955 & 0.4999962\\\\\n",
       "\tF & 0.4999948 & 0.4999959 & 0.4999970 & 0.4999978 & 0.4999985\\\\\n",
       "\tG & 7.92e-07 & 3.75e-07 & 0.00e+00 & 0.00e+00 & 0.00e+00\\\\\n",
       "\tH & 7.5075e-07 & 3.7500e-07 & 0.0000e+00 & 0.0000e+00 & 0.0000e+00\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "1. 3.00075e-06\n",
       "2. 3.751e-06\n",
       "3. 2.292e-06\n",
       "4. 2.292e-06\n",
       "5. 0.4999923325\n",
       "6. 0.499994789\n",
       "7. 7.92e-07\n",
       "8. 7.5075e-07\n",
       "9. 2.625e-06\n",
       "10. 3.375e-06\n",
       "11. 1.875e-06\n",
       "12. 1.875e-06\n",
       "13. 0.499993625\n",
       "14. 0.499995875\n",
       "15. 3.75e-07\n",
       "16. 3.75e-07\n",
       "17. 2.25e-06\n",
       "18. 3e-06\n",
       "19. 1.5e-06\n",
       "20. 1.5e-06\n",
       "21. 0.49999475\n",
       "22. 0.499997\n",
       "23. 0\n",
       "24. 0\n",
       "25. 1.875e-06\n",
       "26. 2.625e-06\n",
       "27. 1.125e-06\n",
       "28. 1.125e-06\n",
       "29. 0.4999955\n",
       "30. 0.49999775\n",
       "31. 0\n",
       "32. 0\n",
       "33. 1.5e-06\n",
       "34. 2.25e-06\n",
       "35. 7.5e-07\n",
       "36. 7.5e-07\n",
       "37. 0.49999625\n",
       "38. 0.4999985\n",
       "39. 0\n",
       "40. 0\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "        markers\n",
       "founders         UNC6 JAX00000010 JAX00240603 JAX00240610 JAX00240613\n",
       "       A 3.000750e-06 0.000002625  0.00000225 0.000001875  0.00000150\n",
       "       B 3.751000e-06 0.000003375  0.00000300 0.000002625  0.00000225\n",
       "       C 2.292000e-06 0.000001875  0.00000150 0.000001125  0.00000075\n",
       "       D 2.292000e-06 0.000001875  0.00000150 0.000001125  0.00000075\n",
       "       E 4.999923e-01 0.499993625  0.49999475 0.499995500  0.49999625\n",
       "       F 4.999948e-01 0.499995875  0.49999700 0.499997750  0.49999850\n",
       "       G 7.920000e-07 0.000000375  0.00000000 0.000000000  0.00000000\n",
       "       H 7.507500e-07 0.000000375  0.00000000 0.000000000  0.00000000"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Check model.probs object\n",
    "model.probs[1,,1:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>120</li>\n",
       "\t<li>8</li>\n",
       "\t<li>76600</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 120\n",
       "\\item 8\n",
       "\\item 76600\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 120\n",
       "2. 8\n",
       "3. 76600\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1]   120     8 76600"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Remove markers with all zeros\n",
    "model.probs = model.probs[, , setdiff(all_markers, markers_zero)]\n",
    "dim(model.probs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>UNC6</th><th scope=col>JAX00000010</th><th scope=col>JAX00240603</th><th scope=col>JAX00240610</th><th scope=col>JAX00240613</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>A</th><td>3.00075e-06</td><td>2.62500e-06</td><td>2.25000e-06</td><td>1.87500e-06</td><td>1.50000e-06</td></tr>\n",
       "\t<tr><th scope=row>B</th><td>3.751e-06</td><td>3.375e-06</td><td>3.000e-06</td><td>2.625e-06</td><td>2.250e-06</td></tr>\n",
       "\t<tr><th scope=row>C</th><td>2.292e-06</td><td>1.875e-06</td><td>1.500e-06</td><td>1.125e-06</td><td>7.500e-07</td></tr>\n",
       "\t<tr><th scope=row>D</th><td>2.292e-06</td><td>1.875e-06</td><td>1.500e-06</td><td>1.125e-06</td><td>7.500e-07</td></tr>\n",
       "\t<tr><th scope=row>E</th><td>0.4999923</td><td>0.4999936</td><td>0.4999947</td><td>0.4999955</td><td>0.4999962</td></tr>\n",
       "\t<tr><th scope=row>F</th><td>0.4999948</td><td>0.4999959</td><td>0.4999970</td><td>0.4999978</td><td>0.4999985</td></tr>\n",
       "\t<tr><th scope=row>G</th><td>7.92e-07</td><td>3.75e-07</td><td>1.00e-20</td><td>1.00e-20</td><td>1.00e-20</td></tr>\n",
       "\t<tr><th scope=row>H</th><td>7.5075e-07</td><td>3.7500e-07</td><td>1.0000e-20</td><td>1.0000e-20</td><td>1.0000e-20</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lllll}\n",
       "  & UNC6 & JAX00000010 & JAX00240603 & JAX00240610 & JAX00240613\\\\\n",
       "\\hline\n",
       "\tA & 3.00075e-06 & 2.62500e-06 & 2.25000e-06 & 1.87500e-06 & 1.50000e-06\\\\\n",
       "\tB & 3.751e-06 & 3.375e-06 & 3.000e-06 & 2.625e-06 & 2.250e-06\\\\\n",
       "\tC & 2.292e-06 & 1.875e-06 & 1.500e-06 & 1.125e-06 & 7.500e-07\\\\\n",
       "\tD & 2.292e-06 & 1.875e-06 & 1.500e-06 & 1.125e-06 & 7.500e-07\\\\\n",
       "\tE & 0.4999923 & 0.4999936 & 0.4999947 & 0.4999955 & 0.4999962\\\\\n",
       "\tF & 0.4999948 & 0.4999959 & 0.4999970 & 0.4999978 & 0.4999985\\\\\n",
       "\tG & 7.92e-07 & 3.75e-07 & 1.00e-20 & 1.00e-20 & 1.00e-20\\\\\n",
       "\tH & 7.5075e-07 & 3.7500e-07 & 1.0000e-20 & 1.0000e-20 & 1.0000e-20\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "1. 3.00075e-06\n",
       "2. 3.751e-06\n",
       "3. 2.292e-06\n",
       "4. 2.292e-06\n",
       "5. 0.4999923325\n",
       "6. 0.499994789\n",
       "7. 7.92e-07\n",
       "8. 7.5075e-07\n",
       "9. 2.625e-06\n",
       "10. 3.375e-06\n",
       "11. 1.875e-06\n",
       "12. 1.875e-06\n",
       "13. 0.499993625\n",
       "14. 0.499995875\n",
       "15. 3.75e-07\n",
       "16. 3.75e-07\n",
       "17. 2.25e-06\n",
       "18. 3e-06\n",
       "19. 1.5e-06\n",
       "20. 1.5e-06\n",
       "21. 0.49999475\n",
       "22. 0.499997\n",
       "23. 1e-20\n",
       "24. 1e-20\n",
       "25. 1.875e-06\n",
       "26. 2.625e-06\n",
       "27. 1.125e-06\n",
       "28. 1.125e-06\n",
       "29. 0.4999955\n",
       "30. 0.49999775\n",
       "31. 1e-20\n",
       "32. 1e-20\n",
       "33. 1.5e-06\n",
       "34. 2.25e-06\n",
       "35. 7.5e-07\n",
       "36. 7.5e-07\n",
       "37. 0.49999625\n",
       "38. 0.4999985\n",
       "39. 1e-20\n",
       "40. 1e-20\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "        markers\n",
       "founders         UNC6 JAX00000010  JAX00240603  JAX00240610  JAX00240613\n",
       "       A 3.000750e-06 0.000002625 2.250000e-06 1.875000e-06 1.500000e-06\n",
       "       B 3.751000e-06 0.000003375 3.000000e-06 2.625000e-06 2.250000e-06\n",
       "       C 2.292000e-06 0.000001875 1.500000e-06 1.125000e-06 7.500000e-07\n",
       "       D 2.292000e-06 0.000001875 1.500000e-06 1.125000e-06 7.500000e-07\n",
       "       E 4.999923e-01 0.499993625 4.999947e-01 4.999955e-01 4.999962e-01\n",
       "       F 4.999948e-01 0.499995875 4.999970e-01 4.999978e-01 4.999985e-01\n",
       "       G 7.920000e-07 0.000000375 1.000000e-20 1.000000e-20 1.000000e-20\n",
       "       H 7.507500e-07 0.000000375 1.000000e-20 1.000000e-20 1.000000e-20"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Make all probabilities non-zero\n",
    "model.probs[model.probs == 0] = 1e-20\n",
    "model.probs[1,,1:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Save genome cache - Females\n",
    "save(model.probs, file=file.path(mapping_dir, 'rix_universal_model_prob_females_27-Jun-2016.rda'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Last Updated: 31-Aug-2016"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.3.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}