Introduction_to_python_week_4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An introduction to solving biological problems with Python\n",
    "\n",
    "## Week 4: Use Cases/Practice\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4.1.1  Filtering TSV Files\n",
    "\n",
    "Work in pairs.\n",
    "\n",
    "Use Case: You need to filter a large file of variants that's too big to open in Excel or R.\n",
    "\n",
    "Write a script that reads in the variant data/test1_example1.tsv file and writes out just the Passed variants that have quality >= 50 to a new tab-separated file. \n",
    "Hint: use the Filter and Qual columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use Case: Combining Multiple Files\n",
    "\n",
    "You need to combine variants from a set of tabular patient files into a single file. You only want variants that have a variant allele frequency (vaf) > 10. But an additional complication is that you want to include variants in the PIK3CA gene that have a vaf > 5 and variants in the ESR1 gene with vaf > 2. The script below shows how the multiple patient files can be combined applying these different filter thresholds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the libraries that we need\n",
    "from glob import glob\n",
    "from csv import DictReader\n",
    "\n",
    "patientIDs = [\"01\", \"10\", \"20\"]\n",
    "dirpath = \"data/patients/variants/\"\n",
    "\n",
    "# Loop through the patient IDs\n",
    "for id in patientIDs:\n",
    "    # Generate the search path\n",
    "    # Note that \"*\" stand for a wildcard (AKA anything)\n",
    "    path = \"\".join([dirpath, \"*\", id, \".csv\"])\n",
    "#     print(path) # The \"search\" term we are using for glob\n",
    "    filespath = glob(path)\n",
    "#     print(filespath) # Prints the full path name to the files\n",
    "    \n",
    "    data = dict() # Initialise an empty dictionary\n",
    "    \n",
    "    # Loop through the files relevant to each patient\n",
    "    for file in filespath:\n",
    "        print(\"Filename:\", file)\n",
    "        with open(file, \"r\") as infile:\n",
    "            fh = DictReader(infile) # fh is the file handle\n",
    "            tmp = next(fh) # Grab the first line to help initialize the dictionary\n",
    "            for k,v in tmp.items():\n",
    "                data.update({k: [v]})\n",
    "            # Iterate through the rest of the file and add in the values to the corresponding key\n",
    "            for line in fh:\n",
    "                for k,v in line.items():\n",
    "                    data[k].append(v)\n",
    "    keep = list() # Initialise an empty list for the indexes of the rows we want to keep\n",
    "    for i in range(len(data[\"gene\"])):\n",
    "        if (data[\"gene\"][i] == \"PIK3CA\" and float(data[\"vaf\"][i]) > 5):\n",
    "            keep.append(i)\n",
    "        elif (data[\"gene\"][i] == \"ESR1\" and float(data[\"vaf\"][i]) > 2):\n",
    "            keep.append(i)\n",
    "        elif (float(data[\"vaf\"][i]) > 10):\n",
    "            keep.append(i)\n",
    "    \n",
    "    # Initialise an empty dictionary to store the results\n",
    "    results = dict()\n",
    "    # Note here the key is the column name and the value is a list\n",
    "    for key, value in data.items(): # Iterate through each key, value\n",
    "        for i in keep: # Iterate through the list of indexes that we want to keep\n",
    "            if key not in results: # Check if the column name is in the dictionary, if not:\n",
    "                results[key] = [value[i]] # Initialise the column name as a key and set the value to be a list of the value\n",
    "            else:\n",
    "                results[key].append(value[i]) # If the column name is already in the dictionary then add to the list\n",
    "    \n",
    "    # Write out the results dictionary to a file\n",
    "    with open(\"\".join([\"patient_\", id, \"_allFiltered.csv\"]), \"w\") as outfile:\n",
    "        outfile.write(\",\".join([\"patientID\", \"gene\", \"vaf\"]) + \"\\n\") # Write out the header\n",
    "        lines = list(zip(results[\"patientID\"], results[\"gene\"], results[\"vaf\"])) # See example below for how zip works\n",
    "        # Write out each row to the output file\n",
    "        for line in lines:\n",
    "            outfile.write(\",\".join(list(line)) + \"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The zip function can be used to pair items, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(zip([\"a\", \"b\", \"c\"], [\"1\", \"2\", \"3\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4.2.1 Combining Multiple Files\n",
    "\n",
    "Work in pairs.\n",
    "\n",
    "Use case: You have 2000 files from TCGA with gene expression counts for patients (1 file per patient). You need to combine them into a single table of counts with a column for each patient. \n",
    "\n",
    "The count files are in data/patients/counts.zip (you will need to unzip the file) and the expected output is in the file data/patients/counts_output/combined_counts.tsv"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}