notebooks/README.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0580641a-d103-4972-b263-9fcc21959431",
   "metadata": {},
   "source": [
    "# Rnalib Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d170c96b-1cc8-4d7b-bb61-fbc3f0ef0e0e",
   "metadata": {
    "tags": []
   },
   "source": [
    "The increasing abundance and complexity of biological data requires programatic approaches for their analysis that allow for the rapid development of custom yet efficient, well-tested and reproducible analysis pipelines. Traditional approaches (e.g., bash command/perl pipelines), that suffer from low readability, maintainability and reproducibility, are increasingly superseded by complex, containerized (e.g., docker/apptainer) analysis pipelines that combine specialized bioinformatics tools with custom scripts developed in higher programming languages such as *Python*, *R* or *Rust*.\n",
    "\n",
    "Despite not being the fastest option, *Python* is a [popular choice by analysts](https://github.blog/2023-03-02-why-python-keeps-growing-explained/), mainly due to the large number of available (bioinformatics) libraries, its simple and well-documented syntax and its broad utility in the emerging data science field. For this reason, an increasing number of python APIs for the handling of genomics data are being developed. \n",
    "While some are mainly wrappers around well-established and efficient bioinformatics tools (e.g., pysam/samtools, pybedtools/BEDTools or cyvcf2/htslib), others are built on popular data wrangling/analysis libraries such as *pandas* or *numpy* (e.g., bioframe/pandas, pyranges/pandas, bionumpy/numpy).\n",
    "\n",
    "Here, we introduce **rnalib**, a Python library for handling transcriptomics data.\n",
    "*Rnalib* implements a **transcriptome** class that represents genomic features (such as genes, transcripts or exons) as Python dataclasses and explicitly models their relationships. It provides a secure way for referencing and annotating such features by structuring them into immutable genomic locations (that can, e.g., safely be used in lookup tables) and mutable, arbitrary annotations. \n",
    "Users can incrementally add annotations by using a comprehensive set of **genomic iterators** that provide access to the most common genomics data formats (e.g., BAM, BED, VCF, etc.). \n",
    "*Genomic iterators* provide a uniform interface to such data sources and are implemented on top of existing, efficient python implementations (pysam, bioframe, pybedtools). \n",
    "\n",
    "*Rnalib* provides **efficient access/querying** of feature sequences and annotations and implements various ways to **export** those, e.g., to pandas dataframes, in order to provide seamless integration with other (genomics) libraries and tools.\n",
    "*Rnalib* is not meant to replace the existing (great) work of others but to complement and integrate it and to fill some missing gaps. \n",
    "\n",
    "*Rnalib*'s target audience are bioinformatics analysts and developers and its main design goal is to enable **fast, readable, reproducible and robust development of novel bioinformatics tools and methods**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63b5cb56-25e2-4b50-b189-fce7cd6d4413",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## Table of Contents\n",
    "To introduce *rnalib*, we have compiled a number of ipython notebooks:\n",
    "\n",
    "This notebook presents the core API of *rnalib*:\n",
    "* Its [genomic intervals](#Genomic-intervals) implementation\n",
    "* Its [transcriptome](#Transcriptome-implementation) implementation\n",
    "* Its [genomic iterators](#Iterators) for efficient iteration over (subregions of) large-scaled genomics datasets based on the pysam library. Iterators keep track of the genomic region of the yielded data, thereby enabling their efficient integration with other data sets\n",
    "* Its strategies for genomic [feature annotation](#Feature-Annotation)\n",
    "\n",
    "A [second (advanced usage) notebook](AdvancedUsage.ipynb) presents:\n",
    "* More advanced usage examples for the topics above\n",
    "* Various *rnalib* utility functions for working with genomics data\n",
    "\n",
    "In [a 'related work' notebook](RelatedWork_performance.ipynb), we compare rnalib to a number of related APIs (pybedtools, bioframe, HTSeq) with regard to \n",
    "* provided features/approaches\n",
    "* performance\n",
    "* common analysis pitfalls that we tried to address in the design of *rnalib*\n",
    "\n",
    "We also provide a set of tutorials for demonstrating rnalib in realistic usage scenarios:\n",
    "* [Tutorial: Transcriptome annotation with genemunge, archs4 and mygene.info](Tutorial_transcriptome_annotation.ipynb): annotation with data from public databases\n",
    "* [Tutorial: SLAM-seq analysis](Tutorial_SLAM-seq.ipynb): Simplified analysis of a SLAM-seq timecourse dataset\n",
    "* [Tutorial: Comparison of gene annotation sets](Tutorial_compare_annotation_sets.ipynb): Comparison of different gene annotation sets (human and fly)\n",
    "* [Tutorial: CTCF analysis with rnalib and bioframe](Tutorial_CTCF_analysis.ipynb): Annotation of genes with CTCF sites\n",
    "* [Tutorial: shRNA analysis](Tutorial_shRNA_analysis.ipynb): a small analysis of shRNA targets\n",
    "\n",
    "For further API examples, please also consider *rnalib*'s [test classes](https://github.com/popitsch/rnalib/tree/main/tests)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf2d6c37-851f-4421-9d72-95c004d9e8f3",
   "metadata": {},
   "source": [
    "## Requirements\n",
    "Before executing this notebook, you need to install all required *rnalib* requirements as well as optional libraries needed by this notebook.\n",
    "It is recommended to do this in a [Python virtual environment](https://rnalib.readthedocs.io/en/latest/readme.html#installation).\n",
    "\n",
    "Note that this notebook as well as *rnalib*'s testing suite use various test resources (genomics data files and indexing structures) that can be created by running the *rnalib* `rnalib_create_testdata` script or by calling the `testdata.create_testdata()` method. Refer to `testdata.py` for details how the test data files were created.\n",
    "There are two separate test resource sets:\n",
    "\n",
    "* *test_resources*: small test datasets that are used by *rnalib*'s test suite\n",
    "* *large_test_resources*: larger test files that are needed to demonstrate *rnalib* under realistic conditions. \n",
    "\n",
    "*Rnalib* knows about the location of the test data via the package-global __RNALIB_TESTDATA__ variable. This variable can either be set via the \"RNALIB_TESTDATA\" environment variable or by 'monkeypatching' (rna.__RNALIB_TESTDATA__ = ...) as shown below. Once this is done, test data resources can be accessed via `rna.get_resource(<resource_id>)`. *Rnalib* will recreate these resources only if they are not found in the provided output folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "97008543-d575-4613-b6b8-e77a13a1e44c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Running rnalib 0.0.3 on python 3.12.1'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "056da0a77125402f912dec484618f63a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Creating testdata:   0%|          | 0/44 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:root:========= All done with 0 errors  ==========\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'Testdata in rnalib_testdata/'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<textarea rows='4' cols='120'>├── bigfiles\n",
       "│   ├── grch38_chr20.fa.gz.gzi (15.73 KB)\n",
       "│   ├── GRCh38.k24.umap.bed.gz (39.06 MB)\n",
       "│   ├── chess3.0.1.gtf.gz (20.45 MB)\n",
       "│   ├── GRCh38.k24.umap.bed.gz.tbi (488.12 KB)\n",
       "│   ├── human_gene_v2.2.h5 (34.94 GB)\n",
       "│   ├── gencode_39.gff3.gz (72.15 MB)\n",
       "│   ├── hgnc_complete_set.txt (15.8 MB)\n",
       "│   ├── gencode_39.gff3.gz.tbi (298.13 KB)\n",
       "│   ├── grch38_chr20.fa.gz (19.87 MB)\n",
       "│   ├── grch38_chr20.fa.gz.fai (23.0 B)\n",
       "│   └── ...\n",
       "├── bed\n",
       "│   ├── test_bed12.bed.gz.tbi (106.0 B)\n",
       "│   ├── test_bed12.bed.gz (115.0 B)\n",
       "│   ├── test.bed.gz (84.0 B)\n",
       "│   ├── test.bedgraph.gz (208.0 B)\n",
       "│   ├── pybedtools_snps.bed.gz (8.14 MB)\n",
       "│   ├── test_nist.b37_chr20_100kbp_at_10mb.bed (156.0 B)\n",
       "│   ├── dmel_randomvalues.bedgraph.gz (185.33 KB)\n",
       "│   ├── GRCh38.k24.umap.ACTB_ex1+2.bedgraph.gz (1.8 KB)\n",
       "│   ├── dmel_randomvalues.bedgraph.gz.tbi (122.0 B)\n",
       "│   ├── test.bedgraph.gz.tbi (103.0 B)\n",
       "│   └── ...\n",
       "├── .DS_Store (6.0 KB)\n",
       "├── div\n",
       "│   ├── hgnc_complete_set.head.txt.gz (1.6 KB)\n",
       "│   ├── slamseq_example_featureCounts_tc_reads.txt (1.86 KB)\n",
       "│   └── slamseq_example_featureCounts_all_reads.txt (1.84 KB)\n",
       "├── slamseq_example\n",
       "│   ├── slamseq_example_tp0+tc.bam.bai (17.12 KB)\n",
       "│   ├── slamseq_example_tp25+tc_tc-only.bam (63.19 KB)\n",
       "│   ├── slamseq_example_tp10.density.tsv.gz (1.28 KB)\n",
       "│   ├── slamseq_example_tp5+tc_tc-only.bam.bai (17.13 KB)\n",
       "│   ├── slamseq_example_tp0.density.tsv.gz (1.54 KB)\n",
       "│   ├── slamseq_example_tp10+tc_tc-only.bam.bai (17.12 KB)\n",
       "│   ├── slamseq_example_tp15+tc.bam (442.16 KB)\n",
       "│   ├── slamseq_example_tp0+tc.bam (529.67 KB)\n",
       "│   ├── slamseq_example_tp5+tc.bam.bai (17.13 KB)\n",
       "│   ├── slamseq_example_tp20+tc_tc-only.bam.bai (17.12 KB)\n",
       "│   └── ...\n",
       "├── bam\n",
       "│   ├── slamseq_example_5.cr0.05.STAR.bam.bai (17.15 KB)\n",
       "│   ├── NA12878_S1.chr20.10_10p1mb.bam (3.74 MB)\n",
       "│   ├── mapt.NA12156.altex.small.bam (1.23 KB)\n",
       "│   ├── slamseq_example_0.cr0.05.STAR.bam.bai (17.12 KB)\n",
       "│   ├── NA12878_S1.chr20.10_10p1mb.bam.bai (5.34 KB)\n",
       "│   ├── slamseq_example_0.cr0.05.STAR.bam (478.66 KB)\n",
       "│   ├── GSM3185042.Actb.bam (5.92 MB)\n",
       "│   ├── slamseq_example_10.cr0.05.STAR.bam.bai (17.12 KB)\n",
       "│   ├── slamseq_example_20.cr0.05.STAR.bam (412.2 KB)\n",
       "│   ├── GSM3185046.pe_example.vcf.gz (2.86 KB)\n",
       "│   └── ...\n",
       "├── fasta\n",
       "│   ├── ACTB+SOX2.fa.gz.fai (40.0 B)\n",
       "│   ├── dmel_r6.36.fa.gz.gzi (34.98 KB)\n",
       "│   ├── ACTB+SOX2.fa.gz (11.74 KB)\n",
       "│   ├── ACTB+SOX2.fa.gz.gzi (8.0 B)\n",
       "│   ├── dmel_r6.36.fa.gz.fai (68.18 KB)\n",
       "│   └── dmel_r6.36.fa.gz (41.41 MB)\n",
       "├── gff\n",
       "│   ├── pybedtools_gff.gff3.gz.tbi (8.54 KB)\n",
       "│   ├── chess3.GRCh38.ACTB+SOX2.gtf.gz.tbi (336.0 B)\n",
       "│   ├── mirgenedb.dme.sorted.gff3.gz.tbi (1.37 KB)\n",
       "│   ├── gencode_44.ACTB+SOX2.gff3.gz (8.98 KB)\n",
       "│   ├── chess3.GRCh38.ACTB+SOX2.gtf.gz (812.0 B)\n",
       "│   ├── pybedtools_gff.gff3.gz (48.03 KB)\n",
       "│   ├── flybase_dmel-2L-r6.36.gtf.gz.tbi (126.0 B)\n",
       "│   ├── UCSC.hg38.ncbiRefSeq.ACTB+SOX2.sorted.gtf.gz (572.0 B)\n",
       "│   ├── slamseq_example_mm10.gff3.gz (918.0 B)\n",
       "│   ├── chess3.GRCh38.ACTB+SOX2.gff3.gz (839.0 B)\n",
       "│   └── ...\n",
       "├── vcf\n",
       "│   ├── dmelanogaster_6_exported_20230523.vcf.gz.tbi (121.0 B)\n",
       "│   ├── test_snps.vcf.gz (227.0 B)\n",
       "│   ├── test.vcf.gz.tbi (157.0 B)\n",
       "│   ├── dmelanogaster_6_exported_20230523.vcf.gz (15.99 KB)\n",
       "│   ├── test.vcf.gz (478.0 B)\n",
       "│   └── test_snps.vcf.gz.tbi (121.0 B)\n",
       "├── .ipynb_checkpoints\n",
       "└── ...</textarea>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# uncomment commands below to install rnalib's requirements \n",
    "#%pip install -r https://raw.githubusercontent.com/popitsch/rnalib/main/requirements.txt\n",
    "# install additional requirements for this notebook\n",
    "#%pip install biopython seaborn\n",
    "\n",
    "# import required libaries \n",
    "import os, pathlib, platform, logging, sys, tempfile, pytest, dataclasses, random, pprint, time\n",
    "from collections import Counter, defaultdict\n",
    "from tqdm.auto import tqdm\n",
    "import itertools\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pysam\n",
    "import pybedtools \n",
    "import bioframe\n",
    "import biotite.sequence as seq\n",
    "from Bio.Align import PairwiseAligner\n",
    "logging.basicConfig(stream=sys.stderr, level=logging.INFO) # set log level\n",
    "\n",
    "# load rnalib\n",
    "import rnalib as rna\n",
    "from rnalib import gi, GI, SEP, display_textarea\n",
    "display(f\"Running rnalib {rna.__version__} on python {platform.python_version()}\")\n",
    "\n",
    "# ensure test data\n",
    "rna.__RNALIB_TESTDATA__ = \"rnalib_testdata/\" # monkeypatch test data dir\n",
    "rna.testdata.create_testdata(rna.__RNALIB_TESTDATA__, rna.testdata.test_resources) # requires additional tools installed\n",
    "display(f\"Testdata in {rna.__RNALIB_TESTDATA__}\")\n",
    "display_textarea('\\n'.join(rna.dir_tree(pathlib.Path(rna.__RNALIB_TESTDATA__))))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "728ecb24-235c-48f4-aadf-7d13776b5c69",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Genomic intervals\n",
    "At the heart of *rnalib* is its representation of genomic regions (intervals), represented by the  *genomic intervals* (GI) class in *rnalib*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "09f43dbe-2db9-45fd-a35e-838ee7b6df35",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<a id=\"help_icon0.16874711614774018\"><img id=\"help_icon\" style=\"cursor: help; float:left; width: 20px\" src=\"help_icon.png\"/>&nbsp;&nbsp;&nbsp;<b>GI docstring</b></a>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "    document.getElementById(\"help_icon0.16874711614774018\").addEventListener(\"click\", function() {\n",
       "        var win = window.open(\"\", \"\", \"toolbar=no, location=no, directories=no, status=no, menubar=no, scrollbars=yes, resizable=yes, width=1000, height=600\");\n",
       "        win.document.body.innerHTML = \"<h1>GI docstring</h1><p style='font-family:Courier New; font-size: 14px;'><br/>    Genomic intervals (GI) in rnalib are inclusive, continuous and 1-based. This model differs from other<br/>    genomics libraries but was chosen to make interpretation of GIs straightforward: start and end coordinates<br/>    represent the first/last included nucleotide as also seen in a genome browser (such as IGV).<br/><br/>    GIs are implemented as (readonly) namedtuples and can safely be used as keys in a dict.<br/>    They are instantiated by the gi() factory method that either parses from a string representation or by passing<br/>    explicitly chrom/start/stop coordinates. Coordinate components can be &#x27;unrestricted&#x27; by setting them to<br/>    None, e.g., gi(&#x27;chr1&#x27;, 1, None) refers to the whole chromosome 1, gi(&#x27;chr1&#x27;, 100000) refers to the section on<br/>    chr1 from (and including) position 100k on, gi(start=100, end=200) refers to positions 100-200 (inclusive) on<br/>    any chromosome.<br/><br/>    Intervals can be stranded. The strand can be &#x27;+&#x27; or &#x27;-&#x27; or None if unstranded. Note that &#x27;.&#x27; will be converted<br/>    to None. Individual genomic positions (points) are represented by intervals with same start and stop<br/>    coordinate. Empty intervals can be represented by start&gt;end coordinates (e.g., gi(&#x27;chr1&#x27;, 1,0).is_empty() -&gt;<br/>    True) but note that this feature is still experimental and may not be supported in all methods.<br/><br/>    GIs can be compared for equality, overlap, envelopment, etc. and can be sorted by chromosome and coordinates.<br/><br/>    GIs are grouped by chromosomes/contigs of a given reference genome and their order is defined by a<br/>    `RefDict` (that extends regular python dicts). RefDicts are typically instantiated automatically<br/>    from indexing data structures (e.g., .tbi or .bai files) and keep track of the order of chromosomes/contigs<br/>    and their lengths. They are used throughout rnalib to assert the compatibility of different genomic datasets<br/>    and to sort GIs by chromosome and coordinate, e.g., via<br/>    `sorted(gis, key=lambda x: (refdict.index(x.chromosome), x))`. Note that the index of chromosome &#x27;None&#x27;<br/>    is always 0, i.e., it is the first chromosome in the sorted list.<br/><br/>    GIs can be iterated over to yield individual positions (points) as GIs, e.g.,<br/>    tuple(gi(&#x27;chr2&#x27;, 1, 3)) -&gt; (gi(&#x27;chr2&#x27;, 1, 1), gi(&#x27;chr2&#x27;, 2, 2), gi(&#x27;chr2&#x27;, 3, 3)).<br/><br/>    GIs should be instantiated by the `gi()` factory method that allows to configure each coordinate component<br/>    individually or can parse from a string representation, e.g., `gi(&#x27;chr1:1-10 (+)&#x27;)`. The factory method<br/>    conducts type and value checking and conversions (e.g., &#x27;.&#x27; to None for strand).<br/><br/>    Attributes<br/>    ----------<br/>    chromosome: str<br/>        Chromosome (default: None)<br/>    start: int<br/>        First included position, 1-based (default: 0).<br/>    end: int<br/>        Last included position, 1-based (default: MAX_INT)<br/>    strand: str<br/>        Strand, either &#x27;+&#x27; or &#x27;-&#x27; or None if unstranded. Note that &#x27;.&#x27; will be converted to None. Default: None<br/><br/>    Examples<br/>    --------<br/>    &gt;&gt;&gt; gi(&#x27;chr3&#x27;)<br/>    &gt;&gt;&gt; gi(&#x27;chr2&#x27;, 1, 10)<br/>    &gt;&gt;&gt; gi(&#x27;1&#x27;, 1, 10, strand=&#x27;+&#x27;)<br/>    &gt;&gt;&gt; gi(&#x27;chr4:1-10 (-)&#x27;)<br/>    &gt;&gt;&gt; assert gi(&#x27;chrX:2-1&#x27;).is_empty()<br/><br/><br/>    Notes<br/>    -----<br/>    Note that there is also a frozen (immutable) dataclass version of this class, `GI_dataclass` that<br/>    shares most functionality with GIs and is used as superclass for genomic features (genes, transcripts,<br/>    etc.) as instantiated by the `Transcriptome` class. The reason for this separation is that dataclass<br/>    instantiation is relatively slow and has performance impacts when instantiating large sets of GIs.<br/><br/>    </p>\";\n",
       "        win.document.title = \"GI docstring\";\n",
       "        win.document.body.style.background = \"#dcf6fa\";\n",
       "        });\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rna.display_help(rna.GI)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "703a79d8-4da1-43de-8dd6-2a4e47181bf4",
   "metadata": {},
   "source": [
    "Genomic intervals should be instantiated via the `gi()` factory method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5c53432c-0db1-405a-bb08-a41d6d725a31",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[chr2:1-100,\n",
       " chr1:5-500 (+),\n",
       " chr3:50-120 (-),\n",
       " chr1:10-20 (-),\n",
       " None:0-1000 (-),\n",
       " None:10-1000 (-)]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# create a list of genomic intervals\n",
    "locs=[gi('chr2', 1, 100), # chr2:1-100, unstranded\n",
    "      gi('chr1:5-500 (+)'),  # parsed from string, strand='+'\n",
    "      gi('chr3:50-120 (-)'), \n",
    "      gi('chr1:10-20 (-)'), \n",
    "      gi(None, None, 1000, '-'), # chromosome and start unbound, end=1000, strand='-'\n",
    "      gi(None, 10, 1000, '-')] # chromosome unbound\n",
    "display(locs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b679cd3f-8313-47fe-a919-640cb0783bc1",
   "metadata": {},
   "source": [
    "Note that this manually created list of genomic intervals is unsorted and below we show how to properly sort such a list by genomic coordinates and defined \n",
    "chromosome order. Note, that typically you don't need to manually sort genomic intervals that were loaded from genomic datasets as  *rnalib* deals mainly with sorted and indexed files (exceptions are discussed). Respective reference dicts (RefDict) that define chromosome order are loaded from index files. For manually created intervals or interval sets derived via other sources that allow access to (potentially) unsorted files (e.g., rnalib), you can use GI.sort() as shown here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f561a885-505f-437c-a296-afe1bfe7d853",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Sorting w/o reference dict leaves order of chromosomes undefined:'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[None:0-1000 (-),\n",
       " chr2:1-100,\n",
       " chr1:5-500 (+),\n",
       " chr3:50-120 (-),\n",
       " chr1:10-20 (-),\n",
       " None:10-1000 (-)]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr2:1-1 > chr1:1-1 = None'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "RefDict (size: 3): dict_keys(['chr1', 'chr2', 'chr3']), dict_values([None, None, None]) name: test "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Properly sorted list:'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[None:0-1000 (-),\n",
       " None:10-1000 (-),\n",
       " chr1:5-500 (+),\n",
       " chr1:10-20 (-),\n",
       " chr2:1-100,\n",
       " chr3:50-120 (-)]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Properly sorting GIs requires knowledge about the order of the chromosomes:\n",
    "display('Sorting w/o reference dict leaves order of chromosomes undefined:', sorted(locs), SEP)\n",
    "\n",
    "# This list is sorted by start coordinates but not chromosomes. \n",
    "# Generally, the order of intervals from different groups (chromosomes) is left undefined in *rnalib*: \n",
    "display(f\"chr2:1-1 > chr1:1-1 = {gi('chr2', 1, 1) > gi('chr1', 1, 1)}\", SEP)\n",
    "\n",
    "# To sort also by chromosome, you can use a reference dict which defines the chromosome order:\n",
    "refdict=rna.RefDict({'chr1':None, 'chr2':None, 'chr3':None}, 'test', None) # manually created RefDict w/o chrom lengths\n",
    "display(refdict, SEP)\n",
    "\n",
    "# And then call the GI.sort(...) method which is basically returning `sorted(locs, key=lambda x: (refdict.index(x.chromosome), x))`\n",
    "# Intervals with no defined chromosome (i.e., chromosome=None) will always be at the beginning of the list (their index is -1). \n",
    "display('Properly sorted list:', GI.sort(locs, refdict))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4552a26-4e84-468e-aa19-3e3a99b94f2e",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Transcriptome implementation\n",
    "*Rnalib* provides a `transcriptome` implementation that enables convenient and efficient access/filtering of genomic annotations via Python. \n",
    "Briefly, a (filtered) input GTF/GFF file is parsed and respective (frozen) dataclasses that derive from the generic `Feature` class are instantiated. \n",
    "The implementation keeps track of parent/child relationships (e.g., between genes, transcripts and their exons) and enables efficient iteration, querying and annotation. Data fields and annotations can conveniently be accessed via 'dot' notation (i.e., `<feature>.<attribute>`). \n",
    "\n",
    "Note that *rnalib* does not implement the full GFF3 format as specified [here](https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md)\n",
    "but rather supports various popular gff 'flavours' as used in popular annotation sets such as _encode_, _ensembl_, _ucsc_, _chess_, _mirgenedb_ and _flybase_.\n",
    "It also supports a 'generic' GFF3 format where all parsed features are interpreted as 'transcript' objects irrespecting the actual feature_type column.\n",
    "Respective parent 'gene' annotations are automatically added, no sub-features (exons, introns, etc.) are created.\n",
    "\n",
    "Here is a brief description of the datamodel depicted below:\n",
    "-   Model contains `genes`, `transcripts` and arbitrary sub-features (e.g., `exons`, `introns`, 3'/5'-`UTR`s, `CDS`) as defined\n",
    "    in the GFF file. Note that `intron` annotations are calculated automatically from exon annotations unless configured otherwise. \n",
    "    Frozen dataclasses (derived from the 'Feature' class) are created for all parsed feature types automatically and users may configure which \n",
    "    GTF/GFF attributes will be parsed and added to those dataclasses as fields.\n",
    "    Gene intervals must envelop transcript intervals which must envelop sub-feature intervals.\n",
    "    \n",
    "-   A transcriptome maintains an `anno` dict that maps (frozen) features to dicts of arbitrary annotation values which supports the (incremental) annotation of annotation features. Annotation values can be accessed like GFF-parsed fields via 'dot' notation: `<feature>.<attribute>`. \n",
    "    \n",
    "-   The implementation exploits the hierarchical relationship between genes and their sub-features to optimize storage and computational requirements where possible.        \n",
    "    Genomic sequences, for example, can be loaded from a reference genome via `load_sequence_data()` but will  be explicitly stored only for gene features. For sub-features (e.g., transcripts, exons, etc.) the respective sequence will be sliced from the gene sequence on demand. \n",
    "    \n",
    "-   The `get_sequence()` method supports several modes:\n",
    "    If `mode='rna'` is passed, the sequence is returned in 5'-3' orientation, i.e., they are reverse-complemented\n",
    "    for minus-strand transcripts. The returned sequence will, however, still use the DNA alphabet (ACTG) to\n",
    "    enable direct alignment/comparison with genomic sequences.\n",
    "    if `mode='spliced'`, the spliced 5'-3' sequence will be returned.\n",
    "    if `mode='translated'`, the spliced 5'-3' CDS sequence will be returned.\n",
    "    \n",
    "-   Genomic range queries via `query()` are supported by a combination of interval and linear search queries.\n",
    "    A transcriptome object maintains one intervaltree per chromosome built from gene annotations.\n",
    "    Overlap/envelop queries will first be applied to the respective intervaltree and the (typically small\n",
    "    result sets) will then be filtered, e.g., for requested sub-feature types.\n",
    "    \n",
    "-   Transcriptome models can also be built from subsets of GFF/GTF files by filtering contained transcripts with a `TranscriptFilter()` that can restrict, e.g., included chromosomes, genomic regions, transcript_ids, gene_types or transcripts with certain 'tag' values (as used in gencode).\n",
    "\n",
    "-   A transcriptome object keeps track of the numbers of filtered GFF entries in a `log` Counter object that is, e.g., useful for debugging. The 'structure' of the parsed features can be accessed via the `get_struct()` method.\n",
    "\n",
    "The following cells demonstrate some use cases."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b09a98f-f140-47d7-acba-f0897b01dbad",
   "metadata": {},
   "source": [
    "![transcriptome datamodel](../docs/_static/rnalib_transcriptome.png \"Transcriptome datamodel\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d91fd797-9af1-4749-bc22-63da3910fe86",
   "metadata": {},
   "source": [
    "Here is a short introduction of how to use *rnalib*'s transcriptome implementation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c6c221a8-a157-4ec5-a47f-850da59dd455",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "    <!-- gh-pages -->\n",
       "    <link rel='stylesheet' href='https://unpkg.com/freezeframe@3.0.10/build/css/freezeframe_styles.min.css'>\n",
       "    <script type='text/javascript' src='https://unpkg.com/freezeframe@3.0.10/build/js/freezeframe.pkgd.min.js'></script>\n",
       "\n",
       "    <script type='text/javascript'> \n",
       "      third = new freezeframe('.my_class_3').capture().setup();\n",
       "      $(function() {\n",
       "        $('.start').click(function(e) {\n",
       "          e.preventDefault();\n",
       "          third.trigger();\n",
       "        });\n",
       "      })\n",
       "    </script>\n",
       "    \n",
       "    <img class=\"my_class_3 freezeframe-responsive\" src=\"../docs/_static/screencasts/introduction.gif\" />\n",
       "    <button class=\"start\">restart</button>\n",
       "    \n",
       "    <script>setTimeout( function() { third.trigger(); }, 1000);</script>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rna.display_animated_gif(\"../docs/_static/screencasts/introduction.gif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "92b69487-7b48-485b-9010-656266ef4e9e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<a id=\"help_icon0.43732405892628956\"><img id=\"help_icon\" style=\"cursor: help; float:left; width: 20px\" src=\"help_icon.png\"/>&nbsp;&nbsp;&nbsp;<b>Transcriptome docstring</b></a>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "    document.getElementById(\"help_icon0.43732405892628956\").addEventListener(\"click\", function() {\n",
       "        var win = window.open(\"\", \"\", \"toolbar=no, location=no, directories=no, status=no, menubar=no, scrollbars=yes, resizable=yes, width=1000, height=600\");\n",
       "        win.document.body.innerHTML = \"<h1>Transcriptome docstring</h1><p style='font-family:Courier New; font-size: 14px;'><br/>    Represents a transcriptome as modelled by a GTF/GFF file.<br/><br/>    *   Model contains genes, transcripts and arbitrary sub-features (e.g., exons, intron, 3&#x27;/5&#x27;-UTRs, CDS) as<br/>        defined in the GFF file. Frozen dataclasses (derived from the &#x27;Feature&#x27; class) are created for all parsed<br/>        feature types automatically and users may configure which GTF/GFF attributes will be added to those (and are<br/>        thus accessible via dot notation, e.g., gene.gene_type).<br/>    *   This `transcriptome` implementation exploits the hierarchical relationship between genes and their<br/>        sub-features to optimize storage and computational requirements, see the `Feature` documentation for<br/>        examples. To enable this, however, parent features *must* envelop (i.e., completely contain) child feature<br/>        locations and this requirement is asserted when building the `transcriptome`.<br/>    *   A `transcriptome` maintains an `anno` dict mapping (frozen) features to dicts of arbitrary annotation<br/>        values. This supports incremental and flexible annotation of `transcriptome` features. Values can directly<br/>        be accessed via dot notation &lt;feature&gt;.&lt;attribute&gt; and can be stored/loaded to/from a (pickled) file.<br/>    *   `Feature` sequences can be added via `load_sequences()` which will extract the sequence of the top-level<br/>        feature (&#x27;gene&#x27;) from the configured reference genome. Sequences can then be accessed via get_sequence().<br/>        For sub-features (e.g., transcripts, exons, etc.) the respective sequence will be sliced from the gene<br/>        sequence. If mode=&#x27;rna&#x27; is passed, the sequence is returned in 5&#x27;-3&#x27; orientation, i.e., they are<br/>        reverse-complemented for minus-strand transcripts. The returned sequence will, however, still use the DNA<br/>        alphabet (ACTG) to enable direct alignment/comparison with genomic sequences.<br/>        if mode=&#x27;spliced&#x27;, the spliced 5&#x27;-3&#x27; sequence will be returned.<br/>        if mode=&#x27;translated&#x27;, the spliced 5&#x27;-3&#x27; CDS sequence will be returned.<br/>    *   Genomic range queries via `query()` are supported by a combination of interval and linear search queries.<br/>        A transcriptome object maintains one intervaltree per chromosome built from gene annotations.<br/>        Overlap/envelop queries will first be applied to the respective intervaltree and the (typically small<br/>        result sets) will then be filtered, e.g., for requested sub-feature types.<br/>    *   When building a transcriptome model from a GFF/GTF file, contained transcripts can be filtered using a<br/>        :func:`TranscriptFilter &lt;rnalib.TranscriptFilter&gt;`.<br/>    * | The current implementation does not implement the full GFF3 format as specified in<br/>      | https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md<br/>      | but currently supports various popular gff &#x27;flavours&#x27; as published in<br/>      | encode, ensembl, ucsc, chess, mirgendb and flybase databases (see<br/>      | :func:`GFF_FLAVOURS &lt;rnalib.constants.GFF_FLAVOURS&gt;`). As such this implementation will likely be extended<br/>      | in the future.<br/><br/>    @see the `README.ipynb &lt;https://github.com/popitsch/rnalib/blob/main/notebooks/README.ipynb&gt;`_ jupyter<br/>    notebook for various querying and iteration examples<br/><br/>    Parameters<br/>    ----------<br/>    annotation_gff : str<br/>        Path to a GFF/GTF file containing the transcriptome annotations.<br/>    annotation_flavour : str<br/>        The annotation flavour of the GFF/GTF file. Currently supported: encode, ensembl, ucsc, chess, mirgendb and<br/>        flybase.<br/>    genome_fa : str<br/>        Path to a FASTA file containing the reference genome. If provided, sequences will be loaded and can be<br/>        accessed via get_sequence().<br/>    gene_name_alias_file : str<br/>        Path to a gene name alias file. If provided, gene symbols will be normalized using the alias file.<br/>    annotation_fun_alias : Callable<br/>        Name of a function that will be used to alias gene names.<br/>        The function must be defined in the global namespace and must accept a gene name and return the alias.<br/>        If provided, gene symbols will be normalized using the alias function.<br/>    copied_fields : tuple<br/>        List of GFF/GTF fields that will be copied to the respective feature annotations.<br/>    calc_introns : bool<br/>        If true, intron features will be added to transcripts that do not have them.<br/>    load_sequence_data : bool<br/>        If true, sequences will be loaded from the genome_fa file.<br/>    disable_progressbar : bool<br/>        If true, progressbars will be disabled.<br/>    genome_offsets : dict<br/>        A dict mapping chromosome names to offsets. If provided, the offset will be added to all coordinates.<br/>    name : str<br/>        An (optional) human-readable name of the transcriptome object. default: &#x27;Transcriptome&#x27;<br/>    feature_filter : TranscriptFilter<br/>        A TranscriptFilter object that will be used to filter transcripts.<br/><br/><br/>    Attributes<br/>    ----------<br/>    genes : List[Gene]<br/>        List of genes in the transcriptome.<br/>    transcripts : List[Transcript]<br/>        List of transcripts in the transcriptome.<br/>    anno : Dict[Feature, Dict[str, Any]]<br/>        Dictionary mapping features (e.g., genes, transcripts) to their annotations.<br/><br/>    Examples<br/>    --------<br/>    &gt;&gt;&gt; fb_gff = ... # path to flybase gff file<br/>    &gt;&gt;&gt; t = rna.Transcriptome(annotation_gff = fb_gff, annotation_flavour =&#x27;flybase&#x27;)<br/>    &gt;&gt;&gt; # download full flybase GTF to temporary directory, use the create_testdata function to download, sort and index<br/>    &gt;&gt;&gt; with tempfile.TemporaryDirectory() as tmpdir: # noqa<br/>    &gt;&gt;&gt;     res = {&quot;flybase_full&quot;: {<br/>    &gt;&gt;&gt;            &quot;uri&quot;: &quot;https://ftp.flybase.net/releases/current/dmel_r6.55/gtf/dmel-all-r6.55.gtf.gz&quot;,<br/>    &gt;&gt;&gt;            &quot;filename&quot;: &quot;dmel-all-r6.55.gtf.gz&quot;}}<br/>    &gt;&gt;&gt;    rna.testdata.create_testdata(tmpdir, res) # download, sort + index; needs bgzip, bedtools, tabix<br/>    &gt;&gt;&gt;    t2 = rna.Transcriptome(annotation_gff = rna.get_resource(&quot;flybase_full&quot;, tmpdir, res), annotation_flavour =&#x27;flybase&#x27;)<br/>    &gt;&gt;&gt;  # instantiate Ensembl with chrom aliasing using a sliced genome FASTA file<br/>    &gt;&gt;&gt; config = {<br/>    &gt;&gt;&gt;    &#x27;genome_fa&#x27;: get_resource(&#x27;ACTB+SOX2_genome&#x27;),  # a sliced genome FASTA file; part of rnalib testdata<br/>    &gt;&gt;&gt;    &#x27;genome_offsets&#x27;: {&#x27;chr3&#x27;: 181711825, &#x27;chr7&#x27;: 5526309}, # genome offsets for the sliced genome<br/>    &gt;&gt;&gt;    &#x27;annotation_gff&#x27;: get_resource(&#x27;ensembl_gff&#x27;), # ensembl GFF file; part of rnalib testdata<br/>    &gt;&gt;&gt;    &#x27;annotation_flavour&#x27;: &#x27;ensembl&#x27;, # ensembl annotation flavour<br/>    &gt;&gt;&gt;    &#x27;annotation_fun_alias&#x27;: &#x27;toggle_chr&#x27;, # aliasing function for chromosome names; will all/remove &#x27;chr&#x27; prefix<br/>    &gt;&gt;&gt;    &#x27;load_sequence_data&#x27;: True # load sequences from genome_fa<br/>    &gt;&gt;&gt;    }<br/>    &gt;&gt;&gt; t3 = Transcriptome(**config)  # instantiate transcriptome with config dict<br/>    </p>\";\n",
       "        win.document.title = \"Transcriptome docstring\";\n",
       "        win.document.body.style.background = \"#dcf6fa\";\n",
       "        });\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rna.display_help(rna.Transcriptome)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e6d1a5e-bd14-44f3-bf51-6ec86c04f82f",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Basic access examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0b23c8c2-a8cf-43b3-959e-e09ef620ca5f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Transcriptome with 14 genes and 41 tx (+seq)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Counter({'parsed_gff_lines': 749})"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# First, let's create a transcriptome containing dmel chromosome 2L annotation from a flybase annotation file.\n",
    "t = rna.Transcriptome(\n",
    "        genome_fa=rna.get_resource('dmel_genome'),      # genome FASTA\n",
    "        annotation_gff=rna.get_resource('flybase_gtf'), # gene annotation GTF/GFF file\n",
    "        annotation_flavour='flybase',                   # flavour of the annotation file.\n",
    "        load_sequence_data=True,                        # load sequences from configured genome FASTA file\n",
    "        disable_progressbar=True,                       # no progressbars\n",
    "        feature_filter={'location': { 'included': { 'chromosomes': ['2L'] } }} # a simple filter that will include only annotations from 2L\n",
    "    ) \n",
    "# show some stats/debugging information\n",
    "display(t) # show some basic info (#genes, #tx)\n",
    "display(t.log) # show stats (such as parsed and filtered lines)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72302fdf-9b16-4b12-85a0-2afc12fc8269",
   "metadata": {},
   "source": [
    "Now, we can access some data. The transcriptome object maintains a list of its gene and transcript objects and enables direct access via featur_ids (gene/transcript ids) or gene names (symbols)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "023aa938-34fa-4bbf-b31b-12df291e5df0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"Gene names: ['CR11023', 'l(2)gl', 'Ir21a', 'asRNA:CR43609', 'Cda5', 'lncRNA:CR46254', 'lncRNA:CR45339', 'lncRNA:CR45340', 'dbr', 'asRNA:CR44987', 'galectin', 'CG11374', 'net', 'Zir']\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"Transcripts of gene 'l(2)gl': ['FBtr0306591', 'FBtr0078169', 'FBtr0078168', 'FBtr0078170', 'FBtr0306592', 'FBtr0306590', 'FBtr0078167', 'FBtr0330655', 'FBtr0078171', 'FBtr0078166', 'FBtr0306589'].\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"Number of exons per Cda5 tx: {'FBtr0309227': 12, 'FBtr0309228': 14, 'FBtr0309226': 14, 'FBtr0309229': 15, 'FBtr0078163': 13, 'FBtr0078164': 14, 'FBtr0113416': 13, 'FBtr0309225': 11}\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"The parent gene of transcript 'FBtr0330655' is l(2)gl\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# show gene names that were parsed from the flybase 'gene_symbol' field\n",
    "display(f\"Gene names: {[g.gene_name for g in t.genes]}\",SEP) \n",
    "\n",
    "# query gene by name and get ids of its transcripts\n",
    "display(f\"Transcripts of gene 'l(2)gl': {[tx.feature_id for tx in t.gene['l(2)gl'].transcript]}.\",SEP) \n",
    "\n",
    "# get number of exons per transcripts\n",
    "display(f\"Number of exons per Cda5 tx: { {tx.feature_id:len(tx.exon) for tx in t.gene['Cda5'].transcript} }\",SEP) \n",
    "\n",
    "# Access transcript by id and show the name field of its gene (=parent)\n",
    "display(f\"The parent gene of transcript 'FBtr0330655' is {t.transcript['FBtr0330655'].parent.gene_name}\",SEP) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46919a79-ae40-4cb5-988b-004c7f581548",
   "metadata": {},
   "source": [
    "Here are some more complex access examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0b8c1364-a2fe-4f94-923e-a19013a38b2b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'The gene l(2)gl is located on chromosome 2L at 2L:9839-21376 (-) '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<textarea rows='7' cols='120'>{'CDS': (CDS@2L:17053-17136,\n",
       "         CDS@2L:14933-15711,\n",
       "         CDS@2L:13683-14874,\n",
       "         CDS@2L:13520-13625,\n",
       "         CDS@2L:12286-12928,\n",
       "         CDS@2L:11779-12221,\n",
       "         CDS@2L:11410-11518,\n",
       "         CDS@2L:11218-11344),\n",
       " 'exon': (exon@2L:21066-21376,\n",
       "          exon@2L:19880-20020,\n",
       "          exon@2L:17053-17212,\n",
       "          exon@2L:14933-15711,\n",
       "          exon@2L:13683-14874,\n",
       "          exon@2L:13520-13625,\n",
       "          exon@2L:12286-12928,\n",
       "          exon@2L:11779-12221,\n",
       "          exon@2L:11410-11518,\n",
       "          exon@2L:9839-11344),\n",
       " 'feature_id': 'FBtr0306591',\n",
       " 'feature_type': 'transcript',\n",
       " 'five_prime_UTR': (five_prime_UTR@2L:21066-21376,\n",
       "                    five_prime_UTR@2L:19880-20020,\n",
       "                    five_prime_UTR@2L:17137-17212),\n",
       " 'gene_type': None,\n",
       " 'gff_feature_type': 'mRNA',\n",
       " 'intron': (intron@2L:20021-21065,\n",
       "            intron@2L:17213-19879,\n",
       "            intron@2L:15712-17052,\n",
       "            intron@2L:14875-14932,\n",
       "            intron@2L:13626-13682,\n",
       "            intron@2L:12929-13519,\n",
       "            intron@2L:12222-12285,\n",
       "            intron@2L:11519-11778,\n",
       "            intron@2L:11345-11409),\n",
       " 'parent': gene@2L:9839-21376,\n",
       " 'source': 'FlyBase',\n",
       " 'subfeature_types': ('five_prime_UTR',\n",
       "                      'intron',\n",
       "                      'CDS',\n",
       "                      'three_prime_UTR',\n",
       "                      'exon'),\n",
       " 'three_prime_UTR': (three_prime_UTR@2L:9839-11214,),\n",
       " 'transcriptome': Transcriptome with 14 genes and 41 tx (+seq)}</textarea>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'transcript_id: FBtr0306591, location: 2L:9839-21376 (-), 1st exon rnk: 1'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Genes can be accessed by name or gene id; show gene name, chromosome and location\n",
    "display(f\"The gene {t['l(2)gl'].gene_name} is located on chromosome {t['l(2)gl'].chromosome} at {t['FBgn0002121'].location} \", SEP) \n",
    "\n",
    "# you can inspect the fields of a transcriptome feature with regular python methods, e.g. vars()\n",
    "# Here we show, e.g., the structure of the 1st transcript\n",
    "tx = t.gene['l(2)gl'].transcript[0]\n",
    "display_textarea(pprint.pformat(vars(tx)), rows=7)\n",
    "\n",
    "# fields as well as dynamic annotations (see below) can be accessed by dot (<feature>.<attribute/field>) notation\n",
    "# this includes calculated annotations such as location and rnk (exon/intron number)\n",
    "display(SEP, f\"transcript_id: {tx.feature_id}, location: {tx.location}, 1st exon rnk: {tx.exon[0].rnk}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0c8efd8-16fc-4a3a-a448-a7216a83c6f2",
   "metadata": {},
   "source": [
    "## Parent/child relationships\n",
    "The hierarchical structure of transcriptome features (i.e., genes, transcripts, exonc, etc.) is modeled via the *parent* attrbute.\n",
    "You can inspect the instantiated structure of the transcriptome by calling the `get_struct()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "ecf4c9c8-284a-46fe-8564-ecb5dfa71469",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'gene': {'transcript'},\n",
       " 'transcript': {'CDS', 'exon', 'five_prime_UTR', 'intron', 'three_prime_UTR'},\n",
       " 'intron': set(),\n",
       " 'exon': set(),\n",
       " 'five_prime_UTR': set(),\n",
       " 'CDS': set(),\n",
       " 'three_prime_UTR': set()}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(t.get_struct()) # Show the hierarchical structure encoded in the GFF"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9aeb55f0-08e4-4316-a68d-dcc4173f7d3d",
   "metadata": {},
   "source": [
    "You can access a gene's transcripts via `<gene>.transcript`, a transcripts exons via `<transcript>.exon` aso."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "a378267e-168b-4e4f-90ff-365e7bf82aea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Selected gene gene@2L:25402-65404 with name Cda5'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Selected tx transcript@2L:25402-65404 with id FBgn0051973'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Selected exon exon@2L:34558-34604 with rank 7'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "g = random.choice(t.genes) # get a random gene from the transcriptome\n",
    "display(f\"Selected gene {g} with name {g.gene_name}\")\n",
    "tx = random.choice(g.transcript) # select a random transcript of this gene\n",
    "display(f\"Selected tx {tx} with id {g.feature_id}\")\n",
    "exon = random.choice(tx.exon) # select a random exon of this transcript\n",
    "display(f\"Selected exon {exon} with rank {exon.rnk}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a1d8f2e-601c-4d5a-998b-88caf6aa4bf9",
   "metadata": {},
   "source": [
    "And you can also traverse the hierachy in child->parent direction via `<transcript>.gene`, aso."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "74e49c7b-ebf2-4e53-93c2-86440df63e7f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'The gene name of the associated gene is Cda5'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(f\"The gene name of the associated gene is {exon.transcript.gene.gene_name}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7b67ea4-326f-4a56-9898-e1721fa601b0",
   "metadata": {},
   "source": [
    "## Sequence access\n",
    "Genomic sequences of transcriptome features are an example for dynamically calculated annotations. \n",
    "First, gene sequences must be loaded from the configured genome FASTA file with `load_sequence_data()` which is done automatically if configured. Sequence data of child features (transcripts, exons, etc.) will be sliced from these sequences on demand.\n",
    "\n",
    "Sequences can be accessed via the `t.get_sequence()` method which supports various sequence conversion modes (e.g., 'dna', 'rna', 'spliced', 'translated', etc.).\n",
    "Sequences can also be accessed via dot notation, i.e., via `<feature>.sequence`, `<feature>.spliced_sequence` or `<feature>.translated_sequence`, which then calls `get_sequence()` with default parameters.\n",
    "\n",
    "Here, we demonstrate some of this functionality by accessing sequence data for a transcript of the 'l(2)gl (FBgn0002121)' gene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2132eed7-d788-42e9-9565-fac7ef977414",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Selected tx: FBtr0078171'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<textarea rows='4' cols='120'>DNA sequence of 1st exon: CTTTATGTATTTCCGTTTGTACGGTTAAGGCGGTGGCCGAGTAATTTTTTGAACTATTTTATTTGCTCACCAACACGCATAAGTTTGTATGTATTCTTATGTTGAGAAACGCCCGCACTGTTCTTCACATCGTTGCTTTAAATGAGAATAAATGAAATAAACGTGTCAACCAAGAACACACAGATAATAATACGACCGGTCAGAGCTTAACCGTGTATCTGTTTTATAAAACGTGAACAATATTATTAGCCCAAAACGATATGCGCGTCATTTTAACACACACAAAGTCGCGATCGGTGGGTCTAGTGTG</textarea>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<textarea rows='4' cols='120'>Spliced sequence: CACACTAGACCCACCGATCGCGACTTTGTGTGTGTTAAAATGACGCGCATATCGTTTTGGGCTAATAATATTGTTCACGTTTTATAAAACAGATACACGGTTAAGCTCTGACCGGTCGTATTATTATCTGTGTGTTCTTGGTTGACACGTTTATTTCATTTATTCTCATTTAAAGCAACGATGTGAAGAACAGTGCGGGCGTTTCTCAACATAAGAATACATACAAACTTATGCGTGTTGGTGAGCAAATAAAATAGTTCAAAAAATTACTCGGCCACCGCCTTAACCGTACAAACGGAAATACATAAAGCATAGAGGCGTAGTCGGTGAGTCAGCCACGCCAACAGTTAACAGTCGTCATTTTGCTGTCTGTCTAACACGCAATTTGTGTGTGTATATGTGAATAAACAGATACCGGTATGCAATCGCAAAGCTGGTCAACTGTAAAGATCGTTTAGAATAAATACGGCCACAATTAGTATACATATTTATTCATCGAAACCAACCCACCACAAGAGCATAACCAATTATGTTAAAGTTTATCAGAGGAAAAGGGCAGCAGCCCAGTGCTGACAGACACCGCCTACAGAAGGACCTTTTTGCTTATCGTAAGACGGCACAGCATGGCTTTCCTCATAAGCCTTCGGCTCTTGCGTATGATCCAGTTTTGAAACTTATGGCAATAGGGACGCAAACAGGGGCTTTAAAAGTTTTCGGTCAACCCGGAGTTGAATTGTACGGTCAGCATACTTTGTTAAACAATTCAGCATCGGAGCTTAATGTACAATTACTTGAATGGGTGTATGGAACTGGTCGCATACTTTCGTTGACGGCAGCGAATCAATTAATTCTATGGGAGCCAGTTGGAGCAACGTTGCTGCCAATCAAAACACTACCGTTTGACGGCAAACTTAAAAAAGTTTCATCGCTGTGCTGTTCTCTCAGTAAGGATCTGCTATGGATTGGAACAGAAGGTGGAAACATCTATCAACTGGATTTACATACATTTACCATTAAGGAGCCTGTAATTTACCATGACGTTGTGCTAGAGCAGGTGCCACCAGCCTACAAGCTAAATCCTGGTGCAATTGAGTCAATCCGCCAACTTCCAAACTCCCCTAGCAAACTTCTAGTTGCATACAATCGCGGCCTTTGTGTTTTGTGGGATTTTGAAAGCGCATCTGTCCAGCGAGCATACATAGCCCCTGGACATGGACAGAGCGTTGGTCTTACAGTGAACTTCGAAGGATCTGAATTTACCTGGTACCACGCTGATGGTTCATACGCCACTTGGAGCATAGATAACCCAGAACCGCCGTCGAATGTTAATTATGTGCCTTATGGACCTGATCCATGCAAAAGCATAAATCGACTGTACAAAGGCAAGCGAAGATCCAACGATGTAATTGTTTTTTCCGGCGGCATGCCACGGTCAGCATATGGTGATCACAATTGTGTGTCCGTTCACGCCAGCGATGGACACAAAGTGTGTCTTGACTTTACGTCTAAAGTGATTGACTTTTTTGTGACCTTTGAAAATAATAGAGATGTCGCTGAAGTTCTTGTTGTACTACTTGAAGAGGAACTCTGCGCTTACGATCTTACTGACCCTAATATTTGTGCTATCAAAGCGCCATATCTTCACTCTGTCCATGCATCAGCTGTCACTTGCAATTACCTTGCTTCTGAAGTCGTACAGTCGGTATATGAAAGTATTTTAAGAGCTGGAGATGAACAAGACATTGACTATAGCAATATTAGCTGGCCTATCACTGGCGGTACTCTCCCGGATAACTTAGAAGAATCTGTAGAAGAGGACGCGACTAAGCTTTATGAGATTTTGTTAACTGGTCACGAAGATGGTTCTGTTAAATTTTGGGACTGCACTGGAGTGTTGCTTAAACCAATTTATAATTTTAAAACTAGCAGCATTTTTGGAAGTGAGTCAGACTTCCGAGATGACGCAGCTGCAGATATGAGTGCGGAACAAGTCGATGAAGGAGAACCGCCATTTCGGAAATCAGGACTTTTTGATCCTTATTCAGATGACCCTCGTTTAGCAGTGAAGAAAATAGCATTCTGCCCAAAAACCGGACAACTTATTGTTGGTGGCACAGCGGGCCAAATAGTTATAGCCGACTTCATAGACTTACCCGAAAAAGTGTCTTTAAAATACATTTCAATGAATTTGGTCAGCGATCGTGATGGATTTGTGTGGAAGGGTCACGATCAGTTAAACGTGCGATCGAACTTATTAGACGGAGAAGCAATTCCTACGACGGAACGTGGTGTAAATATATCGGGAGTACTGCAAGTTTTGCCGCCAGCCAGCATAACATGCATGGCACTCGAAGCAAGCTGGGGCCTAGTATCTGGTGGGACTGCGCACGGCTTAGTTCTCTTTGACTTCAAAAACTTTGTTCCAGTATTTCATCGCTGCACTTTAAACCCAAATGATCTTACTGGAGCAGGAGAGCAGCTGTCTCGTCGAAAGTCTTTTAAGAAATCATTGAGGGAGTCATTTAGAAAGCTTCGCAAGGGTCGATCGACCAGGACCAACCAGAGCAATCAAGTACCAACAACGCTGGAAGCAAGACCCGTCGAGAGGCAAATAGAGGCTCGTTGTGCAGATGACGGGCTAGGATCCATGGTGCGATGTTTACTATTTGCCAAAACTTATGTTACTAATGTCAACATAACGTCGCCAACTTTGTGGTCAGCAACAAATGCCAGTACAGTCTCGGTTTTCCTTCTGCATTTGCCACCAGCGCAGACCGCGGCAACTGCCGTCCCGTCGGCAAGTGGCAATGCACCACCACACATGCCCCGCCGAATTTCTGCGCAGCTTGCTAAAGAAATACAATTAAAACATCGTGCTCCTGTGGTGGGTATTTCTATTTTTGATCAGGCGGGTAGCCCTGTCGATCAGCTGAACGCCGGTGAAAACGGGAGTCCACCGCATCGTGTACTTATTGCTTCCGAGGAACAGTTCAAGGTGTTTTCACTTCCGCAACTAAAGCCGATTAACAAATATAAGCTTACCGCTAACGAAGGTGCTCGGATTCGCCGCATCCATTTTGGTTCGTTTAGTTGTCGCATATCCCCGGAAACACTGCAGAGTATGCACGGTTGTAGCCCAACTAAGTCCACGCGTTCACATGGCGATGGAGAGGCGGATCCTAATATCAGTGGAAGCTTGGCTGTAAGTCGTGGAGATGTATATAACGAAACAGCATTGATATGTTTAACGAATATGGGCGATATCATGGTTTTATCAGTACCTGAATTAAAAAGACAGCTGAATGCCGCAGCAGTGCGACGGGAAGACATTAATGGAGTTTCGTCACTTTGCTTTACAAACTCTGGAGAAGCACTGTATATGATGTCTTCTTCTGAACTGCAGCGTATTGCTTTAGCCACGTCCAGAGTCGTGCAACCCACTGGCGTTGTTCCAGTAGAACCATTAGAAAATGAAGAGTCTGTGTTGGAAGAAAATGATGCAGAGAATAATAAGGAAACCTACGCATGTGATGAAGTTGTGAATACATATGAAATTAAAAATCCATCAGGCATTTCAATATGCACAAGGCCTGCAGAGGAAAACGTTGGAAGAAATAGTGTTCAGCAAGTTAATGGAGTCAACATTTCAAATTCACCTAATCAAGCTAACGAGACTATCAGCAGCTCTATTGGCGATATTACCGTTGACTCGGTGCGCGACCATTTAAATATGACGACCACCACTTTGTGTTCTATTAATACAGAGGAAACCATTGGTCGCCTATCTGTACTTAGCACGCAAACCAACAAAGCCAGTACTACCGTAAACATGAGTGAAATTCCAAATATTAATATTTCTAATTTAGAGGACTTGGAATCGAAAAGAAATACGACGGAAACGAGTACTAGTTCTGTTGTAATTAAATCTATAATTACAAACATTTCTCATGAAAAAACGAACGGAGACAACAAAATAGGAACGCCAAAAACAGCGCCTGAAGAAAGCCAATTTTAACATTGACAGAAGCCGTAACCTACTAATTATTTTATACCTTTACAGAAATACGACGCAAACGAGTACTAGTCCTGATGTAATTTAATCTATAATTACAAACATGTCTCACTAAAAATCGAACGGGGACAGCAAAATAGGAACGCAAAAAAATTAAATTAAATTAAACGCGATTTACATACAAACAGAAATGACAGAATGATAATATAAAATATTCATTTTTTATTTGGCTTAAGCGATGTTGTTGTTCCAAAACCATATAATTATTTATGATTTTATGTAATGTTTTCATGTATTTTCGCGTAGTGACTTTATACCCTTTACCCGTATGATGAAACGAAACATGATATTTGTAGAAAGGTACAAAGGATAAAAGTTGACATTTAAGGCCGATTGTAAAATATTAAAACGCAACCTGGATAACATGCCGATTAATTACTTCTGTTCGTCCGTATAAACACAAATACCGATAAATTGGTCCGTCCCCATTTAATAAATTGTATATTAAAATGGTCTTTATTTTGGATAAATAATTCAAATTATATAGGAACATTTGGGTTTTAAGGGATAGAAACAGCGCTGGTGTCTTTCATTTCCGTTCGTATGACCGTTCAAAGATTATAAAAAAATTGTCGTGCCCACATTTTTTAAATTTCTTTTATTTTCGGTTAAAATATTGTGATACAAAAAAATGCCCTTACAAAATTAGGACAACAACAATTTGCTTTATTATGATGAACAAATTAGGAGTACAAAAAATAGCCTGATGCAAGGATGGGGAGAAAAAAGACCTTGAAGTTCAAAAGAATTTGTTGTATCAAATCCTCAATAAGAAAACTACTTAAAATTAAAAGCGTGTTTCAAATTAATTGTAGAATATTTCGTAAAATGAAATTGCAAACTTTTTTAGCACGAAATTACATGTATTTAATTCTTTGGCTAATACGTCAGTTAAAAAAGGCTACATTTTGATGATACTCACAAATGTAATTTTTGCATTCTACACAGGTGCAATTAAAAAATTATGGTTTTTATGAGTTTTAGAGGGAAGAAAATATGTCAATATGAGAAGTAAGTTAATGAGCCTTCAAGTACTTTAGTACATAATGTAATTTAATTTTAGAAATCCTACGGCTTTATATTCAATTGTAAAGTATTTTAAAATTCCTGACGATTTGTAGTAATTGTTTTGTTCTTAGCTTTTTTGCATTCTCCTCCGTTACTCCGTTACTCGTATAGTAAAAGGTATAGTATATACAATAGTATTGTAGATTCGTTGAGAAGTATGTTACAGGCAGAAGGAAGCTTTTCCGACCATATAAAGTATATAATTGAGTATTTAGCTTAGGATTATTAATAAAACATGTATTTTATTTAAC</textarea>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<textarea rows='4' cols='120'>Spliced sequence showing SJ ('*'): CACACTAGACCCACCGATCGCGACTTTGTGTGTGTTAAAATGACGCGCATATCGTTTTGGGCTAATAATATTGTTCACGTTTTATAAAACAGATACACGGTTAAGCTCTGACCGGTCGTATTATTATCTGTGTGTTCTTGGTTGACACGTTTATTTCATTTATTCTCATTTAAAGCAACGATGTGAAGAACAGTGCGGGCGTTTCTCAACATAAGAATACATACAAACTTATGCGTGTTGGTGAGCAAATAAAATAGTTCAAAAAATTACTCGGCCACCGCCTTAACCGTACAAACGGAAATACATAAAG*CATAGAGGCGTAGTCGGTGAGTCAGCCACGCCAACAGTTAACAGTCGTCATTTTGCTGTCTGTCTAACACGCAATTTGTGTGTGTATATGTGAATAAACAGATACCGGTATGCAATCGCAAAGCTGGTCAACTGTAAAGATCG*TTTAGAATAAATACGGCCACAATTAGTATACATATTTATTCATCGAAACCAACCCACCACAAGAGCATAACCAATTATGTTAAAGTTTATCAGAGGAAAAGGGCAGCAGCCCAGTGCTGACAGACACCGCCTACAGAAGGACCTTTTTGCTTATCGTAAG*ACGGCACAGCATGGCTTTCCTCATAAGCCTTCGGCTCTTGCGTATGATCCAGTTTTGAAACTTATGGCAATAGGGACGCAAACAGGGGCTTTAAAAGTTTTCGGTCAACCCGGAGTTGAATTGTACGGTCAGCATACTTTGTTAAACAATTCAGCATCGGAGCTTAATGTACAATTACTTGAATGGGTGTATGGAACTGGTCGCATACTTTCGTTGACGGCAGCGAATCAATTAATTCTATGGGAGCCAGTTGGAGCAACGTTGCTGCCAATCAAAACACTACCGTTTGACGGCAAACTTAAAAAAGTTTCATCGCTGTGCTGTTCTCTCAGTAAGGATCTGCTATGGATTGGAACAGAAGGTGGAAACATCTATCAACTGGATTTACATACATTTACCATTAAGGAGCCTGTAATTTACCATGACGTTGTGCTAGAGCAGGTGCCACCAGCCTACAAGCTAAATCCTGGTGCAATTGAGTCAATCCGCCAACTTCCAAACTCCCCTAGCAAACTTCTAGTTGCATACAATCGCGGCCTTTGTGTTTTGTGGGATTTTGAAAGCGCATCTGTCCAGCGAGCATACATAGCCCCTGGACATGGACAGAGCGTTGGTCTTACAGTGAACTTCGAAGGATCTGAATTTACCTGGTACCACGCTGATGGTTCATACGCCACTTGGAGCATAGATAACCCAGAACCGCCGTCGAATGTTAATTATGTGCCTTATGGACCTGATCCATGCAAAAGCATAAATCGACTGTACAAAGGCAAGCGAAG*ATCCAACGATGTAATTGTTTTTTCCGGCGGCATGCCACGGTCAGCATATGGTGATCACAATTGTGTGTCCGTTCACGCCAGCGATGGACACAAAGTGTGTCTTGACTTTACGTCTAAAGTGATTGACTTTTTTGTGACCTTTGAAAATAATAGAGATGTCGCTGAAGTTCTTGTTGTACTACTTGAAGAGGAACTCTGCGCTTACGATCTTACTGACCCTAATATTTGTGCTATCAAAGCGCCATATCTTCACTCTGTCCATGCATCAGCTGTCACTTGCAATTACCTTGCTTCTGAAGTCGTACAGTCGGTATATGAAAGTATTTTAAGAGCTGGAGATGAACAAGACATTGACTATAGCAATATTAGCTGGCCTATCACTGGCGGTACTCTCCCGGATAACTTAGAAGAATCTGTAGAAGAGGACGCGACTAAGCTTTATGAGATTTTGTTAACTGGTCACGAAGATGGTTCTGTTAAATTTTGGGACTGCACTGGAGTGTTGCTTAAACCAATTTATAATTTTAAAACTAGCAGCATTTTTGGAAGTGAGTCAGACTTCCGAGATGACGCAGCTGCAGATATGAGTGCGGAACAAGTCGATGAAGGAGAACCGCCATTTCGGAAATCAGGACTTTTTGATCCTTATTCAGATGACCCTCGTTTAGCAGTGAAGAAAATAGCATTCTGCCCAAAAACCGGACAACTTATTGTTGGTGGCACAGCGGGCCAAATAGTTATAGCCGACTTCATAGACTTACCCGAAAAAGTGTCTTTAAAATACATTTCAATGAATTTGGTCAGCGATCGTGATGGATTTGTGTGGAAGGGTCACGATCAGTTAAACGTGCGATCGAACTTATTAGACGGAGAAGCAATTCCTACGACGGAACGTGGTGTAAATATATCGGGAGTACTGCAAGTTTTGCCGCCAGCCAGCATAACATGCATGGCACTCGAAGCAAGCTGGGGCCTAGTATCTGGTGGGACTGCGCACGGCTTAGTTCTCTTTGACTTCAAAAACTTTGTTCCAGTATTTCATCGCTGCACTTTAAACCCAAATGATCTTACTGGAGCAGGAGAGCAGCTGTCTCGTCGAAAGTCTTTTAAGAAATCATTGAGGGAGTCATTTAGAAAGCTTCGCAAGGGTCGATCGACCAGGACCAACCAGAGCAATCAAGTACCAACAACG*CTGGAAGCAAGACCCGTCGAGAGGCAAATAGAGGCTCGTTGTGCAGATGACGGGCTAGGATCCATGGTGCGATGTTTACTATTTGCCAAAACTTATGTTACTAATG*TCAACATAACGTCGCCAACTTTGTGGTCAGCAACAAATGCCAGTACAGTCTCGGTTTTCCTTCTGCATTTGCCACCAGCGCAGACCGCGGCAACTGCCGTCCCGTCGGCAAGTGGCAATGCACCACCACACATGCCCCGCCGAATTTCTGCGCAGCTTGCTAAAGAAATACAATTAAAACATCGTGCTCCTGTGGTGGGTATTTCTATTTTTGATCAGGCGGGTAGCCCTGTCGATCAGCTGAACGCCGGTGAAAACGGGAGTCCACCGCATCGTGTACTTATTGCTTCCGAGGAACAGTTCAAGGTGTTTTCACTTCCGCAACTAAAGCCGATTAACAAATATAAGCTTACCGCTAACGAAGGTGCTCGGATTCGCCGCATCCATTTTGGTTCGTTTAGTTGTCGCATATCCCCGGAAACACTGCAGAGTATGCACGGTTGTAGCCCAACTAAGTCCACGCGTTCACATGGCGATGGAGAGGCGGATCCTAATATCAGTGGAAGCTTGGCTGTAAGTCGTGGAGATGTATATAACGAAACAGCATTGATATGTTTAACGAATATGGGCGATATCATGGTTTTATCAGTACCTGAATTAAAAAGACAGCTGAATGCCGCAGCAGTGCGACGGGAAGACATTAA*TGGAGTTTCGTCACTTTGCTTTACAAACTCTGGAGAAGCACTGTATATGATGTCTTCTTCTGAACTGCAGCGTATTGCTTTAGCCACGTCCAGAGTCGTGCAACCCACTGGCGTTGTTCCAGTAGAACCATTAGAAAATGAAGAGTCTGTGTTGGAAGAAAATGATGCAGAGAATAATAAGGAAACCTACGCATGTGATGAAGTTGTGAATACATATGAAATTAAAAATCCATCAGGCATTTCAATATGCACAAGGCCTGCAGAGGAAAACGTTGGAAGAAATAGTGTTCAGCAAGTTAATGGAGTCAACATTTCAAATTCACCTAATCAAGCTAACGAGACTATCAGCAGCTCTATTGGCGATATTACCGTTGACTCGGTGCGCGACCATTTAAATATGACGACCACCACTTTGTGTTCTATTAATACAGAGGAAACCATTG*GTCGCCTATCTGTACTTAGCACGCAAACCAACAAAGCCAGTACTACCGTAAACATGAGTGAAATTCCAAATATTAATATTTCTAATTTAGAGGACTTGGAATCGAAAAG*AAATACGACGGAAACGAGTACTAGTTCTGTTGTAATTAAATCTATAATTACAAACATTTCTCATGAAAAAACGAACGGAGACAACAAAATAGGAACGCCAAAAACAGCGCCTGAAGAAAGCCAATTTTAACATTGACAGAAGCCGTAACCTACTAATTATTTTATACCTTTACAGAAATACGACGCAAACGAGTACTAGTCCTGATGTAATTTAATCTATAATTACAAACATGTCTCACTAAAAATCGAACGGGGACAGCAAAATAGGAACGCAAAAAAATTAAATTAAATTAAACGCGATTTACATACAAACAGAAATGACAGAATGATAATATAAAATATTCATTTTTTATTTGGCTTAAGCGATGTTGTTGTTCCAAAACCATATAATTATTTATGATTTTATGTAATGTTTTCATGTATTTTCGCGTAGTGACTTTATACCCTTTACCCGTATGATGAAACGAAACATGATATTTGTAGAAAGGTACAAAGGATAAAAGTTGACATTTAAGGCCGATTGTAAAATATTAAAACGCAACCTGGATAACATGCCGATTAATTACTTCTGTTCGTCCGTATAAACACAAATACCGATAAATTGGTCCGTCCCCATTTAATAAATTGTATATTAAAATGGTCTTTATTTTGGATAAATAATTCAAATTATATAGGAACATTTGGGTTTTAAGGGATAGAAACAGCGCTGGTGTCTTTCATTTCCGTTCGTATGACCGTTCAAAGATTATAAAAAAATTGTCGTGCCCACATTTTTTAAATTTCTTTTATTTTCGGTTAAAATATTGTGATACAAAAAAATGCCCTTACAAAATTAGGACAACAACAATTTGCTTTATTATGATGAACAAATTAGGAGTACAAAAAATAGCCTGATGCAAGGATGGGGAGAAAAAAGACCTTGAAGTTCAAAAGAATTTGTTGTATCAAATCCTCAATAAGAAAACTACTTAAAATTAAAAGCGTGTTTCAAATTAATTGTAGAATATTTCGTAAAATGAAATTGCAAACTTTTTTAGCACGAAATTACATGTATTTAATTCTTTGGCTAATACGTCAGTTAAAAAAGGCTACATTTTGATGATACTCACAAATGTAATTTTTGCATTCTACACAGGTGCAATTAAAAAATTATGGTTTTTATGAGTTTTAGAGGGAAGAAAATATGTCAATATGAGAAGTAAGTTAATGAGCCTTCAAGTACTTTAGTACATAATGTAATTTAATTTTAGAAATCCTACGGCTTTATATTCAATTGTAAAGTATTTTAAAATTCCTGACGATTTGTAGTAATTGTTTTGTTCTTAGCTTTTTTGCATTCTCCTCCGTTACTCCGTTACTCGTATAGTAAAAGGTATAGTATATACAATAGTATTGTAGATTCGTTGAGAAGTATGTTACAGGCAGAAGGAAGCTTTTCCGACCATATAAAGTATATAATTGAGTATTTAGCTTAGGATTATTAATAAAACATGTATTTTATTTAAC</textarea>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'tx FBtr0309226, exon 1: dnaseq=CTGATGAGATGTAATTCGGGCAATTTGTACACATTCGTTGCACACCATTTTGCATTGAT, rnaseq=ATCAATGCAAAATGGTGTGCAACGAATGTGTACAAATTGCCCGAATTACATCTCATCAG'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rand_tx = random.sample(t['l(2)gl'].transcript,1)[0] # sample a random transcript of l(2)gl\n",
    "display(f\"Selected tx: {rand_tx.feature_id}\") \n",
    "display_textarea(f\"DNA sequence of 1st exon: {rand_tx.exon[0].sequence}\") # get 'dna' sequence \n",
    "display_textarea(f\"Spliced sequence: {rand_tx.spliced_sequence}\")\n",
    "\n",
    "# SJs can be displayed by calling get_sequence(mode='spliced', show_exon_boundaries=True)\n",
    "display_textarea(f\"Spliced sequence showing SJ ('*'): {t.get_sequence(rand_tx, mode='spliced', show_exon_boundaries=True)}\")\n",
    "\n",
    "# sequences will be reverse complemented for - strand transcripts if mode='rna'\n",
    "# get some random minustrand tx\n",
    "tx=random.sample([tx for tx in t.transcripts if tx.strand=='-'], 1)[0]\n",
    "\n",
    "# show sequence of 1st exon in dna and rna mode (reverse complement as the gene is on the - strand)\n",
    "display(f\"tx {tx.feature_id}, exon 1: dnaseq={t.get_sequence(tx.exon[0], mode='dna')}, rnaseq={t.get_sequence(tx.exon[0], mode='rna')}\", SEP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6075b31a-bbd3-4b77-93e6-2fb8bd81a775",
   "metadata": {},
   "source": [
    "Note that the transcriptome implementation also supports (annotated) CDS annotations. Users can directly access the translated sequence and _in silico_ translate\n",
    "to amino acid sequences as demonstrated here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c0dd24b2-462e-46ac-8879-cb578d3dfdd2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<textarea rows='4' cols='120'>Translated sequence of FBtr0330655: ATGTTAAAGTTTATCAGAGGAAAAGGGCAGCAGCCCAGTGCTGACAGACACCGCCTACAGAAGGACCTTTTTGCTTATCGTAAGACGGCACAGCATGGCTTTCCTCATAAGCCTTCGGCTCTTGCGTATGATCCAGTTTTGAAACTTATGGCAATAGGGACGCAAACAGGGGCTTTAAAAGTTTTCGGTCAACCCGGAGTTGAATTGTACGGTCAGCATACTTTGTTAAACAATTCAGCATCGGAGCTTAATGTACAATTACTTGAATGGGTGTATGGAACTGGTCGCATACTTTCGTTGACGGCAGCGAATCAATTAATTCTATGGGAGCCAGTTGGAGCAACGTTGCTGCCAATCAAAACACTACCGTTTGACGGCAAACTTAAAAAAGTTTCATCGCTGTGCTGTTCTCTCAGTAAGGATCTGCTATGGATTGGAACAGAAGGTGGAAACATCTATCAACTGGATTTACATACATTTACCATTAAGGAGCCTGTAATTTACCATGACGTTGTGCTAGAGCAGGTGCCACCAGCCTACAAGCTAAATCCTGGTGCAATTGAGTCAATCCGCCAACTTCCAAACTCCCCTAGCAAACTTCTAGTTGCATACAATCGCGGCCTTTGTGTTTTGTGGGATTTTGAAAGCGCATCTGTCCAGCGAGCATACATAGCCCCTGGACATGGACAGAGCGTTGGTCTTACAGTGAACTTCGAAGGATCTGAATTTACCTGGTACCACGCTGATGGTTCATACGCCACTTGGAGCATAGATAACCCAGAACCGCCGTCGAATGTTAATTATGTGCCTTATGGACCTGATCCATGCAAAAGCATAAATCGACTGTACAAAGGCAAGCGAAGATCCAACGATGTAATTGTTTTTTCCGGCGGCATGCCACGGTCAGCATATGGTGATCACAATTGTGTGTCCGTTCACGCCAGCGATGGACACAAAGTGTGTCTTGACTTTACGTCTAAAGTGATTGACTTTTTTGTGACCTTTGAAAATAATAGAGATGTCGCTGAAGTTCTTGTTGTACTACTTGAAGAGGAACTCTGCGCTTACGATCTTACTGACCCTAATATTTGTGCTATCAAAGCGCCATATCTTCACTCTGTCCATGCATCAGCTGTCACTTGCAATTACCTTGCTTCTGAAGTCGTACAGTCGGTATATGAAAGTATTTTAAGAGCTGGAGATGAACAAGACATTGACTATAGCAATATTAGCTGGCCTATCACTGGCGGTACTCTCCCGGATAACTTAGAAGAATCTGTAGAAGAGGACGCGACTAAGCTTTATGAGATTTTGTTAACTGGTCACGAAGATGGTTCTGTTAAATTTTGGGACTGCACTGGAGTGTTGCTTAAACCAATTTATAATTTTAAAACTAGCAGCATTTTTGGAAGTGAGTCAGACTTCCGAGATGACGCAGCTGCAGATATGAGTGCGGAACAAGTCGATGAAGGAGAACCGCCATTTCGGAAATCAGGACTTTTTGATCCTTATTCAGATGACCCTCGTTTAGCAGTGAAGAAAATAGCATTCTGCCCAAAAACCGGACAACTTATTGTTGGTGGCACAGCGGGCCAAATAGTTATAGCCGACTTCATAGACTTACCCGAAAAAGTGTCTTTAAAATACATTTCAATGAATTTGGTCAGCGATCGTGATGGATTTGTGTGGAAGGGTCACGATCAGTTAAACGTGCGATCGAACTTATTAGACGGAGAAGCAATTCCTACGACGGAACGTGGTGTAAATATATCGGGAGTACTGCAAGTTTTGCCGCCAGCCAGCATAACATGCATGGCACTCGAAGCAAGCTGGGGCCTAGTATCTGGTGGGACTGCGCACGGCTTAGTTCTCTTTGACTTCAAAAACTTTGTTCCAGTATTTCATCGCTGCACTTTAAACCCAAATGATCTTACTGGAGCAGGAGAGCAGCTGTCTCGTCGAAAGTCTTTTAAGAAATCATTGAGGGAGTCATTTAGAAAGCTTCGCAAGGGTCGATCGACCAGGACCAACCAGAGCAATCAAGTACCAACAACGCTGGAAGCAAGACCCGTCGAGAGGCAAATAGAGGCTCGTTGTGCAGATGACGGGCTAGGATCCATGGTGCGATGTTTACTATTTGCCAAAACTTATGTTACTAATGTCAACATAACGTCGCCAACTTTGTGGTCAGCAACAAATGCCAGTACAGTCTCGGTTTTCCTTCTGCATTTGCCACCAGCGCAGACCGCGGCAACTGCCGTCCCGTCGGCAAGTGGCAATGCACCACCACACATGCCCCGCCGAATTTCTGCGCAGCTTGCTAAAGAAATACAATTAAAACATCGTGCTCCTGTGGTGGGTATTTCTATTTTTGATCAGGCGGGTAGCCCTGTCGATCAGCTGAACGCCGGTGAAAACGGGAGTCCACCGCATCGTGTACTTATTGCTTCCGAGGAACAGTTCAAGGTGTTTTCACTTCCGCAACTAAAGCCGATTAACAAATATAAGCTTACCGCTAACGAAGGTGCTCGGATTCGCCGCATCCATTTTGGTTCGTTTAGTTGTCGCATATCCCCGGAAACACTGCAGAGTATGCACGGTTGTAGCCCAACTAAGTCCACGCGTTCACATGGCGATGGAGAGGCGGATCCTAATATCAGTGGAAGCTTGGCTGTAAGTCGTGGAGATGTATATAACGAAACAGCATTGATATGTTTAACGAATATGGGCGATATCATGGTTTTATCAGTACCTGAATTAAAAAGACAGCTGAATGCCGCAGCAGTGCGACGGGAAGACATTAATGGAGTTTCGTCACTTTGCTTTACAAACTCTGGAGAAGCACTGTATATGATGTCTTCTTCTGAACTGCAGCGTATTGCTTTAGCCACGTCCAGAGTCGTGCAACCCACTGGCGTTGTTCCAGTAGAACCATTAGAAAATGAAGAGTCTGTGTTGGAAGAAAATGATGCAGAGAATAATAAGGAAACCTACGCATGTGATGAAGTTGTGAATACATATGAAATTAAAAATCCATCAGGCATTTCAATATGCACAAGGCCTGCAGAGGAAAACGTTGGAAGAAATAGTGTTCAGCAAGTTAATGGAGTCAACATTTCAAATTCACCTAATCAAGCTAACGAGACTATCAGCAGCTCTATTGGCGATATTACCGTTGACTCGGTGCGCGACCATTTAAATATGACGACCACCACTTTGTGTTCTATTAATACAGAGGAAACCATTGAAATACGACGGAAACGAGTACTAGTTCTGTTG</textarea>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<textarea rows='4' cols='120'>target            0 ------------------------------------------------------------\n",
       "                  0 ------------------------------------------------------------\n",
       "query             0 GCACACTAGACCCACCGATCGCGACTTTGTGTGTGTTAAAATGACGCGCATATCGTTTTG\n",
       "\n",
       "target            0 ------------------------------------------------------------\n",
       "                 60 ------------------------------------------------------------\n",
       "query            60 GGCTAATAATATTGTTCACGTTTTATAAAACAGATACACGGTTAAGCTCTGACCGGTCGT\n",
       "\n",
       "target            0 ------------------------------------------------------------\n",
       "                120 ------------------------------------------------------------\n",
       "query           120 ATTATTATCTGTGTGTTCTTGGTTGACACGTTTATTTCATTTATTCTCATTTAAAGCAAC\n",
       "\n",
       "target            0 ------------------------------------------------------------\n",
       "                180 ------------------------------------------------------------\n",
       "query           180 GATGTGAAGAACAGTGCGGGCGTTTCTCAACATAAGAATACATACAAACTTATGCGTGAT\n",
       "\n",
       "target            0 ------------------------------------------------------------\n",
       "                240 ------------------------------------------------------------\n",
       "query           240 GGTGAGCAAATAAAATAGTTCAAAAAATTACTCGGCCACCGCCTTAACCGTACAAACGGA\n",
       "\n",
       "target            0 ------------------------------------------------------------\n",
       "                300 ------------------------------------------------------------\n",
       "query           300 AATACATAAAGTTTAGAATAAATTCGGCCACAATTAATATACTTATTTATTCATCAAAAC\n",
       "\n",
       "target            0 ------------------------------------------------------------\n",
       "                360 ------------------------------------------------------------\n",
       "query           360 CAAACCACCACAAGAGCATAACCAATCATGTTAAAGTTTATCAGAGGAAAAGGGCAGCAG\n",
       "\n",
       "target            0 ------------------------------------------------------------\n",
       "                420 ------------------------------------------------------------\n",
       "query           420 CCCAGTGCTGACAGACAACGCTACATCGTAAGTTTAGAATAAATACGGCCACAATTAGTA\n",
       "\n",
       "target            0 ------------------------------------------------ATGTTAAAGTTT\n",
       "                480 ------------------------------------------------||||||||||||\n",
       "query           480 TACATATTTATTCATCGAAACCAACCCACCACAAGAGCATAACCAATTATGTTAAAGTTT\n",
       "\n",
       "target           12 ATCAGAGGAAAAGGGCAGCAGCCCAGTGCTGACAGACACCGCCTACAGAAGGACCTTTTT\n",
       "                540 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query           540 ATCAGAGGAAAAGGGCAGCAGCCCAGTGCTGACAGACACCGCCTACAGAAGGACCTTTTT\n",
       "\n",
       "target           72 GCTTATCGTAAGACGGCACAGCATGGCTTTCCTCATAAGCCTTCGGCTCTTGCGTATGAT\n",
       "                600 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query           600 GCTTATCGTAAGACGGCACAGCATGGCTTTCCTCATAAGCCTTCGGCTCTTGCGTATGAT\n",
       "\n",
       "target          132 CCAGTTTTGAAACTTATGGCAATAGGGACGCAAACAGGGGCTTTAAAAGTTTTCGGTCAA\n",
       "                660 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query           660 CCAGTTTTGAAACTTATGGCAATAGGGACGCAAACAGGGGCTTTAAAAGTTTTCGGTCAA\n",
       "\n",
       "target          192 CCCGGAGTTGAATTGTACGGTCAGCATACTTTGTTAAACAATTCAGCATCGGAGCTTAAT\n",
       "                720 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query           720 CCCGGAGTTGAATTGTACGGTCAGCATACTTTGTTAAACAATTCAGCATCGGAGCTTAAT\n",
       "\n",
       "target          252 GTACAATTACTTGAATGGGTGTATGGAACTGGTCGCATACTTTCGTTGACGGCAGCGAAT\n",
       "                780 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query           780 GTACAATTACTTGAATGGGTGTATGGAACTGGTCGCATACTTTCGTTGACGGCAGCGAAT\n",
       "\n",
       "target          312 CAATTAATTCTATGGGAGCCAGTTGGAGCAACGTTGCTGCCAATCAAAACACTACCGTTT\n",
       "                840 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query           840 CAATTAATTCTATGGGAGCCAGTTGGAGCAACGTTGCTGCCAATCAAAACACTACCGTTT\n",
       "\n",
       "target          372 GACGGCAAACTTAAAAAAGTTTCATCGCTGTGCTGTTCTCTCAGTAAGGATCTGCTATGG\n",
       "                900 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query           900 GACGGCAAACTTAAAAAAGTTTCATCGCTGTGCTGTTCTCTCAGTAAGGATCTGCTATGG\n",
       "\n",
       "target          432 ATTGGAACAGAAGGTGGAAACATCTATCAACTGGATTTACATACATTTACCATTAAGGAG\n",
       "                960 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query           960 ATTGGAACAGAAGGTGGAAACATCTATCAACTGGATTTACATACATTTACCATTAAGGAG\n",
       "\n",
       "target          492 CCTGTAATTTACCATGACGTTGTGCTAGAGCAGGTGCCACCAGCCTACAAGCTAAATCCT\n",
       "               1020 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1020 CCTGTAATTTACCATGACGTTGTGCTAGAGCAGGTGCCACCAGCCTACAAGCTAAATCCT\n",
       "\n",
       "target          552 GGTGCAATTGAGTCAATCCGCCAACTTCCAAACTCCCCTAGCAAACTTCTAGTTGCATAC\n",
       "               1080 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1080 GGTGCAATTGAGTCAATCCGCCAACTTCCAAACTCCCCTAGCAAACTTCTAGTTGCATAC\n",
       "\n",
       "target          612 AATCGCGGCCTTTGTGTTTTGTGGGATTTTGAAAGCGCATCTGTCCAGCGAGCATACATA\n",
       "               1140 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1140 AATCGCGGCCTTTGTGTTTTGTGGGATTTTGAAAGCGCATCTGTCCAGCGAGCATACATA\n",
       "\n",
       "target          672 GCCCCTGGACATGGACAGAGCGTTGGTCTTACAGTGAACTTCGAAGGATCTGAATTTACC\n",
       "               1200 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1200 GCCCCTGGACATGGACAGAGCGTTGGTCTTACAGTGAACTTCGAAGGATCTGAATTTACC\n",
       "\n",
       "target          732 TGGTACCACGCTGATGGTTCATACGCCACTTGGAGCATAGATAACCCAGAACCGCCGTCG\n",
       "               1260 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1260 TGGTACCACGCTGATGGTTCATACGCCACTTGGAGCATAGATAACCCAGAACCGCCGTCG\n",
       "\n",
       "target          792 AATGTTAATTATGTGCCTTATGGACCTGATCCATGCAAAAGCATAAATCGACTGTACAAA\n",
       "               1320 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1320 AATGTTAATTATGTGCCTTATGGACCTGATCCATGCAAAAGCATAAATCGACTGTACAAA\n",
       "\n",
       "target          852 GGCAAGCGAAGATCCAACGATGTAATTGTTTTTTCCGGCGGCATGCCACGGTCAGCATAT\n",
       "               1380 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1380 GGCAAGCGAAGATCCAACGATGTAATTGTTTTTTCCGGCGGCATGCCACGGTCAGCATAT\n",
       "\n",
       "target          912 GGTGATCACAATTGTGTGTCCGTTCACGCCAGCGATGGACACAAAGTGTGTCTTGACTTT\n",
       "               1440 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1440 GGTGATCACAATTGTGTGTCCGTTCACGCCAGCGATGGACACAAAGTGTGTCTTGACTTT\n",
       "\n",
       "target          972 ACGTCTAAAGTGATTGACTTTTTTGTGACCTTTGAAAATAATAGAGATGTCGCTGAAGTT\n",
       "               1500 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1500 ACGTCTAAAGTGATTGACTTTTTTGTGACCTTTGAAAATAATAGAGATGTCGCTGAAGTT\n",
       "\n",
       "target         1032 CTTGTTGTACTACTTGAAGAGGAACTCTGCGCTTACGATCTTACTGACCCTAATATTTGT\n",
       "               1560 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1560 CTTGTTGTACTACTTGAAGAGGAACTCTGCGCTTACGATCTTACTGACCCTAATATTTGT\n",
       "\n",
       "target         1092 GCTATCAAAGCGCCATATCTTCACTCTGTCCATGCATCAGCTGTCACTTGCAATTACCTT\n",
       "               1620 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1620 GCTATCAAAGCGCCATATCTTCACTCTGTCCATGCATCAGCTGTCACTTGCAATTACCTT\n",
       "\n",
       "target         1152 GCTTCTGAAGTCGTACAGTCGGTATATGAAAGTATTTTAAGAGCTGGAGATGAACAAGAC\n",
       "               1680 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1680 GCTTCTGAAGTCGTACAGTCGGTATATGAAAGTATTTTAAGAGCTGGAGATGAACAAGAC\n",
       "\n",
       "target         1212 ATTGACTATAGCAATATTAGCTGGCCTATCACTGGCGGTACTCTCCCGGATAACTTAGAA\n",
       "               1740 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1740 ATTGACTATAGCAATATTAGCTGGCCTATCACTGGCGGTACTCTCCCGGATAACTTAGAA\n",
       "\n",
       "target         1272 GAATCTGTAGAAGAGGACGCGACTAAGCTTTATGAGATTTTGTTAACTGGTCACGAAGAT\n",
       "               1800 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1800 GAATCTGTAGAAGAGGACGCGACTAAGCTTTATGAGATTTTGTTAACTGGTCACGAAGAT\n",
       "\n",
       "target         1332 GGTTCTGTTAAATTTTGGGACTGCACTGGAGTGTTGCTTAAACCAATTTATAATTTTAAA\n",
       "               1860 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1860 GGTTCTGTTAAATTTTGGGACTGCACTGGAGTGTTGCTTAAACCAATTTATAATTTTAAA\n",
       "\n",
       "target         1392 ACTAGCAGCATTTTTGGAAGTGAGTCAGACTTCCGAGATGACGCAGCTGCAGATATGAGT\n",
       "               1920 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1920 ACTAGCAGCATTTTTGGAAGTGAGTCAGACTTCCGAGATGACGCAGCTGCAGATATGAGT\n",
       "\n",
       "target         1452 GCGGAACAAGTCGATGAAGGAGAACCGCCATTTCGGAAATCAGGACTTTTTGATCCTTAT\n",
       "               1980 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          1980 GCGGAACAAGTCGATGAAGGAGAACCGCCATTTCGGAAATCAGGACTTTTTGATCCTTAT\n",
       "\n",
       "target         1512 TCAGATGACCCTCGTTTAGCAGTGAAGAAAATAGCATTCTGCCCAAAAACCGGACAACTT\n",
       "               2040 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2040 TCAGATGACCCTCGTTTAGCAGTGAAGAAAATAGCATTCTGCCCAAAAACCGGACAACTT\n",
       "\n",
       "target         1572 ATTGTTGGTGGCACAGCGGGCCAAATAGTTATAGCCGACTTCATAGACTTACCCGAAAAA\n",
       "               2100 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2100 ATTGTTGGTGGCACAGCGGGCCAAATAGTTATAGCCGACTTCATAGACTTACCCGAAAAA\n",
       "\n",
       "target         1632 GTGTCTTTAAAATACATTTCAATGAATTTGGTCAGCGATCGTGATGGATTTGTGTGGAAG\n",
       "               2160 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2160 GTGTCTTTAAAATACATTTCAATGAATTTGGTCAGCGATCGTGATGGATTTGTGTGGAAG\n",
       "\n",
       "target         1692 GGTCACGATCAGTTAAACGTGCGATCGAACTTATTAGACGGAGAAGCAATTCCTACGACG\n",
       "               2220 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2220 GGTCACGATCAGTTAAACGTGCGATCGAACTTATTAGACGGAGAAGCAATTCCTACGACG\n",
       "\n",
       "target         1752 GAACGTGGTGTAAATATATCGGGAGTACTGCAAGTTTTGCCGCCAGCCAGCATAACATGC\n",
       "               2280 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2280 GAACGTGGTGTAAATATATCGGGAGTACTGCAAGTTTTGCCGCCAGCCAGCATAACATGC\n",
       "\n",
       "target         1812 ATGGCACTCGAAGCAAGCTGGGGCCTAGTATCTGGTGGGACTGCGCACGGCTTAGTTCTC\n",
       "               2340 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2340 ATGGCACTCGAAGCAAGCTGGGGCCTAGTATCTGGTGGGACTGCGCACGGCTTAGTTCTC\n",
       "\n",
       "target         1872 TTTGACTTCAAAAACTTTGTTCCAGTATTTCATCGCTGCACTTTAAACCCAAATGATCTT\n",
       "               2400 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2400 TTTGACTTCAAAAACTTTGTTCCAGTATTTCATCGCTGCACTTTAAACCCAAATGATCTT\n",
       "\n",
       "target         1932 ACTGGAGCAGGAGAGCAGCTGTCTCGTCGAAAGTCTTTTAAGAAATCATTGAGGGAGTCA\n",
       "               2460 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2460 ACTGGAGCAGGAGAGCAGCTGTCTCGTCGAAAGTCTTTTAAGAAATCATTGAGGGAGTCA\n",
       "\n",
       "target         1992 TTTAGAAAGCTTCGCAAGGGTCGATCGACCAGGACCAACCAGAGCAATCAAGTACCAACA\n",
       "               2520 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2520 TTTAGAAAGCTTCGCAAGGGTCGATCGACCAGGACCAACCAGAGCAATCAAGTACCAACA\n",
       "\n",
       "target         2052 ACGCTGGAAGCAAGACCCGTCGAGAGGCAAATAGAGGCTCGTTGTGCAGATGACGGGCTA\n",
       "               2580 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2580 ACGCTGGAAGCAAGACCCGTCGAGAGGCAAATAGAGGCTCGTTGTGCAGATGACGGGCTA\n",
       "\n",
       "target         2112 GGATCCATGGTGCGATGTTTACTATTTGCCAAAACTTATGTTACTAATGTCAACATAACG\n",
       "               2640 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2640 GGATCCATGGTGCGATGTTTACTATTTGCCAAAACTTATGTTACTAATGTCAACATAACG\n",
       "\n",
       "target         2172 TCGCCAACTTTGTGGTCAGCAACAAATGCCAGTACAGTCTCGGTTTTCCTTCTGCATTTG\n",
       "               2700 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2700 TCGCCAACTTTGTGGTCAGCAACAAATGCCAGTACAGTCTCGGTTTTCCTTCTGCATTTG\n",
       "\n",
       "target         2232 CCACCAGCGCAGACCGCGGCAACTGCCGTCCCGTCGGCAAGTGGCAATGCACCACCACAC\n",
       "               2760 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2760 CCACCAGCGCAGACCGCGGCAACTGCCGTCCCGTCGGCAAGTGGCAATGCACCACCACAC\n",
       "\n",
       "target         2292 ATGCCCCGCCGAATTTCTGCGCAGCTTGCTAAAGAAATACAATTAAAACATCGTGCTCCT\n",
       "               2820 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2820 ATGCCCCGCCGAATTTCTGCGCAGCTTGCTAAAGAAATACAATTAAAACATCGTGCTCCT\n",
       "\n",
       "target         2352 GTGGTGGGTATTTCTATTTTTGATCAGGCGGGTAGCCCTGTCGATCAGCTGAACGCCGGT\n",
       "               2880 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2880 GTGGTGGGTATTTCTATTTTTGATCAGGCGGGTAGCCCTGTCGATCAGCTGAACGCCGGT\n",
       "\n",
       "target         2412 GAAAACGGGAGTCCACCGCATCGTGTACTTATTGCTTCCGAGGAACAGTTCAAGGTGTTT\n",
       "               2940 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          2940 GAAAACGGGAGTCCACCGCATCGTGTACTTATTGCTTCCGAGGAACAGTTCAAGGTGTTT\n",
       "\n",
       "target         2472 TCACTTCCGCAACTAAAGCCGATTAACAAATATAAGCTTACCGCTAACGAAGGTGCTCGG\n",
       "               3000 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          3000 TCACTTCCGCAACTAAAGCCGATTAACAAATATAAGCTTACCGCTAACGAAGGTGCTCGG\n",
       "\n",
       "target         2532 ATTCGCCGCATCCATTTTGGTTCGTTTAGTTGTCGCATATCCCCGGAAACACTGCAGAGT\n",
       "               3060 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          3060 ATTCGCCGCATCCATTTTGGTTCGTTTAGTTGTCGCATATCCCCGGAAACACTGCAGAGT\n",
       "\n",
       "target         2592 ATGCACGGTTGTAGCCCAACTAAGTCCACGCGTTCACATGGCGATGGAGAGGCGGATCCT\n",
       "               3120 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          3120 ATGCACGGTTGTAGCCCAACTAAGTCCACGCGTTCACATGGCGATGGAGAGGCGGATCCT\n",
       "\n",
       "target         2652 AATATCAGTGGAAGCTTGGCTGTAAGTCGTGGAGATGTATATAACGAAACAGCATTGATA\n",
       "               3180 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          3180 AATATCAGTGGAAGCTTGGCTGTAAGTCGTGGAGATGTATATAACGAAACAGCATTGATA\n",
       "\n",
       "target         2712 TGTTTAACGAATATGGGCGATATCATGGTTTTATCAGTACCTGAATTAAAAAGACAGCTG\n",
       "               3240 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          3240 TGTTTAACGAATATGGGCGATATCATGGTTTTATCAGTACCTGAATTAAAAAGACAGCTG\n",
       "\n",
       "target         2772 AATGCCGCAGCAGTGCGACGGGAAGACATTAATGGAGTTTCGTCACTTTGCTTTACAAAC\n",
       "               3300 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          3300 AATGCCGCAGCAGTGCGACGGGAAGACATTAATGGAGTTTCGTCACTTTGCTTTACAAAC\n",
       "\n",
       "target         2832 TCTGGAGAAGCACTGTATATGATGTCTTCTTCTGAACTGCAGCGTATTGCTTTAGCCACG\n",
       "               3360 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          3360 TCTGGAGAAGCACTGTATATGATGTCTTCTTCTGAACTGCAGCGTATTGCTTTAGCCACG\n",
       "\n",
       "target         2892 TCCAGAGTCGTGCAACCCACTGGCGTTGTTCCAGTAGAACCATTAGAAAATGAAGAGTCT\n",
       "               3420 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          3420 TCCAGAGTCGTGCAACCCACTGGCGTTGTTCCAGTAGAACCATTAGAAAATGAAGAGTCT\n",
       "\n",
       "target         2952 GTGTTGGAAGAAAATGATGCAGAGAATAATAAGGAAACCTACGCATGTGATGAAGTTGTG\n",
       "               3480 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          3480 GTGTTGGAAGAAAATGATGCAGAGAATAATAAGGAAACCTACGCATGTGATGAAGTTGTG\n",
       "\n",
       "target         3012 AATACATATGAAATTAAAAATCCATCAGGCATTTCAATATGCACAAGGCCTGCAGAGGAA\n",
       "               3540 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          3540 AATACATATGAAATTAAAAATCCATCAGGCATTTCAATATGCACAAGGCCTGCAGAGGAA\n",
       "\n",
       "target         3072 AACGTTGGAAGAAATAGTGTTCAGCAAGTTAATGGAGTCAACATTTCAAATTCACCTAAT\n",
       "               3600 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          3600 AACGTTGGAAGAAATAGTGTTCAGCAAGTTAATGGAGTCAACATTTCAAATTCACCTAAT\n",
       "\n",
       "target         3132 CAAGCTAACGAGACTATCAGCAGCTCTATTGGCGATATTACCGTTGACTCGGTGCGCGAC\n",
       "               3660 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          3660 CAAGCTAACGAGACTATCAGCAGCTCTATTGGCGATATTACCGTTGACTCGGTGCGCGAC\n",
       "\n",
       "target         3192 CATTTAAATATGACGACCACCACTTTGTGTTCTATTAATACAGAGGAAACCATTGAAATA\n",
       "               3720 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n",
       "query          3720 CATTTAAATATGACGACCACCACTTTGTGTTCTATTAATACAGAGGAAACCATTGAAATA\n",
       "\n",
       "target         3252 CGACGGAAACGAGTACTAGTTCTGTTG---------------------------------\n",
       "               3780 |||||||||||||||||||||||||||---------------------------------\n",
       "query          3780 CGACGGAAACGAGTACTAGTTCTGTTGTAATTAAATCTATAATTACAAACATTTCTCATG\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               3840 ------------------------------------------------------------\n",
       "query          3840 AAAAAACGAACGGAGACAACAAAATAGGAACGCCAAAAACAGCGCCTGAAGAAAGCCAAT\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               3900 ------------------------------------------------------------\n",
       "query          3900 TTTAACATTGACAGAAGCCGTAACCTACTAATTATTTTATACCTTTACAGAAATACGACG\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               3960 ------------------------------------------------------------\n",
       "query          3960 CAAACGAGTACTAGTCCTGATGTAATTTAATCTATAATTACAAACATGTCTCACTAAAAA\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4020 ------------------------------------------------------------\n",
       "query          4020 TCGAACGGGGACAGCAAAATAGGAACGCAAAAAAATTAAATTAAATTAAACGCGATTTAC\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4080 ------------------------------------------------------------\n",
       "query          4080 ATACAAACAGAAATGACAGAATGATAATATAAAATATTCATTTTTTATTTGGCTTAAGCG\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4140 ------------------------------------------------------------\n",
       "query          4140 ATGTTGTTGTTCCAAAACCATATAATTATTTATGATTTTATGTAATGTTTTCATGTATTT\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4200 ------------------------------------------------------------\n",
       "query          4200 TCGCGTAGTGACTTTATACCCTTTACCCGTATGATGAAACGAAACATGATATTTGTAGAA\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4260 ------------------------------------------------------------\n",
       "query          4260 AGGTACAAAGGATAAAAGTTGACATTTAAGGCCGATTGTAAAATATTAAAACGCAACCTG\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4320 ------------------------------------------------------------\n",
       "query          4320 GATAACATGCCGATTAATTACTTCTGTTCGTCCGTATAAACACAAATACCGATAAATTGG\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4380 ------------------------------------------------------------\n",
       "query          4380 TCCGTCCCCATTTAATAAATTGTATATTAAAATGGTCTTTATTTTGGATAAATAATTCAA\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4440 ------------------------------------------------------------\n",
       "query          4440 ATTATATAGGAACATTTGGGTTTTAAGGGATAGAAACAGCGCTGGTGTCTTTCATTTCCG\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4500 ------------------------------------------------------------\n",
       "query          4500 TTCGTATGACCGTTCAAAGATTATAAAAAAATTGTCGTGCCCACATTTTTTAAATTTCTT\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4560 ------------------------------------------------------------\n",
       "query          4560 TTATTTTCGGTTAAAATATTGTGATACAAAAAAATGCCCTTACAAAATTAGGACAACAAC\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4620 ------------------------------------------------------------\n",
       "query          4620 AATTTGCTTTATTATGATGAACAAATTAGGAGTACAAAAAATAGCCTGATGCAAGGATGG\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4680 ------------------------------------------------------------\n",
       "query          4680 GGAGAAAAAAGACCTTGAAGTTCAAAAGAATTTGTTGTATCAAATCCTCAATAAGAAAAC\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4740 ------------------------------------------------------------\n",
       "query          4740 TACTTAAAATTAAAAGCGTGTTTCAAATTAATTGTAGAATATTTCGTAAAATGAAATTGC\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4800 ------------------------------------------------------------\n",
       "query          4800 AAACTTTTTTAGCACGAAATTACATGTATTTAATTCTTTGGCTAATACGTCAGTTAAAAA\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4860 ------------------------------------------------------------\n",
       "query          4860 AGGCTACATTTTGATGATACTCACAAATGTAATTTTTGCATTCTACACAGGTGCAATTAA\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4920 ------------------------------------------------------------\n",
       "query          4920 AAAATTATGGTTTTTATGAGTTTTAGAGGGAAGAAAATATGTCAATATGAGAAGTAAGTT\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               4980 ------------------------------------------------------------\n",
       "query          4980 AATGAGCCTTCAAGTACTTTAGTACATAATGTAATTTAATTTTAGAAATCCTACGGCTTT\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               5040 ------------------------------------------------------------\n",
       "query          5040 ATATTCAATTGTAAAGTATTTTAAAATTCCTGACGATTTGTAGTAATTGTTTTGTTCTTA\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               5100 ------------------------------------------------------------\n",
       "query          5100 GCTTTTTTGCATTCTCCTCCGTTACTCCGTTACTCGTATAGTAAAAGGTATAGTATATAC\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               5160 ------------------------------------------------------------\n",
       "query          5160 AATAGTATTGTAGATTCGTTGAGAAGTATGTTACAGGCAGAAGGAAGCTTTTCCGACCAT\n",
       "\n",
       "target         3279 ------------------------------------------------------------\n",
       "               5220 ------------------------------------------------------------\n",
       "query          5220 ATAAAGTATATAATTGAGTATTTAGCTTAGGATTATTAATAAAACATGTATTTTATTTAA\n",
       "\n",
       "target         3279 - 3279\n",
       "               5280 - 5281\n",
       "query          5280 C 5281\n",
       "</textarea>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<textarea rows='4' cols='120'>Amino acid sequence of FBtr0078164: MLKFIRGKGQQPSADRHRLQKDLFAYRKTAQHGFPHKPSALAYDPVLKLMAIGTQTGALKVFGQPGVELYGQHTLLNNSASELNVQLLEWVYGTGRILSLTAANQLILWEPVGATLLPIKTLPFDGKLKKVSSLCCSLSKDLLWIGTEGGNIYQLDLHTFTIKEPVIYHDVVLEQVPPAYKLNPGAIESIRQLPNSPSKLLVAYNRGLCVLWDFESASVQRAYIAPGHGQSVGLTVNFEGSEFTWYHADGSYATWSIDNPEPPSNVNYVPYGPDPCKSINRLYKGKRRSNDVIVFSGGMPRSAYGDHNCVSVHASDGHKVCLDFTSKVIDFFVTFENNRDVAEVLVVLLEEELCAYDLTDPNICAIKAPYLHSVHASAVTCNYLASEVVQSVYESILRAGDEQDIDYSNISWPITGGTLPDNLEESVEEDATKLYEILLTGHEDGSVKFWDCTGVLLKPIYNFKTSSIFGSESDFRDDAAADMSAEQVDEGEPPFRKSGLFDPYSDDPRLAVKKIAFCPKTGQLIVGGTAGQIVIADFIDLPEKVSLKYISMNLVSDRDGFVWKGHDQLNVRSNLLDGEAIPTTERGVNISGVLQVLPPASITCMALEASWGLVSGGTAHGLVLFDFKNFVPVFHRCTLNPNDLTGAGEQLSRRKSFKKSLRESFRKLRKGRSTRTNQSNQVPTTLEARPVERQIEARCADDGLGSMVRCLLFAKTYVTNVNITSPTLWSATNASTVSVFLLHLPPAQTAATAVPSASGNAPPHMPRRISAQLAKEIQLKHRAPVVGISIFDQAGSPVDQLNAGENGSPPHRVLIASEEQFKVFSLPQLKPINKYKLTANEGARIRRIHFGSFSCRISPETLQSMHGCSPTKSTRSHGDGEADPNISGSLAVSRGDVYNETALICLTNMGDIMVLSVPELKRQLNAAAVRREDINGVSSLCFTNSGEALYMMSSSELQRIALATSRVVQPTGVVPVEPLENEESVLEENDAENNKETYACDEVVNTYEIKNPSGISICTRPAEENVGRNSVQQVNGVNISNSPNQANETISSSIGDITVDSVRDHLNMTTTTLCSINTEETIEIRRKRVLVLL</textarea>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# The transcriptome implementation also supports (annotated) CDS annotations. \n",
    "# Users can directly access the translated sequence:\n",
    "display_textarea(f\"Translated sequence of FBtr0330655: {t['FBtr0330655'].translated_sequence}\")\n",
    "\n",
    "# align spliced (including 5'/3' UTRs) and translated sequence using biopython\n",
    "display_textarea(next(PairwiseAligner(mode='global', open_gap_score=-3).align(\n",
    "        t['FBtr0330655'].translated_sequence, \n",
    "        t['FBtr0330655'].spliced_sequence)))\n",
    "    \n",
    "# translate with biotite\n",
    "display_textarea(f\"Amino acid sequence of FBtr0078164: {seq.NucleotideSequence(t['FBtr0330655'].translated_sequence).translate(complete=True)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ae15bb8-8a80-422a-89ce-9508cfa8f3b7",
   "metadata": {},
   "source": [
    "Here are some more complex access examples that demonstrate exact kmer search via standard python methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d705390c-ef35-43a9-b751-a1e5f5ce8ce5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"Genes containing ATGC kmer in one of their (spliced) tx: {'asRNA:CR43609', 'lncRNA:CR45339', 'CG11374', 'lncRNA:CR45340', 'dbr', 'asRNA:CR44987', 'lncRNA:CR46254', 'l(2)gl', 'CR11023', 'galectin', 'Cda5', 'net', 'Zir', 'Ir21a'}\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"Genes containing ATGC kmer in a long (>1kb) intron: {'Cda5', 'net', 'galectin', 'l(2)gl'}\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# list all gene names where the kmer ATGC is found in one of the (spliced) transcripts\n",
    "display(f\"Genes containing ATGC kmer in one of their (spliced) tx: { { tx.parent.gene_name for tx in t.transcripts if 'ATGC' in tx.spliced_sequence }}\", SEP)\n",
    "\n",
    "# list all gene names where the kmer ATGC is found in one of its transcripts introns considering only long (>1kb) introns\n",
    "display(f\"Genes containing ATGC kmer in a long (>1kb) intron: { { tx.parent.gene_name for tx in t.transcripts for intron in tx.intron if len(intron)>1000 and 'ATGC' in intron.sequence}}\", SEP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02dcf90d-e276-4310-9a5b-581e8fb73e48",
   "metadata": {},
   "source": [
    "## Querying\n",
    "Transcriptomes can be queried by standard python methods (e.g., list comprehension) or via location (range) queries that are efficiently implemented by intervaltrees.\n",
    "Here we briefly demonstrate both approaches:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "4cfb1374-87bf-4dda-9a02-8d7d59f6fe9b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Cda5', 'l(2)gl', 'net'}"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# find all minus strand genes with transcripts that have >=2 exons and at least one long (>1kb) intron via list comprehension\n",
    "{tx.parent.gene_name for tx in t.transcripts if tx.strand=='-' and len(tx.exon)>=2 and any([len(i)>1000 for i in tx.intron])}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e57ead4-b04d-4cc7-9e73-f568d9467bba",
   "metadata": {},
   "source": [
    "Efficient range queries are implemented by prefiltering the data using per-chromosome \n",
    "[intervaltrees](https://github.com/chaimleib/intervaltree) that are constructed automatically at transcriptome instantiation time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8d068450-f717-4701-9860-c1b2a807deb8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"Genes where an exon overlaps with 2L:10000-20000: {'l(2)gl'}\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"Genes of all transcripts that overlap with gene 'Ir21a': {'Ir21a', 'asRNA:CR43609'}\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"Genes near 'l(2)gl': ['CR11023', 'l(2)gl', 'Ir21a', 'asRNA:CR43609']\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"Genes contained in a +/- 10kb window around 'l(2)gl': ['CR11023', 'l(2)gl', 'asRNA:CR43609']\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"Genes at 2L:20000: ['FBgn0002121']\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# range queries are implemented by intervaltrees:\n",
    "display(f\"Genes where an exon overlaps with 2L:10000-20000: { {ex.parent.parent.gene_name for ex in t.query('2L:10000-20000', 'exon')}}\", SEP) \n",
    "\n",
    "# Query all transcripts that overlap with the 'Ir21a' gene and create a set with their gene names.\n",
    "display(f\"Genes of all transcripts that overlap with gene 'Ir21a': {set(tx.parent.gene_name for tx in t.query(t['Ir21a'], 'transcript'))}\", SEP) \n",
    "\n",
    "# report a coordinate-sorted list of genes in a 10kb window around 'l(2)gl':\n",
    "display(f\"Genes near 'l(2)gl': { [g.gene_name for g in t.query(gi(t['l(2)gl'].chromosome, t['l(2)gl'].start-3000, t['l(2)gl'].end+3000), 'gene')]}\", SEP)\n",
    "\n",
    "# envelop queries (i.e., queries that report only intervals strictly within the query interval) are also supported:\n",
    "display(f\"Genes contained in a +/- 10kb window around 'l(2)gl': { [g.gene_name for g in t.query(gi(t['l(2)gl'].chromosome, t['l(2)gl'].start-3000, t['l(2)gl'].end+3000), 'gene', envelop=True)]}\", SEP)\n",
    "\n",
    "# Direct access to the interval trees is possible via the chr2itree attribute. Here, we query intervals\n",
    "# at position 2L:20000. The respective gene object can be accessed via the data slot of the returned intervaltree Interval objects.\n",
    "display(f\"Genes at 2L:20000: {[iv.data.feature_id for iv in t.chr2itree['2L'].at(20000)]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39c245ce-2479-4381-90e3-4263983993b3",
   "metadata": {},
   "source": [
    "## Summary of access modes\n",
    "This cell demonstrates different ways to access (meta) data associated with transcriptome features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "ae55b3a6-4ee8-4911-b484-39c785ee2bc4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'transcriptome': Transcriptome with 14 genes and 41 tx (+seq),\n",
       " 'feature_id': 'FBtr0078103_exon_2',\n",
       " 'feature_type': 'exon',\n",
       " 'parent': transcript@2L:94739-102086,\n",
       " 'subfeature_types': (),\n",
       " 'source': 'FlyBase',\n",
       " 'gene_type': None,\n",
       " 'gff_feature_type': 'exon'}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "{'my_great_annotation': 24.9}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'FBtr0078103'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'ATCTCCCCGA'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[gene@2L:94739-102086,\n",
       " transcript@2L:94739-102086,\n",
       " CDS@2L:95132-95301,\n",
       " exon@2L:95132-95301,\n",
       " CDS@2L:95132-95301,\n",
       " exon@2L:95132-95301]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "{'Cda5', 'l(2)gl', 'net'}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ex = t['FBtr0078103'].exon[2]         # get the 3rd exon of a transcript with given id\n",
    "t[ex]['my_great_annotation'] = 24.9   # add arbitrary annotations\n",
    "display(vars(ex))                     # show data that is directly associated with this exon\n",
    "display(t[ex])                        # show custom annotations\n",
    "display(ex.rnk)                       # implicitly modeled annotations: rnk is index in exon list\n",
    "display(ex.parent.feature_id)         # access via feature hierarchy (gene->transcript->exon)\n",
    "display(ex.sequence[:10])             # derived annotations: access gene-associated sequence data and slice exon coordinates\n",
    "display(t.query(ex))                  # query for any genomic interval via intervaltrees\n",
    "display({tx.parent.gene_name for tx in t.transcripts if tx.strand=='-' \n",
    "         and len(tx.exon)>=2 and any([len(i)>1000 for i in tx.intron])}, SEP) # query via python list comprehension\n",
    "# and much, much more..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bf0f52e-1c1f-44c3-af27-e2c2d504718f",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Iterators"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4160ea04-a2e5-4ba8-8f3e-49578ecc6969",
   "metadata": {
    "tags": []
   },
   "source": [
    "*rnalib* implements several *genomic iterators* for efficient iteration over (parts of) genomics datasets while keeping track of the [genomic interval](#Genomic-intervals) each yielded item refers to. Most implemented iterators are based on respective [pysam](https://pysam.readthedocs.io/en/latest/api.html) classes but extend them with different filtering/access methods and add sanity checks to avoid common errors when integrating multiple genomics datasets. They also support chunked I/O where feasible and not supported by the underlying (pysam) implementation. The primary design goal of *genomic iterators* is to support the synchronization/integration of different genomics datasets. \n",
    "\n",
    "*Genomic iterators* can iterate whole datasets or specific sub-region (e.g., a certain chromosome or genomic region). Yielded data items are sorted by respective [genomic intervals](#Genomic-intervals) where the order of chromosomes is defined by a [RefDict](#Genomic-intervals) that is derived directly from the iterated dataset (e.g., by querying it's tabix index). When multiple iterators are integrated (e.g., by an AnnotationIterator, see below), their *RefDict*'s are checked for compatibility.\n",
    "\n",
    "Sometimes, genomics datasets from different providers used differing chromosome ids to refer to the same chromosomes which hinders their integration. Location iterators support *chromosome aliasing* as a quick solution for this: users can provide aliasing functions (e.g., for adding/removing 'chr' prefixes from chromosome ids) to dynamically change chromosome ids while keeping the underlying data files untouched.\n",
    "\n",
    "*Genomic iterators* inherit from the *LocationIterator* super class and provide a uniform access interface:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "173f9183-a0b4-4324-80bb-84b8a2040672",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<a id=\"help_icon0.17627878268889163\"><img id=\"help_icon\" style=\"cursor: help; float:left; width: 20px\" src=\"help_icon.png\"/>&nbsp;&nbsp;&nbsp;<b>LocationIterator docstring</b></a>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "    document.getElementById(\"help_icon0.17627878268889163\").addEventListener(\"click\", function() {\n",
       "        var win = window.open(\"\", \"\", \"toolbar=no, location=no, directories=no, status=no, menubar=no, scrollbars=yes, resizable=yes, width=1000, height=600\");\n",
       "        win.document.body.innerHTML = \"<h1>LocationIterator docstring</h1><p style='font-family:Courier New; font-size: 14px;'><br/>    Superclass for genomic iterators (mostly based on pysam) for efficient, indexed iteration over genomic<br/>    datasets. Most rnalib iterables inherit from this superclass.<br/><br/>    A LocationIterator iterates over a genomic dataset and yields `Item`&#x27;s, that are tuples of genomic intervals and<br/>    associated data. The data can be any object (e.g., a string, a dict, a list, etc.). Location iterators can be<br/>    restricted to a genomic region.<br/><br/>    LocationIterators keep track of their current genomic location (self.location) as well as statistics about<br/>    the iterated and yielded items (stats()). They can be consumed and converted to lists (self.to_list()),<br/>    pandas dataframes (self.to_dataframe()), interval trees (self.to_intervaltrees()) or bed files (self.to_bed()).<br/><br/>    LocationIterators can be grouped or tiled, see the group() and tile() methods.<br/><br/>    LocationIterators are iterable, implement the context manager protocol (&#x27;with&#x27; statement) and provide a<br/>    standard interface for region-type filtering (by passing chromosome sets or genomic regions of<br/>    interest). Where feasible and not supported by the underlying (pysam) implementation, they implement chunked<br/>    I/O where for efficient iteration (e.g., FastaIterator).<br/><br/>    The maximum number of iterated items can be queried via the max_items() method which tries to guesstimate this<br/>    number (or an upper bound) from the underlying index data structures.<br/><br/>    Examples<br/>    --------<br/>    &gt;&gt;&gt; # Creates a LocationIterator via factory method and consume its items, splitting into locs and data lists.<br/>    &gt;&gt;&gt; locs, data = zip(*it(&#x27;myfile.bed&#x27;))<br/><br/>    Attributes<br/>    ----------<br/>    _stats : Counter<br/>        A counter that collects statistics (e.g., yielded items per chromosome) for QC and reporting.<br/>    location : GI<br/>        The current genomic location of this iterator<br/><br/><br/>    Parameters<br/>    ----------<br/>    file : str or pysam file object<br/>        A file path or an open file object. In the latter case, the object will not be closed<br/>    region : GI or str<br/>        genomic region to iterate<br/>    file_format : str<br/>        optional, will be determined from filename if omitted<br/>    per_position : bool<br/>        optional, if True, the iterator will yield per-position items (e.g., FASTA files)<br/>    fun_alias : Callable<br/>        optional, if set, the iterator will use this function for aliasing chromosome names<br/>    refdict : RefDict<br/>        optional, if set, the iterator will use the passed reference dict instead of reading it from the file<br/><br/>    Notes<br/>    -----<br/>    TODOs:<br/><br/>    * strand specific iteration<br/>    * url streaming<br/>    * add is_exhausted flag?<br/>    </p>\";\n",
       "        win.document.title = \"LocationIterator docstring\";\n",
       "        win.document.body.style.background = \"#dcf6fa\";\n",
       "        });\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rna.display_help(rna.LocationIterator)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a130a49-555a-47ff-a5e7-68c6c60733fb",
   "metadata": {
    "tags": []
   },
   "source": [
    "Implemented iterators include:\n",
    "* [MemoryIterator](#MemoryIterator): iterates over genomic interval sets stored in various default python data structures (list, dict, set).\n",
    "* [FastaIterator](#FastaIterator): iterates FASTA files per position or genomic window (e.g., Tiling/sliding windows).\n",
    "* [FastqIterator](#FastqIterator): iterates FASTQ entries (Note that this is not a *LocationIterator*)\n",
    "* [GFF3Iterator](#GFF3Iterator): iterates GTF/GFF3 files and parses annotation data. Used for [transcriptome](#transcriptome) building.\n",
    "* [BigWigIterator](#BigWigIterator): iterates BigWig files; based on the pyBigWig library\n",
    "* [ReadIterator](#ReadIterator): iterates SAM/BAM files and implements various filter options. \n",
    "* [FastPileupIterator](#FastPileupIterator): a faster alternative to pysam's pileup method\n",
    "* [VcfIterator](#AnnotationIterator): iterates VCF files and implements additional access and filter methods\n",
    "* [AnnotationIterator](#AnnotationIterator): annotates items from one LocationIterator with items from (multiple) others\n",
    "\n",
    "and many more (see API documentation).\n",
    "\n",
    "Instantiation of an iterator can be done explicitly or via the `it(...)` factory method (recommended).\n",
    "The following sections showcase some of the implemented iterators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "0fca2eae-a424-417a-a00e-03785d02bb31",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "    <!-- gh-pages -->\n",
       "    <link rel='stylesheet' href='https://unpkg.com/freezeframe@3.0.10/build/css/freezeframe_styles.min.css'>\n",
       "    <script type='text/javascript' src='https://unpkg.com/freezeframe@3.0.10/build/js/freezeframe.pkgd.min.js'></script>\n",
       "\n",
       "    <script type='text/javascript'> \n",
       "      third = new freezeframe('.my_class_3').capture().setup();\n",
       "      $(function() {\n",
       "        $('.start').click(function(e) {\n",
       "          e.preventDefault();\n",
       "          third.trigger();\n",
       "        });\n",
       "      })\n",
       "    </script>\n",
       "    \n",
       "    <img class=\"my_class_3 freezeframe-responsive\" src=\"../docs/_static/screencasts/iterator_demo.gif\" />\n",
       "    <button class=\"start\">restart</button>\n",
       "    \n",
       "    <script>setTimeout( function() { third.trigger(); }, 1000);</script>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rna.display_animated_gif(\"../docs/_static/screencasts/iterator_demo.gif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10dc22c5-02ca-4b17-8fa0-90aa1c19d69b",
   "metadata": {},
   "source": [
    "## MemoryIterator\n",
    "Iterates over genomic interval sets stored in various default python data structures (list, dicts)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "341d8bbf-5adc-49a3-8659-916c55cbac8d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Item(location=chr1:10-100, data=2),\n",
       " Item(location=chr1:100-1000, data=0),\n",
       " Item(location=chr1:100-1000, data=1),\n",
       " Item(location=chr2:12-64, data=3)]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[Item(location=chr1:10-100, data='b'),\n",
       " Item(location=chr1:100-1000, data='a2'),\n",
       " Item(location=chr2:12-64, data='c')]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[Item(location=chr1:10-100, data='b'),\n",
       " Item(location=chr1:100-1000, data='a1'),\n",
       " Item(location=chr1:100-1000, data='a2'),\n",
       " Item(location=chr2:12-64, data='c')]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Here we iterate over a list of GIs. GIs will be sorted and their index (in the original input list) will be reported in the data section. \n",
    "# The it(...) factory method returns a MemoryIterator as we passed a list of GIs. We then consume this iterator and convert into a list.\n",
    "display(rna.it([gi('chr1:100-1000'), \n",
    "                gi('chr1:100-1000'), # duplicate interval that will be reported\n",
    "                gi('chr1:10-100'), \n",
    "                gi('chr2:12-64')]).to_list(), SEP )\n",
    "# Here we iterate over a dict and now the dict values will be reported in the data section.\n",
    "display(rna.it({gi('chr1:100-1000'): 'a1', \n",
    "                gi('chr1:100-1000'): 'a2', # This will replace 'a1' in the dict!\n",
    "                gi('chr1:10-100'):   'b', \n",
    "                gi('chr2:12-64'):    'c'}).to_list(), SEP )\n",
    "# We can also pass a (key:GI) dict which allows for iterating duplicate intervals with different IDs.\n",
    "# The MemoryIterator automatically detects whether the keys or the values of a dict are GIs.\n",
    "display(rna.it({'a1': gi('chr1:100-1000'), \n",
    "                'a2': gi('chr1:100-1000'), # a1 and a2 will both be yielded\n",
    "                'b':  gi('chr1:10-100'), \n",
    "                'c':  gi('chr2:12-64')}).to_list() )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "3d50e4d6-fd01-49d6-a725-49fb2857eaf6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<a id=\"help_icon0.024362536134039048\"><img id=\"help_icon\" style=\"cursor: help; float:left; width: 20px\" src=\"help_icon.png\"/>&nbsp;&nbsp;&nbsp;<b>MemoryIterator docstring</b></a>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "    document.getElementById(\"help_icon0.024362536134039048\").addEventListener(\"click\", function() {\n",
       "        var win = window.open(\"\", \"\", \"toolbar=no, location=no, directories=no, status=no, menubar=no, scrollbars=yes, resizable=yes, width=1000, height=600\");\n",
       "        win.document.body.innerHTML = \"<h1>MemoryIterator docstring</h1><p style='font-family:Courier New; font-size: 14px;'><br/>    A location iterator that iterates over intervals retrieved from one of the following datastructures:<br/><br/>    * {str:gi} dict: yields (gi,str); note that the passed strings must be unique<br/>    * {gi:any} dict: yields (gi,any)<br/>    * iterable of gis:  yields (gi, index in input iterable). If a single GI object is passed, its contained<br/>        positions are iterated<br/><br/>    Regions will be sorted according to a refdict that is created from the data automatically.<br/>    For associative datastructures (dicts), the key or value that is not the GI will be yielded as data item.<br/>    For others, the index in the sorted iterable will be yielded as data item.<br/><br/>    Notes<br/>    -----<br/>    * reported stats:<br/>        iterated_items, chromosome: (int, str)<br/>            Number of iterated items<br/>        yielded_items, chromosome: (int, str)<br/>            Number of yielded items (remaining are filtered, e.g., due to region constraints)<br/><br/>    </p>\";\n",
       "        win.document.title = \"MemoryIterator docstring\";\n",
       "        win.document.body.style.background = \"#dcf6fa\";\n",
       "        });\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rna.display_help(rna.MemoryIterator)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "262da942-1ef3-45da-ad9a-3da72463ad8e",
   "metadata": {
    "tags": []
   },
   "source": [
    "## FastaIterator\n",
    "Iterates over a FASTA file but enables also tiling/sliding windows and padding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "2d64ed77-f33f-4360-89c2-b80470284243",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"[(chr7:1-5, 'TTGTG'), (chr7:3-7, 'GTGCC'), (chr7:5-9, 'GCCAT'), (chr7:7-11, 'CATTA'), (chr7:9-13, 'TTACA')]...[(chr7:989-993, 'TCATC'), (chr7:991-995, 'ATCAT'), (chr7:993-997, 'CATTC'), (chr7:995-999, 'TTCAA'), (chr7:997-1001, 'CAAAA')]\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# get 5mer sliding windows (step size 2) with padding from a GRCh38 chr7 region around ACTB; \n",
    "# show first and last 5 results.\n",
    "# item.data returns the data item returned by this iterator (here: a string containing the kmer)\n",
    "# item.location (not used in this example) is a GI describing the genomic location of this sequence\n",
    "with rna.it(rna.get_resource('ACTB+SOX2_genome'), # a sliced version of the reference genome.\n",
    "            region='chr7:1-1001', # iterated region\n",
    "            width=5, # report 5-mers\n",
    "            step=2,  # step size: 2\n",
    "            padding=True) as it: # if padding is True, the yielded windows at the region edges will be padded by 'N' values\n",
    "    kmers = [(loc,kmer) for loc,kmer in it] # consume all items from the iterator\n",
    "display(f\"{kmers[:5]}...{kmers[-5:]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "294b5bcf-2579-41ca-aba0-233c3ed19aba",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<a id=\"help_icon0.4038800888284234\"><img id=\"help_icon\" style=\"cursor: help; float:left; width: 20px\" src=\"help_icon.png\"/>&nbsp;&nbsp;&nbsp;<b>FastaIterator docstring</b></a>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "    document.getElementById(\"help_icon0.4038800888284234\").addEventListener(\"click\", function() {\n",
       "        var win = window.open(\"\", \"\", \"toolbar=no, location=no, directories=no, status=no, menubar=no, scrollbars=yes, resizable=yes, width=1000, height=600\");\n",
       "        win.document.body.innerHTML = \"<h1>FastaIterator docstring</h1><p style='font-family:Courier New; font-size: 14px;'><br/>    Iterates over a FASTA file yielding sequence strings and keeps track of the covered genomic location.<br/><br/>    The yielded sequence of length &lt;width&gt; will be padded with &lt;fillvalue&gt; characters if padding is True.<br/>    This generator will yield every step_size window, for a tiling window approach set step_size = width.<br/>    FASTA files will be automatically indexed if no existing index file is found and will be read in chunked mode.<br/><br/>    Parameters<br/>    ----------<br/>    fasta_file : str or pysam FastaFile<br/>        A file path or an open FastaFile object. In the latter case, the object will not be closed<br/>    region : gi or str<br/>        optional, if set, only features overlapping with this region will be iterated<br/>    width : int<br/>        sequence window size<br/>    step : int<br/>        increment for each yield<br/>    file_format : str<br/>        optional, will be determined from filename if omitted<br/>    chunk_size : int<br/>        optional, chunk size for reading fasta files. default: 1024<br/>    fill_value : str<br/>        optional, fill value for padding. default: &#x27;N&#x27;<br/>    padding : bool<br/>        optional, if True, the sequence will be padded with fill_value characters. default: False<br/>    fun_alias : Callable<br/>        optional, if set, the iterator will use this function for aliasing chromosome names<br/><br/>    Yields<br/>    ------<br/>    sequence: str<br/>        The extracted sequence in including sequence context.<br/>        The core sequence w/o context can be accessed via seq[context_size:context_size+width].<br/><br/>    Notes<br/>    -----<br/>    * reported stats:<br/>        iterated_items, chromosome: (int, str)<br/>            Number of iterated and yielded items<br/>    </p>\";\n",
       "        win.document.title = \"FastaIterator docstring\";\n",
       "        win.document.body.style.background = \"#dcf6fa\";\n",
       "        });\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rna.display_help(rna.FastaIterator)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6d1b71d-897f-4b2a-a2b5-0aa61eca7243",
   "metadata": {},
   "source": [
    "## FastqIterator\n",
    "Iterates over a FASTQ file and returns named tuples (name seq qual) containing read names, sequence \n",
    "and quality strings. Note that this is no Location iterator as unaligned reads have no (known) genomic location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c924ba95-275d-415a-8f4c-91ac318ca8d1",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'@M01558:94:000000000-B4PB4:1:1101:16988:1392 1:N:0:20, @M01558:94:000000000-B4PB4:1:1101:16988:1392 2:N:0:20'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'len r1: 251, len r2: 250'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'@M01558:94:000000000-B4PB4:1:1101:16905:1394 1:N:0:20, @M01558:94:000000000-B4PB4:1:1101:16905:1394 2:N:0:20'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'len r1: 251, len r2: 251'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'@M01558:94:000000000-B4PB4:1:1101:16922:1403 1:N:0:20, @M01558:94:000000000-B4PB4:1:1101:16922:1403 2:N:0:20'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'len r1: 251, len r2: 250'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'@M01558:94:000000000-B4PB4:1:1101:15249:1404 1:N:0:20, @M01558:94:000000000-B4PB4:1:1101:15249:1404 2:N:0:20'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'len r1: 251, len r2: 250'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'@M01558:94:000000000-B4PB4:1:1101:15251:1423 1:N:0:20, @M01558:94:000000000-B4PB4:1:1101:15251:1423 2:N:0:20'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'len r1: 251, len r2: 250'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# iterate over PE reads in two FASTQ files and display their read names and sequence lengths.\n",
    "# We use pythons zip() method to access tuples of read mates\n",
    "for r1,r2 in zip(rna.it(rna.get_resource('small_PE_fastq1')), \n",
    "                 rna.it(rna.get_resource('small_PE_fastq2'))):\n",
    "    display(f\"{r1.name}, {r2.name}\")\n",
    "    display(f\"len r1: {len(r1.seq)}, len r2: {len(r2.seq)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "f9443a95-f5f8-40cf-962b-462da1b68655",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<a id=\"help_icon0.47687518987410604\"><img id=\"help_icon\" style=\"cursor: help; float:left; width: 20px\" src=\"help_icon.png\"/>&nbsp;&nbsp;&nbsp;<b>FastqIterator docstring</b></a>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "    document.getElementById(\"help_icon0.47687518987410604\").addEventListener(\"click\", function() {\n",
       "        var win = window.open(\"\", \"\", \"toolbar=no, location=no, directories=no, status=no, menubar=no, scrollbars=yes, resizable=yes, width=1000, height=600\");\n",
       "        win.document.body.innerHTML = \"<h1>FastqIterator docstring</h1><p style='font-family:Courier New; font-size: 14px;'><br/>    Iterates a fastq file and yields FastqRead objects (named tuples containing sequence name, sequence and quality<br/>    string (omitting line 3 in the FASTQ file)).<br/><br/>    Parameters<br/>    ----------<br/>    fastq_file : str<br/>        The FASTQ file name<br/><br/>    Notes<br/>    -----<br/>    Note that this is not a location iterator.<br/><br/>    reported stats:<br/><br/>    * yielded_items: int<br/>        Number of yielded items<br/>    </p>\";\n",
       "        win.document.title = \"FastqIterator docstring\";\n",
       "        win.document.body.style.background = \"#dcf6fa\";\n",
       "        });\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rna.display_help(rna.FastqIterator)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5404c692-c861-4919-8413-c781ac347856",
   "metadata": {},
   "source": [
    "## VCFIterator\n",
    "Iterates over VCF files and yields VcfRecords that wrap (and partially parse) pysam VcfProxy objects.\n",
    "Info fields and additional format columns are parsed and accessible via dot-notation: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "e214b56d-a71f-4d13-bb87-bb519d15a2cc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "##fileformat=VCFv4.2\n",
      "##contig=<ID=1,length=150>\n",
      "##contig=<ID=2,length=88>\n",
      "##contig=<ID=3,length=12>\n",
      "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n",
      "1\t1\t.\tA\tC\t0\t.\tDP=3;TCF=1.0\n",
      "1\t10\t.\tN\tT\t0\t.\tDP=3;TCF=1.0\n",
      "1\t20\t.\tN\tA\t0\t.\tDP=3;TCF=1.0\n",
      "2\t10\t.\tN\tT\t0\t.\tDP=3;TCF=1.0\n"
     ]
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A C 3 1.0\n",
      "N T 3 1.0\n",
      "N A 3 1.0\n",
      "N T 3 1.0\n"
     ]
    }
   ],
   "source": [
    "rna.print_small_file(rna.get_resource('test_snps_vcf'))\n",
    "display(SEP)\n",
    "with rna.VcfIterator(rna.get_resource('test_snps_vcf')) as it:\n",
    "    for loc, v in it:\n",
    "        print(v.ref, v.alt, v.DP, v.TCF)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7291aac5-912b-4f33-8284-f70f807a64fa",
   "metadata": {},
   "source": [
    "VCFIterators also provide convenience methods to access and filter contained samples and enables filtering based on called \n",
    "genotype, number of calls across samples or zygosity as shown in the [AdvancedUsage notebook](AdvancedUsage.ipynb).\n",
    "In the example below, we demonstrate how to access the names of all samples in a multi-sample VCF and how to access geneotypes by sample name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "8fd4d258-7327-44f0-9a6c-0cb2dac074e4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"Contained samples: 159, e.g., ['DGRP-738', 'DGRP-859', 'DGRP-59']\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2L:4851-4851 G C ./.\n",
      "2L:4860-4860 G A ./.\n",
      "2L:4891-4891 A T 0/0\n",
      "2L:4892-4892 A C 0/0\n",
      "2L:4893-4893 G A 0/0\n"
     ]
    }
   ],
   "source": [
    "#rna.print_small_file(rna.get_resource('dmel_multisample_vcf'))\n",
    "with rna.it(rna.get_resource('dmel_multisample_vcf')) as it:\n",
    "    # show included samples\n",
    "    display(f'Contained samples: {len(it.allsamples)}, e.g., {it.allsamples[:3]}', SEP)\n",
    "    # show some data for the first 5 variants\n",
    "    for loc,vr in itertools.islice(it, 5):\n",
    "        print(loc, vr.ref, vr.alt, vr.GT['DGRP-738'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "4d6e892f-6f7f-4ce0-b03c-f5b5b6678717",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<a id=\"help_icon0.3725770545677406\"><img id=\"help_icon\" style=\"cursor: help; float:left; width: 20px\" src=\"help_icon.png\"/>&nbsp;&nbsp;&nbsp;<b>VcfIterator docstring</b></a>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "    document.getElementById(\"help_icon0.3725770545677406\").addEventListener(\"click\", function() {\n",
       "        var win = window.open(\"\", \"\", \"toolbar=no, location=no, directories=no, status=no, menubar=no, scrollbars=yes, resizable=yes, width=1000, height=600\");\n",
       "        win.document.body.innerHTML = \"<h1>VcfIterator docstring</h1><p style='font-family:Courier New; font-size: 14px;'><br/>    Iterates a VCF file and yields 1-based coordinates and VcfRecord objects (wrapping pysam VcfProxy object)<br/>    @see https://pysam.readthedocs.io/en/latest/api.html#pysam.asVCF<br/><br/>    Besides coordinate-related filtering, users cam also pass a list of sample ids to be considered (samples).<br/>    The iterator will by default report only records with at least 1 call in one of the configured samples.<br/>    To force reporting of all records, set filter_nocalls=False.<br/><br/>    Parameters<br/>    ----------<br/>    vcf_file : str or pysam.VariantFile<br/>        A file path or an open VariantFile object. In the latter case, the object will not be closed<br/>    region : gi or str<br/>        optional, if set, only features overlapping with this region will be iterated<br/>    fun_alias : Callable<br/>        optional, if set, the iterator will use this function for aliasing chromosome names<br/>    samples : list<br/>        optional, if set, only records with calls in these samples will be reported<br/>    filter_nocalls : bool<br/>        if True (default), the iterator will report only records with at least 1 call in one of the configured<br/>        samples.<br/><br/>    Attributes<br/>    ----------<br/>    header:<br/>        pysam VariantFile header<br/>    allsamples: list<br/>        list of all contained samples<br/>    shownsampleindices:<br/>        indices of all configured samples<br/><br/>    Notes<br/>    -----<br/>    * reported stats:<br/>        iterated_items, chromosome: (int, str)<br/>            Number of iterated items<br/>        filtered_nocalls, chromosome: (int, str)<br/>            Number of filtered no-calls<br/>        yielded_items, chromosome: (int, str)<br/>            Number of yielded items (remaining are filtered, e.g., due to region constraints)<br/><br/>    </p>\";\n",
       "        win.document.title = \"VcfIterator docstring\";\n",
       "        win.document.body.style.background = \"#dcf6fa\";\n",
       "        });\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rna.display_help(rna.VcfIterator)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10b2536d-19fe-4f87-80f8-949b2163441e",
   "metadata": {
    "tags": []
   },
   "source": [
    "## GFF3Iterator\n",
    "Iterates over GFF3 files and yields dicts containing parsed annotation fields. This iterator\n",
    "is used for building transcriptome features as discussed below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "fa18afb5-5874-46a4-9fdd-80ac25276c92",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ul>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<li>Item(location=chr3:180989762-181836880 (+), data={'ID': 'ENSG00000242808.12', 'gene_id': 'ENSG00000242808.12', 'gene_type': 'lncRNA', 'gene_name': 'SOX2-OT', 'level': '2', 'hgnc_id': 'HGNC:20209', 'tag': 'ncRNA_host,overlapping_locus', 'havana_gene': 'OTTHUMG00000158186.13', 'feature_type': 'gene', 'source': 'HAVANA', 'score': None, 'phase': None})</li>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<li>Item(location=chr3:180989770-181791029 (+), data={'ID': 'ENST00000626948.3', 'Parent': 'ENSG00000242808.12', 'gene_id': 'ENSG00000242808.12', 'transcript_id': 'ENST00000626948.3', 'gene_type': 'lncRNA', 'gene_name': 'SOX2-OT', 'transcript_type': 'lncRNA', 'transcript_name': 'SOX2-OT-240', 'level': '2', 'transcript_support_level': '5', 'hgnc_id': 'HGNC:20209', 'tag': 'nested_454_RNA_Seq_supported,basic,Ensembl_canonical,TAGENE', 'havana_gene': 'OTTHUMG00000158186.13', 'havana_transcript': 'OTTHUMT00000481527.1', 'feature_type': 'transcript', 'source': 'HAVANA', 'score': None, 'phase': None})</li>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "</ul>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Counter({'exon': 108,\n",
       "         'transcript': 89,\n",
       "         'CDS': 60,\n",
       "         'five_prime_UTR': 33,\n",
       "         'three_prime_UTR': 20,\n",
       "         'start_codon': 17,\n",
       "         'stop_codon': 13,\n",
       "         'gene': 5})"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Counter({'lncRNA': 64,\n",
       "         'protein_coding': 16,\n",
       "         'retained_intron': 4,\n",
       "         'nonsense_mediated_decay': 2,\n",
       "         'protein_coding_CDS_not_defined': 2,\n",
       "         'processed_pseudogene': 1})"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Show the data section first two data items from a GFF3 file\n",
    "rna.display_list(rna.it(rna.get_resource('gencode_gff')).to_list()[:2])\n",
    "\n",
    "# iterate whole GFF3 file and collect stats of contained feature_type annotations\n",
    "display(Counter([d.get('feature_type') for _,d in rna.it(rna.get_resource('gencode_gff'))]))\n",
    "\n",
    "# iterate transcripts and collect stats of contained transcript_type annotations\n",
    "display(Counter([d.get('transcript_type', 'NA') for _,d in \\\n",
    "                 rna.it(rna.get_resource('gencode_gff'), feature_types=('transcript',))]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "8df9e445-ea3d-4fb4-9361-ea30b3089971",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<a id=\"help_icon0.8953289892384203\"><img id=\"help_icon\" style=\"cursor: help; float:left; width: 20px\" src=\"help_icon.png\"/>&nbsp;&nbsp;&nbsp;<b>GFF3Iterator docstring</b></a>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "    document.getElementById(\"help_icon0.8953289892384203\").addEventListener(\"click\", function() {\n",
       "        var win = window.open(\"\", \"\", \"toolbar=no, location=no, directories=no, status=no, menubar=no, scrollbars=yes, resizable=yes, width=1000, height=600\");\n",
       "        win.document.body.innerHTML = \"<h1>GFF3Iterator docstring</h1><p style='font-family:Courier New; font-size: 14px;'><br/>    Iterates a GTF/GFF3 file and yields 1-based coordinates and dicts containing key/value pairs parsed from<br/>    the respective info sections. The feature_type, source, score and phase fields from the GFF/GTF entries are<br/>    copied to this dict (NOTE: attribute fields with the same name will be overloaded).<br/>    @see https://pysam.readthedocs.io/en/latest/api.html#pysam.asGTF<br/><br/>    This iterator is used to build a transcriptome object from a GFF3/GTF file.<br/><br/>    Parameters<br/>    ----------<br/>    gtf_file : str or pysam.TabixFile<br/>        A file path or an open TabixFile object. In the latter case, the object will not be closed<br/>    region : gi or str<br/>        optional, if set, only features overlapping with this region will be iterated<br/>    fun_alias : Callable<br/>        optional, if set, the iterator will use this function for aliasing chromosome names<br/><br/>    Notes<br/>    -----<br/>    * reported stats:<br/>        yielded_items, chromosome: (int, str)<br/>            Number of iterated/yielded items<br/>    </p>\";\n",
       "        win.document.title = \"GFF3Iterator docstring\";\n",
       "        win.document.body.style.background = \"#dcf6fa\";\n",
       "        });\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rna.display_help(rna.GFF3Iterator)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2e8a5ca-8192-4ba5-bf60-6d7e9680509b",
   "metadata": {},
   "source": [
    "## BigWigIterator\n",
    "Wraps PyBigWig to iterates BigWig files (either per position or per interval)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "08dca159-f7ad-46b4-8db2-f450008a1022",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'1:1-1: 0.10000000149011612'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'1:2-2: 0.20000000298023224'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'1:3-3: 0.30000001192092896'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'1:101-150: 1.399999976158142'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'1:151-151: 1.5'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'10:201-300: 2.0'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# iterate a bigwig file and report intervals\n",
    "for loc, val in rna.it(rna.get_resource(\"test_bigwig\")):\n",
    "    display(f\"{loc}: {val}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "01ced1cd-5e24-4001-82a1-17b2e7d9a71a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "63f288c657a74c3c90f4a5cfee4e731e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Chromosome chr1'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr1:1-1: 0.10000000149011612'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr1:2-2: 0.20000000298023224'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr1:3-3: 0.30000001192092896'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr1:4-4: None'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr1:5-5: None'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Chromosome chr10'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr10:1-1: None'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr10:2-2: None'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr10:3-3: None'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr10:4-4: None'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr10:5-5: None'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# use itertools islice to show the first 5 values per chrom when iterating a bigwig file per genomic position (per_position=True)\n",
    "# note that per-position iteration is typically slow.\n",
    "with rna.it(rna.get_resource(\"test_bigwig\"), per_position=True, fun_alias=rna.toggle_chr) as it:\n",
    "    for chrom in tqdm(it.refdict):\n",
    "        display(f\"Chromosome {chrom}\")\n",
    "        it.set_region(gi(chrom, 1, 10000))\n",
    "        for loc, val in itertools.islice(it, 0, 5):\n",
    "            display(f\"{loc}: {val}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a093b8ba-40e3-4e73-81d4-49c978cc3ec3",
   "metadata": {},
   "source": [
    "You can also iterate BigBed files with a BigBedIterator that wraps a pyBigWig iterator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "de068dc6-2e42-4cab-8ca6-8734d75afbb4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"BigBed file with header {'version': 4, 'nLevels': 0, 'nBasesCovered': 24051, 'minVal': 1, 'maxVal': 2, 'sumData': 24268, 'sumSquared': 24702} and chromosomes ['chr1', 'chr11', 'chr15', 'chr2', 'chr4', 'chr7']\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr1:10009335-10009640 (-): 130.0'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr1:10014009-10014289 (-): 136.0'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr1:10014375-10024307 (-): 630.0'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr1:10024990-10028299 (-): 715.0'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'chr1:10028584-10029133 (+): 212.0'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'...'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# iterate a bigbed file and report intervals; note that chroms in this file are sorted 1,11,15,2,4,7\n",
    "with rna.it(rna.get_resource(\"test_bigbed\")) as it:\n",
    "    display(f'BigBed file with header {it.header()} and chromosomes {it.refdict.chromosomes()}')\n",
    "    for loc, val in itertools.islice(it, 5):\n",
    "        display(f\"{loc}: {val.score}\")\n",
    "    display('...')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "8248e442-f486-4a2a-abb5-ee994ccd3f95",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<a id=\"help_icon0.1697000864989857\"><img id=\"help_icon\" style=\"cursor: help; float:left; width: 20px\" src=\"help_icon.png\"/>&nbsp;&nbsp;&nbsp;<b>BigWigIterator docstring</b></a>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "    document.getElementById(\"help_icon0.1697000864989857\").addEventListener(\"click\", function() {\n",
       "        var win = window.open(\"\", \"\", \"toolbar=no, location=no, directories=no, status=no, menubar=no, scrollbars=yes, resizable=yes, width=1000, height=600\");\n",
       "        win.document.body.innerHTML = \"<h1>BigWigIterator docstring</h1><p style='font-family:Courier New; font-size: 14px;'>Iterates over a pyBigWig object.<br/>    If per_position is True, the iterator will yield per-position values, otherwise it will yield interval values.<br/><br/>    Parameters<br/>    ----------<br/>    file : str, PathLike or pyBigWig<br/>        The pyBigWig file to iterate over<br/>    region : gi or str<br/>        optional, if set, only features overlapping with this region will be iterated<br/>    fun_alias : Callable<br/>        optional, if set, the iterator will use this function for aliasing chromosome names<br/>    strand : str<br/>        optional, if set, all yielded intervals will have this strand assigned<br/>    per_position : bool<br/>        optional, if True, the iterator will yield per-position values, otherwise it will yield interval values<br/><br/><br/>    Notes<br/>    -----<br/>    Note that you can also access remote files, but this is slow and not recommended unless files are small.<br/>    Example:<br/>    &gt;&gt;&gt; # iterate a remote bigwig file and report intervals. This file is ~100M<br/>    &gt;&gt;&gt; import itertools<br/>    &gt;&gt;&gt; uri = &#x27;https://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeMapability/wgEncodeCrgMapabilityAlign100mer.bigWig&#x27; # noqa<br/>    &gt;&gt;&gt; for loc, val in itertools.islice(rna.it(uri),0,10):<br/>    &gt;&gt;&gt;     display(f&quot;{loc}: {val}&quot;)<br/>    </p>\";\n",
       "        win.document.title = \"BigWigIterator docstring\";\n",
       "        win.document.body.style.background = \"#dcf6fa\";\n",
       "        });\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rna.display_help(rna.BigWigIterator)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2987a903-200a-4ae4-b536-0c4671bb1757",
   "metadata": {
    "tags": []
   },
   "source": [
    "## ReadIterator\n",
    "ReadIterators iterates over indexed SAM/BAM files and implements several filter options (e.g., filtering by flag, mapping quality or tag values). They keeps track of how many reads were filtered for what reason in the `stats` Counter (e.g.: n_fil_mq: number of reads filtered due to too low mapping quaility).\n",
    "ReadIterators can also yield mismatches wrt. the reference if MD tags are available (MD tags can be added, e.g., by [samtools calmd](http://www.htslib.org/doc/samtools-calmd.html) if not added by the mapper itself). \n",
    "*Rnalib* also provides a `PairedReadIterator` for paired-end datasets, see the [AdvancedUsage notebook](AdvancedUsage.ipynb) for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "15b6f14e-cf57-4f9e-b100-4ed82a0f9944",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'all': Counter({('iterated_items', '1'): 31678,\n",
       "          ('yielded_items', '1'): 31678}),\n",
       " 'def': Counter({('iterated_items', '1'): 31678,\n",
       "          ('yielded_items', '1'): 21932,\n",
       "          ('n_fil_flag', '1'): 9746}),\n",
       " 'mq20': Counter({('iterated_items', '1'): 31678,\n",
       "          ('yielded_items', '1'): 21626,\n",
       "          ('n_fil_flag', '1'): 9746,\n",
       "          ('n_fil_mq', '1'): 306}),\n",
       " 'tag': Counter({('iterated_items', '1'): 31678,\n",
       "          ('yielded_items', '1'): 14544,\n",
       "          ('n_fil_flag', '1'): 9746,\n",
       "          ('n_fil_tag', '1'): 7388})}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# count reads per chromosome in a BAM using different filters\n",
    "stats={x:Counter() for x in ['all', 'def', 'mq20', 'tag']}\n",
    "bam_file = rna.get_resource('small_example_bam')\n",
    "with rna.it(bam_file, flag_filter=0) as it: # all: no filtering\n",
    "    it.to_list() # ignore results, we just want to compare the created stats\n",
    "    stats['all'].update(it.stats)\n",
    "with rna.it(bam_file) as it: # def: default flag filter (as in IGV)\n",
    "    it.to_list()\n",
    "    stats['def'].update(it.stats)\n",
    "with rna.it(bam_file, min_mapping_quality=20) as it: # mq20: default flag filter, mapping quality >= 20\n",
    "    it.to_list()\n",
    "    stats['mq20'].update(it.stats)\n",
    "with rna.it(bam_file, tag_filters=[rna.TagFilter('MD', ['100'], inverse=True)]) as it: # tag: default flag filter, MD tag must be '100'\n",
    "    it.to_list()\n",
    "    stats['tag'].update(it.stats)\n",
    "display(stats)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40faf7bb-d589-45c5-8373-14605e3b99d5",
   "metadata": {},
   "source": [
    "By passing 'include_unmapped=True', you can also iterate all unmapped reads. Note that the yielded location is empty for unmapped reads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d6f8f9e0-eb19-41bf-b5fb-a27323984133",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Item(location=None:0-2147483647, data=<pysam.libcalignedsegment.AlignedSegment object at 0x36a15c520>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Item(location=None:0-2147483647, data=<pysam.libcalignedsegment.AlignedSegment object at 0x347a233a0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Item(location=None:0-2147483647, data=<pysam.libcalignedsegment.AlignedSegment object at 0x36a15c640>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Item(location=None:0-2147483647, data=<pysam.libcalignedsegment.AlignedSegment object at 0x347a233a0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Item(location=None:0-2147483647, data=<pysam.libcalignedsegment.AlignedSegment object at 0x36a15c520>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Item(location=None:0-2147483647, data=<pysam.libcalignedsegment.AlignedSegment object at 0x347a233a0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Item(location=None:0-2147483647, data=<pysam.libcalignedsegment.AlignedSegment object at 0x36a15c4c0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Item(location=None:0-2147483647, data=<pysam.libcalignedsegment.AlignedSegment object at 0x347a233a0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Item(location=None:0-2147483647, data=<pysam.libcalignedsegment.AlignedSegment object at 0x36a15c640>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Item(location=None:0-2147483647, data=<pysam.libcalignedsegment.AlignedSegment object at 0x347a233a0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Counter({'unmapped_reads': 10})"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Here we iterate over a small file containing unmapped reads from the pybedtools test resources.\n",
    "# as pybedtools fdoes not include the respective .bai indices, we first copy the testfile to a temp dir, index it and then iterate the reads\n",
    "with tempfile.TemporaryDirectory() as tmp:\n",
    "    bam = f\"{tmp}/unmapped_reads.bam\"\n",
    "    import shutil\n",
    "    shutil.copyfile(rna.get_resource(\"pybedtools::issue_121.bam\"), bam)\n",
    "    pysam.index(bam)\n",
    "    with rna.it(bam, include_unmapped=True) as it:\n",
    "        for r in it:\n",
    "            display(r)\n",
    "        display(it.stats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d652a6f8-b9a3-4cc2-98ff-67183f3899fd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<a id=\"help_icon0.2987327029441057\"><img id=\"help_icon\" style=\"cursor: help; float:left; width: 20px\" src=\"help_icon.png\"/>&nbsp;&nbsp;&nbsp;<b>ReadIterator docstring</b></a>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "    document.getElementById(\"help_icon0.2987327029441057\").addEventListener(\"click\", function() {\n",
       "        var win = window.open(\"\", \"\", \"toolbar=no, location=no, directories=no, status=no, menubar=no, scrollbars=yes, resizable=yes, width=1000, height=600\");\n",
       "        win.document.body.innerHTML = \"<h1>ReadIterator docstring</h1><p style='font-family:Courier New; font-size: 14px;'>Iterates over a BAM alignment.<br/><br/>    Parameters<br/>    ----------<br/>    bam_file : str<br/>        BAM file name<br/>    region : gi or str<br/>        optional, if set, only features overlapping with this region will be iterated<br/>    file_format : str<br/>        optional, if set, the iterator will assume this file format (e.g., &#x27;bam&#x27; or &#x27;sam&#x27;)<br/>    min_mapping_quality : int<br/>        Minimum mapping quality. If set, reads with lower mapping quality will be filtered.<br/>    flag_filter : int<br/>        Bitwise flag filter. If set, reads with any of the specified flags will be filtered.<br/>    tag_filters : list<br/>        List of TagFilter objects. If set, reads with any of the specified tags will be filtered.<br/>    max_span : int<br/>        Maximum alignment span. If set, reads with a larger alignment span will be filtered.<br/>    report_mismatches : bool<br/>        If set, this iterator will additionally yield (read, mismatches) tuples where &#x27;mismatches&#x27; is a list of<br/>        (read_position, genomic_position, ref_allele, alt_allele) tuples that describe differences wrt. the<br/>        aligned reference sequence. This options requires MD BAM tags to be present.<br/>    min_base_quality : int<br/>        Useful only in combination with report_mismatches; filters mismatches based on minimum per-base quality<br/>    fun_alias : Callable<br/>        optional, if set, the iterator will use this function for aliasing chromosome names<br/><br/>    Notes<br/>    -----<br/>    * Reported stats<br/>        iterated_items, chromosome: (int, str)<br/>            Number of iterated items<br/>        yielded_items, chromosome: (int, str)<br/>            Number of yielded items<br/>        n_fil_flag, chromosome:<br/>            Number of reads filtered due to a FLAG filter<br/>        n_fil_mq, chromosome:<br/>            Number of reads filtered due to low mapping quality<br/>        n_fil_tag chromosome:<br/>            Number of reads filtered due to a TAG filter<br/>        n_fil_max_span, chromosome:<br/>            Number of reads filtered due to exceeding a configured maximum alignment span<br/><br/>    </p>\";\n",
       "        win.document.title = \"ReadIterator docstring\";\n",
       "        win.document.body.style.background = \"#dcf6fa\";\n",
       "        });\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rna.display_help(rna.ReadIterator)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fdf7876-3d0f-475b-a5bf-2ea3bea35a09",
   "metadata": {
    "tags": []
   },
   "source": [
    "## FastPileupIterator\n",
    "*Rnalib* implements a `FastPileupIterator` for fast and lightweight pileup-style access to alignment columns in SAM/BAM files.\n",
    "This implementation is ~10-20X faster than [pysam](https://pysam.readthedocs.io/en/latest/)'s pileup method, see the [AdvancedUsage notebook](AdvancedUsage.ipynb).\n",
    "Note: to instantiate a `FastPileupIterator` via the `it()` factory method, you must pass the  style='pileup' parameter as shown."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "97ff5122-fd16-455d-b3dc-4b7835533bc3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"1:22413312-22413312: Counter({'C': 1396})\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"1:22413313-22413313: Counter({'G': 1379})\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"1:22413314-22413314: Counter({'T': 1374, 'G': 1, 'C': 1})\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"1:22413315-22413315: Counter({'G': 1370})\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"1:22413316-22413316: Counter({'A': 1361, 'C': 1, 'G': 1})\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"1:22413317-22413317: Counter({'C': 1356})\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for loc, dat in rna.it(rna.get_resource('small_example_bam'),  # iterate over a BAM file\n",
    "                       style='pileup',                         # pileup style\n",
    "                       region=gi('1:22413312-22413317')):      # region to pileup\n",
    "    display(f\"{loc}: {dat}\")  # print location and allele counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "2b900962-7993-4844-a43a-fde9bb7fb378",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<a id=\"help_icon0.6950095289382656\"><img id=\"help_icon\" style=\"cursor: help; float:left; width: 20px\" src=\"help_icon.png\"/>&nbsp;&nbsp;&nbsp;<b>FastPileupIterator docstring</b></a>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "    document.getElementById(\"help_icon0.6950095289382656\").addEventListener(\"click\", function() {\n",
       "        var win = window.open(\"\", \"\", \"toolbar=no, location=no, directories=no, status=no, menubar=no, scrollbars=yes, resizable=yes, width=1000, height=600\");\n",
       "        win.document.body.innerHTML = \"<h1>FastPileupIterator docstring</h1><p style='font-family:Courier New; font-size: 14px;'><br/>    Fast pileup iterator that yields a complete pileup (w/o insertions) over a set of genomic positions. This is<br/>    more lightweight and considerably faster than `pysam&#x27;s pileup()` but lacks some features (such as<br/>    `ignore_overlaps` or `ignore_orphans`).<br/>    By default, it basically reports what is seen in the default IGV view (using, e.g., the same read flag filter).<br/><br/>    This iterator uses a ReadIterator to iterate the BAM file and then builds a pileup from the reads that overlap<br/>    the reported positions. The pileup is a dict with genomic positions as keys and a Counter object as values. The<br/>    Counter object contains the counts of each base at the respective position.<br/><br/>    Parameters<br/>    ----------<br/>    bam_file : str<br/>        BAM file name<br/>    chromosome : str<br/>        Chromosome name. If set, only this chromosome will be iterated and reported_positions will be derived from<br/>        the respective parameter. Note, that either chromosome and reported_positions or region must be set.<br/>    reported_positions : range or set<br/>        Range or set of genomic positions for which counts will be reported. The chromosome is derived from the<br/>        respective parameter. Note, that either chromosome and reported_positions or region must be set.<br/>    region : gi<br/>        Genomic region. If set, chromosome and reported_positions will be derived from this region.<br/>        Note, that either chromosome and reported_positions or region must be set.<br/>    file_format : str<br/>        optional, if set, the iterator will assume this file format (e.g., &#x27;bam&#x27; or &#x27;sam&#x27;)<br/>    min_mapping_quality : int<br/>        Filters pileup reads based on minimum mapping quality<br/>    flag_filter : int<br/>        Filters pileup reads based on read flags (see utils.BamFlag for details)<br/>    tag_filters : list<br/>        Filters pileup reads based on BAM tags (see utils.TagFilter for details)<br/>    max_span : int<br/>        Restricts maximum pileup depth.<br/>    min_base_quality : int<br/>        Filters pileup based on minimum per-base quality<br/>    fun_alias : function<br/>        Optional function for aliasing chromosome names<br/><br/>    Returns<br/>    -------<br/>    A &lt;base&gt;:&lt;count&gt; Counter object. Base is &#x27;None&#x27; for deletions.<br/><br/>    Notes<br/>    -----<br/>    * initial performance tests that used a synchronized iterator over a FASTA and 2 BAMs showed ~50kpos/sec<br/>    * reported stats:<br/>        iterated_items, chromosome: (int, str)<br/>            Number of iterated items (reads)<br/>        yielded_items, chromosome: (int, str)<br/>            Number of yielded items (positions)<br/><br/>    Examples<br/>    --------<br/>    &gt;&gt;&gt; pileup_data = rna.it(&quot;test.bam&quot;, style=&quot;pileup&quot;, region=gi(&#x27;chr1:1-100&#x27;)).to_list()<br/>    &gt;&gt;&gt; rna.it(&quot;test.bam&quot;, style=&#x27;pileup&#x27;,<br/>    &gt;&gt;&gt;        region=&#x27;1:22418244-22418244&#x27;, # pileup single position<br/>    &gt;&gt;&gt;        min_mapping_quality=30, flag_filter=0, min_base_quality=40, # read filtering<br/>    &gt;&gt;&gt;        tag_filters=[rna.TagFilter(&#x27;NM&#x27;, [0], True, inverse=True)], # include only reads w/o mismatches<br/>    &gt;&gt;&gt;        ).to_list()<br/><br/>    </p>\";\n",
       "        win.document.title = \"FastPileupIterator docstring\";\n",
       "        win.document.body.style.background = \"#dcf6fa\";\n",
       "        });\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rna.display_help(rna.FastPileupIterator)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f60a2b96-98af-40a4-98f1-f853e7bcaeea",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Feature Annotation\n",
    "\n",
    "Data integration is a central use case for genomic libraries including *rnalib*. This often requires the annotation of genomic features (e.g., genes or exons) with data derived from other, overlapping genomic intervals. One example is to calculate a mean genomic score based on data associated with all overlapping intervals from one or multiple other dataset(s).\n",
    "\n",
    "There are many different approaches to achieve these tasks:\n",
    "- in [pybedtools](https://daler.github.io/pybedtools), which is based on BEDTools, you could:\n",
    "    - [intersect](https://daler.github.io/pybedtools/autodocs/pybedtools.bedtool.BedTool.intersect.html#pybedtools.bedtool.BedTool.intersect) the respective interval sets\n",
    "    - apply a custom annotation method via [map()](https://daler.github.io/pybedtools/autodocs/pybedtools.bedtool.BedTool.map.html#pybedtools.bedtool.BedTool.map) or [each()](https://daler.github.io/pybedtools/autodocs/pybedtools.bedtool.BedTool.each.html#pybedtools.bedtool.BedTool.each) method\n",
    "- In [bioframe](https://bioframe.readthedocs.io/), which is based on pandas dataframes, you could:\n",
    "    - first use [overlap()](https://bioframe.readthedocs.io/en/latest/guide-intervalops.html#overlap) to intersect the interval sets\n",
    "    - then group on the annotated feature location using pandas [groupby](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.groupby.html) method\n",
    "    - and finally apply some [user defined mutate methods](https://pandas.pydata.org/docs/user_guide/gotchas.html#gotchas-udf-mutation) via apply() or aggregate().\n",
    "\n",
    "\n",
    "*Rnalib* transcriptomes contain an `anno` dict (a specialized _defaultdict_) for mapping (frozen) features (e.g., genes or exons) to arbitray annotations, e.g., float values, numpy arrays or custom python objects. Annotations can be added manually or by using the `transcriptome.annotate()` method that implements a generic method to add feature-based annotations based on rnalib's [`AnnotationIterator`](#AnnotationIterator) as shown below. The following section demonstrates some use cases:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "6a911370-0c34-48b3-bce6-c3e3114568eb",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Transcriptome with 5 genes and 89 tx"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# First, let's instantiate another small transcriptome that contains only a few human genes (near ACTB and SOX2) \n",
    "t = rna.Transcriptome(\n",
    "    genome_fa=rna.get_resource('ACTB+SOX2_genome'),       # reference genome FASTA\n",
    "    genome_offsets= {'chr3': 181711825, 'chr7': 5526309}, # This is a sliced reference sequence. To match the gencode annotations we provide respective offset coordinates for the 2 chroms\n",
    "    annotation_gff=rna.get_resource('gencode_gff'),       # gene annotation GFF\n",
    "    annotation_flavour='gencode',                         # this is a gencode flavoured GFF\n",
    "    load_sequence_data=False,                              # load gene sequences from reference genome)\n",
    "    disable_progressbar=True\n",
    ")\n",
    "display(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27d8f063-f3de-42c0-86a0-bb2df9ceb789",
   "metadata": {},
   "source": [
    "## Manual annotation\n",
    "We can now manually set annotation values and access them via 'dot' notation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "3bdcedc9-0c6d-48c4-a239-ec392fd7a526",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 2, 3)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "t.anno[t['ACTB']]['test']=(1,2,3) # add a \"test\" property to gene 'ACTB' and assign a tuple as value\n",
    "display(t['ACTB'].test) # access the annotation via 'dot' notation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cd2d230-ccd0-4c88-a336-520654c4f548",
   "metadata": {
    "tags": []
   },
   "source": [
    "## AnnotationIterator\n",
    "\n",
    "`rnalib` contains a special LocationIterator called `AnnotationIterator` for annotating genomic intervals with data from one or *multiple* other genomic location iterators.\n",
    "This iterator synchronizes genomic locations of the primary iterator (iterating the intervals to be annotated) with all annotating iterators and yields respective locations and data for all overlapping intervals.\n",
    "\n",
    "In the example below, we annotate intervals from a BED file with the sum of scores or all overlapping intervals from a bedgraph file.\n",
    "As all rnalib `LocationIterators` yield their current genomic interval in a uniform format, it is straightforward to use interval arithmetic in the annotation method (in this example we calculate the overlap fraction to calculate the actual score contributions of overlapping bedgraph intervals)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "ccb0020d-9a3d-4db5-a514-86ab4c6a8f66",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rnalib annotation iterator wrapping a Bed- and a BedGraphiterator\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'Annotation 1:6-10 (int1), sum:1.408'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Annotation 1:10-15 (int2), sum:0.3'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Annotation 2:10-150 (int3), sum:0'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Counter({('yielded_items', '1'): 2, ('yielded_items', '2'): 1})"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Annotate all intervals in a BED file with sum of scores from a bedgraph file; \n",
    "# supports bedgraph intervals >1bp: \n",
    "# score contribution is calculated from the interval overlap (this bedgraph file contains, e.g., an interval (1:7-10, 0.3))\n",
    "print('rnalib annotation iterator wrapping a Bed- and a BedGraphiterator')\n",
    "with rna.it(rna.it(rna.get_resource('test_bed')), # iterate this bed file\n",
    "            anno_its=rna.it(rna.get_resource('test_bedgraph')), # annotate it with overlapping intervals from this bedgraph file\n",
    "            labels=['scores'], # adds a 'scores' annotation that contains the values from the bedgraph file\n",
    "            disable_progressbar=True) as it:\n",
    "    for loc, i in it:\n",
    "        # \n",
    "        display(f\"Annotation {i.anno}, sum:{sum([x.data*loc.overlap(x.location) for x in i.scores])}\")\n",
    "    display(it.stats())  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bf7f778-0f91-428d-89c7-47ad4528284c",
   "metadata": {},
   "source": [
    " In this second, more complex example, we iterate over annotations in a genomic subregion of a flybase GTF file and annotate all features with \n",
    "(i) a list of SNPs that were called in 3 particular samples of a multi-sample VCF file and (ii) values from bedgraph file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "79113d3d-12d6-4b1b-a230-ec9f51598e3d",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'gene@2L:7529-9484 (+) NA (gene CR11023)'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'SNPs: [2L:7556T>C, 2L:7902A>G, 2L:8124A>C, 2L:8263G>A, 2L:9204T>C]'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"3'-value: 870, 5'-value: 319\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'pseudogene@2L:7529-9484 (+) FBtr0475186 (gene CR11023)'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'SNPs: [2L:7556T>C, 2L:7902A>G, 2L:8124A>C, 2L:8263G>A, 2L:9204T>C]'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"3'-value: 870, 5'-value: 319\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'exon@2L:7529-8116 (+) FBtr0475186 (gene CR11023)'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'SNPs: [2L:7556T>C, 2L:7902A>G]'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"3'-value: 870, 5'-value: 146\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'exon@2L:8193-9484 (+) FBtr0475186 (gene CR11023)'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'SNPs: [2L:8263G>A, 2L:9204T>C]'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"3'-value: 622, 5'-value: 319\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'exon@2L:9839-11344 (-) FBtr0078171 (gene l(2)gl)'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'SNPs: [2L:10018G>A, 2L:10089G>A, 2L:10610G>T]'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"3'-value: 932, 5'-value: 856\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Annotate flybase GTF annotation with (i) a list of the genotypes of overlapping SNPs for 3 samples and (ii) the maximum value from a bedgraph track\n",
    "with rna.it(\n",
    "    rna.it(rna.get_resource('flybase_gtf'), region='2L:1-30000'), # iterate over 2L:1-30000 in the flybase GTF\n",
    "           anno_its =[ rna.it(rna.get_resource('dmel_multisample_vcf'), # annotate with overlapping SNPs but consider 3 samples only\n",
    "                              samples=['DGRP-208', 'DGRP-325', 'DGRP-721']),\n",
    "                      rna.it(rna.get_resource('dmel_randomvalues')) ], # additionallly annotate with a bedgraph file\n",
    "    disable_progressbar=True) as it:\n",
    "    for loc, (gff_feature, overlapping_snps, rand_values) in it.to_list()[:5]: # here we show just the first 5 features...\n",
    "        # Remember that rnalib iterators return named (location, data) tuples. Here, 'snps' is a list of such \n",
    "        # items and we can access the respective VCFRecord objects via the data field:\n",
    "        overlapping_snps=[x.data for x in overlapping_snps]\n",
    "        # rand_values is a sorted list (by genomic location) of overlapping values from the Bedgraph file.\n",
    "        # We use min/max to get the first/last items from this list and access the value via the \n",
    "        # data field and convert from float (as returned by the bedgraph iterator) to int:        \n",
    "        tp_value, fp_value = int(min(rand_values).data), int(max(rand_values).data) \n",
    "        # now we can display some data\n",
    "        display(f\"{gff_feature['feature_type']}@{loc} {gff_feature.get('transcript_id','NA')} (gene {gff_feature['gene_symbol']})\")\n",
    "        display(f\"SNPs: {overlapping_snps}\")\n",
    "        display(f\"3'-value: {tp_value}, 5'-value: {fp_value}\", SEP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8642fc79-f9c1-4877-a68b-86da460919d6",
   "metadata": {
    "tags": []
   },
   "source": [
    "**Note** that iterating over large interval sets is generally considered slow and is often not needed. \n",
    "In many cases you can use list [comprehension, vectorisation](https://stackoverflow.com/questions/16476924/how-to-iterate-over-rows-in-a-pandas-dataframe/)\n",
    "or mechanisms such as [pandas apply()](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.apply.html) or [pybedtools each()](https://daler.github.io/pybedtools/each.html) as described above. \n",
    "Nevertheless, there are scenarios in which we consider AnnotationIterators  useful, e.g., if you need to access to contextual data (e.g., annotation results from previous data rows) or if you need (fine-grained) access to multiple synchronized datasets in parallel. Performance drawbacks of iteration approaches are also diminished when integrating multiple iterators in parallel. The overall goal of AnnotationIterators is to improve code readability and reduce potential error sources such as [the ones discussed in the 'RelatedWork' notebook](RelatedWork_performance.ipynb#Potential-pifalls). AnnotationIterators are also used by rnalib's [`transcriptome.annotate()`](#Annotation-of-features) method for incremental transcriptome annotation as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "500d6d90-0876-4d35-b12b-98f217f6f008",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<a id=\"help_icon0.9410094463914898\"><img id=\"help_icon\" style=\"cursor: help; float:left; width: 20px\" src=\"help_icon.png\"/>&nbsp;&nbsp;&nbsp;<b>TranscriptomeIterator docstring</b></a>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "\n",
       "    document.getElementById(\"help_icon0.9410094463914898\").addEventListener(\"click\", function() {\n",
       "        var win = window.open(\"\", \"\", \"toolbar=no, location=no, directories=no, status=no, menubar=no, scrollbars=yes, resizable=yes, width=1000, height=600\");\n",
       "        win.document.body.innerHTML = \"<h1>TranscriptomeIterator docstring</h1><p style='font-family:Courier New; font-size: 14px;'><br/>    Iterates over (sorted) features in a transcriptome object. Yielded items contain the respective feature<br/>    (location) and its transcriptome annotation (data).<br/>    Note that no chromosome aliasing is possible with this iterator as returned features are immutable.<br/><br/>    The transcriptome features will be filtered by feature_type and sorted if the transcriptome object<br/>    is not sorted (i.e., if features were added after it was built from a (sorted) GFF3 file).<br/><br/>    A transcriptome iterator can be converted to a pandas dataframe via the to_dataframe() method<br/>    which consumes the remaining items. Column names of the dataframe are automatically estimated<br/>    from the transcriptome model. Data columns may be excluded via the &#x27;excluded_columns&#x27; parameter;<br/>    by default the reserved &#x27;dna_seq&#x27; field is excluded as this field may contain very long sequence strings.<br/>    The default conversion function (fun) extracts values from the iterated location objects (=feature).<br/>    To add custom columns, one needs to provide (i) a list of column names via the &#x27;included_columns&#x27; parameter and<br/>    (ii) a custom extraction methof (fun) that return the respective values. Example:<br/><br/>    Examples<br/>    --------<br/>    &gt;&gt;&gt; # advanced usage: convert to pandas dataframe with a custom conversion function<br/>    &gt;&gt;&gt; # here we add a &#x27;feature length&#x27; column<br/>    &gt;&gt;&gt; def my_fun(tx, item, fun_col, default_value):<br/>    &gt;&gt;&gt;     return [len(tx) if col == &#x27;feature_len&#x27; else tx.get(col, default_value) for col in fun_col] # noqa<br/>    &gt;&gt;&gt; t = Transcriptome(...)<br/>    &gt;&gt;&gt; TranscriptomeIterator(t).to_dataframe(fun=my_fun, included_annotations=[&#x27;feature_len&#x27;]).head()<br/><br/>    Parameters<br/>    ----------<br/>    transcriptome : Transcriptome<br/>        The transcriptome object to iterate over<br/>    region : gi or str<br/>        optional, if set, only features overlapping with this region will be iterated<br/>    feature_types : list<br/>        optional, if set, only features of these types will be iterated<br/>    features : list<br/>        optional, if set, only these features will be iterated<br/>    Notes<br/>    -----<br/>    * reported stats:<br/>        iterated_items, chromosome: (int, str)<br/>            Number of iterated items<br/>        yielded_items, chromosome: (int, str)<br/>            Number of yielded items (remaining are filtered, e.g., due to region constraints)<br/>    </p>\";\n",
       "        win.document.title = \"TranscriptomeIterator docstring\";\n",
       "        win.document.body.style.background = \"#dcf6fa\";\n",
       "        });\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rna.display_help(rna.TranscriptomeIterator)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47b1fd83-a449-446f-bdc7-a23fddf19213",
   "metadata": {},
   "source": [
    "## Transcriptome annotate function\n",
    "`Rnalib` transcriptomes contain an `annotate()` method that creates a `TranscriptomeIterator` over its features and annotates it with a list of LocationIterators, applying a custom annotation method as demonstrated in this example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "7f4904e0-5371-4cff-8925-5a4e9a55a29c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "63ad492068b540d996898968f185d7e3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's annotate exons and introns with mappability scores that are read from a bedgraph file.\n",
    "# We first create a method for calculating the mean mappability for an annotation based on the scores of\n",
    "# all overlapping intervals from the bedgraph file.\n",
    "def calc_mean(item, label='mappability'):\n",
    "    loc, (anno, scores) = item\n",
    "    anno[label]=sum([score*loc.overlap(sloc) for sloc,score in scores])/len(loc)\n",
    "    return calc_mean\n",
    "\n",
    "# now we can call the annotate method, tell it to get the scores from a bedgraph iterator and\n",
    "# annotate exons + introns in our transcriptome. Note that this test bedgraph file covers only ~half of ACTB,\n",
    "# (use a genome browser like IGV to confirm this). So the mappability of the uncovered regions will \n",
    "# be set to zero.\n",
    "t.annotate(anno_its=rna.it(rna.get_resource('human_umap_k24')),\n",
    "           fun_anno=calc_mean,\n",
    "           feature_types=['exon', 'intron'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "994e0c82-4e98-4b17-973f-a686d17de362",
   "metadata": {},
   "source": [
    "Now we can access the annotated data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "e6caaec5-dc5c-4b87-ac01-76dfef306837",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Mappability for a random ACTB transcript:'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[(1, 'ENST00000675515.1_exon_5', 0.0),\n",
       " (2, 'ENST00000675515.1_exon_4', 0.8161705426356589),\n",
       " (3, 'ENST00000675515.1_exon_3', 0.28020416666666664),\n",
       " (4, 'ENST00000675515.1_exon_2', 0.0),\n",
       " (5, 'ENST00000675515.1_exon_1', 0.0),\n",
       " (6, 'ENST00000675515.1_exon_0', 0.0)]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Show the mappability values for exons and introns of one random transcript \n",
    "rnd_tx = random.sample(t.gene['ACTB'].transcript, 1)[0]\n",
    "display(SEP, \"Mappability for a random ACTB transcript:\", [(f.rnk, f.feature_id, f.mappability) for f in rnd_tx.exon])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "c0df59f0-7e7c-4bca-b763-db83dcff933c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Some annotated values:'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'exon mappability value: 0.9702537313432836'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'gene mappability value: None'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Its recommended to use get/getattr to provide sensible default values for missing data:\n",
    "display(SEP, \"Some annotated values:\")\n",
    "display(f\"exon mappability value: {t['ENST00000477812.2'].exon[0].get('mappability',None)}\") # show mappability score for 1st exon \n",
    "display(f\"gene mappability value: {t['ACTB'].get('mappability',None)}\") # only exons/introns were annotated above"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30bc296e-756f-4d0b-a794-1706d5d75a5e",
   "metadata": {},
   "source": [
    "Note that it would be more efficient to annotate genes with mappability score arrays and then calculate mean/median mappability\n",
    "for each feature by slicing these arrays (as done, e.g., for sequences) but this is omitted here for the sake of simplicity.\n",
    "Refer to the [AdvancedUsage notebook](AdvancedUsage.ipynb) for an example. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14a2707e-8f56-4ab4-9018-f543a3239200",
   "metadata": {},
   "source": [
    "We can also convert results to a dataframe (but note that we lose the parent/child relationships here) and further \n",
    "process the data with pandas. Refer to the [AdvancedUsage notebook](AdvancedUsage.ipynb) for more deatils on the `to_dataframe()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "875e655e-3c5b-4716-811c-9f57b95689ca",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\"Dataframe representation, filtered for gene_type=='protein_coding' & mappability>0\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Chromosome</th>\n",
       "      <th>Start</th>\n",
       "      <th>End</th>\n",
       "      <th>Strand</th>\n",
       "      <th>feature_id</th>\n",
       "      <th>feature_type</th>\n",
       "      <th>mappability</th>\n",
       "      <th>test</th>\n",
       "      <th>gene_type</th>\n",
       "      <th>gene_name</th>\n",
       "      <th>source</th>\n",
       "      <th>gff_feature_type</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>119</th>\n",
       "      <td>chr7</td>\n",
       "      <td>5527784</td>\n",
       "      <td>5529570</td>\n",
       "      <td>-</td>\n",
       "      <td>ENST00000676319.1_intron_0</td>\n",
       "      <td>intron</td>\n",
       "      <td>0.373112</td>\n",
       "      <td>None</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>None</td>\n",
       "      <td>HAVANA</td>\n",
       "      <td>exon</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>203</th>\n",
       "      <td>chr7</td>\n",
       "      <td>5528281</td>\n",
       "      <td>5529059</td>\n",
       "      <td>-</td>\n",
       "      <td>ENST00000477812.2_exon_2</td>\n",
       "      <td>exon</td>\n",
       "      <td>0.432606</td>\n",
       "      <td>None</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>None</td>\n",
       "      <td>HAVANA</td>\n",
       "      <td>exon</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>204</th>\n",
       "      <td>chr7</td>\n",
       "      <td>5528281</td>\n",
       "      <td>5529400</td>\n",
       "      <td>-</td>\n",
       "      <td>ENST00000462494.5_exon_2</td>\n",
       "      <td>exon</td>\n",
       "      <td>0.451115</td>\n",
       "      <td>None</td>\n",
       "      <td>protein_coding</td>\n",
       "      <td>None</td>\n",
       "      <td>HAVANA</td>\n",
       "      <td>exon</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Chromosome    Start      End Strand                  feature_id  \\\n",
       "119       chr7  5527784  5529570      -  ENST00000676319.1_intron_0   \n",
       "203       chr7  5528281  5529059      -    ENST00000477812.2_exon_2   \n",
       "204       chr7  5528281  5529400      -    ENST00000462494.5_exon_2   \n",
       "\n",
       "    feature_type  mappability  test       gene_type gene_name  source  \\\n",
       "119       intron     0.373112  None  protein_coding      None  HAVANA   \n",
       "203         exon     0.432606  None  protein_coding      None  HAVANA   \n",
       "204         exon     0.451115  None  protein_coding      None  HAVANA   \n",
       "\n",
       "    gff_feature_type  \n",
       "119             exon  \n",
       "203             exon  \n",
       "204             exon  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "---"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGzCAYAAAA1yP25AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8WgzjOAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1+0lEQVR4nO3deXxU1f3/8feQZSYsGUiALCVA2HepoUIEimJsRKQs+boiBKHFakQlX7XyVUDWACpSJYC1GLCCVCjSKhCEsFg1oKyuZZEtGBJQIWGRJCTn94cP5ueYBJIxyczF1/PxuI8H99wz935yJpA3556bsRljjAAAACyolrcLAAAA8BRBBgAAWBZBBgAAWBZBBgAAWBZBBgAAWBZBBgAAWBZBBgAAWBZBBgAAWBZBBgAAWBZBBkCl2Ww2PfTQQ1fst2jRItlsNh0+fNjVdsMNN+iGG25w7R8+fFg2m02LFi2q+kJ/4ufUDcA3EWQA+Jw1a9bomWee8XYZlXL+/Hk988wz2rx5s7dLAX5RCDIAqs2wYcP0/fffq1mzZuX2adasmb7//nsNGzbM1bZmzRpNmjSpJkosU0Xq/qnz589r0qRJBBmghvl7uwAAVy8/Pz/5+fldto/NZpPD4aihiiqmInX7mpKSEhUWFvrcWALVjRkZoAo988wzstls2rdvn+699145nU41atRI48ePlzFGWVlZGjhwoIKDgxUeHq7nn3/e7fWFhYWaMGGCYmJi5HQ6VadOHfXu3VubNm1y63dpXclzzz2nF154Qc2aNVNQUJD69Omjzz77zK3viBEjVLduXR08eFDx8fGqU6eOIiMjNXnyZBlj3Po+99xzuv766xUaGqqgoCDFxMRoxYoV5X69S5YsUdu2beVwOBQTE6P33nvP7XhF1pr8dI3MiBEjlJqaKumHkHNpM8aoefPmGjhwYKlzXLhwQU6nU/fff3+51/mxVatWqVOnTrLb7erYsaPS09OvWPf27dsVHx+vhg0bKigoSNHR0Ro5cqTra2jUqJEkadKkSa6af3x7bOPGjerdu7fq1Kmj+vXra+DAgfryyy9L1bZ582Z169ZNDodDLVu21Msvv+z6vvqxS+t9lixZoo4dO8put7u+joq+j5fOsXz5cnXo0EFBQUGKjY3Vp59+Kkl6+eWX1apVKzkcDt1www2sGYJvMgCqzMSJE40k07VrV3P33XebefPmmf79+xtJZvbs2aZt27bmgQceMPPmzTM9e/Y0ksyWLVtcrz958qSJiIgwycnJZv78+WbWrFmmbdu2JiAgwOzatcvV79ChQ0aS6dy5s2nevLmZOXOmmTRpkgkJCTGNGjUyOTk5rr6JiYnG4XCY1q1bm2HDhpm5c+ea2267zUgy48ePd6u/SZMm5sEHHzRz5841s2fPNtddd52RZN555x23fpJMp06dTMOGDc3kyZPNzJkzTbNmzUxQUJD59NNPXf3S0tKMJHPo0CFXW58+fUyfPn1KfS1paWnGGGM+/PBDc/PNNxtJ5u9//7trM8aYp556ygQEBJhvv/3WrZ4333zTSDLvvffeZd8fSeaaa64xERERZsqUKWbOnDmmRYsWpnbt2uabb74pt+7c3FzToEED06ZNG/Pss8+aV155xTz11FOmffv2xhhjzp49a+bPn28kmcGDB7tq3rNnjzHGmPXr1xt/f3/Tpk0bM2vWLDNp0iTTsGFD06BBA7ex2blzp7Hb7aZ58+ZmxowZZtq0aSYyMtJcc8015qf/XEsy7du3N40aNTKTJk0yqampru+RyryPXbp0MVFRUWbGjBlmxowZxul0mqZNm5q5c+eaDh06mOeff948/fTTJjAw0Nx4442XHV/AGwgyQBW6FGRGjx7tart48aJp0qSJsdlsZsaMGa72U6dOmaCgIJOYmOjWt6CgwO2cp06dMmFhYWbkyJGutks//IOCgsyxY8dc7du2bTOSzNixY11tiYmJRpIZM2aMq62kpMT079/fBAYGmpMnT7raz58/73btwsJC06lTJ9O3b1+3dklGktm+fbur7ciRI8bhcJjBgwe72jwJMsYYk5SUVOoHtzHG7N2710gy8+fPd2v//e9/b5o3b25KSkpKveandQcGBpoDBw642vbs2WMkmZdeeqncut966y0jyXz88cflnvvkyZNGkpk4cWKpY127djWNGzd2C2B79uwxtWrVMsOHD3e1DRgwwNSuXdt8/fXXrrb9+/cbf3//MoNMrVq1zOeff17qepV5H+12u9v78/LLLxtJJjw83OTn57vax40bV+q9BHwBt5aAavCHP/zB9Wc/Pz9169ZNxhiNGjXK1V6/fn21bdtWBw8edOsbGBgo6Yc1D999950uXryobt26aefOnaWuM2jQIP3qV79y7V933XXq3r271qxZU6rvjx87vnRLobCwUBs2bHC1BwUFuf586tQp5eXlqXfv3mVeOzY2VjExMa79pk2bauDAgVq3bp2Ki4vLH5yfoU2bNurevbuWLFniavvuu++0du1aDR06tNTtl7LExcWpZcuWrv0uXbooODjY7X34qfr160uS3nnnHRUVFVWq5uPHj2v37t0aMWKEQkJC3K578803u96r4uJibdiwQYMGDVJkZKSrX6tWrdSvX78yz92nTx916NChVHtl3sebbrpJzZs3d+13795dkpSQkKB69eqVar/cOAHeQJABqkHTpk3d9p1OpxwOhxo2bFiq/dSpU25tixcvVpcuXeRwOBQaGqpGjRpp9erVysvLK3Wd1q1bl2pr06ZNqbUMtWrVUosWLUr1k+TW95133lGPHj3kcDgUEhKiRo0aaf78+ZW69vnz53Xy5MlSx6rK8OHD9cEHH+jIkSOSpOXLl6uoqMjtqafL+el7I0kNGjQo9T78WJ8+fZSQkKBJkyapYcOGGjhwoNLS0lRQUHDF612qs23btqWOtW/fXt98843OnTunEydO6Pvvv1erVq1K9SurTZKio6PLbK/M+1jW96okRUVFldl+uXECvIEgA1SDsp54Ke8pGPOjBbevv/66RowYoZYtW2rhwoVKT0/X+vXr1bdvX5WUlFRbvZL0n//8R7///e/lcDg0b948rVmzRuvXr9c999xTalGwN911110KCAhwzcq8/vrr6tatW5lBoSwVeR9+ymazacWKFcrMzNRDDz2kr7/+WiNHjlRMTIzOnj1b+S+iivx45uWSyr6P5Y2HJ+MEeAOPXwM+ZMWKFWrRooVWrlzpdptk4sSJZfbfv39/qbZ9+/a53SqQfrhNdfDgQdcszKV+klx9//nPf8rhcGjdunWy2+2ufmlpaZW6du3atV1P8HjqcreIQkJC1L9/fy1ZskRDhw7VBx98oDlz5vys61VUjx491KNHD02bNk1Lly7V0KFDtWzZMv3hD38ot+ZLv4tm7969pY7997//VcOGDVWnTh05HA45HA4dOHCgVL+y2spT2fcRsDpmZAAfcul/wT/+X++2bduUmZlZZv9Vq1bp66+/du1/9NFH2rZtW5lrKubOnev6szFGc+fOVUBAgG666SbXtW02m9v6lsOHD2vVqlVlXjszM9NtzUVWVpb+9a9/6Xe/+93P/h0sderUkSSdPn26zOPDhg3TF198occff1x+fn666667ftb1ruTUqVOlZiK6du0qSa7bS7Vr15ZUuuaIiAh17dpVixcvdjv22Wef6d1339Wtt94q6Yfxj4uL06pVq5Sdne3qd+DAAa1du7bCtVb2fQSsjhkZwIfcdtttWrlypQYPHqz+/fvr0KFDWrBggTp06FDmLYxWrVqpV69eeuCBB1RQUKA5c+YoNDRUTzzxhFs/h8Oh9PR0JSYmqnv37lq7dq1Wr16t//u//3PNnvTv31+zZ8/WLbfconvuuUcnTpxQamqqWrVqpU8++aTUtTt16qT4+Hg9/PDDstvtmjdvniRVyW/kvbSI+OGHH1Z8fHypsNK/f3+FhoZq+fLl6tevnxo3bvyzr3k5ixcv1rx58zR48GC1bNlSZ86c0SuvvKLg4GBXEAkKClKHDh30j3/8Q23atFFISIg6deqkTp066dlnn1W/fv0UGxurUaNG6fvvv9dLL70kp9Pp9rtmnnnmGb377rvq2bOnHnjgARUXF2vu3Lnq1KmTdu/eXaFaK/s+AlbHjAzgQ0aMGKHp06drz549evjhh7Vu3TrXGpCyDB8+XGPGjNHcuXM1bdo0dezYURs3blRERIRbPz8/P6WnpysnJ0ePP/64Pv74Y02cOFFTpkxx9enbt68WLlyonJwcPfroo3rjjTc0c+ZMDR48uMxr9+nTR3PmzNHf//53TZgwQSEhIVq7dq26dOnys8dhyJAhGjNmjNLT0zVs2DDdfffdbscDAwN15513SlKFF/n+HH369FG3bt20bNkyPfzww5o1a5Zat26tjRs3ui24/dvf/qZf/epXGjt2rO6++27XL6GLi4tTenq6QkNDNWHCBD333HPq0aOHPvjgA7fXx8TEaO3atWrQoIHGjx+vhQsXavLkybrpppsq/Bt7K/s+AlZnM6zcAizn8OHDio6O1rPPPqvHHnvssn1HjBihFStWeHVRanUYO3as6wf2pds6V6tBgwbp888/L3NdEvBLx4wMAMu5cOGCXn/9dSUkJFx1Ieb7779329+/f7/WrFmjG264wTsFAT6ONTIALOPEiRPasGGDVqxYoW+//VaPPPKIt0uqci1atNCIESPUokULHTlyRPPnz1dgYGCpdU8AfkCQAWAZX3zxhYYOHarGjRvrxRdfdD05dDW55ZZb9MYbbygnJ0d2u12xsbGaPn16mb+AEABrZAAAgIWxRgYAAFgWQQYAAFjWVb9GpqSkRNnZ2apXr16FPhkXAAB4nzFGZ86cUWRkpGrVKn/e5aoPMtnZ2aU+xRUAAFhDVlaWmjRpUu7xqz7I1KtXT9IPAxEcHOzlagAAQEXk5+crKirK9XO8PFd9kLl0Oyk4OJggAwCAxVxpWQiLfQEAgGURZAAAgGURZAAAgGURZAAAgGURZAAAgGURZAAAgGURZAAAgGURZAAAgGURZAAAgGURZAAAgGURZAAAgGURZAAAgGURZAAAgGURZAAAgGX5e7sAK2v+5Gpvl+CRwzP6e7sEAACqBDMyAADAsggyAADAsggyAADAsggyAADAsrwaZIqLizV+/HhFR0crKChILVu21JQpU2SMcfUxxmjChAmKiIhQUFCQ4uLitH//fi9WDQAAfIVXg8zMmTM1f/58zZ07V19++aVmzpypWbNm6aWXXnL1mTVrll588UUtWLBA27ZtU506dRQfH68LFy54sXIAAOALvPr49YcffqiBAweqf/8fHgdu3ry53njjDX300UeSfpiNmTNnjp5++mkNHDhQkvTaa68pLCxMq1at0l133eW12gEAgPd5dUbm+uuvV0ZGhvbt2ydJ2rNnj95//33169dPknTo0CHl5OQoLi7O9Rqn06nu3bsrMzOzzHMWFBQoPz/fbQMAAFcnr87IPPnkk8rPz1e7du3k5+en4uJiTZs2TUOHDpUk5eTkSJLCwsLcXhcWFuY69lMpKSmaNGlS9RYOAAB8gldnZN58800tWbJES5cu1c6dO7V48WI999xzWrx4scfnHDdunPLy8lxbVlZWFVYMAAB8iVdnZB5//HE9+eSTrrUunTt31pEjR5SSkqLExESFh4dLknJzcxUREeF6XW5urrp27VrmOe12u+x2e7XXDgAAvM+rMzLnz59XrVruJfj5+amkpESSFB0drfDwcGVkZLiO5+fna9u2bYqNja3RWgEAgO/x6ozMgAEDNG3aNDVt2lQdO3bUrl27NHv2bI0cOVKSZLPZ9Oijj2rq1Klq3bq1oqOjNX78eEVGRmrQoEHeLB0AAPgArwaZl156SePHj9eDDz6oEydOKDIyUvfff78mTJjg6vPEE0/o3LlzGj16tE6fPq1evXopPT1dDofDi5UDAABfYDM//jW6V6H8/Hw5nU7l5eUpODi4Ss/d/MnVVXq+mnJ4Rn9vlwAAwGVV9Oc3n7UEAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsy6tBpnnz5rLZbKW2pKQkSdKFCxeUlJSk0NBQ1a1bVwkJCcrNzfVmyQAAwId4Nch8/PHHOn78uGtbv369JOn222+XJI0dO1Zvv/22li9fri1btig7O1tDhgzxZskAAMCH+Hvz4o0aNXLbnzFjhlq2bKk+ffooLy9PCxcu1NKlS9W3b19JUlpamtq3b6+tW7eqR48e3igZAAD4EJ9ZI1NYWKjXX39dI0eOlM1m044dO1RUVKS4uDhXn3bt2qlp06bKzMws9zwFBQXKz8932wAAwNXJZ4LMqlWrdPr0aY0YMUKSlJOTo8DAQNWvX9+tX1hYmHJycso9T0pKipxOp2uLioqqxqoBAIA3+UyQWbhwofr166fIyMifdZ5x48YpLy/PtWVlZVVRhQAAwNd4dY3MJUeOHNGGDRu0cuVKV1t4eLgKCwt1+vRpt1mZ3NxchYeHl3suu90uu91eneUCAAAf4RMzMmlpaWrcuLH69+/vaouJiVFAQIAyMjJcbXv37tXRo0cVGxvrjTIBAICP8fqMTElJidLS0pSYmCh///9fjtPp1KhRo5ScnKyQkBAFBwdrzJgxio2N5YklAAAgyQeCzIYNG3T06FGNHDmy1LEXXnhBtWrVUkJCggoKChQfH6958+Z5oUoAAOCLbMYY4+0iqlN+fr6cTqfy8vIUHBxcpedu/uTqKj1fTTk8o/+VOwEA4EUV/fntE2tkAAAAPOH1W0sAAOAHVpzp9/YsPzMyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsrweZL7++mvde++9Cg0NVVBQkDp37qzt27e7jhtjNGHCBEVERCgoKEhxcXHav3+/FysGAAC+wqtB5tSpU+rZs6cCAgK0du1affHFF3r++efVoEEDV59Zs2bpxRdf1IIFC7Rt2zbVqVNH8fHxunDhghcrBwAAvsDfmxefOXOmoqKilJaW5mqLjo52/dkYozlz5ujpp5/WwIEDJUmvvfaawsLCtGrVKt111101XjMAAPAdXp2R+fe//61u3brp9ttvV+PGjfXrX/9ar7zyiuv4oUOHlJOTo7i4OFeb0+lU9+7dlZmZWeY5CwoKlJ+f77YBAICrk1eDzMGDBzV//ny1bt1a69at0wMPPKCHH35YixcvliTl5ORIksLCwtxeFxYW5jr2UykpKXI6na4tKiqqer8IAADgNV4NMiUlJbr22ms1ffp0/frXv9bo0aP1xz/+UQsWLPD4nOPGjVNeXp5ry8rKqsKKAQCAL/FqkImIiFCHDh3c2tq3b6+jR49KksLDwyVJubm5bn1yc3Ndx37KbrcrODjYbQMAAFcnrwaZnj17au/evW5t+/btU7NmzST9sPA3PDxcGRkZruP5+fnatm2bYmNja7RWAADge7z61NLYsWN1/fXXa/r06brjjjv00Ucf6a9//av++te/SpJsNpseffRRTZ06Va1bt1Z0dLTGjx+vyMhIDRo0yJulAwAAH+DVIPOb3/xGb731lsaNG6fJkycrOjpac+bM0dChQ119nnjiCZ07d06jR4/W6dOn1atXL6Wnp8vhcHixcgAA4Au8GmQk6bbbbtNtt91W7nGbzabJkydr8uTJNVgVAACwAq9/RAEAAICnCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyCDIAAMCyvBpknnnmGdlsNretXbt2ruMXLlxQUlKSQkNDVbduXSUkJCg3N9eLFQMAAF/iUZA5ePBglRXQsWNHHT9+3LW9//77rmNjx47V22+/reXLl2vLli3Kzs7WkCFDquzaAADA2vw9eVGrVq3Up08fjRo1Sv/zP/8jh8PheQH+/goPDy/VnpeXp4ULF2rp0qXq27evJCktLU3t27fX1q1b1aNHD4+vCQAArg4ezcjs3LlTXbp0UXJyssLDw3X//ffro48+8qiA/fv3KzIyUi1atNDQoUN19OhRSdKOHTtUVFSkuLg4V9927dqpadOmyszMLPd8BQUFys/Pd9sAAMDVyaMg07VrV/3lL39Rdna2Xn31VR0/fly9evVSp06dNHv2bJ08ebJC5+nevbsWLVqk9PR0zZ8/X4cOHVLv3r115swZ5eTkKDAwUPXr13d7TVhYmHJycso9Z0pKipxOp2uLiory5EsEAAAW8LMW+/r7+2vIkCFavny5Zs6cqQMHDuixxx5TVFSUhg8fruPHj1/29f369dPtt9+uLl26KD4+XmvWrNHp06f15ptvelzTuHHjlJeX59qysrI8PhcAAPBtPyvIbN++XQ8++KAiIiI0e/ZsPfbYY/rqq6+0fv16ZWdna+DAgZU6X/369dWmTRsdOHBA4eHhKiws1OnTp9365Obmlrmm5hK73a7g4GC3DQAAXJ08CjKzZ89W586ddf311ys7O1uvvfaajhw5oqlTpyo6Olq9e/fWokWLtHPnzkqd9+zZs/rqq68UERGhmJgYBQQEKCMjw3V87969Onr0qGJjYz0pGwAAXGU8empp/vz5GjlypEaMGKGIiIgy+zRu3FgLFy687Hkee+wxDRgwQM2aNVN2drYmTpwoPz8/3X333XI6nRo1apSSk5MVEhKi4OBgjRkzRrGxsTyxBAAAJHkYZPbv33/FPoGBgUpMTLxsn2PHjunuu+/Wt99+q0aNGqlXr17aunWrGjVqJEl64YUXVKtWLSUkJKigoEDx8fGaN2+eJyUDAICrkEdBJi0tTXXr1tXtt9/u1r58+XKdP3/+igHmkmXLll32uMPhUGpqqlJTUz0pEwAAXOU8WiOTkpKihg0blmpv3Lixpk+f/rOLAgAAqAiPgszRo0cVHR1dqr1Zs2auX2gHAABQ3TwKMo0bN9Ynn3xSqn3Pnj0KDQ392UUBAABUhEdB5u6779bDDz+sTZs2qbi4WMXFxdq4caMeeeQR3XXXXVVdIwAAQJk8Wuw7ZcoUHT58WDfddJP8/X84RUlJiYYPH84aGQAAUGM8CjKBgYH6xz/+oSlTpmjPnj0KCgpS586d1axZs6quDwAAoFweBZlL2rRpozZt2lRVLQAAAJXiUZApLi7WokWLlJGRoRMnTqikpMTt+MaNG6ukOAAAgMvxKMg88sgjWrRokfr3769OnTrJZrNVdV0AAABX5FGQWbZsmd58803deuutVV0PAABAhXn0+HVgYKBatWpV1bUAAABUikdB5n//93/1l7/8RcaYqq4HAACgwjy6tfT+++9r06ZNWrt2rTp27KiAgAC34ytXrqyS4gAAAC7HoyBTv359DR48uKprAQAAqBSPgkxaWlpV1wEAAFBpHq2RkaSLFy9qw4YNevnll3XmzBlJUnZ2ts6ePVtlxQEAAFyORzMyR44c0S233KKjR4+qoKBAN998s+rVq6eZM2eqoKBACxYsqOo6AQAASvFoRuaRRx5Rt27ddOrUKQUFBbnaBw8erIyMjCorDgAA4HI8mpH5z3/+ow8//FCBgYFu7c2bN9fXX39dJYUBAABciUczMiUlJSouLi7VfuzYMdWrV+9nFwUAAFARHgWZ3/3ud5ozZ45r32az6ezZs5o4cSIfWwAAAGqMR7eWnn/+ecXHx6tDhw66cOGC7rnnHu3fv18NGzbUG2+8UdU1AgAAlMmjINOkSRPt2bNHy5Yt0yeffKKzZ89q1KhRGjp0qNviXwAAgOrkUZCRJH9/f917771VWQsAAECleBRkXnvttcseHz58uEfFAAAAVIZHQeaRRx5x2y8qKtL58+cVGBio2rVrE2QAAECN8OippVOnTrltZ8+e1d69e9WrVy8W+wIAgBrj8Wct/VTr1q01Y8aMUrM1AAAA1aXKgoz0wwLg7OzsqjwlAABAuTxaI/Pvf//bbd8Yo+PHj2vu3Lnq2bNnlRQGAABwJR4FmUGDBrnt22w2NWrUSH379tXzzz9fFXUBAABckUdBpqSkpKrrAAAAqLQqXSPzc8yYMUM2m02PPvqoq+3ChQtKSkpSaGio6tatq4SEBOXm5nqvSAAA4FM8mpFJTk6ucN/Zs2dfsc/HH3+sl19+WV26dHFrHzt2rFavXq3ly5fL6XTqoYce0pAhQ/TBBx9UumYAAHD18SjI7Nq1S7t27VJRUZHatm0rSdq3b5/8/Px07bXXuvrZbLYrnuvs2bMaOnSoXnnlFU2dOtXVnpeXp4ULF2rp0qXq27evJCktLU3t27fX1q1b1aNHD09KBwAAVxGPbi0NGDBAv/3tb3Xs2DHt3LlTO3fuVFZWlm688Ubddttt2rRpkzZt2qSNGzde8VxJSUnq37+/4uLi3Np37NihoqIit/Z27dqpadOmyszMLPd8BQUFys/Pd9sAAMDVyaMg8/zzzyslJUUNGjRwtTVo0EBTp06t1FNLy5Yt086dO5WSklLqWE5OjgIDA1W/fn239rCwMOXk5JR7zpSUFDmdTtcWFRVV4XoAAIC1eBRk8vPzdfLkyVLtJ0+e1JkzZyp0jqysLD3yyCNasmSJHA6HJ2WUady4ccrLy3NtWVlZVXZuAADgWzwKMoMHD9Z9992nlStX6tixYzp27Jj++c9/atSoURoyZEiFzrFjxw6dOHFC1157rfz9/eXv768tW7boxRdflL+/v8LCwlRYWKjTp0+7vS43N1fh4eHlntdutys4ONhtAwAAVyePFvsuWLBAjz32mO655x4VFRX9cCJ/f40aNUrPPvtshc5x00036dNPP3Vru++++9SuXTv9+c9/VlRUlAICApSRkaGEhARJ0t69e3X06FHFxsZ6UjYAALjKeBRkateurXnz5unZZ5/VV199JUlq2bKl6tSpU+Fz1KtXT506dXJrq1OnjkJDQ13to0aNUnJyskJCQhQcHKwxY8YoNjaWJ5YAAIAkD4PMJcePH9fx48f129/+VkFBQTLGVOiR64p64YUXVKtWLSUkJKigoEDx8fGaN29elZ0fAABYm0dB5ttvv9Udd9yhTZs2yWazaf/+/WrRooVGjRqlBg0aePx5S5s3b3bbdzgcSk1NVWpqqkfnAwAAVzePFvuOHTtWAQEBOnr0qGrXru1qv/POO5Wenl5lxQEAAFyORzMy7777rtatW6cmTZq4tbdu3VpHjhypksIAAACuxKMZmXPnzrnNxFzy3XffyW63/+yiAAAAKsKjINO7d2+99tprrn2bzaaSkhLNmjVLN954Y5UVBwAAcDke3VqaNWuWbrrpJm3fvl2FhYV64okn9Pnnn+u7777jk6kBAECN8WhGplOnTtq3b5969eqlgQMH6ty5cxoyZIh27dqlli1bVnWNAAAAZar0jExRUZFuueUWLViwQE899VR11AQAAFAhlZ6RCQgI0CeffFIdtQAAAFSKR7eW7r33Xi1cuLCqawEAAKgUjxb7Xrx4Ua+++qo2bNigmJiYUp+xNHv27CopDgAA4HIqFWQOHjyo5s2b67PPPtO1114rSdq3b59bn6r8rCUAAIDLqVSQad26tY4fP65NmzZJ+uEjCV588UWFhYVVS3EAAACXU6k1MsYYt/21a9fq3LlzVVoQAABARXm02PeSnwYbAACAmlSpIGOz2UqtgWFNDAAA8JZKrZExxmjEiBGuD4a8cOGC/vSnP5V6amnlypVVVyEAAEA5KhVkEhMT3fbvvffeKi0GAACgMioVZNLS0qqrDgAAgEr7WYt9AQAAvIkgAwAALIsgAwAALIsgAwAALIsgAwAALIsgAwAALIsgAwAALIsgAwAALIsgAwAALIsgAwAALIsgAwAALIsgAwAALIsgAwAALIsgAwAALMurQWb+/Pnq0qWLgoODFRwcrNjYWK1du9Z1/MKFC0pKSlJoaKjq1q2rhIQE5ebmerFiAADgS7waZJo0aaIZM2Zox44d2r59u/r27auBAwfq888/lySNHTtWb7/9tpYvX64tW7YoOztbQ4YM8WbJAADAh/h78+IDBgxw2582bZrmz5+vrVu3qkmTJlq4cKGWLl2qvn37SpLS0tLUvn17bd26VT169CjznAUFBSooKHDt5+fnV98XAAAAvMpn1sgUFxdr2bJlOnfunGJjY7Vjxw4VFRUpLi7O1addu3Zq2rSpMjMzyz1PSkqKnE6na4uKiqqJ8gEAgBd4Pch8+umnqlu3rux2u/70pz/prbfeUocOHZSTk6PAwEDVr1/frX9YWJhycnLKPd+4ceOUl5fn2rKysqr5KwAAAN7i1VtLktS2bVvt3r1beXl5WrFihRITE7VlyxaPz2e322W326uwQgAA4Ku8HmQCAwPVqlUrSVJMTIw+/vhj/eUvf9Gdd96pwsJCnT592m1WJjc3V+Hh4V6qFgAA+BKv31r6qZKSEhUUFCgmJkYBAQHKyMhwHdu7d6+OHj2q2NhYL1YIAAB8hVdnZMaNG6d+/fqpadOmOnPmjJYuXarNmzdr3bp1cjqdGjVqlJKTkxUSEqLg4GCNGTNGsbGx5T6xBAAAflm8GmROnDih4cOH6/jx43I6nerSpYvWrVunm2++WZL0wgsvqFatWkpISFBBQYHi4+M1b948b5YMAAB8iFeDzMKFCy973OFwKDU1VampqTVUEQAAsBKfWyMDAABQUQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWQQZAABgWV4NMikpKfrNb36jevXqqXHjxho0aJD27t3r1ufChQtKSkpSaGio6tatq4SEBOXm5nqpYgAA4Eu8GmS2bNmipKQkbd26VevXr1dRUZF+97vf6dy5c64+Y8eO1dtvv63ly5dry5Ytys7O1pAhQ7xYNQAA8BX+3rx4enq62/6iRYvUuHFj7dixQ7/97W+Vl5enhQsXaunSperbt68kKS0tTe3bt9fWrVvVo0ePUucsKChQQUGBaz8/P796vwgAAOA1PrVGJi8vT5IUEhIiSdqxY4eKiooUFxfn6tOuXTs1bdpUmZmZZZ4jJSVFTqfTtUVFRVV/4QAAwCt8JsiUlJTo0UcfVc+ePdWpUydJUk5OjgIDA1W/fn23vmFhYcrJySnzPOPGjVNeXp5ry8rKqu7SAQCAl3j11tKPJSUl6bPPPtP777//s85jt9tlt9urqCoAAODLfGJG5qGHHtI777yjTZs2qUmTJq728PBwFRYW6vTp0279c3NzFR4eXsNVAgAAX+PVIGOM0UMPPaS33npLGzduVHR0tNvxmJgYBQQEKCMjw9W2d+9eHT16VLGxsTVdLgAA8DFevbWUlJSkpUuX6l//+pfq1avnWvfidDoVFBQkp9OpUaNGKTk5WSEhIQoODtaYMWMUGxtb5hNLAADgl8WrQWb+/PmSpBtuuMGtPS0tTSNGjJAkvfDCC6pVq5YSEhJUUFCg+Ph4zZs3r4YrBQAAvsirQcYYc8U+DodDqampSk1NrYGKAACAlfjEYl8AAABP+Mzj1wAAVKXmT672dgmoAczIAAAAyyLIAAAAyyLIAAAAyyLIAAAAyyLIAAAAyyLIAAAAyyLIAAAAyyLIAAAAyyLIAAAAyyLIAAAAyyLIAAAAyyLIAAAAyyLIAAAAyyLIAAAAyyLIAAAAyyLIAAAAy/L3dgFARTR/crW3S6i0wzP6e7sEALjqMSMDAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsiyADAAAsy6tB5r333tOAAQMUGRkpm82mVatWuR03xmjChAmKiIhQUFCQ4uLitH//fu8UCwAAfI5Xg8y5c+d0zTXXKDU1tczjs2bN0osvvqgFCxZo27ZtqlOnjuLj43XhwoUarhQAAPgif29evF+/furXr1+Zx4wxmjNnjp5++mkNHDhQkvTaa68pLCxMq1at0l133VXm6woKClRQUODaz8/Pr/rCAQCAT/DZNTKHDh1STk6O4uLiXG1Op1Pdu3dXZmZmua9LSUmR0+l0bVFRUTVRLgAA8AKfDTI5OTmSpLCwMLf2sLAw17GyjBs3Tnl5ea4tKyurWusEAADe49VbS9XBbrfLbrd7uwwAAFADfHZGJjw8XJKUm5vr1p6bm+s6BgAAftl8NshER0crPDxcGRkZrrb8/Hxt27ZNsbGxXqwMAAD4Cq/eWjp79qwOHDjg2j906JB2796tkJAQNW3aVI8++qimTp2q1q1bKzo6WuPHj1dkZKQGDRrkvaIBAIDP8GqQ2b59u2688UbXfnJysiQpMTFRixYt0hNPPKFz585p9OjROn36tHr16qX09HQ5HA5vlQwAAHyIV4PMDTfcIGNMucdtNpsmT56syZMn12BVAADAKnx2jQwAAMCVEGQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBlEWQAAIBl+Xu7ANS85k+u9nYJvwiMc804PKO/t0uoNCt+b1hxnPHLwIwMAACwLIIMAACwLIIMAACwLEsEmdTUVDVv3lwOh0Pdu3fXRx995O2SAACAD/D5IPOPf/xDycnJmjhxonbu3KlrrrlG8fHxOnHihLdLAwAAXubzQWb27Nn64x//qPvuu08dOnTQggULVLt2bb366qveLg0AAHiZTz9+XVhYqB07dmjcuHGutlq1aikuLk6ZmZllvqagoEAFBQWu/by8PElSfn5+lddXUnC+ys8JoHKq4+92dbPivx2MM8pTXd8bl85rjLlsP58OMt98842Ki4sVFhbm1h4WFqb//ve/Zb4mJSVFkyZNKtUeFRVVLTUC8C7nHG9X8MvAOKM81f29cebMGTmdznKP+3SQ8cS4ceOUnJzs2i8pKdF3332n0NBQ2Wy2KrtOfn6+oqKilJWVpeDg4Co7L0pjrGsG41wzGOeawTjXjOocZ2OMzpw5o8jIyMv28+kg07BhQ/n5+Sk3N9etPTc3V+Hh4WW+xm63y263u7XVr1+/ukpUcHAwf0lqCGNdMxjnmsE41wzGuWZU1zhfbibmEp9e7BsYGKiYmBhlZGS42kpKSpSRkaHY2FgvVgYAAHyBT8/ISFJycrISExPVrVs3XXfddZozZ47OnTun++67z9ulAQAAL/P5IHPnnXfq5MmTmjBhgnJyctS1a1elp6eXWgBc0+x2uyZOnFjqNhaqHmNdMxjnmsE41wzGuWb4wjjbzJWeawIAAPBRPr1GBgAA4HIIMgAAwLIIMgAAwLIIMgAAwLIIMgAAwLIIMpeRmpqq5s2by+FwqHv37vroo48u23/58uVq166dHA6HOnfurDVr1tRQpdZXmbF+5ZVX1Lt3bzVo0EANGjRQXFzcFd8b/KCy39OXLFu2TDabTYMGDareAq8SlR3n06dPKykpSREREbLb7WrTpg3/flRAZcd5zpw5atu2rYKCghQVFaWxY8fqwoULNVStNb333nsaMGCAIiMjZbPZtGrVqiu+ZvPmzbr22mtlt9vVqlUrLVq0qHqLNCjTsmXLTGBgoHn11VfN559/bv74xz+a+vXrm9zc3DL7f/DBB8bPz8/MmjXLfPHFF+bpp582AQEB5tNPP63hyq2nsmN9zz33mNTUVLNr1y7z5ZdfmhEjRhin02mOHTtWw5VbS2XH+ZJDhw6ZX/3qV6Z3795m4MCBNVOshVV2nAsKCky3bt3Mrbfeat5//31z6NAhs3nzZrN79+4artxaKjvOS5YsMXa73SxZssQcOnTIrFu3zkRERJixY8fWcOXWsmbNGvPUU0+ZlStXGknmrbfeumz/gwcPmtq1a5vk5GTzxRdfmJdeesn4+fmZ9PT0aquRIFOO6667ziQlJbn2i4uLTWRkpElJSSmz/x133GH69+/v1ta9e3dz//33V2udV4PKjvVPXbx40dSrV88sXry4ukq8KngyzhcvXjTXX3+9+dvf/mYSExMJMhVQ2XGeP3++adGihSksLKypEq8KlR3npKQk07dvX7e25ORk07Nnz2qt82pSkSDzxBNPmI4dO7q13XnnnSY+Pr7a6uLWUhkKCwu1Y8cOxcXFudpq1aqluLg4ZWZmlvmazMxMt/6SFB8fX25//MCTsf6p8+fPq6ioSCEhIdVVpuV5Os6TJ09W48aNNWrUqJoo0/I8Ged///vfio2NVVJSksLCwtSpUydNnz5dxcXFNVW25Xgyztdff7127Njhuv108OBBrVmzRrfeemuN1PxL4Y2fhT7/EQXe8M0336i4uLjUxyCEhYXpv//9b5mvycnJKbN/Tk5OtdV5NfBkrH/qz3/+syIjI0v95cH/58k4v//++1q4cKF2795dAxVeHTwZ54MHD2rjxo0aOnSo1qxZowMHDujBBx9UUVGRJk6cWBNlW44n43zPPffom2++Ua9evWSM0cWLF/WnP/1J//d//1cTJf9ilPezMD8/X99//72CgoKq/JrMyMDSZsyYoWXLlumtt96Sw+HwdjlXjTNnzmjYsGF65ZVX1LBhQ2+Xc1UrKSlR48aN9de//lUxMTG688479dRTT2nBggXeLu2qsnnzZk2fPl3z5s3Tzp07tXLlSq1evVpTpkzxdmn4mZiRKUPDhg3l5+en3Nxct/bc3FyFh4eX+Zrw8PBK9ccPPBnrS5577jnNmDFDGzZsUJcuXaqzTMur7Dh/9dVXOnz4sAYMGOBqKykpkST5+/tr7969atmyZfUWbUGefD9HREQoICBAfn5+rrb27dsrJydHhYWFCgwMrNaarciTcR4/fryGDRumP/zhD5Kkzp0769y5cxo9erSeeuop1arF/+urQnk/C4ODg6tlNkZiRqZMgYGBiomJUUZGhqutpKREGRkZio2NLfM1sbGxbv0laf369eX2xw88GWtJmjVrlqZMmaL09HR169atJkq1tMqOc7t27fTpp59q9+7dru33v/+9brzxRu3evVtRUVE1Wb5lePL93LNnTx04cMAVFCVp3759ioiIIMSUw5NxPn/+fKmwcik8Gj47ucp45WdhtS0jtrhly5YZu91uFi1aZL744gszevRoU79+fZOTk2OMMWbYsGHmySefdPX/4IMPjL+/v3nuuefMl19+aSZOnMjj1xVU2bGeMWOGCQwMNCtWrDDHjx93bWfOnPHWl2AJlR3nn+KppYqp7DgfPXrU1KtXzzz00ENm79695p133jGNGzc2U6dO9daXYAmVHeeJEyeaevXqmTfeeMMcPHjQvPvuu6Zly5bmjjvu8NaXYAlnzpwxu3btMrt27TKSzOzZs82uXbvMkSNHjDHGPPnkk2bYsGGu/pcev3788cfNl19+aVJTU3n82pteeukl07RpUxMYGGiuu+46s3XrVtexPn36mMTERLf+b775pmnTpo0JDAw0HTt2NKtXr67hiq2rMmPdrFkzI6nUNnHixJov3GIq+z39YwSZiqvsOH/44Yeme/fuxm63mxYtWphp06aZixcv1nDV1lOZcS4qKjLPPPOMadmypXE4HCYqKso8+OCD5tSpUzVfuIVs2rSpzH9vL41tYmKi6dOnT6nXdO3a1QQGBpoWLVqYtLS0aq3RZgxzagAAwJpYIwMAACyLIAMAACyLIAMAACyLIAMAACyLIAMAACyLIAMAACyLIAMAACyLIAMAACyLIAMAACyLIAMAACyLIAMAACzr/wGOpwQ8TSKXlwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# convert results to a dataframe \n",
    "df = t.iterator().to_dataframe()\n",
    "\n",
    "# Use the pandas query method to show features from protein_coding genes with non-zero mappability...\n",
    "display(SEP, \"Dataframe representation, filtered for gene_type=='protein_coding' & mappability>0\", \n",
    "        df.query('gene_type==\"protein_coding\" & mappability>0').head(3),SEP)\n",
    "\n",
    "# ...and plot the overall mappability distribution...\n",
    "_ = df['mappability'].plot(kind='hist', title='mappability histogram')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1a4fcc3-78ee-4707-837e-01b123a61f8f",
   "metadata": {},
   "source": [
    "# Summary\n",
    "\n",
    "In this notebook, we demonstrated the basic elements of *rnalib*'s API: \n",
    "\n",
    "* how to instantiate a transcriptome from a GFF file,\n",
    "* how to access and query its features and\n",
    "* how to use *genomic iterators* to integrate and annotate it with other genomic datasets. \n",
    "\n",
    "Note that you can also add arbitrary custom features (e.g., genomic intervals representing regulatory regions) to a transcriptome, \n",
    "save and load annotations to/from disk, export them to GFF3 files or convert them to a bioframe compatible dataframe and much more. \n",
    "See the [AdvancedUsage notebook](AdvancedUsage.ipynb) for further details.\n",
    "\n",
    "This concludes our introduction to the basic *rnalib* API. \n",
    "Please see our [AdvancedUsage notebook](AdvancedUsage.ipynb) for more advanced usage examples or \n",
    "refer to our [other notebooks](https://github.com/popitsch/rnalib/tree/main/notebooks) for realistic application scenarios \n",
    "of *rnalib* and to the [test classes](https://github.com/popitsch/rnalib/tree/main/tests) for more usage examples."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "rnalib",
   "language": "python",
   "name": "rnalib"
  },
  "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.12.1"
  },
  "toc-autonumbering": false,
  "toc-showcode": false,
  "toc-showmarkdowntxt": false,
  "toc-showtags": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}