N4_chromVAR.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Motif analysis with chromVAR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scanpy as sc \n",
    "import pandas as pd\n",
    "import anndata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### r2py set-up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rpy2.rinterface_lib.callbacks\n",
    "import logging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ignore R warning messages\n",
    "#Note: this can be commented out to get more verbose R output\n",
    "rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I use `anndata2ri` to handle conversion of the sparse matrix (not too straightforward with `rpy2` alone)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import anndata2ri\n",
    "anndata2ri.activate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Loading the `rpy2` extension enables cell magic to be used. This runs R code in jupyter notebook cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load R packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "library(tidyverse)\n",
    "library(chromVAR)\n",
    "library(chromVARmotifs)\n",
    "library(BSgenome.Hsapiens.UCSC.hg38)\n",
    "library(SummarizedExperiment)\n",
    "library(motifmatchr)\n",
    "# Use cisBP motifs (curated set to remove redudancy)\n",
    "data(\"human_pwms_v2\")\n",
    "# Use ENCODE motifs\n",
    "data(\"encode_pwms\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load dataset\n",
    "\n",
    "AnnData generated in `N2_add_cistopic.ipynb`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "outdir = \"/home/jovyan/data/lung_adult_scATAC/\"\n",
    "adata = sc.read_h5ad(outdir + \"lungAdult__ATAC.wCisTopic.h5ad\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data setup to run chromVAR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import data in R using `anndata2ri`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_chromvar = anndata.AnnData(adata.X)\n",
    "\n",
    "adata_chromvar.var_names = adata.var_names\n",
    "adata_chromvar.obs_names = adata.obs_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "class: SingleCellExperiment \n",
      "dim: 96161 53585 \n",
      "metadata(0):\n",
      "assays(1): X\n",
      "rownames(96161): chr1:816942-817444 chr1:826600-827805 ...\n",
      "  chrY:21574860-21575895 chrY:22366897-22367799\n",
      "rowData names(0):\n",
      "colnames(53585): 59-AAACGAAAGATATGAC-1 59-AAACGAAAGGTAGGAA-1 ...\n",
      "  04-TTTGTGTTCGGATAGG-1 04-TTTGTGTTCTCGTAGA-1\n",
      "colData names(0):\n",
      "reducedDimNames(0):\n",
      "spikeNames(0):\n",
      "altExpNames(0):\n"
     ]
    }
   ],
   "source": [
    "%%R -i adata_chromvar\n",
    "adata_chromvar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transform `SingleCellExperiment` object to `rangedSummarizedExperiment` object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"chr1:816942-817444\" \"chr1:826600-827805\" \"chr1:858056-859079\"\n",
      "[4] \"chr1:869401-870227\" \"chr1:876302-877674\" \"chr1:904112-905616\"\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "## Extract peak ids\n",
    "peaks_ids <- rownames(adata_chromvar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "## String - GRanges conversion\n",
    "## Borrowed from Signac functions \n",
    "## https://satijalab.org/signac/reference/GRangesToString.html\n",
    "StringToGRanges <- function(regions, sep = c(\"-\", \"-\"), ...) {\n",
    "  ranges.df <- data.frame(ranges = regions)\n",
    "  ranges.df <- separate(\n",
    "    data = ranges.df,\n",
    "    col = \"ranges\",\n",
    "    sep = paste0(sep[[1]], \"|\", sep[[2]]),\n",
    "    into = c(\"chr\", \"start\", \"end\")\n",
    "  )\n",
    "  granges <- makeGRangesFromDataFrame(df = ranges.df, ...)\n",
    "  return(granges)\n",
    "}\n",
    "\n",
    "GRangesToString <- function(grange, sep = c(\"-\", \"-\")) {\n",
    "  regions <- paste0(\n",
    "    as.character(x = seqnames(x = grange)),\n",
    "    sep[[1]],\n",
    "    start(x = grange),\n",
    "    sep[[2]],\n",
    "    end(x = grange)\n",
    "  )\n",
    "  return(regions)\n",
    "}\n",
    "\n",
    "# Extend genomicRanges\n",
    "# \n",
    "extend <- function(x, upstream=0, downstream=0)     \n",
    "{\n",
    "    if (any(strand(x) == \"*\"))\n",
    "        warning(\"'*' ranges were treated as '+'\")\n",
    "    on_plus <- strand(x) == \"+\" | strand(x) == \"*\"\n",
    "    new_start <- start(x) - ifelse(on_plus, upstream, downstream)\n",
    "    new_end <- end(x) + ifelse(on_plus, downstream, upstream)\n",
    "    ranges(x) <- IRanges(new_start, new_end)\n",
    "    trim(x)\n",
    "}\n",
    "\n",
    "## Make rangedSummarized experiment object\n",
    "atac_se <- SummarizedExperiment(assays = list(counts=assay(adata_chromvar, \"X\")), rowRanges = StringToGRanges(peaks_ids, sep=c(\":\",\"-\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Match peaks to motifs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "class: RangedSummarizedExperiment \n",
      "dim: 96161 870 \n",
      "metadata(0):\n",
      "assays(1): motifMatches\n",
      "rownames(96161): chr1:816942-817444 chr1:826600-827805 ...\n",
      "  chrY:21574860-21575895 chrY:22366897-22367799\n",
      "rowData names(1): bias\n",
      "colnames(870): ENSG00000008196_LINE2_TFAP2B_D_N1\n",
      "  ENSG00000008197_LINE6_TFAP2D_D ...\n",
      "  ENSG00000112837_LINE19949_TBX18_I_N1\n",
      "  ENSG00000122145_LINE20002_TBX22_I_N1\n",
      "colData names(1): name\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "## Calculate GC bias across genome\n",
    "atac_se <- addGCBias(atac_se, genome = BSgenome.Hsapiens.UCSC.hg38)\n",
    "## Calculate GC bias across genome\n",
    "motif_ix <- matchMotifs(human_pwms_v2, atac_se, genome = BSgenome.Hsapiens.UCSC.hg38, p.cutoff=5e-5)\n",
    "rownames(motif_ix) <- peaks_ids\n",
    "rownames(atac_se) <- peaks_ids\n",
    "\n",
    "motif_ix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute motif deviations (i.e. how much does TF accessibility deviate from expected value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Error in get(as.character(FUN), mode = \"function\", envir = envir) : \n",
      "  object 'as.SimpleList' of mode 'function' was not found\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "bg <- getBackgroundPeaks(object = atac_se)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "head(motif_ix)\n",
    "motif_dev <- computeDeviations(object = atac_se, annotations = motif_ix, background_peaks = bg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "motif_dev"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (emma_env)",
   "language": "python",
   "name": "emma_env"
  },
  "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}