interaction_detection.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# importing all the required packages\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy.stats\n",
    "import math\n",
    "from time import process_time_ns\n",
    "import statsmodels.stats.multitest as ssm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the datafile with SNPs and Phenotype\n",
    "data = pd.read_csv(\"sample10.csv\")\n",
    "phenotype = data.iloc[:,-1] # extracting the phenotype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# defining and initializing required variables\n",
    "\n",
    "# defining the significance level \n",
    "alpha = .05\n",
    "\n",
    "# storing the count of number of samples in the variable n\n",
    "n = data.shape[0]\n",
    "\n",
    "# defining degrees of freedom where n is number of samples\n",
    "df = n-4\n",
    "\n",
    "# defining the critical value\n",
    "critical_value = scipy.stats.t.ppf(q=1-alpha/2,df=df)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is the algoritm implementation for the partial correlation of interaction of genes for Epistasis. We use two different ways of encoding the interaction terms. Firstly, the traditional method of cartesian product and then by using the XOR penetrance function."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Algorithm 1- Algorithm for 2 way epistasis detection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Algorithm for Interaction Coefficient for Pairwise Epistasis\n",
    "\"\"\"Parameters for the method:\n",
    "    snp1 - First snp in the pair to be checked for interaction\n",
    "    snp2 - Second snp in the pair to be checked for interaction\n",
    "    phenotype - Phenotype vector\n",
    "\"\"\"\n",
    "def epistatis(snp1, snp2, phenotype): \n",
    "    try:\n",
    "        # removing the intercept by mean centering\n",
    "        snp1_tilde = snp1 - snp1.mean()\n",
    "        snp2_tilde = snp2 - snp2.mean()\n",
    "        phenotype_tilde = phenotype - phenotype.mean()\n",
    "\n",
    "        # declaring the interaction vector - will contain the interaction term of snp1 and snp2\n",
    "        interaction_vector = pd.Series(dtype='float64')\n",
    "    \n",
    "        # defining the interaction vector - either using cartesian product or XOR\n",
    "        interaction_vector = snp1.mul(snp2) # using cartesian product\n",
    "        # interaction_vector = (snp1%2 + snp2%2)%2 # using XOR penetrance\n",
    "    \n",
    "        # print(interaction_vector)\n",
    "        # mean centering the interaction vector \n",
    "        interaction_vector_tilde = interaction_vector - interaction_vector.mean()\n",
    "   \n",
    "        # computing the dot products as explained in the algorithm\n",
    "        x = (snp1_tilde.dot(interaction_vector_tilde)/snp1_tilde.dot(snp1_tilde)) * (snp1_tilde) # 2nd term in the v variable, using to breakdown the formula\n",
    "        v = interaction_vector_tilde - x\n",
    "        q2 = snp2_tilde - (((snp1_tilde.dot(snp2_tilde)) / (snp1_tilde.dot(snp1_tilde))) * snp1_tilde)\n",
    "        v = v - ((interaction_vector_tilde.dot(q2)/q2.dot(q2))*q2)\n",
    "        b3 = (v.dot(phenotype_tilde)) / (v.dot(v)) # interactionn coefficient, referred as beta_3 in the paper\n",
    "\n",
    "\n",
    "        # residual calculation\n",
    "        residual = phenotype_tilde - (snp1_tilde.dot(phenotype_tilde)/snp1_tilde.dot(snp1_tilde))*snp1_tilde\n",
    "        residual = residual - (phenotype_tilde.dot(q2)/q2.dot(q2))*q2\n",
    "        residual = residual - b3 * v\n",
    "\n",
    "        # print(b3)\n",
    "\n",
    "        v = pd.Series(np.squeeze(np.asarray(v)), dtype='float64')\n",
    "        residual = pd.Series(np.squeeze(np.asarray(residual)), dtype='float64') \n",
    "\n",
    "\n",
    "        t_test = np.sqrt(snp1.shape[0]-4)*np.sqrt(v.dot(v))*b3/(np.sqrt(residual.dot(residual)))\n",
    "\n",
    "        # t_test = np.sqrt(df)*b3/(np.sqrt(1-b3*b3))\n",
    "        p_val = scipy.stats.t.sf(abs(t_test), df) * 2\n",
    "    except Exception as e:\n",
    "        print(\"Error pair detected with error: \", e)\n",
    "        b3 = 0\n",
    "        t_test = 0\n",
    "        p_val = 1\n",
    "\n",
    "    return b3, t_test, p_val # returning the interaction coefficient,, t test value and p value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# applying the above defined method to the dataset\n",
    "\n",
    "p_value_locus = [] # list of tuples containing the p_value and the two interacting loci\n",
    "for i in range(0, data.shape[1]-1):\n",
    "    for j in range(i+1, data.shape[1]-1):\n",
    "        interacting_snp_1 = data.iloc[:,i]\n",
    "        interacting_snp_2 = data.iloc[:,j]\n",
    "        b3, t_test, p_val = epistatis(interacting_snp_1, interacting_snp_2, phenotype)\n",
    "        p_value_locus.append((p_val, data.columns[i], data.columns[j]))\n",
    "\n",
    "        print(\" Interacting SNP 1 = {0} \\t  Interacting SNP 2 ={1} \\t  Beta Coefficient(r) = {2}  \\t t_test = {3} \\t p_val = {4} \".format(data.columns[i], data.columns[j], b3, t_test, p_val))\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# performing fdr correction to correct for multiple tests\n",
    "\n",
    "all_p_values = [] # list of all the p_values of the pairwise combinations\n",
    "for i in range(0, len(p_value_locus)):\n",
    "    all_p_values.append(p_value_locus[i][0])\n",
    "\n",
    "accept, corrected_pvals = ssm.fdrcorrection(np.asarray(all_p_values).flatten(), alpha=0.001)\n",
    "\n",
    "for i in range(0, len(corrected_pvals)):\n",
    "    # printing the significant pairs after p value correction\n",
    "    if corrected_pvals[i] < 0.001:\n",
    "        print(\" p-value = {0} \\t adjusted p-value = {1} \\t loci 1 = {2} \\t loci 2 = {3} \".format( p_value_locus[i][0],  corrected_pvals[i], p_value_locus[i][1], p_value_locus[i][2]))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is the algorithm implementation for 3 way epistasis detection, it can be generalized to k-wise epistasis.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Algorithm 2- Algorithm implementation for 3 way epistasis detection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Algorithm for Interaction Coefficient for 3-way Epistasis\n",
    "\"\"\"Parameters for the method:\n",
    "    i - Index of first snp in the triplet to be checked for interaction\n",
    "    j - Index of second snp in the triplet to be checked for interaction\n",
    "    k - Index of third snp in the triplet to be checked for interaction\n",
    "    phenotype - Phenotype vector\n",
    "    interaction - Type of interaction to be used for the 3-way epistasis; either 'cartesian' or 'xor'\n",
    "\"\"\"\n",
    "\n",
    "def epistasis_three_way(i, j, k, phenotype, interaction = 'cartesian'):\n",
    "    try:\n",
    "        snp = np.matrix(data.iloc[:,:-1]) # reading the SNP data\n",
    "        \n",
    "        phenotype_tilde = phenotype - phenotype.mean() # reading the phenotype data\n",
    "        phenotype_tilde = np.matrix(phenotype_tilde)\n",
    "\n",
    "        x = np.transpose(np.matrix(snp[:,i])) # interacting SNP1\n",
    "        y = np.transpose(np.matrix(snp[:,j])) # interacting SNP2\n",
    "        z = np.transpose(np.matrix(snp[:,k])) # interacting SNP3\n",
    "        q = np.concatenate((x,y,z), axis=0)\n",
    "        # Calculation of the variable q(Q)\n",
    "\n",
    "        if interaction == 'cartesian': \n",
    "          a = np.multiply(x,y)\n",
    "\n",
    "          q = np.concatenate((q,a),axis=0)\n",
    "\n",
    "          a = np.multiply(y,z)\n",
    "\n",
    "          q = np.concatenate((q,a),axis=0)\n",
    "          a = np.multiply(x,z)\n",
    "\n",
    "          q = np.concatenate((q,a),axis=0)\n",
    "          a = np.multiply(x,y)\n",
    "          a = np.multiply(a,z)\n",
    "          q = np.concatenate((q,a),axis=0)\n",
    "        elif interaction == 'xor':\n",
    "          a = (x%2 + y%2)%2\n",
    "          q = np.concatenate((q,a),axis=0)\n",
    "          a = (y%2 + z%2)%2\n",
    "          q = np.concatenate((q,a),axis=0)\n",
    "          a = (x%2 + z%2)%2\n",
    "          q = np.concatenate((q,a),axis=0)\n",
    "          a = (x%2 + y%2)%2\n",
    "          a = (a%2 + z%2)%2\n",
    "          q = np.concatenate((q,a),axis=0)\n",
    "        q = q - q.mean(axis=1)\n",
    "        q = np.concatenate((q,phenotype_tilde),axis=0)\n",
    "        # beta for interaction of 3way epistasis after q is calculated\n",
    "        \n",
    "        q = q.transpose()\n",
    "        beta = 0\n",
    "        for l in range(1,8):\n",
    "          gs_curr = q[:,l]\n",
    "          for m in range(l):\n",
    "            gs_prev = q[:,m]\n",
    "            curr = np.matmul(np.transpose(gs_curr), gs_prev)/np.matmul(np.transpose(gs_prev), gs_prev)\n",
    "            if l == 7 and m == 6:\n",
    "              beta = curr\n",
    "            q[:,l] = q[:,l] - curr[0,0]*gs_prev\n",
    "        v = q[:,6]\n",
    "        resid = q[:,7]\n",
    "        dot_v_v = np.matmul(np.transpose(v),v)[0,0]\n",
    "        dot_resid_resid = np.matmul(np.transpose(resid),resid)[0,0]\n",
    "        t_test = np.sqrt(df)*beta*np.sqrt(dot_v_v)/np.sqrt(dot_resid_resid)\n",
    "        t_test = t_test[0,0]\n",
    "        p_val = 2*scipy.stats.t.sf(abs(t_test), df)\n",
    "\n",
    "    except Exception as e:\n",
    "        print(\"Error pair detected with error: \", e)\n",
    "        t_test = 0\n",
    "        p_val = 1\n",
    "\n",
    "    return t_test, p_val\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# applying the above defined method (3 way Epistasis) to the dataset\n",
    "\n",
    "p_value_locus = [] # list of tuples containing the p_value and the two interacting loci\n",
    "for i in range(0, data.shape[1]-1):\n",
    "    print(i)\n",
    "    for j in range(i+1, data.shape[1]-1):\n",
    "        for k in range(j+1, data.shape[1]-1):\n",
    "            t_test, p_val = epistasis_three_way(i, j, k, phenotype, interaction = 'cartesian') # specifying the interaction model\n",
    "            \n",
    "            p_value_locus.append((p_val, data.columns[i], data.columns[j], data.columns[k]))\n",
    "\n",
    "            print(\" Interacting SNP 1 = {0} \\t  Interacting SNP 2 ={1} \\t Interacting SNP 3 = {2}  \\t t_test = {3} \\t p_val = {4} \".format(data.columns[i], data.columns[j], data.columns[k], t_test, p_val))\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a csv file for the p-values of the three-way interactions\n",
    "pd.DataFrame(p_value_locus).to_csv('p_value_locus_cartesian.csv')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is the implementation of the permutation testing algorithm, which helps to validate the detected epistasis pairs.\n",
    "Algorithm 3- Permutation Testing Algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# defining the significance level \n",
    "alpha = .001\n",
    "\n",
    "# storing the count of number of samples in the variable n\n",
    "n = data.shape[0]\n",
    "\n",
    "# number of snps \n",
    "m = data.shape[1]\n",
    "\n",
    "# defining degrees of freedom where n is number of samples\n",
    "df = n-4\n",
    "\n",
    "# defining the critical value\n",
    "critical_value = scipy.stats.t.ppf(q=1-alpha/2,df=df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "critical_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# making the permuted phenotype matrix. Each column will have a version of the permutated phenotype. \n",
    "# Assuming your phenotype column is named 'y' in your DataFrame 'data\n",
    "\n",
    "phenotype_column = data['bmi_w_tail_res']\n",
    "K = 10 # number of times to be permutated\n",
    "\n",
    "# Create an empty matrix P to store permuted phenotype matrices\n",
    "P = np.zeros((len(phenotype_column), K))\n",
    "\n",
    "# Iterate over each permutation\n",
    "for i in range(K):\n",
    "    # Shuffle the values in the phenotype column\n",
    "    permuted_phenotype = np.random.permutation(phenotype_column)\n",
    "    \n",
    "    # Store the permuted phenotype in the corresponding column of P\n",
    "    P[:, i] = permuted_phenotype\n",
    "P_df = pd.DataFrame(P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# changing data to numpy array from dataframe for easier mathematical manipulations\n",
    "data = data.to_numpy()\n",
    "# centering the data\n",
    "data_tilde = data - data.mean()\n",
    "\n",
    "tilde_exp = data_tilde # defining variable tilde_exp to hold the centered data\n",
    "cols = data.shape[1] # number of columns in the data, i.e number of SNPs in the data\n",
    "\n",
    "tilde_p = phenotype - phenotype.mean() # centering the phenotype\n",
    "\n",
    "permuted_p = P # permutated phenotype matrix\n",
    "num_tests = K # number of permutation tests\n",
    "\n",
    "detected = 0 # to ocunt the number of pairs detetcted\n",
    "results = pd.DataFrame([[0,0,0,0]], columns=['pval','ttest','snp1','snp2'])\n",
    "for i in range(cols):\n",
    "    for j in range(i+1,cols):\n",
    "        print(i,j)\n",
    "        try:\n",
    "           \n",
    "           # Encoding interaction for epistasis\n",
    "           # traditional way to just multiple vectors for loci\n",
    "           #tilde_i = np.multiply(snp[:,i], snp[:,j])\n",
    "           # using xor penetrance function\n",
    "           #tilde_i = (snp[:,i]%2+snp[:,j]%2)%2\n",
    "           tilde_i = (data[:,i]%2+data[:,j]%2)%2\n",
    "           tilde_i = tilde_i - np.mean(tilde_i, axis=0)\n",
    "           print(tilde_i)\n",
    "\n",
    "\n",
    "           # check_epistasis(g_i,g_j, I_ij) - all code below would be this\n",
    "           dot_gi_ti = np.dot(tilde_i, tilde_exp[:,i])\n",
    "           dot_gi_p = np.dot(tilde_p, tilde_exp[:,i])\n",
    "           dot_gi_gi = np.dot(tilde_exp[:,i], tilde_exp[:,i])\n",
    "\n",
    "           v = tilde_i - dot_gi_ti/dot_gi_gi*tilde_exp[:,i]\n",
    "           resid = tilde_p - dot_gi_p/dot_gi_gi*tilde_exp[:,i]\n",
    "\n",
    "           dot_gj_ti = np.dot(tilde_i, tilde_exp[:,j])\n",
    "           dot_gj_gj = np.dot(tilde_exp[:,j], tilde_exp[:,j])\n",
    "           dot_gi_gj = np.dot(tilde_exp[:,i], tilde_exp[:,j])\n",
    "           q_2 = tilde_exp[:,j]-dot_gi_gj/dot_gi_gi*tilde_exp[:,i]\n",
    "\n",
    "           dot_q2_ti = np.dot(tilde_i, q_2)\n",
    "           dot_q2_p = np.dot(tilde_p, q_2)\n",
    "           dot_q2_q2 = np.dot(q_2, q_2)\n",
    "           # residual of I\n",
    "           v = v-dot_q2_ti/dot_q2_q2*q_2 # gives the residual I from the paper\n",
    "           resid = resid - dot_q2_p/dot_q2_q2*q_2 #  residual calculated as per algorithm in the paper\n",
    "\n",
    "           # dot product of I with it's residual\n",
    "           dot_v_v = np.dot(v,v) \n",
    "           dot_v_p = np.dot(v,tilde_p)\n",
    "\n",
    "           #direct-test for beta_3 called r below\n",
    "           r = dot_v_p/dot_v_v\n",
    "           resid = resid - r*v\n",
    "           dot_resid_resid = np.dot(resid,resid)\n",
    "           \n",
    "           t_test = np.sqrt(df)*r*np.sqrt(dot_v_v)/np.sqrt(dot_resid_resid)\n",
    "\n",
    "\n",
    "           #permutation testing for pvals now\n",
    "           # first loci dot products\n",
    "           perm_gi = np.dot(permuted_p.T,tilde_exp[:,i].reshape(n,1))/dot_gi_gi \n",
    "           resid_p = permuted_p - np.outer(tilde_exp[:,i].reshape(n,1),perm_gi) \n",
    "           perm_gj = np.dot(permuted_p.T,q_2)/dot_q2_q2\n",
    "           resid_p = resid_p - np.outer(q_2,perm_gj) \n",
    "           perm_i = np.dot(permuted_p.T,v).reshape(num_tests,1) \n",
    "           print(\"perm beta3s\")\n",
    "           print(\"perm beta3s: \\n\", perm_i/dot_v_v) \n",
    "           resid_p = resid_p - np.outer(tilde_i,perm_i/dot_v_v) \n",
    "           test_betas = np.sum(np.multiply(permuted_p, resid_p),axis=0).reshape(num_tests,1) \n",
    "           test_residp= np.sum(np.multiply(permuted_p, permuted_p),axis=0).reshape(num_tests,1) \n",
    "           test_betas = np.sqrt(dot_v_v*test_betas)\n",
    "           test_betas = np.sqrt(df)*np.divide(perm_i,test_betas)\n",
    "           p_val = sum(sum(abs(test_betas) >= abs(t_test)))/num_tests\n",
    "        except Exception as e:\n",
    "           print(\"Error pair (\",i, \",\",j,\"):\",e)\n",
    "           t_test = 0\n",
    "           p_val = 1\n",
    "        df2 = pd.DataFrame([[p_val,t_test,i,j]], columns=['pval','ttest','snp1','snp2'])\n",
    "        results = pd.concat([results,df2]) \n",
    "        if abs(t_test) > critical_value and p_val < alpha:\n",
    "            detected = detected + 1\n",
    "        print(results)\n",
    "t2 = process_time_ns()\n",
    "results.to_csv(\"perm12_test.csv\")\n",
    "print(\"number of detected pairs:\", detected)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "506c12871c44be1527a6e7c25c95d2dc736aaab0b3457406043d1256ac2febd5"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}