Domain_X_Group_ChiSqr.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "similar-preliminary",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import chisquare\n",
    "import numpy as np\n",
    "import re"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "charming-spyware",
   "metadata": {},
   "outputs": [],
   "source": [
    "def metadata_dommap(f):\n",
    "    file = open(f,\"r\")\n",
    "    raw_data = file.readlines()[21:25]\n",
    "    file.close()\n",
    "    total = int(raw_data[0].rsplit()[-1])\n",
    "    tot_nc = int(raw_data[1].rsplit()[-2])\n",
    "    tot_cp = int(raw_data[2].rsplit()[-2])\n",
    "    tot_is = int(raw_data[3].rsplit()[-2])\n",
    "    tot_cnt = total - tot_nc - tot_cp - tot_is\n",
    "    \n",
    "    return {'NC': tot_nc, 'CP': tot_cp, 'IS': tot_is, 'CNT': tot_cnt, 'Total': total}\n",
    "\n",
    "def read_dommap_tsv(f):\n",
    "    file = open(f,\"r\")\n",
    "    raw_data = file.readlines()[31:]\n",
    "    file.close()\n",
    "    data = [line.rsplit(\"\\t\")[:-1] for line in raw_data]\n",
    "    return data\n",
    "\n",
    "def count_X_grp_properties(data):\n",
    "    X_grp_cnt = {}\n",
    "    for d in data:\n",
    "        uniprot = d[0].split(\"|\")[1]\n",
    "        prop = d[3]\n",
    "        X_grp = d[5]\n",
    "        F_grp = d[7]\n",
    "\n",
    "        if X_grp not in X_grp_cnt.keys():\n",
    "            X_grp_cnt[X_grp] = {\"NC\":0,\"CP\":0,\"IS\":0,\"CNT\":0,\"Total\":0,\"ExpNC\":0,\"ExpCP\":0,\"ExpIS\":0,\"ExpCNT\":0,\n",
    "                                \"NC_Fgrp\":[],\"CP_Fgrp\":[],\"IS_Fgrp\":[]}\n",
    "\n",
    "        if re.findall(\"NC\",prop):\n",
    "            X_grp_cnt[X_grp][\"NC\"] += 1\n",
    "            X_grp_cnt[X_grp][\"NC_Fgrp\"].append([uniprot,F_grp])\n",
    "\n",
    "        if re.findall(\"CP\",prop):\n",
    "            X_grp_cnt[X_grp][\"CP\"] += 1\n",
    "            X_grp_cnt[X_grp][\"CP_Fgrp\"].append([uniprot,F_grp])\n",
    "\n",
    "        if re.findall(\"IS\",prop):\n",
    "            X_grp_cnt[X_grp][\"IS\"] += 1\n",
    "            X_grp_cnt[X_grp][\"IS_Fgrp\"].append([uniprot,F_grp])\n",
    "\n",
    "        if not prop:\n",
    "            X_grp_cnt[X_grp][\"CNT\"] += 1\n",
    "        X_grp_cnt[X_grp][\"Total\"] += 1\n",
    "        \n",
    "    return X_grp_cnt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "nuclear-valley",
   "metadata": {},
   "outputs": [],
   "source": [
    "Mapped_Domains = [\n",
    "\"Yeast.mapped.tsv\",\n",
    "\"Maize.mapped.tsv\",\n",
    "\"FruitFly.mapped.tsv\",\n",
    "\"Methano.mapped.tsv\",\n",
    "\"Ecoli.mapped.tsv\",\n",
    "\"Mouse.mapped.tsv\",\n",
    "\"Human.mapped.tsv\",\n",
    "\"Staph.mapped.tsv\"\n",
    "]\n",
    "\n",
    "Organism = [name.split(\".\")[0] for name in Mapped_Domains]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "attended-garlic",
   "metadata": {},
   "outputs": [],
   "source": [
    "Organism_X_grp_Count = {}\n",
    "Organism_Tot_Dom_Count = {}\n",
    "for MDom,Org in zip(Mapped_Domains, Organism):\n",
    "    \n",
    "    if Org not in Organism_X_grp_Count.keys():\n",
    "        data = read_dommap_tsv(MDom)\n",
    "        Organism_X_grp_Count[Org] = count_X_grp_properties(data)\n",
    "         \n",
    "    if Org not in Organism_Tot_Dom_Count.keys():\n",
    "        Organism_Tot_Dom_Count[Org] = metadata_dommap(MDom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "adjacent-ground",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.seterr(invalid='ignore') # when it finds divide by 0 it will set to 0\n",
    "Chi_Sq_Data = {}\n",
    "for Org in Organism:\n",
    "    # wfile = open(\"{}_Chi_Sq_Stats.tsv\".format(Org),\"w\")\n",
    "    # wfile.write(\"X-Group\\tNC Count\\tCP Count\\tIS Count\\tNormal Count\\tRow Total\\t\\tEst. NC Count\\tEst. CP Count\\tEst. IS Count\\tEst. Normal Count\\t\\tChi Square\\tP-Value\\t\\tNC Obs/Exp\\tCP Obs/Exp\\tIS Obs/Exp\\n\")\n",
    "    Chi_Sq_Data[Org] = {}\n",
    "\n",
    "    Org_X_grp_Data = Organism_X_grp_Count[Org]\n",
    "    Org_Tot_Data = Organism_Tot_Dom_Count[Org]\n",
    "    \n",
    "    exp_NC_ratio = (Org_Tot_Data[\"NC\"]/Org_Tot_Data[\"Total\"])\n",
    "    exp_CP_ratio = (Org_Tot_Data[\"CP\"]/Org_Tot_Data[\"Total\"])\n",
    "    exp_IS_ratio = (Org_Tot_Data[\"IS\"]/Org_Tot_Data[\"Total\"])\n",
    "    exp_CNT_ratio = (Org_Tot_Data[\"CNT\"]/Org_Tot_Data[\"Total\"])\n",
    "\n",
    "    for x_grp in Org_X_grp_Data.keys():\n",
    "        obs_NC = Org_X_grp_Data[x_grp][\"NC\"]\n",
    "        obs_CP = Org_X_grp_Data[x_grp][\"CP\"]\n",
    "        obs_IS = Org_X_grp_Data[x_grp][\"IS\"]\n",
    "        obs_CNT = Org_X_grp_Data[x_grp][\"CNT\"]\n",
    "        \n",
    "        tot = Org_X_grp_Data[x_grp][\"Total\"]\n",
    "        exp_NC = tot*exp_NC_ratio\n",
    "        exp_CP = tot*exp_CP_ratio\n",
    "        exp_IS = tot*exp_IS_ratio\n",
    "        exp_CNT = tot*exp_CNT_ratio\n",
    "        \n",
    "        Org_X_grp_Data[x_grp][\"ExpNC\"] = exp_NC\n",
    "        Org_X_grp_Data[x_grp][\"ExpCP\"] = exp_CP\n",
    "        Org_X_grp_Data[x_grp][\"ExpIS\"] = exp_IS\n",
    "        Org_X_grp_Data[x_grp][\"ExpCNT\"] = exp_CNT\n",
    "\n",
    "        if Org != \"Staph\":\n",
    "            chi_data = chisquare([obs_NC,obs_CP,obs_IS,obs_CNT],[exp_NC,exp_CP,exp_IS,exp_CNT])\n",
    "        else:\n",
    "            chi_data = chisquare([obs_NC,obs_IS,obs_CNT],[exp_NC,exp_IS,exp_CNT])\n",
    "\n",
    "        Chi_Sq_Data[Org][x_grp] = chi_data\n",
    "\n",
    "        if exp_NC > 0:\n",
    "            NC_oe = obs_NC/exp_NC\n",
    "        else:\n",
    "            NC_oe = 0\n",
    "\n",
    "        if exp_CP > 0:\n",
    "            CP_oe = obs_CP/exp_CP\n",
    "        else:\n",
    "            CP_oe = 0\n",
    "\n",
    "        if exp_IS > 0:\n",
    "            IS_oe = obs_IS/exp_IS\n",
    "        else:\n",
    "            IS_oe = 0\n",
    "        \n",
    "        write_str = \"{}\\t{}\\t{}\\t{}\\t{}\\t{}\\t\\t{}\\t{}\\t{}\\t{}\\t\\t{}\\t{}\\t\\t{}\\t{}\\t{}\\n\".format(x_grp,obs_NC,obs_CP,obs_IS,obs_CNT,tot,exp_NC,exp_CP,exp_IS,exp_CNT,chi_data.statistic,chi_data.pvalue,NC_oe,CP_oe,IS_oe)\n",
    "    #     print(chisquare([obs_NC,obs_IS,obs_CNT],[exp_NC,exp_IS,exp_CNT]))\n",
    "    #     wfile.write(write_str)\n",
    "    \n",
    "    # wfile.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "internal-morgan",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "major-elder",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Neurotransmitter-gated ion-channel transmembrane pore ::\n",
      "1.0\n",
      "75.82, 0.00, 0.00, 24.18\n",
      "153 / 0.0E+00\n",
      "Proton glutamate symport protein ::\n",
      "1.0\n",
      "57.14, 0.00, 0.00, 42.86\n",
      "63 / 2.3E-87\n",
      "dicarboxylate/sodium symporter ::\n",
      "1.0\n",
      "45.95, 0.00, 0.00, 54.05\n",
      "37 / 9.6E-32\n",
      "Sodium/Calcium exchanger ::\n",
      "1.1428571428571428\n",
      "20.63, 14.29, 0.00, 79.37\n",
      "63 / 2.8E-73\n",
      "Voltage-gated ion channels ::\n",
      "1.0214909335124245\n",
      "36.27, 2.08, 5.51, 58.29\n",
      "1489 / 0.0E+00\n",
      "Sodium:neurotransmitter symporter family (SNF)-like ::\n",
      "1.0\n",
      "22.78, 0.00, 0.00, 77.22\n",
      "259 / 6.4E-44\n"
     ]
    }
   ],
   "source": [
    "# intresting_xgrp = [\"NAP-like\",\"Class II aaRS and biotin synthetases\",\"HUP domain-like\", \"Alkaline phosphatase-like\",\"DCoH-like\",\"Anticodon-binding domain of a subclass of class I aminoacyl-tRNA synthetases\", \"TIM beta/alpha-barrel\", \"triple barrel\", \"Cysteine proteinases-like\"]\n",
    "intresting_xgrp = [\"Neurotransmitter-gated ion-channel transmembrane pore\", \"Proton glutamate symport protein\", \"dicarboxylate/sodium symporter\",\"Sodium/Calcium exchanger\",\"Voltage-gated ion channels\",\"Sodium:neurotransmitter symporter family (SNF)-like\"]\n",
    "# intresting_xgrp = [\"TIM beta/alpha-barrel\"]\n",
    "\n",
    "# intresting_xgrp = [\"Metal cation-transporting ATPase, ATP-binding domain\",\n",
    "# \"DnaJ/Hsp40 cysteine-rich domain\",\n",
    "# \"FAD-linked reductases, C-terminal domain-like\",\n",
    "# \"Rubredoxin-like\",\n",
    "# \"FKBP-like\",\n",
    "# \"cradle loop barrel\",\n",
    "# \"Flavodoxin-like\"]\n",
    "\n",
    "# intresting_xgrp=[\"Glycogen debranching enzyme (GDE) insertion domain\",\n",
    "# \"NADPH-cytochrome p450 reductase helical insertion domain\",\n",
    "# \"tRNA pseudouridine synthase TruD insertion domain\",\n",
    "# \"Prokaryotic AspRS, insert domain\",\n",
    "# \"Insertion subdomain in DsbA-like\",\n",
    "# \"insertion domain in beta subunit of DNA dependent RNA-polymerase\",\n",
    "# \"Transducin (alpha subunit), insertion domain\"]\n",
    "\n",
    "# intresting_xgrp = [\"Rubredoxin-like\"]\n",
    "\n",
    "Org = \"Human\"\n",
    "# for Org in Chi_Sq_Data.keys():\n",
    "#     print(Org,\"::\")\n",
    "\n",
    "for xgrp in intresting_xgrp:\n",
    "    print(xgrp,\"::\")\n",
    "    Chi_Stat = Chi_Sq_Data[Org][xgrp]\n",
    "    chi_val, p_val = Chi_Stat\n",
    "\n",
    "    Prop_Count = Organism_X_grp_Count[Org][xgrp]\n",
    "    nc_cnt, cp_cnt, is_cnt, norm_cnt, tot, nc_est, cp_est, is_est, norm_est = [Prop_Count[p] for p in Prop_Count][:9]\n",
    "    print(sum([nc_cnt/tot, cp_cnt/tot, is_cnt/tot, norm_cnt/tot]))\n",
    "    print(f\"{100*nc_cnt/tot:3.2f}, {100*cp_cnt/tot:3.2f}, {100*is_cnt/tot:3.2f}, {100*norm_cnt/tot:3.2f}\")\n",
    "    # print(f\"{100*is_cnt/tot:3.2f}, {100*nc_cnt/tot:3.2f}, {100*cp_cnt/tot:3.2f}, {100*norm_cnt/tot:3.2f}\")\n",
    "    print(tot, \"/\",f\"{Chi_Stat.pvalue:3.1E}\")\n",
    "    # print(nc_cnt, cp_cnt, is_cnt, norm_cnt, tot, 0,nc_est, cp_est, is_est, norm_est, 0,chi_val, p_val, nc_cnt/nc_est if nc_est > 0 else 0, cp_cnt/cp_est if cp_est > 0 else 0, is_cnt/is_est if is_est > 0 else 0)\n",
    "    # print((\"{}\\t\"*11).format(*[Prop_Count[p] for p in Prop_Count],Chi_Stat.statistic,Chi_Stat.pvalue))\n",
    "    # print(100*np.array([Prop_Count[p] for p in Prop_Count][0:4])/sum([Prop_Count[p] for p in Prop_Count][0:4]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "macro-string",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "eligible-declaration",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['SRP54', 'ABC_tran_1', 'Septin', 'Beta-Casp', 'CLP1_P', 'MnmE_helical_2nd', 'TIP49_1st', 'Sigma54_activat']\n",
      "['AAA_33', 'ABC_membrane', 'Sigma54_activat', 'Lactamase_B_6', 'MnmE_helical_2nd', 'Fzo_mitofusin', 'AAA_7_N', 'AAA_19']\n",
      "[['SRP54', 'AAA_33'], ['ABC_tran_1', 'ABC_membrane'], ['Septin', 'Sigma54_activat'], ['Beta-Casp', 'Lactamase_B_6'], ['CLP1_P', 'MnmE_helical_2nd'], ['MnmE_helical_2nd', 'Fzo_mitofusin'], ['MnmE_helical_2nd', 'MnmE_helical_2nd'], ['TIP49_1st', 'AAA_7_N'], ['Sigma54_activat', 'AAA_19']] \n",
      "\n",
      "{'AAA_33': ['SRP54'], 'ABC_membrane': ['ABC_tran_1'], 'Sigma54_activat': ['Septin'], 'Lactamase_B_6': ['Beta-Casp'], 'MnmE_helical_2nd': ['CLP1_P', 'MnmE_helical_2nd'], 'Fzo_mitofusin': ['MnmE_helical_2nd'], 'AAA_7_N': ['TIP49_1st'], 'AAA_19': ['Sigma54_activat']} \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Finds the NC F_grp that the IS F_grp of a given X_grp is within\n",
    "\n",
    "hum_data = read_dommap_tsv(\"Human.mapped.tsv\")\n",
    "\n",
    "unique_IS_f_grp = [] # gives all insertional f groups without repeats\n",
    "unique_NC_f_grp = [] # gives all host f groups without repeats\n",
    "pair_f_grp = [] # gives pairs in [IS F group, NC (host) F group] for all insertional domains\n",
    "NC_group_IS_pairs = {} # gives a dict with NC f group keys and values with all IS f groups that insert into it, better organized than pairs\n",
    "for pair in Organism_X_grp_Count[\"Human\"][\"P-loop domains-like\"][\"IS_Fgrp\"]:\n",
    "    \n",
    "    for query in hum_data:\n",
    "        if pair[0] in query[0] and \"NC\" in query[3]:\n",
    "            nc_fgrp = query[7]\n",
    "            if nc_fgrp not in unique_NC_f_grp:\n",
    "                unique_NC_f_grp.append(nc_fgrp)\n",
    "            if [pair[1],nc_fgrp] not in pair_f_grp:\n",
    "                pair_f_grp.append([pair[1],nc_fgrp])\n",
    "            if nc_fgrp not in NC_group_IS_pairs.keys():\n",
    "                NC_group_IS_pairs[nc_fgrp] = []\n",
    "            if pair[1] not in NC_group_IS_pairs[nc_fgrp]:\n",
    "                NC_group_IS_pairs[nc_fgrp].append(pair[1])\n",
    "            \n",
    "\n",
    "    if pair[1] not in unique_IS_f_grp:\n",
    "        unique_IS_f_grp.append(pair[1])\n",
    "\n",
    "\n",
    "print(unique_IS_f_grp)\n",
    "\n",
    "print(unique_NC_f_grp)\n",
    "\n",
    "print(pair_f_grp,\"\\n\")\n",
    "\n",
    "print(NC_group_IS_pairs,\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "20976ab6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7310\n",
      "Yeast N = 7,310 (NC: 10.74%, CP: 0.26%, IS: 2.63%, CON: 86.37%)\n",
      "NC 785 CP 19 IS 192 CON 6,314\n",
      "NC/IS : 4.089\n",
      "119081\n",
      "Maize N = 119,081 (NC: 8.95%, CP: 0.53%, IS: 2.18%, CON: 88.34%)\n",
      "NC 10,654 CP 631 IS 2,598 CON 105,198\n",
      "NC/IS : 4.101\n",
      "44432\n",
      "FruitFly N = 44,432 (NC: 6.13%, CP: 0.33%, IS: 2.44%, CON: 91.10%)\n",
      "NC 2,722 CP 146 IS 1,086 CON 40,478\n",
      "NC/IS : 2.506\n",
      "2256\n",
      "Methano N = 2,256 (NC: 6.47%, CP: 0.44%, IS: 4.26%, CON: 88.83%)\n",
      "NC 146 CP 10 IS 96 CON 2,004\n",
      "NC/IS : 1.521\n",
      "6162\n",
      "Ecoli N = 6,162 (NC: 4.98%, CP: 0.13%, IS: 2.82%, CON: 92.06%)\n",
      "NC 307 CP 8 IS 174 CON 5,673\n",
      "NC/IS : 1.764\n",
      "104048\n",
      "Mouse N = 104,048 (NC: 4.57%, CP: 0.35%, IS: 2.74%, CON: 92.35%)\n",
      "NC 4,752 CP 359 IS 2,848 CON 96,089\n",
      "NC/IS : 1.669\n",
      "155158\n",
      "Human N = 155,158 (NC: 4.52%, CP: 0.40%, IS: 1.72%, CON: 93.36%)\n",
      "NC 7,013 CP 626 IS 2,661 CON 144,858\n",
      "NC/IS : 2.635\n",
      "3582\n",
      "Staph N = 3,582 (NC: 4.05%, CP: 0.14%, IS: 3.85%, CON: 91.96%)\n",
      "NC 145 CP 5 IS 138 CON 3,294\n",
      "NC/IS : 1.051\n"
     ]
    }
   ],
   "source": [
    "for org, data in Organism_Tot_Dom_Count.items():\n",
    "    nc = data[\"NC\"]\n",
    "    cp = data[\"CP\"]\n",
    "    iss = data[\"IS\"]\n",
    "    tot = data[\"Total\"]\n",
    "    cnt = tot - (nc + cp + iss)\n",
    "    pNC = data[\"NC\"]/data[\"Total\"]\n",
    "    pCP = data[\"CP\"]/data[\"Total\"]\n",
    "    pIS = data[\"IS\"]/data[\"Total\"]\n",
    "    pCON = data[\"CNT\"]/data[\"Total\"]\n",
    "    \n",
    "    print(data[\"Total\"])\n",
    "    print(org, f\"N = {tot:,} (NC: {pNC:2.2%}, CP: {pCP:2.2%}, IS: {pIS:2.2%}, CON: {pCON:2.2%})\")\n",
    "    \n",
    "    print(f\"NC {nc:,} CP {cp:,} IS {iss:,} CON {cnt:,}\")\n",
    "\n",
    "    print(f\"NC/IS : {nc/iss:3.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01e20e6e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "7e85ec9bf098c5427e45e2f632dcd4eeff803b007e1abd287d600879388709c1"
  },
  "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}