02_preprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "---\n",
    "output-file: preprocessing.html\n",
    "title: preprocessing\n",
    "\n",
    "---\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# | default_exp preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# | exports\n",
    "\n",
    "import pysam\n",
    "from Bio.Seq import Seq\n",
    "import gzip\n",
    "from LAFITE.utils import Vividict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# | exports\n",
    "\n",
    "def split_bed_line(entry, extends=False):\n",
    "    \"\"\"\n",
    "    split bed line\n",
    "    \"\"\"\n",
    "    cells = entry.strip('\\n').split(\"\\t\")\n",
    "    chrom = start = end = name = score = strand = thick_start = thick_end = item_rgb = block_count = block_sizes = block_starts = None\n",
    "    if len(cells) >= 6:\n",
    "        chrom = cells[0]\n",
    "        start = int(cells[1])\n",
    "        end = int(cells[2])\n",
    "        name = cells[3]\n",
    "        score = cells[4]\n",
    "        strand = cells[5]\n",
    "    if len(cells) >= 12:\n",
    "        thick_start = int(cells[6])\n",
    "        thick_end = int(cells[7])\n",
    "        item_rgb = cells[8]\n",
    "        block_count = int(cells[9])\n",
    "        block_sizes = [int(i) for i in cells[10].rstrip(',').split(',')]\n",
    "        block_starts = [int(i) for i in cells[11].rstrip(',').split(',')]\n",
    "\n",
    "    return chrom, start, end, name, score, strand, thick_start, thick_end, item_rgb, block_count, block_sizes, block_starts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# | exports\n",
    "\n",
    "def bed_block_to_splicing(start, block_count, block_starts, block_sizes):\n",
    "    \"\"\"convert bed block to isoform splicing structure\n",
    "    \"\"\"\n",
    "    read_splicing = []\n",
    "    if block_count > 1:\n",
    "        for i in range(block_count-1):\n",
    "            left_sj = start + block_starts[i] + block_sizes[i]\n",
    "            right_sj = start + block_starts[i+1] + 1\n",
    "            read_splicing.extend([left_sj, right_sj])\n",
    "    return read_splicing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# | exports\n",
    "\n",
    "def read_grouping(bed_file, fasta):\n",
    "    \"\"\"\n",
    "    return junctions with only one long read support\n",
    "    \"\"\"\n",
    "    junction_dict = Vividict()\n",
    "    processed_read = Vividict()\n",
    "    genome = pysam.FastaFile(fasta)\n",
    "    with open(bed_file) as f:\n",
    "        for read in f:\n",
    "            chrom, start, end, name, score, strand, thick_start, thick_end, item_rgb, block_count, block_sizes, block_starts = split_bed_line(\n",
    "                read)  # read the each nanopore read\n",
    "            full_block = [start+1]\n",
    "            if block_count > 1:\n",
    "                for i in range(block_count-1):\n",
    "                    left_sj = start + block_starts[i] + block_sizes[i]\n",
    "                    right_sj = start + block_starts[i+1] + 1\n",
    "                    full_block.extend([left_sj, right_sj])\n",
    "                    tmp_sj = (chrom, strand, left_sj, right_sj)\n",
    "\n",
    "                    if tmp_sj not in junction_dict[(chrom, strand)]:\n",
    "                        if strand == '+':\n",
    "                            intron_motif = f'{genome.fetch(chrom, left_sj, left_sj+2)}-{genome.fetch(chrom, right_sj-3, right_sj-1)}'\n",
    "                        else:\n",
    "                            intron_motif = f'{Seq(genome.fetch(chrom, right_sj-3, right_sj-1)).reverse_complement()}-{Seq(genome.fetch(chrom, left_sj, left_sj+2)).reverse_complement()}'\n",
    "                        if intron_motif in ['GT-AG', 'CT-AC', 'GC-AG', 'CT-GC', 'AT-AC', 'GT-AT']:\n",
    "                            junction_dict[(chrom, strand)][tmp_sj] = [\n",
    "                                1, intron_motif, 'canonical']\n",
    "                        else:\n",
    "                            junction_dict[(chrom, strand)][tmp_sj] = [\n",
    "                                1, intron_motif, 'non-canonical']\n",
    "                    else:\n",
    "                        junction_dict[(chrom, strand)][tmp_sj][0] += 1\n",
    "            full_block.append(end)\n",
    "            name = f'{start+1}_{name}'\n",
    "            processed_read[(chrom, strand)][name] = full_block\n",
    "\n",
    "    return junction_dict, processed_read"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# | exports\n",
    "\n",
    "def polya_signal_import(polyadenylation_event):\n",
    "    \"\"\"\n",
    "    read the polyA information from nanopolish-polya mode\n",
    "    \"\"\"\n",
    "    polya_reads = {}\n",
    "    if polyadenylation_event:\n",
    "        if polyadenylation_event.endswith('.gz'):\n",
    "            OPEN = gzip.open\n",
    "            mode = 'rb'\n",
    "        else:\n",
    "            OPEN = open\n",
    "            mode = 'r'\n",
    "        with OPEN(polyadenylation_event,mode=mode) as f:\n",
    "            next(f)\n",
    "            for line in f:\n",
    "                line = line.rstrip().split(\"\\t\")\n",
    "                if len(line) != 10:\n",
    "                    raise ValueError(\n",
    "                        \"Fatal: input polyA results was not from nanopolish!\")\n",
    "                elif line[9] == 'PASS':\n",
    "                    polya_reads[f'{int(line[2])+1}_{line[0]}'] = True\n",
    "    return polya_reads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# | exports\n",
    "\n",
    "class PolyAFinder():\n",
    "    \"\"\"estimate read polyadenylation event by checking the polya motifs from sequence\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, processed_read, fasta, polyA_motif_file, updis=40,  downdis=10):\n",
    "        self.processed_read = processed_read\n",
    "        self.fasta = fasta\n",
    "        self.polyA_motif_file = polyA_motif_file\n",
    "        self.updis = updis\n",
    "        self.downdis = downdis\n",
    "\n",
    "    def polyA_motif_import(self):\n",
    "        \"\"\"import given polyA motif\n",
    "        \"\"\"\n",
    "        polyA_motif_list = []\n",
    "        for line in open(self.polyA_motif_file):\n",
    "            x = line.strip().upper().replace('U', 'A')\n",
    "            polyA_motif_list.append(x)\n",
    "\n",
    "        return polyA_motif_list\n",
    "\n",
    "    def find_polyA_motif(self, seq, polyA_motif_list):\n",
    "        \"\"\"check the occurrence of given polyA motif in target sequence\n",
    "        \"\"\"\n",
    "        for motif in polyA_motif_list:\n",
    "            i = seq.find(motif)\n",
    "            if i >= 0:\n",
    "                return motif\n",
    "        return None\n",
    "\n",
    "    def polya_estimation(self):\n",
    "        \"\"\"estimate read polyadenylation event\n",
    "        \"\"\"\n",
    "        genome = pysam.FastaFile(self.fasta)\n",
    "        read_polya_dict = {}\n",
    "        polyA_motif_list = self.polyA_motif_import()\n",
    "\n",
    "        for (chrom, strand) in self.processed_read:\n",
    "            for read_id, full_block in self.processed_read[(chrom, strand)].items():\n",
    "\n",
    "                if strand == '-':\n",
    "                    res = full_block[0]\n",
    "                else:\n",
    "                    strand = '+'\n",
    "                    res = full_block[-1]\n",
    "                try:\n",
    "                    if strand == '-':\n",
    "                        seq = f'{Seq(genome.fetch(chrom, res - self.downdis, res + self.updis)).reverse_complement()}'\n",
    "                    else:\n",
    "                        seq = f'{genome.fetch(chrom, res - self.updis, res + self.downdis)}'\n",
    "                    polyA_motif = self.find_polyA_motif(seq, polyA_motif_list)\n",
    "\n",
    "                    if polyA_motif:\n",
    "                        read_polya_dict[read_id] = True\n",
    "                except:\n",
    "                    pass\n",
    "        return read_polya_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "python3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}