{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "## Setting up Environment"
      ],
      "metadata": {
        "id": "EmRILbC1Q_CV"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!apt-get install -y clustalo\n",
        "!apt-get install -y paml"
      ],
      "metadata": {
        "id": "t2a5PjSrB6OB"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install biopython"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "yd-MlGPHLNpe",
        "outputId": "6270c546-7210-4fa4-a8d8-ff9888c81a2e"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Collecting biopython\n",
            "  Downloading biopython-1.81-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (3.1 MB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m3.1/3.1 MB\u001b[0m \u001b[31m12.6 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hRequirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (from biopython) (1.22.4)\n",
            "Installing collected packages: biopython\n",
            "Successfully installed biopython-1.81\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ySEO3ONEZpLi"
      },
      "outputs": [],
      "source": [
        "import os\n",
        "import subprocess\n",
        "\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "from Bio.Align.Applications import ClustalOmegaCommandline\n",
        "from Bio import SeqIO, Entrez\n",
        "from google.colab import drive"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "Entrez.email = ''"
      ],
      "metadata": {
        "id": "LKMazW8QYFaJ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "drive.mount('/content/drive/')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "LaAndAfuaBDy",
        "outputId": "7e74c44e-1dda-4298-f2c9-733fa35d787c"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Mounted at /content/drive/\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "cd /content/drive/MyDrive/Reading_genes_evolution/"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "QjpTLp4xaFFx",
        "outputId": "296169f5-5ce8-4a6a-ed55-db7ae756a38d"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "/content/drive/MyDrive/Reading_genes_evolution\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Create output folders if necessary\n",
        "if not os.path.exists('transcripts'):\n",
        "    os.makedirs('transcripts')\n",
        "\n",
        "if not os.path.exists('proteins'):\n",
        "    os.makedirs('proteins')\n",
        "\n",
        "if not os.path.exists('proteins_alignment'):\n",
        "    os.makedirs('proteins_alignment')\n",
        "\n",
        "if not os.path.exists('cds_alignment'):\n",
        "    os.makedirs('cds_alignment')\n",
        "\n",
        "if not os.path.exists('codeml_data'):\n",
        "    os.makedirs('codeml_data')"
      ],
      "metadata": {
        "id": "4OAceQIZ0iIi"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Path to the proteins folder\n",
        "proteins_dir = '/content/drive/MyDrive/Reading_genes_evolution/proteins'\n",
        "transcripts_dir = '/content/drive/MyDrive/Reading_genes_evolution/transcripts/'\n",
        "\n",
        "protein_aln_dir = '/content/drive/MyDrive/Reading_genes_evolution/proteins_alignment'\n",
        "cds_aln_dir = '/content/drive/MyDrive/Reading_genes_evolution/cds_alignment/'\n",
        "\n",
        "# Define paths to input and output directories for pal2nal\n",
        "codeml_dir = '/content/drive/MyDrive/Reading_genes_evolution/codeml_data/'"
      ],
      "metadata": {
        "id": "4CSyGvMDKNCz"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Analysis of Signals of Positive Selection"
      ],
      "metadata": {
        "id": "NymuR-2ARvaR"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "Hh-KK5DEF6yD"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "## Run Clustal Omega to produce protein alighments\n",
        "## Do not forget to install Clustal Omega in google collab\n",
        "\n",
        "# Loop over all files in the proteins folder\n",
        "for protein_file in os.listdir(proteins_dir):\n",
        "    # Check that the file is a FASTA file\n",
        "    if protein_file.endswith('.fasta'):\n",
        "        # Build the input and output file paths\n",
        "        protein_path = os.path.join(proteins_dir, protein_file)\n",
        "        protein_name = os.path.splitext(protein_file)[0]\n",
        "        out_path = os.path.join(protein_aln_dir, protein_name + '_aln.fasta')\n",
        "\n",
        "        # Build the Clustal Omega commandline\n",
        "        clustalo_cline = ClustalOmegaCommandline(infile=protein_path, outfile=out_path, verbose=True, auto=True)\n",
        "\n",
        "        # Run Clustal Omega\n",
        "        stdout, stderr = clustalo_cline()\n",
        "\n",
        "        # Print the output and error messages\n",
        "        print(stdout)\n",
        "        print(stderr)"
      ],
      "metadata": {
        "id": "RgzslxAB_3Cr"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "## Run pal2nal to produce codon alighments for all genes\n",
        "\n",
        "# Get a list of all files in the protein alignment directory\n",
        "protein_files = [f for f in os.listdir(protein_aln_dir) if os.path.isfile(os.path.join(protein_aln_dir, f))]\n",
        "\n",
        "# Process each protein alignment file\n",
        "for protein_file in protein_files:\n",
        "    if not protein_file.endswith('_protein_aln.fasta'):\n",
        "        continue\n",
        "\n",
        "    # Construct the path to the corresponding transcript file\n",
        "    transcript_file = os.path.join(transcripts_dir, protein_file.replace('_protein_aln.fasta', '_cds.fasta'))\n",
        "\n",
        "    # Construct the path to the output file\n",
        "    output_file = os.path.join(cds_aln_dir, protein_file.replace('_protein_aln.fasta', '_cds_aln.paml'))\n",
        "\n",
        "    # Run the pal2nal command\n",
        "    command = f'perl ./pal2nal.v14/pal2nal.pl {os.path.join(protein_aln_dir, protein_file)} {transcript_file} -nogap -output paml > {output_file}'\n",
        "    subprocess.run(command, shell=True)\n",
        "    #print(command)"
      ],
      "metadata": {
        "id": "u2raE1jaJvqZ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "## Prepare codeml.ctl files for each gene\n",
        "\n",
        "# Loop through all files in the input directory\n",
        "for filename in os.listdir(cds_aln_dir):\n",
        "    # Skip non-fasta files\n",
        "    if not filename.endswith('.paml'):\n",
        "        continue\n",
        "\n",
        "    # Extract the gene symbol from the filename\n",
        "    gene_symbol = filename.split('_')[0]\n",
        "\n",
        "    # Construct the path to the output directory for this gene\n",
        "    gene_output_dir = os.path.join(codeml_dir, gene_symbol)\n",
        "    os.makedirs(gene_output_dir, exist_ok=True)\n",
        "\n",
        "    # Construct the path to the input alignment file\n",
        "    input_cds_aln = os.path.join(cds_aln_dir, filename)\n",
        "\n",
        "    # Construct the path to the output codeml control file\n",
        "    output_codeml = os.path.join(gene_output_dir, 'codeml.ctl')\n",
        "\n",
        "    # Create the codeml control file\n",
        "    codeml_ctl = f\"seqfile = {input_cds_aln}\\n\"\n",
        "    codeml_ctl += \"treefile = /content/drive/MyDrive/Reading_genes_evolution/great_apes.tree\\n\"\n",
        "    codeml_ctl += f\"outfile = {os.path.join(gene_output_dir, filename[:-5])}.out\\n\"\n",
        "    codeml_ctl += \"verbose = 1\\n\"\n",
        "    codeml_ctl += \"runmode = 0\\n\"\n",
        "    codeml_ctl += \"seqtype = 1\\n\"\n",
        "    codeml_ctl += \"CodonFreq = 0\\n\"\n",
        "    codeml_ctl += \"model = 0\\n\"\n",
        "    codeml_ctl += \"NSsites = 0 1 2 7 8\\n\"\n",
        "    codeml_ctl += \"icode = 0\\n\"\n",
        "    codeml_ctl += \"RateAncestor = 0\\n\"\n",
        "    codeml_ctl += \"clock = 0\\n\"\n",
        "    codeml_ctl += \"\\n\"\n",
        "    codeml_ctl += \"method = 0\\n\"\n",
        "    codeml_ctl += \"ndata = 1\\n\"\n",
        "\n",
        "    # Write the codeml control file\n",
        "    with open(output_codeml, 'w') as f:\n",
        "        f.write(codeml_ctl)"
      ],
      "metadata": {
        "id": "TqTlSMyKl18J"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "## Run codeml for each gene\n",
        "\n",
        "# iterate through each gene folder in the directory\n",
        "for gene_folder in os.listdir(codeml_dir):\n",
        "    # check if the current item is a directory\n",
        "    if os.path.isdir(os.path.join(codeml_dir, gene_folder)):\n",
        "        # change directory to the gene folder\n",
        "        os.chdir(os.path.join(codeml_dir, gene_folder))\n",
        "\n",
        "        # run the codeml command in the gene folder\n",
        "        print(gene_folder)\n",
        "        subprocess.run(['/usr/bin/codeml'])\n",
        "\n",
        "        # change back to the original directory\n",
        "        os.chdir(\"..\")"
      ],
      "metadata": {
        "id": "-6D3DKQe9ZnI"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Summary of CODEML results"
      ],
      "metadata": {
        "id": "VpPG9CIDSRja"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!grep --with-filename \"lnL\" ./codeml_data/*/*out | awk '{print $1,$4,$5}' | sed 's/)://g' | tr \" \" \"\\t\" | awk '{ if (NR % 5 == 0) {print $0\"#\"} else {print}}' | tr \"\\n\" \"\\t\" | sed 's/#/\\n/g' | awk '{print $1,$2,$3,$5,$6,$8,$9,$11,$12,$14,$15}' > codeml_data_output.txt"
      ],
      "metadata": {
        "id": "PWWf76Mf0zPT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#M0 vs M1 – Support for neutral selection can be identified if M1 provides a better fit than M0\n",
        "!grep --with-filename \"lnL\"  ./codeml_data/*/*out |awk '{print $1,$4,$5}' |sed 's/)://g' |tr \" \" \"\\t\" |awk '{ if (NR % 5 == 0) {print $0\"#\"} else {print}}'  |tr \"\\n\" \"\\t\" |sed 's/#/\\n/g' |awk '{print $1,$2,$3,$5,$6,$8,$9,$11,$12,$14,$15}' |awk '{print  $1,($4-$2),(2*($5+(-1*$3))) }' |  tr \" \" \"\\t\" >M0M1ChiSquares.tab\n",
        "\n",
        "#M1 vs M2 – Support for positive selection can be identified if M2 provides a better fit than M1\n",
        "!grep --with-filename \"lnL\"  ./codeml_data/*/*out |awk '{print $1,$4,$5}' |sed 's/)://g' |tr \" \" \"\\t\" |awk '{ if (NR % 5 == 0) {print $0\"#\"} else {print}}'  |tr \"\\n\" \"\\t\" |sed 's/#/\\n/g' |awk '{print $1,$2,$3,$5,$6,$8,$9,$11,$12,$14,$15}' |awk '{if ($6>$4) {print  $1,($6-$4),(2*($7+(-1*$5)))}else {print \"FLIP\"$1,($4-$6),(2*($5+(-1*$7)))} }'  |  tr \" \" \"\\t\" >M1M2ChiSquares.tab\n",
        "\n",
        "#M7 vs M8 –The M8–M7 comparison offers a more stringent test for positive selection\n",
        "!grep --with-filename \"lnL\"  ./codeml_data/*/*out |awk '{print $1,$4,$5}' |sed 's/)://g' |tr \" \" \"\\t\" |awk '{ if (NR % 5 == 0) {print $0\"#\"} else {print}}'  |tr \"\\n\" \"\\t\" |sed 's/#/\\n/g' |awk '{print $1,$2,$3,$5,$6,$8,$9,$11,$12,$14,$15}' |awk '{if($10>$8){print  $1, ($10-$8), (2*($11+(-1*$9)))} else {print $1,($8-$10), (2*($9+(-1*$11)))} }' |  tr \" \" \"\\t\" >M7M8ChiSquares.tab"
      ],
      "metadata": {
        "id": "yDXLnWs5kkdZ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "## Tests for neutral selection comparing models M0 and M1a\n",
        "!less M0M1ChiSquares.tab |awk '{print $2}' |while read line; do grep -w \"$line\" chi2lists.txt  ;done |paste M0M1ChiSquares.tab - |awk '$3>$5' |wc\n",
        "\n",
        "!less M0M1ChiSquares.tab |awk '{print $2}' |while read line; do grep -w \"$line\" chi2lists.txt  ;done |paste M0M1ChiSquares.tab - |awk '$3>$5' |sed 's/:/\\t/g' |cut -f 1 |sed 's/FLIP//g' |while read line; do grep --with-filename -A 100 \"Model 1\" $line |grep \"\\*\" - |grep -v \"branch\" - |grep -v \"Positively\" -;done > NeutralSelectedSitesM1.txt\n",
        "\n",
        "!less NeutralSelectedSitesM1.txt |awk '{print $1}' |sort|uniq|wc"
      ],
      "metadata": {
        "id": "fpghOsImF8Ue",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "0966e4aa-f1f4-45e3-ddc0-22f9ee51f737"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "     23     115    1745\n",
            "     22      22     880\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "## Positive selection tests comparing models M1a and M2a\n",
        "\n",
        "!less M1M2ChiSquares.tab |awk '{print $2}' |while read line; do grep -w \"$line\" chi2lists.txt  ;done |paste M1M2ChiSquares.tab - |awk '$3>$5' |wc\n",
        "\n",
        "!less M1M2ChiSquares.tab |awk '{print $2}' |while read line; do grep -w \"$line\" chi2lists.txt  ;done |paste M1M2ChiSquares.tab - |awk '$3>$5' |sed 's/:/\\t/g' |cut -f 1 |sed 's/FLIP//g' |while read line; do grep --with-filename -A 100 \"Model 2\" $line |grep \"\\*\" - |grep -v \"branch\" - |grep -v \"Positively\" -;done > PositivelySelectedSitesM2.txt\n",
        "\n",
        "# Number of orthogroups with significant log ratio threshold and Naive Empirical Bayes (NEB) analysis\n",
        "!less PositivelySelectedSitesM2.txt |awk '{print $1}' |sort |uniq|wc"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "jeJCoLwLCy0T",
        "outputId": "c331e3d1-f080-4470-a1d9-e5043b90966b"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "     17      83    1210\n",
            "      8       8     312\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "gbC5LRGJC9Nu"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "## Positive selection tests comparing models M7 and M8 – a more stringent site model\n",
        "\n",
        "!less M7M8ChiSquares.tab |awk '{print $2}' |while read line; do grep -w \"$line\" chi2lists.txt  ;done |paste M7M8ChiSquares.tab - |awk '$3>$5' |sed 's/:/\\t/g' |cut -f 1 |wc\n",
        "\n",
        "!less M7M8ChiSquares.tab |awk '{print $2}' |while read line; do grep -w \"$line\" chi2lists.txt  ;done |paste M7M8ChiSquares.tab - |awk '$3>$5' |sed 's/:/\\t/g' |cut -f 1 |while read line; do grep --with-filename -A 500 \"Model 8\" $line |grep \"\\*\" - |grep -v \"branch\" - |grep -v \"Positively\" - ;done >PositivelySelectedSitesM8.txt\n",
        "\n",
        "# Number of orthgroups with a significant log ratio threshold and Naive Emperical Bayes (NEB) analysis\n",
        "!less PositivelySelectedSitesM8.txt |awk '{print $1}' |sort|uniq|wc\n",
        "\n",
        "!wc PositivelySelectedSitesM8.txt"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "LhAW4hHeEwYT",
        "outputId": "56b224b0-0e86-4244-f823-f3ac1482a9d1"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "      8       8     302\n",
            "      7       7     273\n",
            "  289  1731 22309 PositivelySelectedSitesM8.txt\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "## How many orthogroups have positive selection under the M0 model – single dN/dS value\n",
        "!grep \"omega\" ./codeml_data/*/*out |awk '$4>1 && $4<998' >PositivelySelectedOrthogroupsM0.txt\n",
        "!wc PositivelySelectedOrthogroupsM0.txt"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "rQWKwfMtFpFv",
        "outputId": "2901a5e1-a783-44d2-a219-5346c94f7cbe"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0 0 0 PositivelySelectedOrthogroupsM0.txt\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "!cat PositivelySelected* |sed 's/_/\\t/g' |awk '{print $1}' |sort |uniq |sed 's/$/_tree.txt/g' |tr \"\\n\" \" \" |sed 's/^/cat /g' |sed 's/$/ >AllPositiveTress.txt/g' >concatenatePositiveSelectionTrees.sh\n",
        "!sh concatenatePositiveSelectionTrees.sh"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "hgvv16d2F7Jf",
        "outputId": "db9bd5a2-8c4b-4877-d8c0-29c6e75f023c"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "cat: ./codeml_tree.txt: No such file or directory\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "!cat PositivelySelected* |sed 's/_/\\t/g' |awk '{print $2}' |sort |uniq |sed 's/$/_tree.txt/g' |tr \"\\n\" \" \" |sed 's/^/cat /g' |sed 's/$/ >AllPositiveTress.txt/g'"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "BzKIw9pTG1vo",
        "outputId": "d2f1f923-84ee-46b2-e994-47033276903b"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "cat data/CEP85L/CEP85L_tree.txt data/ESRRG/ESRRG_tree.txt data/HTT/HTT_tree.txt data/MITF/MITF_tree.txt data/MRPL19/MRPL19_tree.txt data/RBFOX2/RBFOX2_tree.txt data/RFC5/RFC5_tree.txt data/SCAF11/SCAF11_tree.txt  >AllPositiveTress.txt"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "!cat PositivelySelected* |sed 's/_/\\t/g' |awk '{print $2}' |sort |uniq"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Dy9Xy_fzHS2M",
        "outputId": "53a1aef3-200f-4cf9-b42c-face5afa8b1e"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "data/CEP85L/CEP85L\n",
            "data/ESRRG/ESRRG\n",
            "data/HTT/HTT\n",
            "data/MITF/MITF\n",
            "data/MRPL19/MRPL19\n",
            "data/RBFOX2/RBFOX2\n",
            "data/RFC5/RFC5\n",
            "data/SCAF11/SCAF11\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Reading genes evolution"
      ],
      "metadata": {
        "id": "LOw-ZBsZ9Tvj"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "reading_genes = []\n",
        "with open('Reading_genes.txt', 'r') as rg:\n",
        "  for g in rg:\n",
        "    #print(g.strip())\n",
        "    reading_genes.append(g.strip())"
      ],
      "metadata": {
        "id": "5AEobStJQpF4"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "homologene_file = './homologene_homo_sapiens_02.05.2023.txt'\n",
        "\n",
        "xs = []"
      ],
      "metadata": {
        "id": "hLCpe_VabY4-"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "with open(homologene_file, 'r') as hf:\n",
        "  for x in hf:\n",
        "    xs.append(x.strip('\\n'))"
      ],
      "metadata": {
        "id": "OnwqdAJfbiSG"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "xs"
      ],
      "metadata": {
        "id": "Q_enwP8AdYrP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "ys[:20]"
      ],
      "metadata": {
        "id": "b6yQi_mFdaZW"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "genes_age_dict = {}\n",
        "taxon_counts_global = {}\n",
        "taxon_counts_reading = {}\n",
        "\n",
        "for x in xs:\n",
        "  if 'HomoloGene' in x:\n",
        "    hmg_id = x.split()[0].strip(':')\n",
        "    status = x.split()[3]\n",
        "    taxon = x.split()[-1]\n",
        "    #print(hmg_id, status, taxon, x)\n",
        "    if taxon in taxon_counts_global:\n",
        "      taxon_counts_global[taxon] += 1\n",
        "    else:\n",
        "      taxon_counts_global[taxon] = 1\n",
        "  else:\n",
        "    ncbi_id, description, species = x.split('\\t')\n",
        "    if species == 'Homo sapiens' and ncbi_id in reading_genes:\n",
        "      #print(ncbi_id, taxon)\n",
        "      if taxon in taxon_counts_reading:\n",
        "        taxon_counts_reading[taxon] += 1\n",
        "      else:\n",
        "        taxon_counts_reading[taxon] = 1"
      ],
      "metadata": {
        "id": "tU_SxIAydngP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "taxon_counts_global"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "8cVz8WJzdraX",
        "outputId": "214153f5-aeed-49ad-ace3-5f01e4af3736"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "{'Homininae': 465,\n",
              " 'H.sapiens': 52,\n",
              " 'Euarchontoglires': 280,\n",
              " 'Boreoeutheria': 2856,\n",
              " 'Euteleostomi': 6612,\n",
              " 'Tetrapoda': 1470,\n",
              " 'Sordariomycetidae': 4,\n",
              " 'Opisthokonta': 576,\n",
              " 'Catarrhini': 422,\n",
              " 'Bilateria': 2882,\n",
              " 'Amniota': 627,\n",
              " 'Eukaryota': 2476,\n",
              " 'M.oryzae': 1,\n",
              " 'A.thaliana': 1,\n",
              " 'saccharomyceta': 1,\n",
              " 'Diptera': 2,\n",
              " 'Ascomycota': 1,\n",
              " 'Mesangiospermae': 4}"
            ]
          },
          "metadata": {},
          "execution_count": 21
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "taxon_counts_reading"
      ],
      "metadata": {
        "id": "b3Iqq9TLQ1nn",
        "outputId": "f5ee376c-9c3e-4774-8251-7fe5306951d2",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "{'Boreoeutheria': 7,\n",
              " 'Amniota': 7,\n",
              " 'Homininae': 1,\n",
              " 'Euteleostomi': 76,\n",
              " 'Catarrhini': 1,\n",
              " 'Bilateria': 32,\n",
              " 'Eukaryota': 19,\n",
              " 'Opisthokonta': 5,\n",
              " 'Tetrapoda': 13}"
            ]
          },
          "metadata": {},
          "execution_count": 22
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "taxon_counts_reading"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "BUSQsjrwgC0Q",
        "outputId": "9de3d655-0111-410f-a451-18510701c3fe"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "{'Boreoeutheria': 5,\n",
              " 'Catarrhini': 1,\n",
              " 'Bilateria': 16,\n",
              " 'Eukaryota': 13,\n",
              " 'Opisthokonta': 3,\n",
              " 'Amniota': 4,\n",
              " 'Euteleostomi': 50,\n",
              " 'Tetrapoda': 8}"
            ]
          },
          "metadata": {},
          "execution_count": 15
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "ys = []\n",
        "homology_group_linear = []\n",
        "for x in xs:\n",
        "  if 'HomoloGene' in x:\n",
        "    if homology_group_linear:\n",
        "      ys.append(homology_group_linear)\n",
        "    homology_group_linear = []\n",
        "    homology_group_linear.append(x)\n",
        "  else:\n",
        "    homology_group_linear.append(x)\n",
        "ys.append(homology_group_linear)"
      ],
      "metadata": {
        "id": "zsy7QjXHO6gt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "total_list = {}\n",
        "for y in ys:\n",
        "  taxon = y[0].split(' ')[-1]\n",
        "  gene_list_in_taxon = [x.split('\\t')[0] for x in y[1:]]\n",
        "  #print(y[0].split(' ')[-1], [x.split('\\t')[0] for x in y[1:]])\n",
        "  for gene in gene_list_in_taxon:\n",
        "    if gene in reading_genes:\n",
        "      total_list[gene] = taxon\n",
        "      #print(taxon, gene)"
      ],
      "metadata": {
        "id": "9CnSxZStRddw"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "len(ys)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "-2-rZzap8oMX",
        "outputId": "dfbb5d47-3f8d-4f81-d6ad-e3bc704b0195"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "18732"
            ]
          },
          "metadata": {},
          "execution_count": 66
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "xs[:10]"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "p_GSgu4-PVeZ",
        "outputId": "54d37ecf-cd02-4724-bfc5-8556e2b7d699"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "['1: HomoloGene: 141252.Gene conserved in Homininae',\n",
              " 'LOC646021\\tuncharacterized LOC646021\\tHomo sapiens',\n",
              " 'LOC454564\\tuncharacterized LOC454564\\tPan troglodytes',\n",
              " '2: HomoloGene: 141220.Gene exclusive to H.sapiens',\n",
              " 'LOC102724405\\timmunoglobulin superfamily member 3-like\\tHomo sapiens',\n",
              " 'LOC102724844\\tuncharacterized LOC102724844\\tHomo sapiens',\n",
              " 'LOC102724123\\timmunoglobulin superfamily member 3-like\\tHomo sapiens',\n",
              " '3: HomoloGene: 141219.Gene conserved in Homininae',\n",
              " 'LOC102723552\\tankyrin repeat domain-containing protein 20A3-like\\tHomo sapiens',\n",
              " 'LOC474068\\tputative ankyrin repeat domain-containing protein 20A5-like\\tPan troglodytes']"
            ]
          },
          "metadata": {},
          "execution_count": 21
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "GjhD4e4uZRFf"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "f9dYbHPYoA4y"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}