salphagome - short.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: biopython in c:\\users\\drmn\\anaconda3\\lib\\site-packages (1.78)Note: you may need to restart the kernel to use updated packages.\n",
      "Requirement already satisfied: numpy in c:\\users\\drmn\\anaconda3\\lib\\site-packages (from biopython) (1.18.5)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "pip install biopython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Methods in the following modules are needed for a successful run of Salphagome\n",
    "import glob\n",
    "import pandas as pd\n",
    "import re\n",
    "from Bio import Entrez\n",
    "import time\n",
    "import xml.etree.ElementTree as ET\n",
    "import seaborn as sns; sns.set()\n",
    "import copy\n",
    "import numpy as np\n",
    "import subprocess\n",
    "from Bio import SeqIO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Edit the locations below to the location of your results folders\n",
    "files_typhi = glob.glob('typhi_results/*species.txt')\n",
    "files_ptyphiA = glob.glob('ptyphia_results/*species.txt')\n",
    "\n",
    "location_of_fasta = 'databases/just_viruses/just_viruses.fasta'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total Phages: 153\n"
     ]
    }
   ],
   "source": [
    "#List of sequences in the reference FASTA\n",
    "fasta_ids = []\n",
    "\n",
    "with open(location_of_fasta, \"r\") as handle:\n",
    "    for record in SeqIO.parse(handle, \"fasta\") :\n",
    "        fasta_ids.append(record.id)\n",
    "\n",
    "print(\"Total Phages: \" + str(len(fasta_ids)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make Master List\n",
    "master_list = []\n",
    "\n",
    "#Get read hit counts from control folder\n",
    "\n",
    "for file in files_typhi:\n",
    "    \n",
    "    name = str(file)\n",
    "    \n",
    "    list_of_seqs = copy.deepcopy(fasta_ids)\n",
    "    \n",
    "    with open(file, \"r\") as ifile:\n",
    "        for line in ifile:\n",
    "            \n",
    "            data = line.strip().split(' ')\n",
    "            data = data+[name, 'typhi']\n",
    "            master_list.append(data)\n",
    "            \n",
    "    for seq in list_of_seqs:\n",
    "        \n",
    "        data = [0, seq, name, 'typhi']\n",
    "        \n",
    "        master_list.append(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Get read hit counts from disease folder\n",
    "\n",
    "for file in files_ptyphiA:\n",
    "    \n",
    "    list_of_seqs = copy.deepcopy(fasta_ids)\n",
    "    \n",
    "    name = str(file)\n",
    "    \n",
    "    with open(file, \"r\") as ifile:\n",
    "        \n",
    "        for line in ifile:\n",
    "            \n",
    "            data = line.strip().split(' ')\n",
    "            data = data+[name, 'ptyphia']\n",
    "            \n",
    "            list_of_seqs.remove(data[1])\n",
    "            \n",
    "            \n",
    "            master_list.append(data)\n",
    "\n",
    "    for seq in list_of_seqs:\n",
    "        \n",
    "        data = [0, seq, name, 'ptyphia']\n",
    "        \n",
    "        master_list.append(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Name items in our current dataframe.\n",
    "\n",
    "df = pd.DataFrame(columns=['count', 'seq', 'file', 'type'], data= master_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function that gets accession Ids from the SRA.\n",
    "\n",
    "def get_srr_accession(df):\n",
    "    \"\"\"\n",
    "    Extract the SRR number from the file_location field of the dataframe.\n",
    "    \"\"\"\n",
    "    pattern = re.compile(r\"[SED]RR[0-9]+\")\n",
    "\n",
    "    file = df['file']\n",
    "\n",
    "    return re.search(pattern, file).group(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Include accession Ids in our dataframe.\n",
    "\n",
    "df['acc'] = df.apply(get_srr_accession, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert the counts to integer object.\n",
    "\n",
    "df['count'] = df['count'].astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function that gets the spot counts and use for normalization of the data.\n",
    "\n",
    "def get_spot_count(df, email, max_errors):\n",
    "    \"\"\"\n",
    "    Use the Entrez API to get the spot (read) count for that SRR number.\n",
    "    \"\"\"\n",
    "\n",
    "    srr_acc = df['acc']\n",
    "    \n",
    "    # get SRR id\n",
    "    Entrez.email = email\n",
    "\n",
    "    error_count =0\n",
    "\n",
    "    while error_count < max_errors:\n",
    "\n",
    "        try:\n",
    "\n",
    "            handle = Entrez.esearch(db=\"sra\",term=srr_acc)\n",
    "\n",
    "            record = Entrez.read(handle)\n",
    "\n",
    "            handle.close()\n",
    "\n",
    "            srr_id = record[\"IdList\"][0]\n",
    "\n",
    "            break\n",
    "\n",
    "        except:\n",
    "\n",
    "            print ('error occured collecting ID ', srr_acc, error_count)\n",
    "\n",
    "            time.sleep(10)\n",
    "\n",
    "            error_count = error_count +1\n",
    "\n",
    "\n",
    "    # get SRR summary\n",
    "\n",
    "    error_count =0\n",
    "\n",
    "    while error_count < max_errors:\n",
    "\n",
    "        try:\n",
    "\n",
    "            handle = Entrez.esummary(db='sra', id=srr_id)\n",
    "\n",
    "            record = Entrez.read(handle)\n",
    "\n",
    "            handle.close()\n",
    "\n",
    "            my_xml = record[0]['Runs']\n",
    "\n",
    "            # Parse XML\n",
    "            xml_object = ET.fromstringlist([\"<root>\", my_xml, \"</root>\"])\n",
    "\n",
    "            # Get total_spots (reads)\n",
    "            for child in xml_object:\n",
    "\n",
    "                if child.attrib['acc'] == srr_acc:\n",
    "\n",
    "                    #print (srr_acc, child.attrib['total_spots'])\n",
    "\n",
    "                    return int(child.attrib['total_spots'])\n",
    "\n",
    "        except:\n",
    "\n",
    "            print ('error occured collecting spot count', srr_acc, error_count)\n",
    "\n",
    "            time.sleep(10)\n",
    "\n",
    "            error_count = error_count +1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "data2 = my_list = list(df['acc'].unique())\n",
    "\n",
    "spot_df = pd.DataFrame(data2, columns=['acc'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          acc  spot_count\n",
      "0   DRR071000     1137918\n",
      "1   DRR071001      789639\n",
      "2   DRR071002      668472\n",
      "3   DRR071003      566361\n",
      "4  ERR2011049      702860\n"
     ]
    }
   ],
   "source": [
    "#Get spot count for each of the sequences\n",
    "\n",
    "spot_df['spot_count'] = spot_df.apply(get_spot_count, axis=1, args=['visiting.bioinformaticians@gmail.com', 3])\n",
    "print(spot_df.drop_duplicates(subset=['acc'], keep='last'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Write Spot Dataframe to CSV file\n",
    "\n",
    "spot_df.to_csv('spot_map.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "read_count_dict ={}\n",
    "\n",
    "for row in spot_df.iterrows():\n",
    "    read_count_dict[(row[1]['acc'])] = row[1]['spot_count']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function that normalizes the data.\n",
    "\n",
    "def normalise(df):\n",
    "    \n",
    "    srr = df['acc']\n",
    "    \n",
    "    read_count = read_count_dict[srr]\n",
    "    \n",
    "    return read_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normalize the data.\n",
    "\n",
    "df['read_count'] = df.apply(normalise, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "scrolled": true
   },
   "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>count</th>\n",
       "      <th>seq</th>\n",
       "      <th>file</th>\n",
       "      <th>type</th>\n",
       "      <th>acc</th>\n",
       "      <th>read_count</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>20647</td>\n",
       "      <td>NC_010463.1</td>\n",
       "      <td>typhi_results\\typhi_phages.DRR071000.sam.speci...</td>\n",
       "      <td>typhi</td>\n",
       "      <td>DRR071000</td>\n",
       "      <td>1137918</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>12323</td>\n",
       "      <td>NC_019488.1</td>\n",
       "      <td>typhi_results\\typhi_phages.DRR071000.sam.speci...</td>\n",
       "      <td>typhi</td>\n",
       "      <td>DRR071000</td>\n",
       "      <td>1137918</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>7853</td>\n",
       "      <td>NC_028699.1</td>\n",
       "      <td>typhi_results\\typhi_phages.DRR071000.sam.speci...</td>\n",
       "      <td>typhi</td>\n",
       "      <td>DRR071000</td>\n",
       "      <td>1137918</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>6990</td>\n",
       "      <td>NC_010393.1</td>\n",
       "      <td>typhi_results\\typhi_phages.DRR071000.sam.speci...</td>\n",
       "      <td>typhi</td>\n",
       "      <td>DRR071000</td>\n",
       "      <td>1137918</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1712</td>\n",
       "      <td>NC_010392.1</td>\n",
       "      <td>typhi_results\\typhi_phages.DRR071000.sam.speci...</td>\n",
       "      <td>typhi</td>\n",
       "      <td>DRR071000</td>\n",
       "      <td>1137918</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   count          seq                                               file  \\\n",
       "0  20647  NC_010463.1  typhi_results\\typhi_phages.DRR071000.sam.speci...   \n",
       "1  12323  NC_019488.1  typhi_results\\typhi_phages.DRR071000.sam.speci...   \n",
       "2   7853  NC_028699.1  typhi_results\\typhi_phages.DRR071000.sam.speci...   \n",
       "3   6990  NC_010393.1  typhi_results\\typhi_phages.DRR071000.sam.speci...   \n",
       "4   1712  NC_010392.1  typhi_results\\typhi_phages.DRR071000.sam.speci...   \n",
       "\n",
       "    type        acc  read_count  \n",
       "0  typhi  DRR071000     1137918  \n",
       "1  typhi  DRR071000     1137918  \n",
       "2  typhi  DRR071000     1137918  \n",
       "3  typhi  DRR071000     1137918  \n",
       "4  typhi  DRR071000     1137918  "
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Print dataframe's first few lines\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normalise_values(df):\n",
    "    \n",
    "    return (df['count'] / df['read_count']) * 1000000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['normalised_count'] = df.apply(normalise_values, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_csv('typhi.ptyphi.comparison.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1682100db20>"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEJCAYAAAB/pOvWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df1hUdd7/8efA6KTprkozSxFLt5s/Wu+C1Nomt6FMESXUHc3Mbt3b1i6tlk271zIgzF13NW9aKne5u7q3q/vO6l7RLJSlcV31JhPNJK9Qw64uE1LUYcCfKGAw5/tH32bjPgqDMg7g63FdXIfzOZ8z8z5eR15zPmfOORbDMAxERES+IyLcBYiISMejcBAREROFg4iImCgcRETEROEgIiImCgcRETFROIiIiIk13AW0l+PHz+D365INEZFgRERY6Nv36gsu7zLh4PcbCgcRkXaiYSURETFROIiIiInCQURETBQOIiJionAQEREThYOIiJgoHESkU/jkk50sWpTBJ5/sDHcpV4Quc52DiHRtq1a9zYEDX1JfX8fQocPDXU6XpyMHEekU6urqm00ltEIaDvn5+aSkpJCSksLzzz8PQHFxMampqSQlJZGTkxPoW1ZWhtvtZsyYMWRkZNDY2BjK0kREpAUhC4e6ujp+97vfsWLFCvLz89m5cyebNm0iPT2d3NxcCgsL2bNnD0VFRQDMnz+frKws1q9fj2EY5OXlhao0ERFpRcjCoampCb/fT11dHY2NjTQ2NtKrVy/i4uKIjY3FarWSmpqKx+OhsrKS+vp6EhISAHC73Xg8nlCVJiIirQjZCelevXrxxBNPMHbsWHr06MFtt91GVVUVdrs90MfhcOD1ek3tdrsdr9cbqtJERKQVIQuHffv28c4777B582Z69+7Nr3/9a8rLy7FYLIE+hmFgsVjw+/3nbW+LqKhe7Va7iHQ8kZGWwNRu7x3marq+kIXDhx9+iNPpJCoqCvhmqOi1114jMjIy0Mfn8+FwOIiOjsbn8wXaq6urcTgcbXq/mppa3bJbpAtrajICU5/vdJir6fwiIiwtfqgO2TmHwYMHU1xczNmzZzEMg02bNhEfH8+BAweoqKigqamJgoICXC4XMTEx2Gw2SkpKgG++5eRyuUJVmoiItCJkRw4//elP+eyzz3C73XTr1o2bb76ZtLQ0RowYQVpaGg0NDSQmJpKcnAxAdnY2mZmZ1NbWMmTIEGbMmBGq0kREpBUWwzC6xFiMhpVEura5cx/j6NHDREdfx4sv5oa7nE4vbMNKIiLSeSkcRETEROEgIiImCgcRETFROIiIiInCQURETBQOIiJionAQEREThYOIiJgoHERExEThICIiJgoHERExUTiIiIiJwkFEREwUDiIiYqJwEBERk5A9CW7VqlW8+eabgflDhw4xYcIERo0axZIlS2hoaGDs2LHMmzcPgLKyMjIyMjhz5gzDhw9n0aJFWK0hK09ERFoQsiOH+++/n/z8fPLz88nOziYqKopHHnmE9PR0cnNzKSwsZM+ePRQVFQEwf/58srKyWL9+PYZhkJeXF6rSRESkFZdlWOm5555j3rx5HDx4kLi4OGJjY7FaraSmpuLxeKisrKS+vp6EhAQA3G43Ho/ncpQmIiLnEfJwKC4upr6+nrFjx1JVVYXdbg8sczgceL1eU7vdbsfr9Ya6NBERuYCQD+r/5S9/YebMmQD4/X4sFktgmWEYWCyWC7a3RUsPyhaRzi8y0hKY2u29w1xN1xfScDh37hwff/wxS5cuBSA6OhqfzxdY7vP5cDgcpvbq6mocDkeb3qumpha/32ifwkWkw2lqMgJTn+90mKvp/CIiLC1+qA7psNLnn3/ODTfcQM+ePQGIj4/nwIEDVFRU0NTUREFBAS6Xi5iYGGw2GyUlJQDk5+fjcrlCWZqIiLQgpEcOBw8eJDo6OjBvs9lYunQpaWlpNDQ0kJiYSHJyMgDZ2dlkZmZSW1vLkCFDmDFjRihLExGRFlgMw+gSYzEaVhLp2ubOfYyjRw8THX0dL76YG+5yOr2wDiuJiEjnpHAQEREThYOIiJgoHERExEThICIiJgoHERExUTiIiIiJwkFEREwUDiIiYqJwEBERE4WDiIiYKBxERMRE4SAiIiYKBxERMVE4iIiIicJBRERMQhoOmzZtwu12M3bsWBYvXgxAcXExqampJCUlkZOTE+hbVlaG2+1mzJgxZGRk0NjYGMrSRESkBUGHQ2lpqamtuLj4gv0PHjzIwoULyc3NZe3atXz22WcUFRWRnp5Obm4uhYWF7Nmzh6KiIgDmz59PVlYW69evxzAM8vLyLmJzRESkPbQaDp999hl79+7l6aefDvy+d+9ePv30U5577rkLrrdhwwbGjRtHdHQ03bp1Iycnhx49ehAXF0dsbCxWq5XU1FQ8Hg+VlZXU19eTkJAAgNvtxuPxtNtGiohI21hb6/A///M/bN26laqqKn75y1/+Y0WrldGjR19wvYqKCrp168acOXM4cuQId999NwMGDMButwf6OBwOvF4vVVVVzdrtdjter/dit0lERC5Rq+Hw29/+FoCcnBzmzZsX9As3NTWxc+dOVqxYQc+ePXn00Ue56qqrsFgsgT6GYWCxWPD7/edtb4uWHpQtIp1fZKQlMLXbe4e5mq6v1XD41ty5c9mxYwcnT57EMIxAe1JS0nn7X3PNNTidTvr16wfAqFGj8Hg8REZGBvr4fD4cDgfR0dH4fL5Ae3V1NQ6Ho00bUlNTi99vtN5RRDqlpiYjMPX5Toe5ms4vIsLS4ofqoMMhKyuLoqIi4uLiAm0Wi+WC4XDPPffw9NNPc+rUKa6++mq2bNlCcnIyr776KhUVFVx//fUUFBQwadIkYmJisNlslJSUMGzYMPLz83G5XG3YTBERaU9Bh8PWrVspLCykV6/ghm/i4+OZNWsW06ZN4+uvv2bEiBE8+OCD9O/fn7S0NBoaGkhMTCQ5ORmA7OxsMjMzqa2tZciQIcyYMePitkhERC5Z0OFw3XXXBR0M35o8eTKTJ09u1uZ0Olm7dq2p7+DBg1m9enWbXl9EREIj6HAYOnQo8+bN45577uGqq64KtF9oWElERDqvoMNh165dAKxatSrQ1tI5BxER6byCDocVK1aEsg4REelAgg6Hb++N9H9lZma2WzEiItIxBH1vpT59+gR+rr76anbs2BHKukREJIyCPnL47q0zAB555BEeffTRdi9IRETC76Jv2d2rVy+qqqrasxYREekgLuqcg2EY7N27l/79+4ekKBERCa+gw6FPnz7N5sePH8/48ePbvSAREQm/Np1zOHPmDHv37qWxsZFbbrmlzVdMi4hI5xB0OJSWlvLYY49xzTXX0NTUhNfr5ZVXXmHo0KGhrE9ERMIg6HB4/vnnyc7O5o477gBg27ZtLF26VI/zFBHpgoL+ttKZM2cCwQDf3ECvrq4uJEWJiEh4BR0OFouFysrKwPyhQ4eaPbhHRES6jqCHlR5//HEeeOABnE4nFouFDz/8kIULF4ayNhERCZOgw2HUqFH079+f7du34/f7mT17Nj/60Y9CWZuIiIRJ0OGwb98+/vCHP/Dqq6/y+eef8+STT5KTk6ML4URCrO/3u2Ptbgt3GWEXGWkJTO323mGuJvwazzVw/OS5kL1+0OHw3HPP8Ytf/AKAQYMGkZaWxsKFC1u8lff06dM5duwYVus3b/Ob3/yGM2fOsGTJEhoaGhg7dizz5s0DoKysjIyMDM6cOcPw4cNZtGhRYD2RK5m1u42SZbPCXUbYNRz3Bqb694BhT/0ZCF04BH1Cuq6ujtGjRwfmR40aRW1t7QX7G4ZBeXk5+fn5gZ9BgwaRnp5Obm4uhYWF7Nmzh6KiIgDmz59PVlYW69evxzAMfUVWRCSM2vRtpX379gXm9+/fT0TEhVf/8ssvAXj44YcZP348b775JqWlpcTFxREbG4vVaiU1NRWPx0NlZSX19fUkJCQA4Ha78Xg8F7tNIiJyiYIet3niiSeYPn06AwcOxGKxsH//frKzsy/Y/9SpUzidTp599lm+/vprZsyYwaxZs7Db7YE+DocDr9dLVVVVs3a73Y7X623ThkRF6VYeInJlCeW5l6DD4Z577sHj8fDJJ58QGRlJfHw8UVFRAOzYsYPbb7+9Wf9bb72VW2+9NTA/efJkXn75ZYYNGxZoMwwDi8WC3+/HYrGY2tuipqYWv99o0zoinYFOvsqF+HynL3rdiAhLix+q2/Q8h6ioKEaPHs3IkSMDwQCwZMkSU9+dO3eybdu2wLxhGMTExODz+QJtPp8Ph8NBdHR0s/bq6mocDkdbShMRkXZ00Q/7+S7DMH9iP336NMuWLaOhoYHa2lreffddnnzySQ4cOEBFRQVNTU0UFBTgcrmIiYnBZrNRUlICQH5+Pi6Xqz1KExGRi9Au3xU93xDQPffcw6effsrEiRPx+/1MmzaNW2+9laVLl5KWlkZDQwOJiYkkJycDkJ2dTWZmJrW1tQwZMoQZM2a0R2kiInIRQnohwdy5c5k7d26zNqfTydq1a019Bw8ezOrVq0NZjoiIBKldhpVERKRrCdk5BxER6bzaJRx+8pOftMfLiIhIB9HqOYdnnnmmxeVLlixptY+IiHQurR45DBgwgAEDBnD69Gk+//xzBg0axI9//GPKy8tpamq6HDWKiMhl1uqRw8MPPwzAhg0beOutt+jRowcAU6ZM0ddNRUS6qKDPOdTU1NC9e/fAvMVi4fjx4yEpSkREwivo6xycTiezZs3ivvvuwzAM8vPzGTlyZChrExGRMAk6HJ599lneeustNmzYAMDYsWOZOnVqyAoTEZHwCTocrFYrY8aM4YYbbmDEiBFUVVW1+DwHERHpvIL+6/6///u/TJ06lUWLFnHs2DFSUlL4+9//HsraREQkTIIOhz/96U/k5eXxve99D4fDwdtvv83LL78cytpERCRMgg6HpqamZs9YuOmmm9r8QB4REekcgg6HHj16cPjw4UAg7Ny5E5vNFrLCREQkfII+If3rX/+ahx9+GJ/PxwMPPEB5eTnLly8PZW0iIhImQYfDrbfeSl5eHrt27cLv9xMfH0+/fv1CWZuIiIRJ0MNKdXV1fPnllyQmJnL06FGys7M5fPhwq+s9//zzLFiwAIDi4mJSU1NJSkoiJycn0KesrAy3282YMWPIyMigsbHxIjZFRETaS9Dh8Mwzz7Bx40Z2797Nn//8Z6699lqeffbZFtfZtm0b7777LgD19fWkp6eTm5tLYWEhe/bsoaioCID58+eTlZXF+vXrMQyDvLy8S9gkERG5VEGHw8GDB/m3f/s3Nm3axM9+9jPS0tI4ceLEBfufOHGCnJwc5syZA0BpaSlxcXHExsZitVpJTU3F4/FQWVlJfX09CQkJALjdbjwezyVuloiIXIqgzzl8O9Tz4Ycf8vTTT9PU1MTZs2cv2D8rK4t58+Zx5MgRAKqqqrDb7YHlDocDr9drarfb7Xi93jZvSFRUrzavIyLSmdntvUP22kGHQ0JCAuPGjSMyMpKhQ4fy85//nDvvvPO8fVetWsW1116L0+lkzZo1APj9/mbXRRiGgcViuWB7W9XU1OL363Gl0vWE8g+AdG4+3+mLXjciwtLih+qgwyErK4tdu3YxaNAgIiIi+MUvfkFiYuJ5+xYWFuLz+ZgwYQInT57k7NmzVFZWEhkZGejj8/lwOBxER0fj8/kC7dXV1c0uthMRkcuv1XDYtm0bTqeTjRs3Bua/tWHDBpKSkkzrvP7664Hf16xZw44dO1i0aBFJSUlUVFRw/fXXU1BQwKRJk4iJicFms1FSUsKwYcPIz8/H5XK1x7aJiMhFajUc/vrXv+J0OlmxYoVpmcViOW84nI/NZmPp0qWkpaXR0NBAYmIiycnJAGRnZ5OZmUltbS1DhgzRE+ZERMLMYhhGlxio1zkH6ars9t6ULJsV7jLCbtmWg1SfbeSanlaeuis23OWE3bCn/hzecw7ffhX1Ql555ZW2VyUiIh1aq+EwZsyYy1GHiIh0IK2Gw89+9rPzthuGQUVFRbsXJCIi4Rf0V1n/8pe/sGzZMurq6gJt/fr1Y+vWrSEpTEREwifocHj11Vd5/fXX+Y//+A/mzp3L5s2bOXr0aChrExGRMAn63kp9+vQhPj6em266iZqaGh599FE+/vjjUNYmIiJhEnQ4WK1WTp48SVxcHKWlpcA3jw4VEZGuJ+hwmDJlCrNnz+buu+9m5cqVuN1u+vfvH8raREQkTII+5zB58mTGjRtHz549WblyJbt37+auu+4KZW0iIhImQYdDfX09mzdvbvYMh9WrV/PQQw+FpDAREQmfoMNhzpw5nDp1iuuvvz7QZrFYFA4iIl1Q0OHg9XopLCy8qGctiIhI5xL0CemBAwdSXV0dylpERKSDCPrIITk5mbFjxzJw4ECs1n+s9sYbb4SkMBERCZ+gw+FPf/oTs2fP5oc//GEo6xERkQ4g6HDo0aMHjzzySChrERGRDiLocw533nknb731FlVVVZw4cSLw05KXXnqJcePGkZKSEnh0aHFxMampqSQlJZGTkxPoW1ZWhtvtZsyYMWRkZNDY2HiRmyQiIpcq6COH119/nXPnzvHb3/420GaxWCgrKztv/x07drB9+3bWrl1LY2Mj48aNw+l0kp6ezooVK7j22muZPXs2RUVFJCYmMn/+fBYvXkxCQgLp6enk5eUxbdq0S99CERFps6CPHN5++2327dvX7OdCwQBw++2388Ybb2C1WqmpqaGpqYlTp04RFxdHbGwsVquV1NRUPB4PlZWV1NfXk5CQAIDb7cbj8Vz61kmbffLJThYtyuCTT3aGuxQRCaOgw2H+/PltfvFu3brx8ssvk5KSgtPppKqqCrvdHljucDjwer2mdrvdjtfrbfP7yaVbteptysr2smrV2+EuRUTCKOhhpUGDBrFu3TqGDRtGz549A+19+vRpcb1f/epXPPLII8yZM4fy8vJmF9EZhoHFYsHv95+3vS1aelC2BO/cuYbA1G7vHeZqRKQlofw/GnQ4bNy40TTU09I5h/3793Pu3DluuukmevToQVJSEh6Ph8jIyEAfn8+Hw+EgOjoan88XaK+ursbhcLRpQ2pqavH7jTatI2ZNTUZg6vOdDnM1AqH9AyCd26X8H42IsLT4oTroYaXdu3e36ZzDoUOHyMzM5Ny5c5w7d46NGzcydepUDhw4QEVFBU1NTRQUFOByuYiJicFms1FSUgJAfn4+LperDZspIiLtKegjB7/fz2uvvcYHH3xAY2MjI0aMYM6cOc2ulv6uxMRESktLmThxIpGRkSQlJZGSkkK/fv1IS0ujoaGBxMREkpOTAcjOziYzM5Pa2lqGDBnCjBkz2mcLRUSkzYIOhxdeeIF9+/bx85//HL/fz8qVK3n++efJyMi44DppaWmkpaU1a3M6naxdu9bUd/DgwaxevboNpYvIlcRmjWg2ldAKOhy2bNnCO++8Q7du3QC4++67GT9+fMgKExH5rqQb+1J04CSJ//T9cJdyRQg6HAzDCAQDQPfu3ZvNi4iE0k32ntxk79l6R2kXQR+fDR48mN///vd89dVXHDx4kN///vcMHDgwlLWJiEiYBB0OCxcu5NSpUzz44INMmTKF48ePk5WVFcraREQkTIIOhy+++ILKykr69evHNddcw759+/iXf/mXUNYmIiJhEvQ5h6ysLNxuNz/+8Y/1qFARkS4u6HCwWq3MnDkzlLWIiEgHEfSw0oABA/j8889DWYuIiHQQQR85HDx4kEmTJnHddddhs9kC7evWrQtJYSIiEj5Bh8O8efNCWYeIiHQgQYfD7bffHso6RESkA9FNSkRExEThICIiJkEPK3V1vb93FVfZdK+oyEhLYKqHzEB9w9ecPlUf7jJELjuFw/93la0b0556K9xlhF119TdPljpafVr/HsDbyx7iNAoHufJoWElERExCGg5//OMfSUlJISUlhWXLlgFQXFxMamoqSUlJ5OTkBPqWlZXhdrsZM2YMGRkZNDY2hrI0ERFpQcjCobi4mA8//JB3332X9957j71791JQUEB6ejq5ubkUFhayZ88eioqKAJg/fz5ZWVmsX78ewzDIy8sLVWkiItKKkIWD3W5nwYIFgYcC/ehHP6K8vJy4uDhiY2OxWq2kpqbi8XiorKykvr6ehIQEANxuNx6PJ1SliYhIK0IWDgMGDAj8sS8vL+f999/HYrFgt9sDfRwOB16vl6qqqmbtdrsdr9cbqtJERKQVIf+20hdffMHs2bN56qmniIyMpLy8PLDMMAwsFgt+v7/ZbcC/bW+LqKhe7VWySDP6Sq90VKHcN0MaDiUlJfzqV78iPT2dlJQUduzYgc/nCyz3+Xw4HA6io6ObtVdXV+NwONr0XjU1tfj9xkXXqj8AciE+3+mwvr/2TbmQS9k3IyIsLX6oDtmw0pEjR3j88cfJzs4mJSUFgPj4eA4cOEBFRQVNTU0UFBTgcrmIiYnBZrNRUlICQH5+Pi6XK1SliYhIK0J25PDaa6/R0NDA0qVLA21Tp05l6dKlpKWl0dDQQGJiIsnJyQBkZ2eTmZlJbW0tQ4YMYcaMGaEqTUREWhGycMjMzCQzM/O8y9auXWtqGzx4MKtXrw5VOSIi0ga6QlpEREwUDiIiYqJwEBERE4WDNGOJ7NZsKiJXJoWDNNPruqF06xVNr+uGhrsUEQkjPc9BmrF9Pxbb92PDXYaIhJmOHERExEThICIiJgoHERExUTiIiIiJwkFEREwUDiIiYqJwEBERE4WDiIiYKBxERMRE4SAiIiYhDYfa2lruu+8+Dh06BEBxcTGpqakkJSWRk5MT6FdWVobb7WbMmDFkZGTQ2NgYyrJERKQVIQuHTz/9lAcffJDy8nIA6uvrSU9PJzc3l8LCQvbs2UNRUREA8+fPJysri/Xr12MYBnl5eaEqS0REghCycMjLy2PhwoU4HA4ASktLiYuLIzY2FqvVSmpqKh6Ph8rKSurr60lISADA7Xbj8XhCVZaIiAQhZHdl/d3vftdsvqqqCrvdHph3OBx4vV5Tu91ux+v1hqosEREJwmW7Zbff78disQTmDcPAYrFcsL2toqJ6tUudIv+X3d473CWInFco983LFg7R0dH4fL7AvM/nw+FwmNqrq6sDQ1FtUVNTi99vXHR9+gMgF+LznQ7r+2vflAu5lH0zIsLS4ofqy/ZV1vj4eA4cOEBFRQVNTU0UFBTgcrmIiYnBZrNRUlICQH5+Pi6X63KVJSIi53HZjhxsNhtLly4lLS2NhoYGEhMTSU5OBiA7O5vMzExqa2sZMmQIM2bMuFxliYjIeYQ8HDZt2hT43el0snbtWlOfwYMHs3r16lCXIiIiQdIV0iIiYqJwEBERE4WDiIiYKBxERMRE4SAiIiYKBxERMVE4iIiIicJBRERMFA4iImKicBAREROFg4iImCgcRETEROEgIiImCgcRETFROIiIiInCQURETBQOIiJi0qHCYd26dYwbN46kpCTeeuutcJcjInLFumzPkG6N1+slJyeHNWvW0L17d6ZOncpPfvITbrzxxnCXJiJyxekw4VBcXMwdd9xBnz59ABgzZgwej4df/vKXQa0fEWG55Bqu6Xv1Jb+GdD3tsW9dqu7fiwp3CdIBXcq+2dq6HSYcqqqqsNvtgXmHw0FpaWnQ6/dthz/sLz8z8ZJfQ7qeqKhe4S6Bm+c8H+4SpAMK5b7ZYc45+P1+LJZ/JJlhGM3mRUTk8ukw4RAdHY3P5wvM+3w+HA5HGCsSEblydZhwuPPOO9m2bRvHjh2jrq6Ov/3tb7hcrnCXJSJyReow5xx+8IMfMG/ePGbMmMHXX3/N5MmTueWWW8JdlojIFcliGIYR7iJERKRj6TDDSiIi0nEoHERExEThICIiJgoHERExUThcAU6fPs3jjz/e5vWWL1/O8uXLTe0bN27kpZdeao/SRJg+ffpFrbdmzRoWLFhgat+9ezcZGRmXWtYVr8N8lVVC5+TJk5SVlbXb6917773ce++97fZ6cmXbsWNHu77ezTffzM0339yur3klUjhcARYvXkxVVRWPP/44N954I/PmzQNgwYIFuFwuPvjgA2w2G7t37+bMmTM8+uijTJz4zX2mSktLmTp1Kl6vF7fbTVpaGmvWrGHHjh0sXbo0nJslncxHH31Ebm4uVquVQ4cOccstt9CzZ08A7r//fqZMmcL27dt54YUXgG+OXG02Gw0NDRw+fJj9+/dz/PhxHnjgAWbNmgVARUUF06dP5/DhwzidThYvXsxHH33EH//4R1asWMGOHTvIycmhvr6eU6dO8cwzzzBq1Kiw/Rt0JhpWugJkZmbicDhYsGAB69atwzAM6urq2L59e+AI4ODBg6xcuZL//u//ZtmyZYFbmdTU1PDGG2/wzjvv8Nprr1FbWxvOTZFObteuXWRkZODxeGhoaOCGG24AYNWqVYwbN45t27YF9rGCggImTJgAwJ49e3j99ddZs2YNK1euZO/evQAcOXKE5cuX8/777/PBBx/wxRdfNHu/N998k8WLF/Puu++yePFiDYe2gY4criCxsbHExMTw8ccfc/jwYRITE7HZbAC43W66detGdHQ0Q4cOpaSkBIC77rqL7t27069fP/r27cvJkyfDuQnSyd122230798fgAkTJpCXlxdYdvXVV5OYmMiGDRuIjY0lNjaWH/zgBwDcd999XH31N3deHjlyJNu3b6dv374MHz48cJv/H/7whxw/frzZ+/37v/87mzdvxuPx8Omnn3LmzJnLsZldgo4crjCTJk2ioKCAgoIC3G53oD0yMjLwu9/vx2r95nPDt1MAi8WCLqiXS/Hd/cwwjGbz8I/9c926dS3un9/Ot7Z/Tps2jdLSUv75n/+ZOXPmtOu2dHUKhyuA1WqlsbERgOTkZLZt20Z1dTXx8fGBPu+//z6GYVBZWUlpaSnDhg0LV7nShZWUlOD1evH7/bz33nu4XC4iIyMD++fw4cM5evQoH330UbNzA3//+985d+4cJ0+eZPPmzfz0pz9t9b1OnDhBeXk5TzzxBC6Xi40bN9LU1BSybetqNKx0BYiKiuK6665j+vTprFixgoSEBAYOHNuEjCsAAAMlSURBVNisT319PZMmTeLcuXP85je/oW/fvmGqVroyh8PBU089hdfrZcSIEdx///1s3bqVCRMmsGbNGmw2G6NHj+bEiRN07949sJ7NZmPatGnU1tYye/ZsbrzxxlYfBtanTx8mT55MSkoKVquVO+64g/r6es6ePRs4ES4XphvvXUEMw+DMmTM88MAD/Nd//VfgyXsLFizg9ttvb3YYL9LevvstovMxDIOvv/6amTNnkp6ezpAhQwAC19qkpaVdtlpFw0pXlN27dzNy5EimTJnS7JGsIh2Bz+djxIgRxMfHB4JBwkdHDiIiYqIjBxERMVE4iIiIicJBRERMFA4iF+Hhhx/m2LFj4S5DJGQUDiIXYevWreEuQSSkFA4ibfTMM88A39x3atCgQfj9fgDq6upwOp0cO3aMkSNH8sILL+B2uxk9ejRvv/12YP1NmzZx//33M3HiRKZOncquXbvCsh0iLdEV0iJttGTJEtasWcOWLVuYOXMmW7ZsITExkb/+9a84nU769esHfPMcjXfeeQev18vEiRMZNmwYNpuNnJwc3njjDfr27csXX3zBzJkz+dvf/qardqVDUTiIXIKHHnqIvLw8EhMTWblyJU899VRg2bRp07BYLERHR3PXXXexdetWbDYbVVVV/Ou//mugn8Vi4auvvmLw4MFh2AKR81M4iFyC1NRU/vCHP7B9+3bOnj3LbbfdFlj23TuG+v1+IiIi8Pv9OJ1OXnzxxcCyI0eO4HA4LmvdIq3ROQeRi/DtnUR79OjB+PHjSU9PZ+rUqc36vPfeewAcPnyYrVu34nK5cDqdbN26lf379wNQVFTE+PHjqa+vv+zbINIS3T5D5CI8+eST7N27l+XLl+P3+5k6dSpbtmyhd+/ewDcPpLnlllv46quvqK+v57HHHuO+++4Dvrk9+iuvvIJhGFitVtLT0xk+fHg4N0fEROEgcgkMw+A///M/qaysZNGiRYH2kSNH8tJLL+lB99Jp6ZyDyCW49957cTgc5ObmhrsUkXalIwcRETHRCWkRETFROIiIiInCQURETBQOIiJionAQEREThYOIiJj8Pz5ecJXWV6q3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.barplot(x='type', y='normalised_count', data=df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "grouped = df.groupby(['type', 'seq']).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "NOTABLE PHAGE  COVERAGE\n",
      "NC_019488.1 :  89.5397\n",
      "NC_019545.1 :  68.1818\n",
      "NC_010392.1 :  58.427\n",
      "NC_018279.1 :  57.7778\n",
      "NC_010463.1 :  90.5455\n",
      "NC_002371.2 :  58.3333\n",
      "NC_004348.1 :  57.1429\n",
      "NC_017985.1 :  57.1429\n",
      "NC_010393.1 :  80.4511\n",
      "NC_028699.1 :  89.3939\n",
      "NC_010391.1 :  51.6129\n",
      "\n",
      "NOTABLE PHAGE  COVERAGE\n",
      "NC_031129.1 :  54.0984\n",
      "NC_030919.1 :  52.6316\n",
      "NC_019488.1 :  88.9362\n",
      "NC_019545.1 :  72.8814\n",
      "NC_005841.1 :  52.9412\n",
      "NC_010392.1 :  70.1493\n",
      "NC_011802.1 :  68.75\n",
      "NC_014900.1 :  63.1579\n",
      "NC_018279.1 :  62.5\n",
      "NC_010463.1 :  91.6981\n",
      "NC_018275.1 :  52.9412\n",
      "NC_004348.1 :  65.625\n",
      "NC_005340.1 :  55.5556\n",
      "NC_017985.1 :  80\n",
      "NC_010393.1 :  76.6917\n",
      "NC_028699.1 :  89.0052\n",
      "NC_031946.1 :  55\n",
      "\n",
      "NOTABLE PHAGE  COVERAGE\n",
      "NC_019488.1 :  90.9465\n",
      "NC_019545.1 :  69.1176\n",
      "NC_005841.1 :  58.8235\n",
      "NC_010392.1 :  58.3333\n",
      "NC_011802.1 :  57.1429\n",
      "NC_014900.1 :  56.5217\n",
      "NC_018279.1 :  63.4146\n",
      "NC_010463.1 :  92.1642\n",
      "NC_002371.2 :  53.5714\n",
      "NC_004348.1 :  63.8889\n",
      "NC_017985.1 :  53.3333\n",
      "NC_010393.1 :  74.4681\n",
      "NC_028699.1 :  89.0547\n",
      "\n",
      "NOTABLE PHAGE  COVERAGE\n",
      "NC_030919.1 :  57.6923\n",
      "NC_019488.1 :  90.6504\n",
      "NC_019545.1 :  68.2353\n",
      "NC_010392.1 :  70.4545\n",
      "NC_018279.1 :  53.7037\n",
      "NC_010463.1 :  91.8819\n",
      "NC_002371.2 :  60.6061\n",
      "NC_004348.1 :  58.6957\n",
      "NC_017985.1 :  60\n",
      "NC_010393.1 :  76.129\n",
      "NC_028699.1 :  90.8257\n",
      "NC_010391.1 :  51.0638\n",
      "NC_004313.1 :  51.2195\n",
      "NC_031946.1 :  55.5556\n",
      "\n",
      "NOTABLE PHAGE  COVERAGE\n",
      "NC_030919.1 :  52.1739\n",
      "NC_019488.1 :  90.2041\n",
      "NC_019545.1 :  77.7778\n",
      "NC_005841.1 :  52.6316\n",
      "NC_010392.1 :  56.5217\n",
      "NC_018843.1 :  95.5451\n",
      "NC_018279.1 :  52\n",
      "NC_010463.1 :  90.9091\n",
      "NC_017985.1 :  60\n",
      "NC_010393.1 :  76.6423\n",
      "NC_028699.1 :  85.6436\n",
      "NC_031946.1 :  54.5455\n",
      "\n",
      "NOTABLE PHAGE  COVERAGE\n",
      "NC_030919.1 :  78.5714\n",
      "NC_019488.1 :  90.7563\n",
      "NC_019545.1 :  77.9661\n",
      "NC_005841.1 :  66.6667\n",
      "NC_010392.1 :  59.2593\n",
      "NC_011802.1 :  55\n",
      "NC_014900.1 :  60\n",
      "NC_018279.1 :  70.2703\n",
      "NC_010463.1 :  92.0152\n",
      "NC_018275.1 :  61.1111\n",
      "NC_002371.2 :  56\n",
      "NC_004348.1 :  74.1935\n",
      "NC_017985.1 :  80\n",
      "NC_010393.1 :  77.5362\n",
      "NC_028699.1 :  89.3401\n",
      "NC_010391.1 :  53.3333\n",
      "NC_031946.1 :  66.6667\n",
      "\n",
      "NOTABLE PHAGE  COVERAGE\n",
      "NC_019488.1 :  90.3361\n",
      "NC_019545.1 :  67.6471\n",
      "NC_010392.1 :  60.2273\n",
      "NC_018279.1 :  55.3191\n",
      "NC_010463.1 :  90.8088\n",
      "NC_004348.1 :  61.1111\n",
      "NC_017985.1 :  61.5385\n",
      "NC_010393.1 :  78.6765\n",
      "NC_028699.1 :  86.6337\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#\"then calculate the coverage\"\n",
    "\n",
    "COMMAND = '''for f in bedgraphs/*.bedgraph; do awk 'NR>1{if ($4 < 3) {mylo_phage[$1]++;} else {myhi_phage[$1]++;}} END {print \"\"; print \"NOTABLE PHAGE\", \" COVERAGE\"; {for (i in myhi_phage) if (mylo_phage[i] < myhi_phage[i]) print i, \": \", 100*myhi_phage[i]/(myhi_phage[i]+mylo_phage[i]);}}' $f; done'''\n",
    "\n",
    "(o,e) = subprocess.Popen(COMMAND, shell=True, stdout=subprocess.PIPE,stderr=subprocess.PIPE,universal_newlines=True).communicate()\n",
    "\n",
    "print(o)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import rpy2\n",
    "# %load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-56-1e9236fa937f>, line 9)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-56-1e9236fa937f>\"\u001b[1;36m, line \u001b[1;32m9\u001b[0m\n\u001b[1;33m    data$cov[data$cov < 0.5] <- 0\u001b[0m\n\u001b[1;37m        ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "#Calculate percentage of WGS of Salmonella strain that carry each phage\n",
    "%%R\n",
    "\n",
    "library(plyr)\n",
    "library(reshape2)\n",
    "library(reshape)\n",
    "data=read.csv(\"coverage.csv\",header=T,check.names = F)\n",
    "\n",
    "data$cov[data$cov < 0.5] <- 0\n",
    "data$cov[data$cov >= 0.5] <- 1\n",
    "samples_count <- count(data, c('seq','strain'))\n",
    "\n",
    "cov_count=aggregate(cov~seq+strain, data, sum)\n",
    "\n",
    "merged=merge(samples_count, cov_count, all.x = TRUE)\n",
    "merged$percentage=merged$cov*100/merged$freq\n",
    "casted=cast(merged, seq~strain,value=\"percentage\")\n",
    "write.csv(casted,\"new calculation.csv\")"
   ]
  },
  {
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}