{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "setwd('/lustre/scratch117/cellgen/team297/kt16/newcastle_covid')\n",
    "library(ggplot2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "Attaching package: ‘dplyr’\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:stats’:\n",
      "\n",
      "    filter, lag\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:base’:\n",
      "\n",
      "    intersect, setdiff, setequal, union\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(dplyr)\n",
    "makeGeneList <- function(filename){\n",
    "    gl <- readr::read_tsv(filename)\n",
    "    # gl <- readr::read_csv(filename)\n",
    "    y <- grepl('^RPS|^RPL|^MRPL|^MRPS|^MT-', gl$X1)\n",
    "    gl <- gl[!y, ]\n",
    "    gl <- gl %>% dplyr::select(X1, logfoldchanges , pvals)\n",
    "    gl$neglog10pval <- -log10(gl$pvals)\n",
    "    rank <- unlist(gl$neglog10pval*sign(gl$logfoldchanges))\n",
    "    rank[rank == Inf] = 300\n",
    "    rank[rank == -Inf] = -300\n",
    "    names(rank) <- gl$X1\n",
    "    rank <- rev(sort(rank))\n",
    "    gl <- rank\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "celltypes = c('B_exhausted', 'B_immature', 'B_naive', 'B_non-switched_memory', 'B_switched_memory', 'Plasmablast', 'Plasma_cell')\n",
    "groups = c('Critical', 'Severe', 'Moderate', 'Mild', 'Asymptomatic')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n",
      "Warning message:\n",
      "“Missing column names filled in: 'X1' [1]”\n",
      "Parsed with column specification:\n",
      "cols(\n",
      "  X1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  scores = \u001b[32mcol_double()\u001b[39m,\n",
      "  logfoldchanges = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals = \u001b[32mcol_double()\u001b[39m,\n",
      "  pvals_adj = \u001b[32mcol_double()\u001b[39m\n",
      ")\n",
      "\n"
     ]
    }
   ],
   "source": [
    "comparisons = list()\n",
    "for (c in celltypes){\n",
    "    for (g in groups){\n",
    "        comparisons[[c]][[g]] = makeGeneList(paste0('DEG/wilcoxon/',g,'_',c, '_wilcoxon_v2.txt'))\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: Rcpp\n",
      "\n",
      "Warning message:\n",
      "“package ‘Rcpp’ was built under R version 3.6.3”\n"
     ]
    }
   ],
   "source": [
    "library(fgsea)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "h <- as.list(kelvinny::parse_gmt(\"datasets/h.all.v7.2.symbols.gmt\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (25.06% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (24.7% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (22.67% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (26.08% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (30.97% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (21.06% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (19.2% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (18.91% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (20.05% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (23.76% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (15.83% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (15.72% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (16.39% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (16.09% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (16.6% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (23.19% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (22.29% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (20.77% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (20.85% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (27.64% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (19.24% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (17.97% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (17.49% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (18.12% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (21.24% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (25.37% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (23.76% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (23.83% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (25.66% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (34.13% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (20.5% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (19.09% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (19.8% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (20.73% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n",
      "Warning message in fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, :\n",
      "“There are ties in the preranked stats (26.38% of the list).\n",
      "The order of those tied genes will be arbitrary, which may produce unexpected results.”\n"
     ]
    }
   ],
   "source": [
    "res = list()\n",
    "for (c in celltypes){\n",
    "    for (g in groups){\n",
    "        res[[c]][[g]] <- fgsea(pathways = h, stats = comparisons[[c]][[g]], nperm = 10000, minSize = 0, maxSize =1000)\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "for (c in celltypes){\n",
    "    for (g in groups){\n",
    "        res[[c]][[g]] <- as.data.frame(res[[c]][[g]])\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = res\n",
    "\n",
    "# names(result) <- gsub('.csv', \"\",files)\n",
    "for(i in 1:length(comparisons)){\n",
    "    result[[i]] <- lapply(result[[i]], function(x){\n",
    "    x$ranking <- -log10(x$pval)*sign(x$NES) \n",
    "    x <- x[order(x$pathway), ]\n",
    "    return(x)\n",
    "})\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "result <- lapply(result ,function(x){\n",
    "    x[['Critical']]$group = \"Critical\"\n",
    "    x[['Severe']]$group = \"Severe\"\n",
    "    x[['Moderate']]$group = \"Moderate\"\n",
    "    x[['Mild']]$group = \"Mild\"\n",
    "    x[['Asymptomatic']]$group = \"Asymptomatic\"\n",
    "    return(x)\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "result2 <- lapply(result, function(x) {\n",
    "    y <- do.call(rbind, x)\n",
    "    return(y)\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "saveRDS(result2, 'DEG/wilcoxon/Bcell_GSEA_results_v2.RDS')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotGSEA_Hallmark <- function(gsea, group_ref = NULL, cols = NULL, newlabels = NULL, keep_significant_only = FALSE) {\n",
    "    require(ggplot2)\n",
    "    gsea$NES[which(is.na(gsea$NES))] <- 0\n",
    "    gsea$pval[which(is.na(gsea$pval))] <- 1\n",
    "    gsea$padj[which(is.na(gsea$padj))] <- 1\n",
    "    gsea$ranking[which(is.na(gsea$ranking))] <- 0\n",
    "    gsea <- gsea[order(gsea$ranking),]      \n",
    "    gsea_spl <- split(gsea, gsea$group)\n",
    "    if(!is.null(group_ref)){\n",
    "        gsea_spl[[group_ref]] <- gsea_spl[[group_ref]][order(gsea_spl[[group_ref]]$ranking),]\n",
    "        gsea_spl[[group_ref]]$ranking <- gsea_spl[[group_ref]]$ranking*999\n",
    "    } else {\n",
    "        gsea_spl[[2]] <- gsea_spl[[2]]$ranking*999\n",
    "    }\n",
    "    names(gsea_spl) <- NULL\n",
    "\n",
    "    gsea <- do.call(rbind, gsea_spl)\n",
    "    \n",
    "    if (keep_significant_only){\n",
    "        gseax <- split(gsea, gsea$pathway)\n",
    "        for (i in 1:length(gseax)){\n",
    "            if (all(gseax[[i]]$pval >= 0.05)|all(gseax[[i]]$padj >=0.25)){\n",
    "                gseax[[i]] <- NA        \n",
    "            }\n",
    "        }\n",
    "        gseax <- gseax[!is.na(gseax)]\n",
    "        gsea <- do.call(rbind, gseax)\n",
    "    }\n",
    "    gsea <- gsea[order(gsea$ranking), ]\n",
    "    gsea$pathway <- gsub(\"HALLMARK_|\", \"\", gsea$pathway)\n",
    "    gsea$group[which(gsea$pval >= 0.05 & gsea$padj >= 0.25)] <- 'NotSig'\n",
    "    gsea$group[which(gsea$pval < 0.05 & gsea$padj >= 0.25)] <- 'NotSig'\n",
    "    gsea$group[which(gsea$pval >= 0.05)] <- 'NotSig'\n",
    "    gsea$group <- factor(gsea$group, levels = c('NotSig', 'Asymptomatic', 'Mild', 'Moderate', 'Severe', 'Critical'))\n",
    "\n",
    "    x_lim_min <- abs(ceiling(min(-log10(gsea$padj))))\n",
    "    x_lim_max <- abs(ceiling(max(-log10(gsea$padj))))\n",
    "        \n",
    "    if(x_lim_min > x_lim_max){\n",
    "        xval1 <- x_lim_min * -1\n",
    "        xval2 <- x_lim_min\n",
    "    } else {\n",
    "        xval1 <- x_lim_max * -1\n",
    "        xval2 <- x_lim_max\n",
    "    }\n",
    "\n",
    "    if(!is.null(cols)){\n",
    "        gg_color_hue <- function(n) {\n",
    "            hues = seq(15, 375, length = n + 1)\n",
    "            hcl(h = hues, l = 65, c = 100)[1:n]\n",
    "        }\n",
    "        cols. = gg_color_hue(dplyr::n_distinct(gsea$group, na.rm = TRUE))\n",
    "    } else {\n",
    "        cols. = cols\n",
    "    }    \n",
    "    \n",
    "    g <- ggplot(gsea, aes(x = -log10(padj)*sign(NES), y = reorder(pathway, ranking), col = group, size = abs(NES))) + \n",
    "        geom_point() + \n",
    "        labs(x = expression(paste(\"Signed\", \" -log\" [\"10\"], \"adjusted pval\")), y = \"Hallmarks\") +\n",
    "        theme_bw() +\n",
    "        geom_vline(xintercept = 0) +\n",
    "        geom_vline(xintercept = -log10(0.25)) +\n",
    "        geom_vline(xintercept = -log10(0.25)*-1) +\n",
    "        xlim(xval1, xval2) +\n",
    "        scale_size_area(oob = scales::squish, max_size = 3, limits = c(0,2)) +\n",
    "        theme(panel.grid.major = element_blank(), \n",
    "            panel.grid.minor = element_blank(), \n",
    "            panel.background = element_blank(), \n",
    "            axis.line = element_blank(), \n",
    "            axis.ticks = element_blank())\n",
    "    \n",
    "    g$data <- g$data[order(g$data$group, na.last = TRUE), ]\n",
    "    return(g)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "kelvinny::dirCreate('figures/gsea')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Saving 6 x 6.67 in image\n",
      "\n",
      "Saving 6 x 6.67 in image\n",
      "\n",
      "Saving 6 x 6.67 in image\n",
      "\n",
      "Saving 6 x 6.67 in image\n",
      "\n",
      "Saving 6 x 6.67 in image\n",
      "\n",
      "Saving 6 x 6.67 in image\n",
      "\n",
      "Saving 6 x 6.67 in image\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for(i in 1:length(result2)){\n",
    "p <- plotGSEA_Hallmark(result2[[i]], group_ref = \"Moderate\") + scale_color_manual(values = c('#e7e7e7' , '#fed976' ,'#fd8d3c' ,'#e31a1c' ,'#800026' ,'#252525'), drop=FALSE)\n",
    "ggsave(paste0('figures/gsea/',names(result2)[i],\"_v2.pdf\"), plot = p, w = 6, useDingbats = FALSE)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "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>'B_exhausted'</li><li>'B_immature'</li><li>'B_naive'</li><li>'B_non-switched_memory'</li><li>'B_switched_memory'</li><li>'Plasmablast'</li><li>'Plasma_cell'</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 'B\\_exhausted'\n",
       "\\item 'B\\_immature'\n",
       "\\item 'B\\_naive'\n",
       "\\item 'B\\_non-switched\\_memory'\n",
       "\\item 'B\\_switched\\_memory'\n",
       "\\item 'Plasmablast'\n",
       "\\item 'Plasma\\_cell'\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 'B_exhausted'\n",
       "2. 'B_immature'\n",
       "3. 'B_naive'\n",
       "4. 'B_non-switched_memory'\n",
       "5. 'B_switched_memory'\n",
       "6. 'Plasmablast'\n",
       "7. 'Plasma_cell'\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] \"B_exhausted\"           \"B_immature\"            \"B_naive\"              \n",
       "[4] \"B_non-switched_memory\" \"B_switched_memory\"     \"Plasmablast\"          \n",
       "[7] \"Plasma_cell\"          "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "names(result2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp <- c(6, 5, 13, 7, 11, 15, 10)\n",
    "for(i in 1:length(result2)){\n",
    "    p <- plotGSEA_Hallmark(result2[[i]], group_ref = \"Moderate\", keep_significant_only = TRUE) + scale_color_manual(values = c('#e7e7e7','#fed976' ,'#fd8d3c' ,'#e31a1c' ,'#800026' ,'#252525'), drop=FALSE)\n",
    "    ggsave(paste0('figures/gsea/wilcoxon_',names(result2)[i],\"_v3.pdf\"), plot = p, w = 6, h = pp[i]/3, useDingbats = FALSE)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "R version 3.6.1 (2019-07-05)\n",
       "Platform: x86_64-conda_cos6-linux-gnu (64-bit)\n",
       "Running under: Ubuntu 18.04.5 LTS\n",
       "\n",
       "Matrix products: default\n",
       "BLAS/LAPACK: /nfs/users/nfs_k/kt16/miniconda2/envs/dandelion/lib/R/lib/libRblas.so\n",
       "\n",
       "locale:\n",
       " [1] LC_CTYPE=en_GB.utf8       LC_NUMERIC=C             \n",
       " [3] LC_TIME=en_GB.utf8        LC_COLLATE=en_GB.utf8    \n",
       " [5] LC_MONETARY=en_GB.utf8    LC_MESSAGES=en_GB.utf8   \n",
       " [7] LC_PAPER=en_GB.utf8       LC_NAME=C                \n",
       " [9] LC_ADDRESS=C              LC_TELEPHONE=C           \n",
       "[11] LC_MEASUREMENT=en_GB.utf8 LC_IDENTIFICATION=C      \n",
       "\n",
       "attached base packages:\n",
       "[1] stats     graphics  grDevices utils     datasets  methods   base     \n",
       "\n",
       "other attached packages:\n",
       "[1] fgsea_1.12.0  Rcpp_1.0.5    dplyr_1.0.0   ggplot2_3.3.2\n",
       "\n",
       "loaded via a namespace (and not attached):\n",
       " [1] bitops_1.0-6                matrixStats_0.57.0         \n",
       " [3] fs_1.5.0                    doMC_1.3.7                 \n",
       " [5] usethis_1.6.3               devtools_2.3.2             \n",
       " [7] RColorBrewer_1.1-2          rprojroot_1.3-2            \n",
       " [9] GenomeInfoDb_1.22.1         repr_1.1.0                 \n",
       "[11] tools_3.6.1                 backports_1.2.0            \n",
       "[13] R6_2.5.0                    BiocGenerics_0.32.0        \n",
       "[15] colorspace_1.4-1            withr_2.3.0                \n",
       "[17] tidyselect_1.1.0            gridExtra_2.3              \n",
       "[19] prettyunits_1.1.1           processx_3.4.4             \n",
       "[21] compiler_3.6.1              glmnet_4.0-2               \n",
       "[23] cli_2.1.0                   Biobase_2.46.0             \n",
       "[25] desc_1.2.0                  DelayedArray_0.12.3        \n",
       "[27] labeling_0.4.2              scales_1.1.1               \n",
       "[29] readr_1.3.1                 callr_3.5.1                \n",
       "[31] pbdZMQ_0.3-3                digest_0.6.27              \n",
       "[33] XVector_0.26.0              base64enc_0.1-3            \n",
       "[35] pkgconfig_2.0.3             htmltools_0.5.0            \n",
       "[37] sessioninfo_1.1.1           rlang_0.4.10               \n",
       "[39] farver_2.0.3                shape_1.4.5                \n",
       "[41] generics_0.0.2              jsonlite_1.7.1             \n",
       "[43] BiocParallel_1.20.1         RCurl_1.98-1.2             \n",
       "[45] magrittr_1.5                GenomeInfoDbData_1.2.2     \n",
       "[47] Matrix_1.2-18               IRkernel_1.1.1             \n",
       "[49] munsell_0.5.0               S4Vectors_0.24.4           \n",
       "[51] fansi_0.4.1                 lifecycle_0.2.0            \n",
       "[53] kelvinny_2.0.5              SummarizedExperiment_1.16.1\n",
       "[55] zlibbioc_1.32.0             pkgbuild_1.1.0             \n",
       "[57] grid_3.6.1                  parallel_3.6.1             \n",
       "[59] crayon_1.3.4                lattice_0.20-41            \n",
       "[61] IRdisplay_0.7.0             splines_3.6.1              \n",
       "[63] hms_0.5.2                   ps_1.4.0                   \n",
       "[65] pillar_1.4.6                ranger_0.12.1              \n",
       "[67] GenomicRanges_1.38.0        uuid_0.1-4                 \n",
       "[69] codetools_0.2-17            stats4_3.6.1               \n",
       "[71] pkgload_1.1.0               fastmatch_1.1-0            \n",
       "[73] glue_1.4.1                  evaluate_0.14              \n",
       "[75] data.table_1.12.8           remotes_2.2.0              \n",
       "[77] vctrs_0.3.1                 foreach_1.5.0              \n",
       "[79] testthat_3.0.0              gtable_0.3.0               \n",
       "[81] purrr_0.3.4                 assertthat_0.2.1           \n",
       "[83] viridisLite_0.3.0           survival_3.2-7             \n",
       "[85] SingleCellExperiment_1.8.0  tibble_3.0.2               \n",
       "[87] pheatmap_1.0.12             iterators_1.0.12           \n",
       "[89] memoise_1.1.0               IRanges_2.20.2             \n",
       "[91] ellipsis_0.3.1             "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sessionInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
