run_chi_squared.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = 50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas\n",
    "import glob\n",
    "import stats\n",
    "import scipy\n",
    "import numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "train = pandas.read_csv('train.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "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>project_id</th>\n",
       "      <th>run_accession</th>\n",
       "      <th>name</th>\n",
       "      <th>flag</th>\n",
       "      <th>ref_name</th>\n",
       "      <th>ref_pos</th>\n",
       "      <th>map_quality</th>\n",
       "      <th>cigar</th>\n",
       "      <th>next_ref_name</th>\n",
       "      <th>next_ref_pos</th>\n",
       "      <th>length</th>\n",
       "      <th>seq</th>\n",
       "      <th>qual</th>\n",
       "      <th>tags</th>\n",
       "      <th>seq_len</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>PRJNA281708</td>\n",
       "      <td>SRR1982727</td>\n",
       "      <td>SRR1982727.183616</td>\n",
       "      <td>0</td>\n",
       "      <td>15</td>\n",
       "      <td>40555471</td>\n",
       "      <td>60</td>\n",
       "      <td>101M</td>\n",
       "      <td>*</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>AAGTAAATACCAAATTATCAGTGGGTCAGGAGTCTGCATCTTCAAA...</td>\n",
       "      <td>@CCDDFFFHHHHHIIJJJJJIFHCHEDFHJIJFGHIIJJIIIJJJJ...</td>\n",
       "      <td>['AS:i:0', 'XN:i:0', 'XM:i:0', 'XO:i:0', 'XG:i...</td>\n",
       "      <td>101</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    project_id run_accession               name  flag ref_name   ref_pos  \\\n",
       "0  PRJNA281708    SRR1982727  SRR1982727.183616     0       15  40555471   \n",
       "\n",
       "   map_quality cigar next_ref_name  next_ref_pos  length  \\\n",
       "0           60  101M             *             0       0   \n",
       "\n",
       "                                                 seq  \\\n",
       "0  AAGTAAATACCAAATTATCAGTGGGTCAGGAGTCTGCATCTTCAAA...   \n",
       "\n",
       "                                                qual  \\\n",
       "0  @CCDDFFFHHHHHIIJJJJJIFHCHEDFHJIJFGHIIJJIIIJJJJ...   \n",
       "\n",
       "                                                tags  seq_len  \n",
       "0  ['AS:i:0', 'XN:i:0', 'XM:i:0', 'XO:i:0', 'XG:i...      101  "
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train.head(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = pandas.read_csv('test.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clean_seq(sequences, K):\n",
    "    tmp = []\n",
    "    nucleotides = [x*K for x in ['A','C','G','T','N']]\n",
    "    for sequence in sequences:\n",
    "        if any([True if repeated_nucleotide in sequence else False for repeated_nucleotide in nucleotides ]):\n",
    "            continue\n",
    "        else:\n",
    "            tmp.append(sequence)   \n",
    "    return tmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_kmers(sequence, ksize):\n",
    "    kmers = []\n",
    "    n_kmers = len(sequence) - ksize + 1\n",
    "\n",
    "    for i in range(n_kmers):\n",
    "        kmer = sequence[i:i + ksize]\n",
    "        kmers.append(kmer)\n",
    "\n",
    "    return kmers\n",
    "    \n",
    "def read_kmers_from_list(array, ksize):\n",
    "  all_kmers = []\n",
    "  for sequence in array:\n",
    "    kmers = build_kmers(sequence, ksize)\n",
    "    all_kmers += kmers\n",
    "  return all_kmers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_location_data(project_ID):\n",
    "    import requests\n",
    "    import json\n",
    "    import pandas\n",
    "    url = 'https://www.ebi.ac.uk/ena/portal/api/filereport'\n",
    "    params = {'accession':\"{}\".format(project_ID),\n",
    "              'fields':\"center_name,instrument_platform,instrument_model,host\",\n",
    "            'format':'json',\n",
    "            'download':'false',\n",
    "            'result':'read_run'}\n",
    "    \n",
    "    headers = {\"User-Agent\": \"Mozilla/5.0 (Windows NT 10.0; Win64; x64)\"+\n",
    "             \" AppleWebKit/537.36 (KHTML, like Gecko)\"+\n",
    "             \" Chrome/70.0.3538.77 Safari/537.36\"}\n",
    "\n",
    "    print(\"Accessing: \", url)\n",
    "    r = requests.get(url,\n",
    "                  params = params,\n",
    "                  headers = headers\n",
    "                  ) \n",
    "    print(\"\\tMaking Request To: \", r.url) \n",
    "    my_json = json.loads(r.text)\n",
    "    df = pandas.json_normalize(json.loads(r.text))\n",
    "    df = df.assign(project_id=project_ID)\n",
    "    \n",
    "    # Temp Fix if center_name field is empty\n",
    "    df['center_name'] = df['center_name'].replace('', project_ID)\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['PRJNA277616', 'PRJNA281708', 'PRJNA634526']\n"
     ]
    }
   ],
   "source": [
    "#projects  = glob.glob(\"{}/*\".format(bam_files_path), recursive = True)\n",
    "#projects = [x.split(\"/\")[-1] for x in projects]\n",
    "projects = ['PRJNA277616','PRJNA281708', 'PRJNA634526']\n",
    "print(projects)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accessing:  https://www.ebi.ac.uk/ena/portal/api/filereport\n",
      "\tMaking Request To:  https://www.ebi.ac.uk/ena/portal/api/filereport?accession=PRJNA277616&fields=center_name%2Cinstrument_platform%2Cinstrument_model%2Chost&format=json&download=false&result=read_run\n",
      "Accessing:  https://www.ebi.ac.uk/ena/portal/api/filereport\n",
      "\tMaking Request To:  https://www.ebi.ac.uk/ena/portal/api/filereport?accession=PRJNA281708&fields=center_name%2Cinstrument_platform%2Cinstrument_model%2Chost&format=json&download=false&result=read_run\n",
      "Accessing:  https://www.ebi.ac.uk/ena/portal/api/filereport\n",
      "\tMaking Request To:  https://www.ebi.ac.uk/ena/portal/api/filereport?accession=PRJNA634526&fields=center_name%2Cinstrument_platform%2Cinstrument_model%2Chost&format=json&download=false&result=read_run\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>run_accession</th>\n",
       "      <th>sample_accession</th>\n",
       "      <th>center_name</th>\n",
       "      <th>instrument_platform</th>\n",
       "      <th>instrument_model</th>\n",
       "      <th>host</th>\n",
       "      <th>project_id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>SRR1867792</td>\n",
       "      <td>SAMN03394583</td>\n",
       "      <td>Weill Cornell Medical College in Qatar</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina HiSeq 2500</td>\n",
       "      <td></td>\n",
       "      <td>PRJNA277616</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>SRR1909613</td>\n",
       "      <td>SAMN03394592</td>\n",
       "      <td>Weill Cornell Medical College in Qatar</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina HiSeq 2500</td>\n",
       "      <td></td>\n",
       "      <td>PRJNA277616</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>SRR1909637</td>\n",
       "      <td>SAMN03394591</td>\n",
       "      <td>Weill Cornell Medical College in Qatar</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina HiSeq 2500</td>\n",
       "      <td></td>\n",
       "      <td>PRJNA277616</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>SRR1909638</td>\n",
       "      <td>SAMN03394590</td>\n",
       "      <td>Weill Cornell Medical College in Qatar</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina HiSeq 2500</td>\n",
       "      <td></td>\n",
       "      <td>PRJNA277616</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>SRR1909639</td>\n",
       "      <td>SAMN03394589</td>\n",
       "      <td>Weill Cornell Medical College in Qatar</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina HiSeq 2500</td>\n",
       "      <td></td>\n",
       "      <td>PRJNA277616</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>416</th>\n",
       "      <td>SRR14209704</td>\n",
       "      <td>SAMN18699854</td>\n",
       "      <td>SUB9454898</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina NovaSeq 6000</td>\n",
       "      <td>Homo sapiens</td>\n",
       "      <td>PRJNA634526</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>417</th>\n",
       "      <td>SRR14209705</td>\n",
       "      <td>SAMN18699853</td>\n",
       "      <td>SUB9454898</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina NovaSeq 6000</td>\n",
       "      <td>Homo sapiens</td>\n",
       "      <td>PRJNA634526</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>418</th>\n",
       "      <td>SRR14209706</td>\n",
       "      <td>SAMN18699852</td>\n",
       "      <td>SUB9454898</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina NovaSeq 6000</td>\n",
       "      <td>Homo sapiens</td>\n",
       "      <td>PRJNA634526</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>419</th>\n",
       "      <td>SRR14209707</td>\n",
       "      <td>SAMN18699851</td>\n",
       "      <td>SUB9454898</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina NovaSeq 6000</td>\n",
       "      <td>Homo sapiens</td>\n",
       "      <td>PRJNA634526</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>420</th>\n",
       "      <td>SRR14209708</td>\n",
       "      <td>SAMN18699850</td>\n",
       "      <td>SUB9454898</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina NovaSeq 6000</td>\n",
       "      <td>Homo sapiens</td>\n",
       "      <td>PRJNA634526</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>421 rows × 7 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "    run_accession sample_accession                             center_name  \\\n",
       "0      SRR1867792     SAMN03394583  Weill Cornell Medical College in Qatar   \n",
       "1      SRR1909613     SAMN03394592  Weill Cornell Medical College in Qatar   \n",
       "2      SRR1909637     SAMN03394591  Weill Cornell Medical College in Qatar   \n",
       "3      SRR1909638     SAMN03394590  Weill Cornell Medical College in Qatar   \n",
       "4      SRR1909639     SAMN03394589  Weill Cornell Medical College in Qatar   \n",
       "..            ...              ...                                     ...   \n",
       "416   SRR14209704     SAMN18699854                              SUB9454898   \n",
       "417   SRR14209705     SAMN18699853                              SUB9454898   \n",
       "418   SRR14209706     SAMN18699852                              SUB9454898   \n",
       "419   SRR14209707     SAMN18699851                              SUB9454898   \n",
       "420   SRR14209708     SAMN18699850                              SUB9454898   \n",
       "\n",
       "    instrument_platform       instrument_model          host   project_id  \n",
       "0              ILLUMINA    Illumina HiSeq 2500                PRJNA277616  \n",
       "1              ILLUMINA    Illumina HiSeq 2500                PRJNA277616  \n",
       "2              ILLUMINA    Illumina HiSeq 2500                PRJNA277616  \n",
       "3              ILLUMINA    Illumina HiSeq 2500                PRJNA277616  \n",
       "4              ILLUMINA    Illumina HiSeq 2500                PRJNA277616  \n",
       "..                  ...                    ...           ...          ...  \n",
       "416            ILLUMINA  Illumina NovaSeq 6000  Homo sapiens  PRJNA634526  \n",
       "417            ILLUMINA  Illumina NovaSeq 6000  Homo sapiens  PRJNA634526  \n",
       "418            ILLUMINA  Illumina NovaSeq 6000  Homo sapiens  PRJNA634526  \n",
       "419            ILLUMINA  Illumina NovaSeq 6000  Homo sapiens  PRJNA634526  \n",
       "420            ILLUMINA  Illumina NovaSeq 6000  Homo sapiens  PRJNA634526  \n",
       "\n",
       "[421 rows x 7 columns]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pandas.concat([get_location_data(x) for x in projects], ignore_index=True)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "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>project_id</th>\n",
       "      <th>run_accession</th>\n",
       "      <th>name</th>\n",
       "      <th>flag</th>\n",
       "      <th>ref_name</th>\n",
       "      <th>ref_pos</th>\n",
       "      <th>map_quality</th>\n",
       "      <th>cigar</th>\n",
       "      <th>next_ref_name</th>\n",
       "      <th>next_ref_pos</th>\n",
       "      <th>length</th>\n",
       "      <th>seq</th>\n",
       "      <th>qual</th>\n",
       "      <th>tags</th>\n",
       "      <th>seq_len</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>PRJNA281708</td>\n",
       "      <td>SRR1982727</td>\n",
       "      <td>SRR1982727.183616</td>\n",
       "      <td>0</td>\n",
       "      <td>15</td>\n",
       "      <td>40555471</td>\n",
       "      <td>60</td>\n",
       "      <td>101M</td>\n",
       "      <td>*</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>AAGTAAATACCAAATTATCAGTGGGTCAGGAGTCTGCATCTTCAAA...</td>\n",
       "      <td>@CCDDFFFHHHHHIIJJJJJIFHCHEDFHJIJFGHIIJJIIIJJJJ...</td>\n",
       "      <td>['AS:i:0', 'XN:i:0', 'XM:i:0', 'XO:i:0', 'XG:i...</td>\n",
       "      <td>101</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>PRJNA277616</td>\n",
       "      <td>SRR1909638</td>\n",
       "      <td>SRR1909638.13403218</td>\n",
       "      <td>133</td>\n",
       "      <td>MT</td>\n",
       "      <td>2920</td>\n",
       "      <td>0</td>\n",
       "      <td>*</td>\n",
       "      <td>=</td>\n",
       "      <td>2920</td>\n",
       "      <td>0</td>\n",
       "      <td>CTCGTGGAGCCATTCATACAGGTCCCTATGGAGCTTTAATTTATTA...</td>\n",
       "      <td>FFFFFBFBFFFFIIFFFBFFFBFFFFFFFFIBFFFFIIFFIIIIFF...</td>\n",
       "      <td>['YT:Z:UP']</td>\n",
       "      <td>98</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>PRJNA281708</td>\n",
       "      <td>SRR1982826</td>\n",
       "      <td>SRR1982826.8793344</td>\n",
       "      <td>0</td>\n",
       "      <td>MT</td>\n",
       "      <td>2225</td>\n",
       "      <td>60</td>\n",
       "      <td>101M</td>\n",
       "      <td>*</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>CTAAAAAATCCCAAACATATAACTGAACTCCTCACACCCAATTGGA...</td>\n",
       "      <td>CBBCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGG...</td>\n",
       "      <td>['AS:i:0', 'XN:i:0', 'XM:i:0', 'XO:i:0', 'XG:i...</td>\n",
       "      <td>101</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    project_id run_accession                 name  flag ref_name   ref_pos  \\\n",
       "0  PRJNA281708    SRR1982727    SRR1982727.183616     0       15  40555471   \n",
       "1  PRJNA277616    SRR1909638  SRR1909638.13403218   133       MT      2920   \n",
       "2  PRJNA281708    SRR1982826   SRR1982826.8793344     0       MT      2225   \n",
       "\n",
       "   map_quality cigar next_ref_name  next_ref_pos  length  \\\n",
       "0           60  101M             *             0       0   \n",
       "1            0     *             =          2920       0   \n",
       "2           60  101M             *             0       0   \n",
       "\n",
       "                                                 seq  \\\n",
       "0  AAGTAAATACCAAATTATCAGTGGGTCAGGAGTCTGCATCTTCAAA...   \n",
       "1  CTCGTGGAGCCATTCATACAGGTCCCTATGGAGCTTTAATTTATTA...   \n",
       "2  CTAAAAAATCCCAAACATATAACTGAACTCCTCACACCCAATTGGA...   \n",
       "\n",
       "                                                qual  \\\n",
       "0  @CCDDFFFHHHHHIIJJJJJIFHCHEDFHJIJFGHIIJJIIIJJJJ...   \n",
       "1  FFFFFBFBFFFFIIFFFBFFFBFFFFFFFFIBFFFFIIFFIIIIFF...   \n",
       "2  CBBCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGG...   \n",
       "\n",
       "                                                tags  seq_len  \n",
       "0  ['AS:i:0', 'XN:i:0', 'XM:i:0', 'XO:i:0', 'XG:i...      101  \n",
       "1                                        ['YT:Z:UP']       98  \n",
       "2  ['AS:i:0', 'XN:i:0', 'XM:i:0', 'XO:i:0', 'XG:i...      101  "
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "project_id   run_accession\n",
       "PRJNA277616  SRR1867792       3399\n",
       "             SRR1909613       3373\n",
       "             SRR1909637       3423\n",
       "             SRR1909638       3372\n",
       "             SRR1909639       3350\n",
       "             SRR1926132       3378\n",
       "             SRR1926133       3429\n",
       "             SRR1926134       3318\n",
       "             SRR1926135       3418\n",
       "             SRR1926136       3396\n",
       "PRJNA281708  SRR1982584       3464\n",
       "             SRR1982585       3460\n",
       "             SRR1982588       3379\n",
       "             SRR1982612       3391\n",
       "             SRR1982660       3402\n",
       "             SRR1982727       3337\n",
       "             SRR1982783       3330\n",
       "             SRR1982787       3423\n",
       "             SRR1982814       3359\n",
       "             SRR1982826       3481\n",
       "PRJNA634526  SRR14209589      3796\n",
       "             SRR14209591      3767\n",
       "             SRR14209592      3737\n",
       "             SRR14209603      3770\n",
       "             SRR14209614      3778\n",
       "             SRR14209634      3779\n",
       "             SRR14209651      3688\n",
       "             SRR14209655      3826\n",
       "             SRR14209708      3773\n",
       "dtype: int64"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train = train.set_index(['project_id', 'run_accession'])\n",
    "train.groupby(level=['project_id','run_accession']).size()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "train = train.sample(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "train = train.sort_index()\n",
    "indexes = sorted(list(set(train.index)))\n",
    "mydict = {}\n",
    "N = len(indexes)\n",
    "for i in range(N):\n",
    "    key = indexes[i]\n",
    "    seq = train.loc[key]['seq'].to_list()\n",
    "    seq = clean_seq(seq, K)\n",
    "    arr = read_kmers_from_list(seq, K)\n",
    "    mydict[key] = arr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = []\n",
    "for key in mydict.keys():\n",
    "    mydf = pandas.DataFrame()\n",
    "    mydf['kmer_seq'] = mydict[key]\n",
    "    mydf = mydf.assign(project_id=key[0])\n",
    "    mydf = mydf.assign(run_accession=key[1])\n",
    "    mydf = pandas.merge(mydf, df, on=['project_id','run_accession'], how='left')\n",
    "    tmp.append(mydf)\n",
    "kmer_df  = pandas.concat(tmp, ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "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>kmer_seq</th>\n",
       "      <th>project_id</th>\n",
       "      <th>run_accession</th>\n",
       "      <th>sample_accession</th>\n",
       "      <th>center_name</th>\n",
       "      <th>instrument_platform</th>\n",
       "      <th>instrument_model</th>\n",
       "      <th>host</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>GCATGATTTCCCATTCAGAGGCAGGTGCTGCCCTCATATCAGAAAA...</td>\n",
       "      <td>PRJNA277616</td>\n",
       "      <td>SRR1867792</td>\n",
       "      <td>SAMN03394583</td>\n",
       "      <td>Weill Cornell Medical College in Qatar</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina HiSeq 2500</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>CATGATTTCCCATTCAGAGGCAGGTGCTGCCCTCATATCAGAAAAG...</td>\n",
       "      <td>PRJNA277616</td>\n",
       "      <td>SRR1867792</td>\n",
       "      <td>SAMN03394583</td>\n",
       "      <td>Weill Cornell Medical College in Qatar</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina HiSeq 2500</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>ATGATTTCCCATTCAGAGGCAGGTGCTGCCCTCATATCAGAAAAGT...</td>\n",
       "      <td>PRJNA277616</td>\n",
       "      <td>SRR1867792</td>\n",
       "      <td>SAMN03394583</td>\n",
       "      <td>Weill Cornell Medical College in Qatar</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina HiSeq 2500</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>TGATTTCCCATTCAGAGGCAGGTGCTGCCCTCATATCAGAAAAGTA...</td>\n",
       "      <td>PRJNA277616</td>\n",
       "      <td>SRR1867792</td>\n",
       "      <td>SAMN03394583</td>\n",
       "      <td>Weill Cornell Medical College in Qatar</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina HiSeq 2500</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>GATTTCCCATTCAGAGGCAGGTGCTGCCCTCATATCAGAAAAGTAG...</td>\n",
       "      <td>PRJNA277616</td>\n",
       "      <td>SRR1867792</td>\n",
       "      <td>SAMN03394583</td>\n",
       "      <td>Weill Cornell Medical College in Qatar</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina HiSeq 2500</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>62193</th>\n",
       "      <td>CGTGTGTGTGTGTGTGTGTGTGGGTGTGTGTGTGTGTGTGTGTGTG...</td>\n",
       "      <td>PRJNA634526</td>\n",
       "      <td>SRR14209708</td>\n",
       "      <td>SAMN18699850</td>\n",
       "      <td>SUB9454898</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina NovaSeq 6000</td>\n",
       "      <td>Homo sapiens</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>62194</th>\n",
       "      <td>GTGTGTGTGTGTGTGTGTGTGGGTGTGTGTGTGTGTGTGTGTGTGT...</td>\n",
       "      <td>PRJNA634526</td>\n",
       "      <td>SRR14209708</td>\n",
       "      <td>SAMN18699850</td>\n",
       "      <td>SUB9454898</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina NovaSeq 6000</td>\n",
       "      <td>Homo sapiens</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>62195</th>\n",
       "      <td>TGTGTGTGTGTGTGTGTGTGGGTGTGTGTGTGTGTGTGTGTGTGTG...</td>\n",
       "      <td>PRJNA634526</td>\n",
       "      <td>SRR14209708</td>\n",
       "      <td>SAMN18699850</td>\n",
       "      <td>SUB9454898</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina NovaSeq 6000</td>\n",
       "      <td>Homo sapiens</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>62196</th>\n",
       "      <td>GTGTGTGTGTGTGTGTGTGGGTGTGTGTGTGTGTGTGTGTGTGTGT...</td>\n",
       "      <td>PRJNA634526</td>\n",
       "      <td>SRR14209708</td>\n",
       "      <td>SAMN18699850</td>\n",
       "      <td>SUB9454898</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina NovaSeq 6000</td>\n",
       "      <td>Homo sapiens</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>62197</th>\n",
       "      <td>TGTGTGTGTGTGTGTGTGGGTGTGTGTGTGTGTGTGTGTGTGTGTG...</td>\n",
       "      <td>PRJNA634526</td>\n",
       "      <td>SRR14209708</td>\n",
       "      <td>SAMN18699850</td>\n",
       "      <td>SUB9454898</td>\n",
       "      <td>ILLUMINA</td>\n",
       "      <td>Illumina NovaSeq 6000</td>\n",
       "      <td>Homo sapiens</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>62198 rows × 8 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                kmer_seq   project_id  \\\n",
       "0      GCATGATTTCCCATTCAGAGGCAGGTGCTGCCCTCATATCAGAAAA...  PRJNA277616   \n",
       "1      CATGATTTCCCATTCAGAGGCAGGTGCTGCCCTCATATCAGAAAAG...  PRJNA277616   \n",
       "2      ATGATTTCCCATTCAGAGGCAGGTGCTGCCCTCATATCAGAAAAGT...  PRJNA277616   \n",
       "3      TGATTTCCCATTCAGAGGCAGGTGCTGCCCTCATATCAGAAAAGTA...  PRJNA277616   \n",
       "4      GATTTCCCATTCAGAGGCAGGTGCTGCCCTCATATCAGAAAAGTAG...  PRJNA277616   \n",
       "...                                                  ...          ...   \n",
       "62193  CGTGTGTGTGTGTGTGTGTGTGGGTGTGTGTGTGTGTGTGTGTGTG...  PRJNA634526   \n",
       "62194  GTGTGTGTGTGTGTGTGTGTGGGTGTGTGTGTGTGTGTGTGTGTGT...  PRJNA634526   \n",
       "62195  TGTGTGTGTGTGTGTGTGTGGGTGTGTGTGTGTGTGTGTGTGTGTG...  PRJNA634526   \n",
       "62196  GTGTGTGTGTGTGTGTGTGGGTGTGTGTGTGTGTGTGTGTGTGTGT...  PRJNA634526   \n",
       "62197  TGTGTGTGTGTGTGTGTGGGTGTGTGTGTGTGTGTGTGTGTGTGTG...  PRJNA634526   \n",
       "\n",
       "      run_accession sample_accession                             center_name  \\\n",
       "0        SRR1867792     SAMN03394583  Weill Cornell Medical College in Qatar   \n",
       "1        SRR1867792     SAMN03394583  Weill Cornell Medical College in Qatar   \n",
       "2        SRR1867792     SAMN03394583  Weill Cornell Medical College in Qatar   \n",
       "3        SRR1867792     SAMN03394583  Weill Cornell Medical College in Qatar   \n",
       "4        SRR1867792     SAMN03394583  Weill Cornell Medical College in Qatar   \n",
       "...             ...              ...                                     ...   \n",
       "62193   SRR14209708     SAMN18699850                              SUB9454898   \n",
       "62194   SRR14209708     SAMN18699850                              SUB9454898   \n",
       "62195   SRR14209708     SAMN18699850                              SUB9454898   \n",
       "62196   SRR14209708     SAMN18699850                              SUB9454898   \n",
       "62197   SRR14209708     SAMN18699850                              SUB9454898   \n",
       "\n",
       "      instrument_platform       instrument_model          host  \n",
       "0                ILLUMINA    Illumina HiSeq 2500                \n",
       "1                ILLUMINA    Illumina HiSeq 2500                \n",
       "2                ILLUMINA    Illumina HiSeq 2500                \n",
       "3                ILLUMINA    Illumina HiSeq 2500                \n",
       "4                ILLUMINA    Illumina HiSeq 2500                \n",
       "...                   ...                    ...           ...  \n",
       "62193            ILLUMINA  Illumina NovaSeq 6000  Homo sapiens  \n",
       "62194            ILLUMINA  Illumina NovaSeq 6000  Homo sapiens  \n",
       "62195            ILLUMINA  Illumina NovaSeq 6000  Homo sapiens  \n",
       "62196            ILLUMINA  Illumina NovaSeq 6000  Homo sapiens  \n",
       "62197            ILLUMINA  Illumina NovaSeq 6000  Homo sapiens  \n",
       "\n",
       "[62198 rows x 8 columns]"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kmer_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "crosstab = pandas.crosstab(kmer_df.kmer_seq, kmer_df.project_id,margins = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "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>project_id</th>\n",
       "      <th>PRJNA277616</th>\n",
       "      <th>PRJNA281708</th>\n",
       "      <th>PRJNA634526</th>\n",
       "      <th>All</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>kmer_seq</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>AAAAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAA</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AAAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAAC</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACA</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACAA</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACAAA</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTTTT</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTTT</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTT</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTT</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>All</th>\n",
       "      <td>17249</td>\n",
       "      <td>16763</td>\n",
       "      <td>28186</td>\n",
       "      <td>62198</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>54187 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "project_id                                          PRJNA277616  PRJNA281708  \\\n",
       "kmer_seq                                                                       \n",
       "AAAAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAA            1            0   \n",
       "AAAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAAC            1            0   \n",
       "AAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACA            1            0   \n",
       "AAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACAA            1            0   \n",
       "AAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACAAA            1            0   \n",
       "...                                                         ...          ...   \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTTTT            1            0   \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTTT            1            0   \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTT            1            0   \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTT            1            0   \n",
       "All                                                       17249        16763   \n",
       "\n",
       "project_id                                          PRJNA634526    All  \n",
       "kmer_seq                                                                \n",
       "AAAAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAA            0      1  \n",
       "AAAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAAC            0      1  \n",
       "AAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACA            0      1  \n",
       "AAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACAA            0      1  \n",
       "AAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACAAA            0      1  \n",
       "...                                                         ...    ...  \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTTTT            0      1  \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTTT            0      1  \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTT            0      1  \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTT            0      1  \n",
       "All                                                       28186  62198  \n",
       "\n",
       "[54187 rows x 4 columns]"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "crosstab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cross_chi_table(crosstab, row, total, i):\n",
    "\n",
    "    total = total.iloc[0]\n",
    "    matrix = numpy.array([[row[0],total[0]-row[0]],\n",
    "                          [row[1],total[1]-row[1]]])\n",
    "    try:\n",
    "        chi2, p, dof, expected = chi2_contingency(matrix)\n",
    "        crosstab.loc[row.name,\"x_y\"] = p\n",
    "    except:\n",
    "        crosstab.loc[row.name,\"x_y\"] = numpy.NAN\n",
    "        \n",
    "    matrix = numpy.array([[row[1],total[1]-row[1]],\n",
    "                          [row[2],total[2]-row[2]]])\n",
    "    try:\n",
    "        chi2, p, dof, expected = chi2_contingency(matrix)\n",
    "        crosstab.loc[row.name,\"y_z\"] = p\n",
    "    except:\n",
    "        crosstab.loc[row.name,\"y_z\"] = numpy.NAN\n",
    "\n",
    "    matrix = numpy.array([[row[0],total[0]-row[0]],\n",
    "                          [row[2],total[2]-row[2]]])\n",
    "    try:\n",
    "        chi2, p, dof, expected = chi2_contingency(matrix)\n",
    "        crosstab.loc[row.name,\"x_z\"] = p\n",
    "    except:\n",
    "        crosstab.loc[row.name,\"x_z\"] = numpy.NAN\n",
    "    \n",
    "    return crosstab\n",
    "crosstab[\"x_y\"] = 0  \n",
    "crosstab[\"y_z\"] = 0\n",
    "crosstab[\"x_z\"] = 0\n",
    "for i in range(crosstab.shape[0]):\n",
    "    crosstab = cross_chi_table(crosstab,crosstab.iloc[i],crosstab.loc[[\"All\"]],i) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "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>project_id</th>\n",
       "      <th>PRJNA277616</th>\n",
       "      <th>PRJNA281708</th>\n",
       "      <th>PRJNA634526</th>\n",
       "      <th>All</th>\n",
       "      <th>x_y</th>\n",
       "      <th>y_z</th>\n",
       "      <th>x_z</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>kmer_seq</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>AAAAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAA</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AAAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAAC</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACA</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACAA</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACAAA</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTTTT</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTTT</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTT</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTT</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>All</th>\n",
       "      <td>17249</td>\n",
       "      <td>16763</td>\n",
       "      <td>28186</td>\n",
       "      <td>62198</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>54187 rows × 7 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "project_id                                          PRJNA277616  PRJNA281708  \\\n",
       "kmer_seq                                                                       \n",
       "AAAAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAA            1            0   \n",
       "AAAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAAC            1            0   \n",
       "AAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACA            1            0   \n",
       "AAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACAA            1            0   \n",
       "AAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACAAA            1            0   \n",
       "...                                                         ...          ...   \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTTTT            1            0   \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTTT            1            0   \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTT            1            0   \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTT            1            0   \n",
       "All                                                       17249        16763   \n",
       "\n",
       "project_id                                          PRJNA634526    All  x_y  \\\n",
       "kmer_seq                                                                      \n",
       "AAAAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAA            0      1  NaN   \n",
       "AAAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAAC            0      1  NaN   \n",
       "AAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACA            0      1  NaN   \n",
       "AAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACAA            0      1  NaN   \n",
       "AAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACAAA            0      1  NaN   \n",
       "...                                                         ...    ...  ...   \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTTTT            0      1  NaN   \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTTT            0      1  NaN   \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTT            0      1  NaN   \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTT            0      1  NaN   \n",
       "All                                                       28186  62198  NaN   \n",
       "\n",
       "project_id                                          y_z  x_z  \n",
       "kmer_seq                                                      \n",
       "AAAAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAA  NaN  NaN  \n",
       "AAAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAAC  NaN  NaN  \n",
       "AAAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACA  NaN  NaN  \n",
       "AAAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACAA  NaN  NaN  \n",
       "AAAAAAAAAAAAAAAAAAAAAACCCAAAAAAAAAACNCTCTAAGAACAAA  NaN  NaN  \n",
       "...                                                 ...  ...  \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTTTT  NaN  NaN  \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTTT  NaN  NaN  \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTTT  NaN  NaN  \n",
       "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCCCCACTTTT  NaN  NaN  \n",
       "All                                                 NaN  NaN  \n",
       "\n",
       "[54187 rows x 7 columns]"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "crosstab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_val_df = crosstab[(crosstab.x_y < 0.05) | (crosstab.y_z < 0.05) | (crosstab.x_z < 0.05)].\\\n",
    "    sort_values(by=['All'], ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "library = p_val_df.index.to_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"K_{}_library.txt\".format(K), 'w') as file:\n",
    "    file.write('# K={}\\n'.format(K))\n",
    "    file.write('\\n'.join(library))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "ename": "SystemExit",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "An exception has occurred, use %tb to see the full traceback.\n",
      "\u001b[0;31mSystemExit\u001b[0m\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/apps/software/standard/core/anaconda/2020.11-py3.8/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3426: UserWarning: To exit: use 'exit', 'quit', or Ctrl-D.\n",
      "  warn(\"To exit: use 'exit', 'quit', or Ctrl-D.\", stacklevel=1)\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "sys.exit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#https://pythonfordatascienceorg.wordpress.com/chi-square-python/\n",
    "#https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2_contingency.html?highlight=chi2_contingency#scipy.stats.chi2_contingency\n",
    "from scipy.stats import chi2_contingency\n",
    "chi2_contingency(crosstab)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chi2, p, dof, ex = chi2_contingency(crosstab)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"K: \",K)\n",
    "print(\"Chi-Square: \", chi2)\n",
    "print(\"p: \", p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "K = 30\n",
    "(1547445.6979213262,\n",
    " 0.0,\n",
    " 1138764,\n",
    " array([[0.28129705, 0.28784497, 0.43085798],\n",
    "        [0.28129705, 0.28784497, 0.43085798],\n",
    "        [0.5625941 , 0.57568993, 0.86171597],\n",
    "        ...,\n",
    "        [0.28129705, 0.28784497, 0.43085798],\n",
    "        [1.40648525, 1.43922483, 2.15428992],\n",
    "        [0.5625941 , 0.57568993, 0.86171597]]))\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "K = 55\n",
    "(1113199.1051552424,\n",
    " 0.0,\n",
    " 911602,\n",
    " array([[0.25997437, 0.26203298, 0.47799266],\n",
    "        [0.25997437, 0.26203298, 0.47799266],\n",
    "        [0.25997437, 0.26203298, 0.47799266],\n",
    "        ...,\n",
    "        [0.51994873, 0.52406595, 0.95598531],\n",
    "        [0.51994873, 0.52406595, 0.95598531],\n",
    "        [0.51994873, 0.52406595, 0.95598531]]))\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "K = 80\n",
    "(652239.615902771,\n",
    " 0.0,\n",
    " 582422,\n",
    " array([[0.20614538, 0.21375238, 0.58010224],\n",
    "        [0.20614538, 0.21375238, 0.58010224],\n",
    "        [0.20614538, 0.21375238, 0.58010224],\n",
    "        ...,\n",
    "        [0.20614538, 0.21375238, 0.58010224],\n",
    "        [0.20614538, 0.21375238, 0.58010224],\n",
    "        [0.20614538, 0.21375238, 0.58010224]]))\n",
    "\"\"\""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8",
   "language": "python",
   "name": "python-3.8"
  },
  "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}