1-algorithms-that-learn-R.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "source": [
    "Algorithms That Learn To Predict \n",
    "================================\n",
    "\n",
    "Machine learning (ML) refers to the use of algorithms which can learn\n",
    "from data. The inputs to a ML algorithm will generally be some sort of\n",
    "data set---which I will refer to as the *training data*---and the\n",
    "output will usually be another algorithm, which I will call a\n",
    "*fit model*, or sometimes just a *fit* if I'm feeling lazy.\n",
    "-   ML algorithm : training data $\\rightarrow$ fit model\n",
    "\n",
    "The fit model itself also takes data as input, and generally requires\n",
    "that the data provided to it must be very similar in nature to that\n",
    "provided to the ML algorithm as training data: For example, assuming\n",
    "the data sets in question are represented in table form, the data\n",
    "provided to the fit model must usually have all or almost all of the\n",
    "same columns as the training data set did. However, the output from\n",
    "the fit model is usually much simpler, often consisting of a predicted\n",
    "*numeric value* or *categorical label* for each individual\n",
    "sampling unit of the data set.\n",
    "-   Fit model : test data $\\rightarrow$ predicted values\n",
    "\n",
    "It is often convenient to package an ML algorithm up into a\n",
    "`function`; this is the most common strategy employed by ML\n",
    "algorithms implemented in R. A classic example is the `lm` function\n",
    "in R:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "outputs": [],
   "source": [
    "set.seed(123)\n",
    "n = 20\n",
    " ## generate some random data for two variables:\n",
    "predictor1 = rnorm(n)\n",
    "predictor2 = rnorm(n)\n",
    " ## now set a third variable equal to a weighted sum of those\n",
    " ## two variables plus a random error term:\n",
    "output = 2*predictor1 + predictor2 + rnorm(n)\n",
    " ## bundle up the three variables composing our data set into a\n",
    " ## data.frame object:\n",
    "allData = data.frame(p1=predictor1, p2=predictor2, out=output)\n",
    " ## split allData into training and test sets:\n",
    "trainData = allData[1:10, ]\n",
    "testData = allData[11:20, ]  ## should not overlap trainData!\n",
    " ## now train model using only trainData:\n",
    "fitModel = lm(out ~ p1 + p2, data=trainData)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "source": [
    "Now we can use the `fitModel` to make predictions on rows 11--20 of\n",
    "`allData`; in R, this is usually, though not universally, done by\n",
    "calling the `predict` function with\n",
    "1.  the object representing the fit model as first argument and\n",
    "2.  the test data set---usually in the form of a `data.frame`---as\n",
    "    the second argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "outputs": [],
   "source": [
    " ## generate predictions for test data:\n",
    "predictions = predict(fitModel, testData)\n",
    "library(ggplot2); theme_set(theme_bw())\n",
    " ## plot actual values of out column against predicted values\n",
    " ## for the test data using ggplot2::qplot\n",
    "qplot(predictions, testData$out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "source": [
    "This an example of *supervised learning*, in which one of the\n",
    "variables in the training data set (`out` in this case) is treated\n",
    "as an output to be predicted using the others. The test set does not\n",
    "actually need to have this variable present to make predictions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "outputs": [],
   "source": [
    "predictions2 = predict(fitModel, testData[ , c(\"p1\", \"p2\")])\n",
    " ## these predictions should be the same as before:\n",
    "all(predictions2 == predictions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "source": [
    "Thus, in supervised learning approaches the fit model requires only a\n",
    "subset of the variables present in the training data to be present in\n",
    "the test data in order to make predictions.\n",
    "\n",
    "In *unsupervised learning* this is not the case, and we must\n",
    "generally have all variables from the training data also present in\n",
    "any test data that we wish to make predictions on. What is this\n",
    "\"unsupervised learning\", you ask, and what might it be used to\n",
    "predict? Let's consider an example to make things more concrete:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "outputs": [],
   "source": [
    " ## need clue library to make predictions from kmeans fit\n",
    " ## install.packages(\"clue\")  ## uncomment and run if necessary\n",
    "library(clue)\n",
    " ## use k-means clustering algorithm to fit 2 clusters to training data\n",
    "kmeansFit = kmeans(trainData, 2)\n",
    " ## inspect kmeansFit object in R terminal:\n",
    "kmeansFit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "outputs": [],
   "source": [
    " ## use clue::cl_predict instead of ordinary predict with kmeans\n",
    " ## to predict which cluster each test datum is in\n",
    "kmPredictions = cl_predict(kmeansFit, testData)\n",
    "kmPredictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "outputs": [],
   "source": [
    " ## two clusters in this case correspond to low and high values of \"out\":\n",
    "qplot(factor(kmPredictions), testData$out, geom=\"boxplot\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "source": [
    "As seen in this example, unsupervised learning algorithms try to find\n",
    "some latent structure in the training data---such as the carving of\n",
    "the variable space (frequently called *feature space* in ML) into\n",
    "two disjoint clusters done by `kmeans`, about which more will be\n",
    "said later.\n",
    "Many unsupervised learning algorithms, including `kmeans`, produce\n",
    "fit models which can be used to determine how test data would fit into\n",
    "the learned latent structure; for instance, here we were able to\n",
    "assign each test datum to one of the two clusters learned from the\n",
    "training data set. There are some unsupervised learning approaches\n",
    "which generate fit models which are not immediately equipped to make\n",
    "test set predictions, however---hierarchical clustering and tSNE come\n",
    "to mind here---which can limit their utility in some situations.\n\n",
    "Data \n",
    "====\n",
    "\n",
    "Machine learning---perhaps I should lose the qualifier and just say\n",
    "learning---isn't much without data!\n",
    "\n",
    "We're going to see how machine learning algorithms work by applying\n",
    "them to both real and simulated data. It's critical to play with real\n",
    "data in learning machine learning, as it is very difficult to\n",
    "replicate many important features of real data via\n",
    "simulation. Simulation does play an important role in ML as well,\n",
    "however: only with simulated data can we check how our algorithms\n",
    "perform when all of the assumptions that underlie their derivation are\n",
    "truly met. It is also frequently much easier to \"turn the knobs\" on\n",
    "various data set properties of interest---like the number of sampling\n",
    "units $n$, the number of features $m$, the degree of correlation\n",
    "between features, etc.---with simulated data than in the lab or the\n",
    "external world!\n",
    "\n",
    "We will consider two real gene expression data sets:\n",
    "1.  an RNA-seq data set downloaded from Gene Expression Omnibus\n",
    "    (accession\n",
    "    [GSE120430](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE120430))\n",
    "    analyzing transcriptional targets of core promoter factors in\n",
    "    Drosophila neural stem cells ([@neves2019distinct]).\n",
    "2.  a microarray (!) data set from 2006 collected to predict\n",
    "    sensitivity to preoperative chemotherapy using expression levels\n",
    "    measured in fine-needle breast cancer biopsy specimens\n",
    "    ([@hess2006pharmacogenomic]).\n",
    "\n",
    "I'll defer further discussion of the Hess data set until we get to\n",
    "supervised analysis methods.\n",
    "\n",
    "In order to read in the data from file, I'm going to define a\n",
    "convenience function resetting some of the defaults of the\n",
    "`read.table` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "outputs": [],
   "source": [
    "rt = function(f) {\n",
    "    read.table(f, sep=\"\\t\", row.names=1, header=TRUE,\n",
    "               check.names=FALSE, comment.char=\"\", quote=\"\")\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "source": [
    "Now let's use this function to read the Neves data set, along with a\n",
    "file containing Drosophila melanogaster gene annotations, in from the\n",
    "files included here in the github project:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "outputs": [],
   "source": [
    "nevesExpr = log2(rt(\"data/gse120430_deseq_normalized.tsv.gz\") + 1)\n",
    "nevesExpr[1:5, 1:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "outputs": [],
   "source": [
    " ## (note that gene expression matrix files are usually provided\n",
    " ##  using genes-in-rows format)\n",
    " ## simplify nevesExpr by removing genes with no data:\n",
    "nevesExpr = nevesExpr[rowSums(nevesExpr) > 0, ]\n",
    " ## by contrast, sample annotation files generally follow the\n",
    " ## older statistics convention of sampling units-in-rows\n",
    "nevesAnnot = rt(\"data/gse120430_sample_annotation.tsv\")\n",
    "dmGenes = rt(\"data/d_melanogaster_gene_annotations.saf.gz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "source": [
    "Let's take a quick look at `nevesAnnot`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "outputs": [],
   "source": [
    "nevesAnnot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "source": [
    "To minimize the chance of any bugs in our analysis code, it is useful\n",
    "to align the rows of the sample annotation data (and gene annotation\n",
    "data, if we have it) to the columns of the expression matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "outputs": [],
   "source": [
    " ## align sample annotations to expression data:\n",
    "nevesAnnot = nevesAnnot[colnames(nevesExpr), , drop=FALSE]\n",
    " ## align dmGenes to expression data:\n",
    "dmGenes = dmGenes[rownames(nevesExpr), ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "source": [
    "The `group` column indicates whether each sample is in group\n",
    "expressing the control (mCherry) or one of the experimental RNAi\n",
    "transgenes (TAF9, TBP, or TRF2).\n",
    "\n",
    "The sample names in the expression data and sample annotations are\n",
    "Gene Expression Omnibus accession ids; we'll replace these with more\n",
    "descriptive names based on the grouping information in the sample\n",
    "annotations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "outputs": [],
   "source": [
    " ## use more descriptive names for samples\n",
    "betterSampleNames = paste0(nevesAnnot$group, \"-\", 1:3)\n",
    "colnames(nevesExpr) = betterSampleNames\n",
    "rownames(nevesAnnot) = betterSampleNames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "source": [
    "Finally, because the descriptive gene names for the measured\n",
    "Drosophila genes are in one-to-one correspondence with the Flybase\n",
    "gene ids used to label the rows in the file\n",
    "`data/gse120430_deseq_normalized.tsv.gz`, we'll swap them\n",
    "out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "outputs": [],
   "source": [
    " ## use more descriptive names for genes\n",
    "rownames(nevesExpr) = dmGenes$GeneName"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
   "source": [
    "The code shown above for loading in the Neves data set is also\n",
    "contained in the file `load_neves.R`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}