Lecture-7-DataPreprocessing2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "As mentioned in the last notebook, it is necessary to inspect your data before analysis. For example, are you sure you know which file delimiter was used to create the file? Should the data be transposed? Data preprocessing steps also include inspecting the covariates and plotting data to search for outliers that may indicate measurement (or data entry) errors. \n",
    "\n",
    "We also discussed steps that can be used to determine whether the data fit classical statistical assumptions. This notebook continues in this area. In particular, we introduce generalize linear models (GLMs) that can be used to analyze non-normal data, such as count data or percentages.\n",
    "\n",
    "For this notebook, we will use the glucosinolate data that we downloaded earlier.\n",
    "\n",
    "### Load the glucosinolate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rm(list=ls());\n",
    "# open the glucosinolate file in R\n",
    "# same file as before...\n",
    "glucosinolateFileName <- \"data/cmeyer_glucs2015/bmeyer_etal.txt\";  \n",
    "glucs <- read.table(glucosinolateFileName, header=T, sep=\"\\t\", as.is=T, stringsAsFactors=FALSE);  \n",
    "glucs <- glucs[order(glucs[,\"accession_id\"]),];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# what's in the working environment?\n",
    "ls();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What is the statistical model?\n",
    "\n",
    "During today's notebook, we will discuss *Poisson* family analyses. Poisson models (and their derivatives) are used to analyze count data, which are very common.\n",
    "\n",
    "Examples:\n",
    "  - The number of microbial species in each of our gut microbiomes \n",
    "  - The number of people that die from hippo attacks each year\n",
    "  - The number of reads assigned to a gene in an RNA-Seq dataset\n",
    " \n",
    "The glucosinolate data are also count data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "str(glucs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(glucs);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if lme4, ggplot2, and gridExtra aren't installed, install them...\n",
    "if( !require(\"lme4\" )){  \n",
    "    install.packages(\"lme4\");  \n",
    "}\n",
    "\n",
    "if( !require(\"ggplot2\" )){  \n",
    "    install.packages(\"ggplot2\");  \n",
    "}\n",
    "\n",
    "if( !require(\"gridExtra\" )){  \n",
    "    install.packages(\"gridExtra\");  \n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Replicates in the data\n",
    "\n",
    "As noted earlier, the glucosinolate data were generated with *Arabidopsis thaliana*, the plant genetic model species. It is among the main model species used by plant geneticists in part because it is self compatible. This means that replicated inbred lines can be used in experiments to reduce experimental noise. Importantly, the mean phenotype of an inbred line can typically be measured with higher precision than for an outbred line.\n",
    "\n",
    "### Estimate the mean phenotype per genotype\n",
    "\n",
    "The speed of GWAS is affected by the sample size and the number of SNPs, as well as other factors. To speed up analyses, researchers working with repeated measures data (e.g. longitudinal analyses) or inbred lines often reduce the number of observations for each sample or individual. One approach, as illustrated in the previous notebook, is to simply average the phenotypic data per individual. However, this is not ideal (i.e. the errors propagate when analyzing the \"means of means\")."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# recall, there is more than one observation per inbred line\n",
    "counts <- table(glucs$accession_id);\n",
    "cat(\"The distribution of the number of replicates per accession\\n\");\n",
    "table(counts);\n",
    "\n",
    "# earlier, we visualized this as a barplot:\n",
    "# ggplot() + aes(counts) + \n",
    "#         xlab( paste0( \"The number of replicates\" )) + ylab(\"Counts\") +\n",
    "#         geom_bar(stat=\"count\", fill=\"tan1\") + \n",
    "#         geom_vline(aes(xintercept=mean(counts)), linetype=3);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But how should we extract the mean phenotype per individual?\n",
    "\n",
    "Recall that, earlier, we used linear mixed models to fit BLUPs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# previously, we used a mixed model to specify a random effect with the following code:\n",
    "lmer0 <- lmer( G2P ~ 1 + (1|accession_id), data=glucs );\n",
    "linearBlups <- ranef(lmer0)$accession_id; # these are the best-linear unbiased predictors:\n",
    "linearBlups <- data.frame( accession_id=rownames(linearBlups), linear_blup=linearBlups[,1], stringsAsFactors=FALSE );\n",
    "head(linearBlups);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a standard approach with mixed models. However, you'll remember that these are not really *normally distributed* data.\n",
    "\n",
    "Indeed, a lot of the data that biologists work with are non-normal (e.g.):\n",
    "  - number of offspring \n",
    "  - infection rates\n",
    "  - survival/death\n",
    "  - gene expression data\n",
    "  - metabolomic data\n",
    "\n",
    "In general, it is appropriate to analyze such data with generalized linear models (or GLMs). In the case of *count* data, it is typical to use *Poisson* (or Poisson-family) models. An easy rule-of-thumb to remember is that if the data are all non-negative integers, with no natural upper bound, then the Poisson-family models should be used.\n",
    "\n",
    "\n",
    "<!--# however, these are count data (see the output from head above), which suggests \n",
    "# that a Poisson-family model should be used; let's try a simple quasi-Poisson GLM:\n",
    "glm0 <- glm( G2P ~ 1, data=glucs, family=\"quasipoisson\" );\n",
    "glm0.res <- residuals( glm0 );\n",
    "glm0.res[1:10];\n",
    "\n",
    "# note: the names are no longer the accession Ids, but the row names from the glucs data frame. Do you know why?\n",
    "# let's use a workaround:\n",
    "################################################################################\n",
    "## my version of stack, which avoids factor generation\n",
    "################################################################################\n",
    "mstack <- function(arg, newHeaders, setRowNames=T, sorted=TRUE, decreasing=F){\n",
    "    values <- data.frame(names=I(names(arg)), values=as.numeric(arg));\n",
    "\n",
    "    if( setRowNames ){\n",
    "        rownames(values) <- values[,\"names\"];\n",
    "    }\n",
    "    \n",
    "    if( sorted ) {\n",
    "        values <- values[order(values[,\"values\"], decreasing=decreasing),];\t\n",
    "    }\n",
    "    \n",
    "    colnames(values) <- newHeaders;\n",
    "    return(values);\n",
    "}\n",
    "\n",
    "glm0.res <- mstack( glm0.res, newHeaders=c(\"row_id\", \"residual\"), sorted=FALSE );\n",
    "head(glm0.res);-->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the estimates from Meyer et al., however, are areas estimated under the curve, from a QQQ Mass spectrometer.\n",
    "# these can be rounded without a loss of precision \n",
    "glucs[,3:ncol(glucs)] <- round(glucs[,3:ncol(glucs)]);\n",
    "glucs[1:5,];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now, these integers can be fit with GLMs/GLMMs:\n",
    "glmer0 <- glmer( G2P ~ 1 + (1|accession_id), data=glucs, family=\"poisson\" );\n",
    "glmBlups <- ranef(glmer0)$accession_id; # these are the best-linear unbiased predictors:\n",
    "glmBlups <- data.frame( accession_id=rownames(glmBlups), pois_blup=glmBlups[,1], stringsAsFactors=FALSE );\n",
    "head(glmBlups);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For poisson/negative binomial GLMs the distribution of these random effects (or residuals, if a GLM is used) will tend to be skewed, which is to be expected. For larger counts (and larger $\\lambda$) the residuals will tend to normality.\n",
    "\n",
    "### Poisson-family models allow offsets\n",
    "\n",
    "When the counts from a Poisson process should be modeled as a rate, then an offset should be included in the model. Specifying an exposure will thus enable you to express the count data as a function of the *effort* that was necessary to collect the data. For example, if you detect 25,000 counts for G2P in a sample with a high background on a mass spec (the instrument that Meyer et al., used to measure glucosinolates), then you should take that higher *ion count* into **account** when comparing that sample with samples that have less of a background. \n",
    "\n",
    "As another example, if two botanists count the number of species at the botanical gardens and one botanist spends 24 hours counting species and the other 72 hours (assuming the same skill level), then you would typically expect the second botanist to count **more species**, given the extra effort.\n",
    "\n",
    "But how do you specify an offset?\n",
    "<!--log𝜇𝑥 = $\\beta$0 + $\\beta$1𝑥 \n",
    "\n",
    "(where 𝜇𝑥 is the expected count for those with covariate 𝑥), you have\n",
    "\n",
    "log𝜇𝑥𝑡𝑥=𝛽′0+𝛽′1𝑥\n",
    "(where 𝑡𝑥 is the exposure time for those with covariate 𝑥). Now, the last equation could be rewritten\n",
    "\n",
    "log𝜇𝑥=log𝑡𝑥+𝛽′0+𝛽′1𝑥\n",
    "and log𝑡𝑥 plays the role of an offset.-->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# determine the total ion count\n",
    "totalIonCounts <- rowSums( glucs[,-c(1,2)] );\n",
    "head(totalIonCounts);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now include the total ion count as an offset\n",
    "# if neither df has been sorted, you can immediately add it to the model without merging with the original df\n",
    "glm0 <- glm( G2P ~ 1, offset(log(totalIonCounts)), data=glucs, family=\"quasipoisson\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wait a minute, negative weight? The only weighting comes from the offset, which is expressed as a log in Poisson-family models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "which(log(totalIonCounts) < 0);\n",
    "which(totalIonCounts == 0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quality control (data curation) is an ongoing process\n",
    "\n",
    "There are samples in the data that have *no* ion counts. If this wasn't a simple mistake during data collection (e.g. failing to input the data into the file), then the samples could be rerun or dropped. We have a lot of data, and you're a busy person, so let's drop them.\n",
    "\n",
    "<!-- # glm0 <- glm( G2P ~ 1, offset(log(totalIonCounts)), data=glucs2, family=\"quasipoisson\");\n",
    "# summary(glm0);-->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dropouts <- which( totalIonCounts == 0 );\n",
    "glucs2 <- glucs[-dropouts,];\n",
    "dim(glucs2);\n",
    "\n",
    "totalIonCounts <- rowSums( glucs2[,-c(1,2)] );\n",
    "cat(\"The range of ion counts is now:\", range(totalIonCounts), \"\\n\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# 20 is a small number, but the advantage of adding an offset is \n",
    "# that you can include as much data as possible, and let the model weight things appropriately\n",
    "\n",
    "glmer0 <- glmer( G2P ~ 1 + offset(log(totalIonCounts)) + (1|accession_id), data=glucs2, control=glmerControl(optimizer=\"bobyqa\", optCtrl=list(maxfun=2e5)), family=\"poisson\" );\n",
    "glmBlups <- ranef(glmer0)$accession_id; # these are the best-linear unbiased predictors:\n",
    "glmBlups <- data.frame( accession_id=rownames(glmBlups), pois_blup=glmBlups[,1], stringsAsFactors=FALSE );\n",
    "tail(glmBlups);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These BLUPs can be used in GWAS."
   ]
  }
 ],
 "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": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}