pbmc2_media

Back to home

suppressPackageStartupMessages({
  library(CellBench)
  library(scater)
  library(CellMixS)
  library(variancePartition)
  library(purrr)
  library(jcolors)
  library(here)
  library(tidyr)
  library(dplyr)
  library(stringr)
  library(ComplexHeatmap)
  #library(ggtern)
  library(gridExtra)
  library(scran)
  library(cowplot)
  library(CAMERA)
  library(ggrepel)
  library(readr)
})

options(bitmapType='cairo')

Dataset and parameter

sce <- readRDS(params$data)

param <- readRDS(params$param)
celltype <- param[["celltype"]]
batch <- param[["batch"]]
sample <- param[["sample"]]
dataset_name <- param[["dataset_name"]]
dataset_name
## [1] "pbmc2_media"
n_genes <- nrow(sce)

table(colData(sce)[,celltype])
## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12   13 
## 4257 4118 3134 2450 2286 2108 1994 1128  436  390  262  120   87   54
table(colData(sce)[,batch])
## 
##  DMSO fresh MetOH 
## 11463  9004  2357
res_de <- readRDS(params$de)
abund <- readRDS(params$abund)
outputfile <- params$out_file

cols <-c(c(jcolors('pal6'),jcolors('pal8'))[c(1,8,14,5,2:4,6,7,9:13,15:20)],jcolors('pal4'))
names(cols) <- c()

Visualize data

How are sample, celltypes and batches distributed within normalized, but not batch corrected data?

feature_list <- c(batch, celltype, sample)
feature_list <- feature_list[which(!is.na(feature_list))]

lapply(feature_list, function(feature_name){
  visGroup(sce, feature_name, dim_red= "UMAP")
})
## [[1]]

## 
## [[2]]

## 
## [[3]]

Batch strength/size

To compare or describe the severity of a batch effect there are different meassures. In general they can either give an estimate of the relative strength compared to the signal of interest e.g. the celltype signal or an absolut estimate e.g. the number of batch affected genes.

Variance partitioning

How much of the variance within the datasets can we attributed to the batch effect and how much could be explained by the celltype? Which genes are mostly affected?

vp_vars <- c("vp_batch", "vp_celltype", "vp_residuals")
vp <- as_tibble(rowData(sce)[, vp_vars])  %>% dplyr::mutate(gene= rownames(sce)) %>% dplyr::arrange(-vp_batch)
vp_sub <- vp[1:3] %>% set_rownames(vp$gene)
## Warning: Setting row names on a tibble is deprecated.
#plot
plotPercentBars( vp_sub[1:10,] )

plotVarPart( vp_sub )

Variance and gene expression

Are general expression and batch effect related? Does the batch effect or the celltype effect preferable manifest within highly, medium or low expressed genes?

#define expression classes by mean expression quantiles 
th <- quantile(rowMeans(assays(sce)$logcounts), c(.33, .66))
high_th <- th[2]
mid_th <- th[1]

rowData(sce)$expr_class <- ifelse(rowMeans(assays(sce)$logcounts) > high_th, "high",
                                  ifelse(rowMeans(assays(sce)$logcounts) <= high_th &
                                           rowMeans(assays(sce)$logcounts) > mid_th, 
                                         "medium", "low"))
rowData(sce)$mean_expr <- rowMeans(assays(sce)$logcounts)

#plot 
plot_dev <- function(var, var_col){
  ggplot(as.data.frame(rowData(sce)), aes_string(x = "mean_expr", y = var, colour = var_col)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE)
}

#Ternary plots
# ggtern(data=as.data.frame(rowData(sce)),aes(vp_batch, vp_celltype, vp_residuals)) + 
#   stat_density_tern(aes(fill=..level.., alpha=..level..),geom='polygon') +
#   scale_fill_gradient2(high = "red") +
#   guides(color = "none", fill = "none", alpha = "none") +
#   geom_point(size= 0.1, alpha = 0.5)  + 
#   Llab("batch") +
#   Tlab("celltype") +
#   Rlab("other") +
#   theme_bw()
# 
# t1 <- ggtern(data=as.data.frame(rowData(sce)),aes(vp_batch, vp_celltype, vp_residuals)) + 
#   geom_point(size = 0.1) +
#   geom_density_tern() + 
#   Llab("batch") +
#   Tlab("celltype") +
#   Rlab("other") +
#   theme_bw()

## Summarize variance partitioning
# How many genes have a variance component affected by batch with > 1%
n_batch_gene <- vp_sub %>% dplyr::filter(vp_batch > 0.01) %>% nrow()/n_genes
n_batch_gene10 <- vp_sub %>% dplyr::filter(vp_batch > 0.1) %>% nrow()/n_genes
n_celltype_gene <- vp_sub %>% dplyr::filter(vp_celltype> 0.01) %>% nrow()/n_genes
n_rel <- n_batch_gene/n_celltype_gene

# Mean variance that is explained by the batch effect/celltype
m_batch <- mean(vp_sub$vp_batch, na.rm = TRUE)
m_celltype <- mean(vp_sub$vp_celltype, na.rm = TRUE)
m_rel <- m_batch/m_celltype

Scatterplot batch

plot_dev("vp_batch", "vp_batch")
## `geom_smooth()` using formula 'y ~ x'

Scatterplot celltype

plot_dev("vp_celltype", "vp_celltype")
## `geom_smooth()` using formula 'y ~ x'

Ternary plot all genes

#t1

Ternary plot gene expression classes

#t1 + facet_grid(~expr_class)

Cellspecific Mixing score

Overall

#visualize overall cms score
visHist(sce, n_col = 2, prefix = FALSE)

visMetric(sce, metric = "cms_smooth", dim_red = "UMAP")

visGroup(sce, celltype, dim_red = "UMAP")

#summarize
mean_cms <- mean(sce$cms)
n_cms_0.01 <- length(which(sce$cms < 0.01))
cluster_mean_cms <- as_tibble(colData(sce)) %>% group_by_at(celltype) %>% summarize(cms_mean = mean(cms))
var_cms <- var(cluster_mean_cms$cms_mean)

Celltypes cms smooth

#compare by celltypes
visCluster(sce, metric_var = "cms_smooth", cluster_var = celltype)
## Picking joint bandwidth of 0.0261

visCluster(sce, metric_var = "cms_smooth", cluster_var = celltype, violin = TRUE)

Celltypes histogram

#compare histogram by celltype

p <- ggplot(as.data.frame(colData(sce)), 
            aes_string(x = "cms", fill = celltype)) + 
                geom_histogram() + 
                facet_wrap(celltype, scales = "free_y", ncol = 3) +
                scale_fill_manual(values = cols) +
                theme_classic()

p + geom_vline(aes_string(xintercept = "cms_mean", 
                   colour = celltype), 
               cluster_mean_cms, linetype=2) + 
  scale_color_manual(values = cols) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Celltype specificity

Celltype abundance

meta_tib <- as_tibble(colData(sce)) %>% group_by_at(c(batch, celltype)) %>% summarize(n = n()) %>% dplyr::mutate(cell_freq = n / sum(n))

plot_abundance <- function(cluster_var, tib, x_var){
  meta_df <- as.data.frame(eval(tib))
  p <- ggplot(data=meta_df, aes_string(x=x_var, y="cell_freq", fill = cluster_var)) +
    geom_bar(stat="identity") + scale_fill_manual(values=cols, name = "celltype")
  p + coord_flip() + theme_minimal() 
}

