{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get read counts of fastq files\n",
    "\n",
    "In this notebook we are going to count the number of reads assigned to each sample. To do that we are going to create a dataframe that contains all the information.\n",
    "\n",
    "For ease purposes, we are going to create two dataframes: the first one related to the reads generated during the preprocessing step; and the second with the reads mapped to the different profilers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import subprocess\n",
    "import os\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "sys.path.append('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from list_vars import LIST_PROFILERS, DATA_DIR, RESULTS_DIR, POOLS, CONTROLS, DPI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.makedirs(f'{RESULTS_DIR}/counts', exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting counts of fastq files (host mapping)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def count_reads_in_fastq(file_path):\n",
    "    # Path to the Bash script\n",
    "    script_path = \"./sh_funcs/count_reads.sh\"\n",
    "\n",
    "    # Call the script with the file path as an argument\n",
    "    try:\n",
    "        result = subprocess.run(\n",
    "            [script_path, file_path], \n",
    "            text=True, \n",
    "            capture_output=True, \n",
    "            check=True\n",
    "        )\n",
    "        # The script's output is the number of reads\n",
    "        return int(result.stdout.strip())\n",
    "    except subprocess.CalledProcessError as e:\n",
    "        print(f\"Error occurred: {e.stderr}\")\n",
    "        return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(f'{RESULTS_DIR}/counts/mapping_counts.txt', 'w') as file:\n",
    "    file.write(f'SAMPLE\\traw\\t1st_unmapped\\t1st_mapped\\t2nd_unmapped\\t2nd_mapped\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "dict_dict_reads = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dict_reads_ARTIFICIAL = {}\n",
    "\n",
    "dict_reads_ARTIFICIAL['raw'] = count_reads_in_fastq(f'{DATA_DIR}/artificial_v2/artificial_reads_R1.fastq.gz')\n",
    "dict_reads_ARTIFICIAL['1st_unmapped'] = count_reads_in_fastq(f'{RESULTS_DIR}/1stmap/artificial/ARTIFICIAL.unmapped_1.fastq.gz')\n",
    "dict_reads_ARTIFICIAL['1st_mapped'] = dict_reads_ARTIFICIAL['raw'] - dict_reads_ARTIFICIAL['1st_unmapped']\n",
    "dict_reads_ARTIFICIAL['2nd_unmapped'] = count_reads_in_fastq(f'{RESULTS_DIR}/2ndmap/artificial/ARTIFICIAL.unmapped.fastq.1.gz')\n",
    "dict_reads_ARTIFICIAL['2nd_mapped'] = dict_reads_ARTIFICIAL['1st_unmapped'] - dict_reads_ARTIFICIAL['2nd_unmapped']\n",
    "\n",
    "with open(f'{RESULTS_DIR}/counts/mapping_counts.txt', 'a') as file:\n",
    "    file.write(f\"ARTIFICIAL\\t{dict_reads_ARTIFICIAL['raw']}\\t{dict_reads_ARTIFICIAL['1st_unmapped']}\\t{dict_reads_ARTIFICIAL['1st_mapped']}\\t{dict_reads_ARTIFICIAL['2nd_unmapped']}\\t{dict_reads_ARTIFICIAL['2nd_mapped']}\\n\")\n",
    "\n",
    "dict_dict_reads['ARTIFICIAL'] = dict_reads_ARTIFICIAL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for control in CONTROLS:\n",
    "    dict_reads_CONTROLS_i = {}\n",
    "\n",
    "    dict_reads_CONTROLS_i['raw'] = count_reads_in_fastq(f'{DATA_DIR}/EM_EVPools/control_sample/{control}.fastq.gz')\n",
    "    dict_reads_CONTROLS_i['1st_unmapped'] = count_reads_in_fastq(f'{RESULTS_DIR}/1stmap/controls/{control}.unmapped_1.fastq.gz')\n",
    "    dict_reads_CONTROLS_i['1st_mapped'] = dict_reads_CONTROLS_i['raw'] - dict_reads_CONTROLS_i['1st_unmapped']\n",
    "    dict_reads_CONTROLS_i['2nd_unmapped'] = count_reads_in_fastq(f'{RESULTS_DIR}/2ndmap/controls/{control}.unmapped.fastq.1.gz')\n",
    "    dict_reads_CONTROLS_i['2nd_mapped'] = dict_reads_CONTROLS_i['1st_unmapped'] - dict_reads_CONTROLS_i['2nd_unmapped']\n",
    "\n",
    "    with open(f'{RESULTS_DIR}/counts/mapping_counts.txt', 'a') as file:\n",
    "        file.write(f\"{control}\\t{dict_reads_CONTROLS_i['raw']}\\t{dict_reads_CONTROLS_i['1st_unmapped']}\\t{dict_reads_CONTROLS_i['1st_mapped']}\\t{dict_reads_CONTROLS_i['2nd_unmapped']}\\t{dict_reads_CONTROLS_i['2nd_mapped']}\\n\")  \n",
    "\n",
    "    dict_dict_reads[control] = dict_reads_CONTROLS_i  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for pool in POOLS:\n",
    "    dict_reads_POOL_i = {}\n",
    "\n",
    "    dict_reads_POOL_i['raw'] = count_reads_in_fastq(f'{DATA_DIR}/EM_EVPools/control_sample/{pool}.fastq.gz')\n",
    "    dict_reads_POOL_i['1st_unmapped'] = count_reads_in_fastq(f'{RESULTS_DIR}/1stmap/pools/{pool}.unmapped_1.fastq.gz')\n",
    "    dict_reads_POOL_i['1st_mapped'] = dict_reads_POOL_i['raw'] - dict_reads_POOL_i['1st_unmapped']\n",
    "    dict_reads_POOL_i['2nd_unmapped'] = count_reads_in_fastq(f'{RESULTS_DIR}/2ndmap/pools/{pool}.unmapped.fastq.1.gz')\n",
    "    dict_reads_POOL_i['2nd_mapped'] = dict_reads_POOL_i['1st_unmapped'] - dict_reads_POOL_i['2nd_unmapped']\n",
    "\n",
    "    with open(f'{RESULTS_DIR}/counts/mapping_counts.txt', 'a') as file:\n",
    "        file.write(f\"{pool}\\t{dict_reads_POOL_i['raw']}\\t{dict_reads_POOL_i['1st_unmapped']}\\t{dict_reads_POOL_i['1st_mapped']}\\t{dict_reads_POOL_i['2nd_unmapped']}\\t{dict_reads_POOL_i['2nd_mapped']}\\n\")  \n",
    "\n",
    "    dict_dict_reads[pool] = dict_reads_POOL_i                                                                                             "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get counts of profiling results\n",
    "\n",
    "In this case we are going to build a dataframe with the following columns:\n",
    "- SAMPLE\n",
    "- PASS\n",
    "- MODE\n",
    "- PROFILER\n",
    "\n",
    "And for each profiler we are going to read the standarized report adn extract the number of reads mapped in total, mapped to homo sapiens, and lastly extract the reads that were unmapped."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_host_map_info = pd.read_csv(f'{RESULTS_DIR}/counts/mapping_counts.txt', sep='\\t').set_index('SAMPLE')\n",
    "df_host_map_info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_read_info(taxpasta_file, basename, passn):\n",
    "    df = pd.read_csv(taxpasta_file, sep='\\t').set_index('taxonomy_id')\n",
    "\n",
    "    if 9606 in df.index:\n",
    "        counts_human = int(df.loc[9606, 'count'])\n",
    "    else:\n",
    "        counts_human = 0\n",
    "\n",
    "    counts_others = int(df['count'].sum()) - counts_human\n",
    "    if passn == 0:\n",
    "        max_counts = df_host_map_info.loc[basename, 'raw']\n",
    "    elif passn == 2:\n",
    "        max_counts = df_host_map_info.loc[basename, '2nd_unmapped']\n",
    "\n",
    "    counts_unmapped = int(max_counts) - (counts_human + counts_others)\n",
    "\n",
    "    return counts_human, counts_others, counts_unmapped"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_human, counts_others, counts_unmapped = extract_read_info(f'{RESULTS_DIR}/profiling/kraken2/pass0/ARTIFICIAL_mode3/ARTIFICIAL_mode3.report.standardised.species', 'ARTIFICIAL', 0)\n",
    "counts_human, counts_others, counts_unmapped"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = 'ARTIFICIAL'\n",
    "with open(f'{RESULTS_DIR}/counts/profiling_counts_{sample}.txt', 'w') as file:\n",
    "        file.write(f'SAMPLE\\tpass\\tmode\\tprofiler\\tmapped_human\\tmapped_others\\tunmapped\\n')  \n",
    "\n",
    "for profiler in LIST_PROFILERS:\n",
    "    for passn in [0, 2]:\n",
    "        for mode in range(1, 10):\n",
    "            if profiler in ['kaiju']:\n",
    "                taxpasta_file = f'{RESULTS_DIR}/profiling/{profiler}/pass{passn}/{sample}_mode{mode}/{sample}_mode{mode}.results.standardised.species'\n",
    "            elif profiler in ['krakenuniq', 'kraken2']:\n",
    "                taxpasta_file = f'{RESULTS_DIR}/profiling/{profiler}/pass{passn}/{sample}_mode{mode}/{sample}_mode{mode}.report.standardised.species'\n",
    "            elif profiler in ['centrifuge']:\n",
    "                taxpasta_file = f'{RESULTS_DIR}/profiling/{profiler}/pass{passn}/{sample}_mode{mode}/{sample}_mode{mode}.kreport.standardised.species'\n",
    "            elif profiler in ['kmcp']:\n",
    "                taxpasta_file = f'{RESULTS_DIR}/profiling/{profiler}/pass{passn}/{sample}_mode{mode}/{sample}_mode{mode}.profile.standardised.species'\n",
    "            elif profiler in ['ganon']:\n",
    "                taxpasta_file = f'{RESULTS_DIR}/profiling/{profiler}/pass{passn}/{sample}_mode{mode}/{sample}_mode{mode}.report.standardised.species'\n",
    "\n",
    "            if os.path.isfile(taxpasta_file):\n",
    "                counts_human, counts_others, counts_unmapped = extract_read_info(taxpasta_file, sample, passn)\n",
    "\n",
    "                with open(f'{RESULTS_DIR}/counts/profiling_counts_{sample}.txt', 'a') as file:\n",
    "                    file.write(f'{sample}\\t{passn}\\t{mode}\\t{profiler}\\t{counts_human}\\t{counts_others}\\t{counts_unmapped}\\n')\n",
    "            else:\n",
    "                print(f\"{taxpasta_file} not present!\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "for sample in CONTROLS + POOLS:\n",
    "    with open(f'{RESULTS_DIR}/counts/profiling_counts_{sample}.txt', 'w') as file:\n",
    "        file.write(f'SAMPLE\\tpass\\tmode\\tprofiler\\tmapped_human\\tmapped_others\\tunmapped\\n')\n",
    "\n",
    "    for profiler in LIST_PROFILERS:\n",
    "        for passn in [2]:\n",
    "            for mode in [3, 5, 7]:\n",
    "                if profiler in ['kaiju']:\n",
    "                    taxpasta_file = f'{RESULTS_DIR}/profiling/{profiler}/pass{passn}/{sample}_mode{mode}/{sample}_mode{mode}.results.standardised.genus'\n",
    "                elif profiler in ['krakenuniq', 'kraken2']:\n",
    "                    taxpasta_file = f'{RESULTS_DIR}/profiling/{profiler}/pass{passn}/{sample}_mode{mode}/{sample}_mode{mode}.report.standardised.genus'\n",
    "                elif profiler in ['centrifuge']:\n",
    "                    taxpasta_file = f'{RESULTS_DIR}/profiling/{profiler}/pass{passn}/{sample}_mode{mode}/{sample}_mode{mode}.kreport.standardised.genus'\n",
    "                elif profiler in ['kmcp']:\n",
    "                    taxpasta_file = f'{RESULTS_DIR}/profiling/{profiler}/pass{passn}/{sample}_mode{mode}/{sample}_mode{mode}.profile.standardised.genus'\n",
    "                elif profiler in ['ganon']:\n",
    "                    taxpasta_file = f'{RESULTS_DIR}/profiling/{profiler}/pass{passn}/{sample}_mode{mode}/{sample}_mode{mode}.report.standardised.genus'\n",
    "\n",
    "                if os.path.isfile(taxpasta_file):\n",
    "                    counts_human, counts_others, counts_unmapped = extract_read_info(taxpasta_file, sample, passn)\n",
    "\n",
    "                    with open(f'{RESULTS_DIR}/counts/profiling_counts_{sample}.txt', 'a') as file:\n",
    "                        file.write(f'{sample}\\t{passn}\\t{mode}\\t{profiler}\\t{counts_human}\\t{counts_others}\\t{counts_unmapped}\\n')\n",
    "                else:\n",
    "                    print(f\"{taxpasta_file} not present!\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "EVs",
   "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.13.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
