Check_reported_intragenic_IESs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Previously reported intragenic IESs\n",
    "\n",
    "Five genes were previously reported as containing IESs in coding sequences (Hamilton et al., 2016, https://dx.doi.org/10.7554%2FeLife.19090). There are six IESs in total, because one gene overlaps two IESs. Compared to other Tetrahymena IESs, these intragenic IESs were reported to be shorter and precisely excised, with a conserved TTAA junction.\n",
    "\n",
    "We wish to verify that these previously reported IESs were recovered by BleTIES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pybedtools as pbt\n",
    "from Bio import SeqIO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import BletIES predictions and gene annotations\n",
    "# IES predictions from PacBio data\n",
    "bleties_pb = pbt.BedTool('tthe_pb_clr.milraa_subreads.comb.milraa_ies.gff3')\n",
    "# IES predictions from Nanopore data\n",
    "bleties_ont = pbt.BedTool('tthe_ont.milraa_subreads.comb.milraa_ies.gff3')\n",
    "# Gene annotations from ciliate.org\n",
    "genepred = pbt.BedTool('ref/2-upd-Genome-GFF3-latest-2.no_manual.gff3')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Genes containing intragenic IESs\n",
    "# List from Figure 8 of Hamilton et al. 2016, https://dx.doi.org/10.7554/eLife.19090.024\n",
    "# One of the locus IDs in the figure was missing a leading 0: TTHERM_00420400\n",
    "shortlist = ['TTHERM_00420400','TTHERM_00464970', 'TTHERM_00142380', 'TTHERM_00569290', 'TTHERM_00348490']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chr_013\tMILRAA\tinternal_eliminated_sequence_junction\t275467\t275467\t0.9357\t.\t.\tID=BREAK_POINTS_SUBREADS_chr_013_275467_136;IES_length=136;IES_subread_coverage=131;IES_zmw_coverage=131;average_subread_coverage=140;average_zmw_coverage=140;pointer_seq=TTAA;ta_pointer_seq=TAA;ta_pointer_start=275468;ta_pointer_end=275468;\n",
      "chr_015\tMILRAA\tinternal_eliminated_sequence_junction\t409851\t409851\t0.8015\t.\t.\tID=BREAK_POINTS_SUBREADS_chr_015_409851_455;IES_length=455;IES_subread_coverage=105;IES_zmw_coverage=105;average_subread_coverage=131;average_zmw_coverage=131;pointer_seq=TAA;ta_pointer_seq=TAA;ta_pointer_start=409851;ta_pointer_end=409851;\n",
      "chr_015\tMILRAA\tinternal_eliminated_sequence_junction\t412429\t412429\t0.8908\t.\t.\tID=BREAK_POINTS_SUBREADS_chr_015_412429_385;IES_length=385;IES_subread_coverage=106;IES_zmw_coverage=106;average_subread_coverage=119;average_zmw_coverage=119;pointer_seq=TTAA;ta_pointer_seq=TAA;ta_pointer_start=412430;ta_pointer_end=412430;\n",
      "chr_029\tMILRAA\tinternal_eliminated_sequence_junction\t514927\t514927\t0.9216\t.\t.\tID=BREAK_POINTS_SUBREADS_chr_029_514927_102;IES_length=102;IES_subread_coverage=141;IES_zmw_coverage=141;average_subread_coverage=153;average_zmw_coverage=153;\n",
      "chr_033\tMILRAA\tinternal_eliminated_sequence_junction\t780649\t780649\t0.9281\t.\t.\tID=BREAK_POINTS_SUBREADS_chr_033_780649_198;IES_length=198;IES_subread_coverage=129;IES_zmw_coverage=129;average_subread_coverage=139;average_zmw_coverage=139;pointer_seq=TTAA;ta_pointer_seq=TAA;ta_pointer_start=780650;ta_pointer_end=780650;\n",
      "chr_124\tMILRAA\tinternal_eliminated_sequence_junction\t1038139\t1038139\t0.8661\t.\t.\tID=BREAK_POINTS_SUBREADS_chr_124_1038139_327;IES_length=327;IES_subread_coverage=97;IES_zmw_coverage=97;average_subread_coverage=112;average_zmw_coverage=112;pointer_seq=TTAA;ta_pointer_seq=TAA;ta_pointer_start=1038140;ta_pointer_end=1038140;\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Overlaps with PacBio IES predictions\n",
    "print(\n",
    "    bleties_pb.intersect(\n",
    "        genepred.filter(lambda x: x.attrs and 'Name' in x.attrs and x.attrs['Name'] in shortlist)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chr_013\tMILRAA\tinternal_eliminated_sequence_junction\t275467\t275467\t0.8195\t.\t.\tID=BREAK_POINTS_SUBREADS_chr_013_275467_133;IES_length=133;IES_subread_coverage=109;IES_zmw_coverage=109;average_subread_coverage=133;average_zmw_coverage=133;pointer_seq=TTAA;ta_pointer_seq=TAA;ta_pointer_start=275468;ta_pointer_end=275468;\n",
      "chr_015\tMILRAA\tinternal_eliminated_sequence_junction\t409850\t409850\t0.6954\t.\t.\tID=BREAK_POINTS_SUBREADS_chr_015_409850_450;IES_length=450;IES_subread_coverage=105;IES_zmw_coverage=105;average_subread_coverage=151;average_zmw_coverage=151;pointer_seq=TTAA;ta_pointer_seq=TAA;ta_pointer_start=409851;ta_pointer_end=409851;\n",
      "chr_015\tMILRAA\tinternal_eliminated_sequence_junction\t412429\t412429\t0.7241\t.\t.\tID=BREAK_POINTS_SUBREADS_chr_015_412429_381;IES_length=381;IES_subread_coverage=105;IES_zmw_coverage=105;average_subread_coverage=145;average_zmw_coverage=145;pointer_seq=TTAA;ta_pointer_seq=TAA;ta_pointer_start=412430;ta_pointer_end=412430;\n",
      "chr_029\tMILRAA\tinternal_eliminated_sequence_junction\t514921\t514921\t0.7347\t.\t.\tID=BREAK_POINTS_SUBREADS_chr_029_514921_185;IES_length=185;IES_subread_coverage=108;IES_zmw_coverage=108;average_subread_coverage=147;average_zmw_coverage=147;pointer_seq=TTAA;ta_pointer_seq=TAA;ta_pointer_start=514922;ta_pointer_end=514922;\n",
      "chr_033\tMILRAA\tinternal_eliminated_sequence_junction\t780649\t780649\t0.8151\t.\t.\tID=BREAK_POINTS_SUBREADS_chr_033_780649_192;IES_length=192;IES_subread_coverage=97;IES_zmw_coverage=97;average_subread_coverage=119;average_zmw_coverage=119;pointer_seq=TTAA;ta_pointer_seq=TAA;ta_pointer_start=780650;ta_pointer_end=780650;\n",
      "chr_124\tMILRAA\tinternal_eliminated_sequence_junction\t1038139\t1038139\t0.6917\t.\t.\tID=BREAK_POINTS_SUBREADS_chr_124_1038139_325;IES_length=325;IES_subread_coverage=83;IES_zmw_coverage=83;average_subread_coverage=120;average_zmw_coverage=120;pointer_seq=TTAA;ta_pointer_seq=TAA;ta_pointer_start=1038140;ta_pointer_end=1038140;\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Overlaps with Nanopore IES predictions\n",
    "print(\n",
    "    bleties_ont.intersect(\n",
    "        genepred.filter(lambda x: x.attrs and 'Name' in x.attrs and x.attrs['Name'] in shortlist)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All six previously described intragenic IESs were recovered by BleTIES.\n",
    "\n",
    "The annotations from Nanopore data appear to be more accurate, because all six have the conserved TTAA tandem repeats at IES-MDS junctions that were previously reported.\n",
    "\n",
    "BleTIES probably performs better with Nanopore data because of the SPOA aligner/assembler settings."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:bio-env]",
   "language": "python",
   "name": "conda-env-bio-env-py"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}