@@ -250,7 +250,7 @@ In our example, we want to find the genes that change specifically in T cells. F

 ```{r do_pseudobulk}

 # The contrast argument specifies what we want to compare

 # We test the expression difference of stimulated and control T-cells

-de_res <- test_de(fit, contrast = fact(cell = "CD4 T cells", stim = "ctrl") - fact(cell = "CD4 T cells", stim = "stim")) 

+de_res <- test_de(fit, contrast = cond(cell = "CD4 T cells", stim = "ctrl") - cond(cell = "CD4 T cells", stim = "stim")) 

 # Most different genes

 head(de_res[order(de_res$pval), ])

@@ -222,10 +222,17 @@ counts(sce_subset) <- as.matrix(counts(sce_subset))

 sce_subset$cell <- droplevels(sce_subset$cell)

 ```

+In the first step we will aggregate the counts of each patient, condition and cell type and form pseudobulk samples. This ensures that I get reliable p-value by treating each patient as a replicate and not each cell.

+

+```{r}

+sce_reduced <- pseudobulk(sce_subset, group_by = vars(ind, stim, cell))

+```

+

+

 We will identify which genes in CD4 positive T-cells are changed most by the treatment. We will fit a full model including the interaction term `stim:cell`. The interaction term will help us identify cell type specific responses to the treatment:

 ```{r kang_fit}

-fit <- glm_gp(sce_subset, design = ~ cell + stim +  stim:cell - 1,

+fit <- glm_gp(sce_reduced, design = ~ cell + stim +  stim:cell - 1,

               reference_level = "NK cells")

 summary(fit)

 ```

@@ -238,26 +245,18 @@ colnames(fit$Beta)

 ```

-In our example, we want to find the genes that change specifically in T cells. Finding cell type specific responses to a treatment is a big advantage of single cell data over bulk data. To get a proper estimate of the uncertainty (cells from the same donor are **not** independent replicates), we create a pseudobulk for each sample:

+In our example, we want to find the genes that change specifically in T cells. Finding cell type specific responses to a treatment is a big advantage of single cell data over bulk data. 

 ```{r do_pseudobulk}

 # The contrast argument specifies what we want to compare

 # We test the expression difference of stimulated and control T-cells

-#

-# There is no sample label in the colData, so we create it on the fly

-# from `stim` and `ind` columns in colData(fit$data).

-de_res <- test_de(fit, contrast = `stimstim` + `cellCD4 T cells:stimstim`, 

-                  pseudobulk_by = paste0(stim, "-", ind)) 

-

-# The large `lfc` values come from groups were nearly all counts are 0

-# Setting them to Inf makes the plots look nicer

-de_res$lfc <- ifelse(abs(de_res$lfc) > 20, sign(de_res$lfc) * Inf, de_res$lfc)

+de_res <- test_de(fit, contrast = fact(cell = "CD4 T cells", stim = "ctrl") - fact(cell = "CD4 T cells", stim = "stim")) 

 # Most different genes

 head(de_res[order(de_res$pval), ])

 ```

-The test is successful and we identify interesting genes that are differentially expressed in interferon-stimulated T cells:  _IFI6_, _IFIT3_, and _IRF7_ literally stand for _Interferon Induced/Regulated Protein_.

+The test is successful and we identify interesting genes that are differentially expressed in interferon-stimulated T cells:  _IFI6_, _IFIT3_ and _ISG15_ literally stand for _Interferon Induced/Stimulated Protein_.

 To get a more complete overview of the results, we can make a volcano plot that compares the log2-fold change (LFC) vs the logarithmized p-values.

@@ -271,17 +270,19 @@ ggplot(de_res, aes(x = lfc, y = -log10(pval))) +

 Another important task in single cell data analysis is the identification of marker genes for cell clusters. For this we can also use our Gamma-Poisson fit. 

-Let's assume we want to find genes that differ between T cells and the B cells. We can directly compare the corresponding coefficients and find genes that differ in the control condition:

+Let's assume we want to find genes that differ between T cells and the B cells. We can directly compare the corresponding coefficients and find genes that differ in the control condition (this time not accounting for the pseudo-replication structure):

 ```{r find_marker_genes}

-marker_genes <- test_de(fit, `cellCD4 T cells` - `cellB cells`, sort_by = pval)

+fit_full <- glm_gp(sce_subset, design = ~ cell + stim +  stim:cell - 1,

+                   reference_level = "NK cells")

+marker_genes <- test_de(fit_full, `cellCD4 T cells` - `cellB cells`, sort_by = pval)

 head(marker_genes)

 ```

 If we want find genes that differ in the stimulated condition, we just include the additional coefficients in the contrast:

 ```{r do_signif_genes}

