work/1_panel_design/20180716_primer_design_format.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Desired output format:\n",
    "- target and margins - ID Chr Start Stop\n",
    "- confounding SNPs in flanks - ID Chr\tStart\tStop\tref\talt\tSource\ttype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas\n",
    "import pybedtools\n",
    "from collections import defaultdict\n",
    "from Bio import AlignIO\n",
    "# for MSA to variant conversion example\n",
    "from Bio.Alphabet import generic_dna\n",
    "from Bio.Seq import Seq\n",
    "from Bio.SeqRecord import SeqRecord\n",
    "from Bio.Align import MultipleSeqAlignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# input\n",
    "AMPL_FILE = \"data/20180706_comb_ampl.maf\"\n",
    "PRIORITY_FILE = \"data/20180710_priority_ampl_data.xlsx\"\n",
    "PRIORITY_TYPES = ('X', 'Var', 'Cons', 'Other') # sheet names in priority file\n",
    "GENOME_FA = '../../../data/genome/AgamP3/Anopheles-gambiae-PEST_CHROMOSOMES_AgamP3.fa'\n",
    "# output\n",
    "INTERVAL_FILE = \"data/20180716_ampl_intervals.xlsx\"\n",
    "ALL_SNP_FILE = \"data/20180716_ampl_all_snps.xlsx\"\n",
    "HF_SNP_FILE = \"data/20180716_ampl_highfreq_snps.xlsx\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "591"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# alignments to list\n",
    "ampl_seq = list(AlignIO.parse(AMPL_FILE, \"maf\"))\n",
    "len(ampl_seq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chr</th>\n",
       "      <th>start</th>\n",
       "      <th>end</th>\n",
       "      <th>seqid</th>\n",
       "      <th>num_aligned_species</th>\n",
       "      <th>unaligned_species</th>\n",
       "      <th>aligned_len</th>\n",
       "      <th>len_agam</th>\n",
       "      <th>total_indels</th>\n",
       "      <th>total_snvs</th>\n",
       "      <th>...</th>\n",
       "      <th>CpipJ2_identity</th>\n",
       "      <th>CpipJ2_q_span</th>\n",
       "      <th>CpipJ2_s_span</th>\n",
       "      <th>CpipJ2_expect</th>\n",
       "      <th>BDGP6_hits</th>\n",
       "      <th>BDGP6_length</th>\n",
       "      <th>BDGP6_identity</th>\n",
       "      <th>BDGP6_q_span</th>\n",
       "      <th>BDGP6_s_span</th>\n",
       "      <th>BDGP6_expect</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>575</th>\n",
       "      <td>X</td>\n",
       "      <td>876693</td>\n",
       "      <td>876909</td>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>20</td>\n",
       "      <td>AgamM1</td>\n",
       "      <td>277</td>\n",
       "      <td>216</td>\n",
       "      <td>122</td>\n",
       "      <td>32</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>576</th>\n",
       "      <td>X</td>\n",
       "      <td>961289</td>\n",
       "      <td>961516</td>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>21</td>\n",
       "      <td>NaN</td>\n",
       "      <td>231</td>\n",
       "      <td>227</td>\n",
       "      <td>7</td>\n",
       "      <td>42</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>577</th>\n",
       "      <td>X</td>\n",
       "      <td>1026174</td>\n",
       "      <td>1026456</td>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>21</td>\n",
       "      <td>NaN</td>\n",
       "      <td>282</td>\n",
       "      <td>282</td>\n",
       "      <td>0</td>\n",
       "      <td>55</td>\n",
       "      <td>...</td>\n",
       "      <td>0.900709</td>\n",
       "      <td>1-282</td>\n",
       "      <td>supercont3.123:864150-864431</td>\n",
       "      <td>2.274390e-99</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>579</th>\n",
       "      <td>X</td>\n",
       "      <td>1815894</td>\n",
       "      <td>1816086</td>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>18</td>\n",
       "      <td>AgamM1;AsinS1;AdarC2</td>\n",
       "      <td>192</td>\n",
       "      <td>192</td>\n",
       "      <td>5</td>\n",
       "      <td>28</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>582</th>\n",
       "      <td>X</td>\n",
       "      <td>10000121</td>\n",
       "      <td>10000346</td>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>20</td>\n",
       "      <td>AgamM1</td>\n",
       "      <td>225</td>\n",
       "      <td>225</td>\n",
       "      <td>0</td>\n",
       "      <td>32</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>587</th>\n",
       "      <td>X</td>\n",
       "      <td>18742149</td>\n",
       "      <td>18742401</td>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>20</td>\n",
       "      <td>AchrA1</td>\n",
       "      <td>276</td>\n",
       "      <td>252</td>\n",
       "      <td>24</td>\n",
       "      <td>47</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>6 rows × 43 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "    chr     start       end        seqid  num_aligned_species  \\\n",
       "575   X    876693    876909  AgamP3.chrX                   20   \n",
       "576   X    961289    961516  AgamP3.chrX                   21   \n",
       "577   X   1026174   1026456  AgamP3.chrX                   21   \n",
       "579   X   1815894   1816086  AgamP3.chrX                   18   \n",
       "582   X  10000121  10000346  AgamP3.chrX                   20   \n",
       "587   X  18742149  18742401  AgamP3.chrX                   20   \n",
       "\n",
       "        unaligned_species  aligned_len  len_agam  total_indels  total_snvs  \\\n",
       "575                AgamM1          277       216           122          32   \n",
       "576                   NaN          231       227             7          42   \n",
       "577                   NaN          282       282             0          55   \n",
       "579  AgamM1;AsinS1;AdarC2          192       192             5          28   \n",
       "582                AgamM1          225       225             0          32   \n",
       "587                AchrA1          276       252            24          47   \n",
       "\n",
       "        ...       CpipJ2_identity  CpipJ2_q_span  \\\n",
       "575     ...                   NaN            NaN   \n",
       "576     ...                   NaN            NaN   \n",
       "577     ...              0.900709          1-282   \n",
       "579     ...                   NaN            NaN   \n",
       "582     ...                   NaN            NaN   \n",
       "587     ...                   NaN            NaN   \n",
       "\n",
       "                    CpipJ2_s_span  CpipJ2_expect  BDGP6_hits  BDGP6_length  \\\n",
       "575                           NaN            NaN         NaN           NaN   \n",
       "576                           NaN            NaN         NaN           NaN   \n",
       "577  supercont3.123:864150-864431   2.274390e-99         NaN           NaN   \n",
       "579                           NaN            NaN         NaN           NaN   \n",
       "582                           NaN            NaN         NaN           NaN   \n",
       "587                           NaN            NaN         NaN           NaN   \n",
       "\n",
       "     BDGP6_identity  BDGP6_q_span BDGP6_s_span BDGP6_expect  \n",
       "575             NaN           NaN          NaN          NaN  \n",
       "576             NaN           NaN          NaN          NaN  \n",
       "577             NaN           NaN          NaN          NaN  \n",
       "579             NaN           NaN          NaN          NaN  \n",
       "582             NaN           NaN          NaN          NaN  \n",
       "587             NaN           NaN          NaN          NaN  \n",
       "\n",
       "[6 rows x 43 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# metadata to dict of dataframes\n",
    "ampl_data = dict()\n",
    "for p in PRIORITY_TYPES:\n",
    "    ampl_data[p] = pandas.read_excel(PRIORITY_FILE, sheet_name=p)\n",
    "ampl_data['X']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## File with intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chr</th>\n",
       "      <th>start</th>\n",
       "      <th>end</th>\n",
       "      <th>target_start</th>\n",
       "      <th>target_end</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>575</th>\n",
       "      <td>X</td>\n",
       "      <td>876693</td>\n",
       "      <td>876909</td>\n",
       "      <td>876749</td>\n",
       "      <td>876867</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>576</th>\n",
       "      <td>X</td>\n",
       "      <td>961289</td>\n",
       "      <td>961516</td>\n",
       "      <td>961335</td>\n",
       "      <td>961447</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>577</th>\n",
       "      <td>X</td>\n",
       "      <td>1026174</td>\n",
       "      <td>1026456</td>\n",
       "      <td>1026221</td>\n",
       "      <td>1026414</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>579</th>\n",
       "      <td>X</td>\n",
       "      <td>1815894</td>\n",
       "      <td>1816086</td>\n",
       "      <td>1815984</td>\n",
       "      <td>1816045</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>582</th>\n",
       "      <td>X</td>\n",
       "      <td>10000121</td>\n",
       "      <td>10000346</td>\n",
       "      <td>10000169</td>\n",
       "      <td>10000304</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>587</th>\n",
       "      <td>X</td>\n",
       "      <td>18742149</td>\n",
       "      <td>18742401</td>\n",
       "      <td>18742234</td>\n",
       "      <td>18742321</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    chr     start       end  target_start  target_end\n",
       "575   X    876693    876909        876749      876867\n",
       "576   X    961289    961516        961335      961447\n",
       "577   X   1026174   1026456       1026221     1026414\n",
       "579   X   1815894   1816086       1815984     1816045\n",
       "582   X  10000121  10000346      10000169    10000304\n",
       "587   X  18742149  18742401      18742234    18742321"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# example of interval data\n",
    "ampl_data['X'][['chr', 'start', 'end', 'target_start', 'target_end']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "writer = pandas.ExcelWriter(INTERVAL_FILE)\n",
    "for p in PRIORITY_TYPES:\n",
    "    ampl_data[p][['chr', 'start', 'end', 'target_start', 'target_end']].to_excel(writer, p)\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parse variants from multiple sequence alignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check genomic coordinates\n",
    "def coordinate_check(alignment, meta, genome):\n",
    "    '''\n",
    "    Given alignment, series of metadata (chr, start, end used) and corresponding genome,\n",
    "    assert that alignment and metadata have the same underlying sequence\n",
    "    '''\n",
    "    bed_str = '\\t'.join([str(meta['chr']), str(meta['start']), str(meta['end'])])\n",
    "    bed = pybedtools.BedTool(bed_str, from_string=True)\n",
    "    a = bed.sequence(fi=genome)\n",
    "    genome_seq = open(a.seqfn).read().split('\\n')[1]\n",
    "    aln_seq = str(alignment[0].seq).replace(\"-\",\"\")\n",
    "    pybedtools.cleanup(remove_all=True)\n",
    "    if genome_seq != aln_seq:\n",
    "        raise ValueError('Unequal sequences:\\n{}\\n{}\\n{}'.format(genome_seq, aln_seq, alignment[0].seq))\n",
    "    return 0\n",
    "    \n",
    "coordinate_check(ampl_seq[575], ampl_data['X'].loc[575], GENOME_FA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DNAAlphabet() alignment with 3 rows and 7 columns\n",
      "-AAC-T- Alpha.chr1\n",
      "CC-CGTA Beta.chr2\n",
      "GTAGGTA Gamma.chr3\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chr</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt_alleles</th>\n",
       "      <th>alt_genomes</th>\n",
       "      <th>num_genomes</th>\n",
       "      <th>%_alt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Alpha.chr1</td>\n",
       "      <td>100</td>\n",
       "      <td>A</td>\n",
       "      <td>T;C</td>\n",
       "      <td>Gamma;Beta</td>\n",
       "      <td>3</td>\n",
       "      <td>0.666667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Alpha.chr1</td>\n",
       "      <td>101</td>\n",
       "      <td>A</td>\n",
       "      <td>-</td>\n",
       "      <td>Beta</td>\n",
       "      <td>3</td>\n",
       "      <td>0.333333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Alpha.chr1</td>\n",
       "      <td>102</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "      <td>Gamma</td>\n",
       "      <td>3</td>\n",
       "      <td>0.333333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Alpha.chr1</td>\n",
       "      <td>102</td>\n",
       "      <td>C</td>\n",
       "      <td>CG;GG</td>\n",
       "      <td>Beta;Gamma</td>\n",
       "      <td>3</td>\n",
       "      <td>0.666667</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          chr  pos ref alt_alleles alt_genomes  num_genomes     %_alt\n",
       "0  Alpha.chr1  100   A         T;C  Gamma;Beta            3  0.666667\n",
       "1  Alpha.chr1  101   A           -        Beta            3  0.333333\n",
       "2  Alpha.chr1  102   C           G       Gamma            3  0.333333\n",
       "3  Alpha.chr1  102   C       CG;GG  Beta;Gamma            3  0.666667"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def msa_to_var(alignment, ref=0):\n",
    "    '''\n",
    "    Transform MultipleSeqAlignment into list of variants relative to specified reference. \n",
    "    Info for each variant:\n",
    "    chr pos ref alt alt_spp num_aligned alt_freq\n",
    "    \n",
    "    Caveats:\n",
    "    - deletions are per-base denoted as '-', while insertions are relative to previous reference base\n",
    "    - inserts at start and end of reference sequence are not displayed\n",
    "    - substitutions prior to insers are reported as separate variants\n",
    "    '''\n",
    "    \n",
    "    def var_data(alts, seqid, pos, ref_allele, num_genomes):\n",
    "        alt_alleles = ';'.join(alts.keys())\n",
    "        alt_genomes = ';'.join([','.join(alts[k]) for k in alts.keys()])\n",
    "        num_alt = sum((len(v) for v in alts.values()))\n",
    "        percent_alt = num_alt / num_genomes\n",
    "        return {\n",
    "            'chr': seqid,\n",
    "            'pos': pos,\n",
    "            'ref': ref_allele,\n",
    "            'alt_alleles': alt_alleles,\n",
    "            'alt_genomes': alt_genomes,\n",
    "            'num_genomes': num_genomes,\n",
    "            '%_alt': percent_alt\n",
    "        }\n",
    "            \n",
    "    variants = []\n",
    "    \n",
    "    seqid = alignment[ref].id # reference chromosome\n",
    "    pos = alignment[ref].annotations['start'] # position in reference genome\n",
    "    ref_allele = '' # current reference nucleotide\n",
    "    inserts = dict() # insert sequences for non-reference\n",
    "    genomes = [s.id.split('.')[0] for s in alignment] # genome names from seqids \n",
    "    # (expected format <genome>.<chromosome>)\n",
    "    num_genomes = len(alignment)\n",
    "    for i in range(alignment.get_alignment_length()):\n",
    "        # extract alignment column\n",
    "        col = alignment[:, i]\n",
    "        # insertion into reference - only collect insert data\n",
    "        if col[ref] == '-':\n",
    "            if ref_allele == '':\n",
    "                # beginning of reference sequence - skip\n",
    "                continue\n",
    "            # first base of insert - per-genome genotypes with previous base and insert (if any)\n",
    "            if len(inserts) == 0:\n",
    "                prev_col = alignment[:, i - 1]\n",
    "                for (j, gt) in enumerate(col):\n",
    "                    inserts[genomes[j]] = prev_col[j] + gt.replace(\"-\", \"\")\n",
    "            # further inserts - collect insert genotypes\n",
    "            else:\n",
    "                for (j, gt) in enumerate(col):\n",
    "                    inserts[genomes[j]] += gt.replace(\"-\", \"\")\n",
    "            continue\n",
    "            \n",
    "        # insert variants for previous reference base\n",
    "        if len(inserts) > 0:\n",
    "            # set reference allele\n",
    "            ref_insert = inserts[genomes[ref]]\n",
    "            assert len(ref_insert) == 1\n",
    "            # collect alt genotypes from insert dict\n",
    "            alt_inserts = defaultdict(list)\n",
    "            for (k, v) in inserts.items():\n",
    "                if v != ref_insert:\n",
    "                    alt_inserts[v].append(k)\n",
    "            # previous base \n",
    "            variants.append(var_data(alt_inserts, seqid, pos - 1, ref_insert, num_genomes))\n",
    "            inserts = dict()\n",
    "            \n",
    "        # set reference allele\n",
    "        ref_allele = col[ref]\n",
    "        # collect alt genotypes (subst & del)\n",
    "        alts = defaultdict(list)\n",
    "        for (j, gt) in enumerate(col):\n",
    "            if gt != ref_allele:\n",
    "                alts[gt].append(genomes[j])\n",
    "        # write alt genotypes (subst & del)\n",
    "        if len(alts) > 0:\n",
    "            # summary over alt alleles\n",
    "            variants.append(var_data(alts, seqid, pos, ref_allele, num_genomes))\n",
    "        pos += 1\n",
    "        \n",
    "    return pandas.DataFrame(variants, columns=['chr', 'pos', 'ref', 'alt_alleles', \n",
    "                                               'alt_genomes', 'num_genomes', '%_alt'])\n",
    "        \n",
    "# test example\n",
    "a = SeqRecord(Seq(\"-AAC-T-\", generic_dna), id=\"Alpha.chr1\", annotations={'start':100})\n",
    "b = SeqRecord(Seq(\"CC-CGTA\", generic_dna), id=\"Beta.chr2\", annotations={'start':234})\n",
    "c = SeqRecord(Seq(\"GTAGGTA\", generic_dna), id=\"Gamma.chr3\", annotations={'start':345})\n",
    "align = MultipleSeqAlignment([a, b, c])\n",
    "print(align)\n",
    "msa_to_var(align)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATGACCAGCGCAGACTTTCAGGCCATGATACCGTCTCACTTTATCAGCGGCGGCGGGCGGCAGGTAACC---GTGGACAAGTCTGACCACGGCCCGT----------------------------------------------------------CGGTAACGGTGACGGTGCAGAAGGCCGTGCCGGACATGCCGATGTTGCCGCCGGCACCGACCGAGCTCGGCACCGTCACGGAAACGATTACCAAGTCAACGTTCACTGAGACGGTAATGACC\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chr</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt_alleles</th>\n",
       "      <th>alt_genomes</th>\n",
       "      <th>num_genomes</th>\n",
       "      <th>%_alt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>50</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>1026422</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>AfunF1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.047619</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>51</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>1026425</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>AchrA1,AminM1,AculA1,AfunF1,AsteS1,AsteI2,Amac...</td>\n",
       "      <td>21</td>\n",
       "      <td>0.380952</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>52</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>1026431</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>AchrA1,AdarC2</td>\n",
       "      <td>21</td>\n",
       "      <td>0.095238</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>53</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>1026449</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>AculA1,AfunF1,AsteS1,AsteI2</td>\n",
       "      <td>21</td>\n",
       "      <td>0.190476</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>54</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>1026455</td>\n",
       "      <td>C</td>\n",
       "      <td>T;A</td>\n",
       "      <td>AculA1;AchrA1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.095238</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            chr      pos ref alt_alleles  \\\n",
       "50  AgamP3.chrX  1026422   C           T   \n",
       "51  AgamP3.chrX  1026425   G           A   \n",
       "52  AgamP3.chrX  1026431   C           T   \n",
       "53  AgamP3.chrX  1026449   G           A   \n",
       "54  AgamP3.chrX  1026455   C         T;A   \n",
       "\n",
       "                                          alt_genomes  num_genomes     %_alt  \n",
       "50                                             AfunF1           21  0.047619  \n",
       "51  AchrA1,AminM1,AculA1,AfunF1,AsteS1,AsteI2,Amac...           21  0.380952  \n",
       "52                                      AchrA1,AdarC2           21  0.095238  \n",
       "53                        AculA1,AfunF1,AsteS1,AsteI2           21  0.190476  \n",
       "54                                      AculA1;AchrA1           21  0.095238  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# real data examples\n",
    "print(ampl_seq[575][0].seq)\n",
    "msa_to_var(ampl_seq[577]).tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "961289\n",
      "T-GCAGGAAACTAAGCGGCCTTCAAAGGAGATGCAAGTTACTATCGCTCGTCAGCTGGGACTGGAGCCCACGACGGTCGGGAACTTTTTCATGAACGCACGGCGACGCTCGATGGACAAATGGAAGGACGAG---TCGATGAAGGGCCGGAGCAGCAGTAGTGACAATGTTATACTAGAGCACCTGAAGGACGATTCGTCCGACGAGAACCTGGGCTACGATAGTGCGTTC\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chr</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt_alleles</th>\n",
       "      <th>alt_genomes</th>\n",
       "      <th>num_genomes</th>\n",
       "      <th>%_alt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>961289</td>\n",
       "      <td>T</td>\n",
       "      <td>TC</td>\n",
       "      <td>AfarF1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.047619</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>961290</td>\n",
       "      <td>G</td>\n",
       "      <td>T;C;A</td>\n",
       "      <td>AalbS1;AmacM1;AchrA1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.142857</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>961291</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>AfarF1,AdarC2,AalbS1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.142857</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>961311</td>\n",
       "      <td>A</td>\n",
       "      <td>T;C</td>\n",
       "      <td>AquaS1;AdarC2,AalbS1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.142857</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>961317</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>AdarC2,AalbS1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.095238</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>961332</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>AculA1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.047619</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>961449</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>AdarC2,AalbS1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.095238</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>961452</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>AdarC2,AalbS1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.095238</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>41</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>961454</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>AdarC2,AalbS1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.095238</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>961479</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>AdirW1,AatrE1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.095238</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>43</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>961494</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>AchrA1,AsinS1,AdarC2,AalbS1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.190476</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>961509</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>AfarF1,AdirW1,AatrE1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.142857</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>45</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>961512</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>AfarF1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.047619</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>46</th>\n",
       "      <td>AgamP3.chrX</td>\n",
       "      <td>961515</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>AchrA1,AfarF1,AdirW1</td>\n",
       "      <td>21</td>\n",
       "      <td>0.142857</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            chr     pos ref alt_alleles                  alt_genomes  \\\n",
       "0   AgamP3.chrX  961289   T          TC                       AfarF1   \n",
       "1   AgamP3.chrX  961290   G       T;C;A         AalbS1;AmacM1;AchrA1   \n",
       "2   AgamP3.chrX  961291   C           T         AfarF1,AdarC2,AalbS1   \n",
       "3   AgamP3.chrX  961311   A         T;C         AquaS1;AdarC2,AalbS1   \n",
       "4   AgamP3.chrX  961317   G           A                AdarC2,AalbS1   \n",
       "5   AgamP3.chrX  961332   C           T                       AculA1   \n",
       "39  AgamP3.chrX  961449   C           A                AdarC2,AalbS1   \n",
       "40  AgamP3.chrX  961452   T           C                AdarC2,AalbS1   \n",
       "41  AgamP3.chrX  961454   T           C                AdarC2,AalbS1   \n",
       "42  AgamP3.chrX  961479   T           C                AdirW1,AatrE1   \n",
       "43  AgamP3.chrX  961494   C           T  AchrA1,AsinS1,AdarC2,AalbS1   \n",
       "44  AgamP3.chrX  961509   T           C         AfarF1,AdirW1,AatrE1   \n",
       "45  AgamP3.chrX  961512   G           A                       AfarF1   \n",
       "46  AgamP3.chrX  961515   C           T         AchrA1,AfarF1,AdirW1   \n",
       "\n",
       "    num_genomes     %_alt  \n",
       "0            21  0.047619  \n",
       "1            21  0.142857  \n",
       "2            21  0.142857  \n",
       "3            21  0.142857  \n",
       "4            21  0.095238  \n",
       "5            21  0.047619  \n",
       "39           21  0.095238  \n",
       "40           21  0.095238  \n",
       "41           21  0.095238  \n",
       "42           21  0.095238  \n",
       "43           21  0.190476  \n",
       "44           21  0.142857  \n",
       "45           21  0.047619  \n",
       "46           21  0.142857  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def confounding_vars(alignment, meta, genome):\n",
    "    '''\n",
    "    Given alignment, series of metadata (chr, start, end, target_start, target_end used) and corresponding genome,\n",
    "    yield confounding SNPs in flanking sequence.\n",
    "    '''\n",
    "    # does metadata have correct genomic coordinates for alignment? \n",
    "    coordinate_check(alignment, meta, genome)\n",
    "    # get variant\n",
    "    v = msa_to_var(alignment)\n",
    "    # filter flanking variants\n",
    "    return v[(v['pos'] < meta['target_start']) | (v['pos'] >= meta['target_end'])]\n",
    "\n",
    "print(ampl_seq[576][0].annotations['start'])\n",
    "print(ampl_seq[576][0].seq)\n",
    "\n",
    "confounding_vars(ampl_seq[576], ampl_data['X'].loc[576], GENOME_FA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Files with confounding variants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_writer = pandas.ExcelWriter(ALL_SNP_FILE)\n",
    "hf_writer = pandas.ExcelWriter(HF_SNP_FILE)\n",
    "for p in PRIORITY_TYPES:\n",
    "    ampl_vars = []\n",
    "    for i in ampl_data[p].index:\n",
    "        avd = confounding_vars(ampl_seq[i], ampl_data[p].loc[i], GENOME_FA)\n",
    "        avd.index = ['{}_{}'.format(i, j) for j in range(len(avd))]\n",
    "        avd['chr'] = avd['chr'].str.slice(10)\n",
    "        ampl_vars.append(avd)\n",
    "    av_pd = pandas.concat(ampl_vars, axis=0)\n",
    "    av_pd.to_excel(all_writer, p)\n",
    "    # over 1 genome wi\n",
    "    av_pd[(av_pd['alt_genomes'].str.count(';') > 0) |\n",
    "          (av_pd['alt_genomes'].str.count(',') > 0)].to_excel(hf_writer, p)\n",
    "all_writer.save()\n",
    "hf_writer.save()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}