MRE_candidate_gene_discovery_example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "912ec80f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#imports\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn import metrics\n",
    "from sklearn.cluster import DBSCAN\n",
    "import re\n",
    "pd.set_option('display.max_columns', None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b45c7443",
   "metadata": {},
   "outputs": [],
   "source": [
    "#list of all MRE types\n",
    "MREs_list = [\"TGCGCAC\", \"TGCGCGC\", \"TGCGCCC\", \"TGCGCTC\", \"TGCACAC\", \"TGCACGC\", \"TGCACCC\", \"TGCACTC\", \n",
    "             \"GAGTGCA\", \"GGGTGCA\", \"GCGTGCA\", \"GTGTGCA\", \"GAGCGCA\", \"GGGCGCA\", \"GCGCGCA\", \"GTGCGCA\"]\n",
    "\n",
    "#initiate final DataFrame\n",
    "overlapping_ranges_all = pd.DataFrame(columns=[\"df1_row\", \"df2_row\", \"Cluster\", \"MRE_start\",\"MRE_type\",\"MRE_zone_start\",\"MRE_zone_end\",\"MRE_type_str\",\"TrName\",\"Chr\",\"Source\",\"Feature\",\"Start\",\"End\",\"Score\",\"Strand\",\"Frame\",\"Attribute\"])\n",
    "\n",
    "#loop through chromosome FASTAs to discover Metal Responsive Element clusters near genes\n",
    "#filter by MRE cluster types found in the Metallothionein promoter\n",
    "for i in range(1,19):\n",
    "    MRECW_df = pd.DataFrame(columns=[\"MRE_start\", \"MRE_type\", \"Chromosome\"])\n",
    "    \n",
    "    #import gtf of gene predictions for the respective chromosome\n",
    "    df_chr1_gtf = pd.read_csv(f\"./chr_gtf/chr{i}.gtf\", delimiter=\"\\t\", header=None)\n",
    "    df_chr1_names = pd.read_csv(f\"./chr_gtf/uniq_CDS{i}.txt\", delimiter=\"\\t\", header=None)\n",
    "\n",
    "    df_chr1_gtf.columns=['Chr', 'Source', 'Feature', 'Start', 'End', 'Score', 'Strand', 'Frame', 'Attribute']\n",
    "    df_chr1_gtf.insert(0, \"TrName\", df_chr1_names)\n",
    "    \n",
    "    #import chromosome FASTA for the respective chromosme\n",
    "    chr17 = pd.read_csv(f\"./chr_fasta/chr{i}.fa\", names=[\"seq\", \"score\"])\n",
    "    chr17_df = pd.DataFrame([*chr17[\"seq\"][1]], columns=[\"Nucleotide\"])\n",
    "    chr17_df[\"Coordinate\"] = chr17_df.index\n",
    "    chr17_df[\"Score\"] = 1\n",
    "    chr17_df[\"Coordinate\"] = chr17_df[\"Coordinate\"] + 1\n",
    "    \n",
    "    MRE_df = pd.DataFrame(columns=[\"MRE_start\", \"MRE_type\"])\n",
    "    \n",
    "    #find all MRE sites in the chromosome\n",
    "    for MRE in MREs_list:\n",
    "        all_MREs = [m.start() for m in re.finditer(MRE, chr17[\"seq\"][1])]\n",
    "        MRE_df_temp = pd.DataFrame(all_MREs, columns=[\"MRE_start\"])\n",
    "        MRE_df_temp[\"MRE_type\"] = MRE\n",
    "        MRE_df = pd.concat([MRE_df, MRE_df_temp], axis=0)\n",
    "    \n",
    "    #caclulate MRE counts\n",
    "    MRE_freq = []\n",
    "    MRE_frequencies = pd.DataFrame()\n",
    "    for MRE in MREs_list:\n",
    "        MRE_freq.append(MRE_df[MRE_df[\"MRE_type\"]==MRE].shape[0])\n",
    "    MRE_frequencies[\"MRE_type\"] = MREs_list\n",
    "    MRE_frequencies[\"Occurrences\"] = MRE_freq\n",
    "\n",
    "    MRE_df_reset_idx = MRE_df.sort_values(by=\"MRE_start\").reset_index(drop=True)\n",
    "    \n",
    "    #cluster MREs within 350bp of one another \n",
    "    db = DBSCAN(eps=350, min_samples=2).fit(MRE_df_reset_idx[\"MRE_start\"].values.reshape(-1,1))\n",
    "    labels = db.labels_\n",
    "    n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)\n",
    "    n_noise_ = list(labels).count(-1)\n",
    "    \n",
    "    #filter clusters to have exactly 3 MREs only\n",
    "    labels = db.labels_\n",
    "    n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)\n",
    "    n_noise_ = list(labels).count(-1)\n",
    "    MRE_df_reset_idx[\"Cluster\"] = labels\n",
    "    clustered_MREs_df = MRE_df_reset_idx[MRE_df_reset_idx[\"Cluster\"] != -1]\n",
    "    clustered_MREs_df.value_counts(subset=\"Cluster\")\n",
    "    counts = clustered_MREs_df[\"Cluster\"].value_counts(dropna=False)\n",
    "    clustered_MREs_df_over2 = clustered_MREs_df[clustered_MREs_df[\"Cluster\"].map(counts)>2]\n",
    "    clustered_MREs_df_over2_below4 = clustered_MREs_df_over2[clustered_MREs_df_over2[\"Cluster\"].map(counts)<4]\n",
    "    \n",
    "    #sort data structures\n",
    "    df1 = clustered_MREs_df_over2_below4.groupby([\"Cluster\"])[\"MRE_type\"].apply(list).reset_index()\n",
    "    clustered_MREs_df_over2_below4 = clustered_MREs_df_over2_below4.groupby([\"Cluster\"])[\"MRE_start\"].apply(list).reset_index()\n",
    "    clustered_MREs_df_over2_below4[\"MRE_type\"] = df1[\"MRE_type\"]\n",
    "    clustered_MREs_df_over2_below4[\"MRE_type_str\"] = clustered_MREs_df_over2_below4[\"MRE_type\"].astype(str)\n",
    "    clustered_MREs_df_over2_below4[\"MRE_zone_start\"] = clustered_MREs_df_over2_below4[\"MRE_start\"].apply(min)\n",
    "    clustered_MREs_df_over2_below4[\"MRE_zone_end\"] = clustered_MREs_df_over2_below4[\"MRE_start\"].apply(max)\n",
    "    \n",
    "    #filter MREs to be identical to the ones found in the Metallothionein promoter\n",
    "    #chr1:\n",
    "    df1 = clustered_MREs_df_over2_below4[clustered_MREs_df_over2_below4[\"MRE_type_str\"]==\"['GAGCGCA', 'TGCACGC', 'GTGTGCA']\"]\n",
    "    df1_rev = clustered_MREs_df_over2_below4[clustered_MREs_df_over2_below4[\"MRE_type_str\"]==\"['GTGTGCA', 'TGCACGC', 'GAGCGCA']\"]\n",
    "    df2_rc = clustered_MREs_df_over2_below4[clustered_MREs_df_over2_below4[\"MRE_type_str\"]==\"['TGCGCTC', 'GCGTGCA', 'TGCACAC']\"]\n",
    "    df2_rc_rev = clustered_MREs_df_over2_below4[clustered_MREs_df_over2_below4[\"MRE_type_str\"]==\"['TGCACAC', 'GCGTGCA', 'TGCGCTC']\"]\n",
    "    #chr4:\n",
    "    df3 = clustered_MREs_df_over2_below4[clustered_MREs_df_over2_below4[\"MRE_type_str\"]==\"['GCGTGCA', 'TGCACAC', 'GCGTGCA']\"]\n",
    "    df3_rev = clustered_MREs_df_over2_below4[clustered_MREs_df_over2_below4[\"MRE_type_str\"]==\"['GCGTGCA', 'TGCACAC', 'GCGTGCA']\"]\n",
    "    df4_rc = clustered_MREs_df_over2_below4[clustered_MREs_df_over2_below4[\"MRE_type_str\"]==\"['TGCACGC', 'GTGTGCA', 'TGCACGC']\"]\n",
    "    df4_rc_rev = clustered_MREs_df_over2_below4[clustered_MREs_df_over2_below4[\"MRE_type_str\"]==\"['TGCACGC', 'GTGTGCA', 'TGCACGC']\"]\n",
    "    #chr1_CW_variant\n",
    "    df5 = clustered_MREs_df_over2_below4[clustered_MREs_df_over2_below4[\"MRE_type_str\"]==\"['GAGTGCA', 'TGCACGC', 'GTGTGCA']\"]\n",
    "    df5_rev = clustered_MREs_df_over2_below4[clustered_MREs_df_over2_below4[\"MRE_type_str\"]==\"['GTGTGCA', 'TGCACGC', 'GAGTGCA']\"]\n",
    "    df6_rc = clustered_MREs_df_over2_below4[clustered_MREs_df_over2_below4[\"MRE_type_str\"]==\"['TGCACTC', 'GCGTGCA', 'TGCACAC']\"]\n",
    "    df6_rc_rev = clustered_MREs_df_over2_below4[clustered_MREs_df_over2_below4[\"MRE_type_str\"]==\"['TGCACAC', 'GCGTGCA', 'TGCACTC']\"]\n",
    "    \n",
    "    frames = [df1, df1_rev, df2_rc, df2_rc_rev, df3, df3_rev, df4_rc, df4_rc_rev, df5, df5_rev, df6_rc, df6_rc_rev]\n",
    "    result = pd.concat(frames)\n",
    "\n",
    "    #prepare dataframes for gene discovery\n",
    "    df1 = result.sort_values(by='MRE_zone_start').reset_index(drop=True)\n",
    "    df2 = df_chr1_gtf.sort_values(by='Start').reset_index(drop=True)\n",
    "    \n",
    "    #discover genes within 150 bp of MRE clusters\n",
    "    df2[\"Start\"] = df2[\"Start\"] - 150\n",
    "    df2[\"End\"] = df2[\"End\"] + 150\n",
    "    \n",
    "    #create a boolean mask for overlapping ranges of MRE clusters and predicted genes\n",
    "    mask = (df1['MRE_zone_end'].values[:, np.newaxis] >= df2['Start'].values) & (df1['MRE_zone_start'].values[:, np.newaxis] <= df2['End'].values)\n",
    "\n",
    "    # Check if there is any overlap (True in the mask) - discover genes\n",
    "    if mask.any():\n",
    "        overlapping_rows = np.where(mask)\n",
    "        for i, j in zip(*overlapping_rows):\n",
    "            print(f\"Overlap between row {i} in df1 and row {j} in df2.\")\n",
    "\n",
    "    #finalise data structures\n",
    "    overlap_indices = np.argwhere(mask)\n",
    "    overlapping_ranges = pd.DataFrame(overlap_indices, columns=['df1_row', 'df2_row'])\n",
    "    overlapping_ranges = overlapping_ranges.join(df1, on='df1_row')\n",
    "    overlapping_ranges = overlapping_ranges.join(df2, on='df2_row')\n",
    "    overlapping_ranges.drop_duplicates(subset=\"TrName\")\n",
    "    overlapping_ranges[\"Chromosome\"] = i\n",
    "    overlapping_ranges_all = pd.concat([overlapping_ranges_all, overlapping_ranges])\n",
    "    print(i, \"/18\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "791044e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#import UPIMAPI annotations\n",
    "upimapi_df = pd.read_csv(\"/home/Maxim/Bioinformatics/General_software/UPIMAPI/output/CambridgeChrs_adjusted/CambridgeChrs_adjusted/merged/Cambridge_UPIMAPI.tsv\", delimiter=\"\\t\")\n",
    "qseqid_df = upimapi_df[\"qseqid\"].copy()\n",
    "upimapi_df.insert(0, \"TrName\", qseqid_df.str[:-2])\n",
    "test_df_upimapi = upimapi_df.astype(str).groupby(\"TrName\")[['qseqid', 'sseqid', 'pident', 'length', 'mismatch', 'gapopen',\n",
    "       'qstart', 'qend', 'sstart', 'send', 'evalue', 'bitscore', 'Entry',\n",
    "       'Entry Name', 'Gene Names', 'Protein names', 'EC number',\n",
    "       'Function [CC]', 'Pathway', 'Keywords', 'Protein existence',\n",
    "       'Gene Ontology (GO)', 'Protein families', 'Taxonomic lineage',\n",
    "       'Taxonomic lineage (Ids)', 'Organism (ID)', 'BioCyc', 'BRENDA', 'CDD',\n",
    "       'eggNOG', 'Ensembl', 'InterPro', 'KEGG', 'Pfam', 'Reactome', 'RefSeq',\n",
    "       'UniPathway', 'Taxonomic lineage (SUPERKINGDOM)',\n",
    "       'Taxonomic lineage (PHYLUM)', 'Taxonomic lineage (CLASS)',\n",
    "       'Taxonomic lineage (ORDER)', 'Taxonomic lineage (FAMILY)',\n",
    "       'Taxonomic lineage (GENUS)', 'Taxonomic lineage (SPECIES)']].agg('! '.join).apply(lambda x: x[:].str.split('!')).reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "826e35db",
   "metadata": {},
   "outputs": [],
   "source": [
    "#couple discovered candidate genes to their annotations\n",
    "final = overlapping_ranges_all.merge(test_df_upimapi, on=[\"TrName\"], how=\"left\")\n",
    "\n",
    "#gene/MRE cluster coordinates columns\n",
    "final[\"MRE_zone_start_str\"] = final[\"MRE_zone_start\"].astype(str)\n",
    "final[\"MRE_zone_end_str\"] = final[\"MRE_zone_end\"].astype(str)\n",
    "final[\"MRE_coordinates\"] = final[\"MRE_zone_start_str\"] + \" - \" + final[\"MRE_zone_end_str\"]\n",
    "\n",
    "final[\"Start_str\"] = final[\"Start\"].astype(str)\n",
    "final[\"End_str\"] = final[\"End\"].astype(str)\n",
    "final[\"Transcript_coordinates\"] = final[\"Start_str\"] + \" - \" + final[\"End_str\"]\n",
    "\n",
    "#filter redundant columns and save the list of candidantes\n",
    "final_short = final[[\"Chr\",\"MRE_type\", \"MRE_coordinates\", \"TrName\", \"Transcript_coordinates\", \"Entry\"]]\n",
    "final_short.rename(columns={\"TrName\":\"Transcript_name\", \"Entry\":\"UniProt_Entry\", \"Chr\":\"Chromosome\"}).to_csv(\"./MRE_gene_candidates.csv\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}