SH2_mutation_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from CoDIAC import jalviewFunctions, featureTools\n",
    "import pandas as pd\n",
    "\n",
    "# Number of SH2 domain positions that must have a feature to be considered conserved\n",
    "LIGAND_CONSERVATION = 8 \n",
    "\n",
    "# These aren't currently being used, but these might be considered as the cutoff for conserved, based on total numbers\n",
    "PHOSPHOTYROSINE = 5\n",
    "PHOSPHOSERINE = 5\n",
    "UBIQUITINATION = 5\n",
    "PHOSPHOTHREONINE = 3\n",
    "ACETYLATION = 3\n",
    "\n",
    "\n",
    "aln_file = 'Data/Uniprot_Reference/alignment/SH2_IPR000980_promals3d.fasta'\n",
    "aln_mapping_dict, aln = featureTools.return_position_mapping_dict(aln_file)\n",
    "\n",
    "#load the mutations\n",
    "mut_file = 'Data/Features/mutations/SH2_IPR000980_mutations.feature'\n",
    "feature_colors, mut_dict = jalviewFunctions.return_features_from_file(mut_file)\n",
    "mut_ann = featureTools.return_feature_ann_dict(aln, mut_dict)\n",
    "\n",
    "# Load the domain domain interactions \n",
    "domain_file = 'Data/Features/interDomain/SH2_IPR000980_interDomain_phospholipid.feature' #will treat phospholipid as a feature like domain interaction interface for the purposes of interesecting\n",
    "feature_colors, domain_int_dict = jalviewFunctions.return_features_from_file(domain_file)\n",
    "\n",
    "# Will also compare to the ligand binding sites\n",
    "ligand_file = 'Data/Features/ligand/SH2_IPR000980_Ligand.feature'\n",
    "feature_colors, ligand_dict = jalviewFunctions.return_features_from_file(ligand_file)\n",
    "ligand_ann = featureTools.return_feature_ann_dict(aln, ligand_dict)\n",
    "ligand_ann_conserved = ligand_ann['Ligand'].sum() >= LIGAND_CONSERVATION\n",
    "\n",
    "\n",
    "#load all of the PTM features\n",
    "ptm_all_file = 'Data/Features/PTMs_all/PTMs_final/SH2_IPR000980_PTMs_all.feature'\n",
    "feature_colors, ptm_all_dict = jalviewFunctions.return_features_from_file(ptm_all_file)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Produce a mutation-centric information table, based on contacts of ligand and domain-domain interfaces and PTMs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#these are direct intersections between PTMs of type mod and domain and ligand binding sites\n",
    "feature_intersection_domain_dict, headers_no_contacts = featureTools.return_intersection(mut_dict, domain_int_dict)\n",
    "feature_intersection_ligand_dict, headers_no_contacts = featureTools.return_intersection(mut_dict, ligand_dict)\n",
    "feature_intersection_PTMs_dict, headers_no_contacts = featureTools.return_intersection(mut_dict, ptm_all_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = create_report_for_mut(mut_dict, mut_ann['Mutation'], aln_mapping_dict, ligand_ann_conserved, feature_intersection_domain_dict, feature_intersection_PTMs_dict)\n",
    "df.to_csv('Data/Reports/SH2_IPR000980_mutation_summary.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write a readme file that documents what we did, what it means, and what the cutoffs are\n",
    "from datetime import datetime\n",
    "\n",
    "# current date and time\n",
    "now = datetime.now()\n",
    "\n",
    "\n",
    "\n",
    "readme_file = 'Data/Reports/FeatureAnalysis_mutations.md'\n",
    "f = open(readme_file, 'w')\n",
    "f.write('## Feature Analysis\\n')\n",
    "f.write('This document describes the features that were analyzed for the SH2 domain as reported in this directory for mutations\\n')\n",
    "f.write('These reports were generated by the CoDIAC package on: %s\\n'%now.strftime(\"%m/%d/%Y %H:%M:%S\"))\n",
    "f.write('### Analysis Reports (CSV files)\\n')\n",
    "f.write('Mutation Features were analyzed with regard to:\\n')\n",
    "f.write('1. The proximity to a conserved ligand binding sites\\n')\n",
    "f.write('2. The direct intersection with domain-domain interactions\\n')\n",
    "f.write('3. The number of other mutations of that type in the alignment\\n')\n",
    "f.write('4. The direct intersection with a known PTM\\n')\n",
    "f.write('### Feature Cutoffs and files used\\n')\n",
    "f.write('Ligand binding was considered conserved if >= %d residues in that alignment position were in the ligand binding features\\n'%(LIGAND_CONSERVATION))\n",
    "f.write('The following files were used:\\n')\n",
    "f.write('1. %s\\n'%aln_file)\n",
    "f.write('2. %s\\n'%domain_file)\n",
    "f.write('3. %s\\n'%ligand_file)\n",
    "f.write('4. %s\\n'%ptm_all_file)\n",
    "f.write('5. Mutations %s\\n'%mut_file)\n",
    "f.write('### Analysis Report Columns\\n')\n",
    "f.write('1. header: the fasta header of the sequence\\n')\n",
    "f.write('2. Uniprot_ID: the uniprot ID of the sequence\\n')\n",
    "f.write('3. gene_name: the gene name of the sequence\\n')\n",
    "f.write('4. SH2 number: the SH2 domain number (1 if N-terminal and 2 if C-terminal domain in 2 SH2 domain proteins\\n')\n",
    "f.write('5. protein_residue: the residue in the protein sequence numbering\\n')\n",
    "f.write('6. feature_position: the position in the domain, as printed in the feature file, relative to SH2 start\\n')\n",
    "f.write('7. alignment_postion: the position in the alignment\\n')\n",
    "f.write('8. number_mods_at_position: the number of modifications of this type at that position\\n')\n",
    "f.write('9. number_conserved_ligand_residues_nearby: the number of residues in the ligand binding site at that position within +/- 2 amino acids away\\n')\n",
    "f.write('10. number_residues_from_conserved_ligand: the number of residues from the conserved ligand binding site. If 0, it falls directly on a ligand interface.\\n')\n",
    "f.write('11. domain_interaction: the domain interaction at that position, if it intersects directly with a domain-domain interface. We also use this column to indicate if it interacts with a residue identified as interacting with phospholipids\\n')\n",
    "f.write('12. PTM: the modification type(s) at that position, if it intersects directly with a PTM\\n')\n",
    "f.close()\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#def create_report_for_mod(mod, ptm_dict, ptm_annotations_dict, aln_mapping_dict, ligand_ann_conserved, feature_intersection_domain_dict):\n",
    "def create_report_for_mut(mut_dict, mut_ann, aln_mapping_dict, ligand_ann_conserved, feature_intersection_domain_dict, feature_intersection_PTMs_dict):\n",
    "\n",
    "    ptm_features = mut_dict\n",
    "    ptm_annotations = mut_ann\n",
    "    sum_ann = ptm_annotations.sum()\n",
    "    number_dict = {}\n",
    "    ligand_conserved_window = {}\n",
    "    ligand_distance_dict = {}\n",
    "\n",
    "    # Collecting information on PTMs and ligand binding sites\n",
    "    for header in ptm_features:\n",
    "        number_dict[header] = {}\n",
    "        ligand_conserved_window[header] = {}\n",
    "        ligand_distance_dict[header] = {}\n",
    "        aln_map = aln_mapping_dict[header]\n",
    "        for pos in ptm_features[header]:\n",
    "            #print(\"DEBUG: checking %d for %s\"%(int(pos), header))\n",
    "            pos_index = int(pos)-1 #have to translate to 0-based index\n",
    "            aln_pos = aln_map[pos_index] #this is the position in the alignment, 0-based\n",
    "            ptm_number = sum_ann[int(aln_pos)]\n",
    "            number_dict[header][pos] = {'aln_pos':aln_pos+1, 'num_mutations':ptm_number}\n",
    "            ligand_conserved_window[header][pos] = ligand_ann_conserved[aln_pos-2:aln_pos+2].sum() #total number of conserved interaction residues in a window of 5 residues\n",
    "            #let's record the furthest distance it is to a conserved ligand binding sit\n",
    "            ligand_distance_dict[header][pos] = min(abs(aln_pos - np.where(ligand_ann_conserved)[0]))\n",
    "        \n",
    "    df = pd.DataFrame(columns=['header', 'Uniprot_ID', 'gene_name', 'SH2 number', 'protein_residue', 'feature_position', \n",
    "                           'alignment_postion', 'number_mutations_at_position', 'number_conserved_ligand_residues_nearby', 'number_residues_from_conserved_ligand', 'domain_interaction', 'PTMs'])\n",
    "    for header in ptm_features:\n",
    "        uniprot_id, name, number, start, stop = header.split('|')\n",
    "        for pos in ptm_features[header]:\n",
    "            arr_temp = feature_intersection_domain_dict[header][pos]\n",
    "            string_domain = ';'.join(arr_temp)\n",
    "            arr_temp = feature_intersection_PTMs_dict[header][pos]\n",
    "            string_PTM = ';'.join(arr_temp)\n",
    "            dict_temp = {'header':header, 'Uniprot_ID':uniprot_id, 'gene_name':name, 'SH2 number':number, \n",
    "                            'protein_residue':int(start)+int(pos)-1, 'feature_position':pos, 'alignment_postion':number_dict[header][pos]['aln_pos'], 'number_mutations_at_position':number_dict[header][pos]['num_mutations'], \n",
    "                            'number_conserved_ligand_residues_nearby':ligand_conserved_window[header][pos], 'number_residues_from_conserved_ligand':ligand_distance_dict[header][pos], 'domain_interaction':string_domain, 'PTMs':string_PTM}\n",
    "            df_temp = pd.DataFrame.from_dict(dict_temp, orient='index').T\n",
    "            df = pd.concat([df, df_temp])\n",
    "    return df"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "CoDIAC-dev",
   "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}