Differential_accessibility.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# From coverage files to differential accessibility"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check what is your working directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "getwd()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load your sample plan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "splan <-read.csv(\"data/metadata.tsv\", sep=\"\\t\",row.names=1, header=TRUE)\n",
    "splan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read coverage files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list.files('output/peak_coverage')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the path to the individual peak files. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exprs.in <-lapply(rownames(splan), function(x) { paste(\"output/peak_coverage\",'/',x,'_peaks.bed',sep=\"\")}) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then read the peak file.  \n",
    "\n",
    "Column description:\n",
    "1. chr \n",
    "2. start \n",
    "3. end \n",
    "4. peak_ID \n",
    "5. MACS_score  \n",
    "6. MACS_strand \n",
    "7. MACS_fold_enrichment \n",
    "8. MACS_-log10(pval) \n",
    "9. MACS_-log10(qval) \n",
    "10. MACS_summit (0-based offset from chrom start)\n",
    "11. reads_in_peak \n",
    "12. peak_non-zero_bases \n",
    "13. peak_length \n",
    "14. peak_fraction_covered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts.exprs <- lapply(exprs.in, read.csv, sep=\"\\t\", header=FALSE, check.names=FALSE, col.names=c('chr','start','end', 'peak','score','strand','fold_enr','-log10(pval)','-log10(qval)','summit','reads_in_peak','non-zero_bases','peak_lenght','cov_fraction'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combine the data into a single count table, using peak coordinates as an ID.\n",
    "edgeR works on a table of integer read counts, with rows corresponding to regions and columns to independent libraries. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i <- 0\n",
    "for (sample in counts.exprs) {\n",
    "    if (i==0) {\n",
    "        sample$ID <- paste(sample$chr,sample$start,sample$end,sep=\"_\")\n",
    "        count_table <- data.frame(as.numeric(sample$reads_in_peak), row.names=sample$ID)\n",
    "        }\n",
    "    else {count_table <- data.frame(count_table,as.numeric(sample$reads_in_peak))}\n",
    "    i <- i+1\n",
    "}\n",
    "## rename the columns\n",
    "names(count_table) <- rownames(splan)\n",
    "head(count_table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "message(\"size of the table\")\n",
    "message(\"number of ATAC peaks\")\n",
    "message(dim(count_table)[1])\n",
    "message(\"number of samples\")\n",
    "message(dim(count_table)[2])\n",
    "message(\"number of counts per sample\")\n",
    "print(colSums(count_table),na.rm = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Install / load edgeR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## if necessary, install edgeR \n",
    "if (!requireNamespace(\"BiocManager\", quietly = TRUE))\n",
    "    install.packages(\"BiocManager\")\n",
    "\n",
    "BiocManager::install(\"edgeR\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the edgeR package\n",
    "library(edgeR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a DGE object and normalise the counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "edgeR stores data in a simple list-based data object called a DGEList. The groups are defined in the 2nd columns of our sample plan. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y <- DGEList(counts=count_table,group=splan$group)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The calcNormFactors() function normalizes for reads composition by finding a set of scaling factors for the library sizes that minimize the log-fold changes between the samples for most regions. The default method for computing these scale factors uses a trimmed mean of M-values (TMM) between each pair of samples. The effective library size replaces the original library size in all downsteam analyses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y <- calcNormFactors(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check that the design is correct and see the normalisation factors. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y$samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Exploration\n",
    "Before proceeding with the computations for differential accessibility, it is possible to produce a plot showing the sample relations based on multidimensional scaling. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotMDS(y, col=as.numeric(y$samples$group))\n",
    "legend(\"bottomleft\", as.character(unique(y$samples$group)), col=1:3, pch=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimating the Dispersion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y <- estimateCommonDisp(y,verbose=T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y <- estimateTrendedDisp(y)\n",
    "y <- estimateTagwiseDisp(y)\n",
    "plotBCV(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <span style=\"color:red\">TODO, read the doc to check what to do with such dispersion...</span>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differential accessibility\n",
    "Once the dispersions are estimated, we can proceed with testing procedures for determining differential accessibility. The function exactTest() conducts tagwise tests using the exact negative binomial test. The test results for the n most significant tags are conveniently displayed by the topTags() function. By default, Benjamini and Hochberg's algorithm is used to control the false discovery rate (FDR)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "et <- exactTest(y)\n",
    "topTags(et)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of differentially accessible regions at FDR< 0:05 is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Z=summary(de <- decideTestsDGE(et, p.value=0.05, adjust=\"BH\"))\n",
    "Z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function plotSmear generates a plot of the tagwise log-fold-changes against log-cpm. Differential tags are highlighted on the plot. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "detags <- rownames(y)[as.logical(de)]\n",
    "plotSmear(et, ylim=c(-6,6), de.tags=detags)\n",
    "abline(h = c(-2, 2), col = \"blue\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Export the differential analysis results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "et.df = as.data.frame(do.call(rbind, et))\n",
    "write.table(et.df,\"differential_accessibility.tsv\",quote=FALSE,row.names=TRUE,sep=\"\\t\",col.names=NA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <span style=\"color:red\">TODO, read the doc to check if GLM would be more appropriate </span>."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R 4.0.3",
   "language": "R",
   "name": "r-4.0.3"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}