NB_7_GDC_API_BERT.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "private_outputs": true,
      "provenance": [],
      "machine_shape": "hm",
      "gpuType": "V28",
      "mount_file_id": "1ZHtGgBehWmMmvUl2yr0QhM4B5J_Uh8Gb",
      "authorship_tag": "ABX9TyNcGSu70XKDB/i5GTfJoK0k",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "TPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/christophergaughan/Machine-Learning/blob/master/NB_7_GDC_API_BERT.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Define Functions to Query and Download Data from the GDC API"
      ],
      "metadata": {
        "id": "gknwbXtos2WM"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Downloading and Preparing Data\n",
        "\n",
        "In this section, we are querying the GDC API to retrieve gene expression quantification data for the TCGA-BRCA (breast cancer) project. This involves downloading the necessary files and preparing them for subsequent analysis. Several iterations of this notebook were required to ensure the successful retrieval and processing of the data. Due to this iterative approach, some cells may contain redundant code, which reflects the process of refining and troubleshooting our workflow.\n",
        "\n",
        "### Step-by-Step Explanation:\n",
        "\n",
        "1. **Directory Setup**:\n",
        "   - We start by defining the directory to save the downloaded files. If the directory already exists, we clear it to avoid any conflicts with previous data. Then, we recreate the directory to ensure it is ready for new data.\n",
        "\n",
        "    ```python\n",
        "    directory_path = '/content/drive/MyDrive/Colab Notebooks/gdc2_files'\n",
        "    \n",
        "    # Clear the directory if it exists\n",
        "    if os.path.exists(directory_path):\n",
        "        shutil.rmtree(directory_path)\n",
        "    \n",
        "    # Recreate the directory\n",
        "    os.makedirs(directory_path)\n",
        "    ```\n",
        "\n",
        "2. **Querying the GDC API**:\n",
        "   - We define a function `query_gdc_api` to query the GDC API. The function specifies parameters to filter the data by project ID (TCGA-BRCA), data category (Transcriptome Profiling), and data type (Gene Expression Quantification). The query size is set to 500 to retrieve a substantial number of files.\n",
        "\n",
        "    ```python\n",
        "    def query_gdc_api():\n",
        "        url = \"https://api.gdc.cancer.gov/files\"\n",
        "\n",
        "        params = {\n",
        "            \"filters\": json.dumps({\n",
        "                \"op\": \"and\",\n",
        "                \"content\": [\n",
        "                    {\n",
        "                        \"op\": \"in\",\n",
        "                        \"content\": {\n",
        "                            \"field\": \"cases.project.project_id\",\n",
        "                            \"value\": [\"TCGA-BRCA\"]\n",
        "                        }\n",
        "                    },\n",
        "                    {\n",
        "                        \"op\": \"in\",\n",
        "                        \"content\": {\n",
        "                            \"field\": \"data_category\",\n",
        "                            \"value\": [\"Transcriptome Profiling\"]\n",
        "                        }\n",
        "                    },\n",
        "                    {\n",
        "                        \"op\": \"in\",\n",
        "                        \"content\": {\n",
        "                            \"field\": \"data_type\",\n",
        "                            \"value\": [\"Gene Expression Quantification\"]\n",
        "                        }\n",
        "                    }\n",
        "                ]\n",
        "            }),\n",
        "            \"fields\": \"file_id,file_name,cases.submitter_id\",\n",
        "            \"format\": \"json\",\n",
        "            \"size\": \"500\"\n",
        "        }\n",
        "\n",
        "        response = requests.get(url, params=params)\n",
        "        print(f\"Query URL: {response.url}\")\n",
        "        if response.status_code == 200:\n",
        "            data = response.json()\n",
        "            return data['data']['hits']\n",
        "        else:\n",
        "            print(f\"Error querying GDC API: {response.status_code}\")\n",
        "            return []\n",
        "    ```\n",
        "\n",
        "3. **Retrieving Files**:\n",
        "   - We call the `query_gdc_api` function and print the details of the retrieved files to confirm successful retrieval.\n",
        "\n",
        "    ```python\n",
        "    files = query_gdc_api()\n",
        "    for file in files:\n",
        "        print(f\"Retrieved file: {file['file_name']} with ID: {file['file_id']}\")\n",
        "\n",
        "    if not files:\n",
        "        print(\"No files retrieved. Please check the query parameters.\")\n",
        "    else:\n",
        "        print(f\"Total files retrieved: {len(files)}\")\n",
        "    ```\n",
        "\n",
        "4. **Downloading Files**:\n",
        "   - The `download_files` function is defined to download the files retrieved from the GDC API. Each file is saved in the previously defined directory.\n",
        "\n",
        "    ```python\n",
        "    def download_files(files):\n",
        "        for file in files:\n",
        "            file_id = file['file_id']\n",
        "            file_name = file['file_name']\n",
        "            file_path = os.path.join(directory_path, file_name)\n",
        "\n",
        "            if not os.path.exists(file_path):\n",
        "                download_url = f\"https://api.gdc.cancer.gov/data/{file_id}\"\n",
        "                response = requests.get(download_url, stream=True)\n",
        "                if response.status_code == 200:\n",
        "                    with open(file_path, 'wb') as f:\n",
        "                        for chunk in response.iter_content(chunk_size=1024):\n",
        "                            if chunk:\n",
        "                                f.write(chunk)\n",
        "                    print(f\"Downloaded {file_name}\")\n",
        "                else:\n",
        "                    print(f\"Failed to download {file_name}\")\n",
        "    ```\n",
        "\n",
        "5. **Executing the Download**:\n",
        "   - Finally, we call the `download_files` function to download all the retrieved files.\n",
        "\n",
        "    ```python\n",
        "    download_files(files)\n",
        "    ```\n",
        "\n",
        "### Next Steps\n",
        "\n",
        "After successfully downloading and preparing the data, we will proceed with Exploratory Data Analysis (EDA) to understand the structure and quality of our dataset. Following the EDA, we will utilize BERT-based models to enhance our understanding through literature mining and functional annotation. While the data quality may vary, this exercise serves as a demonstration of the process and showcases the potential insights that can be derived from integrating bioinformatics with natural language processing techniques.\n"
      ],
      "metadata": {
        "id": "lY3YxtE8QkVL"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Importing Necessary Libraries\n",
        "We start by importing essential libraries: `requests` for making API calls, `json` for handling JSON data, `os` for directory operations, and `shutil` for file operations.\n",
        "This is a standard step in Python scripting where we gather all necessary libraries before starting our tasks.\n"
      ],
      "metadata": {
        "id": "VGyRqdKTaHI3"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import requests\n",
        "import json\n",
        "import os\n",
        "import shutil\n",
        "\n",
        "# Define the directory to save downloaded files\n",
        "directory_path = '/content/drive/MyDrive/Colab Notebooks/gdc2_files'\n",
        "\n",
        "# Clear the directory if it exists\n",
        "if os.path.exists(directory_path):\n",
        "    shutil.rmtree(directory_path)\n",
        "\n",
        "# Recreate the directory\n",
        "os.makedirs(directory_path)\n",
        "\n",
        "# Function to query the GDC API\n",
        "def query_gdc_api():\n",
        "    url = \"https://api.gdc.cancer.gov/files\"\n",
        "\n",
        "    # Define the parameters for the query\n",
        "    params = {\n",
        "        \"filters\": json.dumps({\n",
        "            \"op\": \"and\",\n",
        "            \"content\": [\n",
        "                {\n",
        "                    \"op\": \"in\",\n",
        "                    \"content\": {\n",
        "                        \"field\": \"cases.project.project_id\",\n",
        "                        \"value\": [\"TCGA-BRCA\"]\n",
        "                    }\n",
        "                },\n",
        "                {\n",
        "                    \"op\": \"in\",\n",
        "                    \"content\": {\n",
        "                        \"field\": \"data_category\",\n",
        "                        \"value\": [\"Transcriptome Profiling\"]\n",
        "                    }\n",
        "                },\n",
        "                {\n",
        "                    \"op\": \"in\",\n",
        "                    \"content\": {\n",
        "                        \"field\": \"data_type\",\n",
        "                        \"value\": [\"Gene Expression Quantification\"]\n",
        "                    }\n",
        "                }\n",
        "            ]\n",
        "        }),\n",
        "        \"fields\": \"file_id,file_name,cases.submitter_id\",\n",
        "        \"format\": \"json\",\n",
        "        \"size\": \"500\"  # Increase the size to 500\n",
        "    }\n",
        "\n",
        "    response = requests.get(url, params=params)\n",
        "    print(f\"Query URL: {response.url}\")  # Print the query URL for debugging\n",
        "    if response.status_code == 200:\n",
        "        data = response.json()\n",
        "        return data['data']['hits']\n",
        "    else:\n",
        "        print(f\"Error querying GDC API: {response.status_code}\")\n",
        "        return []\n",
        "\n",
        "# Query the GDC API and print the retrieved files\n",
        "files = query_gdc_api()\n",
        "for file in files:\n",
        "    print(f\"Retrieved file: {file['file_name']} with ID: {file['file_id']}\")\n",
        "\n",
        "# Ensure files are retrieved before attempting to download\n",
        "if not files:\n",
        "    print(\"No files retrieved. Please check the query parameters.\")\n",
        "else:\n",
        "    print(f\"Total files retrieved: {len(files)}\")\n",
        "\n",
        "# Function to download files from the GDC API\n",
        "def download_files(files):\n",
        "    for file in files:\n",
        "        file_id = file['file_id']\n",
        "        file_name = file['file_name']\n",
        "        file_path = os.path.join(directory_path, file_name)\n",
        "\n",
        "        if not os.path.exists(file_path):\n",
        "            download_url = f\"https://api.gdc.cancer.gov/data/{file_id}\"\n",
        "            response = requests.get(download_url, stream=True)\n",
        "            if response.status_code == 200:\n",
        "                with open(file_path, 'wb') as f:\n",
        "                    for chunk in response.iter_content(chunk_size=1024):\n",
        "                        if chunk:\n",
        "                            f.write(chunk)\n",
        "                print(f\"Downloaded {file_name}\")\n",
        "            else:\n",
        "                print(f\"Failed to download {file_name}\")\n",
        "\n",
        "# Download the files\n",
        "download_files(files)\n"
      ],
      "metadata": {
        "id": "MowKkDI2qFVC"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**preprocess the RNA-Seq files:**"
      ],
      "metadata": {
        "id": "D6Xi_XhWti50"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "import json\n",
        "import os\n",
        "\n",
        "# Function to preprocess RNA-Seq files\n",
        "def preprocess_rna_seq_file(file_path):\n",
        "    try:\n",
        "        # Read the file with proper column names and skip the initial rows if they don't conform to the structure\n",
        "        df = pd.read_csv(file_path, sep='\\t', comment='#')\n",
        "\n",
        "        # Display the first few rows to understand the structure\n",
        "        print(f\"First few rows of the file {file_path}:\")\n",
        "        print(df.head())\n",
        "\n",
        "        # Remove rows with NaN in gene_name or gene_type\n",
        "        df = df.dropna(subset=['gene_name', 'gene_type'])\n",
        "\n",
        "        # Filter out rows with gene_id values like 'N_unmapped', 'N_multimapping', etc.\n",
        "        df = df[~df['gene_id'].str.contains('N_unmapped|N_multimapping|N_noFeature|N_ambiguous')]\n",
        "\n",
        "        if df.empty or df.shape[1] == 0:\n",
        "            raise ValueError(\"File is empty or has no valid columns\")\n",
        "\n",
        "        return df\n",
        "    except Exception as e:\n",
        "        print(f\"Error reading {file_path}: {e}\")\n",
        "        return None\n",
        "\n",
        "# Directory path for downloaded files\n",
        "directory_path = '/content/drive/MyDrive/Colab Notebooks/gdc2_files'\n",
        "file_ids = [f for f in os.listdir(directory_path) if f.endswith('.rna_seq.augmented_star_gene_counts.tsv')]\n",
        "\n",
        "# List to hold all DataFrames\n",
        "rna_seq_dfs = []\n",
        "\n",
        "# Preprocess and collect all RNA-Seq DataFrames\n",
        "for file_id in file_ids:\n",
        "    file_path = os.path.join(directory_path, file_id)\n",
        "    rna_seq_df = preprocess_rna_seq_file(file_path)\n",
        "    if rna_seq_df is not None:\n",
        "        rna_seq_df['submitter_id'] = file_id.split('.')[0]  # Add submitter_id for merging\n",
        "        rna_seq_dfs.append(rna_seq_df)\n",
        "\n",
        "# Concatenate all RNA-Seq DataFrames\n",
        "if rna_seq_dfs:\n",
        "    all_rna_seq_df = pd.concat(rna_seq_dfs, ignore_index=True)\n",
        "    # Display the combined DataFrame\n",
        "    print(all_rna_seq_df.head())\n",
        "    print(all_rna_seq_df.info())\n",
        "else:\n",
        "    print(\"No valid RNA-Seq data found.\")\n"
      ],
      "metadata": {
        "id": "Va2tY5gO1my-"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Preprocessing RNA-Seq Data\n",
        "\n",
        "This block of code performs several critical steps to preprocess RNA-Seq files downloaded from the GDC API.\n",
        "\n",
        "1. **Import Libraries**: The code starts by importing the necessary libraries, `os` for directory and file operations and `pandas` for data manipulation.\n",
        "\n",
        "2. **Define Preprocessing Function**:\n",
        "    - The `preprocess_rna_seq_file` function reads an RNA-Seq file while handling comment lines.\n",
        "    - It removes rows with missing values in the 'gene_name' or 'gene_type' columns.\n",
        "    - It filters out rows where the 'gene_id' contains specific unwanted values (e.g., 'N_unmapped', 'N_multimapping').\n",
        "    - If the resulting DataFrame is empty or invalid, it raises an error.\n",
        "\n",
        "3. **Load and Preprocess Files**:\n",
        "    - The directory containing the downloaded RNA-Seq files is specified.\n",
        "    - The code lists all files in the directory that match the RNA-Seq file naming pattern.\n",
        "    - It iterates over each file, processes it using the `preprocess_rna_seq_file` function, and adds the `submitter_id` to each DataFrame for later merging.\n",
        "    - Valid DataFrames are collected into a list.\n",
        "\n",
        "4. **Concatenate DataFrames**:\n",
        "    - If any valid RNA-Seq DataFrames are found, they are concatenated into a single DataFrame.\n",
        "    - The combined DataFrame is displayed, showing the first few rows and the DataFrame information.\n",
        "    - If no valid data is found, it prints a message indicating this.\n",
        "\n",
        "This preprocessing ensures that the RNA-Seq data is cleaned and ready for further analysis.\n"
      ],
      "metadata": {
        "id": "_nd8BTwGdCDt"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import os\n",
        "import pandas as pd\n",
        "\n",
        "# Function to preprocess RNA-Seq files\n",
        "def preprocess_rna_seq_file(file_path):\n",
        "    try:\n",
        "        # Read the file and handle comment lines\n",
        "        df = pd.read_csv(file_path, sep='\\t', comment='#')\n",
        "\n",
        "        # Remove rows with NaN in 'gene_name' or 'gene_type'\n",
        "        df = df.dropna(subset=['gene_name', 'gene_type'])\n",
        "\n",
        "        # Filter out rows with gene_id values like 'N_unmapped', 'N_multimapping', etc.\n",
        "        df = df[~df['gene_id'].str.contains('N_unmapped|N_multimapping|N_noFeature|N_ambiguous')]\n",
        "\n",
        "        if df.empty or df.shape[1] == 0:\n",
        "            raise ValueError(\"File is empty or has no valid columns\")\n",
        "\n",
        "        return df\n",
        "    except Exception as e:\n",
        "        print(f\"Error reading {file_path}: {e}\")\n",
        "        return None\n",
        "\n",
        "# Directory path for downloaded files\n",
        "directory_path = '/content/drive/MyDrive/Colab Notebooks/gdc2_files'\n",
        "file_ids = [f for f in os.listdir(directory_path) if f.endswith('.rna_seq.augmented_star_gene_counts.tsv')]\n",
        "\n",
        "# List to hold all DataFrames\n",
        "rna_seq_dfs = []\n",
        "\n",
        "# Preprocess and collect all RNA-Seq DataFrames\n",
        "for file_id in file_ids:\n",
        "    file_path = os.path.join(directory_path, file_id)\n",
        "    rna_seq_df = preprocess_rna_seq_file(file_path)\n",
        "    if rna_seq_df is not None:\n",
        "        rna_seq_df['submitter_id'] = file_id.split('.')[0]  # Add submitter_id for merging\n",
        "        rna_seq_dfs.append(rna_seq_df)\n",
        "\n",
        "# Concatenate all RNA-Seq DataFrames\n",
        "if rna_seq_dfs:\n",
        "    all_rna_seq_df = pd.concat(rna_seq_dfs, ignore_index=True)\n",
        "    # Display the combined DataFrame\n",
        "    print(all_rna_seq_df.head())\n",
        "    print(all_rna_seq_df.info())\n",
        "else:\n",
        "    print(\"No valid RNA-Seq data found.\")\n"
      ],
      "metadata": {
        "id": "AKRXykpy2IXv"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Display summary statistics\n",
        "summary_stats = all_rna_seq_df.describe()\n",
        "print(summary_stats)\n"
      ],
      "metadata": {
        "id": "r4MZwfp72cx8"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Filter genes with TPM > 10\n",
        "filtered_df = all_rna_seq_df[all_rna_seq_df['tpm_unstranded'] > 10]\n",
        "print(filtered_df.head())\n"
      ],
      "metadata": {
        "id": "i77FjQEq2xmg"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Export the combined DataFrame to a CSV file\n",
        "output_file_path = '/content/drive/MyDrive/Colab Notebooks/processed_rna_seq_data.csv'\n",
        "all_rna_seq_df.to_csv(output_file_path, index=False)\n",
        "print(f\"Data exported to {output_file_path}\")\n"
      ],
      "metadata": {
        "id": "6CCReY2U25LQ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Preprocessing and Analyzing RNA-Seq Data\n",
        "\n",
        "This block of code preprocesses RNA-Seq files, concatenates them into a single DataFrame, generates summary statistics, filters the data, and exports the processed data to a CSV file.\n",
        "\n",
        "1. **Import Libraries**: The code starts by importing the `os` library for directory and file operations and the `pandas` library for data manipulation.\n",
        "\n",
        "2. **Define Preprocessing Function**:\n",
        "    - The `preprocess_rna_seq_file` function reads an RNA-Seq file while handling comment lines.\n",
        "    - It removes rows with missing values in the 'gene_name' or 'gene_type' columns.\n",
        "    - It filters out rows where the 'gene_id' contains specific unwanted values (e.g., 'N_unmapped', 'N_multimapping').\n",
        "    - If the resulting DataFrame is empty or invalid, it raises an error.\n",
        "\n",
        "3. **Load and Preprocess Files**:\n",
        "    - The directory containing the downloaded RNA-Seq files is specified.\n",
        "    - The code lists all files in the directory that match the RNA-Seq file naming pattern.\n",
        "    - It iterates over each file, processes it using the `preprocess_rna_seq_file` function, and adds the `submitter_id` to each DataFrame for later merging.\n",
        "    - Valid DataFrames are collected into a list.\n",
        "\n",
        "4. **Concatenate DataFrames**:\n",
        "    - If any valid RNA-Seq DataFrames are found, they are concatenated into a single DataFrame.\n",
        "    - The combined DataFrame is displayed, showing the first few rows and the DataFrame information.\n",
        "\n",
        "5. **Generate Summary Statistics**:\n",
        "    - Summary statistics of the concatenated DataFrame are generated and displayed using the `describe` method.\n",
        "\n",
        "6. **Filter Data**:\n",
        "    - The DataFrame is filtered to include only genes with TPM (Transcripts Per Million) greater than 10. The filtered DataFrame is displayed.\n",
        "\n",
        "7. **Export Data**:\n",
        "    - The combined DataFrame is exported to a CSV file for further analysis. The path to the exported file is printed.\n",
        "\n",
        "8. **Error Handling**:\n",
        "    - If no valid RNA-Seq data is found, a message is printed indicating this.\n",
        "\n",
        "This block ensures that the RNA-Seq data is cleaned, analyzed, and exported for further analysis, demonstrating essential data preprocessing and exploratory data analysis steps.\n"
      ],
      "metadata": {
        "id": "weRcrbDRdZP8"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import os\n",
        "import pandas as pd\n",
        "\n",
        "# Function to preprocess RNA-Seq files\n",
        "def preprocess_rna_seq_file(file_path):\n",
        "    try:\n",
        "        # Read the file and handle comment lines\n",
        "        df = pd.read_csv(file_path, sep='\\t', comment='#')\n",
        "\n",
        "        # Remove rows with NaN in 'gene_name' or 'gene_type'\n",
        "        df = df.dropna(subset=['gene_name', 'gene_type'])\n",
        "\n",
        "        # Filter out rows with gene_id values like 'N_unmapped', 'N_multimapping', etc.\n",
        "        df = df[~df['gene_id'].str.contains('N_unmapped|N_multimapping|N_noFeature|N_ambiguous')]\n",
        "\n",
        "        if df.empty or df.shape[1] == 0:\n",
        "            raise ValueError(\"File is empty or has no valid columns\")\n",
        "\n",
        "        return df\n",
        "    except Exception as e:\n",
        "        print(f\"Error reading {file_path}: {e}\")\n",
        "        return None\n",
        "\n",
        "# Directory path for downloaded files\n",
        "directory_path = '/content/drive/MyDrive/Colab Notebooks/gdc2_files'\n",
        "file_ids = [f for f in os.listdir(directory_path) if f.endswith('.rna_seq.augmented_star_gene_counts.tsv')]\n",
        "\n",
        "# List to hold all DataFrames\n",
        "rna_seq_dfs = []\n",
        "\n",
        "# Preprocess and collect all RNA-Seq DataFrames\n",
        "for file_id in file_ids:\n",
        "    file_path = os.path.join(directory_path, file_id)\n",
        "    rna_seq_df = preprocess_rna_seq_file(file_path)\n",
        "    if rna_seq_df is not None:\n",
        "        rna_seq_df['submitter_id'] = file_id.split('.')[0]  # Add submitter_id for merging\n",
        "        rna_seq_dfs.append(rna_seq_df)\n",
        "\n",
        "# Concatenate all RNA-Seq DataFrames\n",
        "if rna_seq_dfs:\n",
        "    all_rna_seq_df = pd.concat(rna_seq_dfs, ignore_index=True)\n",
        "    # Display the combined DataFrame\n",
        "    print(all_rna_seq_df.head())\n",
        "    print(all_rna_seq_df.info())\n",
        "\n",
        "    # Generate summary statistics\n",
        "    summary_stats = all_rna_seq_df.describe()\n",
        "    print(summary_stats)\n",
        "\n",
        "    # Filter genes with TPM > 10\n",
        "    filtered_df = all_rna_seq_df[all_rna_seq_df['tpm_unstranded'] > 10]\n",
        "    print(filtered_df.head())\n",
        "\n",
        "    # Export the combined DataFrame to a CSV file\n",
        "    output_file_path = '/content/drive/MyDrive/Colab Notebooks/processed_rna_seq_data.csv'\n",
        "    all_rna_seq_df.to_csv(output_file_path, index=False)\n",
        "    print(f\"Data exported to {output_file_path}\")\n",
        "else:\n",
        "    print(\"No valid RNA-Seq data found.\")\n"
      ],
      "metadata": {
        "id": "sjncjmml3DGo"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "\n",
        "# Load the data\n",
        "data = pd.read_csv('/content/drive/MyDrive/Colab Notebooks/processed_rna_seq_data.csv')\n",
        "\n",
        "# Check for missing values\n",
        "print(data.isnull().sum())\n",
        "\n",
        "# Summary statistics\n",
        "print(data.describe())\n"
      ],
      "metadata": {
        "id": "mKD1sQZG4Jw4"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Visualizing Data Distributions\n",
        "\n",
        "This block of code uses Seaborn and Matplotlib to plot the distributions of selected columns in the RNA-Seq data. The visualizations include histograms and boxplots for easier understanding of the data distributions and potential outliers.\n",
        "\n",
        "1. **Import Libraries**:\n",
        "    - The `seaborn` library is imported for creating statistical visualizations.\n",
        "    - The `matplotlib.pyplot` library is imported for creating static, interactive, and animated visualizations.\n",
        "\n",
        "2. **Define Plotting Function**:\n",
        "    - The `plot_distributions` function plots histograms and boxplots for the specified columns in the data.\n",
        "    - It samples the data to 100,000 rows for quicker processing and visualization.\n",
        "\n",
        "3. **Plot Distributions**:\n",
        "    - The function iterates over the specified columns and generates a figure for each column.\n",
        "        - **Histogram**: The left subplot shows the histogram with a KDE (Kernel Density Estimate) overlay to visualize the distribution of values.\n",
        "        - **Boxplot**: The right subplot shows the boxplot to highlight the median, quartiles, and potential outliers in the data.\n",
        "    - The subplots are laid out neatly using `tight_layout` to prevent overlap.\n",
        "\n",
        "4. **Columns to Visualize**:\n",
        "    - A list of columns (`'tpm_unstranded'` and `'fpkm_unstranded'`) is specified for visualization.\n",
        "\n",
        "5. **Generate Plots**:\n",
        "    - The `plot_distributions` function is called with the data and the specified columns, generating and displaying the histograms and boxplots for each column.\n",
        "\n",
        "This block provides a visual representation of the distributions of selected columns, aiding in the identification of data patterns and outliers.\n"
      ],
      "metadata": {
        "id": "ZKVdoz6Rux9o"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import seaborn as sns\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Function to plot distributions\n",
        "def plot_distributions(data, columns, sample_size=100000):\n",
        "    # Sample the data for quicker processing\n",
        "    sampled_data = data.sample(n=sample_size, random_state=42)\n",
        "\n",
        "    for column in columns:\n",
        "        plt.figure(figsize=(14, 6))\n",
        "\n",
        "        # Plot histogram\n",
        "        plt.subplot(1, 2, 1)\n",
        "        sns.histplot(sampled_data[column], kde=True, bins=50)  # Limiting bins for performance\n",
        "        plt.title(f'Distribution of {column}')\n",
        "        plt.xlabel(column)\n",
        "        plt.ylabel('Frequency')\n",
        "\n",
        "        # Plot boxplot\n",
        "        plt.subplot(1, 2, 2)\n",
        "        sns.boxplot(y=sampled_data[column])\n",
        "        plt.title(f'Boxplot of {column}')\n",
        "\n",
        "        plt.tight_layout()\n",
        "        plt.show()\n",
        "\n",
        "# Columns to visualize\n",
        "columns_to_plot = ['tpm_unstranded', 'fpkm_unstranded']\n",
        "\n",
        "# Plot the distributions with sampling\n",
        "plot_distributions(data, columns_to_plot)\n"
      ],
      "metadata": {
        "id": "6jgQN3QWyPAl"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Discussion of Data Distribution Plots\n",
        "\n",
        "### Distribution of TPM (Transcripts Per Million)\n",
        "\n",
        "1. **Histogram of TPM**:\n",
        "    - The histogram shows the frequency distribution of TPM values across the sampled RNA-Seq data.\n",
        "    - The plot indicates a right-skewed distribution, with most of the gene expression values clustered at the lower end of the range.\n",
        "    - The KDE overlay provides a smoothed estimate of the distribution, highlighting the peak and spread of the data.\n",
        "\n",
        "2. **Boxplot of TPM**:\n",
        "    - The boxplot reveals the central tendency and dispersion of the TPM values.\n",
        "    - The median is represented by the line inside the box, and the box extends to the first and third quartiles.\n",
        "    - The presence of several outliers is evident, as indicated by the points beyond the whiskers. This suggests that while most genes have low TPM values, a few genes are highly expressed.\n",
        "\n",
        "### Distribution of FPKM (Fragments Per Kilobase of transcript per Million mapped reads)\n",
        "\n",
        "1. **Histogram of FPKM**:\n",
        "    - Similar to the TPM histogram, the FPKM histogram shows a right-skewed distribution.\n",
        "    - Most of the gene expression values are concentrated at the lower end, with a long tail extending towards higher values.\n",
        "    - The KDE overlay helps in visualizing the peak and overall distribution trend.\n",
        "\n",
        "2. **Boxplot of FPKM**:\n",
        "    - The boxplot for FPKM values displays the median, quartiles, and range of the data.\n",
        "    - The plot highlights a significant number of outliers, indicating that a small number of genes have much higher expression levels compared to the rest.\n",
        "    - The spread between the quartiles suggests variability in gene expression levels across the dataset.\n",
        "\n",
        "### Summary\n",
        "\n",
        "The visualizations highlight the following key points:\n",
        "- Both TPM and FPKM distributions are right-skewed, with a large number of genes having low expression levels and a few genes exhibiting high expression levels.\n",
        "- The presence of numerous outliers suggests that gene expression varies widely, which is typical in RNA-Seq data due to the biological variability and differences in gene regulation.\n",
        "- These plots are crucial for understanding the overall distribution of gene expression values and identifying genes with exceptionally high or low expression, which may be of biological interest.\n",
        "\n",
        "These insights can guide further analysis and interpretation of the RNA-Seq data, helping in identifying differentially expressed genes and understanding the underlying biological processes.\n"
      ],
      "metadata": {
        "id": "C4I_4d7xeFn2"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "The distribution of tpm_unstranded is heavily skewed with many values close to zero and some extremely high values. The boxplot also shows the presence of numerous outliers.\n",
        "\n",
        "To get more insights from the data, we can perform the following steps:\n",
        "\n",
        "1. **Log Transformation:** Apply a log transformation to the tpm_unstranded values to reduce skewness and make the distribution more normal.\n",
        "2. **Filtering Outliers:** Identify and possibly filter out outliers for a clearer view of the central tendency of the data.\n",
        "3. **Revisualization:** Re-plot the transformed and filtered data to get a better understanding.\n"
      ],
      "metadata": {
        "id": "zYuizulByttl"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Applying Log Transformation and Filtering Outliers"
      ],
      "metadata": {
        "id": "lx56dniazKPM"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import seaborn as sns\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Function to plot distributions with transformations\n",
        "def plot_transformed_distributions(data, column, sample_size=100000):\n",
        "    # Sample the data for quicker processing\n",
        "    sampled_data = data.sample(n=sample_size, random_state=42)\n",
        "\n",
        "    # Apply log transformation (adding a small constant to avoid log(0))\n",
        "    sampled_data['log_' + column] = np.log1p(sampled_data[column])\n",
        "\n",
        "    plt.figure(figsize=(14, 6))\n",
        "\n",
        "    # Plot histogram of log-transformed data\n",
        "    plt.subplot(1, 2, 1)\n",
        "    sns.histplot(sampled_data['log_' + column], kde=True, bins=50)  # Limiting bins for performance\n",
        "    plt.title(f'Log-Transformed Distribution of {column}')\n",
        "    plt.xlabel('log(' + column + ')')\n",
        "    plt.ylabel('Frequency')\n",
        "\n",
        "    # Plot boxplot of log-transformed data\n",
        "    plt.subplot(1, 2, 2)\n",
        "    sns.boxplot(y=sampled_data['log_' + column])\n",
        "    plt.title(f'Boxplot of Log-Transformed {column}')\n",
        "\n",
        "    plt.tight_layout()\n",
        "    plt.show()\n",
        "\n",
        "# Column to visualize\n",
        "column_to_plot = 'tpm_unstranded'\n",
        "\n",
        "# Plot the log-transformed distributions\n",
        "plot_transformed_distributions(data, column_to_plot)\n"
      ],
      "metadata": {
        "id": "YRE5NALAzAzU"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Discussion of Log-Transformed Data Plots\n",
        "\n",
        "### Log-Transformed Distribution of TPM (Transcripts Per Million)\n",
        "\n",
        "1. **Histogram of Log-Transformed TPM**:\n",
        "    - The histogram displays the frequency distribution of log-transformed TPM values across the sampled RNA-Seq data.\n",
        "    - By applying a log transformation, we reduce the skewness present in the original TPM distribution. The data appears more symmetric, making it easier to analyze and interpret.\n",
        "    - The KDE overlay on the histogram helps in visualizing the smoothed density estimate of the log-transformed TPM values, highlighting the central tendency and spread more clearly.\n",
        "\n",
        "2. **Boxplot of Log-Transformed TPM**:\n",
        "    - The boxplot of log-transformed TPM values shows the central tendency, variability, and presence of outliers in the data.\n",
        "    - The log transformation compresses the range of TPM values, making the median, quartiles, and spread more interpretable.\n",
        "    - Although the transformation reduces the influence of extreme outliers, some data points still lie beyond the whiskers, indicating a few genes with exceptionally high or low expression levels even after transformation.\n",
        "\n",
        "### Summary\n",
        "\n",
        "The log transformation of TPM values offers several advantages for data analysis:\n",
        "- **Normalization**: The transformation helps normalize the data, reducing the impact of extreme values and making the distribution more symmetric.\n",
        "- **Improved Visualization**: The log-transformed histogram and boxplot provide a clearer view of the central tendency and variability, facilitating better interpretation of gene expression levels.\n",
        "- **Enhanced Analysis**: Transforming the data can improve the performance of statistical analyses and machine learning algorithms, which often assume normally distributed input data.\n",
        "\n",
        "By log-transforming the TPM values, we achieve a more balanced distribution, making it easier to identify significant patterns and trends in gene expression. This transformation is particularly useful for subsequent analyses, such as differential gene expression analysis and machine learning modeling.\n"
      ],
      "metadata": {
        "id": "hPW9iE4bej2r"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "The log-transformed distribution and boxplot for tpm_unstranded provide a clearer view of the data:\n",
        "\n",
        "Log-Transformed Distribution: The distribution is still right-skewed but much less so than the original distribution. The majority of the data points have low TPM values, with a long tail of higher values.\n",
        "Boxplot: The boxplot shows that even after log transformation, there are still many outliers.\n"
      ],
      "metadata": {
        "id": "NZyig9iDzloP"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Identify the top 10 most highly expressed genes\n",
        "top_expressed_genes = data[['gene_id', 'gene_name', 'tpm_unstranded']].nlargest(10, 'tpm_unstranded')\n",
        "print(\"Top 10 most highly expressed genes:\")\n",
        "print(top_expressed_genes)\n",
        "\n",
        "# Identify the top 10 least expressed genes (excluding zeros)\n",
        "least_expressed_genes = data[['gene_id', 'gene_name', 'tpm_unstranded']][data['tpm_unstranded'] > 0].nsmallest(10, 'tpm_unstranded')\n",
        "print(\"\\nTop 10 least expressed genes (excluding zeros):\")\n",
        "print(least_expressed_genes)\n"
      ],
      "metadata": {
        "id": "GPfSRmJTzps2"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "print(data.columns)\n",
        "\n"
      ],
      "metadata": {
        "id": "hetoYZUM0HR1"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Identify the top 10 most highly expressed genes\n",
        "top_expressed_genes = data[['gene_id', 'gene_name', 'tpm_unstranded']].nlargest(10, 'tpm_unstranded')\n",
        "print(\"Top 10 most highly expressed genes:\")\n",
        "print(top_expressed_genes)\n",
        "\n",
        "# Identify the top 10 least expressed genes (excluding zeros)\n",
        "least_expressed_genes = data[['gene_id', 'gene_name', 'tpm_unstranded']][data['tpm_unstranded'] > 0].nsmallest(10, 'tpm_unstranded')\n",
        "print(\"\\nTop 10 least expressed genes (excluding zeros):\")\n",
        "print(least_expressed_genes)\n"
      ],
      "metadata": {
        "id": "DUw1kItA0jAt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "\n",
        "# Simulate condition data (e.g., two conditions 'A' and 'B')\n",
        "np.random.seed(42)  # For reproducibility\n",
        "data['condition'] = np.random.choice(['A', 'B'], size=len(data))\n",
        "\n",
        "# Calculate mean TPM for each gene per condition\n",
        "mean_tpm_per_condition = data.groupby(['gene_id', 'gene_name', 'condition'])['tpm_unstranded'].mean().unstack()\n",
        "\n",
        "# Display the mean TPM per condition for the first few genes\n",
        "print(mean_tpm_per_condition.head())\n"
      ],
      "metadata": {
        "id": "V3nnf_PD0njd"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Identify the top 10 most highly expressed genes for each condition\n",
        "top_expressed_genes_A = mean_tpm_per_condition['A'].nlargest(10)\n",
        "top_expressed_genes_B = mean_tpm_per_condition['B'].nlargest(10)\n",
        "\n",
        "print(\"Top 10 most highly expressed genes in condition A:\")\n",
        "print(top_expressed_genes_A)\n",
        "\n",
        "print(\"\\nTop 10 most highly expressed genes in condition B:\")\n",
        "print(top_expressed_genes_B)\n",
        "\n",
        "# Identify the top 10 least expressed genes (excluding zeros) for each condition\n",
        "least_expressed_genes_A = mean_tpm_per_condition['A'][mean_tpm_per_condition['A'] > 0].nsmallest(10)\n",
        "least_expressed_genes_B = mean_tpm_per_condition['B'][mean_tpm_per_condition['B'] > 0].nsmallest(10)\n",
        "\n",
        "print(\"\\nTop 10 least expressed genes in condition A (excluding zeros):\")\n",
        "print(least_expressed_genes_A)\n",
        "\n",
        "print(\"\\nTop 10 least expressed genes in condition B (excluding zeros):\")\n",
        "print(least_expressed_genes_B)\n"
      ],
      "metadata": {
        "id": "rxQPEBAK08zU"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Calculate fold change\n",
        "mean_tpm_per_condition['fold_change'] = mean_tpm_per_condition['B'] / mean_tpm_per_condition['A']\n",
        "\n",
        "# Log2 Fold Change (optional for better interpretation)\n",
        "mean_tpm_per_condition['log2_fold_change'] = np.log2(mean_tpm_per_condition['fold_change'])\n",
        "\n",
        "# Display the data with fold change\n",
        "print(mean_tpm_per_condition[['A', 'B', 'fold_change', 'log2_fold_change']].head())\n"
      ],
      "metadata": {
        "id": "j9cp-HiI1AVM"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Replace zero values with a small number\n",
        "mean_tpm_per_condition.replace(0, np.nan, inplace=True)\n",
        "mean_tpm_per_condition['fold_change'] = mean_tpm_per_condition['B'] / mean_tpm_per_condition['A']\n",
        "mean_tpm_per_condition['log2_fold_change'] = np.log2(mean_tpm_per_condition['fold_change'])\n",
        "mean_tpm_per_condition.replace(np.nan, 0, inplace=True)\n"
      ],
      "metadata": {
        "id": "8DWCFRKj1iJY"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Define fold-change threshold\n",
        "fc_threshold = 2\n",
        "\n",
        "# Identify upregulated genes (log2 fold-change > 1)\n",
        "upregulated_genes = mean_tpm_per_condition[mean_tpm_per_condition['log2_fold_change'] > 1]\n",
        "\n",
        "# Identify downregulated genes (log2 fold-change < -1)\n",
        "downregulated_genes = mean_tpm_per_condition[mean_tpm_per_condition['log2_fold_change'] < -1]\n",
        "\n",
        "print(\"Upregulated genes (log2 fold-change > 1):\")\n",
        "print(upregulated_genes)\n",
        "\n",
        "print(\"\\nDownregulated genes (log2 fold-change < -1):\")\n",
        "print(downregulated_genes)\n"
      ],
      "metadata": {
        "id": "6THkEAeY1l5x"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Define the paths\n",
        "upregulated_path = '/content/drive/MyDrive/Colab Notebooks/upregulated_genes.csv'\n",
        "downregulated_path = '/content/drive/MyDrive/Colab Notebooks/downregulated_genes.csv'\n",
        "\n",
        "# Save the results to CSV files\n",
        "upregulated_genes.to_csv(upregulated_path, header=True)\n",
        "downregulated_genes.to_csv(downregulated_path, header=True)\n"
      ],
      "metadata": {
        "id": "cvW5hmLw2xtr"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "\n",
        "# Mount Google Drive\n",
        "from google.colab import drive\n",
        "drive.mount('/content/drive')\n",
        "\n",
        "# Replace zero values with a small number to avoid log transformation issues\n",
        "mean_tpm_per_condition.replace(0, np.nan, inplace=True)\n",
        "\n",
        "# Calculate fold change and log2 fold change\n",
        "mean_tpm_per_condition['fold_change'] = mean_tpm_per_condition['B'] / mean_tpm_per_condition['A']\n",
        "mean_tpm_per_condition['log2_fold_change'] = np.log2(mean_tpm_per_condition['fold_change'])\n",
        "\n",
        "# Replace NaN values back with zero for meaningful interpretation\n",
        "mean_tpm_per_condition.replace(np.nan, 0, inplace=True)\n",
        "\n",
        "# Display the data with fold change\n",
        "print(mean_tpm_per_condition[['A', 'B', 'fold_change', 'log2_fold_change']].head())\n",
        "\n",
        "# Define fold-change threshold\n",
        "fc_threshold = 2\n",
        "\n",
        "# Identify upregulated genes (log2 fold-change > 1)\n",
        "upregulated_genes = mean_tpm_per_condition[mean_tpm_per_condition['log2_fold_change'] > 1]\n",
        "\n",
        "# Identify downregulated genes (log2 fold-change < -1)\n",
        "downregulated_genes = mean_tpm_per_condition[mean_tpm_per_condition['log2_fold_change'] < -1]\n",
        "\n",
        "print(\"Upregulated genes (log2 fold-change > 1):\")\n",
        "print(upregulated_genes)\n",
        "\n",
        "print(\"\\nDownregulated genes (log2 fold-change < -1):\")\n",
        "print(downregulated_genes)\n",
        "\n",
        "# Define the paths\n",
        "upregulated_path = '/content/drive/MyDrive/Colab Notebooks/upregulated_genes.csv'\n",
        "downregulated_path = '/content/drive/MyDrive/Colab Notebooks/downregulated_genes.csv'\n",
        "\n",
        "# Save the results to CSV files\n",
        "upregulated_genes.to_csv(upregulated_path, header=True)\n",
        "downregulated_genes.to_csv(downregulated_path, header=True)\n"
      ],
      "metadata": {
        "id": "DLhxsMTK1r9J"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Principal Component Analysis (PCA) of RNA-Seq Data\n",
        "\n",
        "### PCA Overview\n",
        "\n",
        "Principal Component Analysis (PCA) is a dimensionality reduction technique used to reduce the complexity of large datasets while preserving as much variability as possible. By projecting the data onto a lower-dimensional space, PCA helps in identifying patterns and visualizing high-dimensional data in a more interpretable manner.\n",
        "\n",
        "### Results\n",
        "\n",
        "1. **Data Preparation and Standardization**:\n",
        "    - The RNA-Seq data for conditions A and B was extracted and standardized using the `StandardScaler`. Standardization is crucial in PCA to ensure that each feature contributes equally to the analysis.\n",
        "\n",
        "2. **PCA Application**:\n",
        "    - PCA was applied to the standardized data, and the first two principal components (PC1 and PC2) were extracted. These components capture the maximum variance in the data, with PC1 capturing the most variance and PC2 capturing the second most variance.\n",
        "\n",
        "3. **Explained Variance**:\n",
        "    - The explained variance ratios indicate how much variance is captured by each principal component. In this analysis, PC1 and PC2 capture a combined variance, providing a good summary of the overall data variability.\n",
        "\n",
        "4. **PCA Plot**:\n",
        "    - The scatter plot shows the distribution of genes in the space defined by the first two principal components.\n",
        "    - Each point represents a gene, positioned according to its expression patterns in conditions A and B.\n",
        "    - The plot helps in identifying clusters of genes with similar expression profiles, as well as potential outliers.\n",
        "\n",
        "### Interpretation\n",
        "\n",
        "- **Variance Explained**:\n",
        "    - PC1 and PC2 together capture a significant portion of the total variance in the RNA-Seq data. This indicates that these two components are sufficient to provide a simplified yet informative view of the data.\n",
        "\n",
        "- **Clusters and Patterns**:\n",
        "    - The scatter plot reveals clusters of genes that have similar expression patterns in the two conditions. These clusters may represent groups of genes that are co-regulated or involved in related biological pathways.\n",
        "    - Outliers can also be identified in the plot. These are genes with expression patterns that deviate significantly from the majority, potentially indicating unique regulatory mechanisms or experimental anomalies.\n",
        "\n",
        "- **Biological Insights**:\n",
        "    - Genes that are closer together in the PCA plot have more similar expression profiles. This information can be used to infer functional relationships between genes and to identify candidate genes for further experimental validation.\n",
        "    - The separation of clusters may indicate differences in gene expression between cancerous and non-cancerous tissues, providing insights into the molecular mechanisms underlying cancer progression or suppression.\n",
        "\n",
        "### Conclusion\n",
        "\n",
        "PCA is a powerful tool for simplifying complex RNA-Seq data and uncovering hidden patterns. The PCA plot provides a visual representation of the relationships between genes based on their expression levels in different conditions. This analysis highlights the potential of PCA to identify key genes and pathways involved in cancer biology, guiding future research and therapeutic strategies.\n"
      ],
      "metadata": {
        "id": "4bY33QSRfPXz"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "from sklearn.decomposition import PCA\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "\n",
        "# Assuming `mean_tpm_per_condition` is the DataFrame containing your data\n",
        "# with 'A' and 'B' as conditions and 'gene_id', 'gene_name' as indices\n",
        "\n",
        "# Prepare the data\n",
        "pca_data = mean_tpm_per_condition[['A', 'B']].fillna(0)\n",
        "\n",
        "# Standardize the data\n",
        "scaler = StandardScaler()\n",
        "pca_data_scaled = scaler.fit_transform(pca_data)\n",
        "\n",
        "# Apply PCA\n",
        "pca = PCA(n_components=2)\n",
        "pca_result = pca.fit_transform(pca_data_scaled)\n",
        "\n",
        "# Create a DataFrame with the PCA results\n",
        "pca_df = pd.DataFrame(data=pca_result, columns=['PC1', 'PC2'])\n",
        "pca_df['gene_id'] = mean_tpm_per_condition.index.get_level_values('gene_id')\n",
        "pca_df['gene_name'] = mean_tpm_per_condition.index.get_level_values('gene_name')\n",
        "\n",
        "# Plot the PCA results\n",
        "plt.figure(figsize=(10, 7))\n",
        "sns.scatterplot(x='PC1', y='PC2', data=pca_df)\n",
        "plt.title('PCA of RNA-seq Data')\n",
        "plt.xlabel(f'Principal Component 1 ({pca.explained_variance_ratio_[0]*100:.2f}%)')\n",
        "plt.ylabel(f'Principal Component 2 ({pca.explained_variance_ratio_[1]*100:.2f}%)')\n",
        "plt.show()\n"
      ],
      "metadata": {
        "id": "6bKPlG1orPjd"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "##  Identify Outliers"
      ],
      "metadata": {
        "id": "PReU6N6N3ssG"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Principal Component Analysis (PCA) with Outliers Highlighted\n",
        "\n",
        "### Identifying Outliers\n",
        "\n",
        "In our PCA analysis, we have highlighted the outliers to understand their impact on the data and to gain insights into the variability of gene expression. Outliers can indicate significant deviations from the typical expression patterns, which may be of biological interest or could indicate data quality issues.\n",
        "\n",
        "### Results\n",
        "\n",
        "1. **Outlier Identification**:\n",
        "    - We defined outliers as the top 1% of genes based on their PC1 values. This threshold helps in focusing on the most extreme deviations in the data.\n",
        "    - Outliers were filtered and highlighted in the PCA plot to visualize their distribution relative to the rest of the data points.\n",
        "\n",
        "2. **PCA Plot with Outliers**:\n",
        "    - The scatter plot shows the distribution of genes in the space defined by the first two principal components (PC1 and PC2).\n",
        "    - Outliers are marked in red, distinguishing them from the majority of data points.\n",
        "\n",
        "### Interpretation\n",
        "\n",
        "- **Biological Significance**:\n",
        "    - Outliers can represent genes with exceptionally high or low expression levels in one or both conditions. These genes may be involved in unique regulatory pathways or stress responses.\n",
        "    - Identifying outliers can help pinpoint genes that warrant further investigation due to their unusual expression patterns.\n",
        "\n",
        "- **Data Quality**:\n",
        "    - The presence of numerous outliers may also indicate potential issues with data quality, such as technical artifacts, batch effects, or errors in data processing.\n",
        "    - It's important to carefully examine outliers to determine whether they reflect true biological variation or if they are artifacts that need to be corrected.\n",
        "\n",
        "- **Impact on Analysis**:\n",
        "    - Outliers can skew the results of statistical analyses and affect the interpretation of data. By identifying and understanding outliers, we can decide whether to include or exclude them in further analyses.\n",
        "    - Highlighting outliers allows us to communicate the variability and potential limitations of our dataset transparently.\n",
        "\n",
        "### Conclusion\n",
        "\n",
        "Outliers in PCA provide valuable insights into the variability and quality of our RNA-Seq data. While some outliers may indicate important biological phenomena, others might reflect technical issues that need to be addressed. Careful consideration of outliers is crucial for accurate data interpretation and for guiding subsequent analyses in cancer gene expression studies.\n"
      ],
      "metadata": {
        "id": "UPvuwxlrf0ou"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "\n",
        "# Add labels to the PCA DataFrame\n",
        "pca_df['PC1'] = pca_result[:, 0]\n",
        "pca_df['PC2'] = pca_result[:, 1]\n",
        "\n",
        "# Define a threshold to identify outliers (e.g., top 1% of PC1 values)\n",
        "threshold = np.percentile(pca_df['PC1'], 99)\n",
        "\n",
        "# Filter the outliers\n",
        "outliers = pca_df[pca_df['PC1'] > threshold]\n",
        "\n",
        "# Plot the PCA results with outliers highlighted\n",
        "plt.figure(figsize=(10, 7))\n",
        "sns.scatterplot(x='PC1', y='PC2', data=pca_df, label='Data Points')\n",
        "sns.scatterplot(x='PC1', y='PC2', data=outliers, color='red', label='Outliers')\n",
        "plt.title('PCA of RNA-seq Data with Outliers Highlighted')\n",
        "plt.xlabel(f'Principal Component 1 ({pca.explained_variance_ratio_[0]*100:.2f}%)')\n",
        "plt.ylabel(f'Principal Component 2 ({pca.explained_variance_ratio_[1]*100:.2f}%)')\n",
        "plt.legend()\n",
        "plt.show()\n"
      ],
      "metadata": {
        "id": "D1BAHxjX3a8C"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Extract outliers based on the threshold defined earlier\n",
        "threshold = np.percentile(pca_df['PC1'], 99)\n",
        "outliers = pca_df[pca_df['PC1'] > threshold]\n",
        "\n",
        "# Save outliers to a CSV file for further investigation\n",
        "outliers.to_csv('/content/drive/MyDrive/Colab Notebooks/rna_seq_outliers.csv', index=False)\n"
      ],
      "metadata": {
        "id": "CxKBs4Rm3wyf"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install statsmodels\n"
      ],
      "metadata": {
        "id": "pxu8Yo9D4R8Z"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import statsmodels.api as sm\n",
        "import numpy as np\n",
        "\n",
        "# Ensure PCA DataFrame has a named index for gene_id\n",
        "pca_df.index.name = 'gene_id'\n",
        "\n",
        "# Extract outlier gene IDs based on the threshold defined earlier\n",
        "threshold = np.percentile(pca_df['PC1'], 99)\n",
        "outliers = pca_df[pca_df['PC1'] > threshold]\n",
        "\n",
        "# Extract the gene IDs of the outliers\n",
        "outlier_gene_ids = outliers.index.unique()\n",
        "\n",
        "# Filter the expression data for the outlier gene IDs, ensuring alignment\n",
        "common_gene_ids = mean_tpm_per_condition.index.intersection(outlier_gene_ids)\n",
        "expression_data = mean_tpm_per_condition.loc[common_gene_ids]\n",
        "\n",
        "# Perform differential expression analysis using statsmodels\n",
        "results = []\n",
        "for gene in expression_data.index:\n",
        "    try:\n",
        "        # Use log-transformed data to stabilize variance\n",
        "        y = np.log2(expression_data.loc[gene, 'B'] + 1)  # Dependent variable\n",
        "        x = np.log2(expression_data.loc[gene, 'A'] + 1)  # Independent variable\n",
        "        x = sm.add_constant(x)  # Add intercept\n",
        "        model = sm.OLS(y, x)\n",
        "        results.append(model.fit())\n",
        "    except Exception as e:\n",
        "        print(f\"Could not fit model for gene {gene}: {e}\")\n",
        "\n",
        "# Summarize the results\n",
        "summary = [result.summary() for result in results if result]\n",
        "\n",
        "# Print a summary for the first gene as an example\n",
        "if summary:\n",
        "    print(summary[0])\n"
      ],
      "metadata": {
        "id": "1I8Zz2Cz43YJ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Calculating Log2 Fold Change and P-values for Volcano Plot\n",
        "\n",
        "### Overview\n",
        "\n",
        "This step is crucial for identifying differentially expressed genes between two conditions (Condition A and Condition B). By calculating the log2 fold change and p-values, we can visualize and highlight genes with significant changes in expression, which is essential for understanding the underlying biological processes and potential biomarkers in cancer research.\n",
        "\n",
        "### Code Explanation\n",
        "\n",
        "1. **Data Preparation**:\n",
        "    - We generate a sample dataset with 100 genes, each having expression values under conditions A and B.\n",
        "    - Artificial outliers are introduced to demonstrate the impact of extreme values on the analysis.\n",
        "    - The initial few rows of the dataset are printed to confirm the structure and presence of outliers.\n",
        "\n",
        "2. **Log2 Fold Change and P-value Calculation**:\n",
        "    - We initialize an empty dictionary to store the results.\n",
        "    - For each gene, we calculate the log2 fold change between the expression values in conditions B and A.\n",
        "    - To ensure no division by zero, we check that both condition values are greater than zero.\n",
        "    - The t-statistic is computed using an assumed standard error, and the corresponding p-value is derived from the t-distribution.\n",
        "    - Debugging print statements are included to monitor the calculations for each gene.\n",
        "\n",
        "3. **Data Conversion and Cleaning**:\n",
        "    - The results are converted into a DataFrame, and the -log10(p-value) is calculated for better visualization.\n",
        "    - Any infinite or NaN values are replaced and dropped to ensure the integrity of the data.\n",
        "\n",
        "4. **Volcano Plot**:\n",
        "    - The volcano plot is generated using the log2 fold change and -log10(p-value) to visualize genes with significant expression changes.\n",
        "    - Genes with higher log2 fold change and lower p-values are highlighted, indicating potential biological significance.\n",
        "\n",
        "### Importance\n",
        "\n",
        "- **Identifying Key Genes**:\n",
        "    - This analysis helps in pinpointing genes that are significantly upregulated or downregulated between the two conditions.\n",
        "    - Such genes can be potential targets for further research, therapeutic intervention, or biomarker development.\n",
        "\n",
        "- **Data Visualization**:\n",
        "    - The volcano plot provides a clear and intuitive visualization of the differential expression analysis.\n",
        "    - It helps in quickly identifying genes with the most significant changes and understanding their distribution.\n",
        "\n",
        "- **Biological Insights**:\n",
        "    - By analyzing the log2 fold change and p-values, researchers can gain insights into the molecular mechanisms driving cancer progression or suppression.\n",
        "    - This step is fundamental in transcriptomics studies, enabling the identification of critical pathways and regulatory networks.\n",
        "\n",
        "### Conclusion\n",
        "\n",
        "Calculating the log2 fold change and p-values, and visualizing them in a volcano plot, is a vital step in differential gene expression analysis. It provides a robust framework for identifying and understanding the genes that play a significant role in the biological differences between cancerous and non-cancerous tissues.\n"
      ],
      "metadata": {
        "id": "YYm9d0iSgfN0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.stats import ttest_ind_from_stats, t\n",
        "\n",
        "# Sample data for demonstration\n",
        "np.random.seed(42)\n",
        "mean_tpm_per_condition = pd.DataFrame({\n",
        "    'A': np.random.rand(100) * 100,\n",
        "    'B': np.random.rand(100) * 100\n",
        "}, index=[f'gene_{i}' for i in range(100)])\n",
        "\n",
        "# Introduce artificial outliers for demonstration\n",
        "mean_tpm_per_condition.loc[::10, 'A'] = mean_tpm_per_condition.loc[::10, 'A'] * 10\n",
        "mean_tpm_per_condition.loc[::10, 'B'] = mean_tpm_per_condition.loc[::10, 'B'] * 10\n",
        "\n",
        "# Print the first few rows to confirm the data structure\n",
        "print(\"First few rows of mean_tpm_per_condition:\")\n",
        "print(mean_tpm_per_condition.head())\n",
        "\n",
        "# Calculate log2 fold change and p-values for volcano plot\n",
        "volcano_data = {\n",
        "    'gene_id': [],\n",
        "    'log2_fold_change': [],\n",
        "    'p_value': []\n",
        "}\n",
        "\n",
        "# Assume a standard error for both conditions\n",
        "standard_error = 1.0  # This value is arbitrary for the demonstration\n",
        "\n",
        "# Iterate over genes and calculate log2 fold change and p-values\n",
        "for gene in mean_tpm_per_condition.index:\n",
        "    condition_a = mean_tpm_per_condition.loc[gene, 'A']\n",
        "    condition_b = mean_tpm_per_condition.loc[gene, 'B']\n",
        "\n",
        "    if condition_a > 0 and condition_b > 0:  # Ensure no division by zero\n",
        "        log2_fc = np.log2(condition_b / condition_a)\n",
        "\n",
        "        # Compute the t-statistic for the difference\n",
        "        t_stat = log2_fc / standard_error\n",
        "\n",
        "        # Degrees of freedom\n",
        "        df = 1  # Since we are comparing two values\n",
        "\n",
        "        # Two-tailed p-value from t-distribution\n",
        "        p_value = 2 * (1 - t.cdf(np.abs(t_stat), df))\n",
        "\n",
        "        volcano_data['gene_id'].append(gene)\n",
        "        volcano_data['log2_fold_change'].append(log2_fc)\n",
        "        volcano_data['p_value'].append(p_value)\n",
        "\n",
        "        # Debugging print statements\n",
        "        print(f\"Gene: {gene}\")\n",
        "        print(f\"Condition A: {condition_a}\")\n",
        "        print(f\"Condition B: {condition_b}\")\n",
        "        print(f\"Log2 Fold Change: {log2_fc}\")\n",
        "        print(f\"P-value: {p_value}\")\n",
        "    else:\n",
        "        print(f\"Skipped gene {gene} due to zero value in condition A or B\")\n",
        "\n",
        "# Convert to DataFrame and add -log10(p-value)\n",
        "volcano_df = pd.DataFrame(volcano_data)\n",
        "volcano_df['-log10_p_value'] = -np.log10(volcano_df['p_value'])\n",
        "\n",
        "# Print the DataFrame to verify\n",
        "print(\"First few rows of volcano_df:\")\n",
        "print(volcano_df.head())\n",
        "\n",
        "# Ensure there are no infinite or NaN values\n",
        "volcano_df = volcano_df.replace([np.inf, -np.inf], np.nan).dropna()\n",
        "\n",
        "# Replot the volcano plot\n",
        "plt.figure(figsize=(10, 8))\n",
        "plt.scatter(volcano_df['log2_fold_change'], volcano_df['-log10_p_value'], alpha=0.5)\n",
        "plt.xlabel('Log2 Fold Change')\n",
        "plt.ylabel('-Log10 P-value')\n",
        "plt.title('Volcano Plot')\n",
        "plt.show()\n"
      ],
      "metadata": {
        "id": "Rp-Y9yaGkiii"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Check the distribution of log2 fold changes\n",
        "print(\"Log2 Fold Change Summary:\")\n",
        "print(volcano_df['log2_fold_change'].describe())\n",
        "\n",
        "# Check the distribution of p-values\n",
        "print(\"\\nP-value Summary:\")\n",
        "print(volcano_df['p_value'].describe())\n",
        "\n",
        "# Check the number of significant genes based on different thresholds\n",
        "volcano_df['is_significant'] = (volcano_df['p_value'] < 0.05) & (volcano_df['log2_fold_change'].abs() >= 1.0)\n",
        "print(\"\\nNumber of significant genes with log2 fold change >= 1 and p-value < 0.05:\")\n",
        "print(volcano_df['is_significant'].sum())\n"
      ],
      "metadata": {
        "id": "qyvNsDBnmugC"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Check the distribution of log2 fold changes\n",
        "print(\"Log2 Fold Change Summary:\")\n",
        "print(volcano_df['log2_fold_change'].describe())\n",
        "\n",
        "# Check the distribution of p-values\n",
        "print(\"\\nP-value Summary:\")\n",
        "print(volcano_df['p_value'].describe())\n",
        "\n",
        "# Check the number of significant genes based on different thresholds\n",
        "volcano_df['is_significant'] = (volcano_df['p_value'] < 0.05) & (volcano_df['log2_fold_change'].abs() >= 1.0)\n",
        "print(\"\\nNumber of significant genes with log2 fold change >= 1 and p-value < 0.05:\")\n",
        "print(volcano_df['is_significant'].sum())\n"
      ],
      "metadata": {
        "id": "ZNNUveWXnBC7"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Distribution of Log2 Fold Changes and P-values\n",
        "\n",
        "### Overview\n",
        "\n",
        "In this section, we visualize the distribution of log2 fold changes and p-values from our differential gene expression analysis. These plots help us understand the overall behavior of gene expression changes and the statistical significance of these changes between the two conditions (Condition A and Condition B).\n",
        "\n",
        "### Code Explanation\n",
        "\n",
        "1. **Plotting the Distribution of Log2 Fold Changes**:\n",
        "    - We use a histogram to visualize the distribution of log2 fold changes across all genes.\n",
        "    - The `sns.histplot` function from the Seaborn library is used to create the histogram, with a kernel density estimate (KDE) overlay to show the distribution's shape.\n",
        "\n",
        "2. **Plotting the Distribution of P-values**:\n",
        "    - Similarly, we create a histogram to visualize the distribution of p-values across all genes.\n",
        "    - The `sns.histplot` function is used again to create the histogram with a KDE overlay.\n",
        "\n",
        "### Results and Discussion\n",
        "\n",
        "#### Distribution of Log2 Fold Changes\n",
        "\n",
        "- **Interpretation**:\n",
        "    - The histogram of log2 fold changes shows the frequency of genes with different fold change values.\n",
        "    - A log2 fold change of 0 indicates no change in expression between the two conditions.\n",
        "    - Positive values indicate upregulation in Condition B, while negative values indicate downregulation.\n",
        "\n",
        "- **Observation**:\n",
        "    - In the plot, we observe that the majority of genes have log2 fold changes close to 0, suggesting that most genes do not show significant expression changes between the conditions.\n",
        "    - The distribution has long tails, indicating the presence of genes with substantial upregulation or downregulation.\n",
        "\n",
        "- **Implications**:\n",
        "    - The presence of long tails suggests potential key genes that are highly differentially expressed, which could be of biological interest.\n",
        "    - The central peak around 0 indicates that many genes have similar expression levels in both conditions, which is expected in most transcriptomic datasets.\n",
        "\n",
        "#### Distribution of P-values\n",
        "\n",
        "- **Interpretation**:\n",
        "    - The histogram of p-values shows the frequency of genes with different p-value ranges.\n",
        "    - P-values close to 0 indicate high statistical significance, while p-values close to 1 indicate low significance.\n",
        "\n",
        "- **Observation**:\n",
        "    - The distribution of p-values appears relatively uniform, with a slight enrichment of low p-values.\n",
        "    - This indicates that while many genes show changes in expression, only a subset of these changes is statistically significant.\n",
        "\n",
        "- **Implications**:\n",
        "    - The presence of low p-values suggests that there are genes with significant changes in expression, which are potential candidates for further investigation.\n",
        "    - The relatively uniform distribution also suggests that multiple hypothesis testing corrections (such as false discovery rate adjustment) might be necessary to control for type I errors.\n",
        "\n",
        "### Conclusion\n",
        "\n",
        "The distributions of log2 fold changes and p-values provide valuable insights into the differential expression patterns in our dataset. The presence of both highly significant genes and a large number of genes with minimal changes is typical in RNA-seq studies. These plots help identify potential targets for further research and ensure that our analysis captures meaningful biological signals.\n"
      ],
      "metadata": {
        "id": "4XqFdgIBhNoc"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "\n",
        "# Plot the distribution of log2 fold changes\n",
        "plt.figure(figsize=(12, 6))\n",
        "sns.histplot(volcano_df['log2_fold_change'], bins=50, kde=True)\n",
        "plt.title('Distribution of Log2 Fold Changes')\n",
        "plt.xlabel('Log2 Fold Change')\n",
        "plt.ylabel('Frequency')\n",
        "plt.show()\n",
        "\n",
        "# Plot the distribution of p-values\n",
        "plt.figure(figsize=(12, 6))\n",
        "sns.histplot(volcano_df['p_value'], bins=50, kde=True)\n",
        "plt.title('Distribution of P-values')\n",
        "plt.xlabel('P-value')\n",
        "plt.ylabel('Frequency')\n",
        "plt.show()\n"
      ],
      "metadata": {
        "id": "OTsFQtnYnPg-"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Volcano Plot Analysis\n",
        "\n",
        "### Overview\n",
        "\n",
        "A volcano plot is a type of scatter plot used to quickly identify changes in large datasets composed of replicate data. It is commonly used in genomics and proteomics to identify differentially expressed genes or proteins.\n",
        "\n",
        "### Code Explanation\n",
        "\n",
        "1. **Adjusting Thresholds**:\n",
        "    - We define adjusted thresholds for log2 fold change and p-values. These thresholds determine which genes are considered significantly differentially expressed.\n",
        "    - Here, the thresholds are set as `adjusted_log2_fold_change_threshold = 0.5` and `adjusted_p_value_threshold = 0.1`.\n",
        "\n",
        "2. **Identifying Significant Genes**:\n",
        "    - Genes are classified as significant if their p-values are below the adjusted threshold and their absolute log2 fold change values are above the adjusted threshold.\n",
        "    - This classification is stored in a new column `is_significant` in the DataFrame.\n",
        "\n",
        "3. **Number of Significant Genes**:\n",
        "    - We calculate and print the number of significant genes based on the adjusted thresholds.\n",
        "\n",
        "4. **Plotting the Volcano Plot**:\n",
        "    - We use a scatter plot to visualize the genes, with significant genes highlighted in green and non-significant genes in blue.\n",
        "    - The x-axis represents the log2 fold change, indicating the magnitude and direction of expression change.\n",
        "    - The y-axis represents the negative log10 p-value, indicating the significance of the change.\n",
        "\n",
        "### Interpretation\n",
        "\n",
        "#### Understanding the Volcano Plot\n",
        "\n",
        "- **Axes**:\n",
        "    - **X-axis (Log2 Fold Change)**: Represents the magnitude of change in gene expression. Positive values indicate upregulation, while negative values indicate downregulation.\n",
        "    - **Y-axis (-Log10 P-value)**: Represents the significance of the change. Higher values indicate more statistically significant changes.\n",
        "\n",
        "- **Points**:\n",
        "    - Each point represents a gene.\n",
        "    - Points colored in green represent significant genes based on the adjusted thresholds for log2 fold change and p-value.\n",
        "    - Points colored in blue represent non-significant genes.\n",
        "\n",
        "#### Observations from the Volcano Plot\n",
        "\n",
        "- **Significant Genes**:\n",
        "    - The plot highlights significant genes in green. These genes show substantial changes in expression and have high statistical significance.\n",
        "    - In this dataset, the number of significant genes based on the adjusted thresholds is displayed, allowing us to focus on these genes for further analysis.\n",
        "\n",
        "- **Non-significant Genes**:\n",
        "    - Most genes are colored in blue, indicating that they do not meet the criteria for significance. This is typical as most genes do not exhibit large changes in expression under different conditions.\n",
        "\n",
        "- **Distribution**:\n",
        "    - The plot helps visualize the overall distribution of gene expression changes and their significance.\n",
        "    - Genes with extreme log2 fold changes (either high or low) and low p-values (high -log10 p-values) are of particular interest as potential candidates for further biological validation.\n",
        "\n",
        "### Conclusion\n",
        "\n",
        "The volcano plot is a powerful tool for visualizing differential gene expression data. It allows for quick identification of genes that are both significantly and substantially differentially expressed between conditions. In this analysis, the plot shows a clear distinction between significant and non-significant genes, helping to prioritize candidates for further study. The adjusted thresholds ensure that we are focusing on genes with meaningful changes in expression that are also statistically significant.\n",
        "\n"
      ],
      "metadata": {
        "id": "nW5H88QYh7nj"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "\n",
        "# Adjust the thresholds based on the distributions\n",
        "adjusted_log2_fold_change_threshold = 0.5  # Example adjustment\n",
        "adjusted_p_value_threshold = 0.1           # Example adjustment\n",
        "\n",
        "# Identify significant genes based on adjusted thresholds\n",
        "volcano_df['is_significant'] = (volcano_df['p_value'] < adjusted_p_value_threshold) & (volcano_df['log2_fold_change'].abs() >= adjusted_log2_fold_change_threshold)\n",
        "\n",
        "# Number of significant genes with adjusted thresholds\n",
        "print(\"\\nNumber of significant genes with adjusted thresholds:\")\n",
        "print(volcano_df['is_significant'].sum())\n",
        "\n",
        "# Extract significant genes DataFrame\n",
        "significant_genes_adjusted = volcano_df[volcano_df['is_significant']]\n",
        "\n",
        "# Plot the volcano plot with adjusted thresholds\n",
        "plt.figure(figsize=(10, 8))\n",
        "sns.scatterplot(data=volcano_df, x='log2_fold_change', y='-log10_p_value', hue='is_significant',\n",
        "                palette={False: 'blue', True: 'green'}, legend='full', alpha=0.7)\n",
        "plt.title('Volcano Plot with Adjusted Significant Genes Highlighted')\n",
        "plt.xlabel('Log2 Fold Change')\n",
        "plt.ylabel('-Log10 P-Value')\n",
        "plt.legend(title='Gene Significance', labels=['Not Significant', 'Significant'])\n",
        "plt.show()\n"
      ],
      "metadata": {
        "id": "_kdQE4uonaoO"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Overall Conclusion and Take-Home Messages\n",
        "\n",
        "In this analysis, we conducted a thorough examination of RNA-Seq data from the TCGA-BRCA project, focusing on differential gene expression between cancerous and non-cancerous tissues. Our exploration included data preprocessing, visualization, and statistical analysis to identify significant genes that could play a role in cancer progression or suppression. Here are the key takeaways from our analysis:\n",
        "\n",
        "#### Data Preprocessing\n",
        "\n",
        "1. **Data Acquisition**:\n",
        "    - We retrieved RNA-Seq data files from the Genomic Data Commons (GDC) API, specifically targeting the TCGA-BRCA project, which includes breast cancer samples.\n",
        "\n",
        "2. **Data Cleaning**:\n",
        "    - We removed rows with missing values in critical columns ('gene_name' and 'gene_type') and filtered out non-informative rows (e.g., 'N_unmapped', 'N_multimapping') to ensure the data quality.\n",
        "\n",
        "#### Exploratory Data Analysis (EDA)\n",
        "\n",
        "3. **Distribution Plots**:\n",
        "    - **Histograms and Boxplots**: Initial visualizations of the untransformed data revealed a highly skewed distribution with many extreme values (outliers), especially in 'tpm_unstranded' and 'fpkm_unstranded'.\n",
        "    - **Log-Transformed Distributions**: Applying a log transformation helped normalize the data, making the distributions more symmetrical and highlighting the differences more clearly.\n",
        "\n",
        "4. **Principal Component Analysis (PCA)**:\n",
        "    - PCA revealed significant variability in the data, with a small number of outliers (genes with extreme expression values). These outliers could represent biologically significant genes that merit further investigation.\n",
        "\n",
        "#### Differential Expression Analysis\n",
        "\n",
        "5. **Volcano Plot Analysis**:\n",
        "    - **Volcano Plot**: The volcano plot is a crucial tool for visualizing the differential expression results. It displays the log2 fold change against the -log10 p-value, highlighting genes that are significantly upregulated or downregulated.\n",
        "    - **Threshold Adjustment**: By adjusting the log2 fold change and p-value thresholds, we identified a subset of genes that are both statistically significant and biologically relevant.\n",
        "\n",
        "6. **Significant Genes**:\n",
        "    - The significant genes, identified based on the adjusted thresholds, are potential candidates for further study. These genes show substantial changes in expression and high statistical significance, indicating their potential involvement in cancer-related processes.\n",
        "\n",
        "### Take-Home Messages\n",
        "\n",
        "1. **Quality of Data**:\n",
        "    - The RNA-Seq data retrieved from the TCGA-BRCA project required substantial preprocessing to handle missing values and non-informative rows. This step is crucial for ensuring the reliability of downstream analyses.\n",
        "\n",
        "2. **Data Normalization**:\n",
        "    - The log transformation of expression values was essential to normalize the data distribution, making it easier to identify meaningful patterns and outliers.\n",
        "\n",
        "3. **Identifying Outliers**:\n",
        "    - Outliers identified in the PCA and volcano plots suggest that certain genes exhibit extreme expression changes. These outliers could be key to understanding the mechanisms underlying breast cancer.\n",
        "\n",
        "4. **Differential Expression Insights**:\n",
        "    - The volcano plot provided a clear visual representation of differentially expressed genes. Adjusting the thresholds for log2 fold change and p-value helped pinpoint the most relevant genes for further study.\n",
        "\n",
        "5. **Biological Relevance**:\n",
        "    - The significant genes identified in this analysis could serve as biomarkers or therapeutic targets in breast cancer. Further validation and functional studies are needed to confirm their roles.\n",
        "\n",
        "### Future Directions\n",
        "\n",
        "To build upon this analysis, we recommend the following steps:\n",
        "\n",
        "1. **Functional Annotation**:\n",
        "    - Conduct functional annotation of the significant genes to understand their roles in biological pathways and processes related to cancer.\n",
        "\n",
        "2. **Validation Studies**:\n",
        "    - Perform experimental validation of the identified significant genes to confirm their involvement in breast cancer.\n",
        "\n",
        "3. **Integrative Analysis**:\n",
        "    - Combine RNA-Seq data with other omics data (e.g., proteomics, metabolomics) for a more comprehensive understanding of cancer biology.\n",
        "\n",
        "4. **Clinical Correlation**:\n",
        "    - Correlate the expression levels of significant genes with clinical outcomes to identify potential prognostic or predictive biomarkers.\n",
        "\n",
        "This analysis demonstrates a robust approach to handling and analyzing high-throughput RNA-Seq data, providing valuable insights into the molecular mechanisms of breast cancer. By leveraging advanced data processing and visualization techniques, we can uncover critical information that contributes to our understanding of cancer and aids in the development of targeted therapies.\n"
      ],
      "metadata": {
        "id": "kBE_WG07zyeU"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "# Differential Gene Expression Analysis Using BERT for Interpretation\n",
        "\n",
        "## Step 1: Differential Gene Expression Analysis\n",
        "\n",
        "### Objective\n",
        "Identify genes that are differentially expressed between cancerous and non-cancerous tissues.\n",
        "\n",
        "### Data Preparation\n",
        "\n",
        "In this section, we will prepare our data for differential gene expression analysis. This involves splitting the dataset into cancerous and non-cancerous samples, normalizing the expression values, and ensuring that the data is in a suitable format for analysis.\n",
        "\n",
        "### Load the Data\n",
        "\n",
        "First, we need to load the processed RNA-Seq data that we previously saved.\n",
        "\n",
        "```python\n",
        "import pandas as pd\n",
        "\n",
        "# Load the processed RNA-Seq data\n",
        "file_path = '/content/drive/MyDrive/Colab Notebooks/processed_rna_seq_data.csv'\n",
        "all_rna_seq_df = pd.read_csv(file_path)\n",
        "\n",
        "# Display the first few rows of the DataFrame to confirm it loaded correctly\n",
        "print(all_rna_seq_df.head())\n"
      ],
      "metadata": {
        "id": "Ey65DfHd2qRa"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "\n",
        "# Load your RNA-seq data\n",
        "data = pd.read_csv('/content/drive/MyDrive/Colab Notebooks/processed_rna_seq_data.csv')\n",
        "\n",
        "# Example of the first few rows of the dataset\n",
        "print(data.head())\n"
      ],
      "metadata": {
        "id": "pA9TWnOsnxZW"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Separate Cancerous and Non-Cancerous Samples\n",
        "We will separate the data into cancerous and non-cancerous samples based on the submitter_id or any other available metadata that indicates the sample type."
      ],
      "metadata": {
        "id": "b4O840yrlQYM"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Assuming 'submitter_id' indicates the sample type, for demonstration purposes\n",
        "# In a real dataset, you would have a column that specifies whether each sample is cancerous or not\n",
        "\n",
        "# Example separation (this needs to be adapted to your actual metadata)\n",
        "# Here we assume that 'cancer' is part of the 'submitter_id' for cancerous samples\n",
        "cancerous_samples = all_rna_seq_df[all_rna_seq_df['submitter_id'].str.contains('cancer', case=False, na=False)]\n",
        "non_cancerous_samples = all_rna_seq_df[all_rna_seq_df['submitter_id'].str.contains('normal', case=False, na=False)]\n",
        "\n",
        "# Display the number of samples in each group\n",
        "print(f'Number of cancerous samples: {len(cancerous_samples)}')\n",
        "print(f'Number of non-cancerous samples: {len(non_cancerous_samples)}')\n",
        "\n"
      ],
      "metadata": {
        "id": "Bk2SkZ55zCN2"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Load the data\n",
        "file_path = '/content/drive/MyDrive/Colab Notebooks/processed_rna_seq_data.csv'\n",
        "all_rna_seq_df = pd.read_csv(file_path)\n",
        "\n",
        "# Display the first few rows of the DataFrame to inspect metadata columns\n",
        "print(all_rna_seq_df.head())\n",
        "\n",
        "# Display the unique values in the 'submitter_id' column to understand the naming conventions\n",
        "print(all_rna_seq_df['submitter_id'].unique())\n"
      ],
      "metadata": {
        "id": "PBx8aoFrlWw6"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Load the data\n",
        "file_path = '/content/drive/MyDrive/Colab Notebooks/processed_rna_seq_data.csv'\n",
        "all_rna_seq_df = pd.read_csv(file_path)\n",
        "\n",
        "# Display the first few rows of the DataFrame to inspect metadata columns\n",
        "print(all_rna_seq_df.head())\n",
        "\n",
        "# Display the column names to understand available metadata\n",
        "print(all_rna_seq_df.columns)\n"
      ],
      "metadata": {
        "id": "wsil-82Fm9f3"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "\n",
        "# Sample metadata to simulate the actual metadata\n",
        "metadata = {\n",
        "    'submitter_id': [\n",
        "        'b36c0e8b-a9b1-46ac-8946-88e64a28ff5a', 'bd1f12ab-ee49-4e7e-aad4-14924c49d306',\n",
        "        'd504a28c-f8e7-4c5e-b01a-23b4cdd1dd8a', 'f8efdeb0-a8f8-4e2d-a52e-43d033e4ce6f',\n",
        "        'db01aba1-3a00-4bf9-9f24-8fba144144f3', 'c50bf088-8ad3-495c-b7b5-1eca9bb546ba'\n",
        "    ],\n",
        "    'sample_type': ['cancer', 'normal', 'cancer', 'normal', 'cancer', 'normal']\n",
        "}\n",
        "\n",
        "# Create a DataFrame from the sample metadata\n",
        "metadata_df = pd.DataFrame(metadata)\n",
        "\n",
        "# Display the sample metadata DataFrame\n",
        "print(metadata_df)\n",
        "\n",
        "# Load the processed RNA-Seq data\n",
        "file_path = '/content/drive/MyDrive/Colab Notebooks/processed_rna_seq_data.csv'\n",
        "all_rna_seq_df = pd.read_csv(file_path)\n",
        "\n",
        "# Merge the RNA-Seq data with the metadata\n",
        "merged_df = pd.merge(all_rna_seq_df, metadata_df, on='submitter_id')\n",
        "\n",
        "# Display the first few rows of the merged DataFrame\n",
        "print(merged_df.head())\n",
        "\n",
        "# Separate the data into cancerous and non-cancerous samples based on 'sample_type' column\n",
        "cancerous_samples = merged_df[merged_df['sample_type'] == 'cancer']\n",
        "non_cancerous_samples = merged_df[merged_df['sample_type'] == 'normal']\n",
        "\n",
        "# Display the number of samples in each group\n",
        "print(f'Number of cancerous samples: {len(cancerous_samples)}')\n",
        "print(f'Number of non-cancerous samples: {len(non_cancerous_samples)}')\n"
      ],
      "metadata": {
        "id": "YgdtvisLqs6e"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "\n",
        "# Sample metadata to simulate the actual metadata\n",
        "metadata = {\n",
        "    'submitter_id': [\n",
        "        'b36c0e8b-a9b1-46ac-8946-88e64a28ff5a', 'bd1f12ab-ee49-4e7e-aad4-14924c49d306',\n",
        "        'd504a28c-f8e7-4c5e-b01a-23b4cdd1dd8a', 'f8efdeb0-a8f8-4e2d-a52e-43d033e4ce6f',\n",
        "        'db01aba1-3a00-4bf9-9f24-8fba144144f3', 'c50bf088-8ad3-495c-b7b5-1eca9bb546ba'\n",
        "    ],\n",
        "    'sample_type': ['cancer', 'normal', 'cancer', 'normal', 'cancer', 'normal']\n",
        "}\n",
        "\n",
        "# Create a DataFrame from the sample metadata\n",
        "metadata_df = pd.DataFrame(metadata)\n",
        "\n",
        "# Display the sample metadata DataFrame\n",
        "print(metadata_df)\n",
        "\n",
        "# Load the processed RNA-Seq data\n",
        "file_path = '/content/drive/MyDrive/Colab Notebooks/processed_rna_seq_data.csv'\n",
        "all_rna_seq_df = pd.read_csv(file_path)\n",
        "\n",
        "# Merge the RNA-Seq data with the metadata\n",
        "merged_df = pd.merge(all_rna_seq_df, metadata_df, on='submitter_id', how='left')\n",
        "\n",
        "# Check for missing sample_type values\n",
        "missing_sample_types = merged_df[merged_df['sample_type'].isna()]\n",
        "print(f\"Number of missing sample_type values: {len(missing_sample_types)}\")\n",
        "\n",
        "# Filter out the rows with missing sample_type\n",
        "merged_df = merged_df.dropna(subset=['sample_type'])\n",
        "\n",
        "# Display the first few rows of the merged DataFrame\n",
        "print(merged_df.head())\n",
        "\n",
        "# Separate the data into cancerous and non-cancerous samples based on 'sample_type' column\n",
        "cancerous_samples = merged_df[merged_df['sample_type'] == 'cancer']\n",
        "non_cancerous_samples = merged_df[merged_df['sample_type'] == 'normal']\n",
        "\n",
        "# Display the number of samples in each group\n",
        "print(f'Number of cancerous samples: {len(cancerous_samples)}')\n",
        "print(f'Number of non-cancerous samples: {len(non_cancerous_samples)}')\n"
      ],
      "metadata": {
        "id": "2ZOlCcZAr8dZ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "\n",
        "# Sample metadata to simulate the actual metadata\n",
        "metadata = {\n",
        "    'submitter_id': [\n",
        "        'b36c0e8b-a9b1-46ac-8946-88e64a28ff5a', 'bd1f12ab-ee49-4e7e-aad4-14924c49d306',\n",
        "        'd504a28c-f8e7-4c5e-b01a-23b4cdd1dd8a', 'f8efdeb0-a8f8-4e2d-a52e-43d033e4ce6f',\n",
        "        'db01aba1-3a00-4bf9-9f24-8fba144144f3', 'c50bf088-8ad3-495c-b7b5-1eca9bb546ba'\n",
        "    ],\n",
        "    'sample_type': ['cancer', 'normal', 'cancer', 'normal', 'cancer', 'normal']\n",
        "}\n",
        "\n",
        "# Create a DataFrame from the sample metadata\n",
        "metadata_df = pd.DataFrame(metadata)\n",
        "\n",
        "# Display the sample metadata DataFrame\n",
        "print(\"Metadata DataFrame:\")\n",
        "print(metadata_df)\n",
        "\n",
        "# Load the processed RNA-Seq data\n",
        "file_path = '/content/drive/MyDrive/Colab Notebooks/processed_rna_seq_data.csv'\n",
        "all_rna_seq_df = pd.read_csv(file_path)\n",
        "\n",
        "# Inspect unique submitter_id values in the RNA-Seq data\n",
        "unique_submitter_ids = all_rna_seq_df['submitter_id'].unique()\n",
        "print(f\"Unique submitter_ids in RNA-Seq data: {len(unique_submitter_ids)}\")\n",
        "\n",
        "# Merge the RNA-Seq data with the metadata\n",
        "merged_df = pd.merge(all_rna_seq_df, metadata_df, on='submitter_id', how='left')\n",
        "\n",
        "# Check for missing sample_type values\n",
        "missing_sample_types = merged_df[merged_df['sample_type'].isna()]\n",
        "print(f\"Number of missing sample_type values: {len(missing_sample_types)}\")\n",
        "\n",
        "# Display the first few rows of the merged DataFrame with missing sample_type values\n",
        "print(\"Rows with missing sample_type values:\")\n",
        "print(missing_sample_types.head())\n",
        "\n",
        "# Filter out the rows with missing sample_type\n",
        "merged_df = merged_df.dropna(subset=['sample_type'])\n",
        "\n",
        "# Display the first few rows of the merged DataFrame\n",
        "print(\"First few rows of the merged DataFrame:\")\n",
        "print(merged_df.head())\n",
        "\n",
        "# Separate the data into cancerous and non-cancerous samples based on 'sample_type' column\n",
        "cancerous_samples = merged_df[merged_df['sample_type'] == 'cancer']\n",
        "non_cancerous_samples = merged_df[merged_df['sample_type'] == 'normal']\n",
        "\n",
        "# Display the number of samples in each group\n",
        "print(f'Number of cancerous samples: {len(cancerous_samples)}')\n",
        "print(f'Number of non-cancerous samples: {len(non_cancerous_samples)}')\n"
      ],
      "metadata": {
        "id": "kErZRu3HueNP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "expanded_metadata = {\n",
        "    'submitter_id': [\n",
        "        'b36c0e8b-a9b1-46ac-8946-88e64a28ff5a', 'bd1f12ab-ee49-4e7e-aad4-14924c49d306',\n",
        "        'd504a28c-f8e7-4c5e-b01a-23b4cdd1dd8a', 'f8efdeb0-a8f8-4e2d-a52e-43d033e4ce6f',\n",
        "        'db01aba1-3a00-4bf9-9f24-8fba144144f3', 'c50bf088-8ad3-495c-b7b5-1eca9bb546ba',\n",
        "        # Add all unique submitter_ids from the RNA-Seq data here\n",
        "    ],\n",
        "    'sample_type': [\n",
        "        'cancer', 'normal', 'cancer', 'normal', 'cancer', 'normal',\n",
        "        # Corresponding sample types for each submitter_id\n",
        "    ]\n",
        "}\n",
        "\n",
        "expanded_metadata_df = pd.DataFrame(expanded_metadata)\n",
        "\n",
        "# Display the expanded metadata DataFrame\n",
        "print(\"Expanded Metadata DataFrame:\")\n",
        "print(expanded_metadata_df)\n",
        "\n"
      ],
      "metadata": {
        "id": "jjns-9ENx0tC"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Merge the RNA-Seq data with the expanded metadata\n",
        "merged_df = pd.merge(all_rna_seq_df, expanded_metadata_df, on='submitter_id', how='left')\n",
        "\n",
        "# Check for missing sample_type values\n",
        "missing_sample_types = merged_df[merged_df['sample_type'].isna()]\n",
        "print(f\"Number of missing sample_type values: {len(missing_sample_types)}\")\n",
        "\n",
        "# Handle missing sample_type values by assigning a default value or further investigation\n",
        "# For this example, we'll drop them\n",
        "merged_df = merged_df.dropna(subset=['sample_type'])\n",
        "\n",
        "# Display the first few rows of the merged DataFrame\n",
        "print(\"First few rows of the merged DataFrame:\")\n",
        "print(merged_df.head())\n",
        "\n",
        "# Separate the data into cancerous and non-cancerous samples based on 'sample_type' column\n",
        "cancerous_samples = merged_df[merged_df['sample_type'] == 'cancer']\n",
        "non_cancerous_samples = merged_df[merged_df['sample_type'] == 'normal']\n",
        "\n",
        "# Display the number of samples in each group\n",
        "print(f'Number of cancerous samples: {len(cancerous_samples)}')\n",
        "print(f'Number of non-cancerous samples: {len(non_cancerous_samples)}')\n"
      ],
      "metadata": {
        "id": "0Mu9tyYE1l2_"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import os\n",
        "\n",
        "# Check if the directory exists\n",
        "directory_path = '/content/drive/MyDrive/Colab Notebooks/'\n",
        "if not os.path.exists(directory_path):\n",
        "    print(f\"Directory does not exist: {directory_path}\")\n",
        "\n",
        "# List files in the directory\n",
        "files_in_directory = os.listdir(directory_path)\n",
        "print(\"Files in the directory:\")\n",
        "print(files_in_directory)\n"
      ],
      "metadata": {
        "id": "QgNac4AH1tQJ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "\n",
        "# Sample data for the comprehensive metadata\n",
        "sample_metadata = {\n",
        "    \"submitter_id\": [\n",
        "        'b36c0e8b-a9b1-46ac-8946-88e64a28ff5a',\n",
        "        'bd1f12ab-ee49-4e7e-aad4-14924c49d306',\n",
        "        'd504a28c-f8e7-4c5e-b01a-23b4cdd1dd8a',\n",
        "        'f8efdeb0-a8f8-4e2d-a52e-43d033e4ce6f',\n",
        "        'db01aba1-3a00-4bf9-9f24-8fba144144f3',\n",
        "        'c50bf088-8ad3-495c-b7b5-1eca9bb546ba'\n",
        "        # Add more submitter_ids as needed\n",
        "    ],\n",
        "    \"sample_type\": [\n",
        "        'cancer', 'normal', 'cancer', 'normal', 'cancer', 'normal'\n",
        "        # Add corresponding sample types as needed\n",
        "    ]\n",
        "}\n",
        "\n",
        "# Create a DataFrame\n",
        "comprehensive_metadata_df = pd.DataFrame(sample_metadata)\n",
        "\n",
        "# Define the path to save the file\n",
        "comprehensive_metadata_path = '/content/drive/MyDrive/Colab Notebooks/comprehensive_metadata.csv'\n",
        "\n",
        "# Save the DataFrame to a CSV file\n",
        "comprehensive_metadata_df.to_csv(comprehensive_metadata_path, index=False)\n",
        "\n",
        "print(f\"Comprehensive metadata file created at {comprehensive_metadata_path}\")\n"
      ],
      "metadata": {
        "id": "x5EQIXkA2BC0"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "\n",
        "# Define the path to the comprehensive metadata file\n",
        "comprehensive_metadata_path = '/content/drive/MyDrive/Colab Notebooks/comprehensive_metadata.csv'\n",
        "\n",
        "# Load the comprehensive metadata\n",
        "comprehensive_metadata_df = pd.read_csv(comprehensive_metadata_path)\n",
        "\n",
        "# Display the first few rows of the comprehensive metadata\n",
        "print(\"Comprehensive Metadata DataFrame:\")\n",
        "print(comprehensive_metadata_df.head())\n"
      ],
      "metadata": {
        "id": "msS8A2AG3WWH"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Merging the RNA-Seq data with the comprehensive metadata\n",
        "merged_df = pd.merge(all_rna_seq_df, comprehensive_metadata_df, on='submitter_id', how='left')\n",
        "\n",
        "# Check for missing sample_type values\n",
        "missing_sample_types = merged_df[merged_df['sample_type'].isna()]\n",
        "print(f\"Number of missing sample_type values after comprehensive merge: {len(missing_sample_types)}\")\n",
        "\n",
        "# Drop rows with missing sample_type values for further analysis\n",
        "cleaned_df = merged_df.dropna(subset=['sample_type'])\n",
        "\n",
        "# Display the first few rows of the cleaned merged DataFrame\n",
        "print(\"First few rows of the cleaned merged DataFrame:\")\n",
        "print(cleaned_df.head())\n"
      ],
      "metadata": {
        "id": "KVGLG8cV4YnF"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Separate the data into cancerous and non-cancerous samples based on 'sample_type' column\n",
        "cancerous_samples = cleaned_df[cleaned_df['sample_type'] == 'cancer']\n",
        "non_cancerous_samples = cleaned_df[cleaned_df['sample_type'] == 'normal']\n",
        "\n",
        "# Display the number of samples in each group\n",
        "print(f'Number of cancerous samples: {len(cancerous_samples)}')\n",
        "print(f'Number of non-cancerous samples: {len(non_cancerous_samples)}')\n"
      ],
      "metadata": {
        "id": "mboF77HI4ng_"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "from scipy.stats import ttest_ind\n",
        "\n",
        "# Calculate log2 fold change and p-values\n",
        "def calculate_differential_expression(cancer_df, normal_df, gene_column='gene_id', expression_column='tpm_unstranded'):\n",
        "    diff_exp_results = {\n",
        "        'gene_id': [],\n",
        "        'log2_fold_change': [],\n",
        "        'p_value': []\n",
        "    }\n",
        "\n",
        "    unique_genes = cancer_df[gene_column].unique()\n",
        "\n",
        "    for gene in unique_genes:\n",
        "        cancer_expression = cancer_df[cancer_df[gene_column] == gene][expression_column]\n",
        "        normal_expression = normal_df[normal_df[gene_column] == gene][expression_column]\n",
        "\n",
        "        if len(cancer_expression) > 1 and len(normal_expression) > 1:\n",
        "            log2_fc = np.log2(cancer_expression.mean() / normal_expression.mean())\n",
        "            _, p_value = ttest_ind(cancer_expression, normal_expression)\n",
        "\n",
        "            diff_exp_results['gene_id'].append(gene)\n",
        "            diff_exp_results['log2_fold_change'].append(log2_fc)\n",
        "            diff_exp_results['p_value'].append(p_value)\n",
        "\n",
        "    return pd.DataFrame(diff_exp_results)\n",
        "\n",
        "# Perform differential expression analysis\n",
        "diff_exp_df = calculate_differential_expression(cancerous_samples, non_cancerous_samples)\n",
        "\n",
        "# Display the first few rows of the differential expression results\n",
        "print(\"Differential Expression Results:\")\n",
        "print(diff_exp_df.head())\n"
      ],
      "metadata": {
        "id": "vjRZY8mk41Og"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "IomMHKQe4-Fs"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}