CarbonSources/ComputeModelDistances.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6adb9ac6-a233-48d2-9ce7-92d9fec9d097",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cobrakbase 0.3.1\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "sys.path.insert(0,\"/Users/chenry/code/ModelSEEDpy\")\n",
    "sys.path.append(\"/Users/chenry/code/cobrakbase\")\n",
    "#sys.path.append(\"/Users/piehl/Desktop/S2022/ModelSEEDpy-main/modelseedpy\")\n",
    "#sys.path.insert(0,\"/Users/piehl/Desktop/S2022/ModelSEEDpy\")\n",
    "#sys.path.append(\"/Users/piehl/Desktop/S2022/cobrakbase-cobra-model\")\n",
    "import os\n",
    "#os.environ[\"HOME\"] = \"/Users/piehl/Desktop/S2022\"\n",
    "import json\n",
    "import cobra\n",
    "import cobrakbase\n",
    "from modelseedpy import MSPackageManager, MSGrowthPhenotypes\n",
    "from cobrakbase.core.kbasefba.fbamodel_from_cobra import CobraModelConverter\n",
    "import pandas as pd\n",
    "from IPython.display import HTML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2e3c88e4-2507-43b5-9d6b-551481b58ff1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "562.8594;not found\n",
      "562.85941;not found\n",
      "562.81369;not found\n",
      "562.81371;not found\n",
      "562.81372;not found\n",
      "562.81363;not found\n",
      "562.85952;not found\n",
      "562.85953;not found\n",
      "562.85942;not found\n",
      "562.85943;not found\n",
      "562.85944;not found\n",
      "562.85945;not found\n",
      "562.85946;not found\n",
      "562.85947;not found\n",
      "562.85948;not found\n",
      "562.85949;not found\n",
      "562.8595;not found\n",
      "562.85951;not found\n",
      "562.55367;no metadata\n",
      "562.55368;no metadata\n",
      "562.55380;no metadata\n",
      "562.55381;no metadata\n",
      "562.55382;no metadata\n",
      "562.55383;no metadata\n",
      "562.55384;no metadata\n",
      "562.55385;no metadata\n",
      "562.55386;no metadata\n",
      "562.55390;no metadata\n",
      "562.55391;no metadata\n",
      "562.55413;no metadata\n",
      "562.55415;no metadata\n",
      "562.55417;no metadata\n",
      "562.55418;no metadata\n",
      "562.55483;no metadata\n",
      "562.55540;no metadata\n",
      "562.55545;no metadata\n",
      "562.55546;no metadata\n",
      "562.55547;no metadata\n",
      "562.55548;no metadata\n",
      "562.55550;no metadata\n",
      "562.55551;no metadata\n",
      "562.55552;no metadata\n",
      "562.55554;no metadata\n",
      "562.55556;no metadata\n",
      "562.55557;no metadata\n",
      "562.55559;no metadata\n",
      "562.55561;no metadata\n",
      "562.55562;no metadata\n",
      "562.55563;no metadata\n",
      "562.55564;no metadata\n",
      "562.55565;no metadata\n",
      "562.55566;no metadata\n",
      "562.55567;no metadata\n",
      "562.55568;no metadata\n",
      "562.55569;no metadata\n",
      "562.55570;no metadata\n",
      "562.55571;no metadata\n",
      "562.55572;no metadata\n",
      "562.55573;no metadata\n",
      "562.55883;no metadata\n",
      "AB_M301_2012_EC11;no translation\n",
      "AB_FA15_FC_1.3;no translation\n",
      "AB_FA15_GC_1.2;no translation\n",
      "AB_FA15_GC_1.1;no translation\n",
      "AB_FA15_FC_1.1;no translation\n",
      "AB_FA15_DD_1.2;no translation\n",
      "AB_FA15_DD_1.3;no translation\n",
      "AB_FA15_DD_1.1;no translation\n",
      "AB_FA15_HF_1.3;no translation\n",
      "AB_FA15_PTC_1.3;no translation\n",
      "AB_SUM15_AL_4.1;no translation\n",
      "AB_FA15_HL_1.2;no translation\n",
      "AB_BQ_SUM18_MCR_2C_9;no translation\n",
      "AB_SUM_2014_141st_7.2;no translation\n",
      "AB_SUM_2014_DD_6.1;no translation\n",
      "AB_SUM15_HF_4.1;no translation\n",
      "AB_SUM16_AL_10.2;no translation\n",
      "AB_SUM16_NB_13.1;no translation\n",
      "AB_SUM16_UM_13.2;no translation\n",
      "AB_SUM16_DD_13.2;no translation\n",
      "AB_SUM15_SB_4.1;no translation\n",
      "AB_FA15_PPC_1.2;no translation\n",
      "AB_SUM16_NB_13.2;no translation\n",
      "AB_SUM15_FC_4.1;no translation\n",
      "AB_SUM15_GC_4.1;no translation\n",
      "AB_SUM16_UM_13.1;no translation\n",
      "AB_SUM16_PTC_10.2;no translation\n",
      "AB_SUM16_DD_13.1;no translation\n",
      "AB_SUM_2014_DD_7.3;no translation\n",
      "AB_FA15_AL_1.1;no translation\n",
      "AB_SUM15_HL_4.1;no translation\n",
      "AB_SUM16_DD_10.5;no translation\n",
      "AB_FA15_AL_1.2;no translation\n",
      "AB_SUM16_AL_13.2;no translation\n",
      "AB_SUM15_PTC_4.1;no translation\n",
      "AB_SUM15_DD_3.1;no translation\n",
      "AB_SUM16_PTC_13.1;no translation\n",
      "AB_SUM16_DD_10.1;no translation\n",
      "AB_SUM16_AL_13.1;no translation\n",
      "AB_SUM16_PTC_13.2;no translation\n"
     ]
    }
   ],
   "source": [
    "#Read in biolog data\n",
    "zpheno_list = []\n",
    "bpheno_list = []\n",
    "data = pd.read_csv(\"BiologMacTestOutput_On-Off.csv\")\n",
    "bdata = {}\n",
    "for index, row in data.iterrows():\n",
    "    media = row[\"Media\"]\n",
    "    if media[0:7] == \"Carbon-\":#Comment this out if you want to use all media\n",
    "        bpheno_list.append(media)\n",
    "        for (columnName, columnData) in data.iteritems():\n",
    "            if columnName != \"Media\":\n",
    "                genome_id = columnName[0:-9]\n",
    "                if genome_id not in bdata:\n",
    "                    bdata[genome_id] = {}    \n",
    "                bdata[genome_id][media] = row[columnName]\n",
    "#Read in metadata sheet and use to translate genome IDs\n",
    "#Deal with trailing zeros and load and save all metadata\n",
    "metadata = {}\n",
    "translation = {}\n",
    "data = pd.read_csv(\"Watershed Metadata 08_23_16-present - PatricMasterSheet 1.csv\")\n",
    "for index, row in data.iterrows():\n",
    "    genome_id = str(row[\"Patric ID\"])\n",
    "    for i in range(3):\n",
    "        if genome_id in bdata:\n",
    "            break\n",
    "        genome_id += \"0\"\n",
    "    if genome_id not in bdata:\n",
    "        print(str(row[\"Patric ID\"])+\";not found\")\n",
    "        next\n",
    "    local_genome = row[\"Local Genome ID from Patric\"]\n",
    "    local_genome = local_genome[17:]\n",
    "    translation[local_genome] = genome_id\n",
    "    metadata[genome_id] = row\n",
    "#Now making sure all genomes from clayton are in metadata\n",
    "for genome_id in bdata:\n",
    "    if genome_id not in metadata:\n",
    "        print(genome_id+\";no metadata\")\n",
    "#Read in biolog predictions from classifiers\n",
    "data = pd.read_csv(\"Aaron_predict_178.csv\")\n",
    "zdata = {}\n",
    "for index, row in data.iterrows():\n",
    "    genome = row[\"genome\"]\n",
    "    if genome[-4:] == \".txt\":\n",
    "        genome = genome[:-4]\n",
    "    if genome not in translation:\n",
    "        print(genome+\";no translation\")\n",
    "    else:\n",
    "        genome = translation[genome]\n",
    "        if genome not in zdata:\n",
    "            zdata[genome] = {}\n",
    "        for (columnName, columnData) in data.iteritems():\n",
    "            if columnName != \"genome\":\n",
    "                if columnName not in zpheno_list:\n",
    "                    zpheno_list.append(columnName)\n",
    "                zdata[genome][columnName] = row[columnName]  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "12bfdfc7-f088-4707-9594-fd540c22eea6",
   "metadata": {},
   "outputs": [],
   "source": [
    "mdlws = 119455\n",
    "#kbase_api = cobrakbase.KBaseAPI()\n",
    "kbase_api = cobrakbase.KBaseCache(token=token,dev=True)\n",
    "models = kbase_api.list_objects(mdlws, object_type=\"KBaseFBA.FBAModel\",include_metadata=True)\n",
    "model_indecies = {}\n",
    "model_names = []\n",
    "modelrxnlists = {}\n",
    "reaction_hash = {}\n",
    "for model in models:\n",
    "    model_indecies[model[1]] = len(model_names)\n",
    "    model_names.append(model[1])\n",
    "    modelobj = kbase_api.get_object(model[1],mdlws)\n",
    "    reactionlist = modelobj['modelreactions']\n",
    "    modelrxnlists[model[1]] = []\n",
    "    for reaction in reactionlist:\n",
    "        if reaction['id'] not in reaction_hash:\n",
    "            reaction_hash[reaction['id']] = reaction\n",
    "        if 'modelReactionProteins' in reaction.keys() and len(reaction['modelReactionProteins']):\n",
    "            modelrxnlists[model[1]].append(reaction['id'])\n",
    "reaction_distances = [[0 for i in range(len(models))] for j in range(len(models))]\n",
    "biolog_distances = [[0 for i in range(len(models))] for j in range(len(models))]\n",
    "classifier_distances = [[0 for i in range(len(models))] for j in range(len(models))]\n",
    "consolidated_distances = [[0 for i in range(len(models))] for j in range(len(models))]\n",
    "count = 0\n",
    "conflicts = {}\n",
    "for modelone in models:\n",
    "    conflicts[modelone[1][0:-9]] = 0\n",
    "    counttwo = 0\n",
    "    for zpheno in zpheno_list:\n",
    "        if modelone[1][0:-9] in zdata and zdata[modelone[1][0:-9]][zpheno] > 0.001:\n",
    "            if modelone[1][0:-9] in bdata and bdata[modelone[1][0:-9]][zpheno] < 0.001:\n",
    "                conflicts[modelone[1][0:-9]] += 1\n",
    "        elif modelone[1][0:-9] in bdata and bdata[modelone[1][0:-9]][zpheno] > 0.001:\n",
    "            conflicts[modelone[1][0:-9]] += 1\n",
    "    for modeltwo in models:\n",
    "        for reaction in reaction_hash.keys():\n",
    "            if reaction in modelrxnlists[modelone[1]]:\n",
    "                if reaction not in modelrxnlists[modeltwo[1]]:\n",
    "                    reaction_distances[count][counttwo] += 1\n",
    "                    reaction_distances[counttwo][count] += 1\n",
    "            if reaction in modelrxnlists[modeltwo[1]]:\n",
    "                if reaction not in modelrxnlists[modelone[1]]:\n",
    "                    reaction_distances[count][counttwo] += 1\n",
    "                    reaction_distances[counttwo][count] += 1\n",
    "        for phenotype in bpheno_list:\n",
    "            if modelone[1][0:-9] in bdata and bdata[modelone[1][0:-9]][phenotype] > 0.001:\n",
    "                if modeltwo[1][0:-9] in bdata and bdata[modeltwo[1][0:-9]][phenotype] < 0.001:\n",
    "                    biolog_distances[count][counttwo] += 1\n",
    "                    biolog_distances[counttwo][count] += 1\n",
    "            if modeltwo[1][0:-9] in bdata and bdata[modeltwo[1][0:-9]][phenotype] > 0.001:\n",
    "                if modelone[1][0:-9] in bdata and bdata[modelone[1][0:-9]][phenotype] < 0.001:\n",
    "                    biolog_distances[count][counttwo] += 1\n",
    "                    biolog_distances[counttwo][count] += 1\n",
    "        for zpheno in zpheno_list:    \n",
    "            if modelone[1][0:-9] in zdata and zdata[modelone[1][0:-9]][zpheno] > 0.001:\n",
    "                if modeltwo[1][0:-9] in zdata and zdata[modeltwo[1][0:-9]][zpheno] < 0.001:\n",
    "                    classifier_distances[count][counttwo] += 1\n",
    "                    classifier_distances[counttwo][count] += 1\n",
    "            if modeltwo[1][0:-9] in zdata and zdata[modeltwo[1][0:-9]][zpheno] > 0.001:\n",
    "                if modelone[1][0:-9] in zdata and zdata[modelone[1][0:-9]][zpheno] < 0.001:\n",
    "                    classifier_distances[count][counttwo] += 1\n",
    "                    classifier_distances[counttwo][count] += 1\n",
    "        counttwo += 1\n",
    "    count += 1\n",
    "with open('ReactionDist.json', 'w') as outfile:\n",
    "    json.dump(reaction_distances, outfile)\n",
    "with open('BiologDist.json', 'w') as outfile:\n",
    "    json.dump(biolog_distances, outfile)\n",
    "with open('ClassifierDist.json', 'w') as outfile:\n",
    "    json.dump(classifier_distances, outfile)   \n",
    "largest_rxn = None\n",
    "largest_class = None\n",
    "largest_biolog = None\n",
    "count = 0\n",
    "for modelone in models:\n",
    "    counttwo = 0\n",
    "    for modeltwo in models:\n",
    "        if not largest_rxn or reaction_distances[count][counttwo] > largest_rxn:\n",
    "            largest_rxn = reaction_distances[count][counttwo]\n",
    "        if not largest_class or classifier_distances[count][counttwo] > largest_class:\n",
    "            largest_class = classifier_distances[count][counttwo]\n",
    "        if not largest_biolog or biolog_distances[count][counttwo] > largest_biolog:\n",
    "            largest_biolog = biolog_distances[count][counttwo]\n",
    "        counttwo += 1\n",
    "    count += 1\n",
    "count = 0\n",
    "for modelone in models:\n",
    "    counttwo = 0\n",
    "    for modeltwo in models:\n",
    "        consolidated_distances[count][counttwo] = reaction_distances[count][counttwo]/largest_rxn + classifier_distances[count][counttwo]/largest_class + biolog_distances[count][counttwo]/largest_biolog\n",
    "#Cluster on consolidated and all individuals\n",
    "#Print cluster size profile for all\n",
    "#Map clusters for each method and compute similarity and show average similarity + distribution\n",
    "#Print cluster representatives for each method\n",
    "#Print metadata for clusters\n",
    "#Compute average and range of consistency values for each cluster\n",
    "#Pick representative genome with the highest inconsistency between classifiers and models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "572e2c09-26b8-47f4-bfb9-36c9675aa5b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('ConsolidatedDist.json', 'w') as outfile:\n",
    "    json.dump(consolidated_distances, outfile)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "943c0bd3-bef9-4c35-93d7-c2d567f9e13d",
   "metadata": {},
   "outputs": [],
   "source": [
    "mdlws = 119455\n",
    "#kbase_api = cobrakbase.KBaseAPI()\n",
    "kbase_api = cobrakbase.KBaseCache()\n",
    "models = kbase_api.list_objects(mdlws, object_type=\"KBaseFBA.FBAModel\",include_metadata=True)\n",
    "model_indecies = {}\n",
    "model_names = []\n",
    "modelrxnlists = {}\n",
    "reaction_hash = {}\n",
    "for model in models:\n",
    "    model_names.append(model[1])\n",
    "with open('ModelIndecies.json', 'w') as outfile:\n",
    "    json.dump(model_names, outfile)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d70cb156-c8fa-498d-a4a2-2b02a4f368c5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}