<commit_before>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Class 21:  joint entropy and the REVEAL algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use the bladder cancer gene expression data to test out the REVEAL algorithm. First, we'll load the data and filter to include only genes for which the median log2 expression level is > 12.  That should give us 164 genes to work with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Dimensions of the original matrix of RNA-seq data: 4473\"\n",
      "[2] \"Dimensions of the original matrix of RNA-seq data: 414\" \n",
      "[1] \"Dimensions of the filtered gene matrix: 164\"\n",
      "[2] \"Dimensions of the filtered gene matrix: 414\"\n"
     ]
    }
   ],
   "source": [
    "gene_matrix_for_network <- read.table(\"shared/bladder_cancer_genes_tcga.txt\",\n",
    "                                     sep=\"\\t\",\n",
    "                                     header=TRUE,\n",
    "                                     row.names=1,\n",
    "                                     stringsAsFactors=FALSE)\n",
    "print(sprintf(\"Dimensions of the original matrix of RNA-seq data: %d\", dim(gene_matrix_for_network)))\n",
    "M <- ncol(gene_matrix_for_network)\n",
    "gene_matrix_filt <- gene_matrix_for_network[apply(gene_matrix_for_network, 1, median) > 12, ]\n",
    "\n",
    "N <- nrow(gene_matrix_filt)\n",
    "\n",
    "print(sprintf(\"Dimensions of the filtered gene matrix: %d\", dim(gene_matrix_filt)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Binarize the gene expression matrix using the mean value as a breakpoint, turning it into a NxM matrix of logicals (true/false)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mean_values <- apply(gene_matrix_filt, 1, mean)\n",
    "gene_matrix_binarized <- (gene_matrix_filt > matrix(rep(mean_values, ncol(gene_matrix_filt)), ncol=ncol(gene_matrix_filt)))\n",
    "M <- ncol(gene_matrix_filt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The core part of the REVEAL algorithm is a function that can compute the joint entropy of a collection of binary (TRUE/FALSE) vectors X1, X2, ..., Xn (where length(X1) = length(Xi) = M).\n",
    "Write a function `entropy_multiple_vecs` that takes as its input a nxM matrix (where n is the number of variables, i.e., genes, and M is the number of samples in which gene expression was measured). The function should use the log2 definition of the Shannon entropy. It should return the joint entropy H(X1, X2, ..., Xn) as a scalar numeric value. I have created a skeleton version of this function for you, in which you can fill in the code. I have also created some test code that you can use to test your function, below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "entropy_multiple_vecs <- function(binary_vecs) {\n",
    "    stopifnot(ncol(binary_vecs) == M)  # sanity check input\n",
    "    stopifnot(is.matrix(binary_vecs))  # sanity check input\n",
    "    n <- nrow(binary_vecs)\n",
    "    \n",
    "    ## convert the nxM matrix of logicals to a Mxn data frame of factors\n",
    "    df <- data.frame(t(binary_vecs))\n",
    "    \n",
    "    ## use \"xtabs\" to tabulate the number of occurrences of each combination of logical\n",
    "    ## values for the n inputs; divide by M and call the resulting matrix \"probmat\"\n",
    "    probmat <- xtabs(~., data.frame(lapply(df, factor)))/M\n",
    "\n",
    "    ## define a matrix \"hmat\" that is the negative of probmat element-wise-times log2 of probmat\n",
    "    hmat <- -probmat*log2(probmat)\n",
    "    \n",
    "    ## some entries of hmat will be NaN (do you know why?).  Set them to zero\n",
    "    hmat[is.nan(hmat)] <- 0\n",
    "    \n",
    "    ## return the sum of entries of hmat\n",
    "    sum(hmat)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are some test cases to run, after you have written your function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Test 1 for your entropy function.  H(X,Y,Z) for three independent RVs X,Y,Z should be very close to 3: 2.993245\"\n",
      "[1] \"Test 2 for your entropy function.  H(X,Y,Z) for X,Y,Z always = FALSE, should be identically zero: 0.000000\"\n",
      "[1] \"Test 3 fo ryour entropy function.  H(X1,X2,X3,X4) for the first four rows of your matrix should be 3.937690: 3.937690\"\n"
     ]
    }
   ],
   "source": [
    "print(sprintf(\"Test 1 for your entropy function.  H(X,Y,Z) for three independent RVs X,Y,Z should be very close to 3: %f\",\n",
    "              entropy_multiple_vecs(t(replicate(3,sample(c(FALSE,TRUE), size=M, replace=TRUE))))))\n",
    "\n",
    "print(sprintf(\"Test 2 for your entropy function.  H(X,Y,Z) for X,Y,Z always = FALSE, should be identically zero: %f\",\n",
    "              entropy_multiple_vecs(matrix(rep(FALSE, 3*M), nrow=3))))\n",
    "\n",
    "print(sprintf(\"Test 3 fo ryour entropy function.  H(X1,X2,X3,X4) for the first four rows of your matrix should be 3.937690: %f\",\n",
    "              entropy_multiple_vecs(gene_matrix_binarized[1:4,])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is an example implementation of the REVEAL algorithm using your function. It will find regulators and store regulators in the list `regulators` (the `regulators` list is initialized to NA for each gene which means that no regulators are known for that gene). Study the code and see if you can figure out exactly how it works. Note the use of `do.call`, `expand.grid`, and `apply` in order to vectorize the code, and the use of `list` and `array` in order to define data structures used in the algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "while(length(genes_to_fit) > 0) {\n",
    "    gene_to_fit <- genes_to_fit[1]\n",
    "    if (last_gene_to_fit > gene_to_fit) {\n",
    "        stage <- stage + 1\n",
    "        if (stage > max_stage) {\n",
    "            break\n",
    "        }\n",
    "        print(sprintf(\"Getting joint entropies for combinations of %d genes\", stage+1))\n",
    "        next_stage_entropies <- array(rep(NA,N^(stage+1)), dim=rep(N,stage+1))\n",
    "        entropies_for_stages[[stage+1]] <- array(pbapply(  expand.grid(rep(list(1:N), (stage+1))), 1, function(gene_set) {\n",
    "            entropy_multiple_vecs(gene_matrix_binarized[gene_set,])\n",
    "            }), dim=rep(N, stage+1))\n",
    "    }\n",
    "    print(sprintf(\"trying regulators for gene %d at stage %d  last_gene_to_fit: %f\",\n",
    "                  gene_to_fit, stage, last_gene_to_fit))\n",
    "    possible_regulators <- setdiff(1:N, gene_to_fit)\n",
    "    if (length(possible_regulators) < stage+1) {\n",
    "        break\n",
    "    }\n",
    "    possible_regulator_combs <- matrix(combn(possible_regulators, stage), nrow=stage)\n",
    "    comb_ratios <- apply(possible_regulator_combs, 2, function(reg_comb) {\n",
    "        HG <- entropies_for_stages[[1]][gene_to_fit]\n",
    "        entropy_G_plus_regulators <- do.call(\"[\",c(list(entropies_for_stages[[stage+1]]),c(gene_to_fit, reg_comb)))\n",
    "        entropy_of_regulators <- do.call(\"[\",c(list(entropies_for_stages[[stage]]),reg_comb))\n",
    "        ratio_res <- (HG - entropy_G_plus_regulators + entropy_of_regulators)/HG\n",
    "        ratio_res\n",
    "    })\n",
    "    good_pred_inds <- which(comb_ratios >= ratio_thresh)\n",
    "    if (length(good_pred_inds) > 0) {\n",
    "        regulators_for_gene <- possible_regulator_combs[,which.max(comb_ratios)]\n",
    "        print(sprintf(\"Regulators for gene %d are: \", gene_to_fit))\n",
    "        print(regulators_for_gene)\n",<commit_msg>fixing dependence on pbapply package<commit_after>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Class 21:  joint entropy and the REVEAL algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use the bladder cancer gene expression data to test out the REVEAL algorithm. First, we'll load the data and filter to include only genes for which the median log2 expression level is > 12.  That should give us 164 genes to work with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Dimensions of the original matrix of RNA-seq data: 4473\"\n",
      "[2] \"Dimensions of the original matrix of RNA-seq data: 414\" \n",
      "[1] \"Dimensions of the filtered gene matrix: 164\"\n",
      "[2] \"Dimensions of the filtered gene matrix: 414\"\n"
     ]
    }
   ],
   "source": [
    "gene_matrix_for_network <- read.table(\"shared/bladder_cancer_genes_tcga.txt\",\n",
    "                                     sep=\"\\t\",\n",
    "                                     header=TRUE,\n",
    "                                     row.names=1,\n",
    "                                     stringsAsFactors=FALSE)\n",
    "print(sprintf(\"Dimensions of the original matrix of RNA-seq data: %d\", dim(gene_matrix_for_network)))\n",
    "M <- ncol(gene_matrix_for_network)\n",
    "gene_matrix_filt <- gene_matrix_for_network[apply(gene_matrix_for_network, 1, median) > 12, ]\n",
    "\n",
    "N <- nrow(gene_matrix_filt)\n",
    "\n",
    "print(sprintf(\"Dimensions of the filtered gene matrix: %d\", dim(gene_matrix_filt)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Binarize the gene expression matrix using the mean value as a breakpoint, turning it into a NxM matrix of logicals (true/false)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mean_values <- apply(gene_matrix_filt, 1, mean)\n",
    "gene_matrix_binarized <- (gene_matrix_filt > matrix(rep(mean_values, ncol(gene_matrix_filt)), ncol=ncol(gene_matrix_filt)))\n",
    "M <- ncol(gene_matrix_filt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The core part of the REVEAL algorithm is a function that can compute the joint entropy of a collection of binary (TRUE/FALSE) vectors X1, X2, ..., Xn (where length(X1) = length(Xi) = M).\n",
    "Write a function `entropy_multiple_vecs` that takes as its input a nxM matrix (where n is the number of variables, i.e., genes, and M is the number of samples in which gene expression was measured). The function should use the log2 definition of the Shannon entropy. It should return the joint entropy H(X1, X2, ..., Xn) as a scalar numeric value. I have created a skeleton version of this function for you, in which you can fill in the code. I have also created some test code that you can use to test your function, below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "entropy_multiple_vecs <- function(binary_vecs) {\n",
    "    stopifnot(ncol(binary_vecs) == M)  # sanity check input\n",
    "    stopifnot(is.matrix(binary_vecs))  # sanity check input\n",
    "    n <- nrow(binary_vecs)\n",
    "    \n",
    "    ## convert the nxM matrix of logicals to a Mxn data frame of factors\n",
    "    df <- data.frame(t(binary_vecs))\n",
    "    \n",
    "    ## use \"xtabs\" to tabulate the number of occurrences of each combination of logical\n",
    "    ## values for the n inputs; divide by M and call the resulting matrix \"probmat\"\n",
    "    probmat <- xtabs(~., data.frame(lapply(df, factor)))/M\n",
    "\n",
    "    ## define a matrix \"hmat\" that is the negative of probmat element-wise-times log2 of probmat\n",
    "    hmat <- -probmat*log2(probmat)\n",
    "    \n",
    "    ## some entries of hmat will be NaN (do you know why?).  Set them to zero\n",
    "    hmat[is.nan(hmat)] <- 0\n",
    "    \n",
    "    ## return the sum of entries of hmat\n",
    "    sum(hmat)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are some test cases to run, after you have written your function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Test 1 for your entropy function.  H(X,Y,Z) for three independent RVs X,Y,Z should be very close to 3: 2.993245\"\n",
      "[1] \"Test 2 for your entropy function.  H(X,Y,Z) for X,Y,Z always = FALSE, should be identically zero: 0.000000\"\n",
      "[1] \"Test 3 fo ryour entropy function.  H(X1,X2,X3,X4) for the first four rows of your matrix should be 3.937690: 3.937690\"\n"
     ]
    }
   ],
   "source": [
    "print(sprintf(\"Test 1 for your entropy function.  H(X,Y,Z) for three independent RVs X,Y,Z should be very close to 3: %f\",\n",
    "              entropy_multiple_vecs(t(replicate(3,sample(c(FALSE,TRUE), size=M, replace=TRUE))))))\n",
    "\n",
    "print(sprintf(\"Test 2 for your entropy function.  H(X,Y,Z) for X,Y,Z always = FALSE, should be identically zero: %f\",\n",
    "              entropy_multiple_vecs(matrix(rep(FALSE, 3*M), nrow=3))))\n",
    "\n",
    "print(sprintf(\"Test 3 fo ryour entropy function.  H(X1,X2,X3,X4) for the first four rows of your matrix should be 3.937690: %f\",\n",
    "              entropy_multiple_vecs(gene_matrix_binarized[1:4,])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is an example implementation of the REVEAL algorithm using your function. It will find regulators and store regulators in the list `regulators` (the `regulators` list is initialized to NA for each gene which means that no regulators are known for that gene). Study the code and see if you can figure out exactly how it works. Note the use of `do.call`, `expand.grid`, and `apply` in order to vectorize the code, and the use of `list` and `array` in order to define data structures used in the algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "while(length(genes_to_fit) > 0) {\n",
    "    gene_to_fit <- genes_to_fit[1]\n",
    "    if (last_gene_to_fit > gene_to_fit) {\n",
    "        stage <- stage + 1\n",
    "        if (stage > max_stage) {\n",
    "            break\n",
    "        }\n",
    "        print(sprintf(\"Getting joint entropies for combinations of %d genes\", stage+1))\n",
    "        entropies_for_stages[[stage+1]] <- array(apply(  expand.grid(rep(list(1:N), (stage+1))), 1, function(gene_set) {\n",
    "            entropy_multiple_vecs(gene_matrix_binarized[gene_set,])\n",
    "            }), dim=rep(N, stage+1))\n",
    "    }\n",
    "    print(sprintf(\"trying regulators for gene %d at stage %d  last_gene_to_fit: %f\",\n",
    "                  gene_to_fit, stage, last_gene_to_fit))\n",
    "    possible_regulators <- setdiff(1:N, gene_to_fit)\n",
    "    if (length(possible_regulators) < stage+1) {\n",
    "        break\n",
    "    }\n",
    "    possible_regulator_combs <- matrix(combn(possible_regulators, stage), nrow=stage)\n",
    "    comb_ratios <- apply(possible_regulator_combs, 2, function(reg_comb) {\n",
    "        HG <- entropies_for_stages[[1]][gene_to_fit]\n",
    "        entropy_G_plus_regulators <- do.call(\"[\",c(list(entropies_for_stages[[stage+1]]),c(gene_to_fit, reg_comb)))\n",
    "        entropy_of_regulators <- do.call(\"[\",c(list(entropies_for_stages[[stage]]),reg_comb))\n",
    "        ratio_res <- (HG - entropy_G_plus_regulators + entropy_of_regulators)/HG\n",
    "        ratio_res\n",
    "    })\n",
    "    good_pred_inds <- which(comb_ratios >= ratio_thresh)\n",
    "    if (length(good_pred_inds) > 0) {\n",
    "        regulators_for_gene <- possible_regulator_combs[,which.max(comb_ratios)]\n",
    "        print(sprintf(\"Regulators for gene %d are: \", gene_to_fit))\n",
    "        print(regulators_for_gene)\n",
{'fraction_non_alphanumeric': 0.15646661532081285, 'fraction_numerical': 0.013981870963919663, 'mean_word_length': 2.6783612987578995, 'pattern_counts': {'":': 110, '<': 56, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5410907', 'n_tokens_mistral': 6127, 'n_tokens_neox': 5772, 'n_words': 1578}