{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Overview \n",
    "\n",
    "\n",
    "\n",
    "This notebook is an example R script on how to prepare the input data prior to building a base GRN.\n",
    "Here, we use Cicero to extract the cis-regulatory connections between scATAC-seq peaks.\n",
    "\n",
    "\n",
    "\n",
    "### Notebook file\n",
    "This notebook is available on CellOracle’s GitHub page as this jupyter notebook (with R kernel) or an R notebook. The notebooks are identical. Please use whichever one you prefer.\n",
    "\n",
    "- Jupyter notebook with R kernel: https://github.com/morris-lab/CellOracle/blob/master/docs/notebooks/01_ATAC-seq_data_processing/option1_scATAC-seq_data_analysis_with_cicero/01_atacdata_analysis_with_cicero_and_monocle3.ipynb\n",
    "- R notebook: https://github.com/morris-lab/CellOracle/blob/master/docs/notebooks/01_ATAC-seq_data_processing/option1_scATAC-seq_data_analysis_with_cicero/01_atacdata_analysis_with_cicero_and_monocle3.Rmd\n",
    "\n",
    "\n",
    "\n",
    "### CAUTION: \n",
    "\n",
    "- This notebook is intended to **demonstrate data preprocessing steps prior to starting a CellOracle analysis**. CellOracle is NOT used in this notebook, and this notebook is not the CellOracle analysis.\n",
    "\n",
    "- Here, we will use `Cicero` to process scATAC-seq data. If you are new to this packages, pelase review the Cicero's documentation to learn the basic process of Cicero in advance.  \n",
    "\n",
    " - `Cicero` documentation: https://cole-trapnell-lab.github.io/cicero-release/docs_m3/ \n",
    "\n",
    "- The R library, cicero and monocle3 is NOT the part of celloracle package. **Please install them yourself if you use this notebook**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 0. Import library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(cicero)\n",
    "library(monocle3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Download data\n",
    "\n",
    "This tutorial uses fetal brain scATAC-seq data from a 10x Genomics database. If you’re using your own scATAC-seq data, you will not need to download this dataset.\n",
    "\n",
    "You can download the demo file with the following command.\n",
    "\n",
    "**Note: If the file download fails, please manually download and unzip the data.**\n",
    "http://cf.10xgenomics.com/samples/cell-atac/1.1.0/atac_v1_E18_brain_fresh_5k/atac_v1_E18_brain_fresh_5k_filtered_peak_bc_matrix.tar.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create folder to store data\n",
    "dir.create(\"data\")\n",
    "\n",
    "# Download demo dataset from 10x genomics \n",
    "download.file(url = \"http://cf.10xgenomics.com/samples/cell-atac/1.1.0/atac_v1_E18_brain_fresh_5k/atac_v1_E18_brain_fresh_5k_filtered_peak_bc_matrix.tar.gz\",\n",
    "              destfile = \"data/matrix.tar.gz\")\n",
    "# Unzip data\n",
    "system(\"tar -xvf data/matrix.tar.gz -C data\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can substitute the data path below to your scATAC-seq data.\n",
    "data_folder <- \"data/filtered_peak_bc_matrix\"\n",
    "\n",
    "# Create a folder to save results\n",
    "output_folder <- \"cicero_output\"\n",
    "dir.create(output_folder)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Load data and make Cell Data Set (CDS) object \n",
    "## 2.1. Process data to make CDS object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in matrix data using the Matrix package\n",
    "indata <- Matrix::readMM(paste0(data_folder, \"/matrix.mtx\")) \n",
    "# Binarize the matrix\n",
    "indata@x[indata@x > 0] <- 1\n",
    "\n",
    "# Format cell info\n",
    "cellinfo <- read.table(paste0(data_folder, \"/barcodes.tsv\"))\n",
    "row.names(cellinfo) <- cellinfo$V1\n",
    "names(cellinfo) <- \"cells\"\n",
    "\n",
    "# Format peak info\n",
    "peakinfo <- read.table(paste0(data_folder, \"/peaks.bed\"))\n",
    "names(peakinfo) <- c(\"chr\", \"bp1\", \"bp2\")\n",
    "peakinfo$site_name <- paste(peakinfo$chr, peakinfo$bp1, peakinfo$bp2, sep=\"_\")\n",
    "row.names(peakinfo) <- peakinfo$site_name\n",
    "\n",
    "row.names(indata) <- row.names(peakinfo)\n",
    "colnames(indata) <- row.names(cellinfo)\n",
    "\n",
    "# Make CDS\n",
    "input_cds <-  suppressWarnings(new_cell_data_set(indata,\n",
    "cell_metadata = cellinfo,\n",
    "gene_metadata = peakinfo))\n",
    "\n",
    "input_cds <- monocle3::detect_genes(input_cds)\n",
    "\n",
    "#Ensure there are no peaks included with zero reads\n",
    "input_cds <- input_cds[Matrix::rowSums(exprs(input_cds)) != 0,] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Qauality check and Filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAADAFBMVEUAAAABAQECAgIDAwME\nBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUW\nFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJyco\nKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6\nOjo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tM\nTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1e\nXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29w\ncHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGC\ngoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OU\nlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWm\npqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4\nuLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnK\nysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc\n3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u\n7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////i\nsF19AAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO3deYAT9d348c/uciyC3MhyI6Ci\nokVBWY4iIioi4oUg3lREtF591B5qKT5aa9U+tdZaba1Hta31qNpf6+NZsZ71Ea0oj4paLSrQ\nxy0Cys3u9zczmSSTY7OZfD9JzM779cfOZHcy35kk791MMhAxAKxJuTcAaA0ICVBASIACQgIU\nEBKggJAABYQEKCAkQAEhAQoICVBASIACQgIUEBKggJAABYQEKCAkQAEhAQoICVBASIACQgIU\nEBKggJAABYQEKCAkQAEhAQoICVBASIACQgIUEBKggJAABYQEKCAkQAEhAQoICVBASIACQgIU\nEBKggJAABYQEKKigkO4QkZdis8NFJhhzt/ON7WXamG2LhrTteFf80o3OlkjVGm++pzt/WkvX\nt932m0U6eTPrr5o8rEPfcV/7e8Gr+ol0bLDYkpTVFPcuedRZe9qm/kWq/haclg0hFeYnbi63\nxi95Icl/u7PvS0lDure7xFTN/KKwNTV0kwssNiR1NSUPyewn9U3Babm0tpC2O9+7pwQbc5BI\n968/F78UC2mhO/vbrCFlbJZSSEvbuA31qHbHPK6wNZ0tbVdYbEjqagrZrfzvsmwh3S9yV3Ba\nLpUc0vt33HFH+m+hUoW0u8h5yUuxkA5xZy/IL6Rs2x6GH9J5IjXXbzDbnhjhDLC0kBV91EZO\ntNiOtNUUslt2ITXuLEMaA9NyqeSQsilVSM4WXJa85IY0QLq4d+RY6VOTR0i2/JAmiJzkXV7h\n/FH6r0JWdInI4wrbY7Eau5DMQpH/F5yWSSWHFD9OaHrw0CG1Qw6+23lOMTN2wOAu1nD5tEGd\n9z8tfvd+MLtXxwP/9mZ9/XjnwiKRYeb+kbsas+32Awa1HzDuBvcA48f19fM3XfiVHfa6sWnL\n5bvUDp27Mjh8cH3+MIviP3NDmi3ypjFba+WYNrGQAqtObFZi4Ni2v1Qd+41wvvPDm435xaJF\nv3BXF5+ad8+b2L1u0vXbMjYgvu+DRBbENuHi0067z5h5IpPcS7eKtGl2jwK3l9ncS/rEfpE/\neerIDnVHP2HSt+qq+vqzv7igX+1u3/3cpNx02Vbjb9at9fUzG3+8Z4fBM5e710lbRepGBu6y\nFCn7/ubRPXYY82A8pODIZrn/VCA+LZNWEFLT0f7x9qgvAvfKX3v63z3Zuyue8S7WXht7gLmP\nhjtFBpst4/yl9t7gPSvbM3b5G5O9Sb+1ydFT1pctpOu9Fx/+R+SaWEjBVaeE5A3sb/uFznf/\naD5s5zy0nGdEo0VGu6uLT+9tH7veHg3pGxC//iSR9td/ktzKtJCy7lHw9jKPi5zlXfMC/5vn\npm+Vs8rxE7wfDV8dvOmyrsbfLGehoSd5P6x9M7ZVwVXkFVLKvj/QyZufFQspZWRj9pbqtcFp\neVRYSEnJkH7qXBp29P7O19PNW4vdx/dz682Kbs7M4P1rna/fdpZZ5z4Ke3eWNomQend3Hw2X\nOt/e7YBeztf/9B9LO/gvg9XVxL7pS13fW88NFDnlucRBuhvSizXO+OZnIs/GQgquOrFZiYH9\nbd84TGRE46nOsO+ZjJDedB5LNSNHOlc9In0D4tf/hbupVWO/dc97sQ3JDClzj4K3l/uU7Nfu\nFR52Lo6ZO8r5em/aVjmrrJKqwe6rGicHb7qsq0mG5Fyrj3udw0zGKlI3MnHbpEjZ9w/cqPp0\n9vaiIW0HjDlT5E/BaXm0gpAOF+8413kw90w+4V4gUu3ctavHOr+yV3mP6ra/M43XVSVCkn63\nL33LjIg9LOeITIs97L653VzuTAa+bVb0FTk6MXja+rIcI705UvYw5hRpszEWUsqq45uVGDj+\ntPQZZ4sujh/fpIV0okgv5xe68/tflqdvQPyv8YIq//YYeqH7NlaWkDL2KOX2Gi/yjnNh61Bn\nQWf6LZEhW1K3ylml7LPCrHWuVr0ssAfZVhMMafr/mYb9RLpkriJtI7MeI6Xsu3Mb7viIafyh\nH1LKyMbc5mxrcFoerSCkvZw7/77PzecPPfTQ1sS9sqvzF8O90lvOw+EBY4b6v72OSYb0vDNt\nuuXmm51HZdMx3iPXedh1d552f+L88Abnh3ODr2ikrS9bSGdJ1Vqzm7MiL6TUVQdDcgdOvqH6\ndW9vxmR7wcn5q3e5Ox3ZqdPt6RuQuP4LCwb7t0jvzMdo1j1Kub36i2x2frLUeYi7Ia51tn1x\n6la5f07cd3sbOopcHdiDbKsJhFTzmTO9K/bIT1tFPiEF9317O//CQbHVpYzsvhcrs0xgWh4V\nFlK2Y6S57n3e9oDvv+y+7urfK1ucB8T93qLOA/D7ZovzrOZh99J9iZC6+Gv96Nf/cajzm9oP\naR/nOw3OpT8b749QIqS09WUNyXnUPLa2Ss4x/osNwVUHQooNnAjh80HOj9oty7K/m5w/C081\ntwGJ6ztWPHjJBPdp2wHZQsrco5Tbq1Y6+LdLwi9Tt8pZ5UBvmCkiXwvsQbbVBEIa7F58wllk\ndcYq8ggpZd/fEf+l/Z/GQkoZ2ZjXRA4MTsujFYS05rjq2CNgtycT98o/nMmL3qKT3b9Fy52L\nr7iXXk6ENNT76fsHedfsGA/JfUblPuweNakhpa0va0jvO0/1nYPu3/ghpaw6EFJs4GQI7hOW\n5FPIAPfx82pzGxAMyduRvZ3f+huzhJS5R8Hba7P/9OhHgZAuTt0qZ5VjvJlTvJfFEnuQbTWB\nkIa5F59MhBRcRR4hpez7Y86FT92ZP8RCSrnDjXk39tsiMS2PVhCS86v/xwe1dW/Y9h/E75XN\nNbFnYN6i3zOfiX8Cz8PJV+3ci1udZwk9z7rnw++2FFLa+rKGZHaSqVeIvB8LKXXVgZC8gZPb\n3rCT86OabOeJrZPAezPpGxC7/tv19fX+Ir8T9zBlnvd3yftxsyGl3F610t79zu+do5DnfO+n\nbpWzykHeCAeLzAvsQbbVNB9ScBVpG5ktpJR9f9O58IY7c5P47yMFRuYvUmjNhLTGue+bzPr7\nnUNQ+XniXhkmMtdd8m3nAXifMT1EvuFenJca0hKJHSQf21JI6evLGtJR0nWa7GRiIaWuuvmQ\n5ni/W/fYnGWH6/yzji5bsODp9A2IXX+1+EdOxtzgPcjOFNndvfT15kNKvb38gxvnQSj/Dgwd\n2Cr3AMd9GDd0EvlRcg+yrqb5kIKrSNvIrMdIwX3f5PwB+oF74XBvH1NHNuZp/+yopws+S0pD\n5YfkniX6B+d7n7d3z7ba7t+6Zzi/UH9rzL/GOU/1P/HerOjg/Em6rSY1JPcdvv9xnsrXtBhS\n2vqyhnSN88xdZvghpa46vlkZIT3k/OBHzsH0Je6FtDdknUdzN+f5za+8X8hpG+Bffz/nZ+f8\n0zmEuqeLyChjrhYvswfaNR9S6u01XuRt58Jm51juW870vRHDh7+SulXuS277rTLrpzurfDu5\nB1lX03xIwVWkbWT8tkmRsu8znOOyp2K/K5yQUkc25naRi0xgWh6VH5J72lubkcdNdp6M1Hxg\njHNMss+Nq80HzgNLdhnvXPJehH7LPRgf2MO9J4Ihfexcbrv/SPcl5FG5Q0pbX9aQnnVX/wM/\npNRVxzcrPaQ1zm/ece7AbZaYjJe/P+zgrsK5IIdmbIB//Wfc92akU29396qc50LPuBfrusb3\nM+sepdxel4jc6Q7pvr42Yf4hnWJHUcGtciuQml2cR733nmtiD7KupvmQAqtI28j4bZMiZd9f\nd2/FgbG3wxrS73B3r7xzgxaU9RyhVhDS8p3E5/5eO1Riiz0d//cFc70TSW70HnHVR6eGFHsB\nSIac6vz+25gzpPT1ZQtpk/u8/Wk/pNRVxzcrPaSTxX0xucGJ5CtbM89suL1NbMiRazI2IP4X\n7Y+D4jvf+7fu5eO8+fZzcoSUcns97p9j1HSy/72xDWlb5VSwz87ejw79LHjTZV1NsyEFV5G2\nkYm7LEXKvnvvIEnVYbGQUu9wM1KqPwtOy6MVhGTW3zBxSIfue8f+bduKY3u16eK+QvevhVMH\n7rjfaU/6135u5qCehzz6UlpIm6/eveO+/7HuldGjR/8qd0hp68sWkhnj/I78Ih5Syqrjm5UW\n0p+c6x3jXLpKvPON0kMySxeM69rnwFu2Z25A4hhr66+OHNm1817Tf7bJu7j92n07dj/y7w+N\nHu2+TpZ9j4K31+ZeUhdb/QOz96gddMjvtqZvlfsa28cn9Nhx8o+897oSe5B1Nc2GFFxF2kYm\n77IUKfv+xCFde0x98AXntlybfoe/J3JwcFomFRSShkfjr9PC5zwpeyznAvEXqy1Wk98qCnS5\ne15gYFomEQnpouOPv8adzheZWu5t+XJZUSNzci6QXwU5V1PMkJqGyc6NgWm5RCSks53jo4ue\nWnym81Tjd+Xeli+Zs6TtP3P9PM8Kcq2mmCE96J8uG5+WS0RCWl/vH55Wfafcm/Jl09BNzs/1\n8zwryLWafFbxbQka1tLiCWNkTFNwWi4RCcls+/1hu3ToNfr0bCe1RVwL/4tQvn9Ocqwmn1X8\n9vignGkHPS1VLwWnZROVkICiIiRAASEBCggJUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJ\nUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJ\nUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJ\nUEBIgAJCAhQQEqCAkAAFhAQoICRAgV1ITQ0rG5U2BKhkFiEtPrFfW5Ga/rMX620OUJkKDmnT\nVJG+Y6ZNqx8gMn2z5iYBlafgkBbK1Ndic8vmyBVamwNUpoJDqh++LT7bNHGczsYAlargkDqf\nmpy/pLPClgAVrOCQxu6+PTE/eazKtgAVq+CQFsnhb8Tmlp8slxe0ipVPhPEkr2jgy6vwV+2m\niQycMOPIiUNEDivsMX56u24hVN9f6KYCRWfxPtLTc/rUiNT0mfVUgSs47bQwS/e+p8BhgOKz\nO7OhcfUqizMbCAmtRjlPESIktBrlPEWIkNBqlPMUIUJCq1HOU4QICa1GOU8RIiS0GuU8RYiQ\n0GqU8xQhQkKrUc5ThAgJrUY5TxEiJLQanCIEKCjOKUIrdh+SMLDz9myLGEJCK2L733E1vrNs\nW+Z3t9x+S8I3ZUsz1yUktBoFh3TZr5wvW6/YQaTd19bkWvB5QkLrV3BIMsn5Ml+6HXvmWNl1\nY44FCQkRYBXS61X7f+rM3imX5ViQkBABViHdLM978+NH51iQkBABViEtlM+9+QUdcyxISIgA\nq5DukqXe/FEDcixISIiAwkPqe8V9L/ea2eTMvtDm2BwLEhIioOCQBlSJ6xFjzm1f+1qOBQkJ\nEVD4G7Iblz5w9elffcKYur3+lms5QkIEKHzQ2Hu5f0xIiIDif2IfISECCAlQQEiAAkICFBQc\nUtcUORYkJERAwSHdtIfIHiPicixISIiAwp/abRguef1XDYSECLA4RvoBIQE+i5AerSUkIIZX\n7QAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCA\nkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCA\nkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCA\nkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCA\nkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCA\nkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCA\nkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCA\nkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCA\nkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqDA\nLqSmhpWNLS1DSIgAi5AWn9ivrUhN/9mLcy5GSIiAgkPaNFWk75hp0+oHiEzfnGNBQkIEFBzS\nQpn6Wmxu2Ry5IseChIQIKDik+uHb4rNNE8flWJCQEAEFh9T51OT8JZ1zLEhIiICCQxq7+/bE\n/OSxORYkJERAwSEtksPfiM0tP1kuz7EgISECCn/VbprIwAkzjpw4ROQwXrVDxFm8j/T0nD41\nIjV9Zj2VczFCQgTYndnQuHoVZzYAnCIEqOAUIUABpwgBCjhFCFDAKUKAAk4RAhRwihCggFOE\nAAWcIgQo4BQhQEGRThH6+ysJtxESWr/inCL0XhsJaO6JHyGh1eAUIUABpwgBCjhFCFDAKUKA\nAk4RAhRwihCggFOEAAWcIgQo4BQhQAH/ixCggE/sAxQQEqCAkAAFhAQoKDikrilyLEhIiICC\nQ7ppD5E9RsTlWJCQEAGFP7XbMLzZf7GXgpAQARbHSD8gJMBnEdKjtYQExPCqHaCAkAAFhAQo\nICRAASEBCggJUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCAkAAFhAQo\nICRAASEBCoIh3bGuGCMQEiIgGJLUHnPvRvURCAkREAzpZwdUS6eT/rRVdwRCQgSkHiOtutFp\nqfsZf2nxY/hCICREQMaLDatunFgtfc5/SW0EQkIEZL5q9/dFO4tjt4eVRiAkREBqSNv+cv4g\nkbozH3v1wk5Vf9EZgZAQAcGQ7j+5m8jQC59vci+8KmfrjEBIiICUl79l70Wvxy+s63mtzgiE\nhAgIhnTd+8UYgZAQAanHSJ+4x0V3vqk6AiEhAoIhbb+oarQzGSzn8j4SEEowpJuk/o/O5K8z\n5JeKIxASIiAY0t7DYg/5pn1GK45ASIiAYEgdz/RnztlRcQRCQgQEQ9ptmj8zfVfFEQgJERAM\n6fSaP3rTx2pOURyBkBABwZA+HSCHXHnrD4+q6vmJ4giEhAhIeR/pgzlV7vmqh6q+kURIiIC0\ns7//9exvnvyn7giEhAjgPz8BFKSEdO/sKTFzFUcgJERAMKRfinTr6dlZcQRCQgQEQ9pz/xVF\nGIGQEAHBkGofK8YIhIQICIbU/8lijEBIiIBgSAtnF2MEQkIEBEPadvIRT33yhUdxBEJCBARD\n6tpF4hRHICREQDCZeUmKIxASIoAzGwAFaSFteuNF7REICRGQEtI/Z7VzDo9umql62iohIQKC\nIa0aKBMOFnNfmz4fK45ASIiAYEjnyK3mbucbL7VfoDgCISECgiENmmi8kMzMXRRHICREQPr/\nIuSFdG5HxREICREQDGn/0X5I40cpjkBIiIBgSFfI5Y1uSD+RbymOQEiIgJRz7SbIsLFy+l6y\np+ZnmxMSIiDlfaQt/zVARHpcsl5zBEJCBKSfIvT5sn8rj0BIiADOtQMUBEM6KeE6xREICRGQ\n+hmyvoFnKI5ASIiAYEibXZs+emjUxA2KIxASIiDbMdL6YRcojkBIiICsLzZc3FdxBEJCBGQN\n6fwOiiMQEiIgS0hNizvvrTgCISECgiF1imkncqfiCISECAiGNN136sOaIxASIoAzGwAFhAQo\nSPlP9FPsrzQCISECgiEt6CdSN6p/lQye4DhIaQRCQgQEQ3q2+qA3nMk7h/X7UHEEQkIEBEM6\nYlDsHLtNQ2cqjkBIiIBgSL1P9WdO7684AiEhAoIhDZzszxzcR3EEQkIEBEOaXf2QN/1T9XTF\nEQgJERAM6R/dqmff9t+3n1Dd/lXFEQgJEZDyhuySA7x/ILv7o5ojEBIiIO3MhqX3/ejXL2xX\nHYGQEAF80BiggA8aAxTwQWOAAj5oDFDAB40BCvigMUABHzQGKOCDxgAFfNBYaN+fEsa0T8u9\nvSgFPmgstPETv5W/C2RJubcXpcAHjYU2/soQC68jpGgIhPTJzc8XYwRCQgQEQlosxxZjBEJC\nBARC2rJnz2IcGEc7pAbpNySEaUXbahRX8Bhp7RGjHv5w/RcuxRGiHdLHct69+bugd9G2GsUV\nDKlup8SHXyqOEPWQfhNi6XsIqVIFk5mXpDhCBYT0ux7dQmgzK8SqCSki4iGdo/lJLikqIKQr\ndwvx9OveDlNCrJqQIiIekpzkfr1N80+RrxJCGh9m6c6EhAypIZ1WhA+nIKT8EVLFIiRDSLBH\nSIaQYI+QDCHBHiEZQoI9QjKEBHuJkAYd79hZjo9RHIGQ8kdIFSsRUirFEQgpf4RUseLJvJJK\ncQRCyh8hVawiHBSlIaT8EVLFIiRDSLBHSIaQYI+QDCHBHiEZQoI9QjKEBHuEZAgJ9gjJEBLs\nEZIhJNgjJENIsEdIhpBgj5AMIcEeIRlCgj1CMoQEe4RkCAn2CMkQEuwRkiEk2CMkQ0iwR0iG\nkGCPkAwhwR4hGUKCPUIyhAR7hGQICfYIyRAS7BGSISTYIyRDSLBHSIaQYI+QDCHBHiEZQoI9\nQjKEBHuEZAgJ9gjJEBLsEZIhJNgjJENIsEdIhpBgj5AMIcEeIRlCgj1CMoQEe3YhNTWsbGxp\nGULKHyFVLIuQFp/Yr61ITf/Zi3MuRkj5I6SKVXBIm6aK9B0zbVr9AJHpm3MsSEj5I6SKVXBI\nC2Xqa7G5ZXPkihwLElL+CKliFRxS/fBt8dmmieNyLEhI+SOkilVwSJ1PTc5f0jnHgoSUP0Kq\nWAWHNHb37Yn5yWNzLEhI+SOkilVwSIvk8Ddic8tPlstzLEhI+SOkilX4q3bTRAZOmHHkxCEi\nh/GqXXMIKSIs3kd6ek6fGpGaPrOeyrkYIeWPkCqW3ZkNjatXcWZDTiFD2mlNGC3e9igZThEy\nX6KQviOhXBBms1FUnCJkvkQhnV/9SgjTQ91+KCpOETJfqpDCbEi42w9FxSlChpBgrzinCK09\nd37CkYSUN0KqWMU5RajhxOMSJktzT/wIKR0hVSxOETKEBHut9BShD3cI9TpyzzAbQkjI1EpP\nEVoiDz+Rv31znb2egZCQqZWeIrRE1oVYegohwVIrPUWIkFBatv8dV+M7y7blXoKQ8kdIFavg\nkC77lfNl6xXOQX27r63JtSAh5Y+QKlbBIckk58t86XbsmWNl1405FiSk/BFSxbIK6fWq/T91\nZu+Uy3IsSEj5I6SKZRXSzfK8Nz9+dI4FCSl/hFSxrEJaKJ978ws65liQkPJHSBXLKqS7ZKk3\nf9SAHAsSUv4IqWIVHlLfK+57udfMJmf2hTbH5liQkPJHSBWr4JAGVHlnqT1izLnta1/LsSAh\n5Y+QKlbhb8huXPrA1ad/9Qlj6vb6W67lCCl/hFSxFD5o7L3cPyak/BFSxWqln9hHSCgtQjKE\nBHuEZAgJ9gjJEBLsEZIhJNgjJENIsEdIhpBgj5AMIcEeIRlCgj1CMoQEe4RkCAn2CMkQEuwR\nkiEk2CMkQ0iwR0iGkGCPkAwhwV7lhNTj0hCfL3ETIaGkKiekmlCfeCSrQ6yakGCrckKqPj/E\nwr+Rj0MsTUiwRUiGkGCPkAwhwR4hGUKCPUIyhAR7hGQICfYIyRAS7BGSISTYIyRDSLBHSIaQ\nYI+QDCHBHiEZQoI9QjKEBHuEZAgJ9gjJEBLsEZIhJNgjJENIsEdIhpBgj5AMIcEeIRlCgj1C\nMoQEe4RkCAn2CMkQEuwRkiEk2CMkQ0iwR0iGkGCPkAwhwR4hGUKCPUIyhAR7hGQICfYIyRAS\n7BGSISTYIyRDSLBHSIaQYI+QDCHBHiEZQoI9QjKEBHuEZAgJ9gjJEBLsEZIhJNgjJENIsEdI\nhpBgj5AMIcEeIRlCgj1CMoQEe4RkCAn2CMkQEuwRkiEk2CMkQ0iwR0iGkGCPkAwhwR4hGUKC\nPUIyhAR7hGQICfYIyRAS7BGSISTYIyRDSLBHSIaQYI+QDCHBHiEZQoI9QjKEBHuEZAgJ9gjJ\nEBLsEZIhJNgjJENIsEdIhpBgj5AMIcEeIZmKDenEI14J4bXGMOtGSIRkKjak3SSUh8KsGyER\nkqnYkIYNWRNCr3vCrBshEZKp3JCGhVm6NyEVEyEZQoI9QjKEBHuEZAgJ9gjJEBLsEZIhJNgj\nJENIsEdIhpBgj5AMIcEeIRlCgj1CMoQEe4RkCAn2CMkQEuwRkiEk2CMkQ0iwR0iGkGCPkAwh\nwR4hGUKCPUIyhAR7hGQICfYIyRAS7BGSISTYIyRDSLBHSIaQYI+QDCHBHiEZQoI9QjKEBHuE\nZAgJ9gjJEBLsEZIhJNgjJENIsEdIhpBgj5AMIcEeIRlCgj1CMoQEe4RkCAn2CMkQEuwRkiEk\n2CMkQ0iwR0iGkGCPkAwhwR4hGUKCPUIyhAR7hGQICfYIyRAS7BGSISTYIyRDSLBHSIaQYI+Q\nDCHBHiEZQoI9QjKEBHuEZCISUu2QUSFM3hZm3SAkVyRCqh5/df7OlXVh1g1CckUjpDC33xJC\nComQDCFlIqSwCMkQUiZCCouQDCFlIqSwCMkQUiZCCouQDCFlIqSwCMkQUiZCCouQDCFlIqSw\nCMkQUiZCCouQDCFlIqSwCMkQUiZCCouQDCFlIqSwCMkQUqYl8vATIXwQZkNaJ0IyhJTp9xLK\nAWE2pHUiJENImcLdfleOD7MhrRMhGULKREhhEZIhpEyEFBYhGULKREhhEZIhpEyEFBYhGULK\nREhhEZIhpEwhQxqzJoTPwmx1xSAkQ0iZwt1+M8K96/SzMJtdKQjJEFKmkLdfx1dCGHlFmM2u\nFHYhNTWsbGxpGULKX6WGFOr2m0BIqRaf2K+tSE3/2YtzLkZI+SOkilVwSJumivQdM21a/QCR\n6ZtzLEhI+YtESLuNnh/C3NrIiU4AAA42SURBVGfDPG1cEWZDVBUc0kKZ+lpsbtkcyfU7hpDy\nF4mQOvc+Ln9HhHsdY2CYDVFVcEj1wxP/zXrTxHE5FiSk/EUjpHC33y9CvLL+8673hvFhmM1u\nQcEhdT41OX9J+g35j17dEnaUrc2s4vR23UKQ9iEW7ihdQizdtirMhlS1DbFwF+kYYun2EmZD\nqqvDLB2F26823N+v+YU++LMoOKSxu29PzE8em/bDxqeT/+jr8bubW8XKMP927IlfPhJi4cd/\nHmbVD94ZZuk7Hwyz9M8fD7HwI78Ms+rf/z7M0tx+GVYW+uDPouCQFsnhb8Tmlp8sl2ttDlCZ\nCn/VbppzaDdhxpETh4gclutVOyACLN5HenpOnxqRmj6zntLbHKAy2Z3Z0Lh6VYtnNgARUPxz\n7YAIICRAASEBCggJUEBIgAJCAhQQEqCAkAAFhAQoICRAASEBCggJUEBIgAJCAhQQEqCAkAAF\nhAQoKGdI9eH+8yQgoWpDGR+42ZQzpBOOCPO/0Wq4qnupR3xG7i71kF85u9Qjnv2VUo94t6wr\n4wM3m3KGFO5/WtVwT+9Sj7hOlpR6yPFXlnrE0n9i3xJCCiCkoiCkciCk4iKkoiCkIEIqCkIq\nB0IqLkIqCkIKIqSiIKRyIKTiIqSiIKQgQioKQioHQiouQioKQgoipKIgpHIoZ0jzNT96MC9/\nGFDqETdUv1HqIQ+8ptQjXnNgqUd8o5pz7ZLWrCn1iNs0P343P++XfMRVJX+MbVhV6hHLcLO2\ngH9GASggJEABIQEKCAlQQEiAAkICFBASoICQAAWEBCggJEABIQEKCAlQQEiAAkICFBASoICQ\nAAXlC2nzf47rPO7yzcUboH/scwsuSx0r+6yGm7rGpi0OpjZufMRS7ena/9ir4y4nf9D8ytWH\nDIxY6rszpPKFdLgMP2VXOaxo699Y1XeS61epY2WfVbBhD/9h3eJgWuPGRyzVnm4YIvULDqnq\n8Eo+46gMGRix1HdnWGUL6Wk5fLvZdqgsLtYAS+WKLGNln7X32A+HS9f8BlMaNzliqfZ0oXzT\n+fqn6r1KtpOBEUt7d4ZXtpDmiPu/grwqJxVrgAfkvixjZZ+1V+s86+ia32BK4yZHLNWejm3v\n/W8QU+RfpdrJwIilvTvDK1tIfWP/oU/ffsUa4Gp5+e7v/eKNtLGyz9rbvHmz/0SrxcGUxk2O\nWKo9/cqh3mSavF2qnQyMWNq7M7xyhdRYM8GbjmnbVKQRTpde7kckLtgaHCv7rM6AI7yHdYuD\nKY4bG7HEe7q6/U7bSrmTsRFLf3eGVK6QVssMbzpNGoo0wldl5uvrn91PrgyOlX1WZ8DYw7rF\nwRTH9UMq6Z6+PURuLelOxkYs/d0ZUrlCWiVHetNpsrJIIzzxiPur6dNuHRsDY2Wf1Rkw9rBu\ncTDFcf2QSrin677Tod2PS7qT/oilvztDKt9Tu4netL6msbgDHSvvBMbKPqszUvypXQuDKY7r\nh+QrwZ4+3FemLTOl3Mn4iHGluztDKtuLDX2GeJOB/Ys8zpmyLDhW9lkV/sO6xcH0xk0Nqfh7\neqkM8V9cLtVOJkf0le7uDKlsIc2Sd52v/yuzi7T+d+vO9abj2m0LjpV9VoX/sG5xML1xYyOW\nbE/vkKPi/3V9iXYyOWLp786QyhbSU3KK8/WE4r1/tlfti87XX8vclLGyz6rwQ2pxML1x/RFL\ntKdNu+34WXy+NDsZHLHkd2dIZQupaaocdOkkObxoA7xU2+bos74qu65JGSv7rAr/Yd3iYHrj\n+iOWaE8/kB5TYj4t0U4GRyz53RlS+c612/S9+s71xTzH8NVj+u+w7yUb08bKPqshfsTS4mBq\n48ZHLM2ePiVxH5doJ1NGLPXdGRL/jAJQQEiAAkICFBASoICQAAWEBCggJEABIQEKCAlQQEiA\nAkICFBASoICQAAWEBCggJEABIQEKCAlQQEiAAkICFBASoICQAAWEBCggJEABIQEKCAlQQEiA\nAkICFBASoICQAAWEBCggJEABIQEKCAlQQEj5WCTyV3/2OyIvqa33JNnufP377N07DD743qYw\n1zzzwsJGLPR6WUwIfn742l5/D0wiiJDy4YR0nj+7a3pID8ndyQsj5ckw6/VC+ll120Pnz+gu\n00OU9FznT8OMk3q9vDYyZbeySgnJXDWmMTCJHkLKxyLZsV/sYf6G7Kgc0lvt+rmfa7/+KLkh\n/yuOKfAPi3e9ooT0RaffBibRQ0j5WCRz5PnYXNXs5kLa8swnpmH1ljDrdUP6idzizf+76tC8\nr/c3eT2/Bd1tyrheXhsZNiRzSn1wEjmElI9Fcm/NBd7cXvUXeSG9dtyA9gNmOg/LKe6HbjeY\neXVrR8jNZp588Wa7Sc7Pt47ovso8KSOd2dhX03DG8E4jb3aPidaes3enfb/pfj63G9J58mBs\nkKt+aMz0Tu7cZjnJmAVdN583vNfRq784c9iOBy51v33XmK49DnjMnTttd+8a266s7zT43FXG\nBMasm/e/Jw0YMNP9I+dvU8b1nI0087puXzSww4hbnYuBaySH93crKLkDbx3Tr99xH3ghJdf9\nmLwamEQOIeVjkSyePMB9brdcrvNCerdLu6PP3K+m+0rz2Hky/47NZt5OM3oe87z3GF0ktxvz\nffcXejCkDwbWTJ0/VC42ZuVgmTB/P9nz81hIP5V9/icxUDCkTlMXvnBt9cjRFz5/Q9uhzpHH\nD6T3nFM617ivevRe4C62ZYKMnj9BBq0wgTHrJnU95d5renVx1ulvU8b1vJC6nDLriUf2k3tT\nrpEc3t+tgOQOPNup+sBTBtYN6p+yTRtqrjLJSeQQUj6ckG7y/hBdJf/wQvqu/NG59FO5K/4c\naJ6M+tTEHqNbRvT49L3aI5xLW1Y3JL6e4F5jS33VSjPfOxi6TC6PhfT5riLjr3ktdggWDEku\ndWZnyNnO15PlfWN22mWTMc/L6cYslTvcxf5LFjlffymzTGDMOvmG8/Xd2smJbcq4nheSTHfm\nVsjslGsEhs94apfYgaZ9qx92DoYmS//gup1DrykmMIkaQsqHE9Lqavc4fdQo44X01K3ui1OP\nyfXJkB42JvYYNS9VnzSl28rUNfyr6hB38vDov2xtu7c7t6VvH/9Vu88uGew8jeo1d4lJC+lN\nZ/Ziec64f2teMVtrhjkLN72y3Hk+JS+6iw0Y6r1CNqH9psCYdW3+5X7zDHnH36bM68VCesqd\n7TYl5Ro5QkruwMtyvDu31AkpsG5jju9hApOoIaR8OCGZAwYZ86FcFQvJ8cUL1+4RDOlDY/yQ\nzDdEfp22hmflSn9uuZzrTY9xFo29j+T8Obh5Zj+puSUtJPfPybfFPXa52gnJTJXh1yzxlr9W\n3Efu5zL2btfB8kZgzLqh3hp/IX+Kb1P69fyQvFchek5JuUaOkJI7cLfc5k3r+gfXbcw5VVsD\nk6ghpHy4If1UXjY/cn5reyF9du7w6qrhhwRDWm9MPKR3peO6tDX8Rm715xb7j8ivO+tyQ1q9\nJvb9Z7vXfuQ/kjfFQnKfEX5b3jN+SGvP6ybS4+sN7pvCq5xvLpO4FwJj1sVeNXtQbolvU/r1\n/JDcDY2FlLxGYPj0kJI7cJ38tzcd1T+4bmMulf8LTKKGkPLhhrSy6ptm3F4mFtJ0OeH+teal\nYEjeAzM2mdHeO7IJelJ+6M/F/yLNlHVeSF338n/wbbnPfyR/lDUkY7Yt/s/dZd8m/y/Lv2OH\nJr74mHWDvIs3OEc0/jalXy8jpOQ1AsOnh5Tcgd+7L2w4+vUPrpu/SGiZG5KZsPPHVYtiIa1v\nO8P99gPZQ7pbrp9V9XzqGv4pR7mTvw6+c2sb79Xwrf17x46RJrZ9L7bIfPmrmd7OfUg+lC2k\n969+1l3qAFlhfh071ukRS/C+m0xgzLpqb3XT5G1/mzKvlx5S8hqB4dNDSu7AEjnBnVte1T+4\nbmPmxA6O5nCMhOZ4IV0vc93jfzekT+Vw57vr93N/Sz/kHTIEQlrdY/T2lZ332OLE0vCZiX89\ntOpx57B8hrOGee57O2ahLIyFdIeMdF9VMH+u7bvBnOoO9PmobCEtl7HOwcj2se02mdflTvca\nl3gv/73U5lgTGLNOjnC+3iMHxLcp83rpISWvERj+If9IKCG5A/tV/9l5/jdN+gfXbcy+k73l\n/EnUEFI+vJA+qpLdjP/UbpJM+e7ZvQ5u2/8W85iMumJDMKRjal415kb31e3g+0hv9mwz49x9\n3Kd1nwyUyWfXJ95HMmdIzb4zZ+wmnV50H76dL7h4l68OzBJS06Ey8vwTBshFzsN5p7Pcwdbv\nKV89b1b7ug9MYMy6/t13nXdQVZdX4iFlXi89pOQ1AsPHdisguQPPdao+5PRhnSb0D67bbGzj\nHfv5k8ghpHx4IZmxconxQ/q/r/XrMulWc+fYRWbTEbXd1wRCukfc18kb92//v6lnNnx80rBO\n+97ivmD92VkjOu7zrfiZDcY8c9TwDn3GXuodsd8xorbuvC/OuC41pGsHLHOu9p1dd+hRf5u7\nhlP28DZr4zf32WHomR+5f04SY9ZNevfIun7HvmsSfyUzrhcIafBcE7xGYPjYbgUld+Cto/vX\nHfPqz05NWfeT4r5+H59EDiFVohdlaXM/qpsU9no5r5G/ufsHJ5FDSBVp/4ua+0nuLLJdTyek\nDZ3vDkyih5Aq0rNdGpr5Se4ssl1PJ6Qf7tcYmEQPIVWmZv+lawtZZLle89e4TpKa/RMYs3an\n1wKTCCIkQAEhAQoICVBASIACQgIUEBKggJAABYQEKCAkQAEhAQoICVBASIACQgIUEBKggJAA\nBYQEKCAkQAEhAQoICVBASIACQgIUEBKggJAABYQEKCAkQAEhAQoICVDw/wFFBEf7rVHpLAAA\nAABJRU5ErkJggg==",
      "text/plain": [
       "Plot with title \"Histogram of Matrix::colSums(exprs(input_cds))\""
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 420
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Visualize peak_count_per_cell\n",
    "hist(Matrix::colSums(exprs(input_cds)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter cells by peak_count\n",
    "# Please set an appropriate threshold values according to your data \n",
    "max_count <-  15000\n",
    "min_count <- 2000\n",
    "input_cds <- input_cds[,Matrix::colSums(exprs(input_cds)) >= min_count] \n",
    "input_cds <- input_cds[,Matrix::colSums(exprs(input_cds)) <= max_count] \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Process Cicero-CDS object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Overlap QC metrics:\n",
      "Cells per bin: 50\n",
      "Maximum shared cells bin-bin: 44\n",
      "Mean shared cells bin-bin: 0.84960828849071\n",
      "Median shared cells bin-bin: 0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Data preprocessing\n",
    "set.seed(2017)\n",
    "\n",
    "input_cds <- detect_genes(input_cds)\n",
    "input_cds <- estimate_size_factors(input_cds)\n",
    "input_cds <- preprocess_cds(input_cds, method = \"LSI\")\n",
    "\n",
    "# Dimensional reduction with umap\n",
    "input_cds <- reduce_dimension(input_cds, reduction_method = 'UMAP', \n",
    "                              preprocess_method = \"LSI\")\n",
    "umap_coords <- reducedDims(input_cds)$UMAP\n",
    "\n",
    "\n",
    "cicero_cds <- make_cicero_cds(input_cds, reduced_coordinates = umap_coords)\n",
    "\n",
    "# Save Cds object (Optional)\n",
    "#saveRDS(cicero_cds, paste0(output_folder, \"/cicero_cds.Rds\"))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Load reference genome information\n",
    "\n",
    "To run Cicero, you need to get a genomic coordinate file that contains the length of each chromosome.\n",
    "You can download the mm10 genomic information with the following command.\n",
    "\n",
    "If your scATAC-seq data was generated with a different reference genome, you will need to get the genome coordinates file for the reference genome you used. See the Cicero documentation for more information.\n",
    "\n",
    "https://cole-trapnell-lab.github.io/cicero-release/docs_m3/#installing-cicero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !!Please make sure that the reference genome information below matches your scATAC-seq reference genome.\n",
    "\n",
    "# If your scATAC-seq was aligned to the mm10 reference genome, you can read the chromosome length file using the following command.\n",
    "download.file(url = \"https://raw.githubusercontent.com/morris-lab/CellOracle/master/docs/demo_data/mm10_chromosome_length.txt\",\n",
    "              destfile = \"./mm10_chromosome_length.txt\")\n",
    "chromosome_length <- read.table(\"./mm10_chromosome_length.txt\")\n",
    "\n",
    "# For mm9 genome, you can use the following command.\n",
    "#data(\"mouse.mm9.genome\")\n",
    "#chromosome_length <- mouse.mm9.genome\n",
    "\n",
    "# For hg19 genome, you can use the following command.\n",
    "#data(\"human.hg19.genome\")\n",
    "#chromosome_length <- mhuman.hg19.genome\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Run Cicero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A data.frame: 6 × 3</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>Peak1</th><th scope=col>Peak2</th><th scope=col>coaccess</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>1</th><td>chr10_100006139_100006389</td><td>chr10_99774288_99774570  </td><td>-0.003546179</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>chr10_100006139_100006389</td><td>chr10_99825945_99826237  </td><td>-0.027536333</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>chr10_100006139_100006389</td><td>chr10_99830012_99830311  </td><td> 0.009588013</td></tr>\n",
       "\t<tr><th scope=row>4</th><td>chr10_100006139_100006389</td><td>chr10_99833211_99833540  </td><td>-0.008067111</td></tr>\n",
       "\t<tr><th scope=row>5</th><td>chr10_100006139_100006389</td><td>chr10_99941805_99941955  </td><td> 0.000000000</td></tr>\n",
       "\t<tr><th scope=row>7</th><td>chr10_100006139_100006389</td><td>chr10_100015291_100017830</td><td>-0.015018099</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 6 × 3\n",
       "\\begin{tabular}{r|lll}\n",
       "  & Peak1 & Peak2 & coaccess\\\\\n",
       "  & <chr> & <fct> & <dbl>\\\\\n",
       "\\hline\n",
       "\t1 & chr10\\_100006139\\_100006389 & chr10\\_99774288\\_99774570   & -0.003546179\\\\\n",
       "\t2 & chr10\\_100006139\\_100006389 & chr10\\_99825945\\_99826237   & -0.027536333\\\\\n",
       "\t3 & chr10\\_100006139\\_100006389 & chr10\\_99830012\\_99830311   &  0.009588013\\\\\n",
       "\t4 & chr10\\_100006139\\_100006389 & chr10\\_99833211\\_99833540   & -0.008067111\\\\\n",
       "\t5 & chr10\\_100006139\\_100006389 & chr10\\_99941805\\_99941955   &  0.000000000\\\\\n",
       "\t7 & chr10\\_100006139\\_100006389 & chr10\\_100015291\\_100017830 & -0.015018099\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 6 × 3\n",
       "\n",
       "| <!--/--> | Peak1 &lt;chr&gt; | Peak2 &lt;fct&gt; | coaccess &lt;dbl&gt; |\n",
       "|---|---|---|---|\n",
       "| 1 | chr10_100006139_100006389 | chr10_99774288_99774570   | -0.003546179 |\n",
       "| 2 | chr10_100006139_100006389 | chr10_99825945_99826237   | -0.027536333 |\n",
       "| 3 | chr10_100006139_100006389 | chr10_99830012_99830311   |  0.009588013 |\n",
       "| 4 | chr10_100006139_100006389 | chr10_99833211_99833540   | -0.008067111 |\n",
       "| 5 | chr10_100006139_100006389 | chr10_99941805_99941955   |  0.000000000 |\n",
       "| 7 | chr10_100006139_100006389 | chr10_100015291_100017830 | -0.015018099 |\n",
       "\n"
      ],
      "text/plain": [
       "  Peak1                     Peak2                     coaccess    \n",
       "1 chr10_100006139_100006389 chr10_99774288_99774570   -0.003546179\n",
       "2 chr10_100006139_100006389 chr10_99825945_99826237   -0.027536333\n",
       "3 chr10_100006139_100006389 chr10_99830012_99830311    0.009588013\n",
       "4 chr10_100006139_100006389 chr10_99833211_99833540   -0.008067111\n",
       "5 chr10_100006139_100006389 chr10_99941805_99941955    0.000000000\n",
       "7 chr10_100006139_100006389 chr10_100015291_100017830 -0.015018099"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Run the main function\n",
    "conns <- run_cicero(cicero_cds, chromosome_length) # Takes a few minutes to run\n",
    "\n",
    "# Save results (Optional)\n",
    "#saveRDS(conns, paste0(output_folder, \"/cicero_connections.Rds\"))\n",
    "\n",
    "# Check results\n",
    "head(conns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. Save results for the next step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_peaks <- row.names(exprs(input_cds))\n",
    "write.csv(x = all_peaks, file = paste0(output_folder, \"/all_peaks.csv\"))\n",
    "write.csv(x = conns, file = paste0(output_folder, \"/cicero_connections.csv\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Please go to next step: TSS annotation**\n",
    "\n",
    "https://morris-lab.github.io/CellOracle.documentation/tutorials/base_grn.html#step2-tss-annotation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "finalized": {
   "timestamp": 1642784300853,
   "trusted": true
  },
  "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}