{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "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": [
    "# set up the model\n",
    "setwd('/lustre/scratch117/cellgen/team297/kt16/newcastle_covid')\n",
    "library(dplyr)\n",
    "\n",
    "bcell.df <- read.csv(\"h5ad/covid_jan_2021_bcells_obs.csv\")\n",
    "clinical_info = read.csv('other_analyses/final_metadata_feb21.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "covid.meta <- bcell.df[,c('sample_id', 'patient_id', 'Status_on_day_collection_summary_v2', 'Age', 'Sex', 'Collection_Day', 'Days_from_onset')] %>% distinct()\n",
    "rownames(covid.meta) <- covid.meta$sample_id\n",
    "covid.meta$Age = clinical_info$Age[match(covid.meta$patient_id, clinical_info$patient_id)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "bcell.meta <- covid.meta[!covid.meta$Status_on_day_collection_summary_v2 %in% c(\"Non_covid\", \"LPS\", 'Healthy'),]\n",
    "bcell.meta$OrderedSeverity <- ordered(bcell.meta$Status_on_day_collection_summary_v2, levels=c(\"Asymptomatic\", \"Mild\", \"Moderate\", \"Severe\", \"Critical\"))\n",
    "bcell.meta['Days_from_onset'][bcell.meta['Days_from_onset'] == 'Not_known'] <- NA\n",
    "bcell.meta['Days_from_onset'][bcell.meta['Days_from_onset'] == 'nan'] <- NA\n",
    "bcell.meta['Days_from_onset'][bcell.meta['Days_from_onset'] == 'Healthy'] <- 0\n",
    "bcell.meta <- bcell.meta[!is.na(bcell.meta['Days_from_onset']),]\n",
    "bcell.df <- bcell.df %>% filter(sample_id %in% row.names(bcell.meta))\n",
    "bcell.model <- model.matrix(~ Sex + Age + as.numeric(Days_from_onset) + OrderedSeverity, data=bcell.meta[bcell.meta$Collection_Day %in% c(\"D0\"), ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# count cells\n",
    "cell.freq.tab <- t(table(bcell.df$sample_id[bcell.df$Collection_Day %in% c(\"D0\") &\n",
    "                                                        !bcell.df$Status_on_day_collection_summary_v2 %in% c(\"LPS\", \"Non_covid\", 'Healthy')],\n",
    "                         bcell.df$celltype_B[bcell.df$Collection_Day %in% c(\"D0\") &\n",
    "                                                        !bcell.df$Status_on_day_collection_summary_v2 %in% c(\"LPS\", \"Non_covid\", 'Healthy')]))\n",
    "cell.freq.tab <- cell.freq.tab[,colSums(cell.freq.tab) != 0]\n",
    "bcell.model <- bcell.model[colnames(cell.freq.tab), ]\n",
    "n.cell.vecc <- table(bcell.df$sample_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: limma\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(edgeR)\n",
    "bcell.dge <- DGEList(cell.freq.tab, lib.size=log(n.cell.vecc[colnames(cell.freq.tab)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#estimate dispersions and fit model\n",
    "bcell.dge <- estimateDisp(bcell.dge, design=bcell.model)\n",
    "bcell.linear.fit <- glmQLFit(bcell.dge, bcell.model, robust=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 9 × 8</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>logFC</th><th scope=col>logCPM</th><th scope=col>F</th><th scope=col>PValue</th><th scope=col>FDR</th><th scope=col>CellType</th><th scope=col>Sig</th><th scope=col>Diff</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>B_exhausted</th><td>-0.022000089</td><td>15.47603</td><td>2.5857822</td><td>0.11149871</td><td>0.1672481</td><td>B_exhausted          </td><td>0</td><td> 0</td></tr>\n",
       "\t<tr><th scope=row>B_immature</th><td>-0.009951020</td><td>16.09104</td><td>0.7288711</td><td>0.39562518</td><td>0.3956252</td><td>B_immature           </td><td>0</td><td> 0</td></tr>\n",
       "\t<tr><th scope=row>B_naive</th><td> 0.004772062</td><td>19.17906</td><td>1.2735821</td><td>0.26223948</td><td>0.3371650</td><td>B_naive              </td><td>0</td><td> 0</td></tr>\n",
       "\t<tr><th scope=row>B_non-switched_memory</th><td> 0.016367280</td><td>15.53889</td><td>2.7378201</td><td>0.10164998</td><td>0.1672481</td><td>B_non-switched_memory</td><td>0</td><td> 0</td></tr>\n",
       "\t<tr><th scope=row>B_switched_memory</th><td> 0.019809624</td><td>16.57053</td><td>4.4268053</td><td>0.03830350</td><td>0.1149105</td><td>B_switched_memory    </td><td>0</td><td> 0</td></tr>\n",
       "\t<tr><th scope=row>Plasma_cell_IgA</th><td>-0.035625748</td><td>15.44095</td><td>6.4128840</td><td>0.01314884</td><td>0.0591698</td><td>Plasma_cell_IgA      </td><td>1</td><td>-1</td></tr>\n",
       "\t<tr><th scope=row>Plasma_cell_IgG</th><td>-0.015731866</td><td>15.62167</td><td>1.0443676</td><td>0.30967893</td><td>0.3483888</td><td>Plasma_cell_IgG      </td><td>0</td><td> 0</td></tr>\n",
       "\t<tr><th scope=row>Plasma_cell_IgM</th><td>-0.027554743</td><td>14.21102</td><td>3.4250695</td><td>0.06765268</td><td>0.1522185</td><td>Plasma_cell_IgM      </td><td>0</td><td> 0</td></tr>\n",
       "\t<tr><th scope=row>Plasmablast</th><td>-0.050900382</td><td>15.66965</td><td>6.8703183</td><td>0.01036463</td><td>0.0591698</td><td>Plasmablast          </td><td>1</td><td>-1</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 9 × 8\n",
       "\\begin{tabular}{r|llllllll}\n",
       "  & logFC & logCPM & F & PValue & FDR & CellType & Sig & Diff\\\\\n",
       "  & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <chr> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\tB\\_exhausted & -0.022000089 & 15.47603 & 2.5857822 & 0.11149871 & 0.1672481 & B\\_exhausted           & 0 &  0\\\\\n",
       "\tB\\_immature & -0.009951020 & 16.09104 & 0.7288711 & 0.39562518 & 0.3956252 & B\\_immature            & 0 &  0\\\\\n",
       "\tB\\_naive &  0.004772062 & 19.17906 & 1.2735821 & 0.26223948 & 0.3371650 & B\\_naive               & 0 &  0\\\\\n",
       "\tB\\_non-switched\\_memory &  0.016367280 & 15.53889 & 2.7378201 & 0.10164998 & 0.1672481 & B\\_non-switched\\_memory & 0 &  0\\\\\n",
       "\tB\\_switched\\_memory &  0.019809624 & 16.57053 & 4.4268053 & 0.03830350 & 0.1149105 & B\\_switched\\_memory     & 0 &  0\\\\\n",
       "\tPlasma\\_cell\\_IgA & -0.035625748 & 15.44095 & 6.4128840 & 0.01314884 & 0.0591698 & Plasma\\_cell\\_IgA       & 1 & -1\\\\\n",
       "\tPlasma\\_cell\\_IgG & -0.015731866 & 15.62167 & 1.0443676 & 0.30967893 & 0.3483888 & Plasma\\_cell\\_IgG       & 0 &  0\\\\\n",
       "\tPlasma\\_cell\\_IgM & -0.027554743 & 14.21102 & 3.4250695 & 0.06765268 & 0.1522185 & Plasma\\_cell\\_IgM       & 0 &  0\\\\\n",
       "\tPlasmablast & -0.050900382 & 15.66965 & 6.8703183 & 0.01036463 & 0.0591698 & Plasmablast           & 1 & -1\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 9 × 8\n",
       "\n",
       "| <!--/--> | logFC &lt;dbl&gt; | logCPM &lt;dbl&gt; | F &lt;dbl&gt; | PValue &lt;dbl&gt; | FDR &lt;dbl&gt; | CellType &lt;chr&gt; | Sig &lt;dbl&gt; | Diff &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|\n",
       "| B_exhausted | -0.022000089 | 15.47603 | 2.5857822 | 0.11149871 | 0.1672481 | B_exhausted           | 0 |  0 |\n",
       "| B_immature | -0.009951020 | 16.09104 | 0.7288711 | 0.39562518 | 0.3956252 | B_immature            | 0 |  0 |\n",
       "| B_naive |  0.004772062 | 19.17906 | 1.2735821 | 0.26223948 | 0.3371650 | B_naive               | 0 |  0 |\n",
       "| B_non-switched_memory |  0.016367280 | 15.53889 | 2.7378201 | 0.10164998 | 0.1672481 | B_non-switched_memory | 0 |  0 |\n",
       "| B_switched_memory |  0.019809624 | 16.57053 | 4.4268053 | 0.03830350 | 0.1149105 | B_switched_memory     | 0 |  0 |\n",
       "| Plasma_cell_IgA | -0.035625748 | 15.44095 | 6.4128840 | 0.01314884 | 0.0591698 | Plasma_cell_IgA       | 1 | -1 |\n",
       "| Plasma_cell_IgG | -0.015731866 | 15.62167 | 1.0443676 | 0.30967893 | 0.3483888 | Plasma_cell_IgG       | 0 |  0 |\n",
       "| Plasma_cell_IgM | -0.027554743 | 14.21102 | 3.4250695 | 0.06765268 | 0.1522185 | Plasma_cell_IgM       | 0 |  0 |\n",
       "| Plasmablast | -0.050900382 | 15.66965 | 6.8703183 | 0.01036463 | 0.0591698 | Plasmablast           | 1 | -1 |\n",
       "\n"
      ],
      "text/plain": [
       "                      logFC        logCPM   F         PValue     FDR      \n",
       "B_exhausted           -0.022000089 15.47603 2.5857822 0.11149871 0.1672481\n",
       "B_immature            -0.009951020 16.09104 0.7288711 0.39562518 0.3956252\n",
       "B_naive                0.004772062 19.17906 1.2735821 0.26223948 0.3371650\n",
       "B_non-switched_memory  0.016367280 15.53889 2.7378201 0.10164998 0.1672481\n",
       "B_switched_memory      0.019809624 16.57053 4.4268053 0.03830350 0.1149105\n",
       "Plasma_cell_IgA       -0.035625748 15.44095 6.4128840 0.01314884 0.0591698\n",
       "Plasma_cell_IgG       -0.015731866 15.62167 1.0443676 0.30967893 0.3483888\n",
       "Plasma_cell_IgM       -0.027554743 14.21102 3.4250695 0.06765268 0.1522185\n",
       "Plasmablast           -0.050900382 15.66965 6.8703183 0.01036463 0.0591698\n",
       "                      CellType              Sig Diff\n",
       "B_exhausted           B_exhausted           0    0  \n",
       "B_immature            B_immature            0    0  \n",
       "B_naive               B_naive               0    0  \n",
       "B_non-switched_memory B_non-switched_memory 0    0  \n",
       "B_switched_memory     B_switched_memory     0    0  \n",
       "Plasma_cell_IgA       Plasma_cell_IgA       1   -1  \n",
       "Plasma_cell_IgG       Plasma_cell_IgG       0    0  \n",
       "Plasma_cell_IgM       Plasma_cell_IgM       0    0  \n",
       "Plasmablast           Plasmablast           1   -1  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# linear changes\n",
    "bcell.res <- as.data.frame(topTags(glmQLFTest(bcell.linear.fit, coef=4), sort.by='none', n=Inf))\n",
    "bcell.res$CellType <- rownames(bcell.res)\n",
    "bcell.res$Sig <- as.numeric(bcell.res$FDR < 0.1 & bcell.res$PValue < 0.05)\n",
    "bcell.res$Diff <- sign(bcell.res$logFC)\n",
    "bcell.res$Diff[bcell.res$FDR >= 0.1 | bcell.res$PValue >= 0.05] <- 0\n",
    "bcell.res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 9 × 8</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>logFC</th><th scope=col>logCPM</th><th scope=col>F</th><th scope=col>PValue</th><th scope=col>FDR</th><th scope=col>CellType</th><th scope=col>Sig</th><th scope=col>Diff</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>B_exhausted</th><td> 0.50638457</td><td>15.47603</td><td> 2.13745812</td><td>0.147390797</td><td>0.223010196</td><td>B_exhausted          </td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>B_immature</th><td>-0.03605731</td><td>16.09104</td><td> 0.01579397</td><td>0.900283778</td><td>0.900283778</td><td>B_immature           </td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>B_naive</th><td>-0.15756913</td><td>19.17906</td><td> 1.88402647</td><td>0.173452375</td><td>0.223010196</td><td>B_naive              </td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>B_non-switched_memory</th><td>-0.25783179</td><td>15.53889</td><td> 0.95509224</td><td>0.331172309</td><td>0.372568847</td><td>B_non-switched_memory</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>B_switched_memory</th><td>-0.45695532</td><td>16.57053</td><td> 3.32389167</td><td>0.071760524</td><td>0.161461178</td><td>B_switched_memory    </td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>Plasma_cell_IgA</th><td> 0.50037281</td><td>15.44095</td><td> 2.00559435</td><td>0.160337879</td><td>0.223010196</td><td>Plasma_cell_IgA      </td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>Plasma_cell_IgG</th><td> 0.80302329</td><td>15.62167</td><td> 4.25417278</td><td>0.042176302</td><td>0.126528905</td><td>Plasma_cell_IgG      </td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>Plasma_cell_IgM</th><td> 1.57655146</td><td>14.21102</td><td>14.74752590</td><td>0.000234591</td><td>0.002111319</td><td>Plasma_cell_IgM      </td><td>1</td><td>1</td></tr>\n",
       "\t<tr><th scope=row>Plasmablast</th><td> 1.53639973</td><td>15.66965</td><td>10.68238382</td><td>0.001556629</td><td>0.007004830</td><td>Plasmablast          </td><td>1</td><td>1</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 9 × 8\n",
       "\\begin{tabular}{r|llllllll}\n",
       "  & logFC & logCPM & F & PValue & FDR & CellType & Sig & Diff\\\\\n",
       "  & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <chr> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\tB\\_exhausted &  0.50638457 & 15.47603 &  2.13745812 & 0.147390797 & 0.223010196 & B\\_exhausted           & 0 & 0\\\\\n",
       "\tB\\_immature & -0.03605731 & 16.09104 &  0.01579397 & 0.900283778 & 0.900283778 & B\\_immature            & 0 & 0\\\\\n",
       "\tB\\_naive & -0.15756913 & 19.17906 &  1.88402647 & 0.173452375 & 0.223010196 & B\\_naive               & 0 & 0\\\\\n",
       "\tB\\_non-switched\\_memory & -0.25783179 & 15.53889 &  0.95509224 & 0.331172309 & 0.372568847 & B\\_non-switched\\_memory & 0 & 0\\\\\n",
       "\tB\\_switched\\_memory & -0.45695532 & 16.57053 &  3.32389167 & 0.071760524 & 0.161461178 & B\\_switched\\_memory     & 0 & 0\\\\\n",
       "\tPlasma\\_cell\\_IgA &  0.50037281 & 15.44095 &  2.00559435 & 0.160337879 & 0.223010196 & Plasma\\_cell\\_IgA       & 0 & 0\\\\\n",
       "\tPlasma\\_cell\\_IgG &  0.80302329 & 15.62167 &  4.25417278 & 0.042176302 & 0.126528905 & Plasma\\_cell\\_IgG       & 0 & 0\\\\\n",
       "\tPlasma\\_cell\\_IgM &  1.57655146 & 14.21102 & 14.74752590 & 0.000234591 & 0.002111319 & Plasma\\_cell\\_IgM       & 1 & 1\\\\\n",
       "\tPlasmablast &  1.53639973 & 15.66965 & 10.68238382 & 0.001556629 & 0.007004830 & Plasmablast           & 1 & 1\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 9 × 8\n",
       "\n",
       "| <!--/--> | logFC &lt;dbl&gt; | logCPM &lt;dbl&gt; | F &lt;dbl&gt; | PValue &lt;dbl&gt; | FDR &lt;dbl&gt; | CellType &lt;chr&gt; | Sig &lt;dbl&gt; | Diff &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|\n",
       "| B_exhausted |  0.50638457 | 15.47603 |  2.13745812 | 0.147390797 | 0.223010196 | B_exhausted           | 0 | 0 |\n",
       "| B_immature | -0.03605731 | 16.09104 |  0.01579397 | 0.900283778 | 0.900283778 | B_immature            | 0 | 0 |\n",
       "| B_naive | -0.15756913 | 19.17906 |  1.88402647 | 0.173452375 | 0.223010196 | B_naive               | 0 | 0 |\n",
       "| B_non-switched_memory | -0.25783179 | 15.53889 |  0.95509224 | 0.331172309 | 0.372568847 | B_non-switched_memory | 0 | 0 |\n",
       "| B_switched_memory | -0.45695532 | 16.57053 |  3.32389167 | 0.071760524 | 0.161461178 | B_switched_memory     | 0 | 0 |\n",
       "| Plasma_cell_IgA |  0.50037281 | 15.44095 |  2.00559435 | 0.160337879 | 0.223010196 | Plasma_cell_IgA       | 0 | 0 |\n",
       "| Plasma_cell_IgG |  0.80302329 | 15.62167 |  4.25417278 | 0.042176302 | 0.126528905 | Plasma_cell_IgG       | 0 | 0 |\n",
       "| Plasma_cell_IgM |  1.57655146 | 14.21102 | 14.74752590 | 0.000234591 | 0.002111319 | Plasma_cell_IgM       | 1 | 1 |\n",
       "| Plasmablast |  1.53639973 | 15.66965 | 10.68238382 | 0.001556629 | 0.007004830 | Plasmablast           | 1 | 1 |\n",
       "\n"
      ],
      "text/plain": [
       "                      logFC       logCPM   F           PValue      FDR        \n",
       "B_exhausted            0.50638457 15.47603  2.13745812 0.147390797 0.223010196\n",
       "B_immature            -0.03605731 16.09104  0.01579397 0.900283778 0.900283778\n",
       "B_naive               -0.15756913 19.17906  1.88402647 0.173452375 0.223010196\n",
       "B_non-switched_memory -0.25783179 15.53889  0.95509224 0.331172309 0.372568847\n",
       "B_switched_memory     -0.45695532 16.57053  3.32389167 0.071760524 0.161461178\n",
       "Plasma_cell_IgA        0.50037281 15.44095  2.00559435 0.160337879 0.223010196\n",
       "Plasma_cell_IgG        0.80302329 15.62167  4.25417278 0.042176302 0.126528905\n",
       "Plasma_cell_IgM        1.57655146 14.21102 14.74752590 0.000234591 0.002111319\n",
       "Plasmablast            1.53639973 15.66965 10.68238382 0.001556629 0.007004830\n",
       "                      CellType              Sig Diff\n",
       "B_exhausted           B_exhausted           0   0   \n",
       "B_immature            B_immature            0   0   \n",
       "B_naive               B_naive               0   0   \n",
       "B_non-switched_memory B_non-switched_memory 0   0   \n",
       "B_switched_memory     B_switched_memory     0   0   \n",
       "Plasma_cell_IgA       Plasma_cell_IgA       0   0   \n",
       "Plasma_cell_IgG       Plasma_cell_IgG       0   0   \n",
       "Plasma_cell_IgM       Plasma_cell_IgM       1   1   \n",
       "Plasmablast           Plasmablast           1   1   "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Quadratic changes.\n",
    "bcell.quad.res <- as.data.frame(topTags(glmQLFTest(bcell.linear.fit, coef=5), sort.by='none', n=Inf))\n",
    "bcell.quad.res$CellType <- rownames(bcell.quad.res)\n",
    "bcell.quad.res$Sig <- as.numeric(bcell.quad.res$FDR < 0.1 & bcell.quad.res$PValue < 0.05)\n",
    "bcell.quad.res$Diff <- sign(bcell.quad.res$logFC)\n",
    "bcell.quad.res$Diff[bcell.quad.res$FDR >= 0.1 | bcell.quad.res$PValue >= 0.05] <- 0\n",
    "bcell.quad.res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## remove Critical group to check for robustness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "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": [
    "# set up the model\n",
    "setwd('/lustre/scratch117/cellgen/team297/kt16/newcastle_covid')\n",
    "library(dplyr)\n",
    "\n",
    "bcell.df <- read.csv(\"h5ad/covid_jan_2021_bcells_obs.csv\")\n",
    "clinical_info = read.csv('other_analyses/final_metadata_feb21.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "covid.meta <- bcell.df[,c('sample_id', 'patient_id', 'Status_on_day_collection_summary_v2', 'Age', 'Sex', 'Collection_Day', 'Days_from_onset')] %>% distinct()\n",
    "rownames(covid.meta) <- covid.meta$sample_id\n",
    "covid.meta$Age = clinical_info$Age[match(covid.meta$patient_id, clinical_info$patient_id)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "bcell.meta <- covid.meta[!covid.meta$Status_on_day_collection_summary_v2 %in% c(\"Non_covid\", \"LPS\", 'Healthy', 'Critical'),]\n",
    "bcell.meta$OrderedSeverity <- ordered(bcell.meta$Status_on_day_collection_summary_v2, levels=c(\"Asymptomatic\", \"Mild\", \"Moderate\", \"Severe\"))\n",
    "bcell.meta['Days_from_onset'][bcell.meta['Days_from_onset'] == 'Not_known'] <- NA\n",
    "bcell.meta['Days_from_onset'][bcell.meta['Days_from_onset'] == 'nan'] <- NA\n",
    "bcell.meta['Days_from_onset'][bcell.meta['Days_from_onset'] == 'Healthy'] <- 0\n",
    "bcell.meta <- bcell.meta[!is.na(bcell.meta['Days_from_onset']),]\n",
    "bcell.df <- bcell.df %>% filter(sample_id %in% row.names(bcell.meta))\n",
    "bcell.model <- model.matrix(~ Sex + Age + as.numeric(Days_from_onset) + OrderedSeverity, data=bcell.meta[bcell.meta$Collection_Day %in% c(\"D0\"), ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# count cells\n",
    "cell.freq.tab <- t(table(bcell.df$sample_id[bcell.df$Collection_Day %in% c(\"D0\") &\n",
    "                                                        !bcell.df$Status_on_day_collection_summary_v2 %in% c(\"LPS\", \"Non_covid\", 'Healthy')],\n",
    "                         bcell.df$celltype_B[bcell.df$Collection_Day %in% c(\"D0\") &\n",
    "                                                        !bcell.df$Status_on_day_collection_summary_v2 %in% c(\"LPS\", \"Non_covid\", 'Healthy')]))\n",
    "cell.freq.tab <- cell.freq.tab[,colSums(cell.freq.tab) != 0]\n",
    "bcell.model <- bcell.model[colnames(cell.freq.tab), ]\n",
    "n.cell.vecc <- table(bcell.df$sample_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: limma\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(edgeR)\n",
    "bcell.dge <- DGEList(cell.freq.tab, lib.size=log(n.cell.vecc[colnames(cell.freq.tab)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#estimate dispersions and fit model\n",
    "bcell.dge <- estimateDisp(bcell.dge, design=bcell.model)\n",
    "bcell.linear.fit <- glmQLFit(bcell.dge, bcell.model, robust=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 9 × 8</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>logFC</th><th scope=col>logCPM</th><th scope=col>F</th><th scope=col>PValue</th><th scope=col>FDR</th><th scope=col>CellType</th><th scope=col>Sig</th><th scope=col>Diff</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>B_exhausted</th><td>-0.0334312670</td><td>15.45495</td><td>3.282996270</td><td>0.07390033</td><td>0.3145728</td><td>B_exhausted          </td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>B_immature</th><td> 0.0027339951</td><td>16.07439</td><td>0.031530621</td><td>0.85952810</td><td>0.9228529</td><td>B_immature           </td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>B_naive</th><td> 0.0005444977</td><td>19.17061</td><td>0.009439796</td><td>0.92285292</td><td>0.9228529</td><td>B_naive              </td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>B_non-switched_memory</th><td> 0.0184035709</td><td>15.53814</td><td>2.039204173</td><td>0.15733381</td><td>0.3145728</td><td>B_non-switched_memory</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>B_switched_memory</th><td> 0.0194959951</td><td>16.59030</td><td>2.816523550</td><td>0.09735324</td><td>0.3145728</td><td>B_switched_memory    </td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>Plasma_cell_IgA</th><td>-0.0235493759</td><td>15.46497</td><td>1.674859422</td><td>0.19947747</td><td>0.3145728</td><td>Plasma_cell_IgA      </td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>Plasma_cell_IgG</th><td>-0.0023611299</td><td>15.66264</td><td>0.015098097</td><td>0.90252682</td><td>0.9228529</td><td>Plasma_cell_IgG      </td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>Plasma_cell_IgM</th><td>-0.0274229775</td><td>14.24536</td><td>1.769025640</td><td>0.18742947</td><td>0.3145728</td><td>Plasma_cell_IgM      </td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>Plasmablast</th><td>-0.0320698858</td><td>15.68467</td><td>1.600023704</td><td>0.20971518</td><td>0.3145728</td><td>Plasmablast          </td><td>0</td><td>0</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 9 × 8\n",
       "\\begin{tabular}{r|llllllll}\n",
       "  & logFC & logCPM & F & PValue & FDR & CellType & Sig & Diff\\\\\n",
       "  & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <chr> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\tB\\_exhausted & -0.0334312670 & 15.45495 & 3.282996270 & 0.07390033 & 0.3145728 & B\\_exhausted           & 0 & 0\\\\\n",
       "\tB\\_immature &  0.0027339951 & 16.07439 & 0.031530621 & 0.85952810 & 0.9228529 & B\\_immature            & 0 & 0\\\\\n",
       "\tB\\_naive &  0.0005444977 & 19.17061 & 0.009439796 & 0.92285292 & 0.9228529 & B\\_naive               & 0 & 0\\\\\n",
       "\tB\\_non-switched\\_memory &  0.0184035709 & 15.53814 & 2.039204173 & 0.15733381 & 0.3145728 & B\\_non-switched\\_memory & 0 & 0\\\\\n",
       "\tB\\_switched\\_memory &  0.0194959951 & 16.59030 & 2.816523550 & 0.09735324 & 0.3145728 & B\\_switched\\_memory     & 0 & 0\\\\\n",
       "\tPlasma\\_cell\\_IgA & -0.0235493759 & 15.46497 & 1.674859422 & 0.19947747 & 0.3145728 & Plasma\\_cell\\_IgA       & 0 & 0\\\\\n",
       "\tPlasma\\_cell\\_IgG & -0.0023611299 & 15.66264 & 0.015098097 & 0.90252682 & 0.9228529 & Plasma\\_cell\\_IgG       & 0 & 0\\\\\n",
       "\tPlasma\\_cell\\_IgM & -0.0274229775 & 14.24536 & 1.769025640 & 0.18742947 & 0.3145728 & Plasma\\_cell\\_IgM       & 0 & 0\\\\\n",
       "\tPlasmablast & -0.0320698858 & 15.68467 & 1.600023704 & 0.20971518 & 0.3145728 & Plasmablast           & 0 & 0\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 9 × 8\n",
       "\n",
       "| <!--/--> | logFC &lt;dbl&gt; | logCPM &lt;dbl&gt; | F &lt;dbl&gt; | PValue &lt;dbl&gt; | FDR &lt;dbl&gt; | CellType &lt;chr&gt; | Sig &lt;dbl&gt; | Diff &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|\n",
       "| B_exhausted | -0.0334312670 | 15.45495 | 3.282996270 | 0.07390033 | 0.3145728 | B_exhausted           | 0 | 0 |\n",
       "| B_immature |  0.0027339951 | 16.07439 | 0.031530621 | 0.85952810 | 0.9228529 | B_immature            | 0 | 0 |\n",
       "| B_naive |  0.0005444977 | 19.17061 | 0.009439796 | 0.92285292 | 0.9228529 | B_naive               | 0 | 0 |\n",
       "| B_non-switched_memory |  0.0184035709 | 15.53814 | 2.039204173 | 0.15733381 | 0.3145728 | B_non-switched_memory | 0 | 0 |\n",
       "| B_switched_memory |  0.0194959951 | 16.59030 | 2.816523550 | 0.09735324 | 0.3145728 | B_switched_memory     | 0 | 0 |\n",
       "| Plasma_cell_IgA | -0.0235493759 | 15.46497 | 1.674859422 | 0.19947747 | 0.3145728 | Plasma_cell_IgA       | 0 | 0 |\n",
       "| Plasma_cell_IgG | -0.0023611299 | 15.66264 | 0.015098097 | 0.90252682 | 0.9228529 | Plasma_cell_IgG       | 0 | 0 |\n",
       "| Plasma_cell_IgM | -0.0274229775 | 14.24536 | 1.769025640 | 0.18742947 | 0.3145728 | Plasma_cell_IgM       | 0 | 0 |\n",
       "| Plasmablast | -0.0320698858 | 15.68467 | 1.600023704 | 0.20971518 | 0.3145728 | Plasmablast           | 0 | 0 |\n",
       "\n"
      ],
      "text/plain": [
       "                      logFC         logCPM   F           PValue     FDR      \n",
       "B_exhausted           -0.0334312670 15.45495 3.282996270 0.07390033 0.3145728\n",
       "B_immature             0.0027339951 16.07439 0.031530621 0.85952810 0.9228529\n",
       "B_naive                0.0005444977 19.17061 0.009439796 0.92285292 0.9228529\n",
       "B_non-switched_memory  0.0184035709 15.53814 2.039204173 0.15733381 0.3145728\n",
       "B_switched_memory      0.0194959951 16.59030 2.816523550 0.09735324 0.3145728\n",
       "Plasma_cell_IgA       -0.0235493759 15.46497 1.674859422 0.19947747 0.3145728\n",
       "Plasma_cell_IgG       -0.0023611299 15.66264 0.015098097 0.90252682 0.9228529\n",
       "Plasma_cell_IgM       -0.0274229775 14.24536 1.769025640 0.18742947 0.3145728\n",
       "Plasmablast           -0.0320698858 15.68467 1.600023704 0.20971518 0.3145728\n",
       "                      CellType              Sig Diff\n",
       "B_exhausted           B_exhausted           0   0   \n",
       "B_immature            B_immature            0   0   \n",
       "B_naive               B_naive               0   0   \n",
       "B_non-switched_memory B_non-switched_memory 0   0   \n",
       "B_switched_memory     B_switched_memory     0   0   \n",
       "Plasma_cell_IgA       Plasma_cell_IgA       0   0   \n",
       "Plasma_cell_IgG       Plasma_cell_IgG       0   0   \n",
       "Plasma_cell_IgM       Plasma_cell_IgM       0   0   \n",
       "Plasmablast           Plasmablast           0   0   "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# linear changes\n",
    "bcell.res <- as.data.frame(topTags(glmQLFTest(bcell.linear.fit, coef=4), sort.by='none', n=Inf))\n",
    "bcell.res$CellType <- rownames(bcell.res)\n",
    "bcell.res$Sig <- as.numeric(bcell.res$FDR < 0.1 & bcell.res$PValue < 0.05)\n",
    "bcell.res$Diff <- sign(bcell.res$logFC)\n",
    "bcell.res$Diff[bcell.res$FDR >= 0.1 | bcell.res$PValue >= 0.05] <- 0\n",
    "bcell.res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 9 × 8</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>logFC</th><th scope=col>logCPM</th><th scope=col>F</th><th scope=col>PValue</th><th scope=col>FDR</th><th scope=col>CellType</th><th scope=col>Sig</th><th scope=col>Diff</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>B_exhausted</th><td> 0.1616904</td><td>15.45495</td><td> 0.2760204</td><td>6.008314e-01</td><td>6.008314e-01</td><td>B_exhausted          </td><td>0</td><td> 0</td></tr>\n",
       "\t<tr><th scope=row>B_immature</th><td>-0.3098771</td><td>16.07439</td><td> 1.4774286</td><td>2.278914e-01</td><td>3.418371e-01</td><td>B_immature           </td><td>0</td><td> 0</td></tr>\n",
       "\t<tr><th scope=row>B_naive</th><td>-0.2113975</td><td>19.17061</td><td> 4.3488903</td><td>4.034672e-02</td><td>7.262410e-02</td><td>B_naive              </td><td>1</td><td>-1</td></tr>\n",
       "\t<tr><th scope=row>B_non-switched_memory</th><td>-0.1353764</td><td>15.53814</td><td> 0.3330929</td><td>5.655274e-01</td><td>6.008314e-01</td><td>B_non-switched_memory</td><td>0</td><td> 0</td></tr>\n",
       "\t<tr><th scope=row>B_switched_memory</th><td>-0.1482451</td><td>16.59030</td><td> 0.4909614</td><td>4.856110e-01</td><td>6.008314e-01</td><td>B_switched_memory    </td><td>0</td><td> 0</td></tr>\n",
       "\t<tr><th scope=row>Plasma_cell_IgA</th><td> 0.7192348</td><td>15.46497</td><td> 4.8448680</td><td>3.072471e-02</td><td>6.913061e-02</td><td>Plasma_cell_IgA      </td><td>1</td><td> 1</td></tr>\n",
       "\t<tr><th scope=row>Plasma_cell_IgG</th><td> 1.2398394</td><td>15.66264</td><td>11.7834434</td><td>9.646690e-04</td><td>2.894007e-03</td><td>Plasma_cell_IgG      </td><td>1</td><td> 1</td></tr>\n",
       "\t<tr><th scope=row>Plasma_cell_IgM</th><td> 1.8654506</td><td>14.24536</td><td>23.3519016</td><td>6.738117e-06</td><td>6.064305e-05</td><td>Plasma_cell_IgM      </td><td>1</td><td> 1</td></tr>\n",
       "\t<tr><th scope=row>Plasmablast</th><td> 1.5772031</td><td>15.68467</td><td>13.0170670</td><td>5.464777e-04</td><td>2.459150e-03</td><td>Plasmablast          </td><td>1</td><td> 1</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 9 × 8\n",
       "\\begin{tabular}{r|llllllll}\n",
       "  & logFC & logCPM & F & PValue & FDR & CellType & Sig & Diff\\\\\n",
       "  & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <chr> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\tB\\_exhausted &  0.1616904 & 15.45495 &  0.2760204 & 6.008314e-01 & 6.008314e-01 & B\\_exhausted           & 0 &  0\\\\\n",
       "\tB\\_immature & -0.3098771 & 16.07439 &  1.4774286 & 2.278914e-01 & 3.418371e-01 & B\\_immature            & 0 &  0\\\\\n",
       "\tB\\_naive & -0.2113975 & 19.17061 &  4.3488903 & 4.034672e-02 & 7.262410e-02 & B\\_naive               & 1 & -1\\\\\n",
       "\tB\\_non-switched\\_memory & -0.1353764 & 15.53814 &  0.3330929 & 5.655274e-01 & 6.008314e-01 & B\\_non-switched\\_memory & 0 &  0\\\\\n",
       "\tB\\_switched\\_memory & -0.1482451 & 16.59030 &  0.4909614 & 4.856110e-01 & 6.008314e-01 & B\\_switched\\_memory     & 0 &  0\\\\\n",
       "\tPlasma\\_cell\\_IgA &  0.7192348 & 15.46497 &  4.8448680 & 3.072471e-02 & 6.913061e-02 & Plasma\\_cell\\_IgA       & 1 &  1\\\\\n",
       "\tPlasma\\_cell\\_IgG &  1.2398394 & 15.66264 & 11.7834434 & 9.646690e-04 & 2.894007e-03 & Plasma\\_cell\\_IgG       & 1 &  1\\\\\n",
       "\tPlasma\\_cell\\_IgM &  1.8654506 & 14.24536 & 23.3519016 & 6.738117e-06 & 6.064305e-05 & Plasma\\_cell\\_IgM       & 1 &  1\\\\\n",
       "\tPlasmablast &  1.5772031 & 15.68467 & 13.0170670 & 5.464777e-04 & 2.459150e-03 & Plasmablast           & 1 &  1\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 9 × 8\n",
       "\n",
       "| <!--/--> | logFC &lt;dbl&gt; | logCPM &lt;dbl&gt; | F &lt;dbl&gt; | PValue &lt;dbl&gt; | FDR &lt;dbl&gt; | CellType &lt;chr&gt; | Sig &lt;dbl&gt; | Diff &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|\n",
       "| B_exhausted |  0.1616904 | 15.45495 |  0.2760204 | 6.008314e-01 | 6.008314e-01 | B_exhausted           | 0 |  0 |\n",
       "| B_immature | -0.3098771 | 16.07439 |  1.4774286 | 2.278914e-01 | 3.418371e-01 | B_immature            | 0 |  0 |\n",
       "| B_naive | -0.2113975 | 19.17061 |  4.3488903 | 4.034672e-02 | 7.262410e-02 | B_naive               | 1 | -1 |\n",
       "| B_non-switched_memory | -0.1353764 | 15.53814 |  0.3330929 | 5.655274e-01 | 6.008314e-01 | B_non-switched_memory | 0 |  0 |\n",
       "| B_switched_memory | -0.1482451 | 16.59030 |  0.4909614 | 4.856110e-01 | 6.008314e-01 | B_switched_memory     | 0 |  0 |\n",
       "| Plasma_cell_IgA |  0.7192348 | 15.46497 |  4.8448680 | 3.072471e-02 | 6.913061e-02 | Plasma_cell_IgA       | 1 |  1 |\n",
       "| Plasma_cell_IgG |  1.2398394 | 15.66264 | 11.7834434 | 9.646690e-04 | 2.894007e-03 | Plasma_cell_IgG       | 1 |  1 |\n",
       "| Plasma_cell_IgM |  1.8654506 | 14.24536 | 23.3519016 | 6.738117e-06 | 6.064305e-05 | Plasma_cell_IgM       | 1 |  1 |\n",
       "| Plasmablast |  1.5772031 | 15.68467 | 13.0170670 | 5.464777e-04 | 2.459150e-03 | Plasmablast           | 1 |  1 |\n",
       "\n"
      ],
      "text/plain": [
       "                      logFC      logCPM   F          PValue       FDR         \n",
       "B_exhausted            0.1616904 15.45495  0.2760204 6.008314e-01 6.008314e-01\n",
       "B_immature            -0.3098771 16.07439  1.4774286 2.278914e-01 3.418371e-01\n",
       "B_naive               -0.2113975 19.17061  4.3488903 4.034672e-02 7.262410e-02\n",
       "B_non-switched_memory -0.1353764 15.53814  0.3330929 5.655274e-01 6.008314e-01\n",
       "B_switched_memory     -0.1482451 16.59030  0.4909614 4.856110e-01 6.008314e-01\n",
       "Plasma_cell_IgA        0.7192348 15.46497  4.8448680 3.072471e-02 6.913061e-02\n",
       "Plasma_cell_IgG        1.2398394 15.66264 11.7834434 9.646690e-04 2.894007e-03\n",
       "Plasma_cell_IgM        1.8654506 14.24536 23.3519016 6.738117e-06 6.064305e-05\n",
       "Plasmablast            1.5772031 15.68467 13.0170670 5.464777e-04 2.459150e-03\n",
       "                      CellType              Sig Diff\n",
       "B_exhausted           B_exhausted           0    0  \n",
       "B_immature            B_immature            0    0  \n",
       "B_naive               B_naive               1   -1  \n",
       "B_non-switched_memory B_non-switched_memory 0    0  \n",
       "B_switched_memory     B_switched_memory     0    0  \n",
       "Plasma_cell_IgA       Plasma_cell_IgA       1    1  \n",
       "Plasma_cell_IgG       Plasma_cell_IgG       1    1  \n",
       "Plasma_cell_IgM       Plasma_cell_IgM       1    1  \n",
       "Plasmablast           Plasmablast           1    1  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Quadratic changes.\n",
    "bcell.quad.res <- as.data.frame(topTags(glmQLFTest(bcell.linear.fit, coef=5), sort.by='none', n=Inf))\n",
    "bcell.quad.res$CellType <- rownames(bcell.quad.res)\n",
    "bcell.quad.res$Sig <- as.numeric(bcell.quad.res$FDR < 0.1 & bcell.quad.res$PValue < 0.05)\n",
    "bcell.quad.res$Diff <- sign(bcell.quad.res$logFC)\n",
    "bcell.quad.res$Diff[bcell.quad.res$FDR >= 0.1 | bcell.quad.res$PValue >= 0.05] <- 0\n",
    "bcell.quad.res"
   ]
  },
  {
   "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": "4.0.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
