counting.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def necessary_information(filename):\n",
    "    lengths = np.empty(0)\n",
    "    longest = ''\n",
    "    current = ''\n",
    "    with open(filename) as fd:\n",
    "        lines = fd.readlines()\n",
    "    for line in lines:\n",
    "        if line[0] == '>':\n",
    "            if (lengths.size > 0) and (lengths[-1] > len(longest)):\n",
    "                longest = current\n",
    "            lengths = np.append(lengths, 0)\n",
    "            current = ''\n",
    "        else:\n",
    "            line = line.replace(' ', '')\n",
    "            line = line.replace('\\n', '')\n",
    "            current += line\n",
    "            lengths[-1] += len(line)\n",
    "    if (lengths.size > 0) and (lengths[-1] > len(longest)):\n",
    "        longest = current\n",
    "    return longest, lengths\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_stats(lengths, name):\n",
    "    lengths = np.sort(lengths)[::-1]\n",
    "    sum = lengths.sum()\n",
    "    print(f'Statistics for {name}:')\n",
    "    print(f'Total number: {lengths.size}')\n",
    "    print(f'Total length: {sum}')\n",
    "    print(f'Length of the longest: {lengths[0]}')\n",
    "    prefices_sum = np.cumsum(lengths)\n",
    "    print(f'N50: {lengths[np.argmax(prefices_sum >= sum / 2)]}')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_gap_stats(sequence, name):\n",
    "    total_number = sequence.count('AN') + sequence.count('TN') + sequence.count('GN') + sequence.count('CN') \n",
    "    total_length = sequence.count('N')\n",
    "    print(f'Gap statistics for {name}:')\n",
    "    print(f'Total number: {total_number}')\n",
    "    print(f'Total length of gaps: {total_length}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Statistics for contigs:\n",
      "Total number: 594\n",
      "Total length: 3922956.0\n",
      "Length of the longest: 179304.0\n",
      "N50: 48055.0\n",
      "\n",
      "Statistics for scaffolds:\n",
      "Total number: 67\n",
      "Total length: 3874616.0\n",
      "Length of the longest: 3834998.0\n",
      "N50: 3834998.0\n",
      "\n",
      "Gap statistics for longest scaffold before gapclose:\n",
      "Total number: 61\n",
      "Total length of gaps: 7380\n",
      "\n",
      "Gap statistics for longest scaffold after gapclose:\n",
      "Total number: 7\n",
      "Total length of gaps: 2037\n"
     ]
    }
   ],
   "source": [
    "longest_contig, contig_length = necessary_information('contig/out_contig.fa')\n",
    "longest_scaffold, scaffold_length = necessary_information('scaffold/Poil_scaffold.fa')\n",
    "longest_scaffold_gap, scaffold_gap_length = necessary_information('gapclose/Poil_gapClosed.fa')\n",
    "\n",
    "#print(len(longest_contig))\n",
    "#print(len(longest_scaffold))\n",
    "#print(len(longest_scaffold_gap))\n",
    "\n",
    "print_stats(contig_length, 'contigs')\n",
    "print()\n",
    "print_stats(scaffold_length, 'scaffolds')\n",
    "print()\n",
    "print_gap_stats(longest_scaffold, 'longest scaffold before gapclose')\n",
    "print()\n",
    "print_gap_stats(longest_scaffold_gap, 'longest scaffold after gapclose')"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
  },
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "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.8.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}