04_GWAS_Case_Control_Association_Testing_practice.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyOvkZfYEWU3kHJpVhq8Ombs",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "ir",
      "display_name": "R"
    },
    "language_info": {
      "name": "R"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/luuloi/GWAS_Introduction_2023/blob/main/04_GWAS_Case_Control_Association_Testing_practice.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### **Task 1: Data Inspection and Basic Descriptive Statistics**\n",
        "\n",
        "1. After loading your data:\n",
        "    - Check the last few rows of the `LHON` dataframe.\n",
        "    - Check the unique values in the `GENO` and `PHENO` columns.\n",
        "    - Calculate the mean, median, and mode for any numerical column present in the `LHON` dataframe.\n",
        "    \n",
        "### **Task 2: Data Visualization**\n",
        "\n",
        "1. Create a pie chart for the `PHENO` column showcasing the distribution of cases and controls.\n",
        "2. Plot the distribution (histogram) of any numeric columns present in the dataframe.\n",
        "3. Create a scatter plot between two numeric columns (if available) and use color to differentiate between the two `PHENO` groups.\n",
        "\n",
        "### **Task 3: Statistical Analysis**\n",
        "\n",
        "1. Choose any two categorical columns from the `LHON` dataframe (other than `GENO` and `PHENO`) and create a contingency table.\n",
        "2. Perform the Chi-square test on this new table and interpret the results.\n",
        "\n",
        "### **Task 4: Logistic Regression and Interpretation**\n",
        "\n",
        "1. If you have more categorical predictors, encode them similarly to `genoadd` (using numerical encoding) and add them to your logistic regression model.\n",
        "2. Interpret the significance of each predictor using the summary results.\n",
        "3. Predict the probabilities for the training set using your model and classify them with a different cutoff (say 0.7).\n",
        "4. Construct a confusion matrix using this new cutoff.\n",
        "\n",
        "### **Task 5: ROC and Precision-Recall Curve Analysis**\n",
        "\n",
        "1. Calculate the area under the precision-recall curve.\n",
        "2. Compare the AUC from the ROC and the precision-recall curve. What can you conclude from the differences or similarities in the AUC values?\n",
        "3. Change the logistic regression model (e.g., add or remove predictors) and observe how it affects the ROC curve and precision-recall curve."
      ],
      "metadata": {
        "id": "hNv8Bv8nucCI"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "## Read data\n",
        "LHON = read.table(\"https://raw.githubusercontent.com/luuloi/GWAS_Introduction_2023/main/data/LHON.txt\", header=TRUE, stringsAsFactors=TRUE)"
      ],
      "metadata": {
        "id": "tslxYxhJu5hv"
      },
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "## Check first few rows\n",
        "head(LHON)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 286
        },
        "id": "a4_xRc7MvAjj",
        "outputId": "a4fba167-fb47-47a8-ba6f-52b7c389425f"
      },
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "<table class=\"dataframe\">\n",
              "<caption>A data.frame: 6 × 3</caption>\n",
              "<thead>\n",
              "\t<tr><th></th><th scope=col>IID</th><th scope=col>GENO</th><th scope=col>PHENO</th></tr>\n",
              "\t<tr><th></th><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;fct&gt;</th></tr>\n",
              "</thead>\n",
              "<tbody>\n",
              "\t<tr><th scope=row>1</th><td>ID1</td><td>TT</td><td>CONTROL</td></tr>\n",
              "\t<tr><th scope=row>2</th><td>ID2</td><td>CT</td><td>CONTROL</td></tr>\n",
              "\t<tr><th scope=row>3</th><td>ID3</td><td>TT</td><td>CASE   </td></tr>\n",
              "\t<tr><th scope=row>4</th><td>ID4</td><td>CT</td><td>CONTROL</td></tr>\n",
              "\t<tr><th scope=row>5</th><td>ID5</td><td>TT</td><td>CONTROL</td></tr>\n",
              "\t<tr><th scope=row>6</th><td>ID6</td><td>TT</td><td>CONTROL</td></tr>\n",
              "</tbody>\n",
              "</table>\n"
            ],
            "text/markdown": "\nA data.frame: 6 × 3\n\n| <!--/--> | IID &lt;fct&gt; | GENO &lt;fct&gt; | PHENO &lt;fct&gt; |\n|---|---|---|---|\n| 1 | ID1 | TT | CONTROL |\n| 2 | ID2 | CT | CONTROL |\n| 3 | ID3 | TT | CASE    |\n| 4 | ID4 | CT | CONTROL |\n| 5 | ID5 | TT | CONTROL |\n| 6 | ID6 | TT | CONTROL |\n\n",
            "text/latex": "A data.frame: 6 × 3\n\\begin{tabular}{r|lll}\n  & IID & GENO & PHENO\\\\\n  & <fct> & <fct> & <fct>\\\\\n\\hline\n\t1 & ID1 & TT & CONTROL\\\\\n\t2 & ID2 & CT & CONTROL\\\\\n\t3 & ID3 & TT & CASE   \\\\\n\t4 & ID4 & CT & CONTROL\\\\\n\t5 & ID5 & TT & CONTROL\\\\\n\t6 & ID6 & TT & CONTROL\\\\\n\\end{tabular}\n",
            "text/plain": [
              "  IID GENO PHENO  \n",
              "1 ID1 TT   CONTROL\n",
              "2 ID2 CT   CONTROL\n",
              "3 ID3 TT   CASE   \n",
              "4 ID4 CT   CONTROL\n",
              "5 ID5 TT   CONTROL\n",
              "6 ID6 TT   CONTROL"
            ]
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# First create a genotype variable with an additive coding based on the counts of the number of T alleles\n",
        "LHON$genoadd <- with(LHON, 0 + 1*(GENO==\"CT\") + 2*(GENO==\"TT\"))\n",
        "head(LHON)"
      ],
      "metadata": {
        "id": "iGKbdoMmvqEC"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 1.1 Check the last few rows of the LHON dataframe"
      ],
      "metadata": {
        "id": "NGz_AZLDvr2f"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "tail()"
      ],
      "metadata": {
        "id": "PpjUGRvpvz6g"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 1.2 Check the unique values in the GENO and PHENO columns"
      ],
      "metadata": {
        "id": "CF9z5pZsvvZi"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "unique(LHON$)\n",
        "unique(LHON$)"
      ],
      "metadata": {
        "id": "rKdBfqE7v2m2"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 1.3 Calculate the mean, median, and mode"
      ],
      "metadata": {
        "id": "FZc5j7ghv9VQ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "mean(LHON$genoadd)\n",
        "median(LHON$genoadd)\n",
        "\n",
        "get_mode <- function(v) {\n",
        "   uniqv <- unique(v)\n",
        "   uniqv[which.max(tabulate(match(v, uniqv)))]\n",
        "}\n",
        "mode_genoadd <- get_mode(LHON$)"
      ],
      "metadata": {
        "id": "ynPKRPYTurHn"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 2.1 Pie chart for the PHENO column"
      ],
      "metadata": {
        "id": "lcfgLOhTwW-c"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "pie(table(LHON$), main=\"Distribution of PHENO\", col=c(\"lightblue\", \"lightcoral\"))"
      ],
      "metadata": {
        "id": "H3dMJDEhwRH1"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 2.2 Histogram for numeric column"
      ],
      "metadata": {
        "id": "zucATOIGwb8A"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "hist(LHON$, main=\"Distribution of genoadd\", col=\"lightgreen\")"
      ],
      "metadata": {
        "id": "equLPDUqwbJX"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 2.3 Scatter plot between two numeric columns (assuming we have another numeric column 'numeric_column_2')"
      ],
      "metadata": {
        "id": "_ttCmSB9wigG"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# For this demonstration, let's use genoadd vs. genoadd (just to showcase the command)\n",
        "library(ggplot2)\n",
        "ggplot(LHON, aes(x=genoadd, y=genoadd, color=)) + geom_point() + theme_minimal()"
      ],
      "metadata": {
        "id": "l2XoKBz2whbl"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 3.1 Contingency table"
      ],
      "metadata": {
        "id": "N5cjHVY9xMBK"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Assuming we use GENO and newpheno for the contingency table\n",
        "cont_table_new <- table(LHON$, LHON$)"
      ],
      "metadata": {
        "id": "ckqxbPn5xJ96"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 3.2 Chi-square test"
      ],
      "metadata": {
        "id": "YH-1IE4pxOuz"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "chi_test_new <- chisq.test(cont_table_new)\n",
        "print(chi_test_new)"
      ],
      "metadata": {
        "id": "AcwTBqY-xP5y"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 4.1 Adding genoadd to the model"
      ],
      "metadata": {
        "id": "MRh4IVX9xbYt"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "logistmod4 <- glm(newpheno ~ genoadd + , family=binomial(link=\"logit\"), data=LHON)"
      ],
      "metadata": {
        "id": "g_NTtUOAxWhL"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 4.2 Interpret the summary"
      ],
      "metadata": {
        "id": "X88yVs-mxin4"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "summary()"
      ],
      "metadata": {
        "id": "yXCcMVyXxgBX"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 4.3 Predict using the model"
      ],
      "metadata": {
        "id": "a7x2KaJAxpZZ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "LHON$pred_new <- predict(, type=\"response\")\n",
        "\n",
        "# Classification with 0.7 as the cutoff\n",
        "LHON$predicted_pheno <- ifelse(LHON$pred_new > 0.7, \"CASE\", \"CONTROL\")"
      ],
      "metadata": {
        "id": "hJTlLqlexoRL"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 4.4 Confusion matrix with new cutoff"
      ],
      "metadata": {
        "id": "7SaujGRjxywr"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "confusionMatrix(table(LHON$predicted_pheno, LHON$))"
      ],
      "metadata": {
        "id": "LB4uloB0xxze"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 5.1 Area under the precision-recall curve"
      ],
      "metadata": {
        "id": "EIqQtw6n0lG2"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Install precrec package\n",
        "pacman::p_load(\"precrec\")\n",
        "\n",
        "# loading library\n",
        "library(precrec)\n",
        "\n",
        "ROCnPR_new <- evalmod(scores = LHON$pred_new, labels = LHON$)\n",
        "auc_new <- auc(ROCnPR_new)\n",
        "auc_new[which(auc_new$type == \"PRC\"),]"
      ],
      "metadata": {
        "id": "7VdjrXR90VUn"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 5.2 Compare AUC from ROC and precision-recall curve"
      ],
      "metadata": {
        "id": "Bcro9U5P00Gi"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "print(auc_new)"
      ],
      "metadata": {
        "id": "YXBhEnjh0YRq"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 5.3 Change logistic regression model (here, we remove GENO)"
      ],
      "metadata": {
        "id": "OstqKVCp0yMC"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "logistmod5 <- glm(newpheno ~ genoadd, family=binomial(link=\"logit\"), data=LHON)\n",
        "LHON$pred5 <- predict(logistmod5, type=\"response\")\n",
        "ROCnPR5 <- evalmod(scores = LHON$pred5, labels = LHON$)\n",
        "autoplot(ROCnPR5, \"ROC\")\n",
        "autoplot(ROCnPR5, \"PRC\")"
      ],
      "metadata": {
        "id": "oQpu0Fmh0ZfE"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}