docs/source/parsing_gff3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parsing GFF3 files\n",
    "\n",
    "Without specification, the default GFF3 parsing function will be used. This function will call `gffutils` and instruct it to build a new database in memory, with the flag `merge_strategy=\"create_unique\"`. These can be adjusted by modfying the arguments to `GffutilsParseArgs`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from inscripta.biocantor.io.gff3.parser import parse_standard_gff3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gff3 = \"tests/data/INSC1006_chrI.gff3\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# GFF3 parsing requires that the file be passed directly instead as an open handle\n",
    "annotation_records = list(parse_standard_gff3(gff3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After parsing, there will be one `AnnotationRecord` built for every sequence in the GFF3 file. This dataclass is a wrapper that contains two objects, a `SeqRecord` object and a `AnnotationCollectionModel`. Because this parser is not given any FASTA information, the `SeqRecord` is null. The `AnnotationCollectionModel` objects are `marshmallow_dataclass` objects, and so can be dumped to and loaded directly from JSON."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OrderedDict([('feature_collections', None),\n",
       "             ('genes',\n",
       "              [OrderedDict([('transcripts',\n",
       "                             [OrderedDict([('exon_starts', [16174]),\n",
       "                                           ('exon_ends', [18079]),\n",
       "                                           ('strand', 'MINUS'),\n",
       "                                           ('cds_starts', None),\n",
       "                                           ('cds_ends', None),\n",
       "                                           ('cds_frames', None),\n",
       "                                           ('qualifiers',\n",
       "                                            {'ncrna_class': ['other'],\n",
       "                                             'note': ['CAT transcript id: T0000001; CAT alignment id: IsoSeq-PB.2586.1; CAT novel prediction: IsoSeq']}),\n",
       "                                           ('is_primary_tx', False),\n",
       "                                           ('transcript_id', 'GI526_G0000001'),\n",
       "                                           ('protein_id', None),\n",
       "                                           ('product', None),\n",
       "                                           ('transcript_symbol',\n",
       "                                            'GI526_G0000001'),\n",
       "                                           ('transcript_type', 'ncRNA'),\n",
       "                                           ('sequence_name', 'CM021111.1'),\n",
       "                                           ('sequence_guid', None),\n",
       "                                           ('transcript_interval_guid', None),\n",
       "                                           ('transcript_guid', None)])]),\n",
       "                            ('gene_id',\n",
       "                             '8ad3f444-384e-35e0-e560-aef88bd2863f'),\n",
       "                            ('gene_symbol', None),\n",
       "                            ('gene_type', 'ncRNA'),\n",
       "                            ('locus_tag', 'GI526_G0000001'),\n",
       "                            ('qualifiers', None),\n",
       "                            ('sequence_name', 'CM021111.1'),\n",
       "                            ('sequence_guid', None),\n",
       "                            ('gene_guid', None)]),\n",
       "               OrderedDict([('transcripts',\n",
       "                             [OrderedDict([('exon_starts', [37461]),\n",
       "                                           ('exon_ends', [39103]),\n",
       "                                           ('strand', 'PLUS'),\n",
       "                                           ('cds_starts', [37637]),\n",
       "                                           ('cds_ends', [39011]),\n",
       "                                           ('cds_frames', ['ZERO']),\n",
       "                                           ('qualifiers',\n",
       "                                            {'codon_start': ['1'],\n",
       "                                             'gene': ['GDH3'],\n",
       "                                             'note': ['CAT transcript id: T0000002; CAT alignment id: NM_001178204.1-0; CAT source transcript id: NM_001178204.1; CAT source transcript biotype: protein_coding'],\n",
       "                                             'translation': ['MTSEPEFQQAYDEIVSSVEDSKIFEKFPQYKKVLPIVSVPERIIQFRVTWENDNGEQEVAQGYRVQFNSAKGPYKGGLRFHPSVNLSILKFLGFEQIFKNALTGLDMGGGKGGLCVDLKGKSDNEIRRICYAFMRELSRHIGKDTDVPAGDIGVGGREIGYLFGAYRSYKNSWEGVLTGKGLNWGGSLIRPEATGFGLVYYTQAMIDYATNGKESFEGKRVTISGSGNVAQYAALKVIELGGIVVSLSDSKGCIISETGITSEQIHDIASAKIRFKSLEEIVDEYSTFSESKMKYVAGARPWTHVSNVDIALPCATQNEVSGDEAKALVASGVKFVAEGANMGSTPEAISVFETARSTATNAKDAVWFGPPKAANLGGVAVSGLEMAQNSQKVTWTAERVDQELKKIMINCFNDCIQAAQEYSTEKNTNTLPSLVKGANIASFVMVADAMLDQGDVF']}),\n",
       "                                           ('is_primary_tx', False),\n",
       "                                           ('transcript_id', 'GI526_G0000002'),\n",
       "                                           ('protein_id', 'KAF1903245.1'),\n",
       "                                           ('product', 'GDH3 isoform 1'),\n",
       "                                           ('transcript_symbol', 'GDH3'),\n",
       "                                           ('transcript_type',\n",
       "                                            'protein_coding'),\n",
       "                                           ('sequence_name', 'CM021111.1'),\n",
       "                                           ('sequence_guid', None),\n",
       "                                           ('transcript_interval_guid', None),\n",
       "                                           ('transcript_guid', None)])]),\n",
       "                            ('gene_id',\n",
       "                             'a1b669f1-57f6-ae9b-8f4f-a27a6e84d15a'),\n",
       "                            ('gene_symbol', 'GDH3'),\n",
       "                            ('gene_type', 'protein_coding'),\n",
       "                            ('locus_tag', 'GI526_G0000002'),\n",
       "                            ('qualifiers', {'gene': ['GDH3']}),\n",
       "                            ('sequence_name', 'CM021111.1'),\n",
       "                            ('sequence_guid', None),\n",
       "                            ('gene_guid', None)]),\n",
       "               OrderedDict([('transcripts',\n",
       "                             [OrderedDict([('exon_starts', [39518]),\n",
       "                                           ('exon_ends', [40772]),\n",
       "                                           ('strand', 'PLUS'),\n",
       "                                           ('cds_starts', [39518]),\n",
       "                                           ('cds_ends', [40772]),\n",
       "                                           ('cds_frames', ['ZERO']),\n",
       "                                           ('qualifiers',\n",
       "                                            {'codon_start': ['1'],\n",
       "                                             'gene': ['BDH2'],\n",
       "                                             'translation': ['MRALAYFGKGNIRFTNHLKEPHIVAPDELVIDIAWCGICGTDLHEYTDGPIFFPEDGHTHEISHNPLPQAMGHEMAGTVLEVGPSVKNLKVGDKVVVEPTGTCRDRYRWPLSPKVDKEWCAACKKGYYNICSYLGLCGAGVQSGGFAEGVVMNESHCYKVPDFVPLDVAALIQPLAVCWHAIRVCEFKAGSTALIIGAGPIGLGTILALNAAGCKDIVVSEPAKVRRELAEKMGARVYDPTAHAAKESIDYLRSIADGGDGFDYTFDCSGLEVTLNAAIQCLTFRGTAVNLAMWGHHKIQFSPMDITLHERKYTGSMCYTHHDFETVIEALEEGRIDIDRARHMITGRVNIEDGLDGAIMKLINEKESTIKIILTPNNHGELNREADNEKKEISELSSRKDQERLRESINEAKLRHT']}),\n",
       "                                           ('is_primary_tx', False),\n",
       "                                           ('transcript_id', 'GI526_G0000003'),\n",
       "                                           ('protein_id', 'KAF1903246.1'),\n",
       "                                           ('product', 'BDH2 isoform 1'),\n",
       "                                           ('transcript_symbol', 'BDH2'),\n",
       "                                           ('transcript_type',\n",
       "                                            'protein_coding'),\n",
       "                                           ('sequence_name', 'CM021111.1'),\n",
       "                                           ('sequence_guid', None),\n",
       "                                           ('transcript_interval_guid', None),\n",
       "                                           ('transcript_guid', None)])]),\n",
       "                            ('gene_id',\n",
       "                             '4967ade5-6d91-faeb-79ed-e57093e4e5f2'),\n",
       "                            ('gene_symbol', 'BDH2'),\n",
       "                            ('gene_type', 'protein_coding'),\n",
       "                            ('locus_tag', 'GI526_G0000003'),\n",
       "                            ('qualifiers', {'gene': ['BDH2']}),\n",
       "                            ('sequence_name', 'CM021111.1'),\n",
       "                            ('sequence_guid', None),\n",
       "                            ('gene_guid', None)]),\n",
       "               OrderedDict([('transcripts',\n",
       "                             [OrderedDict([('exon_starts', [41085]),\n",
       "                                           ('exon_ends', [42503]),\n",
       "                                           ('strand', 'PLUS'),\n",
       "                                           ('cds_starts', None),\n",
       "                                           ('cds_ends', None),\n",
       "                                           ('cds_frames', None),\n",
       "                                           ('qualifiers',\n",
       "                                            {'gene': ['BDH1'],\n",
       "                                             'note': ['CAT transcript id: T0000004; CAT alignment id: NM_001178202.2-0; CAT source transcript id: NM_001178202.2; CAT source transcript biotype: protein_coding']}),\n",
       "                                           ('is_primary_tx', False),\n",
       "                                           ('transcript_id', 'GI526_G0000004'),\n",
       "                                           ('protein_id', None),\n",
       "                                           ('product', None),\n",
       "                                           ('transcript_symbol', 'BDH1'),\n",
       "                                           ('transcript_type',\n",
       "                                            'protein_coding'),\n",
       "                                           ('sequence_name', 'CM021111.1'),\n",
       "                                           ('sequence_guid', None),\n",
       "                                           ('transcript_interval_guid', None),\n",
       "                                           ('transcript_guid', None)])]),\n",
       "                            ('gene_id',\n",
       "                             '278a932c-a0a7-e31b-3156-5860ca4a4021'),\n",
       "                            ('gene_symbol', 'BDH1'),\n",
       "                            ('gene_type', 'protein_coding'),\n",
       "                            ('locus_tag', 'GI526_G0000004'),\n",
       "                            ('qualifiers', {'gene': ['BDH1']}),\n",
       "                            ('sequence_name', 'CM021111.1'),\n",
       "                            ('sequence_guid', None),\n",
       "                            ('gene_guid', None)]),\n",
       "               OrderedDict([('transcripts',\n",
       "                             [OrderedDict([('exon_starts', [42579]),\n",
       "                                           ('exon_ends', [43218]),\n",
       "                                           ('strand', 'PLUS'),\n",
       "                                           ('cds_starts', None),\n",
       "                                           ('cds_ends', None),\n",
       "                                           ('cds_frames', None),\n",
       "                                           ('qualifiers', {'gene': ['ECM1']}),\n",
       "                                           ('is_primary_tx', False),\n",
       "                                           ('transcript_id', 'GI526_G0000005'),\n",
       "                                           ('protein_id', None),\n",
       "                                           ('product', None),\n",
       "                                           ('transcript_symbol', 'ECM1'),\n",
       "                                           ('transcript_type', 'ncRNA'),\n",
       "                                           ('sequence_name', 'CM021111.1'),\n",
       "                                           ('sequence_guid', None),\n",
       "                                           ('transcript_interval_guid', None),\n",
       "                                           ('transcript_guid', None)])]),\n",
       "                            ('gene_id',\n",
       "                             'e8b37537-588b-43a2-eb26-c88883f06014'),\n",
       "                            ('gene_symbol', 'ECM1'),\n",
       "                            ('gene_type', 'ncRNA'),\n",
       "                            ('locus_tag', 'GI526_G0000005'),\n",
       "                            ('qualifiers', {'gene': ['ECM1']}),\n",
       "                            ('sequence_name', 'CM021111.1'),\n",
       "                            ('sequence_guid', None),\n",
       "                            ('gene_guid', None)])]),\n",
       "             ('name', None),\n",
       "             ('id', None),\n",
       "             ('sequence_name', 'CM021111.1'),\n",
       "             ('sequence_guid', None),\n",
       "             ('sequence_path', None),\n",
       "             ('qualifiers', None),\n",
       "             ('start', None),\n",
       "             ('end', None),\n",
       "             ('completely_within', None),\n",
       "             ('parent_or_seq_chunk_parent', None)])"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "annotation_collection_model = annotation_records[0].annotation\n",
    "annotation_collection_model.Schema().dump(annotation_collection_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Converting models to BioCantor data structures\n",
    "After loading an `AnnotationCollectionModel`, this object can be directly converted in to an `AnnotationCollection` with sequence information.\n",
    "\n",
    "`AnnotationCollection` objects are the core data structure, and contain a set of genes and features as children."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "annotation_collection = annotation_collection_model.to_annotation_collection()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GeneInterval(identifiers={'8ad3f444-384e-35e0-e560-aef88bd2863f', 'GI526_G0000001'}, Intervals:TranscriptInterval((16174-18079:-), cds=[None], symbol=GI526_G0000001))\n",
      "GeneInterval(identifiers={'GI526_G0000002', 'a1b669f1-57f6-ae9b-8f4f-a27a6e84d15a', 'GDH3'}, Intervals:TranscriptInterval((37461-39103:+), cds=[CDS((37637-39011:+), (CDSFrame.ZERO)], symbol=GDH3))\n",
      "GeneInterval(identifiers={'4967ade5-6d91-faeb-79ed-e57093e4e5f2', 'BDH2', 'GI526_G0000003'}, Intervals:TranscriptInterval((39518-40772:+), cds=[CDS((39518-40772:+), (CDSFrame.ZERO)], symbol=BDH2))\n",
      "GeneInterval(identifiers={'BDH1', 'GI526_G0000004', '278a932c-a0a7-e31b-3156-5860ca4a4021'}, Intervals:TranscriptInterval((41085-42503:+), cds=[None], symbol=BDH1))\n",
      "GeneInterval(identifiers={'ECM1', 'e8b37537-588b-43a2-eb26-c88883f06014', 'GI526_G0000005'}, Intervals:TranscriptInterval((42579-43218:+), cds=[None], symbol=ECM1))\n"
     ]
    }
   ],
   "source": [
    "# this example dataset has 4 genes and 0 features\n",
    "for child in annotation_collection:\n",
    "    print(child)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene1 = annotation_collection.genes[0]\n",
    "tx1 = gene1.transcripts[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18078"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# convert mRNA coordinates to genomic coordinates\n",
    "tx1.transcript_pos_to_sequence(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "NoncodingTranscriptError",
     "evalue": "No CDS positions on non-coding transcript",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNoncodingTranscriptError\u001b[0m                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-12-0276456c0756>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# NoncodingTranscriptError is raised when trying to convert CDS coordinates on a non-coding transcript\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mtx1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcds_pos_to_sequence\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/Users/ian.fiddes/repos/inscripta_content_store/dependencies/biocantor/inscripta/biocantor/gene/transcript.py\u001b[0m in \u001b[0;36mcds_pos_to_sequence\u001b[0;34m(self, pos)\u001b[0m\n\u001b[1;32m    549\u001b[0m         \u001b[0;34m\"\"\"Converts a relative position along the CDS to sequence coordinate.\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    550\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_coding\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 551\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mNoncodingTranscriptError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"No CDS positions on non-coding transcript\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    552\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcds\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mchromosome_location\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrelative_to_parent_pos\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpos\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    553\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNoncodingTranscriptError\u001b[0m: No CDS positions on non-coding transcript"
     ]
    }
   ],
   "source": [
    "# NoncodingTranscriptError is raised when trying to convert CDS coordinates on a non-coding transcript\n",
    "tx1.cds_pos_to_sequence(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Primary transcripts\n",
    "\n",
    "It is often useful to have an understanding of what isoform of a gene is the 'most important'. An input dataset can provide this information based on the parser implementation used. If this information is not provided, then this value is inferred by the simple heuristic of:\n",
    "\n",
    "1. Longest CDS isoform.\n",
    "2. Longest isoform (if no coding isoforms)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gene1.get_primary_transcript() == tx1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Incorporating sequence information\n",
    "\n",
    "The above parsing call treated the input GFF3 as a standard GFF3 without sequence. However, this file actually has a FASTA suffix, and so the sequence information can be loaded either from that file, or from a separate FASTA file.\n",
    "\n",
    "When parsers are used that return sequence information, the returned `ParsedAnnotationRecord` now has a non-null `seqrecord` member. This object can then be used to instantiate an `AnnotationCollection` with sequence information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Sequence;\n",
       "  Alphabet=NT_EXTENDED_GAPPED;\n",
       "  Length=50040;\n",
       "  Parent=None;\n",
       "  Type=None>"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from inscripta.biocantor.io.gff3.parser import parse_gff3_embedded_fasta\n",
    "\n",
    "parsed_with_sequence = list(parse_gff3_embedded_fasta(gff3))\n",
    "annotation_collection_with_sequence = parsed_with_sequence[0].to_annotation_collection()\n",
    "annotation_collection_with_sequence.get_reference_sequence()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GI526_G0000001: mRNA: GCGGCGCTCT...\n",
      "GDH3: mRNA: AAACAGTTAA... protein: MTSEPEFQQA...\n",
      "BDH2: mRNA: ATGAGAGCCT... protein: MRALAYFGKG...\n",
      "BDH1: mRNA: GGGGCAGATA...\n",
      "ECM1: mRNA: ATGTGGGAAC...\n"
     ]
    }
   ],
   "source": [
    "for gene in annotation_collection_with_sequence:\n",
    "    for tx in gene.transcripts:\n",
    "        if tx.is_coding:\n",
    "            print(f\"{tx.transcript_symbol}: mRNA: {tx.get_spliced_sequence()[:10]}... protein: {tx.get_protein_sequence()[:10]}...\")\n",
    "        else:\n",
    "            print(f\"{tx.transcript_symbol}: mRNA: {tx.get_spliced_sequence()[:10]}...\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sequence information can also come from a separate FASTA file. In this case, a different function call is used, and it gets provided the FASTA directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from inscripta.biocantor.io.gff3.parser import parse_gff3_fasta\n",
    "\n",
    "fasta = \"tests/data/INSC1006_chrI.fasta\"\n",
    "parsed_with_sequence_from_fasta = list(parse_gff3_embedded_fasta(gff3))\n",
    "annotation_collection_with_sequence_from_fasta = parsed_with_sequence_from_fasta[0].to_annotation_collection()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Sequence;\n",
       "  Alphabet=NT_EXTENDED_GAPPED;\n",
       "  Length=50040;\n",
       "  Parent=None;\n",
       "  Type=None>"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "annotation_collection_with_sequence_from_fasta.get_reference_sequence()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Querying the collection\n",
    "\n",
    "`AnnotationCollections` have the ability to be subsetted. These range queries can be performed in two modes, controlled by the flag `completely_within`. When `completely_within = True`, the positions in the query are exact bounds. When `completely_within = False`, any constituent object that overlaps the range query will be retained.\n",
    "\n",
    "`start` and `end` are not required to be set, and are inferred to be `0` and `len(sequence)` respectively if not used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['GDH3', 'BDH2', 'BDH1', 'ECM1']\n"
     ]
    }
   ],
   "source": [
    "# remove GI526_G0000001 by moving the start position to within its bounds, when strict boundaries are required\n",
    "subset1 = annotation_collection.query_by_position(start=16175, completely_within=True)\n",
    "print([x.gene_symbol for x in subset1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['BDH2', 'BDH1']\n"
     ]
    }
   ],
   "source": [
    "# select BDH1 and BDH2\n",
    "subset2 = annotation_collection.query_by_position(start=40000, end=42000, completely_within=False)\n",
    "print([x.gene_symbol for x in subset2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modifying the parser\n",
    "\n",
    "Each of the GFF3 wrapper parser functions, `parse_standard_gff3`, `parse_gff3_embedded_fasta`, and `parse_gff3_fasta`, all have arguments that allow you to modify how the GFF3 file is interpreted. These include:\n",
    "\n",
    "1. `gffutil_parse_args=GffutilsParseArgs()`: This dataclass contains arguments that are passed to gffutils.\n",
    "2. `parse_func`: The default implementation of this function can be written bespoke to your annotation data.\n",
    "3. `gffutil_transform_func`: This function is passed straight to the same argument in `gffutils`.\n",
    "4. `db_fn`: Change this value to retain the `gffutils` sqlite database on disk."
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}