-marker_genes2 <- test_de(fit, (`cellCD4 T cells` + `cellCD4 T cells:stimstim`) - 

+marker_genes2 <- test_de(fit_full, (`cellCD4 T cells` + `cellCD4 T cells:stimstim`) - 

                                (`cellB cells` + `cellB cells:stimstim`), 

                         sort_by = pval)

@@ -4,7 +4,7 @@ output: github_document

 <!-- README.md is generated from README.Rmd. Please edit that file -->

-```{r, include = FALSE}

+```{r init_chunk, include = FALSE}

 knitr::opts_chunk$set(

   collapse = TRUE,

   comment = "#>",

@@ -59,14 +59,14 @@ If you use this package in a scientific publication, please cite:

 Load the glmGamPoi package

-```{r}

+```{r load_glmGamPoi}

 library(glmGamPoi)

 ```

 To fit a single Gamma-Poisson GLM do:

-```{r}

+```{r single_gp_fit}

 # overdispersion = 1/size

 counts <- rnbinom(n = 10, mu = 5, size = 1/0.7)

@@ -83,12 +83,12 @@ The `glm_gp()` function returns a list with the results of the fit. Most importa

 Fitting repeated Gamma-Poisson GLMs for each gene of a single cell dataset is just as easy:

 I will first load an example dataset using the `TENxPBMCData` package. The dataset has 33,000 genes and 4340 cells. It takes roughly 1.5 minutes to fit the Gamma-Poisson model on the full dataset. For demonstration purposes, I will subset the dataset to 300 genes, but keep the 4340 cells:

-```{r, warning=FALSE, message = FALSE}

+```{r load_additional_packages, warning=FALSE, message = FALSE}

 library(SummarizedExperiment)

 library(DelayedMatrixStats)

 ```

-```{r}

+```{r load_pbmc_data}

 # The full dataset with 33,000 genes and 4340 cells

 # The first time this is run, it will download the data

 pbmcs <- TENxPBMCData::TENxPBMCData("pbmc4k")

@@ -101,7 +101,7 @@ pbmcs_subset

 I call `glm_gp()` to fit one GLM model for each gene and force the calculation to happen in memory.

-```{r}

+```{r simple_fit}

 fit <- glm_gp(pbmcs_subset, on_disk = FALSE)

 summary(fit)

 ```

@@ -114,7 +114,7 @@ I compare my method (in-memory and on-disk) with `r BiocStyle::Biocpkg("DESeq2")

-```{r, warning=FALSE}

+```{r run_bench_mark, warning=FALSE}

 # Explicitly realize count matrix in memory so that it is a fair comparison

 pbmcs_subset <- as.matrix(assay(pbmcs_subset))

 model_matrix <- matrix(1, nrow = ncol(pbmcs_subset))

@@ -145,7 +145,7 @@ On this dataset, `glmGamPoi` is more than 5 times faster than `edgeR` and more t

 To demonstrate that the method does not sacrifice accuracy, I compare the parameters that each method estimates. The means and β coefficients are identical, but that the overdispersion estimates from `glmGamPoi` are more reliable:

-```{r message=FALSE, warning=FALSE}

+```{r compare_with_deseq2_and_edger, message=FALSE, warning=FALSE}

 # Results with my method

 fit <- glm_gp(pbmcs_subset, design = model_matrix, on_disk = FALSE)

@@ -201,14 +201,14 @@ The method scales linearly, with the number of rows and columns in the dataset.

 `glmGamPoi` provides an interface to do quasi-likelihood ratio testing to identify differentially expressed genes. To demonstrate this feature, we will use the data from [Kang _et al._ (2018)](https://www.ncbi.nlm.nih.gov/pubmed/29227470) provided by the `MuscData` package. This is a single cell dataset of 8 Lupus patients for which 10x droplet-based scRNA-seq was performed before and after treatment with interferon beta. The `SingleCellExperiment` object conveniently provides the patient id (`ind`), treatment status (`stim`) and cell type (`cell`):

-```{r}

+```{r load_kang_data}

 sce <- muscData::Kang18_8vs8()

 colData(sce)

 ```

 For demonstration purpose, I will work on a subset of the genes and cells:

-```{r}

+```{r subset_kang_data}

 set.seed(1)

 # Take highly expressed genes and proper cells:

 sce_subset <- sce[rowSums(counts(sce)) > 100, 

@@ -224,7 +224,7 @@ sce_subset$cell <- droplevels(sce_subset$cell)

 We will identify which genes in CD4 positive T-cells are changed most by the treatment. We will fit a full model including the interaction term `stim:cell`. The interaction term will help us identify cell type specific responses to the treatment:

-```{r}

+```{r kang_fit}

 fit <- glm_gp(sce_subset, design = ~ cell + stim +  stim:cell - 1,

               reference_level = "NK cells")

 summary(fit)

@@ -233,14 +233,14 @@ summary(fit)

 To see how the coefficient of our model are called, we look at the `colnames(fit$Beta)`:

-```{r}

+```{r show_beta_colnames}

 colnames(fit$Beta)

 ```

 In our example, we want to find the genes that change specifically in T cells. Finding cell type specific responses to a treatment is a big advantage of single cell data over bulk data. To get a proper estimate of the uncertainty (cells from the same donor are **not** independent replicates), we create a pseudobulk for each sample:

-```{r}

+```{r do_pseudobulk}

 # The contrast argument specifies what we want to compare

 # We test the expression difference of stimulated and control T-cells

 #

@@ -261,7 +261,7 @@ The test is successful and we identify interesting genes that are differentially

 To get a more complete overview of the results, we can make a volcano plot that compares the log2-fold change (LFC) vs the logarithmized p-values.

-```{r}

+```{r make_volcano_plot}

 library(ggplot2)

 ggplot(de_res, aes(x = lfc, y = -log10(pval))) +

   geom_point(size = 0.6, aes(color = adj_pval < 0.1)) +

@@ -273,14 +273,14 @@ Another important task in single cell data analysis is the identification of mar

 Let's assume we want to find genes that differ between T cells and the B cells. We can directly compare the corresponding coefficients and find genes that differ in the control condition:

-```{r}

+```{r find_marker_genes}

 marker_genes <- test_de(fit, `cellCD4 T cells` - `cellB cells`, sort_by = pval)

 head(marker_genes)

 ```

 If we want find genes that differ in the stimulated condition, we just include the additional coefficients in the contrast:

-```{r}

+```{r do_signif_genes}

 marker_genes2 <- test_de(fit, (`cellCD4 T cells` + `cellCD4 T cells:stimstim`) - 

                                (`cellB cells` + `cellB cells:stimstim`), 

                         sort_by = pval)

@@ -292,7 +292,7 @@ We identify many genes related to the human leukocyte antigen (HLA) system that

 A note of caution: applying `test_de()` to single cell data without the pseudobulk gives overly optimistic p-values. This is due to the fact that cells from the same sample are not independent replicates! It can still be fine to use the method for identifying marker genes, as long as one is aware of the difficulties interpreting the results.

-```{r}

+```{r plot_marker_genes}

 # Create a data.frame with the expression values, gene names, and cell types

 tmp <- data.frame(gene = rep(marker_genes$name[1:6], times = ncol(sce_subset)),

                   expression = c(counts(sce_subset)[marker_genes$name[1:6], ]),

@@ -316,7 +316,7 @@ This work was supported by the EMBL International PhD Programme and the European

 # Session Info

-```{r}

+```{r print_sessionInfo}

 sessionInfo()

 ```

@@ -307,7 +307,9 @@ ggplot(tmp, aes(x = celltype, y = expression)) +

+# Acknowlegments

+This work was supported by the EMBL International PhD Programme and the European Research Council Synergy grant DECODE under grant agreement No. 810296.

@@ -53,7 +53,7 @@ If you use this package in a scientific publication, please cite:

 > glmGamPoi: Fitting Gamma-Poisson Generalized Linear Models on Single Cell Count Data  

 > Constantin Ahlmann-Eltze, Wolfgang Huber  

-> bioRxiv 2020.08.13.249623; doi: https://doi.org/10.1101/2020.08.13.249623

+> Bioinformatics; 2020-12-09; doi: https://doi.org/10.1093/bioinformatics/btaa1009

 # Example

@@ -245,13 +245,13 @@ In our example, we want to find the genes that change specifically in T cells. F

 # We test the expression difference of stimulated and control T-cells

 #

 # There is no sample label in the colData, so we create it on the fly

-# (paste0(stim, "-", ind)) to aggregate the counts by sample

+# from `stim` and `ind` columns in colData(fit$data).

 de_res <- test_de(fit, contrast = `stimstim` + `cellCD4 T cells:stimstim`, 

                   pseudobulk_by = paste0(stim, "-", ind)) 

 # The large `lfc` values come from groups were nearly all counts are 0

 # Setting them to Inf makes the plots look nicer

-de_res$lfc = ifelse(abs(de_res$lfc) > 20, sign(de_res$lfc) * Inf, de_res$lfc)

+de_res$lfc <- ifelse(abs(de_res$lfc) > 20, sign(de_res$lfc) * Inf, de_res$lfc)

 # Most different genes

 head(de_res[order(de_res$pval), ])

@@ -265,7 +265,7 @@ To get a more complete overview of the results, we can make a volcano plot that

 library(ggplot2)

 ggplot(de_res, aes(x = lfc, y = -log10(pval))) +

   geom_point(size = 0.6, aes(color = adj_pval < 0.1)) +

-  ggtitle("Volcano Plot", "Genes that change most through\ninterferon-beta treatment in T cells")

+  ggtitle("Volcano Plot", "Genes that change most through interferon-beta treatment in T cells")

 ```

@@ -288,7 +288,7 @@ marker_genes2 <- test_de(fit, (`cellCD4 T cells` + `cellCD4 T cells:stimstim`) -

 head(marker_genes2)

 ```

-Unsurprisingly, we identify many genes related to the human leukocyte antigen (HLA) system that is important for antigen presenting cells like B-cells, but are not expressed by T helper cells. The plot below shows the expression differences.

+We identify many genes related to the human leukocyte antigen (HLA) system that is important for antigen presenting cells like B-cells, but are not expressed by T helper cells. The plot below shows the expression differences.

 A note of caution: applying `test_de()` to single cell data without the pseudobulk gives overly optimistic p-values. This is due to the fact that cells from the same sample are not independent replicates! It can still be fine to use the method for identifying marker genes, as long as one is aware of the difficulties interpreting the results.

@@ -199,58 +199,118 @@ The method scales linearly, with the number of rows and columns in the dataset.

 ## Differential expression analysis

-`glmGamPoi` provides an interface to do quasi-likelihood ratio testing to identify differentially expressed genes:

+`glmGamPoi` provides an interface to do quasi-likelihood ratio testing to identify differentially expressed genes. To demonstrate this feature, we will use the data from [Kang _et al._ (2018)](https://www.ncbi.nlm.nih.gov/pubmed/29227470) provided by the `MuscData` package. This is a single cell dataset of 8 Lupus patients for which 10x droplet-based scRNA-seq was performed before and after treatment with interferon beta. The `SingleCellExperiment` object conveniently provides the patient id (`ind`), treatment status (`stim`) and cell type (`cell`):

 ```{r}

-# Create random categorical assignment to demonstrate DE

-group <- sample(c("Group1", "Group2"), size = ncol(pbmcs_subset), replace = TRUE)

-

-# Fit model with group vector as design

-fit <- glm_gp(pbmcs_subset, design = group)

-# Compare against model without group 

-res <- test_de(fit, reduced_design = ~ 1)

-# Look at first 6 genes

-head(res)

+sce <- muscData::Kang18_8vs8()

+colData(sce)

 ```

-The p-values agree well with the ones that `edgeR` is calculating. This is because `glmGamPoi` uses the same framework of quasi-likelihood ratio tests that was invented by `edgeR` and is described in [Lund et al. (2012)](https://doi.org/10.1515/1544-6115.1826).

+For demonstration purpose, I will work on a subset of the genes and cells:

-```{r fig.height=3, fig.width=3, warning=FALSE, message=FALSE}

-model_matrix <- model.matrix(~ group, data = data.frame(group = group))

-edgeR_data <- edgeR::DGEList(pbmcs_subset)

-edgeR_data <- edgeR::calcNormFactors(edgeR_data)

-edgeR_data <- edgeR::estimateDisp(edgeR_data, design = model_matrix)

-edgeR_fit <- edgeR::glmQLFit(edgeR_data, design = model_matrix)

-edgeR_test <- edgeR::glmQLFTest(edgeR_fit, coef = 2)

-edgeR_res <- edgeR::topTags(edgeR_test, sort.by = "none", n = nrow(pbmcs_subset))

+```{r}

+set.seed(1)

+# Take highly expressed genes and proper cells:

+sce_subset <- sce[rowSums(counts(sce)) > 100, 

+                  sample(which(sce$multiplets == "singlet" & 

+                              ! is.na(sce$cell) &

+                              sce$cell %in% c("CD4 T cells", "B cells", "NK cells")), 

+                         1000)]

+# Convert counts to dense matrix

+counts(sce_subset) <- as.matrix(counts(sce_subset))

+# Remove empty levels because glm_gp() will complain otherwise

+sce_subset$cell <- droplevels(sce_subset$cell)

 ```

+We will identify which genes in CD4 positive T-cells are changed most by the treatment. We will fit a full model including the interaction term `stim:cell`. The interaction term will help us identify cell type specific responses to the treatment:

-```{r fig.height=4, fig.width=3, warning=FALSE, message=FALSE, echo = FALSE}

-par(cex.main = 2, cex.lab = 1.5)

-plot(res$pval, edgeR_res$table$PValue, pch = 16, log = "xy",

-     main = "p-values", xlab = "glmGamPoi", ylab = "edgeR")

-abline(0,1)

+```{r}

+fit <- glm_gp(sce_subset, design = ~ cell + stim +  stim:cell - 1,

+              reference_level = "NK cells")

+summary(fit)

+```

+

+

+To see how the coefficient of our model are called, we look at the `colnames(fit$Beta)`:

+

+```{r}

+colnames(fit$Beta)

+```

+

+

+In our example, we want to find the genes that change specifically in T cells. Finding cell type specific responses to a treatment is a big advantage of single cell data over bulk data. To get a proper estimate of the uncertainty (cells from the same donor are **not** independent replicates), we create a pseudobulk for each sample:

+

+```{r}

+# The contrast argument specifies what we want to compare

+# We test the expression difference of stimulated and control T-cells

+#

+# There is no sample label in the colData, so we create it on the fly

+# (paste0(stim, "-", ind)) to aggregate the counts by sample

+de_res <- test_de(fit, contrast = `stimstim` + `cellCD4 T cells:stimstim`, 

+                  pseudobulk_by = paste0(stim, "-", ind)) 

+

+# The large `lfc` values come from groups were nearly all counts are 0

+# Setting them to Inf makes the plots look nicer

+de_res$lfc = ifelse(abs(de_res$lfc) > 20, sign(de_res$lfc) * Inf, de_res$lfc)

+

+# Most different genes

+head(de_res[order(de_res$pval), ])

+```

+

+The test is successful and we identify interesting genes that are differentially expressed in interferon-stimulated T cells:  _IFI6_, _IFIT3_, and _IRF7_ literally stand for _Interferon Induced/Regulated Protein_.

+

+To get a more complete overview of the results, we can make a volcano plot that compares the log2-fold change (LFC) vs the logarithmized p-values.

+

+```{r}

+library(ggplot2)

+ggplot(de_res, aes(x = lfc, y = -log10(pval))) +

+  geom_point(size = 0.6, aes(color = adj_pval < 0.1)) +

+  ggtitle("Volcano Plot", "Genes that change most through\ninterferon-beta treatment in T cells")

 ```

-#### Pseudobulk

+Another important task in single cell data analysis is the identification of marker genes for cell clusters. For this we can also use our Gamma-Poisson fit. 

+

+Let's assume we want to find genes that differ between T cells and the B cells. We can directly compare the corresponding coefficients and find genes that differ in the control condition:

-Be very careful how you interpret the p-values of a single cell experiment. Cells that come from one individual are not independent replicates. That means that you cannot turn your RNA-seq experiment with 3 treated and 3 control samples into a 3000 vs 3000 experiment by measuring 1000 cells per sample. The actual unit of replication are still the 3 samples in each condition. 

+```{r}

+marker_genes <- test_de(fit, `cellCD4 T cells` - `cellB cells`, sort_by = pval)

+head(marker_genes)

+```

-Nonetheless, single cell data is valuable because it allows you to compare the effect of a treatment on specific cell types. The simplest way to do such a test is called pseudobulk. This means that the data is subset to the cells of a specific cell type. Then the counts of cells from the same sample are combined to form a "pseudobulk" sample. The `test_de()` function of glmGamPoi supports this feature directly through the `pseudobulk_by` and `subset_to` parameters:

+If we want find genes that differ in the stimulated condition, we just include the additional coefficients in the contrast:

 ```{r}

-# say we have cell type labels for each cell and know from which sample they come originally

-sample_labels <- rep(paste0("sample_", 1:6), length = ncol(pbmcs_subset))

-cell_type_labels <- sample(c("T-cells", "B-cells", "Macrophages"), ncol(pbmcs_subset), replace = TRUE)

-

-test_de(fit, contrast = Group1 - Group2,

-        pseudobulk_by = sample_labels, 

-        subset_to = cell_type_labels == "T-cells",

-        n_max = 4, sort_by = pval, decreasing = FALSE)

+marker_genes2 <- test_de(fit, (`cellCD4 T cells` + `cellCD4 T cells:stimstim`) - 

+                               (`cellB cells` + `cellB cells:stimstim`), 

+                        sort_by = pval)

+

+head(marker_genes2)

 ```

+Unsurprisingly, we identify many genes related to the human leukocyte antigen (HLA) system that is important for antigen presenting cells like B-cells, but are not expressed by T helper cells. The plot below shows the expression differences.

+

+A note of caution: applying `test_de()` to single cell data without the pseudobulk gives overly optimistic p-values. This is due to the fact that cells from the same sample are not independent replicates! It can still be fine to use the method for identifying marker genes, as long as one is aware of the difficulties interpreting the results.

+

+```{r}

+# Create a data.frame with the expression values, gene names, and cell types

+tmp <- data.frame(gene = rep(marker_genes$name[1:6], times = ncol(sce_subset)),

+                  expression = c(counts(sce_subset)[marker_genes$name[1:6], ]),

+                  celltype = rep(sce_subset$cell, each = 6))

+

+ggplot(tmp, aes(x = celltype, y = expression)) +

+  geom_jitter(height = 0.1) +

+  stat_summary(geom = "crossbar", fun = "mean", color = "red") +

+  facet_wrap(~ gene, scales = "free_y") +

+  ggtitle("Marker genes of B vs. T cells")

+```

+

+

+

+

+

+

+

 # Session Info

@@ -49,6 +49,11 @@ BiocManager::install("glmGamPoi")

 For the latest developments, see the `r BiocStyle::Githubpkg("const-ae/glmGamPoi", "GitHub")` repo.

+If you use this package in a scientific publication, please cite:

+

+> glmGamPoi: Fitting Gamma-Poisson Generalized Linear Models on Single Cell Count Data  

+> Constantin Ahlmann-Eltze, Wolfgang Huber  

+> bioRxiv 2020.08.13.249623; doi: https://doi.org/10.1101/2020.08.13.249623

 # Example

@@ -24,12 +24,12 @@ set.seed(2)

 Pronounciation: [`dʒi əl əm ɡam ˈpwɑ`](http://ipa-reader.xyz/?text=d%CA%92i%20%C9%99l%20%C9%99m%20%C9%A1am%20%CB%88pw%C9%91)

-The core design aims of `gmlGamPoi` are:

+The core design aims of `glmGamPoi` are:

-* Fit the Gamma-Poisson models on arbitrarily large or small datasets

+* Fit Gamma-Poisson models on arbitrarily large or small datasets

 * Be faster than alternative methods, such as `DESeq2` or `edgeR`

 * Calculate exact or approximate results based on user preference

-* Support in memory or on-disk data

+* Support in-memory or on-disk data

 * Follow established conventions around tools for RNA-seq analysis

 * Present a simple user-interface

 * Avoid unnecessary dependencies

@@ -64,6 +64,7 @@ To fit a single Gamma-Poisson GLM do:

 ```{r}

 # overdispersion = 1/size

 counts <- rnbinom(n = 10, mu = 5, size = 1/0.7)

+

 # design = ~ 1 means that an intercept-only model is fit

 fit <- glm_gp(counts, design = ~ 1)

 fit

@@ -86,6 +87,7 @@ library(DelayedMatrixStats)

 # The full dataset with 33,000 genes and 4340 cells

 # The first time this is run, it will download the data

 pbmcs <- TENxPBMCData::TENxPBMCData("pbmc4k")

+

 # I want genes where at least some counts are non-zero

 non_empty_rows <- which(rowSums2(assay(pbmcs)) > 0)

 pbmcs_subset <- pbmcs[sample(non_empty_rows, 300), ]

@@ -187,7 +189,7 @@ I am comparing the gene-wise estimates of the coefficients from all three method

 ## Scalability

-The method scales linearly, with the number of rows and columns in the dataset. For example: fitting the full `pbmc4k` dataset with subsampling on a modern MacBook Pro in-memory takes ~1 minute and on-disk a little over 4 minutes. Fitting the `pbmc68k` (17x the size) takes ~73 minutes (17x the time) on-disk. Fitting that dataset in-memory is not possible because it is just too big: the maximum in-memory matrix size is `2^31-1 ≈ 2.1e9` is elements, the `pbmc68k` dataset however has nearly 100 million elements more than that.

+The method scales linearly, with the number of rows and columns in the dataset. For example: fitting the full `pbmc4k` dataset with subsampling on a modern MacBook Pro in-memory takes ~1 minute and on-disk a little over 4 minutes. Fitting the `pbmc68k` (17x the size) takes ~73 minutes (17x the time) on-disk.

 ## Differential expression analysis

@@ -229,7 +231,7 @@ abline(0,1)

 #### Pseudobulk

-Be very careful how you interpret the p-values of a single cell experiment. Cells that come from one individual are not independent replicates. That means that you cannot turn your RNA-seq experiment with 3 treated and 3 control samples into a 3000 vs 3000 experiment through single cell experiment with 1000 cells per sample. The actual unit of replication are still the 3 samples in each condition. 

+Be very careful how you interpret the p-values of a single cell experiment. Cells that come from one individual are not independent replicates. That means that you cannot turn your RNA-seq experiment with 3 treated and 3 control samples into a 3000 vs 3000 experiment by measuring 1000 cells per sample. The actual unit of replication are still the 3 samples in each condition. 

 Nonetheless, single cell data is valuable because it allows you to compare the effect of a treatment on specific cell types. The simplest way to do such a test is called pseudobulk. This means that the data is subset to the cells of a specific cell type. Then the counts of cells from the same sample are combined to form a "pseudobulk" sample. The `test_de()` function of glmGamPoi supports this feature directly through the `pseudobulk_by` and `subset_to` parameters:

@@ -13,15 +13,17 @@ knitr::opts_chunk$set(

 )

 set.seed(2)

 ```

-# glmGamPoi

+

+# glmGamPoi <a href='https://github.com/const-ae/glmGamPoi'><img src='man/figures/logo.svg' align="right" height="139" /></a>

 <!-- badges: start -->

 [![codecov](https://codecov.io/gh/const-ae/glmGamPoi/branch/master/graph/badge.svg)](https://codecov.io/gh/const-ae/glmGamPoi)

 <!-- badges: end -->

-

 > Fit Gamma-Poisson Generalized Linear Models Reliably.

+Pronounciation: [`dʒi əl əm ɡam ˈpwɑ`](http://ipa-reader.xyz/?text=d%CA%92i%20%C9%99l%20%C9%99m%20%C9%A1am%20%CB%88pw%C9%91)

+

 The core design aims of `gmlGamPoi` are:

 * Fit the Gamma-Poisson models on arbitrarily large or small datasets

@@ -132,7 +132,7 @@ bench::mark(

 )

 ```

-On this dataset, `glmGamPoi` is more than 6 times faster than `edgeR` and more than 20 times faster than `DESeq2`. `glmGamPoi` does **not** use approximations to achieve this performance increase. The performance comes from an optimized algorithm for inferring the overdispersion for each gene. It is tuned for datasets typically encountered in single RNA-seq with many samples and many small counts, by avoiding duplicate calculations.

+On this dataset, `glmGamPoi` is more than 5 times faster than `edgeR` and more than 18 times faster than `DESeq2`. `glmGamPoi` does **not** use approximations to achieve this performance increase. The performance comes from an optimized algorithm for inferring the overdispersion for each gene. It is tuned for datasets typically encountered in single RNA-seq with many samples and many small counts, by avoiding duplicate calculations.

 To demonstrate that the method does not sacrifice accuracy, I compare the parameters that each method estimates. The means and β coefficients are identical, but that the overdispersion estimates from `glmGamPoi` are more reliable:

@@ -144,13 +144,12 @@ fit <- glm_gp(pbmcs_subset, design = model_matrix, on_disk = FALSE)

 dds <- DESeq2::DESeqDataSetFromMatrix(pbmcs_subset, 

                         colData = data.frame(name = seq_len(4340)),

                         design = ~ 1)

-dds <- DESeq2::estimateSizeFactors(dds, "poscounts")

+sizeFactors(dds)  <- fit$size_factors

 dds <- DESeq2::estimateDispersions(dds, quiet = TRUE)

 dds <- DESeq2::nbinomWaldTest(dds, minmu = 1e-6)

 #edgeR

-edgeR_data <- edgeR::DGEList(pbmcs_subset)

-edgeR_data <- edgeR::calcNormFactors(edgeR_data)

+edgeR_data <- edgeR::DGEList(pbmcs_subset, lib.size = fit$size_factors)

 edgeR_data <- edgeR::estimateDisp(edgeR_data, model_matrix)

 edgeR_fit <- edgeR::glmFit(edgeR_data, design = model_matrix)

 ```

@@ -165,7 +164,7 @@ plot(fit$Beta[,1], edgeR_fit$unshrunk.coefficients[,1], pch = 16,

      main = "Beta Coefficients", xlab = "glmGamPoi", ylab = "edgeR")

 abline(0,1)

-plot(fit$Mu[,1], rowData(dds)$baseMean, pch = 16, log="xy",

+plot(fit$Mu[,1], assay(dds, "mu")[,1], pch = 16, log="xy",

      main = "Gene Mean", xlab = "glmGamPoi", ylab = "DESeq2")

 abline(0,1)

 plot(fit$Mu[,1], edgeR_fit$fitted.values[,1], pch = 16, log="xy",

@@ -62,7 +62,7 @@ To fit a single Gamma-Poisson GLM do:

 ```{r}

 # overdispersion = 1/size

 counts <- rnbinom(n = 10, mu = 5, size = 1/0.7)

-# size_factors = FALSE, because only a single GLM is fitted

+# design = ~ 1 means that an intercept-only model is fit

 fit <- glm_gp(counts, design = ~ 1)

 fit

@@ -226,6 +226,24 @@ abline(0,1)

 ```

+#### Pseudobulk

+

+Be very careful how you interpret the p-values of a single cell experiment. Cells that come from one individual are not independent replicates. That means that you cannot turn your RNA-seq experiment with 3 treated and 3 control samples into a 3000 vs 3000 experiment through single cell experiment with 1000 cells per sample. The actual unit of replication are still the 3 samples in each condition. 

+

+Nonetheless, single cell data is valuable because it allows you to compare the effect of a treatment on specific cell types. The simplest way to do such a test is called pseudobulk. This means that the data is subset to the cells of a specific cell type. Then the counts of cells from the same sample are combined to form a "pseudobulk" sample. The `test_de()` function of glmGamPoi supports this feature directly through the `pseudobulk_by` and `subset_to` parameters:

+

+```{r}

+# say we have cell type labels for each cell and know from which sample they come originally

+sample_labels <- rep(paste0("sample_", 1:6), length = ncol(pbmcs_subset))

+cell_type_labels <- sample(c("T-cells", "B-cells", "Macrophages"), ncol(pbmcs_subset), replace = TRUE)

+

+test_de(fit, contrast = Group1 - Group2,

+        pseudobulk_by = sample_labels, 

+        subset_to = cell_type_labels == "T-cells",

+        n_max = 4, sort_by = pval, decreasing = FALSE)

+```

+

+

 # Session Info

 ```{r}

@@ -50,13 +50,20 @@ For the latest developments, see the `r BiocStyle::Githubpkg("const-ae/glmGamPoi

 # Example

+Load the glmGamPoi package

+

+```{r}

+library(glmGamPoi)

+```

+

+

 To fit a single Gamma-Poisson GLM do:

 ```{r}

 # overdispersion = 1/size

 counts <- rnbinom(n = 10, mu = 5, size = 1/0.7)

 # size_factors = FALSE, because only a single GLM is fitted

-fit <- glmGamPoi::glm_gp(counts, design = ~ 1)

+fit <- glm_gp(counts, design = ~ 1)

 fit

 # Internally fit is just a list:

@@ -86,7 +93,7 @@ pbmcs_subset

 I call `glm_gp()` to fit one GLM model for each gene and force the calculation to happen in memory.

 ```{r}

-fit <- glmGamPoi::glm_gp(pbmcs_subset, on_disk = FALSE)

+fit <- glm_gp(pbmcs_subset, on_disk = FALSE)

 summary(fit)

 ```

@@ -106,9 +113,9 @@ model_matrix <- matrix(1, nrow = ncol(pbmcs_subset))

 bench::mark(

   glmGamPoi_in_memory = {

-    glmGamPoi::glm_gp(pbmcs_subset, design = model_matrix, on_disk = FALSE)

+    glm_gp(pbmcs_subset, design = model_matrix, on_disk = FALSE)

   }, glmGamPoi_on_disk = {

-    glmGamPoi::glm_gp(pbmcs_subset, design = model_matrix, on_disk = TRUE)

+    glm_gp(pbmcs_subset, design = model_matrix, on_disk = TRUE)

   }, DESeq2 = suppressMessages({

     dds <- DESeq2::DESeqDataSetFromMatrix(pbmcs_subset,

                         colData = data.frame(name = seq_len(4340)),

@@ -131,7 +138,7 @@ To demonstrate that the method does not sacrifice accuracy, I compare the parame

 ```{r message=FALSE, warning=FALSE}

 # Results with my method

-fit <- glmGamPoi::glm_gp(pbmcs_subset, design = model_matrix, on_disk = FALSE)

+fit <- glm_gp(pbmcs_subset, design = model_matrix, on_disk = FALSE)

 # DESeq2

 dds <- DESeq2::DESeqDataSetFromMatrix(pbmcs_subset, 

@@ -191,9 +198,9 @@ The method scales linearly, with the number of rows and columns in the dataset.

 group <- sample(c("Group1", "Group2"), size = ncol(pbmcs_subset), replace = TRUE)

 # Fit model with group vector as design

-fit <- glmGamPoi::glm_gp(pbmcs_subset, design = group)

+fit <- glm_gp(pbmcs_subset, design = group)

 # Compare against model without group 

-res <- glmGamPoi::gampoi_test_qlr(pbmcs_subset, fit, reduced_design = ~ 1)

+res <- test_de(fit, reduced_design = ~ 1)

 # Look at first 6 genes

 head(res)

 ```

@@ -56,11 +56,11 @@ To fit a single Gamma-Poisson GLM do:

 # overdispersion = 1/size

 counts <- rnbinom(n = 10, mu = 5, size = 1/0.7)

 # size_factors = FALSE, because only a single GLM is fitted

-fit <- glmGamPoi::glm_gp(counts, design = ~ 1, size_factors = FALSE, overdispersion_shrinkage = FALSE)

+fit <- glmGamPoi::glm_gp(counts, design = ~ 1)

 fit

 # Internally fit is just a list:

-as.list(fit)

+as.list(fit)[1:2]

 ```

 The `glm_gp()` function returns a list with the results of the fit. Most importantly, it contains the estimates for the coefficients β and the overdispersion.

@@ -127,7 +127,7 @@ bench::mark(

 On this dataset, `glmGamPoi` is more than 6 times faster than `edgeR` and more than 20 times faster than `DESeq2`. `glmGamPoi` does **not** use approximations to achieve this performance increase. The performance comes from an optimized algorithm for inferring the overdispersion for each gene. It is tuned for datasets typically encountered in single RNA-seq with many samples and many small counts, by avoiding duplicate calculations.

-To demonstrate that the method is not sacrificing accuracy, I compare the parameters that each method estimates. I find that means and β coefficients are identical, but that the estimates of the overdispersion estimates from `glmGamPoi` are more reliable:

+To demonstrate that the method does not sacrifice accuracy, I compare the parameters that each method estimates. The means and β coefficients are identical, but that the overdispersion estimates from `glmGamPoi` are more reliable:

 ```{r message=FALSE, warning=FALSE}

 # Results with my method

@@ -182,6 +182,43 @@ I am comparing the gene-wise estimates of the coefficients from all three method

 The method scales linearly, with the number of rows and columns in the dataset. For example: fitting the full `pbmc4k` dataset with subsampling on a modern MacBook Pro in-memory takes ~1 minute and on-disk a little over 4 minutes. Fitting the `pbmc68k` (17x the size) takes ~73 minutes (17x the time) on-disk. Fitting that dataset in-memory is not possible because it is just too big: the maximum in-memory matrix size is `2^31-1 ≈ 2.1e9` is elements, the `pbmc68k` dataset however has nearly 100 million elements more than that.

+## Differential expression analysis

+

+`glmGamPoi` provides an interface to do quasi-likelihood ratio testing to identify differentially expressed genes:

+

+```{r}

+# Create random categorical assignment to demonstrate DE

+group <- sample(c("Group1", "Group2"), size = ncol(pbmcs_subset), replace = TRUE)

+

+# Fit model with group vector as design

+fit <- glmGamPoi::glm_gp(pbmcs_subset, design = group)

+# Compare against model without group 

+res <- glmGamPoi::gampoi_test_qlr(pbmcs_subset, fit, reduced_design = ~ 1)

+# Look at first 6 genes

+head(res)

+```

+

+The p-values agree well with the ones that `edgeR` is calculating. This is because `glmGamPoi` uses the same framework of quasi-likelihood ratio tests that was invented by `edgeR` and is described in [Lund et al. (2012)](https://doi.org/10.1515/1544-6115.1826).

+

+```{r fig.height=3, fig.width=3, warning=FALSE, message=FALSE}

+model_matrix <- model.matrix(~ group, data = data.frame(group = group))

+edgeR_data <- edgeR::DGEList(pbmcs_subset)

+edgeR_data <- edgeR::calcNormFactors(edgeR_data)

+edgeR_data <- edgeR::estimateDisp(edgeR_data, design = model_matrix)

+edgeR_fit <- edgeR::glmQLFit(edgeR_data, design = model_matrix)

+edgeR_test <- edgeR::glmQLFTest(edgeR_fit, coef = 2)

+edgeR_res <- edgeR::topTags(edgeR_test, sort.by = "none", n = nrow(pbmcs_subset))

+```

+

+

+```{r fig.height=4, fig.width=3, warning=FALSE, message=FALSE, echo = FALSE}

+par(cex.main = 2, cex.lab = 1.5)

+plot(res$pval, edgeR_res$table$PValue, pch = 16, log = "xy",

+     main = "p-values", xlab = "glmGamPoi", ylab = "edgeR")

+abline(0,1)

+```

+

+

 # Session Info

 ```{r}

@@ -56,11 +56,11 @@ To fit a single Gamma-Poisson GLM do:

 # overdispersion = 1/size

 counts <- rnbinom(n = 10, mu = 5, size = 1/0.7)

 # size_factors = FALSE, because only a single GLM is fitted

-fit <- glmGamPoi::glm_gp(counts, design = ~ 1, size_factors = FALSE)

+fit <- glmGamPoi::glm_gp(counts, design = ~ 1, size_factors = FALSE, overdispersion_shrinkage = FALSE)

 fit

 # Internally fit is just a list:

-c(fit)

+as.list(fit)

 ```

 The `glm_gp()` function returns a list with the results of the fit. Most importantly, it contains the estimates for the coefficients β and the overdispersion.

@@ -96,8 +96,10 @@ summary(fit)

 I compare my method (in-memory and on-disk) with `r BiocStyle::Biocpkg("DESeq2")` and `r BiocStyle::Biocpkg("edgeR")`. Both are classical methods for analyzing RNA-Seq datasets and have been around for almost 10 years. Note that both tools can do a lot more than just fitting the Gamma-Poisson model, so this benchmark only serves to give a general impression of the performance.

+

+

 ```{r, warning=FALSE}

-# Explicitly realize count matrix in memory

+# Explicitly realize count matrix in memory so that it is a fair comparison

 pbmcs_subset <- as.matrix(assay(pbmcs_subset))

 model_matrix <- matrix(1, nrow = ncol(pbmcs_subset))

@@ -119,7 +121,7 @@ bench::mark(

     edgeR_data <- edgeR::calcNormFactors(edgeR_data)

     edgeR_data <- edgeR::estimateDisp(edgeR_data, model_matrix)

     edgeR_fit <- edgeR::glmFit(edgeR_data, design = model_matrix)

-  }, check = FALSE

+  }, check = FALSE, min_iterations = 3

 )

 ```

@@ -16,7 +16,6 @@ set.seed(2)

 # glmGamPoi

 <!-- badges: start -->

-[![R build status](https://github.com/const-ae/glmGamPoi/workflows/R-CMD-check/badge.svg)](https://github.com/const-ae/glmGamPoi/actions)

 [![codecov](https://codecov.io/gh/const-ae/glmGamPoi/branch/master/graph/badge.svg)](https://codecov.io/gh/const-ae/glmGamPoi)

 <!-- badges: end -->

@@ -37,7 +37,7 @@ The core design aims of `gmlGamPoi` are:

 # Installation

-You can install the release version of `r BiocStyle::Biocpkg("glmGamPoi")` from BioConductor:

+You can install the release version of *[glmGamPoi](https://bioconductor.org/packages/glmGamPoi)* from BioConductor:

 ``` r

 if (!requireNamespace("BiocManager", quietly = TRUE))

@@ -11,42 +11,180 @@ knitr::opts_chunk$set(

   fig.path = "man/figures/README-"

   # out.width = "100%"

 )

-set.seed(1)

+set.seed(2)

 ```

 # glmGamPoi

 <!-- badges: start -->

-[![R build status](https://github.com/const-ae/glmGamPoi/workflows/R-CMD-check/badge.svg)](https://github.com/const-ae/glmGamPoi)

+[![R build status](https://github.com/const-ae/glmGamPoi/workflows/R-CMD-check/badge.svg)](https://github.com/const-ae/glmGamPoi/actions)

 [![codecov](https://codecov.io/gh/const-ae/glmGamPoi/branch/master/graph/badge.svg)](https://codecov.io/gh/const-ae/glmGamPoi)

 <!-- badges: end -->

-```{r child = "man/fragments/paragraph-intro.Rmd"}

-```

+> Fit Gamma-Poisson Generalized Linear Models Reliably.

+

+The core design aims of `gmlGamPoi` are:

+

+* Fit the Gamma-Poisson models on arbitrarily large or small datasets

+* Be faster than alternative methods, such as `DESeq2` or `edgeR`

+* Calculate exact or approximate results based on user preference

+* Support in memory or on-disk data

+* Follow established conventions around tools for RNA-seq analysis

+* Present a simple user-interface

+* Avoid unnecessary dependencies

+* Make integration into other tools easy

+

+

+# Installation

+You can install the release version of `r BiocStyle::Biocpkg("glmGamPoi")` from BioConductor:

-## Installation

+``` r

+if (!requireNamespace("BiocManager", quietly = TRUE))

+    install.packages("BiocManager")

-```{r child = "man/fragments/paragraph-installation.Rmd"}

+BiocManager::install("glmGamPoi")

 ```

+For the latest developments, see the `r BiocStyle::Githubpkg("const-ae/glmGamPoi", "GitHub")` repo.

-## Example

-```{r child = "man/fragments/paragraph-simple_example.Rmd"}

+# Example

+

+To fit a single Gamma-Poisson GLM do:

+

+```{r}

+# overdispersion = 1/size

+counts <- rnbinom(n = 10, mu = 5, size = 1/0.7)

+# size_factors = FALSE, because only a single GLM is fitted

+fit <- glmGamPoi::glm_gp(counts, design = ~ 1, size_factors = FALSE)

+fit

+

+# Internally fit is just a list:

+c(fit)

 ```

+The `glm_gp()` function returns a list with the results of the fit. Most importantly, it contains the estimates for the coefficients β and the overdispersion.

 Fitting repeated Gamma-Poisson GLMs for each gene of a single cell dataset is just as easy:

-```{r child = "man/fragments/paragraph-pbmc4k_example.Rmd"}

+I will first load an example dataset using the `TENxPBMCData` package. The dataset has 33,000 genes and 4340 cells. It takes roughly 1.5 minutes to fit the Gamma-Poisson model on the full dataset. For demonstration purposes, I will subset the dataset to 300 genes, but keep the 4340 cells:

+```{r, warning=FALSE, message = FALSE}

+library(SummarizedExperiment)

+library(DelayedMatrixStats)

+```

+

+```{r}

+# The full dataset with 33,000 genes and 4340 cells

+# The first time this is run, it will download the data