COSMIC_Data_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# COSMIC DATA ANALYSIS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`COSMIC_csv_reader.py\n",
    "COSMIC_sample_printer.py\n",
    "pyfaidx_extraction.py\n",
    "context_finder.py\n",
    "threemer_count_jellyfish.sh\n",
    "create_spectra.sh\n",
    "create_uhc_dendrogram_heatmap.py`<br>\n",
    "<br>\n",
    "The objective of this project was to analyze the sequence contexts of point mutation data downloaded from the Catalogue of Somatic Mutations in Cancer (COSMIC) database (https://cancer.sanger.ac.uk/cosmic/download) and compare it to sequence context data from UVB exposure experiments conducted in the Essigmann lab."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Descriptions and Methodology"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filtering and outputting COSMIC point mutation data into separate .csv files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`COSMIC_csv_reader.py`<br>\n",
    "Reads the downloaded .csv file from COSMIC, V85 and filters by certain parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if row[0] == \"GENE_NAME\": continue # Skips the header of the .csv\n",
    "if row[PRIMARY_SITE] != \"skin\": continue\n",
    "if row[PRIMARY_HISTOLOGY] != \"malignant_melanoma\": continue\n",
    "if \"del\" in row[MUTATION_CDS] or \"ins\" in row[MUTATION_CDS]: continue # Ignore indels\n",
    "if row[GRCH] == \"null\" or int(row[GRCH]) != 38: continue\n",
    "if \"null\" in row[MUTATION_STRAND]: continue\n",
    "if row[MUTATION_SOMATIC_STATUS] != \"Confirmed somatic variant\": continue"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It then creates .csv files for each sample group in a new directory and outputs the following information into the following format:\n",
    "\n",
    "| Position    | Sample Name | Strand:Substitution | Count |\n",
    "|:-----------:|:-----------:|:-------------------:|:-----:|\n",
    "| 12:71994489 | 1_RESISTANT | +:A>G               |   2   |\n",
    "| ...         | ...         | ...                 | ...   |\n",
    "\n",
    "Notes:\n",
    "- Position indicates the genomic position, in the Chromosome number:Base position format.\n",
    "- A '+' strand indicates that the mutation position is located on the forward strand.\n",
    "- A '-' strand indicates it is located on the other strand, thus the complementary base must be searched (In this case: T) and the reverse complement sequence must be written to output (ex: ATA where T is the lesion site).\n",
    "- Substitution is formatted as original base>mutation.\n",
    "- The count of 2 in the example represents that a A>G mutation occurred twice at 12:71994489 for the 1_RESISTANT sample group.\n",
    "- Any tandem base substitutions that are read are split up into individual entries occurring at their respective positions.\n",
    "\n",
    "Execute the script using the following line as an example: `python COSMIC_csv_reader.py -i V85_38_MUTANT.csv -o ./samples/`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`COSMIC_sample_printer.py`<br>\n",
    "A small portion of COSMIC_csv_reader.py that filters by certain parameters and prints the samples to console. Used primarily for searching for details of particular samples.<br>\n",
    "<br>\n",
    "Execute the script using the following line as an example: `python COSMIC_sample_printer.py -i V85_38_MUTANT.csv\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finding sequence contexts within the human genome"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`context_finder.py`<br>\n",
    "Reads the .csv files outputted by COSMIC_csv_reader.py and uses pyfaidx to locate the sequence context of the mutation position. The program calculates and prints the percent of successful matches of the original base to genomic position. As the program counts the number of occurrences of a sequence context occurring at a certain genomic position within an individual sample, a great amount of memory is required and it is recommended that the script run with about 4 GB of memory. The sequence contexts for each sample are outputted into .csv files within a new directory. If a base does not match to the position, it is not written to the outputted to its respective .csv file. The format of each outputted file is different than those outputted by COSMIC_csv_reader.py to allow create_spectra.sh (specifically mutpos_manual.py) to read and create spectra.<br>\n",
    "<br>\n",
    "\n",
    "The format of the output .csv files should look like this:\n",
    "\n",
    "| Position    | Context | Substitution | Count |\n",
    "|:-----------:|:-------:|:------------:|:-----:|\n",
    "| 12:71994489 | CAG     | A>G          | 2     |\n",
    "| ...         | ...     | ...          | ...   |\n",
    "\n",
    "Execute the script using the following line as an example: `python context_finder.py -f hg38.fa -i ./samples/ -o ./samples_contexted/`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Copying chromosome sequences to a .fasta file for threemer count comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`pyfaidx_extraction.py`<br>\n",
    "A script that uses pyfaidx to quickly read chromosome sequences (>chr1, >chr2, ..., >chrX, >chrY) from the .fasta file containing the human genome, hg38.fa, and output them to a new .fasta file called hg38_chromosomes.fa. This script was made to see if the outputted reference counts .txt of the threemer_count_jellyfish.sh script was impacted by .fasta entries within hg38.fa that were not relevant to the point mutation data. While there is some difference, proportionality for each threemer count was maintained.<br>\n",
    "<br>\n",
    "Execute the script using the following line: `python pyfaidx_extraction.py`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Counting of threemers in the human genome"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`threemer_count_jellyfish.sh`<br>\n",
    "Uses jellyfish to quickly count the number of threemers within the human genome (hg38.fa or hg38_chromosomes.fa) and outputs the counts into mer_counts.txt file. The .txt file is formatted into two tab-separated columns so that it can be read by the create_spectra.sh script. Jellyfish will also output a histogram (which I admittedly do not know how to open). The script was run for both hg38.fa and hg38_chromosomes.fa. Although the sequence counts were all higher for hg38.fa, the difference between counts for both were proportional.<br> \n",
    "<br>\n",
    "Run this script on ROUS using the following line as an example: `qsub threemer_count_jellyfish hg38.fa`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creation of Spectra for Each Sample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`create_spectra.sh`<br>\n",
    "Counts the number of .csv entries (unique mutations) within each contexted sample .csv file and only creates spectra for sample .csv files that contain 50 or more entries. The header of the .csv file does not count as an entry. The spectra are specifically created by the lab's in-house script, mutpos_manual.py. For each sample group with >= 50 unique mutations, the script will output a mutation spectrum and .csv file into a new directory.\n",
    "\n",
    "Execute the script using the following line as an example: `./create_spectra.sh ./samples_contexted ../ref/hg38_chromosomes_ref_counts.txt ./sample_contexted_figures`<br>\n",
    "Where the format is ./create_spectra.sh {input contexted .csv dir} {reference count .txt made from threemer_count_jellyfish.sh} {output spectra and .csv files}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creation of UHC Dendrograms and Heatmaps for Visual Comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`create_uhc_dendrogram_heatmap.py`<br>\n",
    "Performs unsupervised hierarchial clustering upon .csv sample files placed in a directory. The script uses functions from the in-house script, cospec.py, to create a dendrogram and heatmap of cosine similarity among samples. For the heatmap, it is intended to take an input of 10 or so samples. For the dendrogram, it is recommended to change the function in cospec.py to match the following code if you are working with hundreds of samples:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_uhc_dendrogram(spec_list, cluster_names):\n",
    "    spectra = []\n",
    "    for sig in spec_list:\n",
    "        spectra.append(list(spec_list[sig].values()))\n",
    "    data = spectra\n",
    "    labels = cluster_names\n",
    "    linkages = uhc_cluster(spec_list)\n",
    "    with plt.rc_context({'lines.linewidth':0.1}): #Helpeful to make the tree more interpretable\n",
    "        dendrogram(linkage(pd.DataFrame(data, [labels, ['']*len(labels)]),\n",
    "                                 method='ward',\n",
    "                                 metric='cosine'),\n",
    "                         no_labels=False,\n",
    "                         labels=labels,\n",
    "                         leaf_rotation=90,\n",
    "                         leaf_font_size=1, #Needs to be very small to make labels readable\n",
    "                         link_color_func=lambda x:'k')\n",
    "    plt.show()\n",
    "    plt.savefig('file.pdf', dpi=450) #File is saved as file.pdf in the same directory containing cospec.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Execute the script using the following line as an example: `python create_uhc_dendrogram_heatmap.py -i ../../muts_csv/`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Author"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scripts `cospec.py` and `mutpos_manual.py` were written by klkim [at] mit.edu<br>\n",
    "All other scripts featured were written by jggatter [at] mit.edu<br>\n",
    "pyfaidx (https://pythonhosted.org/pyfaidx/) <br>\n",
    "jellyfish (http://www.cbcb.umd.edu/software/jellyfish/)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "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
}