dge.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "766704fd-538b-45de-982c-714e0d6d952b",
   "metadata": {},
   "source": [
    "# Differential gene expression for samples from BeatAML study*\n",
    "\n",
    "The analysis is based on two datasets from a study by Tyner et al. (2018). The first one can be obtained from http://www.vizome.org/aml (_Sample data_ -> _All sample data_). It contains rich data on tumor samples from patients suffering from AML (Acute Myeloid Leukemia). See the original paper as well as _Data dictionary_ on the provided website. The dataset is already saved in the input directory as _`BeatAML_RNASeq_rawcounts_2018_10_24.csv.gz`_.\n",
    "\n",
    "The second dataset can be obtained from http://vizome.org/additional_figures_BeatAML.html as _Full Cohort RNASeq Counts_. It contains raw read counts for samples subjected to RNA-Seq analysis. The dataset is already saved in the input directory as _`BeatAML_all_sample_data.csv.gz`_.\n",
    "\n",
    "Utilise those two datasets by taking subsequent steps in order to perform differential gene expression (DGE) between two groups: samples from patients with or without a mutation in NPM1 gene. The logic behind the analysis is concisely explained step by step.\n",
    "\n",
    "---\n",
    "\n",
    "<font size=\"1\">*Tyner JW, Tognon CE, Bottomly D, Wilmot B, Kurtz SE, Savage SL, Long N, Schultz AR, Traer E, Abel M, Agarwal A. Functional genomic landscape of acute myeloid leukaemia. Nature. 2018 Oct 25;562(7728):526-31.</font>\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc834f29-bcc5-4767-8b29-fb5f2f3bd07c",
   "metadata": {},
   "source": [
    "- Import `DESeq2` library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f104d5d3-5ee1-4f0f-a965-6cb06a86e695",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: S4Vectors\n",
      "\n",
      "Loading required package: stats4\n",
      "\n",
      "Loading required package: BiocGenerics\n",
      "\n",
      "\n",
      "Attaching package: ‘BiocGenerics’\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:stats’:\n",
      "\n",
      "    IQR, mad, sd, var, xtabs\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:base’:\n",
      "\n",
      "    anyDuplicated, aperm, append, as.data.frame, basename, cbind,\n",
      "    colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,\n",
      "    get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,\n",
      "    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,\n",
      "    Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,\n",
      "    table, tapply, union, unique, unsplit, which.max, which.min\n",
      "\n",
      "\n",
      "\n",
      "Attaching package: ‘S4Vectors’\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:base’:\n",
      "\n",
      "    expand.grid, I, unname\n",
      "\n",
      "\n",
      "Loading required package: IRanges\n",
      "\n",
      "Loading required package: GenomicRanges\n",
      "\n",
      "Loading required package: GenomeInfoDb\n",
      "\n",
      "Loading required package: SummarizedExperiment\n",
      "\n",
      "Loading required package: MatrixGenerics\n",
      "\n",
      "Loading required package: matrixStats\n",
      "\n",
      "\n",
      "Attaching package: ‘MatrixGenerics’\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:matrixStats’:\n",
      "\n",
      "    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,\n",
      "    colCounts, colCummaxs, colCummins, colCumprods, colCumsums,\n",
      "    colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,\n",
      "    colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,\n",
      "    colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,\n",
      "    colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,\n",
      "    colWeightedMeans, colWeightedMedians, colWeightedSds,\n",
      "    colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,\n",
      "    rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,\n",
      "    rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,\n",
      "    rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,\n",
      "    rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,\n",
      "    rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,\n",
      "    rowWeightedMads, rowWeightedMeans, rowWeightedMedians,\n",
      "    rowWeightedSds, rowWeightedVars\n",
      "\n",
      "\n",
      "Loading required package: Biobase\n",
      "\n",
      "Welcome to Bioconductor\n",
      "\n",
      "    Vignettes contain introductory material; view with\n",
      "    'browseVignettes()'. To cite Bioconductor, see\n",
      "    'citation(\"Biobase\")', and for packages 'citation(\"pkgname\")'.\n",
      "\n",
      "\n",
      "\n",
      "Attaching package: ‘Biobase’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:MatrixGenerics’:\n",
      "\n",
      "    rowMedians\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:matrixStats’:\n",
      "\n",
      "    anyMissing, rowMedians\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(\"DESeq2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af2ca8fd-1597-47e7-9c43-62a39e3fafa7",
   "metadata": {},
   "source": [
    "- Load BeatAML sample data on all samples to `smplData` data frame. Select only those rows for which `RNASeq` column values are `TRUE`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3d9b9575-75da-40a5-ac13-89f15bea7e93",
   "metadata": {},
   "outputs": [],
   "source": [
    "smplData <- read.csv(\"input/BeatAML_all_sample_data.csv.gz\", header=TRUE, check.names=FALSE, skip=3)\n",
    "smplData <- smplData[smplData[, \"RNAseq\"],]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc5918f6-abc8-439b-a401-4940626e30a3",
   "metadata": {},
   "source": [
    "- Get sample IDs (`LLS_SampleID` column) and check for which the column `SpecificDxAtInclusion` (diagnosis at inclusion) indicates samples from patients with AML diagnosis with mutated NPM1 gene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ae8bcbdb-2035-410c-a05b-16b754702341",
   "metadata": {},
   "outputs": [],
   "source": [
    "smpls   <- smplData[, \"LLS_SampleID\"]\n",
    "mutated <- smplData[, \"SpecificDxAtInclusion\"] == \"AML with mutated NPM1\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3522483-ccaf-4184-a3e2-44f60b30bef6",
   "metadata": {},
   "source": [
    "- Create a data frame `metaData` describing samples for DGE. Row names denote sample IDs, the `groups` column contains the following labels: `NPM1norm` for samples without the mutation diagnosed, `NPM1mut` for samples with the mutation diagnosed. The type of this column is set to `factor` and `NPM1norm` as the reference group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "db11c9ff-4a0d-4273-8fba-68b95fd622c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "groups   <- rep(\"NPM1norm\", length(smpls))\n",
    "groups[mutated] <- \"NPM1mut\"\n",
    "groups   <- factor(groups)\n",
    "groups   <- relevel(groups, ref=\"NPM1norm\")\n",
    "metaData <- data.frame(row.names=smpls, groups)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20fa8e66-756b-4bf8-8849-6b775db7b17b",
   "metadata": {},
   "source": [
    "- Look up the dimensions and first 5 rows of the resulting metaData data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "110e76b2-8267-4e93-ba53-83ec83bb884e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>507</li><li>1</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 507\n",
       "\\item 1\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 507\n",
       "2. 1\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 507   1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 6 × 1</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>groups</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;fct&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>15-00246</th><td>NPM1mut </td></tr>\n",
       "\t<tr><th scope=row>17-00040</th><td>NPM1norm</td></tr>\n",
       "\t<tr><th scope=row>14-00514</th><td>NPM1norm</td></tr>\n",
       "\t<tr><th scope=row>15-00981</th><td>NPM1norm</td></tr>\n",
       "\t<tr><th scope=row>13-00150</th><td>NPM1norm</td></tr>\n",
       "\t<tr><th scope=row>14-00376</th><td>NPM1norm</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 6 × 1\n",
       "\\begin{tabular}{r|l}\n",
       "  & groups\\\\\n",
       "  & <fct>\\\\\n",
       "\\hline\n",
       "\t15-00246 & NPM1mut \\\\\n",
       "\t17-00040 & NPM1norm\\\\\n",
       "\t14-00514 & NPM1norm\\\\\n",
       "\t15-00981 & NPM1norm\\\\\n",
       "\t13-00150 & NPM1norm\\\\\n",
       "\t14-00376 & NPM1norm\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 6 × 1\n",
       "\n",
       "| <!--/--> | groups &lt;fct&gt; |\n",
       "|---|---|\n",
       "| 15-00246 | NPM1mut  |\n",
       "| 17-00040 | NPM1norm |\n",
       "| 14-00514 | NPM1norm |\n",
       "| 15-00981 | NPM1norm |\n",
       "| 13-00150 | NPM1norm |\n",
       "| 14-00376 | NPM1norm |\n",
       "\n"
      ],
      "text/plain": [
       "         groups  \n",
       "15-00246 NPM1mut \n",
       "17-00040 NPM1norm\n",
       "14-00514 NPM1norm\n",
       "15-00981 NPM1norm\n",
       "13-00150 NPM1norm\n",
       "14-00376 NPM1norm"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dim(metaData)\n",
    "head(metaData)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce8a64ed-cf06-413d-a2e6-e2867eee6346",
   "metadata": {},
   "source": [
    "- Look the number of samples in each group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "508eb195-ceef-413f-9fbe-4b78331a0fe0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "NPM1norm  NPM1mut \n",
       "     400      107 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "table(metaData[, \"groups\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1c0189f-b159-4534-8077-1cc37d936437",
   "metadata": {},
   "source": [
    "- Load BeatAML raw gene counts to the countData data frame, look up the dimensions and the first 5 rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ea854300-af99-40dc-b7c3-46921c266249",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>63677</li><li>502</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 63677\n",
       "\\item 502\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 63677\n",
       "2. 502\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 63677   502"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 6 × 502</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>Symbol</th><th scope=col>Chr</th><th scope=col>Exon_Start</th><th scope=col>Exon_End</th><th scope=col>Strand</th><th scope=col>Length</th><th scope=col>GeneStart</th><th scope=col>GeneEnd</th><th scope=col>12-00023</th><th scope=col>12-00051</th><th scope=col>⋯</th><th scope=col>17-00043</th><th scope=col>17-00044</th><th scope=col>17-00045</th><th scope=col>17-00046</th><th scope=col>17-00047</th><th scope=col>17-00049</th><th scope=col>17-00051</th><th scope=col>17-00052</th><th scope=col>17-00053</th><th scope=col>17-00055</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>⋯</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>ENSG00000000003</th><td>TSPAN6  </td><td>X </td><td>99883667;99885756;99887482;99888402;99888928;99890175;99890555;99891188;99894942                                                                                                                                                                                                                                                         </td><td>99884983;99885863;99887565;99888536;99889026;99890249;99890743;99892101;99894988                                                                                                                                                                                                                                                         </td><td>-</td><td>2968</td><td> 99883667</td><td> 99894988</td><td>    8</td><td> 52</td><td>⋯</td><td>   41</td><td>   16</td><td>   18</td><td>   28</td><td>   35</td><td>   34</td><td>   35</td><td>   27</td><td>  34</td><td>   20</td></tr>\n",
       "\t<tr><th scope=row>ENSG00000000005</th><td>TNMD    </td><td>X </td><td>99839799;99840228;99848621;99849258;99852501;99854013;99854505                                                                                                                                                                                                                                                                           </td><td>99840063;99840359;99849032;99849359;99852654;99854179;99854882                                                                                                                                                                                                                                                                           </td><td>+</td><td>1610</td><td> 99839799</td><td> 99854882</td><td>    0</td><td>  0</td><td>⋯</td><td>    0</td><td>    1</td><td>    0</td><td>    1</td><td>    0</td><td>    0</td><td>    0</td><td>    0</td><td>   0</td><td>    0</td></tr>\n",
       "\t<tr><th scope=row>ENSG00000000419</th><td>DPM1    </td><td>20</td><td>49551404;49552685;49557402;49557642;49558568;49562274;49562384;49565166;49571723;49574900                                                                                                                                                                                                                                                </td><td>49551773;49552799;49557492;49557746;49558663;49562299;49562460;49565199;49571822;49575092                                                                                                                                                                                                                                                </td><td>-</td><td>1207</td><td> 49551404</td><td> 49575092</td><td> 1080</td><td>920</td><td>⋯</td><td> 1087</td><td> 1208</td><td> 1058</td><td> 1297</td><td> 1069</td><td> 1163</td><td> 1041</td><td> 1366</td><td>1389</td><td> 1198</td></tr>\n",
       "\t<tr><th scope=row>ENSG00000000457</th><td>SCYL3   </td><td>1 </td><td>169818772;169823411;169824937;169828182;169831754;169833510;169836037;169838069;169839396;169842837;169845119;169847775;169857817;169862929;169863148                                                                                                                                                                                    </td><td>169822913;169824105;169825098;169828353;169831938;169833649;169836114;169838269;169839498;169842893;169845232;169847960;169858031;169863093;169863408                                                                                                                                                                                    </td><td>-</td><td>6876</td><td>169818772</td><td>169863408</td><td>  701</td><td>587</td><td>⋯</td><td>  881</td><td>  990</td><td>  927</td><td> 1008</td><td>  706</td><td>  996</td><td>  791</td><td>  865</td><td> 651</td><td>  751</td></tr>\n",
       "\t<tr><th scope=row>ENSG00000000460</th><td>C1orf112</td><td>1 </td><td>169631245;169652610;169652897;169752952;169754018;169763871;169764181;169764550;169767998;169770024;169771762;169772310;169773216;169775145;169776932;169790820;169792549;169795234;169796192;169796852;169798404;169799400;169799837;169801558;169806068;169811558;169812815;169816831;169818641;169819406;169819595;169820958;169821931</td><td>169631761;169652766;169653073;169753069;169754054;169764046;169764354;169765124;169768099;169770112;169771866;169772450;169773521;169775229;169777070;169790900;169792613;169795258;169796340;169796981;169798955;169799482;169799888;169801647;169806253;169811680;169812910;169816944;169818745;169819486;169819707;169821077;169823221</td><td>+</td><td>6354</td><td>169631245</td><td>169823221</td><td>  249</td><td> 82</td><td>⋯</td><td>  560</td><td>  548</td><td>  511</td><td>  569</td><td>  459</td><td>  421</td><td>  625</td><td>  444</td><td> 667</td><td>  416</td></tr>\n",
       "\t<tr><th scope=row>ENSG00000000938</th><td>FGR     </td><td>1 </td><td>27938575;27939730;27940941;27941361;27941945;27942200;27943368;27943704;27948070;27949553;27949895;27951600;27952557;27952971;27961576                                                                                                                                                                                                   </td><td>27939633;27939861;27941094;27941437;27942124;27942355;27943517;27943807;27948168;27949655;27950573;27951662;27952751;27953080;27961788                                                                                                                                                                                                   </td><td>-</td><td>3474</td><td> 27938575</td><td> 27961788</td><td>29052</td><td>541</td><td>⋯</td><td>21732</td><td>26652</td><td>24034</td><td>18079</td><td>25757</td><td>21527</td><td>20409</td><td>27210</td><td>1604</td><td>22325</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 6 × 502\n",
       "\\begin{tabular}{r|lllllllllllllllllllll}\n",
       "  & Symbol & Chr & Exon\\_Start & Exon\\_End & Strand & Length & GeneStart & GeneEnd & 12-00023 & 12-00051 & ⋯ & 17-00043 & 17-00044 & 17-00045 & 17-00046 & 17-00047 & 17-00049 & 17-00051 & 17-00052 & 17-00053 & 17-00055\\\\\n",
       "  & <chr> & <chr> & <chr> & <chr> & <chr> & <int> & <int> & <int> & <int> & <int> & ⋯ & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int>\\\\\n",
       "\\hline\n",
       "\tENSG00000000003 & TSPAN6   & X  & 99883667;99885756;99887482;99888402;99888928;99890175;99890555;99891188;99894942                                                                                                                                                                                                                                                          & 99884983;99885863;99887565;99888536;99889026;99890249;99890743;99892101;99894988                                                                                                                                                                                                                                                          & - & 2968 &  99883667 &  99894988 &     8 &  52 & ⋯ &    41 &    16 &    18 &    28 &    35 &    34 &    35 &    27 &   34 &    20\\\\\n",
       "\tENSG00000000005 & TNMD     & X  & 99839799;99840228;99848621;99849258;99852501;99854013;99854505                                                                                                                                                                                                                                                                            & 99840063;99840359;99849032;99849359;99852654;99854179;99854882                                                                                                                                                                                                                                                                            & + & 1610 &  99839799 &  99854882 &     0 &   0 & ⋯ &     0 &     1 &     0 &     1 &     0 &     0 &     0 &     0 &    0 &     0\\\\\n",
       "\tENSG00000000419 & DPM1     & 20 & 49551404;49552685;49557402;49557642;49558568;49562274;49562384;49565166;49571723;49574900                                                                                                                                                                                                                                                 & 49551773;49552799;49557492;49557746;49558663;49562299;49562460;49565199;49571822;49575092                                                                                                                                                                                                                                                 & - & 1207 &  49551404 &  49575092 &  1080 & 920 & ⋯ &  1087 &  1208 &  1058 &  1297 &  1069 &  1163 &  1041 &  1366 & 1389 &  1198\\\\\n",
       "\tENSG00000000457 & SCYL3    & 1  & 169818772;169823411;169824937;169828182;169831754;169833510;169836037;169838069;169839396;169842837;169845119;169847775;169857817;169862929;169863148                                                                                                                                                                                     & 169822913;169824105;169825098;169828353;169831938;169833649;169836114;169838269;169839498;169842893;169845232;169847960;169858031;169863093;169863408                                                                                                                                                                                     & - & 6876 & 169818772 & 169863408 &   701 & 587 & ⋯ &   881 &   990 &   927 &  1008 &   706 &   996 &   791 &   865 &  651 &   751\\\\\n",
       "\tENSG00000000460 & C1orf112 & 1  & 169631245;169652610;169652897;169752952;169754018;169763871;169764181;169764550;169767998;169770024;169771762;169772310;169773216;169775145;169776932;169790820;169792549;169795234;169796192;169796852;169798404;169799400;169799837;169801558;169806068;169811558;169812815;169816831;169818641;169819406;169819595;169820958;169821931 & 169631761;169652766;169653073;169753069;169754054;169764046;169764354;169765124;169768099;169770112;169771866;169772450;169773521;169775229;169777070;169790900;169792613;169795258;169796340;169796981;169798955;169799482;169799888;169801647;169806253;169811680;169812910;169816944;169818745;169819486;169819707;169821077;169823221 & + & 6354 & 169631245 & 169823221 &   249 &  82 & ⋯ &   560 &   548 &   511 &   569 &   459 &   421 &   625 &   444 &  667 &   416\\\\\n",
       "\tENSG00000000938 & FGR      & 1  & 27938575;27939730;27940941;27941361;27941945;27942200;27943368;27943704;27948070;27949553;27949895;27951600;27952557;27952971;27961576                                                                                                                                                                                                    & 27939633;27939861;27941094;27941437;27942124;27942355;27943517;27943807;27948168;27949655;27950573;27951662;27952751;27953080;27961788                                                                                                                                                                                                    & - & 3474 &  27938575 &  27961788 & 29052 & 541 & ⋯ & 21732 & 26652 & 24034 & 18079 & 25757 & 21527 & 20409 & 27210 & 1604 & 22325\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 6 × 502\n",
       "\n",
       "| <!--/--> | Symbol &lt;chr&gt; | Chr &lt;chr&gt; | Exon_Start &lt;chr&gt; | Exon_End &lt;chr&gt; | Strand &lt;chr&gt; | Length &lt;int&gt; | GeneStart &lt;int&gt; | GeneEnd &lt;int&gt; | 12-00023 &lt;int&gt; | 12-00051 &lt;int&gt; | ⋯ ⋯ | 17-00043 &lt;int&gt; | 17-00044 &lt;int&gt; | 17-00045 &lt;int&gt; | 17-00046 &lt;int&gt; | 17-00047 &lt;int&gt; | 17-00049 &lt;int&gt; | 17-00051 &lt;int&gt; | 17-00052 &lt;int&gt; | 17-00053 &lt;int&gt; | 17-00055 &lt;int&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n",
       "| ENSG00000000003 | TSPAN6   | X  | 99883667;99885756;99887482;99888402;99888928;99890175;99890555;99891188;99894942                                                                                                                                                                                                                                                          | 99884983;99885863;99887565;99888536;99889026;99890249;99890743;99892101;99894988                                                                                                                                                                                                                                                          | - | 2968 |  99883667 |  99894988 |     8 |  52 | ⋯ |    41 |    16 |    18 |    28 |    35 |    34 |    35 |    27 |   34 |    20 |\n",
       "| ENSG00000000005 | TNMD     | X  | 99839799;99840228;99848621;99849258;99852501;99854013;99854505                                                                                                                                                                                                                                                                            | 99840063;99840359;99849032;99849359;99852654;99854179;99854882                                                                                                                                                                                                                                                                            | + | 1610 |  99839799 |  99854882 |     0 |   0 | ⋯ |     0 |     1 |     0 |     1 |     0 |     0 |     0 |     0 |    0 |     0 |\n",
       "| ENSG00000000419 | DPM1     | 20 | 49551404;49552685;49557402;49557642;49558568;49562274;49562384;49565166;49571723;49574900                                                                                                                                                                                                                                                 | 49551773;49552799;49557492;49557746;49558663;49562299;49562460;49565199;49571822;49575092                                                                                                                                                                                                                                                 | - | 1207 |  49551404 |  49575092 |  1080 | 920 | ⋯ |  1087 |  1208 |  1058 |  1297 |  1069 |  1163 |  1041 |  1366 | 1389 |  1198 |\n",
       "| ENSG00000000457 | SCYL3    | 1  | 169818772;169823411;169824937;169828182;169831754;169833510;169836037;169838069;169839396;169842837;169845119;169847775;169857817;169862929;169863148                                                                                                                                                                                     | 169822913;169824105;169825098;169828353;169831938;169833649;169836114;169838269;169839498;169842893;169845232;169847960;169858031;169863093;169863408                                                                                                                                                                                     | - | 6876 | 169818772 | 169863408 |   701 | 587 | ⋯ |   881 |   990 |   927 |  1008 |   706 |   996 |   791 |   865 |  651 |   751 |\n",
       "| ENSG00000000460 | C1orf112 | 1  | 169631245;169652610;169652897;169752952;169754018;169763871;169764181;169764550;169767998;169770024;169771762;169772310;169773216;169775145;169776932;169790820;169792549;169795234;169796192;169796852;169798404;169799400;169799837;169801558;169806068;169811558;169812815;169816831;169818641;169819406;169819595;169820958;169821931 | 169631761;169652766;169653073;169753069;169754054;169764046;169764354;169765124;169768099;169770112;169771866;169772450;169773521;169775229;169777070;169790900;169792613;169795258;169796340;169796981;169798955;169799482;169799888;169801647;169806253;169811680;169812910;169816944;169818745;169819486;169819707;169821077;169823221 | + | 6354 | 169631245 | 169823221 |   249 |  82 | ⋯ |   560 |   548 |   511 |   569 |   459 |   421 |   625 |   444 |  667 |   416 |\n",
       "| ENSG00000000938 | FGR      | 1  | 27938575;27939730;27940941;27941361;27941945;27942200;27943368;27943704;27948070;27949553;27949895;27951600;27952557;27952971;27961576                                                                                                                                                                                                    | 27939633;27939861;27941094;27941437;27942124;27942355;27943517;27943807;27948168;27949655;27950573;27951662;27952751;27953080;27961788                                                                                                                                                                                                    | - | 3474 |  27938575 |  27961788 | 29052 | 541 | ⋯ | 21732 | 26652 | 24034 | 18079 | 25757 | 21527 | 20409 | 27210 | 1604 | 22325 |\n",
       "\n"
      ],
      "text/plain": [
       "                Symbol   Chr\n",
       "ENSG00000000003 TSPAN6   X  \n",
       "ENSG00000000005 TNMD     X  \n",
       "ENSG00000000419 DPM1     20 \n",
       "ENSG00000000457 SCYL3    1  \n",
       "ENSG00000000460 C1orf112 1  \n",
       "ENSG00000000938 FGR      1  \n",
       "                Exon_Start                                                                                                                                                                                                                                                                                                                               \n",
       "ENSG00000000003 99883667;99885756;99887482;99888402;99888928;99890175;99890555;99891188;99894942                                                                                                                                                                                                                                                         \n",
       "ENSG00000000005 99839799;99840228;99848621;99849258;99852501;99854013;99854505                                                                                                                                                                                                                                                                           \n",
       "ENSG00000000419 49551404;49552685;49557402;49557642;49558568;49562274;49562384;49565166;49571723;49574900                                                                                                                                                                                                                                                \n",
       "ENSG00000000457 169818772;169823411;169824937;169828182;169831754;169833510;169836037;169838069;169839396;169842837;169845119;169847775;169857817;169862929;169863148                                                                                                                                                                                    \n",
       "ENSG00000000460 169631245;169652610;169652897;169752952;169754018;169763871;169764181;169764550;169767998;169770024;169771762;169772310;169773216;169775145;169776932;169790820;169792549;169795234;169796192;169796852;169798404;169799400;169799837;169801558;169806068;169811558;169812815;169816831;169818641;169819406;169819595;169820958;169821931\n",
       "ENSG00000000938 27938575;27939730;27940941;27941361;27941945;27942200;27943368;27943704;27948070;27949553;27949895;27951600;27952557;27952971;27961576                                                                                                                                                                                                   \n",
       "                Exon_End                                                                                                                                                                                                                                                                                                                                 \n",
       "ENSG00000000003 99884983;99885863;99887565;99888536;99889026;99890249;99890743;99892101;99894988                                                                                                                                                                                                                                                         \n",
       "ENSG00000000005 99840063;99840359;99849032;99849359;99852654;99854179;99854882                                                                                                                                                                                                                                                                           \n",
       "ENSG00000000419 49551773;49552799;49557492;49557746;49558663;49562299;49562460;49565199;49571822;49575092                                                                                                                                                                                                                                                \n",
       "ENSG00000000457 169822913;169824105;169825098;169828353;169831938;169833649;169836114;169838269;169839498;169842893;169845232;169847960;169858031;169863093;169863408                                                                                                                                                                                    \n",
       "ENSG00000000460 169631761;169652766;169653073;169753069;169754054;169764046;169764354;169765124;169768099;169770112;169771866;169772450;169773521;169775229;169777070;169790900;169792613;169795258;169796340;169796981;169798955;169799482;169799888;169801647;169806253;169811680;169812910;169816944;169818745;169819486;169819707;169821077;169823221\n",
       "ENSG00000000938 27939633;27939861;27941094;27941437;27942124;27942355;27943517;27943807;27948168;27949655;27950573;27951662;27952751;27953080;27961788                                                                                                                                                                                                   \n",
       "                Strand Length GeneStart GeneEnd   12-00023 12-00051 ⋯ 17-00043\n",
       "ENSG00000000003 -      2968    99883667  99894988     8     52      ⋯    41   \n",
       "ENSG00000000005 +      1610    99839799  99854882     0      0      ⋯     0   \n",
       "ENSG00000000419 -      1207    49551404  49575092  1080    920      ⋯  1087   \n",
       "ENSG00000000457 -      6876   169818772 169863408   701    587      ⋯   881   \n",
       "ENSG00000000460 +      6354   169631245 169823221   249     82      ⋯   560   \n",
       "ENSG00000000938 -      3474    27938575  27961788 29052    541      ⋯ 21732   \n",
       "                17-00044 17-00045 17-00046 17-00047 17-00049 17-00051 17-00052\n",
       "ENSG00000000003    16       18       28       35       34       35       27   \n",
       "ENSG00000000005     1        0        1        0        0        0        0   \n",
       "ENSG00000000419  1208     1058     1297     1069     1163     1041     1366   \n",
       "ENSG00000000457   990      927     1008      706      996      791      865   \n",
       "ENSG00000000460   548      511      569      459      421      625      444   \n",
       "ENSG00000000938 26652    24034    18079    25757    21527    20409    27210   \n",
       "                17-00053 17-00055\n",
       "ENSG00000000003   34        20   \n",
       "ENSG00000000005    0         0   \n",
       "ENSG00000000419 1389      1198   \n",
       "ENSG00000000457  651       751   \n",
       "ENSG00000000460  667       416   \n",
       "ENSG00000000938 1604     22325   "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "countData <- read.csv(\"input/BeatAML_RNASeq_rawcounts_2018_10_24.csv.gz\", header=TRUE, check.names=FALSE, row.names=1)\n",
    "dim(countData)\n",
    "head(countData)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2db26fdd-bad1-42f0-94e5-feecebb292fc",
   "metadata": {},
   "source": [
    "- Select columns the names of which match sample ID pattern, look up the dimensions and the first 5 rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "79bb7a09-1565-4544-b9d3-529cd9f0e9ce",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>63677</li><li>494</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 63677\n",
       "\\item 494\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 63677\n",
       "2. 494\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 63677   494"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 6 × 494</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>12-00023</th><th scope=col>12-00051</th><th scope=col>12-00066</th><th scope=col>12-00150</th><th scope=col>12-00211</th><th scope=col>12-00258</th><th scope=col>12-00294</th><th scope=col>12-00301</th><th scope=col>12-00372</th><th scope=col>12-00423</th><th scope=col>⋯</th><th scope=col>17-00043</th><th scope=col>17-00044</th><th scope=col>17-00045</th><th scope=col>17-00046</th><th scope=col>17-00047</th><th scope=col>17-00049</th><th scope=col>17-00051</th><th scope=col>17-00052</th><th scope=col>17-00053</th><th scope=col>17-00055</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>⋯</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>ENSG00000000003</th><td>    8</td><td> 52</td><td>    2</td><td>  60</td><td>   0</td><td>    0</td><td>   0</td><td>    1</td><td>  32</td><td>  18</td><td>⋯</td><td>   41</td><td>   16</td><td>   18</td><td>   28</td><td>   35</td><td>   34</td><td>   35</td><td>   27</td><td>  34</td><td>   20</td></tr>\n",
       "\t<tr><th scope=row>ENSG00000000005</th><td>    0</td><td>  0</td><td>    0</td><td>   0</td><td>   0</td><td>    0</td><td>   0</td><td>    0</td><td>   0</td><td>   0</td><td>⋯</td><td>    0</td><td>    1</td><td>    0</td><td>    1</td><td>    0</td><td>    0</td><td>    0</td><td>    0</td><td>   0</td><td>    0</td></tr>\n",
       "\t<tr><th scope=row>ENSG00000000419</th><td> 1080</td><td>920</td><td>  664</td><td>1128</td><td> 747</td><td>  921</td><td> 955</td><td>  324</td><td> 797</td><td>1019</td><td>⋯</td><td> 1087</td><td> 1208</td><td> 1058</td><td> 1297</td><td> 1069</td><td> 1163</td><td> 1041</td><td> 1366</td><td>1389</td><td> 1198</td></tr>\n",
       "\t<tr><th scope=row>ENSG00000000457</th><td>  701</td><td>587</td><td>  300</td><td> 624</td><td>1120</td><td>  491</td><td>1198</td><td>  279</td><td>1172</td><td> 760</td><td>⋯</td><td>  881</td><td>  990</td><td>  927</td><td> 1008</td><td>  706</td><td>  996</td><td>  791</td><td>  865</td><td> 651</td><td>  751</td></tr>\n",
       "\t<tr><th scope=row>ENSG00000000460</th><td>  249</td><td> 82</td><td>  317</td><td> 395</td><td> 373</td><td>  257</td><td> 577</td><td>  133</td><td> 386</td><td> 487</td><td>⋯</td><td>  560</td><td>  548</td><td>  511</td><td>  569</td><td>  459</td><td>  421</td><td>  625</td><td>  444</td><td> 667</td><td>  416</td></tr>\n",
       "\t<tr><th scope=row>ENSG00000000938</th><td>29052</td><td>541</td><td>25829</td><td>1720</td><td>3326</td><td>31518</td><td>8188</td><td>13510</td><td>5928</td><td>9186</td><td>⋯</td><td>21732</td><td>26652</td><td>24034</td><td>18079</td><td>25757</td><td>21527</td><td>20409</td><td>27210</td><td>1604</td><td>22325</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 6 × 494\n",
       "\\begin{tabular}{r|lllllllllllllllllllll}\n",
       "  & 12-00023 & 12-00051 & 12-00066 & 12-00150 & 12-00211 & 12-00258 & 12-00294 & 12-00301 & 12-00372 & 12-00423 & ⋯ & 17-00043 & 17-00044 & 17-00045 & 17-00046 & 17-00047 & 17-00049 & 17-00051 & 17-00052 & 17-00053 & 17-00055\\\\\n",
       "  & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int> & ⋯ & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int>\\\\\n",
       "\\hline\n",
       "\tENSG00000000003 &     8 &  52 &     2 &   60 &    0 &     0 &    0 &     1 &   32 &   18 & ⋯ &    41 &    16 &    18 &    28 &    35 &    34 &    35 &    27 &   34 &    20\\\\\n",
       "\tENSG00000000005 &     0 &   0 &     0 &    0 &    0 &     0 &    0 &     0 &    0 &    0 & ⋯ &     0 &     1 &     0 &     1 &     0 &     0 &     0 &     0 &    0 &     0\\\\\n",
       "\tENSG00000000419 &  1080 & 920 &   664 & 1128 &  747 &   921 &  955 &   324 &  797 & 1019 & ⋯ &  1087 &  1208 &  1058 &  1297 &  1069 &  1163 &  1041 &  1366 & 1389 &  1198\\\\\n",
       "\tENSG00000000457 &   701 & 587 &   300 &  624 & 1120 &   491 & 1198 &   279 & 1172 &  760 & ⋯ &   881 &   990 &   927 &  1008 &   706 &   996 &   791 &   865 &  651 &   751\\\\\n",
       "\tENSG00000000460 &   249 &  82 &   317 &  395 &  373 &   257 &  577 &   133 &  386 &  487 & ⋯ &   560 &   548 &   511 &   569 &   459 &   421 &   625 &   444 &  667 &   416\\\\\n",
       "\tENSG00000000938 & 29052 & 541 & 25829 & 1720 & 3326 & 31518 & 8188 & 13510 & 5928 & 9186 & ⋯ & 21732 & 26652 & 24034 & 18079 & 25757 & 21527 & 20409 & 27210 & 1604 & 22325\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 6 × 494\n",
       "\n",
       "| <!--/--> | 12-00023 &lt;int&gt; | 12-00051 &lt;int&gt; | 12-00066 &lt;int&gt; | 12-00150 &lt;int&gt; | 12-00211 &lt;int&gt; | 12-00258 &lt;int&gt; | 12-00294 &lt;int&gt; | 12-00301 &lt;int&gt; | 12-00372 &lt;int&gt; | 12-00423 &lt;int&gt; | ⋯ ⋯ | 17-00043 &lt;int&gt; | 17-00044 &lt;int&gt; | 17-00045 &lt;int&gt; | 17-00046 &lt;int&gt; | 17-00047 &lt;int&gt; | 17-00049 &lt;int&gt; | 17-00051 &lt;int&gt; | 17-00052 &lt;int&gt; | 17-00053 &lt;int&gt; | 17-00055 &lt;int&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n",
       "| ENSG00000000003 |     8 |  52 |     2 |   60 |    0 |     0 |    0 |     1 |   32 |   18 | ⋯ |    41 |    16 |    18 |    28 |    35 |    34 |    35 |    27 |   34 |    20 |\n",
       "| ENSG00000000005 |     0 |   0 |     0 |    0 |    0 |     0 |    0 |     0 |    0 |    0 | ⋯ |     0 |     1 |     0 |     1 |     0 |     0 |     0 |     0 |    0 |     0 |\n",
       "| ENSG00000000419 |  1080 | 920 |   664 | 1128 |  747 |   921 |  955 |   324 |  797 | 1019 | ⋯ |  1087 |  1208 |  1058 |  1297 |  1069 |  1163 |  1041 |  1366 | 1389 |  1198 |\n",
       "| ENSG00000000457 |   701 | 587 |   300 |  624 | 1120 |   491 | 1198 |   279 | 1172 |  760 | ⋯ |   881 |   990 |   927 |  1008 |   706 |   996 |   791 |   865 |  651 |   751 |\n",
       "| ENSG00000000460 |   249 |  82 |   317 |  395 |  373 |   257 |  577 |   133 |  386 |  487 | ⋯ |   560 |   548 |   511 |   569 |   459 |   421 |   625 |   444 |  667 |   416 |\n",
       "| ENSG00000000938 | 29052 | 541 | 25829 | 1720 | 3326 | 31518 | 8188 | 13510 | 5928 | 9186 | ⋯ | 21732 | 26652 | 24034 | 18079 | 25757 | 21527 | 20409 | 27210 | 1604 | 22325 |\n",
       "\n"
      ],
      "text/plain": [
       "                12-00023 12-00051 12-00066 12-00150 12-00211 12-00258 12-00294\n",
       "ENSG00000000003     8     52          2      60        0         0       0    \n",
       "ENSG00000000005     0      0          0       0        0         0       0    \n",
       "ENSG00000000419  1080    920        664    1128      747       921     955    \n",
       "ENSG00000000457   701    587        300     624     1120       491    1198    \n",
       "ENSG00000000460   249     82        317     395      373       257     577    \n",
       "ENSG00000000938 29052    541      25829    1720     3326     31518    8188    \n",
       "                12-00301 12-00372 12-00423 ⋯ 17-00043 17-00044 17-00045\n",
       "ENSG00000000003     1      32       18     ⋯    41       16       18   \n",
       "ENSG00000000005     0       0        0     ⋯     0        1        0   \n",
       "ENSG00000000419   324     797     1019     ⋯  1087     1208     1058   \n",
       "ENSG00000000457   279    1172      760     ⋯   881      990      927   \n",
       "ENSG00000000460   133     386      487     ⋯   560      548      511   \n",
       "ENSG00000000938 13510    5928     9186     ⋯ 21732    26652    24034   \n",
       "                17-00046 17-00047 17-00049 17-00051 17-00052 17-00053 17-00055\n",
       "ENSG00000000003    28       35       34       35       27      34        20   \n",
       "ENSG00000000005     1        0        0        0        0       0         0   \n",
       "ENSG00000000419  1297     1069     1163     1041     1366    1389      1198   \n",
       "ENSG00000000457  1008      706      996      791      865     651       751   \n",
       "ENSG00000000460   569      459      421      625      444     667       416   \n",
       "ENSG00000000938 18079    25757    21527    20409    27210    1604     22325   "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "countData <- countData[, grepl(\"^[0-9]+\\\\-[0-9]+$\", colnames(countData))]\n",
    "dim(countData)\n",
    "head(countData)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b14329b-f10e-4be9-811a-f4b7ce2b3e03",
   "metadata": {},
   "source": [
    "- Find the intersection of sample IDs from sample data (`metaData`) and count data (`countData`), look up the number of obtained sample IDs and the first 5 of them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3ebe3bd9-00ae-4322-afb8-fad1d3c72697",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "493"
      ],
      "text/latex": [
       "493"
      ],
      "text/markdown": [
       "493"
      ],
      "text/plain": [
       "[1] 493"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>'15-00246'</li><li>'17-00040'</li><li>'14-00514'</li><li>'15-00981'</li><li>'13-00150'</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item '15-00246'\n",
       "\\item '17-00040'\n",
       "\\item '14-00514'\n",
       "\\item '15-00981'\n",
       "\\item '13-00150'\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. '15-00246'\n",
       "2. '17-00040'\n",
       "3. '14-00514'\n",
       "4. '15-00981'\n",
       "5. '13-00150'\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] \"15-00246\" \"17-00040\" \"14-00514\" \"15-00981\" \"13-00150\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "common <- intersect(rownames(metaData), colnames(countData))\n",
    "length(common)\n",
    "common[1:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48cef808-23a1-49ae-a432-fbe716f1f1f4",
   "metadata": {},
   "source": [
    "- From the `metaData` and `countData` data frames select respectively rows and columns that correspond to sample IDs common for both data frames, look up dimensions of the resulting data frames."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "db331edd-4a4c-42dc-b802-bff8243d437a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>493</li><li>1</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 493\n",
       "\\item 1\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 493\n",
       "2. 1\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 493   1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>63677</li><li>493</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 63677\n",
       "\\item 493\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 63677\n",
       "2. 493\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 63677   493"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "metaData  <- metaData[common, , drop=FALSE]\n",
    "countData <- countData[, common, drop=FALSE]\n",
    "dim(metaData)\n",
    "dim(countData)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6491e8e-f642-4bff-81e6-ffd3b4778115",
   "metadata": {},
   "source": [
    "- Using `metaData` and `countData` data frames, create a DESEq2 dataset (`dss`). Get rid of rows with less than total 10 counts and perform DGE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "bc8680de-dd15-42c3-8e88-29f263bf3bae",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "estimating size factors\n",
      "\n",
      "estimating dispersions\n",
      "\n",
      "gene-wise dispersion estimates\n",
      "\n",
      "mean-dispersion relationship\n",
      "\n",
      "final dispersion estimates\n",
      "\n",
      "fitting model and testing\n",
      "\n",
      "-- replacing outliers and refitting for 4541 genes\n",
      "-- DESeq argument 'minReplicatesForReplace' = 7 \n",
      "-- original counts are preserved in counts(dds)\n",
      "\n",
      "estimating dispersions\n",
      "\n",
      "fitting model and testing\n",
      "\n"
     ]
    }
   ],
   "source": [
    "dds <- DESeqDataSetFromMatrix(countData=countData, colData=metaData, design=~groups)\n",
    "dds <- dds[ rowSums( counts(dds) ) >= 10, ]\n",
    "dds <- DESeq(dds)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1ad9fcf-7a6f-428a-857c-88a2ae6e3768",
   "metadata": {},
   "source": [
    "- Obtain results with the FDR (false discovery rate, `alpha`) threshold set to `0.01` and perform `log2FoldChange` shirnkage using _apegml_ method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "84e2e37d-ada8-42d3-8a32-25792d3c74f2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "using 'apeglm' for LFC shrinkage. If used in published research, please cite:\n",
      "    Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for\n",
      "    sequence count data: removing the noise and preserving large differences.\n",
      "    Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895\n",
      "\n"
     ]
    }
   ],
   "source": [
    "res  <- results(dds, alpha=0.01)\n",
    "coef <- resultsNames(dds)[2]\n",
    "res  <- lfcShrink(dds, res=res, coef=coef, type=\"apeglm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7937551a-2da8-46d5-909f-b9efbec5ce5f",
   "metadata": {},
   "source": [
    "- Look up the result summary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f140896d-efa3-4ae3-8e00-bd7119b78187",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "out of 42142 with nonzero total read count\n",
      "adjusted p-value < 0.01\n",
      "LFC > 0 (up)       : 2747, 6.5%\n",
      "LFC < 0 (down)     : 5022, 12%\n",
      "outliers [1]       : 0, 0%\n",
      "low counts [2]     : 9811, 23%\n",
      "(mean count < 0)\n",
      "[1] see 'cooksCutoff' argument of ?results\n",
      "[2] see 'independentFiltering' argument of ?results\n",
      "\n"
     ]
    }
   ],
   "source": [
    "summary(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9499d85f-517e-420e-9432-f023e75d90dc",
   "metadata": {},
   "source": [
    "- Check which rows (gene IDs) were skipped prior performing the DGE analysis, add them to final results and fill with NA values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "3ecc9da5-040b-4934-b3ca-3881eb4f675d",
   "metadata": {},
   "outputs": [],
   "source": [
    "skipped <- setdiff(rownames(countData), rownames(res))\n",
    "res[skipped,] <- NA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5406ef7-6f88-4f72-b1d5-6e70f15a6ed6",
   "metadata": {},
   "source": [
    "- Create a data frame with gene IDs as the first column and DGE results as remianing. Save the complete and final results to a TSV file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "2ec0f059-073f-413f-a73b-9f40afc93438",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message in dir.create(\"output\", recursive = TRUE):\n",
      "“'output' already exists”\n"
     ]
    }
   ],
   "source": [
    "res_df <- cbind(rownames(res), data.frame(res, row.names=NULL))\n",
    "colnames(res_df)[1] <- \"gene_id\"\n",
    "\n",
    "dir.create(\"output\", recursive=TRUE)\n",
    "write.table(res_df, file=\"output/DGE_results.tsv\", sep=\"\\t\", row.names=FALSE, quote=FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "248f664c-71df-4b4a-94d5-6ec9c4d23a06",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.2.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}