mirpipe_pipeline.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Open source RNA-Seq pipeline for identification of novel-mirs and their gene regulatory networks<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#MiRPipe-Flowchart\" data-toc-modified-id=\"MiRPipe-Flowchart-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>MiRPipe Flowchart</a></span></li><li><span><a href=\"#FastQ-Files\" data-toc-modified-id=\"FastQ-Files-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>FastQ Files</a></span></li><li><span><a href=\"#Pre-processing\" data-toc-modified-id=\"Pre-processing-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Pre-processing</a></span><ul class=\"toc-item\"><li><span><a href=\"#Quality-Check-of-Fastq-Files\" data-toc-modified-id=\"Quality-Check-of-Fastq-Files-3.1\"><span class=\"toc-item-num\">3.1&nbsp;&nbsp;</span>Quality Check of Fastq Files</a></span></li><li><span><a href=\"#Adaptor-Timming-and-Fastq-Splitting\" data-toc-modified-id=\"Adaptor-Timming-and-Fastq-Splitting-3.2\"><span class=\"toc-item-num\">3.2&nbsp;&nbsp;</span>Adaptor Timming and Fastq Splitting</a></span></li></ul></li><li><span><a href=\"#Sequence-Alignment\" data-toc-modified-id=\"Sequence-Alignment-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Sequence Alignment</a></span><ul class=\"toc-item\"><li><span><a href=\"#Downloading-the-Mirdeep*-aligner-and-mirbase\" data-toc-modified-id=\"Downloading-the-Mirdeep*-aligner-and-mirbase-4.1\"><span class=\"toc-item-num\">4.1&nbsp;&nbsp;</span><strong>Downloading the Mirdeep* aligner and mirbase</strong></a></span></li><li><span><a href=\"#piRNA-Pipeline\" data-toc-modified-id=\"piRNA-Pipeline-4.2\"><span class=\"toc-item-num\">4.2&nbsp;&nbsp;</span>piRNA Pipeline</a></span></li><li><span><a href=\"#miRNA-sequence-Alignment\" data-toc-modified-id=\"miRNA-sequence-Alignment-4.3\"><span class=\"toc-item-num\">4.3&nbsp;&nbsp;</span>miRNA sequence Alignment</a></span></li></ul></li><li><span><a href=\"#Post-Alignment-Analysis\" data-toc-modified-id=\"Post-Alignment-Analysis-5\"><span class=\"toc-item-num\">5&nbsp;&nbsp;</span>Post-Alignment Analysis</a></span><ul class=\"toc-item\"><li><span><a href=\"#Processing-of-raw-counts-from-Mirdeep*-results\" data-toc-modified-id=\"Processing-of-raw-counts-from-Mirdeep*-results-5.1\"><span class=\"toc-item-num\">5.1&nbsp;&nbsp;</span><strong>Processing of raw counts from Mirdeep* results</strong></a></span></li></ul></li><li><span><a href=\"#Blast-Search\" data-toc-modified-id=\"Blast-Search-6\"><span class=\"toc-item-num\">6&nbsp;&nbsp;</span>Blast Search</a></span><ul class=\"toc-item\"><li><span><a href=\"#Search-for-Novel-miRNA-sequence-in-DashR-Database\" data-toc-modified-id=\"Search-for-Novel-miRNA-sequence-in-DashR-Database-6.1\"><span class=\"toc-item-num\">6.1&nbsp;&nbsp;</span>Search for Novel miRNA sequence in DashR Database</a></span></li><li><span><a href=\"#DashR-Results-post-processing\" data-toc-modified-id=\"DashR-Results-post-processing-6.2\"><span class=\"toc-item-num\">6.2&nbsp;&nbsp;</span>DashR Results post-processing</a></span></li></ul></li><li><span><a href=\"#Seed-based-Sequence-Clustering\" data-toc-modified-id=\"Seed-based-Sequence-Clustering-7\"><span class=\"toc-item-num\">7&nbsp;&nbsp;</span>Seed based Sequence Clustering</a></span><ul class=\"toc-item\"><li><span><a href=\"#Dictionary-Preparation\" data-toc-modified-id=\"Dictionary-Preparation-7.1\"><span class=\"toc-item-num\">7.1&nbsp;&nbsp;</span><strong>Dictionary Preparation</strong></a></span></li><li><span><a href=\"#CD-HIT\" data-toc-modified-id=\"CD-HIT-7.2\"><span class=\"toc-item-num\">7.2&nbsp;&nbsp;</span>CD-HIT</a></span></li><li><span><a href=\"#Seed-Clustering-Algorithm\" data-toc-modified-id=\"Seed-Clustering-Algorithm-7.3\"><span class=\"toc-item-num\">7.3&nbsp;&nbsp;</span>Seed Clustering Algorithm</a></span></li></ul></li><li><span><a href=\"#miRNA-Differential-expression-Analysis-using-DESeq2\" data-toc-modified-id=\"miRNA-Differential-expression-Analysis-using-DESeq2-8\"><span class=\"toc-item-num\">8&nbsp;&nbsp;</span>miRNA Differential expression Analysis using DESeq2</a></span></li><li><span><a href=\"#Finding-significantly-altered-differentially-expressed-miRNAs\" data-toc-modified-id=\"Finding-significantly-altered-differentially-expressed-miRNAs-9\"><span class=\"toc-item-num\">9&nbsp;&nbsp;</span>Finding significantly altered differentially expressed miRNAs</a></span></li><li><span><a href=\"#Naming-of-Novel-miRNAs\" data-toc-modified-id=\"Naming-of-Novel-miRNAs-10\"><span class=\"toc-item-num\">10&nbsp;&nbsp;</span>Naming of Novel miRNAs</a></span><ul class=\"toc-item\"><li><span><a href=\"#Rfam-Database-Integration\" data-toc-modified-id=\"Rfam-Database-Integration-10.1\"><span class=\"toc-item-num\">10.1&nbsp;&nbsp;</span>Rfam Database Integration</a></span></li><li><span><a href=\"#Renaming-of-novel-miRNAs\" data-toc-modified-id=\"Renaming-of-novel-miRNAs-10.2\"><span class=\"toc-item-num\">10.2&nbsp;&nbsp;</span>Renaming of novel miRNAs</a></span></li></ul></li><li><span><a href=\"#piRNA-Counts-Analysis\" data-toc-modified-id=\"piRNA-Counts-Analysis-11\"><span class=\"toc-item-num\">11&nbsp;&nbsp;</span>piRNA Counts Analysis</a></span></li><li><span><a href=\"#piRNA-Differential-expression-Analysis-using-DESeq2\" data-toc-modified-id=\"piRNA-Differential-expression-Analysis-using-DESeq2-12\"><span class=\"toc-item-num\">12&nbsp;&nbsp;</span>piRNA Differential expression Analysis using DESeq2</a></span></li><li><span><a href=\"#Visualizations\" data-toc-modified-id=\"Visualizations-13\"><span class=\"toc-item-num\">13&nbsp;&nbsp;</span>Visualizations</a></span><ul class=\"toc-item\"><li><span><a href=\"#Network-for-miRNA-Gene-Interaction-using-miRNet\" data-toc-modified-id=\"Network-for-miRNA-Gene-Interaction-using-miRNet-13.1\"><span class=\"toc-item-num\">13.1&nbsp;&nbsp;</span>Network for miRNA-Gene Interaction using miRNet</a></span><ul class=\"toc-item\"><li><span><a href=\"#miRNA-Gene-Interaction\" data-toc-modified-id=\"miRNA-Gene-Interaction-13.1.1\"><span class=\"toc-item-num\">13.1.1&nbsp;&nbsp;</span>miRNA-Gene Interaction</a></span></li><li><span><a href=\"#2D-interactive-network-visualization-using-Plotly\" data-toc-modified-id=\"2D-interactive-network-visualization-using-Plotly-13.1.2\"><span class=\"toc-item-num\">13.1.2&nbsp;&nbsp;</span>2D interactive network visualization using Plotly</a></span></li></ul></li></ul></li><li><span><a href=\"#Output\" data-toc-modified-id=\"Output-14\"><span class=\"toc-item-num\">14&nbsp;&nbsp;</span>Output</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MiRPipe Flowchart"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](Figures/miRPipe_Flowchart.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FastQ Files\n",
    "**Loading libraries**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-20T12:51:49.057874Z",
     "start_time": "2019-05-20T12:51:43.698812Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import pickle\n",
    "import pandas as pd\n",
    "import subprocess\n",
    "import glob\n",
    "from splinter import Browser\n",
    "import time\n",
    "import networkx as nx\n",
    "from operator import add\n",
    "from os import path\n",
    "import threading\n",
    "import subprocess\n",
    "import tempfile\n",
    "from rna_tools.Seq import RNASequence\n",
    "from rna_tools.rna_tools_config import RFAM_DB_PATH\n",
    "from threading import Semaphore\n",
    "screenlock = Semaphore(value=1)\n",
    "from tqdm import tqdm\n",
    "from Bio.Seq import Seq\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import shutil"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Declaring ENV Variables**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-20T12:51:49.151397Z",
     "start_time": "2019-05-20T12:51:49.135554Z"
    }
   },
   "outputs": [],
   "source": [
    "# setting env variables:\n",
    "os.environ['HOME_DIR'] = os.getcwd()\n",
    "os.environ['SEQ_DIR'] = os.path.join(os.environ['HOME_DIR'],'data')\n",
    "os.environ['Tools_DIR'] = os.path.join(os.environ['HOME_DIR'],'Tools')    \n",
    "os.environ['REF_DIR'] = os.path.join(os.environ['HOME_DIR'], 'Tools',\n",
    "                                     'MDS_command_line_v38/MDS_command_line')\n",
    "os.mkdir(os.path.join(os.environ['SEQ_DIR'],'output'))\n",
    "print('Please be sure that all the fastq files are present in %s. All the results will be saved in the same folder also.' %os.environ['SEQ_DIR'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ensure that path is correct and fastq files are present in this path\n",
    "For Smooth Functioning of Pipeline, please make sure that only input fastq files are present in the data folder.\n",
    "Delete all the clutter before re-running the pipeline for trouble-free execution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-20T12:51:49.634658Z",
     "start_time": "2019-05-20T12:51:49.156447Z"
    }
   },
   "outputs": [],
   "source": [
    "!echo $HOME_DIR\n",
    "!echo $SEQ_DIR\n",
    "!echo $REF_DIR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pre-processing\n",
    "## Quality Check of Fastq Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "    Following script perform the following task:\n",
    "    1. Perform quality checking of rach samples using FastQC.\n",
    "    2. Prepare the better info-graphic and detailed report of all quality checking\n",
    "       reports from all samples using multiqc\n",
    "'''\n",
    "!Rscript $HOME_DIR/scripts/FastQC.R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adaptor Timming and Fastq Splitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print('The default adaptor sequence used for adaptor trimming is \\\n",
    "TCGTATGCCGTCTTCTGCTTG. If your adaptor sequence is different then please \\\n",
    "enter your data specific adaptor sequence. Please selec appropriate option.')\n",
    "print('1. Default adaptor sequence: TCGTATGCCGTCTTCTGCTTG')\n",
    "print('2. User specific adaptor sequence')\n",
    "user_choice =int(input('Please select your option:'))\n",
    "\n",
    "if user_choice == 1:\n",
    "    adaptor = 'TCGTATGCCGTCTTCTGCTTG'\n",
    "    os.environ['adaptor'] = adaptor\n",
    "elif user_choice == 2:\n",
    "    adaptor = input('Please enter the adaptor sequence')\n",
    "    os.environ['adaptor'] = adaptor    \n",
    "\n",
    "print(f'The choosen adaptor sequence is:{adaptor}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add script for adaptor trimming here using trim_galore\n",
    "!bash scripts/adaptor_trimming.sh\n",
    "\n",
    "# Add script for read length based spliting (into 3 parts) here using bbduk\n",
    "!bash scripts/fastq_split.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sequence Alignment "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Downloading the Mirdeep* aligner and mirbase**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Please enter your choice: ')\n",
    "print(\"Option 1: hg19 based alignment using Mirdeep* with miRBase v19 \")\n",
    "print(\"Option 2: hg19 based alignment using Mirdeep* with miRBase v20\")\n",
    "print(\"Option 3: hg38 based alignment using Mirdeep* with miRBase 21\")\n",
    "print(\"Option 4: hg38 based alignment using Mirdeep* with miRBase 22 (Default Condition)\")\n",
    "print(\"Please enter either 1, 2, 3 or 4\")\n",
    "choice = int(input(\"Please enter your choice:  \"))\n",
    "\n",
    "if choice == 1:\n",
    "    print(\"You chose hg19 and miRBase v19 based alignment using Mirdeep*\")\n",
    "    with open('data/mirdeep.conf',\"w\") as handle:\n",
    "        handle.write(\"Human Genome = hg19\\nmiRBase = 19\")   \n",
    "    print(\"Downloading Mirdeep*....\")\n",
    "    command = \"wget --content-disposition  http://sourceforge.net/projects/mirdeepstar/files/MDS_command_line_v37.zip -O $Tools_DIR/MDS_command_line_v37.zip && \"\n",
    "    command += \"unzip -o $Tools_DIR/MDS_command_line_v37.zip -d $Tools_DIR\"\n",
    "    subprocess.call(command, shell=True)    \n",
    "    print(\"Downloading Mirdeep* is complete. Now downloading hg19 human reference genome....\")\n",
    "    ref_var = 'hg19'    \n",
    "    os.environ['REF_DIR'] = os.path.join(os.environ['HOME_DIR'], 'Tools',\n",
    "                                     'MDS_command_line_v37/MDS_command_line')\n",
    "    print(\"\")\n",
    "    command = \"\"\n",
    "    command += \"mkdir $HOME_DIR/refs/hg19 && \"\n",
    "    command += \"wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz -P $HOME_DIR/refs/hg19/ && \"\n",
    "    command += \"pigz -p 5 -d $HOME_DIR/refs/hg19/hg19.fa.gz && echo Download complete. Now building bowtie indexes && \"\n",
    "    command += \"$HOME_DIR/Tools/bowtie-1.2.3-linux-x86_64/bowtie-build --threads 8 $HOME_DIR/refs/hg19/hg19.fa $HOME_DIR/refs/hg19/hg19\"\n",
    "    subprocess.call(command, shell=True)  \n",
    "    print(\"Human refenrence genome hg19 has been downloaded and index files generation complete.\")\n",
    "    \n",
    "elif choice == 2:\n",
    "    print(\"You chose hg19 and miRBase v20 based alignment using Mirdeep*\")\n",
    "    with open('data/mirdeep.conf',\"w\") as handle:\n",
    "        handle.write(\"Human Genome = hg19\\nmiRBase = 20\")   \n",
    "    command = \"\"\n",
    "    command = \"wget --content-disposition  http://sourceforge.net/projects/mirdeepstar/files/MDS_command_line_v37.zip -O $Tools_DIR/MDS_command_line_v37.zip && \"\n",
    "    command += \"unzip -o $Tools_DIR/MDS_command_line_v37.zip -d $Tools_DIR && \"\n",
    "    command += \"rm $HOME_DIR/Tools/MDS_command_line_v37/MDS_command_line/genome/hg19/miRBase/* && \"\n",
    "    command += \"wget --content-disposition https://www.mirbase.org/ftp/20/hairpin.fa.gz -P $HOME_DIR/Tools/MDS_command_line_v37/MDS_command_line/genome/hg19/miRBase && \"\n",
    "    command += \"pigz -p 5 -d $HOME_DIR/Tools/MDS_command_line_v37/MDS_command_line/genome/hg19/miRBase/hairpin.fa.gz && \"\n",
    "    command += \"wget --content-disposition  https://www.mirbase.org/ftp/20/mature.fa.gz -P $HOME_DIR/Tools/MDS_command_line_v37/MDS_command_line/genome/hg19/miRBase && \"\n",
    "    command += \"pigz -p 5 -d $HOME_DIR/Tools/MDS_command_line_v37/MDS_command_line/genome/hg19/miRBase/mature.fa.gz &&  \"\n",
    "    command += \"wget --content-disposition  https://www.mirbase.org/ftp/20/genomes/hsa.gff3 -P $HOME_DIR/Tools/MDS_command_line_v37/MDS_command_line/genome/hg19/miRBase && \"\n",
    "    command += \"mv $HOME_DIR/Tools/MDS_command_line_v37/MDS_command_line/genome/hg19/miRBase/hsa.gff3 $HOME_DIR/Tools/MDS_command_line_v37/MDS_command_line/genome/hg19/miRBase/knownMiR.gff3\"                \n",
    "    subprocess.call(command, shell=True)   \n",
    "    print(\"Downloading Mirdeep* is complete. Now downloading hg19 human reference genome....\")\n",
    "    ref_var = 'hg19'    \n",
    "    os.environ['REF_DIR'] = os.path.join(os.environ['HOME_DIR'], 'Tools',\n",
    "                                     'MDS_command_line_v37/MDS_command_line')\n",
    "    print(\"\")\n",
    "    command = \"\"\n",
    "    command += \"mkdir $HOME_DIR/refs/hg19 && \"\n",
    "    command += \"wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz -P $HOME_DIR/refs/hg19 && \"\n",
    "    command += \"pigz -p 5 -d $HOME_DIR/refs/hg19/hg19.fa.gz && echo Download complete. Now building bowtie indexes && \"\n",
    "    command += \"$HOME_DIR/Tools/bowtie-1.2.3-linux-x86_64/bowtie-build --threads 8 $HOME_DIR/refs/hg19/hg19.fa $HOME_DIR/refs/hg19/hg19\"\n",
    "    subprocess.call(command, shell=True)    \n",
    "    print(\"Human refenrence genome hg19 has been downloaded and index files generation complete.\")\n",
    "\n",
    "elif choice == 3:\n",
    "    print(\"You chose hg38 and miRBase v21 based genome alignment using Mirdeep*\")\n",
    "    with open('data/mirdeep.conf',\"w\") as handle:\n",
    "        handle.write(\"Human Genome = hg38\\nmiRBase = 21\")   \n",
    "    print(\"Downloading Mirdeep*....\")\n",
    "    command = \"\"\n",
    "    command += \"rm -rf $Tools_DIR/MDS_command_line_v38/MDS_command_line/genome/hg38 && \" \n",
    "    command += \"wget --content-disposition  http://sourceforge.net/projects/mirdeepstar/files/Index_files/hg38.zip/download -P $Tools_DIR/ && \"\n",
    "    command += \"unzip -o $Tools_DIR/hg38.zip -d $Tools_DIR/MDS_command_line_v38/MDS_command_line/genome/ && \"\n",
    "    command += \"rm -r $Tools_DIR/MDS_command_line_v38/MDS_command_line/genome/hg19\"\n",
    "    subprocess.call(command, shell=True)   \n",
    "    print(\"Downloading Mirdeep* is complete.\")\n",
    "    ref_var = 'hg38'   \n",
    "\n",
    "elif choice == 4:\n",
    "    print(\"You chose hg38 and miRBase v22 based genome alignment using Mirdeep*(Default Condition)\")\n",
    "    with open('data/mirdeep.conf',\"w\") as handle:\n",
    "        handle.write(\"Human Genome = hg38\\nmiRBase = 22\")   \n",
    "    ref_var = 'hg38'\n",
    "\n",
    "else:\n",
    "    print(\"Please enter valid option\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## piRNA Pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def piRNA_pipeline():\n",
    "    cmd = 'bash scripts/piRNA_pipeline.sh'\n",
    "    os.system(cmd)\n",
    "\n",
    "threading.Thread(target=piRNA_pipeline).start()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## miRNA sequence Alignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uncomment the following command if you want to perform sequence alignment sequentially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Sequential sample processing\n",
    "# !bash $HOME_DIR/scripts/Mirdeep_star.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Preparing batches for miRNA sequence alignment for multi-thread processing**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def isTreadAlive():\n",
    "    for t in threads:\n",
    "        if t.is_alive():\n",
    "            return 1\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mirdeep_star(batch,batch_id):    \n",
    "    os.chdir(os.environ['REF_DIR'])\n",
    "    files = os.listdir(os.getcwd())    \n",
    "    ref_id = os.path.join(os.getcwd(),'genome')\n",
    "    if ref_var == 'hg19':\n",
    "        ref_idx = 'hg19'\n",
    "    elif ref_var == 'hg38':\n",
    "        ref_idx = 'hg38'\n",
    "    else:\n",
    "        print('Reference index are missing in ', ref_id)\n",
    "        sys.exit(2)\n",
    "    \n",
    "    batch = pd.read_csv(os.path.join(os.environ['SEQ_DIR'],batch),index_col=0)\n",
    "    for idx in range(batch.shape[0]):\n",
    "        basename = batch.iloc[idx,0].split('.')[0] +'_trimmed.fastq'\n",
    "        path_new = os.path.join(os.environ['SEQ_DIR'],'fastq_17_24',basename)\n",
    "        \n",
    "        # Constraining miRDeep* to take only 4 gb in each thread\n",
    "        try:\n",
    "            command = 'java -Xmx4096m -jar MD.jar -g '+ ref_idx + ' -a ' + os.environ['adaptor'] + ' -t 17 -l 24 -s -20 -r 5 -p 20 -m 101 ' + path_new + ' && rm ' + path_new\n",
    "            os.system(command)\n",
    "        except:\n",
    "            screenlock.acquire()\n",
    "            print('Sample ',basename,' is not present in the directory.')\n",
    "            screenlock.release()\n",
    "    \n",
    "    screenlock.acquire()\n",
    "    print('---------------------------')\n",
    "    print(batch_id, ' is complete.')\n",
    "    print('---------------------------')\n",
    "    screenlock.release()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Collecting the sample_list provided by user\n",
    "os.chdir(os.environ['HOME_DIR'])\n",
    "sample_list = pd.read_csv('data/sample_list.csv',index_col=0)\n",
    "\n",
    "# Generating batches\n",
    "if sample_list.shape[0] < 10:\n",
    "    no_of_batch = 5\n",
    "elif sample_list.shape[0] > 10 and sample_list.shape[0] <= 30:\n",
    "    no_of_batch = 10\n",
    "elif sample_list.shape[0] > 30 and sample_list.shape[0] <= 50:\n",
    "    no_of_batch = 15\n",
    "elif sample_list.shape[0] > 50 and sample_list.shape[0] <= 100:\n",
    "    no_of_batch = 20\n",
    "elif sample_list.shape[0] > 100:\n",
    "    no_of_batch = 25\n",
    "    \n",
    "file_list_splitted = np.array_split(sample_list, no_of_batch)\n",
    "for batch_idx in range(len(file_list_splitted)):\n",
    "    file_list_splitted[batch_idx].to_csv('data/batch_'+str(batch_idx+1)+'.csv',encoding='utf-8',index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calling mirdeep* for each batch in individual threads.\n",
    "threads = []\n",
    "for i in range(1,no_of_batch+1):\n",
    "    print('Running batch_'+ str(i)+ ' on thread-'+ str(i))\n",
    "    t = threading.Thread(target=mirdeep_star,args=('batch_'+str(i)+'.csv', 'batch_'+str(i),))\n",
    "    threads.append(t)\n",
    "    t.start()\n",
    "\n",
    "flag =1\n",
    "while (flag):\n",
    "    time.sleep(5.0)\n",
    "    flag = isTreadAlive()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Post-Alignment Analysis\n",
    "## **Processing of raw counts from Mirdeep* results**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(os.getcwd())\n",
    "files = [] \n",
    "[files.append(i) for i in os.listdir(\"data/fastq_17_24\") if (\".result\" in i and not \"known\" in i)]\n",
    "result_data = {}\n",
    "print('Preparing master look up table for result data')\n",
    "for file in tqdm(files):    \n",
    "    result_data['_'.join(file.split(\"_\")[:-1])] = {}\n",
    "    readfile = open(os.path.join(\"data/fastq_17_24\", file), 'r').readlines()\n",
    "    header = readfile[0].split(\"\\t\")\n",
    "    for line in readfile[1:]:\n",
    "        result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]] = {}\n",
    "        result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[1]] = line.split(\"\\t\")[1]\n",
    "        result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[2]] = line.split(\"\\t\")[2]\n",
    "        result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[3]] = line.split(\"\\t\")[3]\n",
    "        result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[4]] = line.split(\"\\t\")[4]\n",
    "        result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[5]] = int(line.split(\"\\t\")[5])\n",
    "        result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[6]] = line.split(\"\\t\")[6]\n",
    "        result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[7]] = line.split(\"\\t\")[7]\n",
    "        result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[8]] = line.split(\"\\t\")[8]\n",
    "        result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[9]] = line.split(\"\\t\")[9]\n",
    "        result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[10]] = line.split(\"\\t\")[10]\n",
    "\n",
    "    known_miR_filename = file.split('.')[0] + '.known_miR.result'\n",
    "    readfile = open(os.path.join(\"data/fastq_17_24\", known_miR_filename), 'r').readlines()\n",
    "    for line in readfile[1:]:\n",
    "        if not line.split(\"\\t\")[0] in result_data['_'.join(file.split(\"_\")[:-1])].keys():\n",
    "            result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]] = {}\n",
    "            result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[1]] = line.split(\"\\t\")[1]\n",
    "            result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[2]] = line.split(\"\\t\")[2]\n",
    "            result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[3]] = line.split(\"\\t\")[3]\n",
    "            result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[4]] = line.split(\"\\t\")[4]\n",
    "            result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[5]] = int(line.split(\"\\t\")[5])\n",
    "            result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[6]] = line.split(\"\\t\")[6]\n",
    "            result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[7]] = line.split(\"\\t\")[7]\n",
    "            result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[8]] = line.split(\"\\t\")[8]\n",
    "            result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[9]] = line.split(\"\\t\")[9].split('(')[0].lower()\n",
    "            result_data['_'.join(file.split(\"_\")[:-1])][line.split(\"\\t\")[0]][header[10]] = line.split(\"\\t\")[10]\n",
    "        \n",
    "print('Master lookup table is generated. Now collecting all unique miRs from master lookup table')\n",
    "mir_dict = {}\n",
    "for file in tqdm(result_data.keys()):\n",
    "    # print('Working on ',file)\n",
    "    for k,v in result_data[file].items():\n",
    "        i = 1\n",
    "        if not k in mir_dict.keys():\n",
    "            mir_dict[k] = [v['chr'],v['mature_loci'].split('-')[0],\n",
    "                           v['mature_loci'].split('-')[1],v['mature miR'],v['sequence']]\n",
    "        else:\n",
    "            mir_dict1 = mir_dict.copy()\n",
    "            for k1,v1 in mir_dict1.items():\n",
    "                new_loc_start = int(v['mature_loci'].split('-')[0])\n",
    "                new_loc_end = int(v['mature_loci'].split('-')[0])\n",
    "                new_mir_seq = v['mature miR']\n",
    "                if k1 == k:                    \n",
    "                    if v1[0] == v['chr']:                          \n",
    "                        if (new_loc_start == int(v1[1])) and (new_loc_end == int(v1[2])):\n",
    "                            if new_mir_seq == v1[3]:\n",
    "                                pass\n",
    "                    else:\n",
    "                        mir_dict[k+'_'+str(i)] = [v['chr'],v['mature_loci'].split('-')[0],\n",
    "                                       v['mature_loci'].split('-')[1],v['mature miR'],v['sequence']]\n",
    "                        i += 1\n",
    "\n",
    "\n",
    "print('Prepating dataframe for all unique miRs and samples....')\n",
    "counts = pd.DataFrame(columns=list(result_data.keys()))\n",
    "mir_name = list(mir_dict.keys())\n",
    "index = 0\n",
    "\n",
    "for k,v in tqdm(mir_dict.items()):\n",
    "    for file in result_data.keys():\n",
    "        for k1,v1 in result_data[file].items():\n",
    "            if v1['mature miR'] == v[3]:\n",
    "                loc_start = int(v1['mature_loci'].split('-')[0])\n",
    "                loc_end = int(v1['mature_loci'].split('-')[1])\n",
    "                loc_chr = v1['chr']\n",
    "                if loc_chr == v[0]:\n",
    "                    if loc_start == int(v[1]) and loc_end == int(v[2]):   \n",
    "                        counts.loc[k,file] = int(v1['expression(number of mature reads)'])\n",
    "                   \n",
    "counts['mir_ID'] = mir_name\n",
    "counts = counts.set_index('mir_ID')\n",
    "counts['chr'] = [v[0] for v in mir_dict.values()]\n",
    "counts['chr_start'] = [v[1] for v in mir_dict.values()]\n",
    "counts['chr_end'] = [v[2] for v in mir_dict.values()]\n",
    "counts['miRNA_sequence'] = [v[3] for v in mir_dict.values()]\n",
    "counts['miRNA_precursor_sequence'] = [v[4].lower() for v in mir_dict.values()]\n",
    "# rearranging columns\n",
    "cols = counts.columns.tolist()\n",
    "cols = cols[-5:]+ cols[:-5]\n",
    "counts = counts[cols]\n",
    "\n",
    "counts = counts.fillna(0)\n",
    "print(counts.shape)\n",
    "counts.to_csv(\"data/counts.csv\",encoding='utf-8',index=True)\n",
    "counts.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_known = counts[counts.index.str.contains('hsa')]\n",
    "counts_unknown = counts[~counts.index.str.contains('hsa')]\n",
    "print('All %d known and %d novel miRNA are successfully separated' %(counts_known.shape[0],counts_unknown.shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Blast Search\n",
    "## Search for Novel miRNA sequence in DashR Database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "   Run this cell only when running this notebook first time. After successfully running this, \n",
    "   pickle file will be saved which can be used further without generating the \n",
    "   results again and will save lot of time.\n",
    "   This code has been tested with Firefox version 66.0.4 (64-bit).\n",
    "'''\n",
    "\n",
    "novel_mirs = list(counts_unknown.index.values)\n",
    "gecko = os.path.abspath('geckodriver')\n",
    "dashR_Results = []\n",
    "browser = Browser('firefox',executable_path=gecko,headless=True)\n",
    "browser.visit('http://dashr2.lisanwanglab.org/search.php#')\n",
    "\n",
    "if ref_var == 'hg19':\n",
    "    #  Select DASHR2 GEO HG19 as reference\n",
    "    xpath = '/html/body/div/div[1]/div/div/select/option[3]' \n",
    "else:\n",
    "    #  Select DASHR2 GEO HG38 as reference\n",
    "    xpath = '/html/body/div/div[1]/div/div/select/option[4]'\n",
    "browser.find_by_xpath(xpath).click()\n",
    "time.sleep(3)\n",
    "\n",
    "browser.find_by_text('Search by sequence ').click()\n",
    "time.sleep(3)\n",
    "n_idx = len(list(counts_unknown.loc[:,'miRNA_sequence']))\n",
    "n = 0\n",
    "for idx in tqdm(range(n_idx)):\n",
    "    name = list(counts_unknown.index.values)[idx]\n",
    "    seq = list(counts_unknown.loc[:,'miRNA_sequence'])[idx]\n",
    "    dashR_Results.append(name)\n",
    "    browser.fill('querySeq', seq)\n",
    "    time.sleep(2)\n",
    "    xpath = '//*[@id=\"search\"]/div[4]/div/button'    \n",
    "    browser.find_by_xpath(xpath).click()\n",
    "    if not browser.is_text_present('No matches found'):\n",
    "        n += 1\n",
    "        xpath = '//*[@id=\"sequence-results\"]/pre/table' \n",
    "        results = browser.find_by_xpath(xpath)\n",
    "        i=0\n",
    "        for search_result in results:\n",
    "            title = search_result.text.encode('utf8') \n",
    "            link = search_result[\"href\"]             \n",
    "            dashR_Results.append((title, link)) \n",
    "            i += 1   \n",
    "    else:\n",
    "        seq = str(Seq(seq).reverse_complement())\n",
    "        browser.fill('querySeq', seq)\n",
    "        time.sleep(2)\n",
    "        xpath = '//*[@id=\"search\"]/div[4]/div/button'    \n",
    "        browser.find_by_xpath(xpath).click()\n",
    "        if not browser.is_text_present('No matches found'):\n",
    "            n += 1\n",
    "            xpath = '//*[@id=\"sequence-results\"]/pre/table' \n",
    "            results = browser.find_by_xpath(xpath)\n",
    "            i=0\n",
    "            for search_result in results:\n",
    "                title = search_result.text.encode('utf8') \n",
    "                link = search_result[\"href\"]             \n",
    "                dashR_Results.append((title, link)) \n",
    "                i += 1\n",
    "        else:\n",
    "            dashR_Results.append('miRNA NOT FOUND')\n",
    "    \n",
    "\n",
    "\n",
    "browser.quit()\n",
    "with open(\"data/dashR_Results.pickle\", 'wb') as handle:\n",
    "    pickle.dump(dashR_Results, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uncomment the following cell if you want to load the pickle file obtained from dashR search module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# '''\n",
    "#    Get DashR database searching results from saved pickle file.   \n",
    "# '''\n",
    "# import pickle\n",
    "# print('Loading saved results from DashR Database....')\n",
    "# with open(\"data/dashR_Results.pickle\", 'rb') as handle:\n",
    "#     dashR_Results = pickle.load(handle)\n",
    "# len(dashR_Results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DashR Results post-processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loc_check(rna_name, df1,idx):\n",
    "    \n",
    "    rna_name1 = rna_name.split(' ')[1]\n",
    "    if '3p' in rna_name1 or '5p' in rna_name1:\n",
    "        rna_name1 = rna_name1[:-3]\n",
    "    rna_name1 = rna_name1.lower()\n",
    "    \n",
    "    rna_chr1 = counts_unknown.loc[dashR_Results[idx],'chr']\n",
    "    rna_chr_start1 = int(counts_unknown.loc[dashR_Results[idx],'chr_start'])\n",
    "    rna_chr_end1 = int(counts_unknown.loc[dashR_Results[idx],'chr_end'])\n",
    "\n",
    "    if not counts_unknown.shape[0] == df1.shape[0]:\n",
    "        try:\n",
    "            rna_chr2 = df1.loc[rna_name1,'chr']\n",
    "            rna_chr_start2 = int(df1.loc[rna_name1,'chr_start'])\n",
    "            rna_chr_end2 = int(df1.loc[rna_name1,'chr_end'])\n",
    "        except:\n",
    "            rna_chr2 = df1.loc[dashR_Results[idx],'chr']\n",
    "            rna_chr_start2 = int(df1.loc[dashR_Results[idx],'chr_start'])\n",
    "            rna_chr_end2 = int(df1.loc[dashR_Results[idx],'chr_end'])\n",
    "\n",
    "    else:\n",
    "        rna_chr2 = rna_name.split(' ')[3].split(':')[0]\n",
    "        rna_chr_start2 = int(rna_name.split(' ')[3].split(':')[1].split('[')[0].split('-')[0])\n",
    "        rna_chr_end2 = int(rna_name.split(' ')[3].split(':')[1].split('[')[0].split('-')[1])\n",
    "    \n",
    "    if (rna_chr1 == rna_chr2) and np.abs(rna_chr_start1 - rna_chr_start2)<=2 and np.abs(rna_chr_end1 - rna_chr_end2)<=2:\n",
    "        flag = True\n",
    "        return flag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_unknown2 = counts_unknown.copy()\n",
    "counts_known2 = counts_known.copy()\n",
    "counts_unknown2 = counts_unknown2.drop(['miRNA_sequence','chr','chr_start','chr_end','miRNA_precursor_sequence'],axis=1)\n",
    "counts_known2 = counts_known2.drop(['miRNA_sequence','chr','chr_start','chr_end','miRNA_precursor_sequence'],axis=1)\n",
    "novel_mir_old_index = counts_unknown.index.tolist()\n",
    "paralogue_dict = {}\n",
    "\n",
    "for idx in (range(0,len(dashR_Results),2)):\n",
    "    results = dashR_Results[idx+1]    \n",
    "    if isinstance((dashR_Results[idx+1]),tuple):\n",
    "        results = list(dashR_Results[idx+1])[0].decode(\"utf-8\")        \n",
    "        rna_type_results = results.split('\\n')[-1]\n",
    "        rna_type1 = rna_type_results.split(' ')[1].lower()\n",
    "        char_no = 1\n",
    "        if 'hsa-mir' in rna_type1:                            \n",
    "            idx2 = list(counts_unknown2.index.values).index(dashR_Results[idx])\n",
    "            count_to_add = list(map(int,list(counts_unknown2.iloc[idx2,:].values)))\n",
    "            if '3p' in rna_type1 or '5p' in rna_type1:\n",
    "                rna_type1 = rna_type1[:-3]\n",
    "            if rna_type1.lower() in list(counts_known.index.values):    \n",
    "                counts_add_flag = loc_check(rna_type_results,counts_known,idx)\n",
    "                if counts_add_flag:                        \n",
    "                    count_increament = list(counts_known2.loc[rna_type1.lower(),:].values)                        \n",
    "                    count_increament = list(map(add,count_increament,count_to_add))               \n",
    "                    counts_known2.iloc[list(counts_known2.index.values).index(rna_type1.lower()),:] = count_increament\n",
    "                    counts_unknown2 = counts_unknown2.drop([str(dashR_Results[idx])],axis = 0)                    \n",
    "                    novel_mir_old_index.remove(novel_mir_old_index[novel_mir_old_index.index(dashR_Results[idx])])\n",
    "                else:                    \n",
    "                    if rna_type1 in paralogue_dict.keys():\n",
    "                        char_string = paralogue_dict[rna_type1][-1].split('_')[-1]\n",
    "                        char_no = int(char_string) + 1\n",
    "                    else:\n",
    "                        char_no = 1\n",
    "\n",
    "                    new_name = rna_type1 + '_' + str(char_no)\n",
    "                    unique_flag = False\n",
    "                    while unique_flag == False:\n",
    "                        if new_name in counts_known2.index.tolist():\n",
    "                            new_name = new_name[:-2]+'_' + str(int(new_name[-1])+1)\n",
    "                        else:\n",
    "                            unique_flag = True\n",
    "\n",
    "                    old_index = counts_unknown2.index.tolist()                        \n",
    "                    old_index[idx2] = new_name\n",
    "                    counts_unknown2.index = old_index\n",
    "                    if rna_type1 in paralogue_dict.keys():\n",
    "                        paralogue_dict[rna_type1].append(new_name)\n",
    "                    else:\n",
    "                        paralogue_dict[rna_type1] = []\n",
    "                        paralogue_dict[rna_type1].append(new_name)\n",
    "            else:\n",
    "                mir_add_flag = loc_check(rna_type_results,counts_unknown,idx)\n",
    "                if mir_add_flag:\n",
    "                    char_no += 1\n",
    "                    old_index = counts_unknown2.index.tolist()\n",
    "                    # Novel mir name will be changed to known mir name\n",
    "                    if old_index[idx2] in counts_known2.index.tolist():\n",
    "                        old_index[idx2] = old_index[idx2][:-2] + '_'+ str(char_no)\n",
    "                        char_no += 1\n",
    "                    \n",
    "                    old_index[idx2] = rna_type1[:-2] + '_'+ str(char_no)\n",
    "                    counts_unknown2.index = old_index\n",
    "                else:\n",
    "                    # Novel miRNA is the paralogues of known mirna\n",
    "                    if rna_type1 in paralogue_dict.keys():\n",
    "                        char_string = paralogue_dict[rna_type1][-1].split('_')[-1]\n",
    "                        char_no = int(char_string) + 1\n",
    "                    else:\n",
    "                        char_no = 1\n",
    "\n",
    "                    new_name = rna_type1 + '_' + str(char_no)\n",
    "                    unique_flag = False\n",
    "                    while unique_flag == False:\n",
    "                        if new_name in counts_known2.index.tolist():\n",
    "                            new_name = new_name[:-2]+'_' + str(int(new_name[-1])+1)\n",
    "                        else:\n",
    "                            unique_flag = True\n",
    "                    old_index = counts_unknown2.index.tolist()                        \n",
    "                    old_index[idx2] = new_name\n",
    "                    counts_unknown2.index = old_index \n",
    "                    if rna_type1 in paralogue_dict.keys():\n",
    "                        paralogue_dict[rna_type1].append(new_name)\n",
    "                    else:\n",
    "                        paralogue_dict[rna_type1] = []\n",
    "                        paralogue_dict[rna_type1].append(new_name)\n",
    "            \n",
    "\n",
    "counts_known2['chr'] = counts_known[[\"chr\"]]\n",
    "counts_known2['chr_start'] = counts_known[[\"chr_start\"]]\n",
    "counts_known2['chr_end'] = counts_known[[\"chr_end\"]]\n",
    "counts_known2['miRNA_sequence'] = counts_known[[\"miRNA_sequence\"]]\n",
    "counts_known2['miRNA_precursor_sequence'] = counts_known[[\"miRNA_precursor_sequence\"]]\n",
    "\n",
    "new_mir_list = counts_unknown2.index.tolist()\n",
    "counts_unknown2.index = novel_mir_old_index\n",
    "counts_unknown2['chr'] = counts_unknown[[\"chr\"]]\n",
    "counts_unknown2['chr_start'] = counts_unknown[[\"chr_start\"]]\n",
    "counts_unknown2['chr_end'] = counts_unknown[[\"chr_end\"]]\n",
    "counts_unknown2['miRNA_sequence'] = counts_unknown[[\"miRNA_sequence\"]]\n",
    "counts_unknown2['miRNA_precursor_sequence'] = counts_unknown[[\"miRNA_precursor_sequence\"]]\n",
    "counts_unknown2.index = new_mir_list\n",
    "\n",
    "# rearranging columns\n",
    "cols = counts_known2.columns.tolist()\n",
    "cols = cols[-5:]+ cols[:-5]\n",
    "counts_known2 = counts_known2[cols]\n",
    "\n",
    "cols = counts_unknown2.columns.tolist()\n",
    "cols = cols[-5:]+ cols[:-5]\n",
    "counts_unknown2 = counts_unknown2[cols]\n",
    "\n",
    "# add 2 dfs\n",
    "new_mirs_df = counts_unknown2[counts_unknown2.index.str.contains('hsa')]\n",
    "frames = [counts_known2, new_mirs_df]\n",
    "counts_known2 = pd.concat(frames)\n",
    "counts_unknown2 = counts_unknown2[~counts_unknown2.index.str.contains('hsa')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "known_mirna_paralogue = counts_known2[['miRNA_sequence']]\n",
    "known_mirna_paralogue = pd.DataFrame(known_mirna_paralogue.groupby('miRNA_sequence').apply(lambda x: x.index.tolist()).reset_index(name='[mir_name]'))\n",
    "known_mirna_paralogue.to_csv(\"data/list_of_known_mirna_paralogues.csv\",index=False,encoding='utf-8')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Seed based Sequence Clustering\n",
    "## **Dictionary Preparation**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_known2 = counts_known2.drop_duplicates(subset=['miRNA_sequence'], keep='first')\n",
    "seed_dict = {}\n",
    "Xseed_dict = {}\n",
    "unknown_mir_dict = {}\n",
    "for index in counts_unknown2.index:\n",
    "    seed = counts_unknown2.loc[index,'miRNA_sequence'][1:7]\n",
    "    Xseed = counts_unknown2.loc[index,'miRNA_sequence'][7:]\n",
    "    unknown_mir_dict[index] = [seed,counts_unknown2.loc[index,'chr'],counts_unknown2.loc[index,'chr_start'],\n",
    "                              counts_unknown2.loc[index,'chr_end'],counts_unknown2.loc[index,'miRNA_sequence']]\n",
    "    seed_dict[index] = seed\n",
    "    Xseed_dict[index] = Xseed\n",
    "\n",
    "for index in counts_known2.index:\n",
    "    seed = counts_known2.loc[index,'miRNA_sequence'][1:7]\n",
    "    Xseed = counts_known2.loc[index,'miRNA_sequence'][7:]\n",
    "    unknown_mir_dict[index] = [seed,counts_known2.loc[index,'chr'],counts_known2.loc[index,'chr_start'],\n",
    "                              counts_known2.loc[index,'chr_end'],counts_known2.loc[index,'miRNA_sequence']]\n",
    "    seed_dict[index] = seed\n",
    "    Xseed_dict[index] = Xseed\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import defaultdict \n",
    "\n",
    "exact_seed_cluster = defaultdict(list) \n",
    "for key,values in seed_dict.items():\n",
    "    exact_seed_cluster[values].append(key)\n",
    "    \n",
    "exact_seed_cluster = dict(exact_seed_cluster)\n",
    "\n",
    "exact_Xseed_cluster = defaultdict(list) \n",
    "for key,values in Xseed_dict.items():\n",
    "    exact_Xseed_cluster[values].append(key)\n",
    "    \n",
    "exact_Xseed_cluster = dict(exact_Xseed_cluster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Fasta File**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "count = 1\n",
    "writer = open(\"data/seed_seqs.fasta\", 'w')\n",
    "for key in unknown_mir_dict.keys():\n",
    "    line = \">\" + str(count) + \"\\t\" + key + \"\\n\"\n",
    "    line += str(unknown_mir_dict[key][0]) + \"\\n\"\n",
    "    writer.write(line)\n",
    "    count+=1\n",
    "writer.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "count = 1\n",
    "writer = open(\"data/seed_outer_seqs.fasta\", 'w')\n",
    "for key in unknown_mir_dict.keys():\n",
    "    line = \">\" + str(count) + \"\\t\" + key + \"\\n\"\n",
    "    line += str(unknown_mir_dict[key][-1][7:]) + \"\\n\"\n",
    "    writer.write(line)\n",
    "    count+=1\n",
    "writer.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CD-HIT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute CD-HIT command\n",
    "\n",
    "!cd Tools/cd-hit-v4.6.8-2017-1208/; echo \"***** Clustering miRs Seed Sequences Only**********\"; \\\n",
    "                                                           ./cd-hit -l 5 -n 2 -c 1 \\\n",
    "                                                                    -i ../../data/seed_seqs.fasta \\\n",
    "                                                                    -o ../../data/cdhit_seed;     \\\n",
    "                                                           echo \" \"; \\\n",
    "                                                           echo \"***** Clustering miRs Rest Sequences Only**********\"; \\\n",
    "                                                           ./cd-hit -l 8 -i ../../data/seed_outer_seqs.fasta \\\n",
    "                                                                    -o ../../data/cdhit_seed_outer;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Seed Clustering Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Preparing lookup table for Xseed Clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Cluster dictionary preparation\n",
    "Cluster : cluster_no:[mir1,mir2]\n",
    "\"\"\"\n",
    "def fetch_id(xseed_id_fasta, xid):\n",
    "    for line in xseed_id_fasta:\n",
    "        if line:\n",
    "            if '>' in line[0]:   \n",
    "                if xid == line.split('\\t')[0].split('>')[1]:\n",
    "                    return line.split('\\t')[1]\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "def cluster_dict_preparation(cluster_file, fasta_file):\n",
    "    cluster_file = open(cluster_file,'r').read().split('\\n')\n",
    "    fasta_id = open(fasta_file,'r').read().split('\\n')\n",
    "    cluster_dict = {}    \n",
    "    for line in cluster_file:\n",
    "        if line:\n",
    "            if '>' in line[0]:\n",
    "                line_no = cluster_file.index(line)\n",
    "                flag = True    \n",
    "                k = line[1:].replace(' ','_')\n",
    "                xid = cluster_file[line_no+1].split('>')[1].split('.')[0]\n",
    "                cluster_dict[k] = [fetch_id(fasta_id,xid)]\n",
    "                i = 2\n",
    "                while flag:\n",
    "                    if cluster_file[line_no+i]:\n",
    "                        if '>' in cluster_file[line_no+i][0]:\n",
    "                            flag = False\n",
    "                            pass\n",
    "                        else:\n",
    "                            xid = cluster_file[line_no+i].split('>')[1].split('.')[0]\n",
    "                            cluster_dict[k].append(fetch_id(fasta_id,xid))\n",
    "                            i += 1\n",
    "                    else:\n",
    "                        flag = False\n",
    "                        \n",
    "    return cluster_dict\n",
    "\n",
    "\n",
    "seed_cluster = cluster_dict_preparation('data/cdhit_seed.clstr','data/seed_seqs.fasta')\n",
    "Xseed_cluster = cluster_dict_preparation('data/cdhit_seed_outer.clstr','data/seed_outer_seqs.fasta')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Only for outer seed region clusters\n",
    "rev_dict : mir_name:[cluster1,cluster2]\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "def rev_dictionary(cluster):\n",
    "    rev_dict = defaultdict(list)\n",
    "    mir_list = []\n",
    "    for k in list(unknown_mir_dict.keys()):\n",
    "        mir_list.append(k)\n",
    "        \n",
    "    for mir in mir_list:\n",
    "        for k,v in cluster.items():\n",
    "            for mir1 in v:\n",
    "                if mir == mir1:\n",
    "                    rev_dict[mir].append(k)\n",
    "                    \n",
    "    return dict(rev_dict)\n",
    "\n",
    "\n",
    "exact_seed_rev_dict = rev_dictionary(exact_seed_cluster)\n",
    "seed_rev_dict = rev_dictionary(seed_cluster)\n",
    "Xseed_rev_dict = rev_dictionary(Xseed_cluster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Function to update and merge the counts if all conditions (seed,seq and loc) are satisfied\n",
    "Input : df1,df2,mir1 (from df1), mir2 (from df2)\n",
    "Output : df1_modified (with updated counts of mir1) and df2_modified (with mir2 removed)\n",
    "\"\"\"\n",
    "\n",
    "def update_counts(df1,df2,mir1,mir2):\n",
    "    mir1_counts = list(df1.loc[mir1,:].values)[5:]\n",
    "    mir2_counts = list(df2.loc[mir2,:].values)[5:]\n",
    "    mir2_counts = list(map(add,mir1_counts,mir2_counts))\n",
    "    df1.loc[mir2,:] = list(df1.loc[mir2,:].values)[:5] + mir2_counts\n",
    "    df2 = df2.drop(mir1,axis = 0)\n",
    "    return df1, df2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "novel_miRNA2 = counts_unknown2.copy()\n",
    "mir_reannotation = pd.DataFrame()\n",
    "mir_old_name = []\n",
    "mir_new_name = []\n",
    "novel_mir_seq = []\n",
    "family_name = []\n",
    "attribute = []\n",
    "paralogue_dict = {}\n",
    "\n",
    "for k,v in tqdm(seed_cluster.items()):    \n",
    "    char_no = 0\n",
    "    novel_mir_name = seed_cluster[k]\n",
    "    novel_mir_name = [i for i in novel_mir_name if not 'hsa' in i]\n",
    "    for n_mir in novel_mir_name:   \n",
    "        novel_mir_name = [i for i in novel_mir_name if i is not None]                        \n",
    "        if len(novel_mir_name)>1:  \n",
    "            for n_mir in novel_mir_name:\n",
    "                if not n_mir == None:\n",
    "                    char_no = 1\n",
    "                    novel_mir_name1 = novel_mir_name.copy()\n",
    "                    novel_mir_name1.remove(n_mir)\n",
    "                    n_chr_loc = unknown_mir_dict[n_mir][1]\n",
    "                    n_chr_loc_start = int(unknown_mir_dict[n_mir][2])\n",
    "                    n_chr_loc_end = int(unknown_mir_dict[n_mir][3])\n",
    "                    seed_common_mir = list(set(novel_mir_name1) & set(exact_seed_cluster[exact_seed_rev_dict[n_mir][0]]))\n",
    "                    if seed_common_mir:                    \n",
    "                        Xseed_common_mir = list(set(seed_common_mir) & set(Xseed_cluster[Xseed_rev_dict[n_mir][0]]))\n",
    "                        if Xseed_common_mir:\n",
    "                            for c_mir1 in Xseed_common_mir:\n",
    "                                c_chr_loc = unknown_mir_dict[c_mir1][1]\n",
    "                                c_chr_loc_start = int(unknown_mir_dict[c_mir1][2])\n",
    "                                c_chr_loc_end = int(unknown_mir_dict[c_mir1][3])\n",
    "                                if (n_chr_loc == c_chr_loc) and np.abs(n_chr_loc_start - c_chr_loc_start)<=2 and np.abs(n_chr_loc_end - c_chr_loc_end)<=2:\n",
    "                                    _,novel_miRNA2 = update_counts(novel_miRNA2,novel_miRNA2,c_mir1,n_mir)\n",
    "                                    novel_mir_name[novel_mir_name.index(c_mir1)] = None\n",
    "                                    mir_old_name.append(c_mir1)\n",
    "                                    mir_new_name.append(n_mir)\n",
    "                                    family_name.append(n_mir)\n",
    "                                    novel_mir_seq.append(unknown_mir_dict[n_mir][-1])\n",
    "                                    attribute.append('Same Seed, Xseed and loc')\n",
    "                                else:\n",
    "                                    if n_mir in paralogue_dict.keys():\n",
    "                                        char_string = paralogue_dict[n_mir][-1].split('_')[-1]\n",
    "                                        char_no = int(char_string) + 1\n",
    "                                    else:\n",
    "                                        char_no = 1\n",
    "                                        \n",
    "                                    new_name = n_mir+'_'+ str(char_no)\n",
    "                                    unique_flag = False\n",
    "                                    while unique_flag == False:\n",
    "                                        if new_name in novel_miRNA2.index.tolist():\n",
    "                                            new_name = new_name[:-2]+'_' + str(int(new_name[-1])+1)\n",
    "                                        else:\n",
    "                                            unique_flag = True\n",
    "                                    novel_miRNA2 = novel_miRNA2.rename(index={c_mir1 : new_name})\n",
    "                                    \"\"\"\n",
    "                                        The incoming miR shared same seed, altered xseed with different genomic location with cluster head.\n",
    "                                        We need to check whether the incoming miR shares the genomic location with the existing \n",
    "                                        members of the same cluster or not. If not, then it'll be added to the cluster with \n",
    "                                        different name. Otherwise, it'll be merged to the member with whom it shares the exact \n",
    "                                        seed, xseed and genomic location.\n",
    "                                    \"\"\"\n",
    "                                    \n",
    "                                    if n_mir in paralogue_dict.keys():\n",
    "                                        if len(paralogue_dict[n_mir]) >= 1:\n",
    "                                            for m in paralogue_dict[n_mir]:\n",
    "                                                if m in novel_miRNA2.index.tolist():\n",
    "                                                    if (novel_miRNA2.loc[m,'chr'] == c_chr_loc) and np.abs(int(novel_miRNA2.loc[m,'chr_start']) - int(novel_miRNA2.loc[new_name,'chr_start']))<=2 and np.abs(int(novel_miRNA2.loc[m,'chr_end']) - int(novel_miRNA2.loc[new_name,'chr_end']))<=2:\n",
    "                                                        _,novel_miRNA2 = update_counts(novel_miRNA2,novel_miRNA2,new_name,m)\n",
    "                                                        mir_old_name.append(c_mir1)\n",
    "                                                        mir_new_name.append(m)\n",
    "                                                        family_name.append(m)\n",
    "                                                        novel_mir_seq.append(unknown_mir_dict[n_mir][-1])\n",
    "                                                        attribute.append('Same Seed, Xseed and loc')\n",
    "                                                        break\n",
    "                                    novel_mir_name[novel_mir_name.index(c_mir1)] = None\n",
    "                                    mir_old_name.append(c_mir1)\n",
    "                                    mir_new_name.append(new_name)\n",
    "                                    family_name.append(n_mir)\n",
    "                                    novel_mir_seq.append(unknown_mir_dict[n_mir][-1])\n",
    "                                    attribute.append('Paralogues')                                    \n",
    "                                    if n_mir in paralogue_dict.keys():\n",
    "                                        paralogue_dict[n_mir].append(new_name)\n",
    "                                    else:\n",
    "                                        paralogue_dict[n_mir] = []\n",
    "                                        paralogue_dict[n_mir].append(new_name)\n",
    "                        else:\n",
    "                            for c_mir1 in seed_common_mir: \n",
    "                                c_chr_loc = unknown_mir_dict[c_mir1][1]\n",
    "                                c_chr_loc_start = int(unknown_mir_dict[c_mir1][2])\n",
    "                                c_chr_loc_end = int(unknown_mir_dict[c_mir1][3])                               \n",
    "                                if n_mir in paralogue_dict.keys():\n",
    "                                        char_string = paralogue_dict[n_mir][-1].split('_')[-1]\n",
    "                                        char_no = int(char_string) + 1\n",
    "                                else:\n",
    "                                    char_no = 1\n",
    "                                    \n",
    "                                new_name = n_mir+'_'+ str(char_no)\n",
    "                                unique_flag = False\n",
    "                                while unique_flag == False:\n",
    "                                    if new_name in novel_miRNA2.index.tolist():\n",
    "                                        new_name = new_name[:-2]+'_' + str(int(new_name[-1])+1)\n",
    "                                    else:\n",
    "                                        unique_flag = True\n",
    "                                novel_miRNA2 = novel_miRNA2.rename(index={c_mir1 : new_name})\n",
    "                                if n_mir in paralogue_dict.keys():\n",
    "                                        if len(paralogue_dict[n_mir]) >= 1:\n",
    "                                            for m in paralogue_dict[n_mir]:\n",
    "                                                if m in novel_miRNA2.index.tolist():\n",
    "                                                    if (novel_miRNA2.loc[m,'chr'] == c_chr_loc) and np.abs(int(novel_miRNA2.loc[m,'chr_start']) - int(novel_miRNA2.loc[new_name,'chr_start']))<=2 and np.abs(int(novel_miRNA2.loc[m,'chr_end']) - int(novel_miRNA2.loc[new_name,'chr_end']))<=2:\n",
    "                                                        _,novel_miRNA2 = update_counts(novel_miRNA2,novel_miRNA2,new_name,m)\n",
    "                                                        mir_old_name.append(c_mir1)\n",
    "                                                        mir_new_name.append(m)\n",
    "                                                        family_name.append(m)\n",
    "                                                        novel_mir_seq.append(unknown_mir_dict[n_mir][-1])\n",
    "                                                        attribute.append('Same Seed, Xseed and loc')\n",
    "                                                        break\n",
    "                                novel_mir_name[novel_mir_name.index(c_mir1)] = None\n",
    "                                mir_old_name.append(c_mir1)\n",
    "                                mir_new_name.append(new_name)\n",
    "                                family_name.append(n_mir)\n",
    "                                novel_mir_seq.append(unknown_mir_dict[n_mir][-1])\n",
    "                                attribute.append('Paralogues')                                    \n",
    "                                if n_mir in paralogue_dict.keys():\n",
    "                                    paralogue_dict[n_mir].append(new_name)\n",
    "                                else:\n",
    "                                    paralogue_dict[n_mir] = []\n",
    "                                    paralogue_dict[n_mir].append(new_name)\n",
    "\n",
    "                                \n",
    "\n",
    "            \n",
    "            \n",
    "mir_reannotation['Old_miR_ID'] = mir_old_name\n",
    "mir_reannotation['New_miR_ID'] = mir_new_name\n",
    "mir_reannotation['Family'] = family_name\n",
    "mir_reannotation['Relation'] = attribute\n",
    "mir_reannotation['sequence'] = novel_mir_seq\n",
    "mir_reannotation.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed_dict2 = {k:v for k,v in seed_dict.items() if 'hsa' in k}\n",
    "for idx in tqdm(range(novel_miRNA2.shape[0])):\n",
    "    novel_seq = novel_miRNA2.iloc[idx,3]\n",
    "    novel_seed = novel_seq[1:7]\n",
    "    suffix = 1\n",
    "    for k,v in seed_dict2.items():\n",
    "        if novel_seed == v:\n",
    "            if k+'_'+str(suffix) in novel_miRNA2.index:\n",
    "                paralog_index = novel_miRNA2[novel_miRNA2.index.str.contains(k)].index.tolist()\n",
    "                suffix = sorted([int(i.split('_')[-1]) for i in paralog_index])[-1]\n",
    "                new_name = k + '_' + str(suffix+1)\n",
    "            elif not '_' in k:\n",
    "                suffix = 1\n",
    "                new_name = k + '_' + str(suffix)\n",
    "            novel_miRNA2 = novel_miRNA2.rename(index={novel_miRNA2.index[idx] : new_name})\n",
    "            break    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "novel_miRNA3 = novel_miRNA2[['miRNA_sequence']]\n",
    "novel_miRNA3 = pd.DataFrame(novel_miRNA3.groupby('miRNA_sequence').apply(lambda x: x.index.tolist()).reset_index(name='[mir_name]'))\n",
    "novel_miRNA3.to_csv(\"data/list_of_novel_mir_paralogues.csv\",index=False,encoding='utf-8')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frame = [counts_known2,novel_miRNA2]\n",
    "counts_final = pd.concat(frame)\n",
    "counts_final = counts_final.drop_duplicates(subset=['miRNA_sequence'], keep='first')\n",
    "counts_final.to_csv(\"data/counts_final.csv\",encoding='utf-8',index=True)\n",
    "print(counts_final.shape)\n",
    "\n",
    "# Preparing the final count matrix for differential expression analysis\n",
    "counts_final = counts_final.reindex(sorted(counts_final.columns), axis=1)\n",
    "counts_final = counts_final.drop(['chr','chr_start','chr_end','miRNA_sequence','miRNA_precursor_sequence'],axis=1)\n",
    "counts_final.to_csv(\"data/counts_final_1.csv\",encoding='utf-8',index=True) #For DESeq2\n",
    "counts_final.to_csv(\"data/output/miRNA_expression_counts.csv\",encoding='utf-8',index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# miRNA Differential expression Analysis using DESeq2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-05-22T10:24:40.850010Z",
     "start_time": "2019-05-22T10:24:30.691078Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    Before reunning this module, please confirm that the order of sample \n",
    "    written in sample_list.csv and order of columns in counts_final_1.csv\n",
    "    is same.e.g.\n",
    "    In sample_list.csv the samples are:\n",
    "    Sample        File        Condition\n",
    "    SampleA   sampleA.fastq    treated\n",
    "    SampleB   sampleB.fastq    treated\n",
    "    SampleC   sampleC.fastq    treated\n",
    "    SampleD   sampleD.fastq    treated\n",
    "       .           .              .   \n",
    "       .           .              .   \n",
    "       .           .              .   \n",
    "       .           .              .   \n",
    "       \n",
    "    The column order in counts_final_1.csv should be:\n",
    "    \n",
    "       .    sampleA    sampleB    sampleC    sampleD    \n",
    "       .      25         102         10         5\n",
    "       .       .          .           .         .\n",
    "       .       .          .           .         .\n",
    "       .       .          .           .         .\n",
    "\"\"\"\n",
    "# !Rscript $HOME_DIR/scripts/install_packages.R\n",
    "\n",
    "!Rscript scripts/norm_diff_exp.R miRNA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Naming of sequences that are discovered by miRPipe as miRs but not found in miRBase"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Disclaimer**\n",
    "\n",
    "We have renamed the identified novel miRNAs according to the miRBase nomenclature system. The rules for miRNA nomenclature are as follows:\n",
    "\n",
    "    1. The paralogs of a miRNA will be renamed by adding \"_1\" or \"_2\" at the end of the actual miRNA name.\n",
    "    \n",
    "    2. If the novel miRNA sequence has been already annotated in another organism, then the same identifier will be used to other organism also. For this, we have implemented Rfam database search for all novel miRNAs with E-value less than or equal to 0.01 and then renamed them using the same identifier.\n",
    "    \n",
    "    3. If above conditions are not met for any novel miRNAs then the renaming was done sequentially ().\n",
    "\n",
    "\n",
    "**Note: The suffix \" * \" is added after each name that represents that these names are putative names only. Users are requested to check names according to the miRBase and Rfam latest version.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rfam Database Integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reference: https://rna-tools.readthedocs.io/en/latest/_modules/rna_tools/RfamSearch.html\n",
    "class RfamSearchError(Exception):\n",
    "    pass\n",
    "\n",
    "class RfamSearch():\n",
    "    \"\"\"RfamSearch (local).\n",
    "\n",
    "    Rfam is a collection of multiple sequence alignments and covariance models representing non-coding RNA families. Rfam is available on the web http://rfam.xfam.org/. The website allow the user to search a query sequence against a library of covariance models, and view multiple sequence alignments and family annotation. The database can also be downloaded in flatfile form and searched locally using the INFERNAL package (http://infernal.wustl.edu/). The first release of Rfam (1.0) contains 25 families, which annotate over 50 000 non-coding RNA genes in the taxonomic divisions of the EMBL nucleotide database.\n",
    "\n",
    "    Infernal (\"INFERence of RNA ALignment\") is for searching DNA sequence databases for RNA structure and sequence similarities. It is an implementation of a special case of profile stochastic context-free grammars called covariance models (CMs). A CM is like a sequence profile, but it scores a combination of sequence consensus and RNA secondary structure consensus, so in many cases, it is more capable of identifying RNA homologs that conserve their secondary structure more than their primary sequence.\n",
    "\n",
    "    Infernal `cmscan` is used to search the CM-format Rfam database.\n",
    "\n",
    "    Setup:\n",
    "\n",
    "    - download the database from ftp://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT (file: Rfam.cm.gz, ~30mb)\n",
    "    - install http://eddylab.org/infernal/\n",
    "    - set up ``RFAM_DB_PATH`` in the config file of rna-tools.\n",
    "    - compress Rfam.cm\n",
    "    \n",
    "    Example of compressing the database::\n",
    "\n",
    "         $ cmpress Rfam.cm\n",
    "         Working...    done.\n",
    "         Pressed and indexed 3016 CMs and p7 HMM filters (3016 names and 3016 accessions).\n",
    "         Covariance models and p7 filters pressed into binary file:  Rfam.cm.i1m\n",
    "         SSI index for binary covariance model file:                 Rfam.cm.i1i\n",
    "         Optimized p7 filter profiles (MSV part)  pressed into:      Rfam.cm.i1f\n",
    "         Optimized p7 filter profiles (remainder) pressed into:      Rfam.cm.i1p\n",
    "\n",
    "    Cite: Nawrocki and S. R. Eddy, Infernal 1.1: 100-fold faster RNA homology searches, Bioinformatics 29:2933-2935 (2013). \"\"\"  # noqa\n",
    "\n",
    "    def __init__(self):\n",
    "        pass\n",
    "\n",
    "\n",
    "    def cmscan(self, seq, verbose=False):\n",
    "        \"\"\"Run cmscan on the seq.\n",
    "\n",
    "        Usage::\n",
    "\n",
    "           >>> seq = RNASequence(\"GGCGCGGCACCGUCCGCGGAACAAACGG\")\n",
    "           >>> rs = RfamSearch()\n",
    "           >>> hit = rs.cmscan(seq)\n",
    "           >>> print(hit)  #doctest: +ELLIPSIS\n",
    "           # cmscan :: search sequence(s) against a CM database...\n",
    "\n",
    "        :param seq: string\n",
    "        :returns: result\n",
    "        :rtype: string \"\"\"\n",
    "\n",
    "        # make tmp file\n",
    "        tf = tempfile.NamedTemporaryFile(delete=False)\n",
    "        tf.name += '.fa'\n",
    "        with open(tf.name, 'w') as f:\n",
    "            f.write('>target\\n')\n",
    "            f.write(seq.seq + '\\n')\n",
    "\n",
    "        # make output file\n",
    "        of = tempfile.NamedTemporaryFile(delete=False)\n",
    "\n",
    "        # run cmscan\n",
    "        cmd = 'cmscan -E 0.01 refs/Rfam.cm' + ' ' + tf.name + '  > ' + of.name\n",
    "        o = subprocess.Popen(\n",
    "            cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "        err = o.stderr.read().strip()\n",
    "        if err:\n",
    "            raise RfamSearchError(err)\n",
    "        self.output = open(of.name).read()\n",
    "        return self.output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RfamBlastSearchResult(seq):\n",
    "    seq = RNASequence(seq)\n",
    "    rs = RfamSearch()\n",
    "    hit = rs.cmscan(seq)\n",
    "    result = hit.split('\\n')\n",
    "    if not 'No hits detected' in result[16]:\n",
    "        if not len(str.strip(result[16]))==0:\n",
    "            try:\n",
    "                str_list = list(filter(None, result[16].split(' ')))\n",
    "                if float(str_list[2]) <= 0.1:\n",
    "                    if 'mir' in str_list[5] or 'let' in str_list[5]:\n",
    "                        mir_name = str_list[5]\n",
    "                    else:\n",
    "                        mir_name = 'None'\n",
    "            except:\n",
    "                str_list = list(filter(None, result[15].split(' ')))\n",
    "                if float(str_list[2]) <= 0.1:\n",
    "                    if 'mir' in str_list[5] or 'let' in str_list[5]:\n",
    "                        mir_name = str_list[5]\n",
    "                    else:\n",
    "                        mir_name = 'None'\n",
    "        else:\n",
    "            str_list = list(filter(None, result[15].split(' ')))\n",
    "            if float(str_list[2]) <= 0.1:\n",
    "                if 'mir' in str_list[5] or 'let' in str_list[5]:\n",
    "                    mir_name = str_list[5]\n",
    "                else:\n",
    "                    mir_name = 'None'\n",
    "    else:\n",
    "        mir_name = 'None'\n",
    "        \n",
    "    return mir_name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Renaming of dysregualted novel miRNAs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = pd.read_csv('data/miRNA_significantly_DE_mir.csv', sep=',',header=None)\n",
    "df2 = pd.read_csv('data/counts_final.csv', sep=',',header=0)\n",
    "try:    \n",
    "    df2 = df2.set_index(['mirID'])\n",
    "except:\n",
    "    a = df2[df2.columns[df2.columns.str.startswith('Unnamed:')]]\n",
    "    df2['mirID'] = list(a.loc[:,list(a.columns.values)[0]].values)\n",
    "\n",
    "    # Reomove Unnamed column\n",
    "    df2 = df2[df2.columns[~df2.columns.str.startswith('Unnamed:')]]\n",
    "    df2 = df2.set_index(['mirID'])\n",
    "df3 = pd.DataFrame()\n",
    "for idx in range(df1.shape[0]):\n",
    "    idx2 = list(df2.index.values).index(df1.iloc[idx,0])   \n",
    "    df3 = df3.append(df2.iloc[idx2,:])[df2.columns.tolist()]\n",
    "df3 = df3.reindex(sorted(df3.index), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df6 = df3[df3.index.str.contains('novel')]\n",
    "\n",
    "new_name_dict = {}\n",
    "\"\"\"\n",
    "    First we'll each novel miRNA sequence in Rfam database and if\n",
    "    it matched to any miRNA of any other organism then the novel\n",
    "    miRNA will be annotated using the same identifier.\n",
    "\"\"\"\n",
    "print('Searching novel miRNA sequence using blast search at Rfam database....')\n",
    "for s in tqdm(df6['miRNA_sequence'].to_list()):\n",
    "    name2 = RfamBlastSearchResult(s)\n",
    "    if not 'None' in name2:\n",
    "        if 'let' in name2:\n",
    "            if not s in new_name_dict.keys():\n",
    "                new_name_dict[s] = []\n",
    "                new_name_dict[s].append('hsa-let-'+name2.split('-')[1]+str(j)+'*')\n",
    "            else:\n",
    "                j = len(new_name_dict[s]) + 1\n",
    "                new_name_dict[s].append('hsa-let-'+name2.split('-')[1]+'-'+str(j)+'*')\n",
    "        else:\n",
    "            if not s in new_name_dict.keys():\n",
    "                new_name_dict[s] = []\n",
    "                new_name_dict[s].append('hsa-miR-'+name2.split('-')[1]+str(j)+'*')\n",
    "            else:\n",
    "                j = len(new_name_dict[s]) + 1\n",
    "                new_name_dict[s].append('hsa-miR-'+name2.split('-')[1]+'-'+str(j)+'*')\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "    Now, we'll rename the remaining miRNAs sequentially.\n",
    "\n",
    "\"\"\"\n",
    "print('Blast search is complete. Now renaming the remaining novel miRNAs sequentially.')\n",
    "last_mir_number = 12136 # As per miRBase version-22, the hairpin precursor number ends at 12136\n",
    "uniq_seq = [v for v in list(set(df6['miRNA_sequence'].to_list())) if not v in new_name_dict.keys()]\n",
    "name1 = 'hsa-mir-'\n",
    "new_name = [name1 + str(last_mir_number+i+1) for i in range(len(uniq_seq))]\n",
    "mir_name_dict = dict(zip(uniq_seq,new_name))\n",
    "\n",
    "new_name = []\n",
    "new_name2 = []\n",
    "j = 1\n",
    "\n",
    "for i in tqdm(range(df6.shape[0])):\n",
    "    if df6.iloc[i,3] in new_name_dict.keys():\n",
    "        new_name.append(new_name_dict[df6.iloc[i,3]])\n",
    "        if len(new_name_dict[df6.iloc[i,3]]) == 1:\n",
    "            new_name2.append(new_name_dict[df6.iloc[i,3]])\n",
    "        else:\n",
    "            last_element = [v for v in new_name2 if new_name_dict[df6.iloc[i,3] in v][-1]]\n",
    "            last_number = int(last_element.split('_')[-1])\n",
    "            new_name2.append(new_name_dict[df6.iloc[i,3]] + '_'+ str(j+1))\n",
    "    else:\n",
    "        if not mir_name_dict[df6.iloc[i,3]] in new_name:\n",
    "            new_name.append(mir_name_dict[df6.iloc[i,3]])\n",
    "            new_name2.append(mir_name_dict[df6.iloc[i,3]]+'*')\n",
    "        else:\n",
    "            new_name1 = mir_name_dict[df6.iloc[i,3]] + '_'+ str(j) + '*'\n",
    "            j += 1\n",
    "            new_name.append(mir_name_dict[df6.iloc[i,3]])\n",
    "            new_name2.append(new_name1)\n",
    "\n",
    "df6['new_name'] = new_name2\n",
    "\n",
    "cols = df6.columns.tolist()\n",
    "cols = [cols[-1]] + cols[:-1]\n",
    "df6 = df6[cols]\n",
    "df6.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Final output files preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_names = []\n",
    "for n in df3.index.tolist():\n",
    "    if 'hsa' in n:\n",
    "        new_names.append('-')\n",
    "for n in df6['new_name'].values.tolist():\n",
    "    new_names.append(n)\n",
    "\n",
    "df3['New name assigned to novel miRNAs'] = new_names\n",
    "cols = df3.columns.tolist()\n",
    "cols = [cols[-1]] + cols[:-1]\n",
    "df3 = df3[cols]\n",
    "\n",
    "mir_para = []\n",
    "for idx in tqdm(range(df3.shape[0])):\n",
    "    seq = df3.iloc[idx,4]\n",
    "    try:        \n",
    "        known_para = known_mirna_paralogue.iloc[known_mirna_paralogue['miRNA_sequence'].values.tolist().index(seq),1]\n",
    "        if len(known_para) == 1:\n",
    "            mir_para.append(known_para[0])\n",
    "        else:\n",
    "            mir_para.append(';'.join(known_para))\n",
    "    except:        \n",
    "        novel_para = novel_miRNA3.iloc[novel_miRNA3['miRNA_sequence'].values.tolist().index(seq),1]\n",
    "        if len(novel_para) == 1:\n",
    "            mir_para.append(novel_para[0])\n",
    "        else:\n",
    "            mir_para.append(';'.join(novel_para))\n",
    "\n",
    "new_names2 = new_names.copy()\n",
    "for mir in mir_para:\n",
    "    mir_idx = mir_para.index(mir)\n",
    "    old_name = new_names2[mir_idx]\n",
    "    if not 'hsa' in mir and ';' in mir:\n",
    "        count = 1\n",
    "        mir_list = mir.split(';')\n",
    "        para_name_list = []\n",
    "        for idx in range(1,len(mir_list)):\n",
    "            para_name = old_name[:-1] + '_' + str(count) + '*'\n",
    "            count += 1\n",
    "            para_name_list.append(para_name)\n",
    "    \n",
    "        new_names2[mir_idx] = old_name +';'+ ';'.join(para_name_list)\n",
    "\n",
    "df3['New name assigned to novel miRNAs'] = new_names2\n",
    "old_index = df3.index.tolist()\n",
    "df3.index = mir_para\n",
    "df3 = df3.drop(['miRNA_precursor_sequence'],axis=1)\n",
    "df3.to_csv('data/output/diff_exp_miRNAs_expression_counts.csv',index=True,encoding='utf-8')\n",
    "df4 = pd.read_csv(\"data/miRNA_DESeq_expression_data_controlIsUnreated.csv\",sep=',',header=0,index_col=0)\n",
    "df4 = df4.set_index(['cts...1.'])\n",
    "\n",
    "p_adj = []\n",
    "up_down_expr = []\n",
    "fold_change = []\n",
    "SE_mir = list(old_index)\n",
    "for mir in SE_mir:    \n",
    "    idx = list(df4.index.values).index(mir)\n",
    "    p_adj_val = \"{0:.5f}\".format(df4.loc[str(mir),'padj'])\n",
    "    p_adj.append(p_adj_val)\n",
    "    if df4.loc[str(mir),'log2FoldChange']>0:\n",
    "        up_down_expr_val = 'up'\n",
    "    else:\n",
    "        up_down_expr_val = 'down'\n",
    "        \n",
    "    up_down_expr.append(up_down_expr_val)\n",
    "    fold_change.append(df4.loc[str(mir),'log2FoldChange'])\n",
    "        \n",
    "final_df = pd.DataFrame()\n",
    "final_df['Differentially expressed miRNA'] = mir_para\n",
    "new_names = []\n",
    "for n in df3.index.tolist():\n",
    "    if 'hsa' in n:\n",
    "        new_names.append('-')\n",
    "for n in df6['new_name'].values.tolist():\n",
    "    new_names.append(n)\n",
    "\n",
    "final_df['New name assigned to novel miRNAs'] = new_names2\n",
    "final_df['mature_sequence'] = df3['miRNA_sequence'].values.tolist()\n",
    "final_df['Fold Change'] = fold_change\n",
    "final_df['up/down regulation'] = up_down_expr\n",
    "final_df['p_adj_value'] = p_adj\n",
    "final_df = final_df.sort_values(by=['Differentially expressed miRNA'], ascending=True)\n",
    "final_df.to_csv('data/output/final_diff_exp_miRNAs.csv',index=False,encoding='utf-8')\n",
    "final_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# piRNA Counts Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.environ['COUNT_DIR'] = os.path.join(os.environ['HOME_DIR'],'data','piRNA','pirna_counts')   \n",
    "# os.chdir(os.environ['COUNT_DIR'])\n",
    "\n",
    "raw_counts = pd.DataFrame()\n",
    "for file in tqdm(os.listdir(os.environ['COUNT_DIR'])):\n",
    "    if file.endswith('_counts.txt'):\n",
    "        file_open = open(os.path.join(os.environ['COUNT_DIR'],file),\"r\").read().split(\"\\n\")        \n",
    "        pi_name = []\n",
    "        pi_count = []\n",
    "        for line in file_open:\n",
    "            if len(line) > 1:\n",
    "                pi_name.append(line.split(\"\\t\")[-2])\n",
    "                pi_count.append(line.split(\"\\t\")[-1])\n",
    "        raw_counts['_'.join(file.split(\"_\")[:-1])] = pi_count\n",
    "        \n",
    "\n",
    "raw_counts[\"piRNA_ID\"] = pi_name\n",
    "raw_counts = raw_counts.set_index([\"piRNA_ID\"])\n",
    "\n",
    "raw_counts = raw_counts.apply(pd.to_numeric)\n",
    "raw_counts = raw_counts.groupby(['piRNA_ID']).sum()\n",
    "\n",
    "\n",
    "# Count filtering : Removing those piRNA's which are not expressed in any sample\n",
    "raw_counts1 = pd.DataFrame(columns = raw_counts.columns)\n",
    "for index in list(raw_counts.index.values):\n",
    "    if sum([int(val) for val in raw_counts.loc[index,:].values]) == 0:\n",
    "        next\n",
    "    else:\n",
    "        raw_counts1.loc[index,:] = [int(val) for val in raw_counts.loc[index,:].values]\n",
    "        \n",
    "# Sorting columns\n",
    "raw_counts1 = raw_counts1.reindex(sorted(raw_counts1.columns), axis=1)\n",
    "# Saving Results\n",
    "raw_counts1.to_csv('data/piRNA/piRNA_raw_counts.csv',encoding='utf-8',index=True) # FOr DESeq2\n",
    "raw_counts1.to_csv('data/output/piRNA_raw_counts.csv',encoding='utf-8',index=True)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# piRNA Differential expression Analysis using DESeq2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!Rscript $HOME_DIR/scripts/norm_diff_exp.R piRNA\n",
    "shutil.copy2('data/piRNA/significantly_DE_piRNA.csv', 'data/output/significantly_DE_piRNA.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Deleting temporary files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.remove('data/counts.csv')\n",
    "for file in glob.glob('data/batch_*.csv'):\n",
    "    os.remove(file)\n",
    "\n",
    "for file in glob.glob('data/cdhit_*'):\n",
    "    os.remove(file)\n",
    "\n",
    "for file in glob.glob('data/seed_*'):\n",
    "    os.remove(file)\n",
    "\n",
    "os.remove('data/miRNA_significantly_DE_mir.csv')\n",
    "os.remove('data/counts_final.csv')\n",
    "os.remove('data/counts_final_1.csv')\n",
    "os.remove('data/list_of_known_mirna_paralogues.csv')\n",
    "os.remove('data/list_of_novel_mir_paralogues.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Output\n",
    "\n",
    "You can check your results output directory in the following files:\n",
    "\n",
    "    1. final_diff_exp_miRNAs.csv : This file contains the list of significantly differentially expressed miRNAs along with fold change, new names (in case of novel miRNAs), adjusted p-value, and their genomic locations.\n",
    "    \n",
    "    2. diff_exp_miRNAs_expression_counts.csv  : This file contains expression counts of significantly differentially expressed miRNAs.\n",
    "    \n",
    "    3. miRNA_expression_counts.csv : This file contains the raw expression counts of miRNAs in all samples before differential expression analysis.\n",
    "    \n",
    "    4. piRNA_raw_counts.csv :  This file contains raw counts of all piRNAs ontained from all samples before differential expression analysis.\n",
    "    \n",
    "    5. significantly_DE_piRNA.csv : This file contains the list  of all significantly differentially expressed piRNAs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "nbTranslate": {
   "displayLangs": [
    "*"
   ],
   "hotkey": "alt-t",
   "langInMainMenu": true,
   "sourceLang": "en",
   "targetLang": "fr",
   "useGoogleTranslate": true
  },
  "notify_time": "5",
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Open source RNA-Seq pipeline for identification of novel-mirs and their gene regulatory networks",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "181.767px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}