enterococcus_metadata_cleaning.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import requests\n",
    "import numpy as np\n",
    "import bs4\n",
    "import re\n",
    "from pathlib import Path\n",
    "from collections import Counter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gouliouris et. al. 2018 UK Enterococcus Paper Metadata\n",
    "\n",
    "## Add sample and project accessions to Gouliouris 2018 metadata\n",
    "\n",
    "Starting from table S2 in https://mbio.asm.org/content/9/6/e01780-18 which contains `run_accessions` for deposited raw sequencing data in ENA, use the ENA API to automatically add `study_accession` and `sample_accession`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# table S2 from https://mbio.asm.org/content/9/6/e01780-18\n",
    "gouliouris_2018_uk_metadata = pd.read_csv('datasets/inline-supplementary-material-3.csv', sep='\\t', skiprows=1)\n",
    "gouliouris_2018_uk_metadata = gouliouris_2018_uk_metadata.rename(columns={'Accession number': 'Run_accession', \n",
    "                                          'Isolate ID': 'Isolate_name',\n",
    "                                          'ST': 'Sequence_Type'})\n",
    "\n",
    "# break into 500 accession chunks to more efficiently query the API\n",
    "api_responses = \"\"\n",
    "\n",
    "run_accessions = gouliouris_2018_uk_metadata.loc[gouliouris_2018_uk_metadata['Run_accession'].str.startswith('ERR'), 'Run_accession'].values\n",
    "for run_accession_chunk in [run_accessions[i:i + 500] for i in range(0, len(run_accessions), 500)]:\n",
    "    run_accession_chunk = \",\".join(run_accession_chunk)\n",
    "    api_response = requests.get(f\"https://www.ebi.ac.uk/ena/browser/api/xml/{run_accession_chunk}\")    \n",
    "    api_responses += api_response.text\n",
    "\n",
    "# parse combined xml records into a soup for easier traversal \n",
    "api_response_soup = bs4.BeautifulSoup(api_responses, 'lxml')\n",
    "\n",
    "\n",
    "def add_study_and_sample_metadata(row, api_response_soup):\n",
    "    \"\"\"\n",
    "    Use the ENA API to get the study and sample accessions\n",
    "    \"\"\"\n",
    "    run_accession = row['Run_accession']\n",
    "        \n",
    "    run_soup = api_response_soup.find(accession=run_accession)\n",
    "    if run_soup:\n",
    "        for xref in run_soup.find_all('xref_link'):\n",
    "            if xref.db.text == 'ENA-STUDY':\n",
    "                row['Study_accession'] = xref.id.text.strip()\n",
    "            elif xref.db.text == 'ENA-SAMPLE':\n",
    "                row['Sample_accession'] = xref.id.text.strip()\n",
    "    else:\n",
    "        row['Study_accession'] = np.nan\n",
    "        row['Sample_accession'] = np.nan\n",
    "    return row\n",
    "\n",
    "gouliouris_2018_uk_metadata = gouliouris_2018_uk_metadata.apply(lambda x: add_study_and_sample_metadata(x, api_response_soup), axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's grab the collection date to be consistent with Raven 2016"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_accessions = gouliouris_2018_uk_metadata.loc[~gouliouris_2018_uk_metadata['Sample_accession'].isna(), 'Sample_accession'].values\n",
    "for sample_accession_chunk in [sample_accessions[i:i + 500] for i in range(0, len(sample_accessions), 500)]:\n",
    "    sample_accession_chunk_text = \",\".join(sample_accession_chunk)\n",
    "    api_response = requests.get(f\"https://www.ebi.ac.uk/ena/browser/api/xml/{sample_accession_chunk_text}\")    \n",
    "    api_responses += api_response.text\n",
    "\n",
    "# parse combined xml records into a soup for easier traversal \n",
    "api_response_soup = bs4.BeautifulSoup(api_responses, 'lxml')\n",
    "\n",
    "collection_date_dict = {}\n",
    "for sample in api_response_soup.find_all('sample'):\n",
    "    sample_acc = sample['accession']\n",
    "    for tag in sample.find_all('tag'):\n",
    "        if tag.text == 'collection_date':\n",
    "            collection_date_dict[sample_acc] = tag.findNext('value').text\n",
    "            \n",
    "            \n",
    "gouliouris_2018_uk_metadata['Isolation_date'] = gouliouris_2018_uk_metadata['Sample_accession'].apply(lambda x: collection_date_dict[x] if x in collection_date_dict else np.nan)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then let's tidy up the UK metadata sheet by renaming some fields and adding some extra metadata to make our life easier merging in with the alberta data later. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# tidy up UK metadata for merging\n",
    "gouliouris_2018_uk_metadata = gouliouris_2018_uk_metadata.rename(columns={'Origin': 'Origin',\n",
    "                                          'BAPS group': 'BAPS_group',\n",
    "                                          'Location': 'Location',\n",
    "                                          'Ampicillin resistance': 'Ampicillin',\n",
    "                                          'Vancomycin resistance': 'Vancomycin'})\n",
    "\n",
    "gouliouris_2018_uk_metadata.loc[gouliouris_2018_uk_metadata['Removed in deduplication'] == 'Removed', 'Metadata_status'] = 'Removed for deduplication in original paper (10.1128/mBio.01780-18)'\n",
    "gouliouris_2018_uk_metadata = gouliouris_2018_uk_metadata.drop('Removed in deduplication', axis=1)\n",
    "\n",
    "# all in this paper are E. faecium\n",
    "gouliouris_2018_uk_metadata['Species'] = 'Enterococcus faecium'\n",
    "gouliouris_2018_uk_metadata['Country/Province'] = 'United Kingdom'\n",
    "\n",
    "# drop reference strains\n",
    "gouliouris_2018_uk_metadata = gouliouris_2018_uk_metadata.loc[gouliouris_2018_uk_metadata['Origin'] != 'Reference strain']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add read data for Gouliouris 2018 \n",
    "\n",
    "Then parse the read data for all the re-downloaded reads by FM (instead of those from HS and archived on `deivos.research.cs.dal.ca`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_reads(read_dir):\n",
    "    reads_data = {'Run_accession': [], 'R1': [], 'R2': []}\n",
    "    for read_dir in Path(read_dir).glob(\"*/\"):\n",
    "        run_accession = str(read_dir.name)\n",
    "        reads_data['Run_accession'].append(run_accession)\n",
    "\n",
    "        r1 = read_dir / f\"{run_accession}_1.fastq.gz\"\n",
    "        if r1.exists():\n",
    "            reads_data['R1'].append(str(r1))\n",
    "        else:\n",
    "            raise ValueError(f\"{r1} is missing\")\n",
    "\n",
    "        r2 = read_dir / f\"{run_accession}_2.fastq.gz\"\n",
    "        if r2.exists():\n",
    "            reads_data['R2'].append(str(r2))\n",
    "        else:\n",
    "            raise ValueError(f\"{r2} is missing\")\n",
    "    return pd.DataFrame(reads_data)  \n",
    "\n",
    "gouliouris_2018_uk_reads = get_reads('genomes/gouliouris_2018_uk')\n",
    "\n",
    "gouliouris_2018_uk_metadata_and_reads = pd.merge(gouliouris_2018_uk_metadata, \n",
    "                                                 gouliouris_2018_uk_reads, \n",
    "                                                 on='Run_accession', \n",
    "                                                 validate='one_to_one', \n",
    "                                                 how='outer')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Raven et. al. 2016 UK Hospital Enterococcus Paper Metadata\n",
    "\n",
    "## Add sample and project accessions to Raven 2016 metadata\n",
    "\n",
    "Starting from table S1 in https://www.nature.com/articles/nmicrobiol201533#Sec11 which \n",
    "Starting from table S2 in https://mbio.asm.org/content/9/6/e01780-18 which contains `run_accessions` for deposited raw sequencing data in ENA, use the ENA API to automatically add `study_accession` and `sample_accession`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "raven_2016_uk_hospital_metadata = pd.read_excel('datasets/41564_2016_BFnmicrobiol201533_MOESM202_ESM.xlsx', \n",
    "                                                skiprows=1,\n",
    "                                                skipfooter=4)\n",
    "\n",
    "# we are redoing genomics so just want phenotypic vanco resistance\n",
    "raven_2016_uk_hospital_metadata['Vancomycin'] = raven_2016_uk_hospital_metadata['Phenotypic vancomycin resistance'].replace({'VanA': 'R', 'VSE**': 'R', 'VanA***': 'R', 'VanB':'R', 'VSE': 'S', 'Unknown': np.nan}) \n",
    "raven_2016_uk_hospital_metadata['Location'] = raven_2016_uk_hospital_metadata['Hospital (HA), community (CA) or health-care associated (HCA)'].replace({'HA': 'Hospital', 'CA': 'Community', 'HCA': 'Health-Care Associated', 'Unknown': np.nan}) \n",
    "\n",
    "raven_2016_uk_hospital_metadata['Isolation_date'] = raven_2016_uk_hospital_metadata['Year of isolation'].astype(str).str.split('-').str.get(0)\n",
    "\n",
    "# drop uneeded columns (especially as we are doing our own genomics)\n",
    "raven_2016_uk_hospital_metadata = raven_2016_uk_hospital_metadata.drop(['Phenotypic vancomycin resistance', \n",
    "                                                                        'Hospital (HA), community (CA) or health-care associated (HCA)',\n",
    "                                                                        'Van transposon', 'Collection', 'Year of isolation'], axis=1)\n",
    "\n",
    "# rename columns to align better with gouil\n",
    "raven_2016_uk_hospital_metadata = raven_2016_uk_hospital_metadata.rename(columns={'Sequence ID': 'Run_accession',\n",
    "                                          'Source': 'Origin',\n",
    "                                          'ID in collection': 'Isolate_name',\n",
    "                                          'MLST type': 'Sequence_Type'})\n",
    "\n",
    "raven_2016_uk_hospital_metadata['Country/Province'] = \"United Kingdom\"\n",
    "raven_2016_uk_hospital_metadata['Species'] = \"Enterococcus faecalis\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# break into 500 accession chunks to more efficiently query the API\n",
    "api_responses = \"\"\n",
    "\n",
    "run_accessions = raven_2016_uk_hospital_metadata.loc[:,'Run_accession'].values\n",
    "for run_accession_chunk in [run_accessions[i:i + 500] for i in range(0, len(run_accessions), 500)]:\n",
    "    run_accession_chunk = \",\".join(run_accession_chunk)\n",
    "    api_response = requests.get(f\"https://www.ebi.ac.uk/ena/browser/api/xml/{run_accession_chunk}\")    \n",
    "    api_responses += api_response.text\n",
    "\n",
    "# parse combined xml records into a soup for easier traversal \n",
    "api_response_soup = bs4.BeautifulSoup(api_responses, 'lxml')\n",
    "\n",
    "\n",
    "def add_study_and_sample_metadata(row, api_response_soup):\n",
    "    \"\"\"\n",
    "    Use the ENA API to get the study and sample accessions\n",
    "    \"\"\"\n",
    "    run_accession = row['Run_accession']\n",
    "        \n",
    "    run_soup = api_response_soup.find(accession=run_accession)\n",
    "    if run_soup:\n",
    "        for xref in run_soup.find_all('xref_link'):\n",
    "            if xref.db.text == 'ENA-STUDY':\n",
    "                row['Study_accession'] = xref.id.text.strip()\n",
    "            elif xref.db.text == 'ENA-SAMPLE':\n",
    "                row['Sample_accession'] = xref.id.text.strip()\n",
    "    else:\n",
    "        row['Study_accession'] = np.nan\n",
    "        row['Sample_accession'] = np.nan\n",
    "    return row\n",
    "\n",
    "raven_2016_uk_hospital_metadata = raven_2016_uk_hospital_metadata.apply(lambda x: add_study_and_sample_metadata(x, api_response_soup), axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add read data for Raven 2016\n",
    "\n",
    "Data downloaded from ENA by FM, checksums verified (see `genomes/`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "raven_2016_uk_hospital_reads = get_reads('genomes/raven_2016_uk_hospital')\n",
    "raven_2016_uk_hospital_metadata_and_reads = pd.merge(raven_2016_uk_hospital_metadata, \n",
    "                                                     raven_2016_uk_hospital_reads, \n",
    "                                                     on='Run_accession', \n",
    "                                                     validate='one_to_one', \n",
    "                                                     how='outer')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Collate All UK Data\n",
    "\n",
    "All UK data can be collated into one table now that download is reliable instead of original protocol used by HS:\n",
    "1. get all accessions\n",
    "2. try to download\n",
    "3. keep the subset that successfully downloaded as a \"random sample with similar habitat distribution as all Alberta E. faecium metadata\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "uk_metadata_and_reads = pd.concat([gouliouris_2018_uk_metadata_and_reads, \n",
    "                                   raven_2016_uk_hospital_metadata_and_reads])\n",
    "uk_metadata_and_reads.loc[uk_metadata_and_reads['R1'].isna(), 'Read_status'] = \"Removed for deduplication in original paper\"\n",
    "uk_metadata_and_reads.loc[~uk_metadata_and_reads['R1'].isna(), 'Read_status'] = \"Validated download\"\n",
    "\n",
    "# # drop extraneous columns\n",
    "# uk_metadata_and_reads = uk_metadata_and_reads.drop(['Strain_name', 'Species_reads'], axis=1)\n",
    "\n",
    "# # make species information from metadata the official species for the UK data as this seems reliable for this dataset\n",
    "# uk_metadata_and_reads = uk_metadata_and_reads.rename(columns={'Species_metadata': 'Species'})\n",
    "\n",
    "# # re-order columns into natural groupings\n",
    "uk_metadata_and_reads = uk_metadata_and_reads[\n",
    "                      ['Isolate_name', 'Study_accession', 'Sample_accession', \n",
    "                       'Run_accession',  \n",
    "                       'Metadata_status',\n",
    "                       'Species', 'BAPS_group', 'Sequence_Type',\n",
    "                       'Country/Province', 'Origin', 'Location',\n",
    "                      'Ampicillin', 'Vancomycin', 'Read_status', \n",
    "                      'R1', \"R2\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# AAFC Enterococcus Paper Metadata\n",
    "\n",
    "## Collect metadata from NCBI\n",
    "\n",
    "https://www.nature.com/articles/s41598-020-61002-5\n",
    "\n",
    "Raw sequencing data is now deposited for this paper under `BioProject` accession (`PRJNA604849`) i.e., `study_accession`.\n",
    "\n",
    "So first let's get all the metadata from NCBI and link the reads before trying to merge in the metadata from the paper with variable names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('datasets/PRJNA604849_full_biosample_list.xml') as fh:\n",
    "    PRJNA604849_xml = bs4.BeautifulSoup(fh, 'lxml')\n",
    "\n",
    "parsed_xml_data = {}\n",
    "for biosample in PRJNA604849_xml.find_all('biosample'):\n",
    "    biosample_data = {}\n",
    "    biosample_data['Sample_name'] = biosample.find(db_label='Sample name').text\n",
    "    biosample_data['Species'] = biosample.find('organismname').text\n",
    "    biosample_data['Strain_name'] = biosample.find(attribute_name=\"strain\").text\n",
    "    biosample_data['Study_accession'] = biosample.find(target=\"bioproject\")['label']\n",
    "    parsed_xml_data[biosample['accession']] = biosample_data\n",
    "    \n",
    "alberta_ncbi = pd.DataFrame(parsed_xml_data).T.reset_index().rename(columns={'index': 'Sample_accession'})\n",
    "\n",
    "# tidy up as we don't use sample_name for matching because strain_name was a closer match to the isolate name\n",
    "# in the paper metadata\n",
    "alberta_ncbi = alberta_ncbi.drop('Sample_name', axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add read data to Alberta metadata from NCBI\n",
    "\n",
    "Despite reconciling this at various stages the constant changes to metadata has meant this pairing has broken again, I'm applying the same approach as used for merging in the paper metadata to NCBI again. \n",
    "\n",
    "Having persuaded AAFC to actually upload the reads I can link biosamples based on SRA accessions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "alberta_sra_data = pd.read_csv('datasets/PRJNA604849_sra_run_info.csv')\n",
    "alberta_sra_data = alberta_sra_data[['Run','BioSample']]\n",
    "alberta_sra_data = alberta_sra_data.rename(columns={'BioSample': 'Sample_accession', \n",
    "                                                    'Run': 'Run_accession'})\n",
    "alberta_read_data = get_reads('genomes/zaheer_2020_alberta')\n",
    "alberta_read_data = pd.merge(alberta_sra_data, \n",
    "                              alberta_read_data, \n",
    "                              on='Run_accession', \n",
    "                              validate='one_to_one', \n",
    "                              how='outer')\n",
    "\n",
    "alberta_ncbi_metadata_and_reads = pd.merge(alberta_ncbi,\n",
    "                                     alberta_read_data,\n",
    "                                     on='Sample_accession',\n",
    "                                     how='outer',\n",
    "                                     validate='one_to_many', suffixes=['', '_sra'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are duplicates where the same sample has multiple sequences (111) of them but there is no way to tell which is the resequencing with how AAFC uploaded the raw data as well as some samples with no reads associated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "alberta_ncbi_metadata_and_reads.loc[alberta_ncbi_metadata_and_reads['Sample_accession'].duplicated(), 'Metadata_status'] = 'Two sets of reads for this sample without time metadata to reconcile'\n",
    "alberta_ncbi_metadata_and_reads.loc[alberta_ncbi_metadata_and_reads['R1'].isna(), 'Read_status'] = 'No reads associated with sample'\n",
    "alberta_ncbi_metadata_and_reads.loc[~alberta_ncbi_metadata_and_reads['R1'].isna(), 'Read_status'] = 'Validated download'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merge in metadata from paper reconciling name problems\n",
    "\n",
    "Now we need to parse and tidy up the metadata from https://www.nature.com/articles/s41598-020-61002-5  `41598_2020_61002_MOESM2_ESM.csv` as this seems to have some disconnections with the metadata via NCBI e.g., species assignments we will use the NCBI data when in doubt. Given that NCBI confirm species assignments this seems a prudent choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "alberta_metadata = pd.read_csv('datasets/41598_2020_61002_MOESM2_ESM.csv', sep='\\t')\n",
    "\n",
    "# drop the 2 inexplicably duplicated rows \n",
    "alberta_metadata = alberta_metadata.drop_duplicates()\n",
    "\n",
    "# get rid of the identical rows apart from _ vs - \n",
    "alberta_metadata = alberta_metadata[alberta_metadata['ISOLATE'] != 'SWEntR-0393']\n",
    "# identical row apart from incorrect species (when compared to NCBI assembly)\n",
    "alberta_metadata = alberta_metadata[alberta_metadata['ISOLATE'] != 'SWEntR 0262']\n",
    "\n",
    "# tidy column names to help with merging later\n",
    "alberta_metadata = alberta_metadata.rename(columns={'ISOLATION SOURCE': 'Origin',\n",
    "                                                    'SPECIFIC LOCATION': 'Location',\n",
    "                                                    'VNCO': 'Vancomycin',\n",
    "                                                    'AMPI': 'Ampicillin',\n",
    "                                                    'TEIC': 'Teicoplanin',\n",
    "                                                    'DOXY': 'Doxycycline',\n",
    "                                                    'ERTH': 'Erythromycin',\n",
    "                                                    'GENT': 'Gentamicin',\n",
    "                                                    'LNZD': 'Linezolid',\n",
    "                                                    'LVFL': 'Levofloxacin',\n",
    "                                                    'QUIN': 'Quinolone',\n",
    "                                                    'STEP': 'Streptomycin',\n",
    "                                                    'NTRO': 'Nitrofurantoin',\n",
    "                                                    'TGC': 'Tigecycline',\n",
    "                                                    'SPECIATION': 'Species',\n",
    "                                                    'SOURCE CODE': 'Source_code',\n",
    "                                                    'ISOLATE': 'Isolate_name_paper'})\n",
    "\n",
    "# get rid of useless columns and add extra column to help with merging UK and AB metadata\n",
    "alberta_metadata = alberta_metadata.drop(['Unnamed: 18', \"Resistance count\", 'LOCATION'], axis=1)\n",
    "alberta_metadata['Country/Province'] = 'Canada/Alberta'\n",
    "\n",
    "# Remove trailing spaces that were left in the paper metadata\n",
    "alberta_metadata['Species'] = alberta_metadata['Species'].str.strip()\n",
    "\n",
    "# To ensure we merge the correct identifiers we are going to use the Source_code\n",
    "# information AS WELL as the isolate_name, therefore let's create a new identifier\n",
    "# out of the source code and isolate_name (and then remove spaces/underscores/hyphens etc)\n",
    "def combine_source_and_isolate_name(row):\n",
    "    \"\"\"\n",
    "    Try to combine isolate name with the source code field\n",
    "    if it isn't already prefixed by the source code information\n",
    "    \"\"\"\n",
    "    # Spaces, hyphens and dashes are a major source of disconnect so just remove them and try to map\n",
    "    metadata_isolate_name = row['Isolate_name_paper'].replace('_', '').replace('-', '').replace(' ', '')\n",
    "    if metadata_isolate_name.startswith(row['Source_code']):\n",
    "        return metadata_isolate_name\n",
    "    else:\n",
    "        return row['Source_code'] +  metadata_isolate_name\n",
    "\n",
    "\n",
    "alberta_metadata['Source_code_and_isolate_name'] = alberta_metadata.apply(\\\n",
    "                                                            combine_source_and_isolate_name, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to merge in the accession data from the bioproject.\n",
    "\n",
    "The NCBI biosamples have both a distinct `Strain_name` and a distinct `Sample_name` whereas the paper metadata just has an `Isolate_name_paper` (original `ISOLATE`). \n",
    "\n",
    "`Isolate_name_paper` seems to correspond to either one of these NCBI identifiers with no apparent pattern. However, `Strain_name` does match the paper metadata `Isolate_name_paper` more often and seems to be generally closer.\n",
    "\n",
    "Therefore, we are going to try and identify mappings between the NCBI `Strain_name` and the paper metadata `Isolate_name_paper`.  Then we are going to treat the NCBI `Strain_name` as the true `Isolate_name`.\n",
    "\n",
    "Fortunately, it is usually fairly obvious what the correct mapping is as they usually only differ in hyphens/dashes/spacing/number of leading 0's. \n",
    "The `Strain_name` in NCBI often contain part of what is the `Source_code` in the paper metadata (i.e., source), therefore I combined `Source_code` with `Isolate_name_paper` as `Source_code_and_isolate_name` and then searched for mappings using that.   \n",
    "\n",
    "This means I can match `Strain_name` to `Source_code_and_isolate_name` by finding the `Source_code_and_isolate_name` that shares the longest suffix with each `Strain_name` in NCBI.\n",
    "\n",
    "For extra security to prevent mis-assignments, I also added the following filter conditions: \n",
    "\n",
    "1. All the numerical portions of `Strain_name` and `Source_code_and_isolate_name` had to match to be valid\n",
    "2. The relationship between the two sets of names had to be one-to-one i.e., each `Strain_name` was assigned to one and only one `Source_code_and_isolate_name`.\n",
    "\n",
    "Any remaining `Strain_name` in the NCBI data that wasn't assigned a `Source_code_and_isolate_name` can then be manually reviewed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def attempt_to_reconcile_name_sets(set1, set2):\n",
    "    \"\"\"\n",
    "    Try and find the closest match between strings in set1 to set2\n",
    "    First find the closest matching string by longest shared suffix\n",
    "    Then as a safety move extract all digits and confirm they match between\n",
    "    the biosample identifier and the closest paper identifier\n",
    "    \n",
    "    returns: dictionary mapping likely paper identifiers from set1 to set2\n",
    "    \"\"\"  \n",
    "    closest_match_strings_by_suffix_length = []\n",
    "    # for each name in set1 clean it up then and revese it (to make the suffix a prefix)\n",
    "    for set1_name in set1:\n",
    "        set1_name_clean = set1_name.replace('-', '').replace('_', '').lower()[::-1]\n",
    "        \n",
    "        # compare to each cleaned name in set2\n",
    "        distances_between_set1_name_and_all_set2 = []\n",
    "        longest_suffix = {'suffix_length': 0, 'set2_name': '', 'set1_name': set1_name}\n",
    "        for set2_name in set2:\n",
    "            set2_name_clean = set2_name.replace('-', '').replace('_', '').lower()[::-1]\n",
    "            \n",
    "            # and recover the string that has the longest shared suffix \n",
    "            suffix_length = 0\n",
    "            for set1_name_char, set2_name_char in zip(set1_name_clean, set2_name_clean):\n",
    "                if set1_name_char == set2_name_char:\n",
    "                    suffix_length += 1\n",
    "                else:\n",
    "                    break\n",
    "                    \n",
    "            if suffix_length > longest_suffix['suffix_length']:\n",
    "                longest_suffix['suffix_length'] = suffix_length\n",
    "                longest_suffix['set2_name'] = set2_name\n",
    "\n",
    "            closest_match_strings_by_suffix_length.append(longest_suffix)\n",
    "    \n",
    "    # then check that the best matchs ALSO share all the same numerical components\n",
    "    closest_match_strings_by_suffix_size_and_numerical = []\n",
    "    for closest_pairing in closest_match_strings_by_suffix_length:\n",
    "        \n",
    "        # handle mapping long date names (failing due to numbers in the date portion)\n",
    "        if closest_pairing['set1_name'].startswith('ES-'):\n",
    "            set1_numerical = \"\".join(re.findall(r'\\d+', closest_pairing['set1_name'].split('-')[-1]))\n",
    "        else:\n",
    "            set1_numerical = \"\".join(re.findall(r'\\d+', closest_pairing['set1_name']))\n",
    "\n",
    "        if closest_pairing['set2_name'].startswith('ES-'):\n",
    "            set2_numerical =  \"\".join(re.findall(r'\\d+', closest_pairing['set2_name'].split('-')[-1]))\n",
    "        else:\n",
    "            set2_numerical =  \"\".join(re.findall(r'\\d+', closest_pairing['set2_name']))\n",
    "        \n",
    "        if set1_numerical == set2_numerical:\n",
    "            closest_match_strings_by_suffix_size_and_numerical.append((closest_pairing['set1_name'], \n",
    "                                                                       closest_pairing['set2_name']))\n",
    "        else:\n",
    "            # check if there is an issue with different numbers of leading 0s\n",
    "            # as this is a common issue\n",
    "            if set1_numerical.lstrip('0') == set2_numerical.lstrip('0'):\n",
    "                 closest_match_strings_by_suffix_size_and_numerical.append((closest_pairing['set1_name'], \n",
    "                                                                            closest_pairing['set2_name']))\n",
    "    \n",
    "    # create a dictionary from the closest hits and extract any in set1 without a mapping to set2\n",
    "    closest_match_strings_by_suffix_size_and_numerical = dict(closest_match_strings_by_suffix_size_and_numerical)\n",
    "    \n",
    "    # check for any non-one to one mappings in the dicionary i.e., different set1_name -> the same set2_name\n",
    "    # delete and manually resolve \n",
    "    set2_name_counter = Counter(closest_match_strings_by_suffix_size_and_numerical.values())\n",
    "    duplicates = []\n",
    "    closest_match_strings_by_suffix_size_and_numerical_no_duplicates = {}\n",
    "    for set1_name, set2_name in closest_match_strings_by_suffix_size_and_numerical.items():\n",
    "        # i.e. drop all those who don't have one-to-one mapping\n",
    "        if set2_name_counter[set2_name] == 1:\n",
    "            closest_match_strings_by_suffix_size_and_numerical_no_duplicates[set1_name] = set2_name\n",
    "    \n",
    "    set1_names_without_matches = set(set1) - set(closest_match_strings_by_suffix_size_and_numerical_no_duplicates.keys())\n",
    "    return closest_match_strings_by_suffix_size_and_numerical_no_duplicates, set1_names_without_matches\n",
    "\n",
    "alberta_ncbi_to_metadata_match, alberta_ncbi_without_metadata_match = attempt_to_reconcile_name_sets(alberta_ncbi_metadata_and_reads['Strain_name'].values, \n",
    "                                                                                                    alberta_metadata['Source_code_and_isolate_name'].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which NCBI strain_names didn't get a metadata match:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'CB_0150',\n",
       " 'CB_0182',\n",
       " 'CB_0383',\n",
       " 'ES-C-ST002-07DEC15-0142B',\n",
       " 'FC_0142B',\n",
       " 'HC_NS0026',\n",
       " 'HC_NS0078',\n",
       " 'HC_NS0150',\n",
       " 'HC_NS0238',\n",
       " 'HC_NS0383',\n",
       " 'HC_NS0854',\n",
       " 'HC_NS1042',\n",
       " 'HC_NS1090',\n",
       " 'HC_NS1104',\n",
       " 'HC_NS210',\n",
       " 'HC_SS0002',\n",
       " 'HC_SS0025',\n",
       " 'HC_SS0026',\n",
       " 'HC_VRE0078',\n",
       " 'SW_0002',\n",
       " 'SW_0025',\n",
       " 'SW_0182',\n",
       " 'SW_0238',\n",
       " 'WW_0060M',\n",
       " 'WW_0089I'}"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alberta_ncbi_without_metadata_match"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can try and manually fix these and identify cases of missing accessions or unresolvably ambiguous assignments (i.e., >1 perfectly valid appearing mapping from `Strain_name` to `Source_code_and_isolate_name` was possible)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "missing = {'HC_NS0026',  # ambiguous \n",
    "           'HC_NS0078', # ambiguous\n",
    "           'HC_NS0854', # missing\n",
    "           'HC_NS1042', # missing\n",
    "           'HC_NS1090', # missing\n",
    "           'HC_NS1104', # missing\n",
    "           'HC_VRE0078', # missing\n",
    "           'HC_SS0026', #missing\n",
    "           'WW_0060M', # missing\n",
    "           'WW_0089I'} # missing\n",
    "\n",
    "manual_fixes = {'CB_0150': 'CBEntR0150',\n",
    "                 'CB_0182': 'CBEntR0182',\n",
    "                 'CB_0383': 'CBSWEntR0383',\n",
    "                 'ES-C-ST002-07DEC15-0142B': 'WW0142B',\n",
    "                 'FC_0142B': 'FC0142B',\n",
    "                 'HC_NS0150': 'NSSNS0150',\n",
    "                 'HC_NS0238': 'NSSNS0238',\n",
    "                 'HC_NS0383': 'NSSNS0383',\n",
    "                 'HC_NS210': 'NSSNS0210',\n",
    "                 'HC_SS0002': 'SS0002',\n",
    "                 'HC_SS0025': 'SS0025',\n",
    "                 'SW_0002': 'NWSEnt0002',\n",
    "                 'SW_0025': 'NWSEnt0025',\n",
    "                 'SW_0182': 'NWSEnt0182',\n",
    "                 'SW_0238': 'NWSEnt0238'}\n",
    "\n",
    "alberta_ncbi_to_metadata_match.update(manual_fixes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately, this still leaves us with 10 isolates that are in the deposited NCBI data but don't seem to have any supplied metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "unresolved_ncbi = alberta_ncbi_metadata_and_reads[~alberta_ncbi_metadata_and_reads['Strain_name'].isin(alberta_ncbi_to_metadata_match.keys())].sort_values('Strain_name')['Strain_name'].values\n",
    "# add a status to the NCBI data indicating missing\n",
    "alberta_ncbi_metadata_and_reads.loc[alberta_ncbi_metadata_and_reads['Strain_name'].isin(unresolved_ncbi), 'Metadata_status'] = 'No linkable metadata in original paper (10.1038/s41598-020-61002-5)'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using our mapping from `Strain_name` to `Source_code_and_isolate_name` let's add a `Source_code_and_isolate_name` to the NCBI metadata sheet and merge using that (ensuring valid one-to-one merging).  \n",
    "\n",
    "We are doing a `left` merge because we only want to keep the metadata from the paper for isolates that correspond to the NCBI bioproject."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# translate the incorrect isolate_names in the paper metadata to the correct ones deposited in NCBI\n",
    "alberta_ncbi_metadata_and_reads['Source_code_and_isolate_name'] = alberta_ncbi_metadata_and_reads['Strain_name'].apply(lambda x: alberta_ncbi_to_metadata_match[x] if x in alberta_ncbi_to_metadata_match else f\"UNMATCHED: {x}\")\n",
    "\n",
    "# add the 10 unresolved identifiers to the paper metadata before merging\n",
    "alberta_metadata = pd.concat([alberta_metadata, pd.DataFrame({'Source_code_and_isolate_name': unresolved_ncbi})])\n",
    "\n",
    "alberta_merged = pd.merge(alberta_ncbi_metadata_and_reads, alberta_metadata, \n",
    "                          how='left',\n",
    "                          on='Source_code_and_isolate_name', suffixes=['_ncbi', '_paper'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we want to tidy up the alberta metadata to make our ultimate merging with the UK data easier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rename NCBI species information as official species information\n",
    "# and change the NCBI `Strain_name` to `Isolate_name` to match better with other datasets (even though what\n",
    "# NCBI calls things should remain the master)\n",
    "alberta_merged = alberta_merged.rename(columns={'Species_ncbi': 'Species',\n",
    "                                                'Strain_name': 'Isolate_name'})\n",
    "\n",
    "# drop the superflous and likely erroneous species information from paper metadata sheet\n",
    "alberta_merged = alberta_merged.drop(['Species_paper', 'Source_code_and_isolate_name'], axis=1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Final tidy up before merging of UK and AAFC data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Merging the two datasets\n",
    "\n",
    "Combine the two datasets and tidy the names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_data = pd.concat([alberta_merged.reset_index(drop=True), uk_metadata_and_reads.reset_index(drop=True)], \n",
    "                     join='outer', ignore_index=True)\n",
    "\n",
    "all_data = all_data[['Study_accession', 'Sample_accession', 'Run_accession', 'Isolate_name',\n",
    "                     'Isolate_name_paper', \n",
    "                     'Metadata_status', 'Species', 'BAPS_group','Sequence_Type', 'Country/Province', \n",
    "                     'Origin', 'Location', 'Source_code', \n",
    "                     'Ampicillin', 'Vancomycin', 'Teicoplanin', 'Doxycycline',\n",
    "                     'Erythromycin', 'Nitrofurantoin', 'Gentamicin', 'Linezolid',\n",
    "                     'Levofloxacin', 'Quinolone', 'Streptomycin', 'Tigecycline', 'Read_status',\n",
    "                     'R1', 'R2']]\n",
    "\n",
    "all_data.to_csv('all_combined_enterococcus_metadata.tsv', index=False, sep='\\t')"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}