reannotation_ppacificus_all.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os.path\n",
    "from collections import defaultdict\n",
    "import statistics\n",
    "\n",
    "import matplotlib\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from orthologue_analysis.orthogroups import init_orthogroup_df\n",
    "from orthologue_analysis.species import AltSourceMixin, Species, SpeciesList\n",
    "from orthologue_analysis.utils import SequenceIDMapping, orthofinder_paths\n",
    "from reannotation.analysis import (\n",
    "    interpro_accessions_frequently_missed_by_all_tools,\n",
    "    interpro_accessions_in_novel_transcripts,\n",
    "    interpro_accessions_in_missed_transcripts,\n",
    "    missed_transcripts_with_significantly_more_frequent_accessions\n",
    ")\n",
    "from reannotation.pipelines import interpro_accession_pipeline, suspicious_orthologue_pipeline, novel_orthologue_pipeline\n",
    "from reannotation.statistics import fisher_exact_for_two_lists_of_accessions\n",
    "from reannotation.utils import extract_accessions_from_transcript\n",
    "from utils.esm import extract_esm_means\n",
    "from utils.gffutils import init_db\n",
    "\n",
    "matplotlib.use(\"TkAgg\")\n",
    "\n",
    "\n",
    "class Pristionchus(Species):\n",
    "    abbr = \"P\"\n",
    "    genus = \"pristionchus\"\n",
    "    clade = 0\n",
    "\n",
    "\n",
    "class PristionchusFromTool(AltSourceMixin, Pristionchus):\n",
    "    pass\n",
    "\n",
    "\n",
    "results_label = \"Results_Aug21\"\n",
    "wbps_ann_path = \"data/from_WBPS/pristionchus_pacificus.PRJNA12644.WBPS19.annotations.gff3\"\n",
    "braker_path = \"data/from_MARS/Pristionchus_pacificus_braker3_full.gff3\"\n",
    "helixer_path = \"data/from_MARS/Pristionchus_pacificus_helixer_full.gff3\"\n",
    "anno_path = \"data/from_EBI/pristionchus_pacificus_gca000180635v4.gff3\"\n",
    "db = init_db(wbps_ann_path, \"db/Ppac_wbps.db\")\n",
    "of = orthofinder_paths(results_label)\n",
    "\n",
    "wbps_col = \"Ppac_LT\"\n",
    "braker_col = \"Ppac_braker3_LT\"\n",
    "helixer_col = \"Ppac_helixer_LT\"\n",
    "anno_col = \"Ppac_anno_LT\"\n",
    "\n",
    "hog_df = init_orthogroup_df(of[\"orthogroups\"])\n",
    "seq_id_map = SequenceIDMapping(of[\"wd\"])\n",
    "mars_data_dir = os.path.join(\"data\", \"from_MARS\", \"\")\n",
    "ebi_data_dir = os.path.join(\"data\", \"from_EBI\", \"\")\n",
    "\n",
    "species_list = SpeciesList([\n",
    "    PristionchusFromTool(\"pacificus\", data_dir=mars_data_dir, data_label=\"Ppac_LT\", prot_filename_suffix=\".fa\"),\n",
    "    PristionchusFromTool(\"pacificus_braker3_reann\", data_dir=mars_data_dir, data_label=\"Ppac_braker3_LT\", prot_filename_suffix=\".fa\"),\n",
    "    PristionchusFromTool(\"pacificus_helixer_reann\", data_dir=mars_data_dir, data_label=\"Ppac_helixer_LT\", prot_filename_suffix=\".fa\"),\n",
    "    PristionchusFromTool(\"pacificus_anno_reann\", data_dir=ebi_data_dir, data_label=\"Ppac_anno_LT\", prot_filename_suffix=\".fa\")],\n",
    "    wd_path=of[\"wd\"],\n",
    "    load_blast=True\n",
    ")\n",
    "\n",
    "wbps_species = species_list.get_species_with_data_label(wbps_col)\n",
    "\n",
    "min_freq = 10\n",
    "\n",
    "interproscan_dir = \"data/from_MARS/interproscan/ppac\"\n",
    "\n",
    "acc_product = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# General Statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Shared orthologues with WBPS:\")\n",
    "print(\"WBPS: {}\".format(len(hog_df[~hog_df[wbps_col].isna()])))\n",
    "print(\"BRAKER3: {}\".format(len(hog_df[~hog_df[wbps_col].isna() & ~hog_df[braker_col].isna()])))\n",
    "print(\"Helixer: {}\".format(len(hog_df[~hog_df[wbps_col].isna() & ~hog_df[helixer_col].isna()])))\n",
    "print(\"Anno: {}\".format(len(hog_df[~hog_df[wbps_col].isna() & ~hog_df[anno_col].isna()])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assessing merged/split genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "def pickle_cache_suspicious_orthologue_pipeline(tool, *args, **kwargs):\n",
    "    merged_path = os.path.join(\"data\", \"tmp\", \"ppac_{}_merged.pickle\".format(tool))\n",
    "    split_path = os.path.join(\"data\", \"tmp\", \"ppac_{}_split.pickle\".format(tool))\n",
    "    if os.path.isfile(merged_path) and os.path.isfile(split_path):\n",
    "        with open(merged_path, \"rb\") as f:\n",
    "            merged = pickle.load(f)\n",
    "        with open(split_path, \"rb\") as f:\n",
    "            split = pickle.load(f)\n",
    "    else:\n",
    "        merged, split = suspicious_orthologue_pipeline(*args, **kwargs)\n",
    "        with open(merged_path, 'wb') as f:\n",
    "            pickle.dump(merged, f, protocol=pickle.HIGHEST_PROTOCOL)\n",
    "        with open(split_path, 'wb') as f:\n",
    "            pickle.dump(split, f, protocol=pickle.HIGHEST_PROTOCOL)\n",
    "    return merged, split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "braker_merged, braker_split = pickle_cache_suspicious_orthologue_pipeline(\"braker\", hog_df, wbps_col, braker_col, species_list, seq_id_map, wbps_prefix=\"Transcript\")\n",
    "anno_merged, anno_split = pickle_cache_suspicious_orthologue_pipeline(\"anno\", hog_df, wbps_col, anno_col, species_list, seq_id_map, wbps_prefix=\"Transcript\")\n",
    "helixer_merged, helixer_split = pickle_cache_suspicious_orthologue_pipeline(\"helixer\", hog_df, wbps_col, helixer_col, species_list, seq_id_map, wbps_prefix=\"Transcript\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_genes = len(list(species_list.get_species_with_data_label(\"Ppac_braker3_LT\").db.all_features(featuretype=\"gene\")))\n",
    "print(f\"BRAKER3: merged={len(braker_merged)}, split={len(braker_split)}, total={round(100*(len(braker_split) + len(braker_merged)*2)/num_genes, 2)}\")\n",
    "num_genes = len(list(species_list.get_species_with_data_label(\"Ppac_helixer_LT\").db.all_features(featuretype=\"gene\")))\n",
    "print(f\"Helixer: merged={len(helixer_merged)}, split={len(helixer_split)}, total={round(100*(len(helixer_split) + len(helixer_merged)*2)/num_genes, 2)}\")\n",
    "num_genes = len(list(species_list.get_species_with_data_label(\"Ppac_anno_LT\").db.all_features(featuretype=\"gene\")))\n",
    "print(f\"Anno: merged={len(anno_merged)}, split={len(anno_split)}, total={round(100*(len(anno_split) + len(anno_merged)*2)/num_genes, 2)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anno_merged, anno_split = suspicious_orthologue_pipeline(hog_df, wbps_col, anno_col, species_list, seq_id_map, wbps_prefix=\"Transcript\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "helixer_merged, helixer_split = suspicious_orthologue_pipeline(hog_df, wbps_col, helixer_col, species_list, seq_id_map, wbps_prefix=\"Transcript\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# InterPro accession investigation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### BRAKER3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acc_product, acc_tally_shared, acc_tally_missed_braker3, acc_tally_novel_braker3, missed_transcripts = interpro_accession_pipeline(db, hog_df, wbps_col, braker_col, interproscan_dir, acc_product, prefix=\"Transcript\")\n",
    "\n",
    "# Find InterPro accessions occurring with significantly different frequency than in control (acc_tally_shared)\n",
    "braker3_novel_results = fisher_exact_for_two_lists_of_accessions(acc_tally_novel_braker3, acc_tally_shared)\n",
    "braker3_missed_results = fisher_exact_for_two_lists_of_accessions(acc_tally_missed_braker3, acc_tally_shared)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interpro_accessions_in_missed_transcripts(acc_product, acc_tally_missed_braker3, acc_tally_novel_braker3, braker3_missed_results, braker3_novel_results, min_freq=75)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### There are no IPR accessions in the annotation... Will have to write brand new pipelines for dealing with PFAM domains (which are present)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Novel orthologues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Novel orthogroups\")\n",
    "total_orthogroups = hog_df.shape[0]\n",
    "braker3_novel_ogs = hog_df[(hog_df[wbps_col].isna()) & ~(hog_df[braker_col].isna())].shape[0]\n",
    "helixer_novel_ogs = hog_df[(hog_df[wbps_col].isna()) & ~(hog_df[helixer_col].isna())].shape[0]\n",
    "anno_novel_ogs = hog_df[(hog_df[wbps_col].isna()) & ~(hog_df[anno_col].isna())].shape[0]\n",
    "print(f\"BRAKER3: {braker3_novel_ogs} ({round(100*braker3_novel_ogs/total_orthogroups, 2)}%)\")\n",
    "print(f\"Helixer: {helixer_novel_ogs} ({round(100*helixer_novel_ogs/total_orthogroups, 2)}%)\")\n",
    "print(f\"Anno: {anno_novel_ogs} ({round(100*anno_novel_ogs/total_orthogroups, 2)}%)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# novel_orthologue_pipeline(hog_df, wbps_col, anno_col, species_list, out_dir=\"data/novel_orthologue_sequences/ppac/anno/\")\n",
    "# novel_orthologue_pipeline(hog_df, wbps_col, braker_col, species_list, out_dir=\"data/novel_orthologue_sequences/ppac/braker3/\")\n",
    "# novel_orthologue_pipeline(hog_df, wbps_col, helixer_col, species_list, out_dir=\"data/novel_orthologue_sequences/ppac/helixer/\")\n",
    "anno_esm_means = extract_esm_means(\"data/from_MARS/Ppac_esm_pLDDTs_anno.txt\").values()\n",
    "braker3_esm_means = extract_esm_means(\"data/from_MARS/Ppac_esm_pLDDTs_braker3.txt\").values()\n",
    "helixer_esm_means = extract_esm_means(\"data/from_MARS/Ppac_esm_pLDDTs_helixer.txt\").values()\n",
    "\n",
    "print(statistics.mean(map(float, anno_esm_means)))\n",
    "print(statistics.mean(map(float, braker3_esm_means)))\n",
    "print(statistics.mean(map(float, helixer_esm_means)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "cols = (\n",
    "    \"fn\",\n",
    "    \"mean\",\n",
    "    \"median\",\n",
    "    \"stdev\",\n",
    "    \"var\",\n",
    "    \"max\",\n",
    "    \"min\",\n",
    "    \"perc_confident\"\n",
    ")\n",
    "df = pd.read_csv(\"data/from_MARS/pLDDT_ppac.csv\", names=cols)\n",
    "print(f\"Mean of means: {df['mean'].mean()}\")\n",
    "print(f\"% that are \\\"Confident\\\": {100*df[df['mean'] >= 70].shape[0]/df.shape[0]}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_anno = pd.read_csv(\"data/from_MARS/pLDDT_ppac_anno.csv\", names=cols)\n",
    "df_braker3 = pd.read_csv(\"data/from_MARS/pLDDT_ppac_braker3.csv\", names=cols)\n",
    "df_helixer = pd.read_csv(\"data/from_MARS/pLDDT_ppac_helixer.csv\", names=cols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Mean of WBPS % that are \\\"Confident\\\" residues: {df['perc_confident'].mean()}\")\n",
    "print(f\"Mean of BRAKER3 % that are \\\"Confident\\\" residues: {df_braker3['perc_confident'].mean()}\")\n",
    "print(f\"Mean of Anno % that are \\\"Confident\\\" residues: {df_anno['perc_confident'].mean()}\")\n",
    "print(f\"Mean of Helixer % that are \\\"Confident\\\" residues: {df_helixer['perc_confident'].mean()}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}