plot_abundance(cluster_var = celltype, tib = meta_tib, x_var = batch)

#summarize diff abundance
mean_rel_abund_diff <- mean(unlist(abund))
min_rel_abund_diff <- min(unlist(abund))
max_rel_abund_diff <- max(unlist(abund))

Batch and celltype specific count distributions

Do the overall count distribution vary between batches? Are count distributions celltype depended

#batch level
bids <- levels(as.factor(colData(sce)[, batch]))
names(bids) <- bids
cids <- levels(as.factor(colData(sce)[, celltype]))
names(cids) <- cids

#mean gene expression by batch and cluster
mean_list <- lapply(bids, function(batch_var){
  mean_cluster <- lapply(cids, function(cluster_var){
    counts_sc <- as.matrix(logcounts(
      sce[, colData(sce)[, batch] %in% batch_var & 
            colData(sce)[, celltype] %in% cluster_var]))
  })
  mean_c <- mean_cluster %>% map(rowMeans) %>% bind_rows %>%
    dplyr::mutate(gene=rownames(sce)) %>% 
    gather(cluster, logcounts, cids)
})
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(cids)` instead of `cids` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
mean_expr <- mean_list %>% bind_rows(.id= "batch")

ggplot(mean_expr, aes(x=logcounts, colour=batch)) + geom_density(alpha=.3) +
  theme_classic() +
  facet_wrap( ~ cluster, ncol = 3) +
  scale_colour_manual(values = cols) +
  scale_x_continuous(limits =  c(0, 7))
## Warning: Removed 36 rows containing non-finite values (stat_density).

Batch to batch comparisons of expression distributions

Differentially expressed genes

Upset plot

## Upset plot\
cont <- param[["cont"]]
cs <- names(cont)
names(cs) <- cs

# Filter DEG by pvalue
FilterDEGs <- function (degDF = df, filter = c(FDR = 5)){
  rownames(degDF) <- degDF$gene
  #pval <- degDF[, grep("adj.P.Val$", colnames(degDF)), drop = FALSE]
  pval <- degDF[, grep("PValue$", colnames(degDF)), drop = FALSE]
  pf <- pval <= filter["FDR"]/100
  pf[is.na(pf)] <- FALSE
  DEGlistUPorDOWN <- sapply(colnames(pf), function(x) rownames(pf[pf[, x, drop = FALSE], , drop = FALSE]), simplify = FALSE)
}

result <- list()
m2 <- list()

for(jj in 1:length(cs)){
  result[[jj]] <- sapply(res_de[[1]][[names(cs)[jj]]], function(x) FilterDEGs(x))
  names(result[[jj]]) <- cids
  m2[[jj]] = make_comb_mat(result[[jj]], mode = "intersect")
}

names(result) <- names(cs)
names(m2) <- names(cs)

lapply(m2, function(x) UpSet(x))
## $`DMSO-fresh`

## 
## $`MetOH-Fresh`

## 
## $`DMSO-MetOH`

Logfold_change and GSEA

# DE genes (per cluster and mean)
res <- res_de[["table"]]
#n_de <- lapply(res, function(y) vapply(y, function(x) sum(x$adj.P.Val < 0.05), numeric(1)))
n_de <- lapply(res, function(y) vapply(y, function(x) sum(x$adj.PValue < 0.05), numeric(1)))
n_genes_lfc1 <- lapply(res, function(y) vapply(y, function(x) sum(abs(x$logFC) > 1), numeric(1)))
mean_n_genes_lfc1 <- mean(unlist(n_genes_lfc1))/n_genes

# plot DE for all comparison and check gene sets
#get geneset
gs <- read_delim(params$gs, delim = "\n", col_names = "cat")
## Parsed with column specification:
## cols(
##   cat = col_character()
## )
cats <- sapply(gs$cat, function(u) strsplit(u, "\t")[[1]][-2], 
               USE.NAMES = FALSE)
names(cats) <- sapply(cats, .subset, 1)
cats <- lapply(cats, function(u) u[-1])

plotDE <- function(cont_var){
  #res_s <- res[[cont_var]] %>% map(filter, adj.P.Val < .05) %>% map(filter, abs(logFC) > 1)
  res_s <- res[[cont_var]] %>% map(dplyr::filter, PValue < .05) %>% map(dplyr::filter, abs(logFC) > 1)
  
  #plot
  lapply(names(res[[cont_var]]), function(ct){
    ct_de <- res[[cont_var]][[ct]]
    ct_de$gene <- gsub('[A-z0-9]*\\.', '', ct_de$gene)
    res_s[[ct]]$gene <- gsub('[A-z0-9]*\\.', '', res_s[[ct]]$gene)
    
    #p <- ggplot(ct_de, aes(x = AveExpr, y = logFC, colour = abs(logFC) > 1, label = gene)) + 
    p <- ggplot(ct_de, aes(x = logCPM, y = logFC, colour = abs(logFC) > 2, label = gene)) + 
      geom_point(size = 2, alpha = .5) +
      geom_label_repel(data = res_s[[ct]]) +
      ggtitle(paste0(ct,": ", cont_var)) +
      theme_classic()
    print(p)
    
    cat("Cluster:", ct, "Contrast:", cont_var, 
        "Num genes:", nrow(ct_de), "Num DE:", nrow(res_s[[ct]]), "\n" )

    # run 'camera' for this comparison
    inds <- ids2indices(cats, ct_de$gene)
    #cm <- cameraPR(ct_de$t, inds) #LR
    cm <- cameraPR(ct_de$F, inds)
    print(cm %>% rownames_to_column("category") %>% 
      filter(FDR < .05 & NGenes >= 5) %>% head(8))
  })
}

if( length(names(res)) <= 3 ){
    pathways <- lapply(names(res), plotDE)
}

## Cluster: 0 Contrast: DMSO-fresh Num genes: 8331 Num DE: 6 
##                                category NGenes Direction       PValue
## 1 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 1.838334e-63
## 2       GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up 1.770386e-22
## 3                       GO_RRNA_BINDING     47        Up 2.825708e-16
## 4                       GO_MRNA_BINDING    106        Up 5.995141e-10
## 5                     GO_OXYGEN_BINDING      6        Up 4.316126e-05
##            FDR
## 1 1.623249e-60
## 2 7.816255e-20
## 3 8.317001e-14
## 4 1.323427e-07
## 5 7.622279e-03

## Cluster: 1 Contrast: DMSO-fresh Num genes: 8331 Num DE: 6 
##                                category NGenes Direction       PValue
## 1 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 6.962559e-53
## 2       GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up 8.577222e-19
## 3                       GO_RRNA_BINDING     47        Up 9.249362e-14
## 4                       GO_MRNA_BINDING    106        Up 1.200581e-07
##            FDR
## 1 6.147939e-50
## 2 3.786844e-16
## 3 2.722396e-11
## 4 2.650283e-05

## Cluster: 2 Contrast: DMSO-fresh Num genes: 8331 Num DE: 8 
##                                category NGenes Direction       PValue
## 1 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 2.936690e-37
## 2       GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up 3.597659e-14
## 3                       GO_RRNA_BINDING     47        Up 1.671247e-09
## 4                       GO_MRNA_BINDING    106        Up 4.343016e-06
##            FDR
## 1 2.593097e-34
## 2 1.588367e-11
## 3 4.919036e-07
## 4 9.587208e-04

## Cluster: 3 Contrast: DMSO-fresh Num genes: 8331 Num DE: 33 
##                                category NGenes Direction       PValue
## 1 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 8.874833e-13
## 2       GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up 1.047358e-05
## 3                GO_TROPOMYOSIN_BINDING      5        Up 1.204128e-05
## 4               GO_S100_PROTEIN_BINDING      9        Up 4.909950e-05
##            FDR
## 1 7.836478e-10
## 2 3.544150e-03
## 3 3.544150e-03
## 4 1.083871e-02

## Cluster: 4 Contrast: DMSO-fresh Num genes: 8331 Num DE: 5 
##                                category NGenes Direction       PValue
## 1 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 6.261523e-27
## 2       GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up 9.461716e-10
## 3                       GO_RRNA_BINDING     47        Up 1.715529e-06
## 4                     GO_OXYGEN_BINDING      6        Up 1.785109e-05
## 5                       GO_MRNA_BINDING    106        Up 4.337498e-05
##            FDR
## 1 5.528925e-24
## 2 4.177348e-07
## 3 5.049375e-04
## 4 3.940629e-03
## 5 7.660022e-03

## Cluster: 5 Contrast: DMSO-fresh Num genes: 8331 Num DE: 3 
##                                category NGenes Direction       PValue
## 1 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 1.360843e-32
## 2       GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up 3.849157e-12
## 3                       GO_MRNA_BINDING    106        Up 9.314689e-09
## 4                       GO_RRNA_BINDING     47        Up 1.255129e-05
##            FDR
## 1 1.201624e-29
## 2 1.699403e-09
## 3 2.741623e-06
## 4 2.770698e-03

## Cluster: 6 Contrast: DMSO-fresh Num genes: 8331 Num DE: 11 
##                                category NGenes Direction       PValue
## 1 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 1.492729e-17
## 2       GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up 2.993152e-07
##            FDR
## 1 1.318080e-14
## 2 1.321477e-04

## Cluster: 7 Contrast: DMSO-fresh Num genes: 8331 Num DE: 6 
##                                category NGenes Direction       PValue
## 1 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 2.336582e-30
## 2       GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up 5.517514e-11
## 3                       GO_RRNA_BINDING     47        Up 2.961837e-06
## 4                       GO_MRNA_BINDING    106        Up 5.691083e-05
##            FDR
## 1 2.063202e-27
## 2 2.435982e-08
## 3 8.717674e-04
## 4 1.256307e-02

## Cluster: 8 Contrast: DMSO-fresh Num genes: 8331 Num DE: 11 
##                                category NGenes Direction       PValue
## 1            GO_PEPTIDE_ANTIGEN_BINDING     18        Up 8.245751e-11
## 2 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 1.541058e-09
## 3     GO_MHC_CLASS_II_RECEPTOR_ACTIVITY      8        Up 2.265937e-07
## 4                    GO_ANTIGEN_BINDING     66        Up 1.709366e-05
## 5                 GO_CHEMOKINE_ACTIVITY     17        Up 3.208435e-05
## 6     GO_CCR_CHEMOKINE_RECEPTOR_BINDING     14        Up 6.619015e-05
## 7         GO_CHEMOKINE_RECEPTOR_BINDING     24        Up 2.055289e-04
##            FDR
## 1 7.280998e-08
## 2 6.803771e-07
## 3 6.669409e-05
## 4 3.773426e-03
## 5 5.666097e-03
## 6 9.740984e-03
## 7 2.592601e-02

## Cluster: 9 Contrast: DMSO-fresh Num genes: 8331 Num DE: 30 
##                                                    category NGenes Direction
## 1                                         GO_OXYGEN_BINDING      6        Up
## 2 GO_OXIDOREDUCTASE_ACTIVITY_ACTING_ON_PEROXIDE_AS_ACCEPTOR     23        Up
## 3                                   GO_TETRAPYRROLE_BINDING     35        Up
## 4                                       GO_IRON_ION_BINDING     45        Up
## 5                        GO_MAP_KINASE_PHOSPHATASE_ACTIVITY     11        Up
## 6                       GO_SOLUTE_PROTON_SYMPORTER_ACTIVITY      8        Up
## 7                       GO_SOLUTE_CATION_SYMPORTER_ACTIVITY     20        Up
##          PValue          FDR
## 1 2.308877e-102 2.038739e-99
## 2  6.165724e-21 2.722167e-18
## 3  3.186288e-13 9.378309e-11
## 4  2.509152e-09 5.538954e-07
## 5  2.726391e-07 4.814807e-05
## 6  4.675847e-06 6.881288e-04
## 7  4.186948e-04 4.621344e-02

## Cluster: 10 Contrast: DMSO-fresh Num genes: 8331 Num DE: 23 
##                                category NGenes Direction       PValue
## 1 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 1.085872e-11
## 2       GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up 4.043869e-05
##            FDR
## 1 9.588252e-09
## 2 1.785368e-02

## Cluster: 11 Contrast: DMSO-fresh Num genes: 8331 Num DE: 399 
##                                category NGenes Direction       PValue
## 1 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 3.146508e-13
## 2       GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up 4.042991e-05
## 3                       GO_RRNA_BINDING     47        Up 7.796956e-05
##            FDR
## 1 2.778367e-10
## 2 1.784980e-02
## 3 2.294904e-02

## Cluster: 12 Contrast: DMSO-fresh Num genes: 8331 Num DE: 126 
##                                category NGenes Direction       PValue
## 1 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 2.121785e-15
## 2       GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up 2.570222e-05
## 3                       GO_RRNA_BINDING     47        Up 6.452382e-05
##            FDR
## 1 1.873536e-12
## 2 1.134753e-02
## 3 1.899151e-02

## Cluster: 13 Contrast: DMSO-fresh Num genes: 8331 Num DE: 1021 
##                                                          category NGenes
## 1                 GO_OXIDOREDUCTASE_ACTIVITY_OXIDIZING_METAL_IONS      6
## 2                                   GO_CYTOKINE_RECEPTOR_ACTIVITY     36
## 3 GO_OXIDOREDUCTASE_ACTIVITY_ACTING_ON_NAD_P_H_OXYGEN_AS_ACCEPTOR      5
## 4                GO_PHOSPHATIDYLINOSITOL_3_4_BISPHOSPHATE_BINDING     11
## 5                                           GO_COMPLEMENT_BINDING     10
## 6                                    GO_PEPTIDE_RECEPTOR_ACTIVITY     21
## 7                                         GO_S100_PROTEIN_BINDING      9
## 8                                             GO_IRON_ION_BINDING     45
##   Direction       PValue          FDR
## 1        Up 2.298351e-08 1.199705e-05
## 2        Up 8.232470e-07 2.058644e-04
## 3        Up 9.325682e-07 2.058644e-04
## 4        Up 3.312243e-06 5.849422e-04
## 5        Up 1.018324e-05 1.498634e-03
## 6        Up 3.783337e-05 4.772409e-03
## 7        Up 6.009060e-05 6.632500e-03
## 8        Up 6.828096e-05 6.699121e-03

## Cluster: 0 Contrast: MetOH-Fresh Num genes: 8331 Num DE: 25 
## [1] category  NGenes    Direction PValue    FDR      
## <0 rows> (or 0-length row.names)

## Cluster: 1 Contrast: MetOH-Fresh Num genes: 8331 Num DE: 28 
## [1] category  NGenes    Direction PValue    FDR      
## <0 rows> (or 0-length row.names)

## Cluster: 2 Contrast: MetOH-Fresh Num genes: 8331 Num DE: 30 
##                                                    category NGenes Direction
## 1                                         GO_OXYGEN_BINDING      6        Up
## 2 GO_OXIDOREDUCTASE_ACTIVITY_ACTING_ON_PEROXIDE_AS_ACCEPTOR     23        Up
## 3                                   GO_TETRAPYRROLE_BINDING     35        Up
## 4                                       GO_IRON_ION_BINDING     45        Up
##         PValue          FDR
## 1 1.138330e-77 1.005145e-74
## 2 1.582975e-15 6.988836e-13
## 3 5.859002e-12 1.724500e-09
## 4 2.429449e-08 5.363008e-06

## Cluster: 3 Contrast: MetOH-Fresh Num genes: 8331 Num DE: 122 
##                           category NGenes Direction       PValue          FDR
## 1         GO_RAGE_RECEPTOR_BINDING      9        Up 6.659767e-12 5.880574e-09
## 2 GO_LONG_CHAIN_FATTY_ACID_BINDING      6        Up 4.008544e-05 1.769772e-02
## 3           GO_TROPOMYOSIN_BINDING      5        Up 1.548790e-04 4.558604e-02

## Cluster: 4 Contrast: MetOH-Fresh Num genes: 8331 Num DE: 19 
## [1] category  NGenes    Direction PValue    FDR      
## <0 rows> (or 0-length row.names)

## Cluster: 5 Contrast: MetOH-Fresh Num genes: 8331 Num DE: 22 
## [1] category  NGenes    Direction PValue    FDR      
## <0 rows> (or 0-length row.names)

## Cluster: 6 Contrast: MetOH-Fresh Num genes: 8331 Num DE: 43 
##                   category NGenes Direction       PValue        FDR
## 1 GO_RAGE_RECEPTOR_BINDING      9        Up 1.976667e-05 0.01745397

## Cluster: 7 Contrast: MetOH-Fresh Num genes: 8331 Num DE: 39 
##                                                    category NGenes Direction
## 1                                         GO_OXYGEN_BINDING      6        Up
## 2 GO_OXIDOREDUCTASE_ACTIVITY_ACTING_ON_PEROXIDE_AS_ACCEPTOR     23        Up
## 3                                   GO_TETRAPYRROLE_BINDING     35        Up
## 4                                       GO_IRON_ION_BINDING     45        Up
##         PValue          FDR
## 1 1.973522e-66 1.742620e-63
## 2 2.051177e-16 9.055944e-14
## 3 5.612312e-12 1.651891e-09
## 4 3.175164e-07 7.009175e-05

## Cluster: 8 Contrast: MetOH-Fresh Num genes: 8331 Num DE: 119 
##                                                    category NGenes Direction
## 1                                         GO_OXYGEN_BINDING      6        Up
## 2 GO_OXIDOREDUCTASE_ACTIVITY_ACTING_ON_PEROXIDE_AS_ACCEPTOR     23        Up
## 3                                   GO_TETRAPYRROLE_BINDING     35        Up
## 4                                       GO_IRON_ION_BINDING     45        Up
##         PValue          FDR
## 1 3.281320e-69 2.897406e-66
## 2 2.451515e-13 1.082344e-10
## 3 5.660405e-08 1.666046e-05
## 4 1.185453e-06 2.616887e-04

## Cluster: 9 Contrast: MetOH-Fresh Num genes: 8331 Num DE: 68 
## [1] category  NGenes    Direction PValue    FDR      
## <0 rows> (or 0-length row.names)

## Cluster: 10 Contrast: MetOH-Fresh Num genes: 8331 Num DE: 434 
##                                                    category NGenes Direction
## 1                                         GO_OXYGEN_BINDING      6        Up
## 2 GO_OXIDOREDUCTASE_ACTIVITY_ACTING_ON_PEROXIDE_AS_ACCEPTOR     23        Up
## 3                                   GO_TETRAPYRROLE_BINDING     35        Up
## 4                                       GO_IRON_ION_BINDING     45        Up
##          PValue           FDR
## 1 2.570759e-282 2.269980e-279
## 2  3.014323e-59  1.330824e-56
## 3  1.259411e-34  3.706867e-32
## 4  7.641491e-29  1.686859e-26

## Cluster: 11 Contrast: MetOH-Fresh Num genes: 8331 Num DE: 471 
##                                category NGenes Direction       PValue
## 1 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 5.930564e-24
## 2                       GO_RRNA_BINDING     47        Up 1.261220e-11
## 3       GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up 2.380511e-10
## 4   GO_FIBROBLAST_GROWTH_FACTOR_BINDING     11        Up 1.451591e-05
## 5    GO_TRNA_METHYLTRANSFERASE_ACTIVITY      7        Up 4.160682e-05
## 6            GO_PEPTIDE_ANTIGEN_BINDING     18        Up 8.955881e-05
##            FDR
## 1 5.236688e-21
## 2 5.568285e-09
## 3 7.006638e-08
## 4 3.204388e-03
## 5 7.347764e-03
## 6 1.318007e-02

## Cluster: 12 Contrast: MetOH-Fresh Num genes: 8331 Num DE: 221 
##                                                               category NGenes
## 1                                                    GO_OXYGEN_BINDING      6
## 2            GO_OXIDOREDUCTASE_ACTIVITY_ACTING_ON_PEROXIDE_AS_ACCEPTOR     23
## 3                                              GO_TETRAPYRROLE_BINDING     35
## 4                                                  GO_IRON_ION_BINDING     45
## 5 GO_NUCLEOBASE_CONTAINING_COMPOUND_TRANSMEMBRANE_TRANSPORTER_ACTIVITY     16
##   Direction        PValue           FDR
## 1        Up 7.675202e-175 6.777203e-172
## 2        Up  7.142274e-40  3.153314e-37
## 3        Up  4.089366e-21  1.203637e-18
## 4        Up  6.526673e-16  1.440763e-13
## 5        Up  2.289766e-04  4.043727e-02

## Cluster: 13 Contrast: MetOH-Fresh Num genes: 8331 Num DE: 426 
##                                     category NGenes Direction       PValue
## 1                GO_PHOSPHOLIPASE_C_ACTIVITY     10        Up 2.355024e-08
## 2                        GO_COLLAGEN_BINDING     16        Up 6.890368e-06
## 3                     GO_FIBRONECTIN_BINDING      8        Up 1.416769e-05
## 4   GO_PHOSPHORIC_DIESTER_HYDROLASE_ACTIVITY     30        Up 2.905674e-05
## 5                        GO_INTEGRIN_BINDING     40        Up 4.270163e-05
## 6 GO_TRANSFORMING_GROWTH_FACTOR_BETA_BINDING      8        Up 9.839959e-05
## 7        GO_FIBROBLAST_GROWTH_FACTOR_BINDING     11        Up 2.532854e-04
##            FDR
## 1 1.039743e-05
## 2 2.028065e-03
## 3 3.127517e-03
## 4 5.131420e-03
## 5 5.386505e-03
## 6 1.086085e-02
## 7 2.485011e-02

## Cluster: 0 Contrast: DMSO-MetOH Num genes: 8331 Num DE: 25 
##                                category NGenes Direction       PValue
## 1                     GO_OXYGEN_BINDING      6        Up 2.154031e-13
## 2 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 4.046929e-09
## 3           GO_PHOSPHOLIPASE_C_ACTIVITY     10        Up 3.052416e-08
##            FDR
## 1 1.902009e-10
## 2 1.786719e-06
## 3 8.984278e-06

## Cluster: 1 Contrast: DMSO-MetOH Num genes: 8331 Num DE: 32 
##                                category NGenes Direction       PValue
## 1 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 6.723456e-09
## 2           GO_PHOSPHOLIPASE_C_ACTIVITY     10        Up 7.099730e-07
## 3                     GO_OXYGEN_BINDING      6        Up 2.030438e-05
##            FDR
## 1 5.936812e-06
## 2 3.134531e-04
## 3 5.976256e-03

## Cluster: 2 Contrast: DMSO-MetOH Num genes: 8331 Num DE: 28 
##                                                    category NGenes Direction
## 1                                         GO_OXYGEN_BINDING      6        Up
## 2 GO_OXIDOREDUCTASE_ACTIVITY_ACTING_ON_PEROXIDE_AS_ACCEPTOR     23        Up
## 3                                   GO_TETRAPYRROLE_BINDING     35        Up
## 4                                       GO_IRON_ION_BINDING     45        Up
## 5                     GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up
## 6                               GO_PHOSPHOLIPASE_C_ACTIVITY     10        Up
## 7                           GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up
## 8                                 GO_PHOSPHOLIPASE_ACTIVITY     29        Up
##          PValue           FDR
## 1 2.336705e-196 2.063310e-193
## 2  3.925038e-37  1.732904e-34
## 3  4.162066e-28  1.225035e-25
## 4  1.358048e-18  2.997891e-16
## 5  3.184490e-13  5.623809e-11
## 6  1.437682e-08  2.115789e-06
## 7  3.298359e-05  3.640563e-03
## 8  2.213209e-04  2.171404e-02

## Cluster: 3 Contrast: DMSO-MetOH Num genes: 8331 Num DE: 419 
##                                category NGenes Direction       PValue
## 1 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 4.958307e-36
## 2              GO_RAGE_RECEPTOR_BINDING      9        Up 1.280768e-31
## 3                GO_TROPOMYOSIN_BINDING      5        Up 4.494973e-27
## 4               GO_S100_PROTEIN_BINDING      9        Up 2.673340e-23
## 5       GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up 4.234791e-15
## 6  GO_CALCIUM_DEPENDENT_PROTEIN_BINDING     29        Up 3.068241e-14
## 7      GO_LONG_CHAIN_FATTY_ACID_BINDING      6        Up 1.085093e-08
## 8                       GO_RRNA_BINDING     47        Up 9.970032e-08
##            FDR
## 1 4.378185e-33
## 2 5.654593e-29
## 3 1.323020e-24
## 4 5.901398e-21
## 5 7.478641e-13
## 6 4.515429e-12
## 7 1.197672e-06
## 8 9.781709e-06

## Cluster: 4 Contrast: DMSO-MetOH Num genes: 8331 Num DE: 20 
##                      category NGenes Direction       PValue          FDR
## 1 GO_PHOSPHOLIPASE_C_ACTIVITY     10        Up 1.612137e-08 1.423517e-05

## Cluster: 5 Contrast: DMSO-MetOH Num genes: 8331 Num DE: 27 
## [1] category  NGenes    Direction PValue    FDR      
## <0 rows> (or 0-length row.names)

## Cluster: 6 Contrast: DMSO-MetOH Num genes: 8331 Num DE: 37 
##                                                    category NGenes Direction
## 1                               GO_PHOSPHOLIPASE_C_ACTIVITY     10        Up
## 2                                  GO_RAGE_RECEPTOR_BINDING      9        Up
## 3                                  GO_PEPTIDOGLYCAN_BINDING      5        Up
## 4                                         GO_OXYGEN_BINDING      6        Up
## 5                                 GO_PHOSPHOLIPASE_ACTIVITY     29        Up
## 6 GO_MONOCARBOXYLIC_ACID_TRANSMEMBRANE_TRANSPORTER_ACTIVITY      8        Up
## 7                     GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up
## 8                                     GO_SYMPORTER_ACTIVITY     31        Up
##         PValue          FDR
## 1 5.013858e-12 4.427237e-09
## 2 8.585676e-07 3.790576e-04
## 3 2.799116e-05 8.238733e-03
## 4 7.663526e-05 1.552254e-02
## 5 8.789659e-05 1.552254e-02
## 6 1.059544e-04 1.559296e-02
## 7 1.331528e-04 1.679628e-02
## 8 2.823210e-04 3.116118e-02

## Cluster: 7 Contrast: DMSO-MetOH Num genes: 8331 Num DE: 40 
##                                                    category NGenes Direction
## 1                                         GO_OXYGEN_BINDING      6        Up
## 2 GO_OXIDOREDUCTASE_ACTIVITY_ACTING_ON_PEROXIDE_AS_ACCEPTOR     23        Up
## 3                                   GO_TETRAPYRROLE_BINDING     35        Up
## 4                                       GO_IRON_ION_BINDING     45        Up
## 5                     GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up
## 6                           GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up
##          PValue           FDR
## 1 3.770330e-190 3.329201e-187
## 2  9.601924e-42  4.239249e-39
## 3  2.880389e-26  8.477946e-24
## 4  1.342378e-18  2.963300e-16
## 5  1.865617e-14  3.294679e-12
## 6  1.820519e-05  2.679197e-03

## Cluster: 8 Contrast: DMSO-MetOH Num genes: 8331 Num DE: 79 
##                                                    category NGenes Direction
## 1                                         GO_OXYGEN_BINDING      6        Up
## 2 GO_OXIDOREDUCTASE_ACTIVITY_ACTING_ON_PEROXIDE_AS_ACCEPTOR     23        Up
## 3                                   GO_TETRAPYRROLE_BINDING     35        Up
## 4                                       GO_IRON_ION_BINDING     45        Up
## 5                                          GO_NADPH_BINDING      7        Up
##          PValue           FDR
## 1 3.710664e-162 3.276516e-159
## 2  1.542843e-34  6.811652e-32
## 3  4.036629e-19  1.188115e-16
## 4  5.001726e-15  1.104131e-12
## 5  1.111019e-04  1.962060e-02

## Cluster: 9 Contrast: DMSO-MetOH Num genes: 8331 Num DE: 28 
##                                                    category NGenes Direction
## 1                                         GO_OXYGEN_BINDING      6        Up
## 2 GO_OXIDOREDUCTASE_ACTIVITY_ACTING_ON_PEROXIDE_AS_ACCEPTOR     23        Up
## 3                                   GO_TETRAPYRROLE_BINDING     35        Up
## 4                                       GO_IRON_ION_BINDING     45        Up
## 5      GO_MONOSACCHARIDE_TRANSMEMBRANE_TRANSPORTER_ACTIVITY      6        Up
## 6                                GO_DEATH_RECEPTOR_ACTIVITY     14        Up
##          PValue           FDR
## 1 5.293795e-186 4.674421e-183
## 2  2.802013e-39  1.237089e-36
## 3  3.688863e-24  1.085755e-21
## 4  2.521594e-20  5.566419e-18
## 5  9.648539e-05  1.703932e-02
## 6  1.291875e-04  1.901209e-02

## Cluster: 10 Contrast: DMSO-MetOH Num genes: 8331 Num DE: 418 
##                                                    category NGenes Direction
## 1                                         GO_OXYGEN_BINDING      6        Up
## 2 GO_OXIDOREDUCTASE_ACTIVITY_ACTING_ON_PEROXIDE_AS_ACCEPTOR     23        Up
## 3                                   GO_TETRAPYRROLE_BINDING     35        Up
## 4                                       GO_IRON_ION_BINDING     45        Up
##         PValue          FDR
## 1 0.000000e+00 0.000000e+00
## 2 2.942990e-84 1.299330e-81
## 3 2.440511e-50 7.183237e-48
## 4 4.992880e-40 1.102178e-37

## Cluster: 11 Contrast: DMSO-MetOH Num genes: 8331 Num DE: 655 
##                                category NGenes Direction       PValue
## 1 GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME    169        Up 1.810356e-27
## 2       GO_STRUCTURAL_MOLECULE_ACTIVITY    314        Up 1.324484e-11
## 3                       GO_RRNA_BINDING     47        Up 8.771444e-09
## 4               GO_GLUCOSIDASE_ACTIVITY      5        Up 2.108795e-05
## 5   GO_SOLUTE_PROTON_SYMPORTER_ACTIVITY      8        Up 6.173424e-05
##            FDR
## 1 1.598544e-24
## 2 5.847596e-09
## 3 2.581728e-06
## 4 4.655166e-03
## 5 1.090227e-02

## Cluster: 12 Contrast: DMSO-MetOH Num genes: 8331 Num DE: 228 
##                                                    category NGenes Direction
## 1                                         GO_OXYGEN_BINDING      6        Up
## 2 GO_OXIDOREDUCTASE_ACTIVITY_ACTING_ON_PEROXIDE_AS_ACCEPTOR     23        Up
## 3                                   GO_TETRAPYRROLE_BINDING     35        Up
## 4                                       GO_IRON_ION_BINDING     45        Up
## 5                              GO_TBP_CLASS_PROTEIN_BINDING     18        Up
## 6                          GO_LONG_CHAIN_FATTY_ACID_BINDING      6        Up
##          PValue           FDR
## 1 5.842930e-260 5.159308e-257
## 2  1.132019e-55  4.997866e-53
## 3  2.850146e-30  8.388930e-28
## 4  2.481850e-23  5.478685e-21
## 5  5.632416e-05  9.946846e-03
## 6  1.453885e-04  2.139634e-02

## Cluster: 13 Contrast: DMSO-MetOH Num genes: 8331 Num DE: 412 
##                                           category NGenes Direction
## 1                           GO_FIBRONECTIN_BINDING      8        Up
## 2                      GO_PHOSPHOLIPASE_C_ACTIVITY     10        Up
## 3  GO_OXIDOREDUCTASE_ACTIVITY_OXIDIZING_METAL_IONS      6        Up
## 4    GO_PHOSPHOLIPID_TRANSLOCATING_ATPASE_ACTIVITY      6        Up
## 5                            GO_COMPLEMENT_BINDING     10        Up
## 6                              GO_COLLAGEN_BINDING     16        Up
## 7      GO_LOW_DENSITY_LIPOPROTEIN_PARTICLE_BINDING      7        Up
## 8 GO_PHOSPHATIDYLINOSITOL_3_4_BISPHOSPHATE_BINDING     11        Up
##         PValue          FDR
## 1 1.269442e-12 5.604587e-10
## 2 2.174672e-09 6.400786e-07
## 3 2.526340e-07 5.576895e-05
## 4 8.149231e-07 1.439154e-04
## 5 2.428015e-06 3.573229e-04
## 6 1.205113e-05 1.520164e-03
## 7 2.242339e-05 2.474982e-03
## 8 4.878419e-05 4.786271e-03

Summarize differential expression analysis

# DE genes (per cluster and mean)
#n_de <- lapply(res, function(y) vapply(y, function(x) sum(x$adj.P.Val < 0.05), numeric(1)))
n_de <- lapply(res, function(y) vapply(y, function(x) sum(x$PValue < 0.05), numeric(1)))
n_de_cl <- lapply(res, function(y) vapply(y, function(x) nrow(x), numeric(1)))
mean_n_de <- lapply(n_de, function(x) mean(x))
mean_mean_n_de <- mean(unlist(mean_n_de))/n_genes
min_mean_n_de <- min(unlist(mean_n_de))/n_genes
max_mean_n_de <- max(unlist(mean_n_de))/n_genes

# Genes with lfc > 1
n_genes_lfc1 <- lapply(res, function(y) vapply(y, function(x) sum(abs(x$logFC) > 1), numeric(1)))
mean_n_genes_lfc1 <- mean(unlist(n_genes_lfc1))/n_genes
min_n_genes_lfc1 <- min(unlist(n_genes_lfc1))/n_genes
max_n_genes_lfc1 <- max(unlist(n_genes_lfc1))/n_genes

# DE genes overlap between celltypes (celltype specific de genes)
# Genes are "overlapping" if they are present in all clusters with at least 10% of all cells
de_overlap <- lapply(result, function(x){
  result2 <- x[table(colData(sce)[, celltype]) > ncol(sce) * 0.1]
  de_overlap <- length(Reduce(intersect, result2))
  de_overlap
})

mean_de_overlap <- mean(unlist(de_overlap))/n_genes
min_de_overlap <- min(unlist(de_overlap))/n_genes
max_de_overlap <- max(unlist(de_overlap))/n_genes

#Genes unique to single celltypes
unique_genes_matrix <- NULL
unique_genes <- NULL
cb <- length(names(result[[1]]))
unique_genes <- lapply(result,function(x){
  for( i in 1:cb ){
    unique_genes[i] <-as.numeric(length(setdiff(unlist(x[i]),unlist(x[-i]))))
  }
  unique_genes_matrix <- cbind(unique_genes_matrix, unique_genes)
  unique_genes_matrix
})

unique_genes <- Reduce('cbind', unique_genes)
colnames(unique_genes) <- names(result)
rownames(unique_genes) <- names(result[[1]])

# Relative cluster specificity (unique/overlapping)
rel_spec1 <- NULL
for( i in 1:dim(unique_genes)[2] ){
  rel_spec <- unique_genes[,i]/de_overlap[[i]]
  rel_spec1 <- cbind(rel_spec1,rel_spec)
}

mean_rel_spec <- mean(rel_spec1)
min_rel_spec <- min(rel_spec1)
max_rel_spec <- max(rel_spec1)

Celltype specific DE distributions

How similar is the batch effect between celltypes. Do we have similar logFC distributions or different?

combine_folds <- function(cont_var){
  #extract the contrast of interest and change log2fold colums names to be unique
  B <- res[[cont_var]] 
  new_name <- function(p){
    colnames(B[[p]])[3] <- paste0("logFC_", p)
    return(B[[p]][,c(1,3)])
    }
  B_new_names <- lapply(names(B),new_name)
  names(B_new_names) <- names(B)
  #combine log2fold colums
  Folds <- Reduce(function(...){inner_join(..., by="gene")}, B_new_names)
}

all_folds <- lapply(cs, combine_folds)

#define pannels for pairs() function
panel.cor <- function(x, y, digits = 2, cex.cor){
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  txt <- format(c(r, 0.123456789), digits=digits)[1]
  test <- cor.test(x,y)
  Signif <- ifelse(round(test$p.value, 3) < 0.001,
                   "p<0.001",
                   paste("p=",round(test$p.value,3)))  
  text(0.5, 0.25, paste("r=",txt), cex = 3)
  text(.5, .75, Signif, cex = 3)
}

panel.smooth <- function (x, y, col = "blue", bg = NA, pch = 18, cex = 1.5, 
                          col.smooth = "red", span = 2/3, iter = 3, ...){
  points(x, y, pch = pch, col = col, bg = bg, cex = cex)
  ok <- is.finite(x) & is.finite(y)
  if( any(ok) ) 
    lines(stats::lowess(x[ok], y[ok], f = span, iter = iter), 
          col = col.smooth, ...)
}

panel.hist <- function(x, ...){
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(usr[1:2], 0, 1.5) )
  h <- hist(x, plot = FALSE)
  breaks <- h$breaks
  nB <- length(breaks)
  y <- h$counts
  y <- y/max(y)
  rect(breaks[-nB], 0, breaks[-1], y, col="cyan", ...)
}
#plot correlations
lapply(names(all_folds), function(x) pairs(all_folds[[x]][,-1],
                                           lower.panel = panel.smooth, 
                                           upper.panel = panel.cor, 
                                           diag.panel = panel.hist, main = x))

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
#extract correlation coefficients
# correlation coefficients from celltype specific gege logFC

lfc_cor_list <-lapply(names(all_folds), function(com){
  exclude <- which(table(colData(sce)[,celltype]) < 100)
  r <- cor(all_folds[[com]][, -c(1, (exclude + 1))])
  mean_r <- (sum(r) - ncol(r))/ (ncol(r)^2 - ncol(r))
})

mean_lfc_cor <- mean(unlist(lfc_cor_list))

Batch categorization

How does the batch effect manifest? Can we describe it by “simple” mean shifts of expression levels for some genes for all the cells in a given celltype and batch? Can we “remove” the batch effcet using a linear model with batch, batch and celltype or batch and celltype interacting?

#Visualize different models
vis_type <- function(dim_red){
  g <- visGroup(sce, batch, dim_red = dim_red) + 
    ggtitle("unadjusted")
  g1 <- visGroup(sce, batch, dim_red = paste0(dim_red, "_Xadj1")) + 
    ggtitle("constant batch effect")
  g2 <- visGroup(sce, batch, dim_red = paste0(dim_red, "_Xadj2")) + 
    ggtitle("constant batch effect, different ct composition")
  g3 <- visGroup(sce, batch, dim_red = paste0(dim_red, "_Xadj3")) + 
    ggtitle("celltype and batch effect interact")
  do.call("grid.arrange", c(list(g, g1, g2, g3), ncol = 2))
}

PCA

vis_type("PCA")

UMAP

vis_type("UMAP")

Cellspecific Mixing Score

# #Cellspecific Mixing score (Batch effect strength after "removal")
visHist(sce, metric = c("cms", "cms.Xadj1", "cms.Xadj2", "cms.Xadj3"), prefix = FALSE)

visIntegration(sce, metric = "cms", metric_name = "cms")
## Picking joint bandwidth of 0.0246

Simulation parameter

Extract parameter to use as input into simualation

#percentage of batch affected genes
cond <- gsub("-.*", "", names(n_de))
cond <- c(cond, unique(gsub(".*-", "", names(n_de))))
cond <- unique(cond)
de_be_tab <- n_de %>% bind_cols()
de_cl_tab <- n_de_cl %>% bind_cols()

de_be <- cond %>% map(function(x){
  de_tab <- de_be_tab[, grep(x, colnames(de_be_tab))]
  de_be <- rowMeans(de_tab)
}) %>% bind_cols() %>% set_colnames(cond)

n_cl <- cond %>% map(function(x){
  cl_tab <- de_cl_tab[, grep(x, colnames(de_cl_tab))]
  de_cl <- rowMeans(cl_tab)
}) %>% bind_cols() %>% set_colnames(cond)


p_be <- de_be/n_cl
mean_p_be <- mean(colMeans(p_be))
min_p_be <- min(colMins(as.matrix(p_be)))
max_p_be <- max(colMaxs(as.matrix(p_be)))
sd_p_be <- mean(colSds(as.matrix(p_be)))
if(is.na(sd_p_be)){ sd_p_be <- 0 }


#### Percentage of celltype specific genes "p_ct"
n_de_unique <- lapply(result,function(x){
  de_genes <- unlist(x) %>% unique() %>% length()
  de_genes <- de_genes/length(x)
}) %>% bind_cols()


rel_spec2 <- NULL
for(i in 1:length(de_overlap)){
  rel_spec <- de_overlap[[i]]/mean(n_de[[i]][table(colData(sce)[, celltype]) > dim(expr)[2] * 0.1])
  rel_spec2 <- cbind(rel_spec2, rel_spec)
}

mean_p_ct <- 1 - mean(rel_spec2)
max_p_ct <- 1 - min(rel_spec2)
min_p_ct <- 1 - max(rel_spec2)
sd_p_ct <- sd(rel_spec2)
if(is.na(sd_p_ct)){ sd_p_ct <- 0 }

# Logfold change
#logFoldchange batch effect distribution
mean_lfc_cl <- lapply(res, function(y) vapply(y, function(x){
  #de_genes <- which(x$adj.P.Val < 0.05)
  de_genes <- which(x$adj.PValue < 0.05)
  mean_de <- mean(abs(x[, "logFC"]))}
  , numeric(1))) %>% bind_cols()

mean_lfc_be <- mean(colMeans(mean_lfc_cl, na.rm = TRUE))
min_lfc_be <- min(colMins(as.matrix(mean_lfc_cl), na.rm = TRUE))
max_lfc_be <- max(colMaxs(as.matrix(mean_lfc_cl), na.rm = TRUE))

Summarize batch effect

  • Batch size
  • Celltype specificity
  • “Batch genes”
  • batch type
#Size? How much of the variance can be attributed to the batch effect?
size <- data.frame("batch_genes_1per" = n_batch_gene,       # 1.variance partition
                   "batch_genes_10per" = n_batch_gene10,
                   "celltype_gene_1per" = n_celltype_gene,
                   "relative_batch_celltype" = n_rel,
                   "mean_var_batch" = m_batch,
                   "mean_var_celltype" = m_celltype,
                   "rel_mean_ct_batch" = m_rel,
                   "mean_cms" = mean_cms,                    #2.cms
                   "n_cells_cms_0.01" = n_cms_0.01,
                   "mean_mean_n_de_genes" = mean_mean_n_de,  #3.de genes
                   "max_mean_n_de_genes" = max_mean_n_de,
                   "min_mean_n_de_genes" = min_mean_n_de,
                   "mean_n_genes_lfc1" = mean_n_genes_lfc1,
                   "min_n_genes_lfc1" = min_n_genes_lfc1,
                   "max_n_genes_lfc1" = max_n_genes_lfc1,
                   "n_cells_total" = ncol(sce),              #4.general
                   "n_genes_total" = nrow(sce))


#Celltype-specificity? How celltype/cluster specific are batch effects? 
# Differences in size, distribution or abundance? Do we find correlations between lfcs,
# overlap in de genes, pathways? Interaction between ct and be?
celltype <- data.frame('mean_rel_abund_diff' = mean_rel_abund_diff, #1.abundance
                       'min_rel_abund_diff' = min_rel_abund_diff,
                       'max_rel_abund_diff' = max_rel_abund_diff,
                       "celltype_var_cms" = var_cms,                 #2.size/strength
                       "mean_de_overlap" = mean_de_overlap,
                       "min_de_overlap" = min_de_overlap,
                       "max_de_overlap" = max_de_overlap,
                       "mean_rel_cluster_spec"= mean_rel_spec,
                       "min_rel_cluster_spec"= min_rel_spec,
                       "max_rel_cluster_spec"= max_rel_spec,
                       "mean_lfc_cor" = mean_lfc_cor )


sim <- data.frame("mean_p_be" = mean_p_be,
                  "max_p_be" = max_p_be,
                  "min_p_be" = min_p_be,
                  "sd_p_be" = sd_p_be,
                  "mean_lfc_be" = mean_lfc_be,
                  "min_lfc_be" = min_lfc_be,
                  "max_lfc_be" = max_lfc_be,
                  "mean_p_ct"= mean_p_ct,
                  "min_p_ct"= min_p_ct,
                  "max_p_ct"= max_p_ct,
                  "sd_p_ct" = sd_p_ct)


summary <- cbind(size, celltype, sim) %>% set_rownames(dataset_name)

### -------------- save summary object ----------------------###
saveRDS(summary, file = outputfile)
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.6 LTS
## 
## Matrix products: default
## BLAS:   /home/aluetg/R/lib/R/lib/libRblas.so
## LAPACK: /home/aluetg/R/lib/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] readr_1.3.1                 ggrepel_0.8.2              
##  [3] CAMERA_1.42.0               xcms_3.8.2                 
##  [5] MSnbase_2.12.0              ProtGenerics_1.18.0        
##  [7] mzR_2.20.0                  Rcpp_1.0.3                 
##  [9] cowplot_1.0.0               scran_1.14.6               
## [11] gridExtra_2.3               ComplexHeatmap_2.2.0       
## [13] stringr_1.4.0               dplyr_0.8.5                
## [15] tidyr_1.0.2                 here_0.1                   
## [17] jcolors_0.0.4               purrr_0.3.3                
## [19] variancePartition_1.16.1    scales_1.1.0               
## [21] foreach_1.4.8               limma_3.42.2               
## [23] CellMixS_1.2.4              kSamples_1.2-9             
## [25] SuppDists_1.1-9.5           scater_1.14.6              
## [27] ggplot2_3.3.0               CellBench_1.2.0            
## [29] tibble_2.1.3                magrittr_1.5               
## [31] SingleCellExperiment_1.8.0  SummarizedExperiment_1.16.1
## [33] DelayedArray_0.12.2         BiocParallel_1.20.1        
## [35] matrixStats_0.55.0          Biobase_2.46.0             
## [37] GenomicRanges_1.38.0        GenomeInfoDb_1.22.0        
## [39] IRanges_2.20.2              S4Vectors_0.24.3           
## [41] BiocGenerics_0.32.0        
## 
## loaded via a namespace (and not attached):
##   [1] tidyselect_1.0.0         lme4_1.1-21              RSQLite_2.2.0           
##   [4] htmlwidgets_1.5.1        munsell_0.5.0            codetools_0.2-16        
##   [7] preprocessCore_1.48.0    statmod_1.4.34           withr_2.1.2             
##  [10] colorspace_1.4-1         knitr_1.28               rstudioapi_0.11         
##  [13] robustbase_0.93-5        mzID_1.24.0              labeling_0.3            
##  [16] GenomeInfoDbData_1.2.2   farver_2.0.3             bit64_0.9-7             
##  [19] rprojroot_1.3-2          vctrs_0.2.4              xfun_0.12               
##  [22] BiocFileCache_1.10.2     R6_2.4.1                 doParallel_1.0.15       
##  [25] ggbeeswarm_0.6.0         clue_0.3-57              rsvd_1.0.3              
##  [28] locfit_1.5-9.1           bitops_1.0-6             assertthat_0.2.1        
##  [31] nnet_7.3-13              beeswarm_0.2.3           gtable_0.3.0            
##  [34] affy_1.64.0              rlang_0.4.5              GlobalOptions_0.1.1     
##  [37] splines_3.6.1            acepack_1.4.1            impute_1.60.0           
##  [40] checkmate_2.0.0          BiocManager_1.30.10      yaml_2.2.1              
##  [43] reshape2_1.4.3           backports_1.1.5          Hmisc_4.3-1             
##  [46] MassSpecWavelet_1.52.0   RBGL_1.62.1              tools_3.6.1             
##  [49] ellipsis_0.3.0           affyio_1.56.0            gplots_3.0.3            
##  [52] RColorBrewer_1.1-2       ggridges_0.5.2           plyr_1.8.6              
##  [55] base64enc_0.1-3          progress_1.2.2           zlibbioc_1.32.0         
##  [58] RCurl_1.98-1.1           prettyunits_1.1.1        rpart_4.1-15            
##  [61] GetoptLong_0.1.8         viridis_0.5.1            cluster_2.1.0           
##  [64] colorRamps_2.3           data.table_1.12.8        circlize_0.4.8          
##  [67] RANN_2.6.1               pcaMethods_1.78.0        packrat_0.5.0           
##  [70] hms_0.5.3                evaluate_0.14            pbkrtest_0.4-8.6        
##  [73] XML_3.99-0.3             jpeg_0.1-8.1             shape_1.4.4             
##  [76] compiler_3.6.1           KernSmooth_2.23-16       ncdf4_1.17              
##  [79] crayon_1.3.4             minqa_1.2.4              htmltools_0.4.0         
##  [82] mgcv_1.8-31              Formula_1.2-3            lubridate_1.7.4         
##  [85] DBI_1.1.0                dbplyr_1.4.2             MASS_7.3-51.5           
##  [88] rappdirs_0.3.1           boot_1.3-24              Matrix_1.2-18           
##  [91] cli_2.0.2                vsn_3.54.0               gdata_2.18.0            
##  [94] igraph_1.2.4.2           pkgconfig_2.0.3          foreign_0.8-76          
##  [97] MALDIquant_1.19.3        vipor_0.4.5              dqrng_0.2.1             
## [100] multtest_2.42.0          XVector_0.26.0           digest_0.6.25           
## [103] graph_1.64.0             rmarkdown_2.1            htmlTable_1.13.3        
## [106] edgeR_3.28.1             DelayedMatrixStats_1.8.0 listarrays_0.3.1        
## [109] curl_4.3                 gtools_3.8.1             rjson_0.2.20            
## [112] nloptr_1.2.2.1           lifecycle_0.2.0          nlme_3.1-145            
## [115] BiocNeighbors_1.4.2      fansi_0.4.1              viridisLite_0.3.0       
## [118] pillar_1.4.3             lattice_0.20-40          httr_1.4.1              
## [121] DEoptimR_1.0-8           survival_3.1-11          glue_1.3.1              
## [124] png_0.1-7                iterators_1.0.12         bit_1.1-15.2            
## [127] stringi_1.4.6            blob_1.2.1               BiocSingular_1.2.2      
## [130] latticeExtra_0.6-29      caTools_1.18.0           memoise_1.1.0           
## [133] irlba_2.3.3