statistical_analysis_figure_(observed_vs_expected).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4dcc4d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The first part of this script consists of processing the data and obtaining random sampling with replacement counts\n",
    "# The second part of the script was used for plotting the counts.\n",
    "# The files generated by this script can be found in the supplementary data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3d8651d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e74f64d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_path = 'PDB_metadata.csv'\n",
    "PDB_df = pd.read_csv(file_path, sep=',',low_memory=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d22bf339",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Processing phage taxonomy column for bacteria (Prophage-DB)\n",
    "all_bacteria_df = PDB_df[PDB_df['host_domain'] == 'Bacteria'].copy()\n",
    "lineage_split = all_bacteria_df['lineage_genomad'].str.split(';')\n",
    "\n",
    "realm_list = []\n",
    "kingdom_list = []\n",
    "phylum_list = []\n",
    "class_list = []\n",
    "order_list = []\n",
    "family_list = []\n",
    "\n",
    "# Iterate through each row in the DataFrame\n",
    "for row in lineage_split:\n",
    "    # Initialize variables to store lineage information for the current row\n",
    "    realm = None\n",
    "    kingdom = None\n",
    "    phylum = None\n",
    "    class_ = None\n",
    "    order_ = None\n",
    "    family = None\n",
    "    \n",
    "    # Iterate through each element in the row\n",
    "    for element in row:\n",
    "        element = element.strip()\n",
    "        if element.endswith('viria'):\n",
    "            realm = element\n",
    "        elif element.endswith('virae'):\n",
    "            kingdom = element\n",
    "        elif element.endswith('viricota'):\n",
    "            phylum = element\n",
    "        elif element.endswith('viricetes'):\n",
    "            class_ = element\n",
    "        elif element.endswith('virales'):\n",
    "            order_ = element\n",
    "        elif element.endswith('viridae'):\n",
    "            family = element\n",
    "    \n",
    "    # Append the lineage information for the current row to the lists\n",
    "    realm_list.append(realm)\n",
    "    kingdom_list.append(kingdom)\n",
    "    phylum_list.append(phylum)\n",
    "    class_list.append(class_)\n",
    "    order_list.append(order_)\n",
    "    family_list.append(family)\n",
    "\n",
    "# Assign the lists to the DataFrame as new columns\n",
    "all_bacteria_df['realm_phage'] = realm_list\n",
    "all_bacteria_df['kingdom_phage'] = kingdom_list\n",
    "all_bacteria_df['phylum_phage'] = phylum_list\n",
    "all_bacteria_df['class_phage'] = class_list\n",
    "all_bacteria_df['order_phage'] = order_list\n",
    "all_bacteria_df['family_phage'] = family_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbac8a81",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Processing phage taxonomy column for archaea (Prophage-DB)\n",
    "all_archaea_df = PDB_df[PDB_df['host_domain'] == 'Archaea'].copy()\n",
    "lineage_split = all_archaea_df['lineage_genomad'].str.split(';')\n",
    "\n",
    "realm_list = []\n",
    "kingdom_list = []\n",
    "phylum_list = []\n",
    "class_list = []\n",
    "order_list = []\n",
    "family_list = []\n",
    "\n",
    "# Iterate through each row in the DataFrame\n",
    "for row in lineage_split:\n",
    "    # Initialize variables to store lineage information for the current row\n",
    "    realm = None\n",
    "    kingdom = None\n",
    "    phylum = None\n",
    "    class_ = None\n",
    "    order_ = None\n",
    "    family = None\n",
    "    \n",
    "    # Iterate through each element in the row\n",
    "    for element in row:\n",
    "        element = element.strip()\n",
    "        if element.endswith('viria'):\n",
    "            realm = element\n",
    "        elif element.endswith('virae'):\n",
    "            kingdom = element\n",
    "        elif element.endswith('viricota'):\n",
    "            phylum = element\n",
    "        elif element.endswith('viricetes'):\n",
    "            class_ = element\n",
    "        elif element.endswith('virales'):\n",
    "            order_ = element\n",
    "        elif element.endswith('viridae'):\n",
    "            family = element\n",
    "    \n",
    "    # Append the lineage information for the current row to the lists\n",
    "    realm_list.append(realm)\n",
    "    kingdom_list.append(kingdom)\n",
    "    phylum_list.append(phylum)\n",
    "    class_list.append(class_)\n",
    "    order_list.append(order_)\n",
    "    family_list.append(family)\n",
    "\n",
    "# Assign the lists to the DataFrame as new columns\n",
    "all_archaea_df['realm_phage'] = realm_list\n",
    "all_archaea_df['kingdom_phage'] = kingdom_list\n",
    "all_archaea_df['phylum_phage'] = phylum_list\n",
    "all_archaea_df['class_phage'] = class_list\n",
    "all_archaea_df['order_phage'] = order_list\n",
    "all_archaea_df['family_phage'] = family_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba8a9a81",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to extract information for each taxonomic group (bacterial phages)\n",
    "def extract_taxonomic_info(row, group):\n",
    "    parts = row.split(';')\n",
    "    for part in parts:\n",
    "        if part.startswith(group):\n",
    "            return part.split('__')[1]\n",
    "    return None\n",
    "\n",
    "# Select the columns of interest only if they exist in the DataFrame\n",
    "columns_of_interest = ['realm_phage','kingdom_phage','phylum_phage','class_phage','order_phage','family_phage']\n",
    "bacteria_result_df = all_bacteria_df[columns_of_interest].copy()  # Create a copy\n",
    "\n",
    "# Update specified columns with 'Unclassified' only for the selected row\n",
    "columns_to_write_unclassified = ['realm_phage','kingdom_phage','phylum_phage','class_phage','order_phage','family_phage']\n",
    "bacteria_result_df = bacteria_result_df.fillna('Unclassified')\n",
    "\n",
    "columns_to_process =  ['realm_phage','kingdom_phage','phylum_phage','class_phage','order_phage','family_phage']\n",
    "# Define the desired sample size\n",
    "desired_sample_size = 10000\n",
    "\n",
    "# DataFrame to store sampled data and word counts for each column\n",
    "column_results_df = pd.DataFrame()\n",
    "\n",
    "# Iterate through the columns\n",
    "for column in columns_to_process:\n",
    "    column_without_na = bacteria_result_df[column].dropna()\n",
    "    # Perform random sampling without replacement on the current column\n",
    "    sampled_data = column_without_na.sample(n=desired_sample_size, replace=True)\n",
    "    # Reset the index to avoid duplicate labels\n",
    "    sampled_data = sampled_data.reset_index(drop=True)\n",
    "    # Get the count of each unique word in the sampled data\n",
    "    word_counts = sampled_data.value_counts()\n",
    "    # Add results to the DataFrame\n",
    "    column_results_df[column + '_sampled_data'] = sampled_data\n",
    "    column_results_df[column + '_word_counts'] = word_counts\n",
    "\n",
    "# Exclude 'Taxonomic classification' column from the count\n",
    "columns_to_count = column_results_df.columns.str.replace('_sampled_data', '')\n",
    "\n",
    "column_counts = {}\n",
    "\n",
    "# Order of taxonomic groups\n",
    "order_of_groups =  ['realm_phage','kingdom_phage','phylum_phage','class_phage','order_phage','family_phage']\n",
    "\n",
    "for column in order_of_groups:\n",
    "    if column in columns_to_count:\n",
    "        counts = column_results_df[column + '_sampled_data'].value_counts()\n",
    "        column_counts[column] = pd.DataFrame({'Name': counts.index, 'Count': counts.values})\n",
    "\n",
    "# Create a new DataFrame by concactenating the counts\n",
    "bacteria_taxonomic_counts_df = pd.concat(column_counts.values(), keys=column_counts.keys(), axis=1)\n",
    "# Iterate over each level 'r__', 'k__', 'p__', 'c__', 'o__', 'f__', 'g__', 's__'\n",
    "for level in  ['realm_phage','kingdom_phage','phylum_phage','class_phage','order_phage','family_phage']:\n",
    "    # Get the counts for the current level\n",
    "    counts = bacteria_taxonomic_counts_df[level, 'Count']\n",
    "    # Calculate the total count for the current level\n",
    "    total_count = counts.sum()\n",
    "    # Calculate the percentage for each count within the current level\n",
    "    percentage = counts / total_count\n",
    "    # Add a new column for the percentage to the original DataFrame using .loc\n",
    "    bacteria_taxonomic_counts_df.loc[:, (level, 'Percentage')] = percentage\n",
    "    # Add a new column for the logarithm of the percentage using .loc\n",
    "    bacteria_taxonomic_counts_df.loc[:, (level, 'Logarithm')] = percentage.apply(lambda x: np.log10(x))\n",
    "\n",
    "# Reorganize the columns to have 'Name', 'Count', 'Percentage', and 'Logarithm' for each level\n",
    "new_columns = []\n",
    "for level in  ['realm_phage','kingdom_phage','phylum_phage','class_phage','order_phage','family_phage']:\n",
    "    new_columns.extend([(level, 'Name'), (level, 'Count'), (level, 'Percentage'), (level, 'Logarithm')])\n",
    "\n",
    "# Reorder the DataFrame columns\n",
    "bacteria_taxonomic_counts_df = bacteria_taxonomic_counts_df[new_columns]\n",
    "\n",
    "bacteria_taxonomic_counts_df.to_csv('counts_for_bacterial_phages_PDB.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25d002df",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to extract information for each taxonomic group (archaeal phages)\n",
    "def extract_taxonomic_info(row, group):\n",
    "    parts = row.split(';')\n",
    "    for part in parts:\n",
    "        if part.startswith(group):\n",
    "            return part.split('__')[1]\n",
    "    return None\n",
    "\n",
    "# Select the columns of interest only if they exist in the DataFrame\n",
    "columns_of_interest = ['realm_phage','kingdom_phage','phylum_phage','class_phage','order_phage','family_phage']\n",
    "archaea_result_df = all_archaea_df[columns_of_interest].copy()  # Create a copy\n",
    "\n",
    "# Update specified columns with 'Unclassified' only for the selected row\n",
    "columns_to_write_unclassified = ['realm_phage','kingdom_phage','phylum_phage','class_phage','order_phage','family_phage']\n",
    "archaea_result_df = archaea_result_df.fillna('Unclassified')\n",
    "\n",
    "columns_to_process =  ['realm_phage','kingdom_phage','phylum_phage','class_phage','order_phage','family_phage']\n",
    "# Define the desired sample size\n",
    "desired_sample_size = 10000\n",
    "\n",
    "# DataFrame to store sampled data and word counts for each column\n",
    "column_results_df = pd.DataFrame()\n",
    "\n",
    "# Iterate through the columns\n",
    "for column in columns_to_process:\n",
    "    column_without_na = archaea_result_df[column].dropna()\n",
    "    # Perform random sampling without replacement on the current column\n",
    "    sampled_data = column_without_na.sample(n=desired_sample_size, replace=True)\n",
    "    # Reset the index to avoid duplicate labels\n",
    "    sampled_data = sampled_data.reset_index(drop=True)\n",
    "    # Get the count of each unique word in the sampled data\n",
    "    word_counts = sampled_data.value_counts()\n",
    "    # Add results to the DataFrame\n",
    "    column_results_df[column + '_sampled_data'] = sampled_data\n",
    "    column_results_df[column + '_word_counts'] = word_counts\n",
    "\n",
    "# Exclude 'Taxonomic classification' column from the count\n",
    "columns_to_count = column_results_df.columns.str.replace('_sampled_data', '')\n",
    "\n",
    "column_counts = {}\n",
    "\n",
    "# Order of taxonomic groups\n",
    "order_of_groups =  ['realm_phage','kingdom_phage','phylum_phage','class_phage','order_phage','family_phage']\n",
    "\n",
    "for column in order_of_groups:\n",
    "    if column in columns_to_count:\n",
    "        counts = column_results_df[column + '_sampled_data'].value_counts()\n",
    "        column_counts[column] = pd.DataFrame({'Name': counts.index, 'Count': counts.values})\n",
    "\n",
    "# Create a new DataFrame by concactenating the counts\n",
    "archaea_taxonomic_counts_df = pd.concat(column_counts.values(), keys=column_counts.keys(), axis=1)\n",
    "# Iterate over each level 'r__', 'k__', 'p__', 'c__', 'o__', 'f__', 'g__', 's__'\n",
    "for level in  ['realm_phage','kingdom_phage','phylum_phage','class_phage','order_phage','family_phage']:\n",
    "    # Get the counts for the current level\n",
    "    counts = archaea_taxonomic_counts_df[level, 'Count']\n",
    "    # Calculate the total count for the current level\n",
    "    total_count = counts.sum()\n",
    "    # Calculate the percentage for each count within the current level\n",
    "    percentage = counts / total_count\n",
    "    # Add a new column for the percentage to the original DataFrame using .loc\n",
    "    archaea_taxonomic_counts_df.loc[:, (level, 'Percentage')] = percentage\n",
    "    # Add a new column for the logarithm of the percentage using .loc\n",
    "    archaea_taxonomic_counts_df.loc[:, (level, 'Logarithm')] = percentage.apply(lambda x: np.log10(x))\n",
    "\n",
    "# Reorganize the columns to have 'Name', 'Count', 'Percentage', and 'Logarithm' for each level\n",
    "new_columns = []\n",
    "for level in  ['realm_phage','kingdom_phage','phylum_phage','class_phage','order_phage','family_phage']:\n",
    "    new_columns.extend([(level, 'Name'), (level, 'Count'), (level, 'Percentage'), (level, 'Logarithm')])\n",
    "\n",
    "# Reorder the DataFrame columns\n",
    "archaea_taxonomic_counts_df = archaea_taxonomic_counts_df[new_columns]\n",
    "\n",
    "archaea_taxonomic_counts_df.to_csv('counts_for_archaeal_phages_PDB.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27451257",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This file contains all metadata in IMG/VR 4 that is classified as provirus\n",
    "file_path = 'imgvr_provirus.csv'\n",
    "# Read the CSV file into a DataFrame\n",
    "imgvr_df = pd.read_csv(file_path, sep=',',low_memory=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb758d57",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Processing phage taxonomy column for bacteria (IMGVR)\n",
    "def extract_taxonomic_info(row, group):\n",
    "    parts = row.split(';')\n",
    "    for part in parts:\n",
    "        if part.startswith(group):\n",
    "            return part.split('__')[1]\n",
    "    return None\n",
    "\n",
    "all_bacteria_df_copy = all_bacteria_df.copy()\n",
    "\n",
    "# Apply the function to create new columns for each taxonomic group\n",
    "for group in ['r', 'k', 'p', 'c', 'o', 'f', 'g', 's']:\n",
    "    all_bacteria_df_copy[group + '__'] = all_bacteria_df_copy['Taxonomic classification'].apply(lambda x: extract_taxonomic_info(x, group))\n",
    "\n",
    "# Select the columns of interest only if they exist in the DataFrame\n",
    "columns_of_interest = ['Taxonomic classification'] + [group + '__' for group in ['r', 'k', 'p', 'c', 'o', 'f', 'g', 's']]\n",
    "all_bacteria_df_result_df = all_bacteria_df_copy[columns_of_interest].copy()  # Create a copy\n",
    "\n",
    "# Add 'Unclassified' column based on 'Taxonomic classification'\n",
    "all_bacteria_df_result_df['Unclassified'] = all_bacteria_df_result_df['Taxonomic classification'].apply(lambda x: 'Unclassified' if ';;;;;;;' in x else '')\n",
    "\n",
    "# Move 'Unclassified' column next to 's__'\n",
    "all_bacteria_df_result_df.insert(\n",
    "    all_bacteria_df_result_df.columns.get_loc('s__') + 1,\n",
    "    'Unclassified',\n",
    "    all_bacteria_df_result_df.pop('Unclassified')\n",
    ")\n",
    "\n",
    "# Add Unclassified to DataFrame\n",
    "mask = all_bacteria_df_result_df['Taxonomic classification'] == ';;;;;;;'\n",
    "columns_to_write_unclassified = ['r__', 'k__', 'p__', 'c__', 'o__', 'f__', 'g__', 's__']\n",
    "all_bacteria_df_result_df.loc[mask, columns_to_write_unclassified] = 'Unclassified'\n",
    "all_bacteria_df_result_df = all_bacteria_df_result_df.fillna('Unclassified')\n",
    "\n",
    "# Random Sampling With Replacement\n",
    "columns_to_process = ['r__', 'k__', 'p__', 'c__', 'o__', 'f__', 'g__', 's__']\n",
    "\n",
    "# Define the desired sample size\n",
    "desired_sample_size = 10000\n",
    "\n",
    "# DataFrame to store sampled data and word counts for each column\n",
    "column_results_df = pd.DataFrame()\n",
    "\n",
    "# Iterate through the columns\n",
    "for column in columns_to_process:\n",
    "\n",
    "    column_without_na = all_bacteria_df_result_df[column].dropna()\n",
    "    # Perform random sampling without replacement on the current column\n",
    "    sampled_data = column_without_na.sample(n=desired_sample_size, replace=True)\n",
    "    \n",
    "    # Reset the index to avoid duplicate labels\n",
    "    sampled_data = sampled_data.reset_index(drop=True)\n",
    "    \n",
    "    # Get the count of each unique word in the sampled data\n",
    "    word_counts = sampled_data.value_counts()\n",
    "    \n",
    "    # Add results to the DataFrame\n",
    "    column_results_df[column + '_sampled_data'] = sampled_data\n",
    "    column_results_df[column + '_word_counts'] = word_counts\n",
    "\n",
    "# Exclude 'Taxonomic classification' column from the count\n",
    "columns_to_count = column_results_df.columns.str.replace('_sampled_data', '')\n",
    "\n",
    "column_counts = {}\n",
    "\n",
    "# Order of taxonomic groups\n",
    "order_of_groups = ['r__', 'k__', 'p__', 'c__', 'o__', 'f__', 'g__', 's__']\n",
    "\n",
    "for column in order_of_groups:\n",
    "    if column in columns_to_count:\n",
    "        counts = column_results_df[column + '_sampled_data'].value_counts()\n",
    "        column_counts[column] = pd.DataFrame({'Name': counts.index, 'Count': counts.values})\n",
    "\n",
    "# Create a new DataFrame by concatenating the counts\n",
    "all_bacteria_taxonomic_counts_df = pd.concat(column_counts.values(), keys=column_counts.keys(), axis=1)\n",
    "\n",
    "# Get proportions (percentage) and log (base10)\n",
    "for level in ['r__', 'k__', 'p__', 'c__', 'o__', 'f__', 'g__', 's__']:\n",
    "    # Get the counts for the current level\n",
    "    counts = all_bacteria_taxonomic_counts_df[level, 'Count']\n",
    "    # Calculate the total count for the current level\n",
    "    total_count = counts.sum()\n",
    "    # Calculate the percentage for each count within the current level\n",
    "    percentage = counts / total_count\n",
    "    # Add a new column for the percentage to the original DataFrame using .loc\n",
    "    all_bacteria_taxonomic_counts_df.loc[:, (level, 'Percentage')] = percentage\n",
    "    # Add a new column for the logarithm of the percentage using .loc\n",
    "    all_bacteria_taxonomic_counts_df.loc[:, (level, 'Logarithm')] = percentage.apply(lambda x: np.log10(x))\n",
    "\n",
    "# Reorganize the columns to have 'Name', 'Count', 'Percentage', and 'Logarithm' for each level\n",
    "new_columns = []\n",
    "for level in ['r__', 'k__', 'p__', 'c__', 'o__', 'f__', 'g__', 's__']:\n",
    "    new_columns.extend([(level, 'Name'), (level, 'Count'), (level, 'Percentage'), (level, 'Logarithm')])\n",
    "\n",
    "# Reorder the DataFrame columns\n",
    "all_bacteria_taxonomic_counts_df = all_bacteria_taxonomic_counts_df[new_columns]\n",
    "\n",
    "all_bacteria_taxonomic_counts_df.to_csv('counts_for_bacterial_phages_IMGVR.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fe9ca78",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Processing phage taxonomy column for archaea (IMGVR)\n",
    "def extract_taxonomic_info(row, group):\n",
    "    parts = row.split(';')\n",
    "    for part in parts:\n",
    "        if part.startswith(group):\n",
    "            return part.split('__')[1]\n",
    "    return None\n",
    "\n",
    "# Create a copy of the DataFrame to avoid SettingWithCopyWarning\n",
    "all_archaea_df_copy = all_archaea_df.copy()\n",
    "\n",
    "# Apply the function to create new columns for each taxonomic group\n",
    "for group in ['r', 'k', 'p', 'c', 'o', 'f', 'g', 's']:\n",
    "    all_archaea_df_copy[group + '__'] = all_archaea_df_copy['Taxonomic classification'].apply(lambda x: extract_taxonomic_info(x, group))\n",
    "\n",
    "# Select the columns of interest only if they exist in the DataFrame\n",
    "columns_of_interest = ['Taxonomic classification'] + [group + '__' for group in ['r', 'k', 'p', 'c', 'o', 'f', 'g', 's']]\n",
    "all_archaea_df_result_df = all_archaea_df_copy[columns_of_interest].copy()  # Create a copy\n",
    "\n",
    "# Add 'Unclassified' column based on 'Taxonomic classification'\n",
    "all_archaea_df_result_df['Unclassified'] = all_archaea_df_result_df['Taxonomic classification'].apply(lambda x: 'Unclassified' if ';;;;;;;' in x else '')\n",
    "\n",
    "# Move 'Unclassified' column next to 's__'\n",
    "all_archaea_df_result_df.insert(\n",
    "    all_archaea_df_result_df.columns.get_loc('s__') + 1,\n",
    "    'Unclassified',\n",
    "    all_archaea_df_result_df.pop('Unclassified')\n",
    ")\n",
    "\n",
    "# Add Unclassified to DataFrame\n",
    "mask = all_archaea_df_result_df['Taxonomic classification'] == ';;;;;;;'\n",
    "columns_to_write_unclassified = ['r__', 'k__', 'p__', 'c__', 'o__', 'f__', 'g__', 's__']\n",
    "all_archaea_df_result_df.loc[mask, columns_to_write_unclassified] = 'Unclassified'\n",
    "all_archaea_df_result_df = all_archaea_df_result_df.fillna('Unclassified')\n",
    "\n",
    "# Random Sampling With Replacement\n",
    "columns_to_process = ['r__', 'k__', 'p__', 'c__', 'o__', 'f__', 'g__', 's__']\n",
    "\n",
    "# Define the desired sample size\n",
    "desired_sample_size = 10000\n",
    "\n",
    "# DataFrame to store sampled data and word counts for each column\n",
    "column_results_df = pd.DataFrame()\n",
    "\n",
    "# Iterate through the columns\n",
    "for column in columns_to_process:\n",
    "\n",
    "    column_without_na = all_archaea_df_result_df[column].dropna()\n",
    "    # Perform random sampling without replacement on the current column\n",
    "    sampled_data = column_without_na.sample(n=desired_sample_size, replace=True)\n",
    "    \n",
    "    # Reset the index to avoid duplicate labels\n",
    "    sampled_data = sampled_data.reset_index(drop=True)\n",
    "    \n",
    "    # Get the count of each unique word in the sampled data\n",
    "    word_counts = sampled_data.value_counts()\n",
    "    \n",
    "    # Add results to the DataFrame\n",
    "    column_results_df[column + '_sampled_data'] = sampled_data\n",
    "    column_results_df[column + '_word_counts'] = word_counts\n",
    "\n",
    "# Exclude 'Taxonomic classification' column from the count\n",
    "columns_to_count = column_results_df.columns.str.replace('_sampled_data', '')\n",
    "\n",
    "column_counts = {}\n",
    "\n",
    "# Order of taxonomic groups\n",
    "order_of_groups = ['r__', 'k__', 'p__', 'c__', 'o__', 'f__', 'g__', 's__']\n",
    "\n",
    "for column in order_of_groups:\n",
    "    if column in columns_to_count:\n",
    "        counts = column_results_df[column + '_sampled_data'].value_counts()\n",
    "        column_counts[column] = pd.DataFrame({'Name': counts.index, 'Count': counts.values})\n",
    "\n",
    "# Create a new DataFrame by concatenating the counts\n",
    "all_archaea_taxonomic_counts_df = pd.concat(column_counts.values(), keys=column_counts.keys(), axis=1)\n",
    "\n",
    "# Get proportions (percentage) and log (base10)\n",
    "for level in ['r__', 'k__', 'p__', 'c__', 'o__', 'f__', 'g__', 's__']:\n",
    "    # Get the counts for the current level\n",
    "    counts = all_archaea_taxonomic_counts_df[level, 'Count']\n",
    "    # Calculate the total count for the current level\n",
    "    total_count = counts.sum()\n",
    "    # Calculate the percentage for each count within the current level\n",
    "    percentage = counts / total_count\n",
    "    # Add a new column for the percentage to the original DataFrame using .loc\n",
    "    all_archaea_taxonomic_counts_df.loc[:, (level, 'Percentage')] = percentage\n",
    "    # Add a new column for the logarithm of the percentage using .loc\n",
    "    all_archaea_taxonomic_counts_df.loc[:, (level, 'Logarithm')] = percentage.apply(lambda x: np.log10(x))\n",
    "\n",
    "# Reorganize the columns to have 'Name', 'Count', 'Percentage', and 'Logarithm' for each level\n",
    "new_columns = []\n",
    "for level in ['r__', 'k__', 'p__', 'c__', 'o__', 'f__', 'g__', 's__']:\n",
    "    new_columns.extend([(level, 'Name'), (level, 'Count'), (level, 'Percentage'), (level, 'Logarithm')])\n",
    "\n",
    "# Reorder the DataFrame columns\n",
    "all_archaea_taxonomic_counts_df = all_archaea_taxonomic_counts_df[new_columns]\n",
    "\n",
    "# Display the updated DataFrame\n",
    "all_archaea_taxonomic_counts_df\n",
    "all_archaea_taxonomic_counts_df.to_csv('counts_for_archaeal_phages_IMGVR.csv', index=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ff92bbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting bacterial phages\n",
    "def read_and_convert_data(file1, file2, pdb_column, imgvr_column):\n",
    "    # Read data from files\n",
    "    pdb_data = pd.read_csv(file1)\n",
    "    imgvr_data = pd.read_csv(file2)\n",
    "    \n",
    "    # Extract required columns\n",
    "    pdb_data = pdb_data[[pdb_column + '_PDB', pdb_column + '_percentage_PDB']]\n",
    "    imgvr_data = imgvr_data[[imgvr_column + '_IMGVR', imgvr_column + '_percentage_IMGVR']]\n",
    "    \n",
    "    # Remove rows with missing values (NaN)\n",
    "    pdb_data.dropna(inplace=True)\n",
    "    imgvr_data.dropna(inplace=True)\n",
    "    \n",
    "    return pdb_data, imgvr_data\n",
    "\n",
    "def plot_combined_counts(ax, data, border_width=1.5):\n",
    "    y_pos = 0\n",
    "    y_ticks = []\n",
    "    y_labels = []\n",
    "    label_colors = []\n",
    "\n",
    "    # Colors for observed and expected counts\n",
    "    observed_color = '#8ac926'\n",
    "    expected_color = '#008bf8'\n",
    "\n",
    "    # Colors for y-axis labels based on columns\n",
    "    column_colors = {\n",
    "        'Realm': '#ff6b6b',\n",
    "        'Kingdom': '#6A5ACD',\n",
    "        'Phylum': '#FFA500',\n",
    "        'Class': '#A0B4FF',\n",
    "        'Order': '#6bd425',\n",
    "        'Family': '#00BFFF'\n",
    "    }\n",
    "\n",
    "    # Plot each column's data\n",
    "    for column_name, pdb_data, imgvr_data in data:\n",
    "        # Merge dataframes on specified column with an outer join to ensure all names are included\n",
    "        merged_data = pd.merge(pdb_data, imgvr_data, how='outer', left_on=column_name + '_PDB', right_on=column_name + '_IMGVR')\n",
    "        \n",
    "        # Fill missing percentages with 0\n",
    "        merged_data[column_name + '_percentage_PDB'].fillna(0, inplace=True)\n",
    "        merged_data[column_name + '_percentage_IMGVR'].fillna(0, inplace=True)\n",
    "        \n",
    "        # Filter out rows where category values are NaN\n",
    "        merged_data = merged_data[~merged_data[column_name + '_PDB'].isna()]\n",
    "        \n",
    "        # Extract observed and expected counts from the merged dataframe\n",
    "        observed_col = column_name + '_percentage_PDB'\n",
    "        expected_col = column_name + '_percentage_IMGVR'\n",
    "        observed_counts = merged_data[observed_col].values\n",
    "        expected_counts = merged_data[expected_col].values\n",
    "        \n",
    "        # Calculate the standard deviation of each group\n",
    "        observed_std = np.std(observed_counts, ddof=1)\n",
    "        expected_std = np.std(expected_counts, ddof=1)\n",
    "\n",
    "        # Calculate the sample sizes of each group\n",
    "        observed_n = len(observed_counts)\n",
    "        expected_n = len(expected_counts)\n",
    "\n",
    "        # Calculate the standard error for each group\n",
    "        observed_se = observed_std / np.sqrt(observed_n)\n",
    "        expected_se = expected_std / np.sqrt(expected_n)\n",
    "\n",
    "        # Extract column names\n",
    "        column_names = merged_data[column_name + '_PDB'].values\n",
    "\n",
    "        # Create a horizontal bar graph with error bars and outlines\n",
    "        y = np.arange(y_pos, y_pos + len(observed_counts) * 2, 2)\n",
    "        height = 0.8  # Adjust the height to make the bars slightly smaller than the gap\n",
    "        ax.barh(y - height/2, observed_counts, height, edgecolor='black', color=observed_color, linewidth=border_width, label='Observed')\n",
    "        ax.barh(y + height/2, expected_counts, height, edgecolor='black', color=expected_color, linewidth=border_width, label='Expected')\n",
    "        ax.errorbar(observed_counts, y - height/2, xerr=observed_se, fmt='none', capsize=3, ecolor='black')\n",
    "        ax.errorbar(expected_counts, y + height/2, xerr=expected_se, fmt='none', capsize=3, ecolor='black')\n",
    "\n",
    "        # Append the current y positions and labels for setting y-ticks later\n",
    "        y_ticks.extend(y)\n",
    "        y_labels.extend(column_names)\n",
    "        label_colors.extend([column_colors[column_name]] * len(column_names))\n",
    "        y_pos += len(observed_counts) * 2  # Increase y_pos by double the number of counts to create separation\n",
    "\n",
    "    # Set y-ticks without labels\n",
    "    ax.set_yticks(y_ticks)\n",
    "    ax.set_yticklabels([])\n",
    "\n",
    "    # Add colored labels manually\n",
    "    for y, label, color in zip(y_ticks, y_labels, label_colors):\n",
    "        ax.text(-0.02, y, label, ha='right', va='center', fontsize=20, fontweight='bold', color=color, transform=ax.get_yaxis_transform())\n",
    "\n",
    "    # Set y-axis limits to fit the bars without extra space\n",
    "    ax.set_ylim(min(y_ticks) - 1, max(y_ticks) + 1)\n",
    "\n",
    "    # Customize x-axis ticks to remove the -0.2 value\n",
    "    x_ticks = ax.get_xticks()\n",
    "    x_ticks = x_ticks[x_ticks >= 0]\n",
    "    x_ticks = x_ticks[x_ticks <= 1.0]\n",
    "    ax.set_xticks(x_ticks)\n",
    "\n",
    "    # Unique legend handles and labels\n",
    "    handles, labels = ax.get_legend_handles_labels()\n",
    "    by_label = dict(zip(labels, handles))\n",
    "    ax.legend(by_label.values(), by_label.keys(), prop={'size': 20, 'weight': 'bold'}, loc='upper right')\n",
    "    \n",
    "    ax.tick_params(axis='x', labelsize=20)\n",
    "    for label in ax.get_xticklabels():\n",
    "        label.set_fontweight('bold')\n",
    "\n",
    "file1 = \"counts_for_bacterial_phages_PDB.csv\"\n",
    "file2 = \"counts_for_bacterial_phages_IMGVR.csv\"\n",
    "\n",
    "# Assuming the column names in your dataframes are 'Kingdom_PDB', 'Kingdom_IMGVR', 'Realm_PDB', 'Realm_IMGVR', etc.\n",
    "columns = ['Realm', 'Kingdom', 'Phylum', 'Class', 'Order', 'Family']\n",
    "\n",
    "data = []\n",
    "for column in columns:\n",
    "    pdb_data, imgvr_data = read_and_convert_data(file1, file2, pdb_column=column, imgvr_column=column)\n",
    "    data.append((column, pdb_data, imgvr_data))\n",
    "\n",
    "# Create a single figure with one set of axes\n",
    "fig, ax = plt.subplots(figsize=(10, 25))\n",
    "\n",
    "# Plot combined data\n",
    "plot_combined_counts(ax, data, border_width=0.01)\n",
    "\n",
    "# Add a title to the figure, centered above the bars\n",
    "ax.text(0.5, 1.02, 'Bacterial phages', fontsize=30, fontweight='bold', ha='center', transform=ax.transAxes)\n",
    "\n",
    "# Adjust the layout \n",
    "plt.tight_layout(rect=[0, 0, 1, 0.96])\n",
    "plt.subplots_adjust(top=0.95, bottom=0.05)\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cd0a336",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting archaeal phages\n",
    "\n",
    "def read_and_convert_data(file1, file2, pdb_column, imgvr_column):\n",
    "    # Read data from files\n",
    "    pdb_data = pd.read_csv(file1)\n",
    "    imgvr_data = pd.read_csv(file2)\n",
    "    \n",
    "    # Extract required columns\n",
    "    pdb_data = pdb_data[[pdb_column + '_PDB', pdb_column + '_percentage_PDB']]\n",
    "    imgvr_data = imgvr_data[[imgvr_column + '_IMGVR', imgvr_column + '_percentage_IMGVR']]\n",
    "    \n",
    "    # Remove rows with missing values (NaN)\n",
    "    pdb_data.dropna(inplace=True)\n",
    "    imgvr_data.dropna(inplace=True)\n",
    "    \n",
    "    return pdb_data, imgvr_data\n",
    "\n",
    "def plot_combined_counts(ax, data, border_width=1.5):\n",
    "    y_pos = 0\n",
    "    y_ticks = []\n",
    "    y_labels = []\n",
    "    label_colors = []\n",
    "\n",
    "    # Colors for observed and expected counts\n",
    "    observed_color = '#ff6b6b'\n",
    "    expected_color = '#008bf8'\n",
    "\n",
    "    # Colors for y-axis labels based on columns\n",
    "    column_colors = {\n",
    "        'Realm': '#ff6b6b',\n",
    "        'Kingdom': '#6A5ACD',\n",
    "        'Phylum': '#FFA500',\n",
    "        'Class': '#A0B4FF',\n",
    "        'Order': '#6bd425',\n",
    "        'Family': '#00BFFF'\n",
    "    }\n",
    "\n",
    "    # Plot each column's data\n",
    "    for column_name, pdb_data, imgvr_data in data:\n",
    "        # Merge dataframes on specified column with an outer join to ensure all names are included\n",
    "        merged_data = pd.merge(pdb_data, imgvr_data, how='outer', left_on=column_name + '_PDB', right_on=column_name + '_IMGVR')\n",
    "        \n",
    "        # Fill missing percentages with 0\n",
    "        merged_data[column_name + '_percentage_PDB'].fillna(0, inplace=True)\n",
    "        merged_data[column_name + '_percentage_IMGVR'].fillna(0, inplace=True)\n",
    "        \n",
    "        # Filter out rows where category values are NaN\n",
    "        merged_data = merged_data[~merged_data[column_name + '_PDB'].isna()]\n",
    "        \n",
    "        # Extract observed and expected counts from the merged dataframe\n",
    "        observed_col = column_name + '_percentage_PDB'\n",
    "        expected_col = column_name + '_percentage_IMGVR'\n",
    "        observed_counts = merged_data[observed_col].values\n",
    "        expected_counts = merged_data[expected_col].values\n",
    "        \n",
    "        # Calculate the standard deviation of each group\n",
    "        observed_std = np.std(observed_counts, ddof=1)\n",
    "        expected_std = np.std(expected_counts, ddof=1)\n",
    "\n",
    "        # Calculate the sample sizes of each group\n",
    "        observed_n = len(observed_counts)\n",
    "        expected_n = len(expected_counts)\n",
    "\n",
    "        # Calculate the standard error for each group\n",
    "        observed_se = observed_std / np.sqrt(observed_n)\n",
    "        expected_se = expected_std / np.sqrt(expected_n)\n",
    "\n",
    "        # Extract column names\n",
    "        column_names = merged_data[column_name + '_PDB'].values\n",
    "\n",
    "        # Create a horizontal bar graph with error bars and outlines\n",
    "        y = np.arange(y_pos, y_pos + len(observed_counts) * 2, 2)\n",
    "        height = 0.8  # Adjust the height to make the bars slightly smaller than the gap\n",
    "        ax.barh(y - height/2, observed_counts, height, edgecolor='black', color=observed_color, linewidth=border_width, label='Observed')\n",
    "        ax.barh(y + height/2, expected_counts, height, edgecolor='black', color=expected_color, linewidth=border_width, label='Expected')\n",
    "        ax.errorbar(observed_counts, y - height/2, xerr=observed_se, fmt='none', capsize=3, ecolor='black')\n",
    "        ax.errorbar(expected_counts, y + height/2, xerr=expected_se, fmt='none', capsize=3, ecolor='black')\n",
    "\n",
    "        # Append the current y positions and labels for setting y-ticks later\n",
    "        y_ticks.extend(y)\n",
    "        y_labels.extend(column_names)\n",
    "        label_colors.extend([column_colors[column_name]] * len(column_names))\n",
    "        y_pos += len(observed_counts) * 2  # Increase y_pos by double the number of counts to create separation\n",
    "\n",
    "    # Set y-ticks without labels\n",
    "    ax.set_yticks(y_ticks)\n",
    "    ax.set_yticklabels([])\n",
    "\n",
    "    # Add colored labels manually\n",
    "    for y, label, color in zip(y_ticks, y_labels, label_colors):\n",
    "        ax.text(-0.02, y, label, ha='right', va='center', fontsize=20, fontweight='bold', color=color, transform=ax.get_yaxis_transform())\n",
    "\n",
    "    # Set y-axis limits to fit the bars without extra space\n",
    "    ax.set_ylim(min(y_ticks) - 1, max(y_ticks) + 1)\n",
    "\n",
    "    # Customize x-axis ticks to remove the -0.2 value\n",
    "    x_ticks = ax.get_xticks()\n",
    "    x_ticks = x_ticks[x_ticks >= 0]\n",
    "    x_ticks = x_ticks[x_ticks <= 1.0]\n",
    "    ax.set_xticks(x_ticks)\n",
    "\n",
    "    # Get unique legend handles and labels\n",
    "    handles, labels = ax.get_legend_handles_labels()\n",
    "    by_label = dict(zip(labels, handles))\n",
    "    ax.legend(by_label.values(), by_label.keys(), prop={'size': 20, 'weight': 'bold'}, loc='upper right')\n",
    "    \n",
    "    ax.tick_params(axis='x', labelsize=20)\n",
    "    for label in ax.get_xticklabels():\n",
    "        label.set_fontweight('bold')\n",
    "\n",
    "file1 = \"counts_for_archaeal_phages_PDB.csv\"\n",
    "file2 = \"counts_for_archaeal_phages_IMGVR.csv\"\n",
    "\n",
    "# Assuming the column names in your dataframes are 'Kingdom_PDB', 'Kingdom_IMGVR', 'Realm_PDB', 'Realm_IMGVR', etc.\n",
    "columns = ['Realm', 'Kingdom', 'Phylum', 'Class', 'Order', 'Family']\n",
    "\n",
    "data = []\n",
    "for column in columns:\n",
    "    pdb_data, imgvr_data = read_and_convert_data(file1, file2, pdb_column=column, imgvr_column=column)\n",
    "    data.append((column, pdb_data, imgvr_data))\n",
    "\n",
    "# Create a single figure with one set of axes\n",
    "fig, ax = plt.subplots(figsize=(10, 25))\n",
    "\n",
    "# Plot combined data\n",
    "plot_combined_counts(ax, data, border_width=0.01)\n",
    "\n",
    "# Add a title to the figure, centered above the bars\n",
    "ax.text(0.5, 1.02, 'Archaeal phages', fontsize=30, fontweight='bold', ha='center', transform=ax.transAxes)\n",
    "\n",
    "# Adjust the layout \n",
    "plt.tight_layout(rect=[0, 0, 1, 0.96])\n",
    "plt.subplots_adjust(top=0.95, bottom=0.05)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}