families_to_Ka_Ks_values.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Families to Ka/Ks values\n",
    "\n",
    "This notebook needs as input a file containing sequences ID and their respective family ID. Such file can be obtained withe the pipeline FTAG Finder v3 on Galaxy http://stat.genopole.cnrs.fr/galaxy.\n",
    "\n",
    "It will perform the Ka/Ks computation through the Paml programm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from Bio import Seq\n",
    "from Bio import SeqIO\n",
    "import subprocess\n",
    "import re\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_seq_prot = list(SeqIO.parse(\"source_data/Theobroma_cacao.Theobroma_cacao_20110822.pep.all.fa\", \"fasta\"))\n",
    "raw_seq_cds = list(SeqIO.parse(\"source_data/Theobroma_cacao.Theobroma_cacao_20110822.cds.all.fa\", \"fasta\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "ceb = pd.read_table(\"results_from_galaxy_FTAG_finder/families_filter30-70.tabular\", sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# get unique family\n",
    "families_set = set(ceb.family)\n",
    "gene_id_by_family = list()\n",
    "\n",
    "for family in families_set:\n",
    "    a = list(ceb[ceb.family == family][\"geneName\"])\n",
    "    gene_id_by_family.append(a)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In order to reduce the amount of computation only one couple of sequence per family are considered. \n",
    "gene_id_by_family = [ x[0:2] for x in gene_id_by_family]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check family length\n",
    "# [len(x) for x in gene_id_by_family]\n",
    "# seq_ID_to_write = list(ceb[ceb[\"family\"] == target_family][\"geneName\"])\n",
    "# seq_ID_to_write = list(ceb[0:500][\"geneName\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_ka = list()\n",
    "all_ks = list()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Main loop\n",
    "\n",
    "**First level**\n",
    "\n",
    "First loop run over the list *gene_id_by_family*. Each element of the list (here called *IDs* ) is also a list containing sequences ID for a given family. \n",
    "In each run of the loop, the focus is put on one family, whose raw pep and cds sequences are written in target_pep/cds.fa files. \n",
    "\n",
    "**Second level**\n",
    "\n",
    "Fasta files created in the first loop are read and each sequence inside is written one at a time in a tmp fasta file.\n",
    "Of course this step is useless in the case where only one couple of sequence per family is used. However, it will be needed if the full Ka/Ks analysis need to be done.\n",
    "\n",
    "\n",
    "The Paml script is then launched through a bash script. Than tmp files are erased in order not to crowd current folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "for IDs in gene_id_by_family:\n",
    "    with open(\"target_pep.fa\", \"w\") as OUT_pep:\n",
    "        for seq in raw_seq_prot:\n",
    "            if seq.id in IDs:\n",
    "                SeqIO.write(seq, OUT_pep, \"fasta\")\n",
    "\n",
    "    with open(\"target_cds.fa\", \"w\") as OUT_cds:\n",
    "        for seq in raw_seq_cds:\n",
    "            if seq.id in IDs:\n",
    "                SeqIO.write(seq, OUT_cds, \"fasta\")\n",
    "        \n",
    "    pep_list = list(SeqIO.parse(\"target_pep.fa\", \"fasta\"))\n",
    "    cds_list = list(SeqIO.parse(\"target_cds.fa\", \"fasta\"))\n",
    "    nb_seq = len(pep_list)\n",
    "    \n",
    "    for i in range(0, nb_seq-1):\n",
    "        for j in range(i+1, nb_seq):\n",
    "            with open(\"tmp_pep.fa\", \"w\") as OUT_pep:\n",
    "                SeqIO.write(pep_list[i], OUT_pep, \"fasta\")\n",
    "                SeqIO.write(pep_list[j], OUT_pep, \"fasta\")\n",
    "            with open(\"tmp_cds.fa\", \"w\") as OUT_cds:\n",
    "                SeqIO.write(cds_list[i], OUT_cds, \"fasta\")\n",
    "                SeqIO.write(cds_list[j], OUT_cds, \"fasta\")\n",
    "\n",
    "            cmd = \"./paml_run.sh tmp_pep.fa tmp_cds.fa\"\n",
    "            subprocess.run(cmd, stdout = subprocess.PIPE, stderr = subprocess.PIPE, shell = True)\n",
    "\n",
    "            # Line with Ks and Ka values in yn result file is the only one starting with \n",
    "            res_line = \"\".join([line for line in open('yn') if re.match('^\\s+[0-9]',line)])\n",
    "            res_line = re.split(\"\\s+\", res_line)\n",
    "            all_ka.append(res_line[8])\n",
    "            all_ks.append(res_line[11])\n",
    "            subprocess.run(\"rm 2YN.* cds.ali.phy prot.ali.aln res rst rst1 rub tmp_pep.dnd yn tmp_cds.fa tmp_pep.fa \", \n",
    "                           stdout = subprocess.PIPE, stderr = subprocess.PIPE, shell = True)\n",
    "    subprocess.run(\"rm target_pep.fa target_cds.fa \", stdout = subprocess.PIPE, stderr = subprocess.PIPE, shell = True)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save all output (Ka, Ks and ratio of both) in a csv file\n",
    "OUT_ka_ks = open(\"ka_ks_values\", \"w\")\n",
    "OUT_ka_ks.write(\"ka;ks;ratio\\n\")\n",
    "for i in range(len(all_ka)):\n",
    "        if float(all_ks[i]) == 0:\n",
    "            ratio = 0\n",
    "        else :\n",
    "            ratio = float(all_ka[i]) / float(all_ks[i])\n",
    "        OUT_ka_ks.write(str(all_ka[i]) + \";\" + str(all_ks[i]) + \";\" + str(ratio) + \"\\n\" )\n",
    "OUT_ka_ks.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}