COVID_combat_minor.Rmd

---
title: "COVID COMBAT study"
subtitle: '[COvid-19 Multi-omics Blood ATlas (COMBAT) Consortium, Cell, 2022](https://doi.org/10.1016/j.cell.2022.01.012)'
author: "Developed by [Gabriel Hoffman](http://gabrielhoffman.github.io/)"
date: "Run on `r Sys.time()`"
output: 
  html_document:
    toc: true
    smart: true
    self_contained: true
---

<!---


cd /sc/arion/projects/CommonMind/hoffman/dreamlet_analysis/COVID_combat/
R --vanilla


system("git pull"); rmarkdown::render("COVID_combat_minor.Rmd");


# https://hoffmg01.hpc.mssm.edu/dreamlet_analysis/COVID_combat

# resolve issue reading in 'raw'
conda activate /sc/arion/work/hoffmg01/R.cache/R/basilisk/1.14.3/zellkonverter/1.13.3/zellkonverterAnnDataEnv-0.10.6

adata = ad.read(h5ad("local.h5ad")
adata.X = adata.raw.X
adata.write("covid_counts.h5ad")



# cd 

--->


```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  warning=FALSE,
  message=FALSE,
  error = FALSE,
  tidy = FALSE,
  dev = c("png", "pdf"),
  package.startup.message = FALSE,
  cache = TRUE,
  cache.lazy = FALSE)
```


```{r load, cache=FALSE}
library(SingleCellExperiment)
library(zellkonverter)
library(dreamlet)
library(zenith)
library(GSEABase)
library(qvalue)
library(RColorBrewer) 
library(cowplot)
library(kableExtra)
library(tidyverse) 
```

```{r read_data}
# read single cell RNA-seq
file = "/sc/arion/projects/CommonMind/hoffman/scRNAseq_data/covid19_combat/covid_counts.h5ad"
sce = readH5AD(file, use_hdf5=TRUE)
counts(sce) = assay(sce, 'X')

# use HGNC names
rownames(sce) = rowData(sce)$feature_name

# remove doublets
sce = sce[,sce$GEX_region != "E: Doublets"]

# remove nan cells
sce = sce[,sce$minor_subset != "nan"]
```


# UMAP provided with data
```{r plots.UMAP, fig.widht=8, fig.height=6}
plotProjection(sce, "X_umap", "minor_subset")
```



```{r merge.data}
# read metadata
file = "/sc/arion/projects/CommonMind/hoffman/scRNAseq_data/covid19_combat/CBD-KEY-CLINVAR/COMBAT_CLINVAR_for_processed.txt"
df = read.table(file, header=TRUE)

# filter and merge metadata
df = df[df$scRNASeq_sample_ID %in% sce$scRNASeq_sample_ID,]
idx = match(sce$scRNASeq_sample_ID, df$scRNASeq_sample_ID)
colData(sce) = cbind(colData(sce), df[idx,])

# For each donor, select the most severe sample
res = colData(sce) %>%
	as_tibble %>%
	group_by(donor_id) %>%
	select(scRNASeq_sample_ID, Source) %>%
	distinct %>%
	summarize(scRNASeq_sample_ID, Source, i = which.max(Source),
		use_sample_id = scRNASeq_sample_ID[which.max(Source)])

# subset to one sample per donor
sceSub = sce[,sce$scRNASeq_sample_ID %in% droplevels(res$use_sample_id)]

# Keep only Oxford samples since Source and Institute are confounded:
# Flu was the only one collected at St_Georges
sceSub = sceSub[,sceSub$Institute == "Oxford"]

# create pseudobulk
pb <- aggregateToPseudoBulk(sceSub,
    assay = "counts",     
    cluster_id = "minor_subset",  
    sample_id = "scRNASeq_sample_ID",
    verbose = FALSE)
```


## Properties of dataset
```{r properties}
df = data.frame(table(pb$Source))
colnames(df) = c("Disease status", "Count")

df %>%
  kbl(row.names=FALSE) %>% 
  kable_classic(full_width = FALSE) 

df = list()
df$`# Samples` = ncol(pb)
df$`# Cells` = ncol(sce)
df = unlist(df)
df = data.frame(Property=names(df), count = df)

df %>%  
  kbl(row.names=FALSE) %>% 
  kable_classic(full_width = FALSE) 
```


```{r proc}
form = ~ Age + (1|sex) + (1|Source) 
res.proc = processAssays(pb, form,  
                          min.count = 3,
                          min.cells = 3,
                          min.samples = 8,
                          BPPARAM = SnowParam(6))

form = ~ Age + (1|sex) + (1|Source)
res.vp = fitVarPart(res.proc, form, BPPARAM=SnowParam(6) )

form = ~ Source + Age + sex 
fit = dreamlet(res.proc, form, BPPARAM=SnowParam(6))
```


### Show details of processing
```{r details1}
details(res.proc) 
```


### Show voom-style mean-variance trends
```{r voom.plot, fig.height=14, fig.width=9}
# show voom plot for each cell clusters  
plotVoom( res.proc, ncol=4) 
```

```{r details2}
details(fit) 
```


# Variance partitioning results
```{r vp.plot, fig.height=20, fig.width=12}
plotVarPart(sortCols(res.vp), label.angle=60, ncol=4)
```




```{r volcano.plot, fig.height=28, fig.width=12}
plotVolcano(fit, coef="SourceCOVID_MILD", ncol=4) 
```

```{r examples, eval=FALSE}
coef_array = c("SourceCOVID_HCW_MILD", "SourceCOVID_MILD", "SourceCOVID_SEV", "SourceCOVID_CRIT", "SourceSepsis")

df = coef_array %>%
  map_df(function(x) 
    topTable(fit, coef = x, number=Inf) %>%
    as_tibble %>%
    mutate(coef = factor(x, coef_array))) %>%
    mutate(se = logFC / t) %>%
    group_by(coef,assay)

genes = c("IFI27", "PPARG", "ZNF217", "PIM3")

df %>% 
  filter(ID %in% genes) %>%
  filter(assay %in% c("cMono", 'ncMono')) %>%
  ggplot(aes(logFC, coef, color=coef)) +
    geom_point() +
    geom_errorbar(aes(xmin = logFC - 1.96*se, xmax = logFC + 1.96*se), width=0) +
    theme_classic() +
    theme(aspect.ratio=1) +
    facet_grid(ID ~assay)

genes = df %>%
  filter(assay %in% c("cMono", 'ncMono')) %>%
  arrange(-abs(logFC)) %>%
  ungroup %>%
  select(ID) %>%
  distinct %>%
  head(100) %>%
  c

tab = res.gsa %>%
  as_tibble %>%
  mutate(coef = factor(coef, coef_array)) %>%
  filter(assay %in% c("cMono", 'ncMono')) %>%
  arrange(PValue) 



gs = "interferon|interlukin|activiation|inflamm|extracell|chemokin"

tab2 = tab %>%
  filter(grepl(gs, Geneset)) %>%
  filter(assay == "cMono") %>%
  mutate(tstat = delta/ se)


nrow = length(unique(tab2$coef))
ncol = length(unique(tab2$Geneset))
zmax = max(abs(tab2$tstat), na.rm=TRUE)

fig = tab2 %>%
  ggplot(aes(coef, Geneset, fill=tstat, label=ifelse(FDR < 0.05, '*', ''))) +
    geom_tile() +    
    theme_classic() + 
    scale_fill_gradient2("t-statistic", low="blue", mid="white", high="red", limits=c(-zmax, zmax)) + 
    theme(aspect.ratio=ncol/nrow, axis.text.x = element_text(angle = 45, vjust = 1, hjust=1), plot.title = element_text(hjust = 0.5)) + 
    geom_text(vjust=1, hjust=0.5) 
fig

# this showed enrichment for type I and II interferon pathways in the less severe hospitalized COVID-19 patients across cell types.


gs = c("GO0002544: chronic inflammatory response", 
"GO0031012: extracellular matrix",
"GO0035456: response to interferon−beta",
"GO0071346: cellular response to type II interferon")

genes = geneIds(go.gs[[gs[4]]])

tab3 = df %>%
  filter(assay == c( "ncMono"), ID %in% genes) 

nrow = length(unique(tab3$coef))
ncol = length(unique(tab3$ID))
zmax = max(abs(tab3$logFC), na.rm=TRUE)

fig = tab3  %>%
  ggplot(aes(coef, ID, fill=logFC, label=ifelse(adj.P.Val < 0.05, '*', ''))) +
    geom_tile() +    
    theme_classic() + 
    scale_fill_gradient2("logFC", low="blue", mid="white", high="red", limits=c(-zmax, zmax)) + 
    theme(aspect.ratio=ncol/nrow, axis.text.x = element_text(angle = 45, vjust = 1, hjust=1), plot.title = element_text(hjust = 0.5)) + 
    geom_text(vjust=1, hjust=0.5) +
    facet_wrap(~ assay)
    
ggsave(fig, file="~/www/test.pdf", width=15, height=18)



```


### Summarize differential expression 
```{r plot.pi, fig.height=16, fig.width=15}
file = "./topTable_COVID_2022.tsv"

coef_array = c("SourceCOVID_HCW_MILD", "SourceCOVID_MILD", "SourceCOVID_SEV", "SourceCOVID_CRIT", "SourceSepsis")

coef_array %>%
  map_df(function(x) 
    topTable(fit, coef = x, number=Inf) %>%
    as_tibble %>%
    mutate(coef = x)) %>%
    write.table(file, quote=FALSE, sep='\t', row.names=FALSE)

# Summarize differential expression for each coef and assay
df = coef_array %>%
  map_df(function(x) 
    topTable(fit, coef = x, number=Inf) %>%
    as_tibble %>%
    mutate(coef = factor(x, coef_array))) %>%
  group_by(coef,assay) %>%
  summarize( nDE = sum(adj.P.Val < 0.05), 
            pi1 = 1 - pi0est(P.Value)$pi0,
            nGenes = length(adj.P.Val)) %>%
  mutate(assay = factor(assay, assayNames(pb))) 


ymax = 1.05*max(df$nGenes)
fig1 = ggplot(df, aes(nGenes, assay, fill=assay)) + 
    geom_bar(stat="identity") + 
    theme_classic() +
    theme(aspect.ratio=3, legend.position="none") +
    scale_x_continuous(limits=c(0,ymax), expand=c(0,0)) +
    xlab("# genes expressed") +
    ylab("Cell type") +
    facet_wrap(~ coef, nrow=1)

ymax = 1.05*max(df$nDE)
fig2 = ggplot(df, aes(nDE, assay, fill=assay)) + 
    geom_bar(stat="identity") + 
    theme_classic() +
    theme(aspect.ratio=3, legend.position="none") +
    scale_x_continuous(limits=c(0,ymax), expand=c(0,0)) +
    xlab("# genes with FDR < 5%") +
    ylab("Cell type") +
    facet_wrap(~ coef, nrow=1)

fig3 = ggplot(df, aes(pi1, assay, fill=assay)) + 
    geom_bar(stat="identity") + 
    theme_classic() +
    theme(aspect.ratio=3, legend.position="none") +
    scale_x_continuous(limits=c(0,1), expand=c(0,0)) +
    xlab(bquote(pi[1])) +
    ylab("Cell type") +
    facet_wrap(~ coef, nrow=1)

plot_grid(fig1, fig2, fig3, labels=LETTERS[1:3], nrow=3, axis="tblr", align="hv")
```    


## Gene set analysis
```{r zenith}
# Load Gene Ontology database 
go.gs = get_GeneOntology(to="SYMBOL")

# remove set with > 1000 genes
go.gs = go.gs[sapply(go.gs, function(x) length(geneIds(x))) < 1000]

# zenith analysis
res.gsa = zenith_gsa(fit, go.gs, coefs=coef_array, progressbar=FALSE)
```

```{r zenith.plot, fig.height=80, fig.width=20}
plotZenithResults(res.gsa, 5, 3)
```


```{r zenith.plot2, fig.height=40, fig.width=20}
i = grep("Mono|cDC", res.gsa$assay)
plotZenithResults(res.gsa[i,], 5, 3) + facet_wrap(~coef, nrow=1)
```



### MSigDB
```{r zenith.MSigDB}
# Load Gene Ontology database 
msigdb.gs = get_MSigDB("H", to="SYMBOL")

# zenith analysis
res.msigdb = zenith_gsa(fit, msigdb.gs, coefs=coef_array, progressbar=FALSE)
```

```{r zenith.plot.MSigDB, fig.height=40, fig.width=20}
plotZenithResults(res.msigdb, 5, 3)
```


```{r zenith.plot2.MSigDB, fig.height=20, fig.width=20}
i = grep("Mono|cDC", res.msigdb$assay)
fig = plotZenithResults(res.msigdb[i,], 5, 3) + facet_wrap(~coef, nrow=1)
fig

gs.array = c("M5936_HALLMARK_OXIDATIVE_PHOSPHORYLATION",
'M5892_HALLMARK_CHOLESTEROL_HOMEOSTASIS',
'M5936_HALLMARK_OXIDATIVE_PHOSPHORYLATION',
'M5890_HALLMARK_TNFA_SIGNALING_VIA_NFKB',
'M5921_HALLMARK_COMPLEMENT',
'M5913_HALLMARK_INTERFERON_GAMMA_RESPONSE')

gs1 = 'M5890_HALLMARK_TNFA_SIGNALING_VIA_NFKB'
genes = sort(geneIds(msigdb.gs[[gs1]]))
# Summarize differential expression for each coef and assay
df = coef_array %>%
  map_df(function(x) 
    topTable(fit, coef = x, number=Inf) %>%
    as_tibble %>%
    mutate(coef = factor(x, coef_array))) %>%
  mutate(assay = factor(assay, assayNames(pb))) %>%
  filter(grepl("Mono", assay), ID %in% genes) %>%
  droplevels 

# heatmap of significant genes
genes = df %>% 
  filter(assay == 'cMono') %>%
  group_by(ID) %>%
  summarize(FDR = min(adj.P.Val)) %>%
  filter(FDR < 0.001) %>%
  pull(ID) %>%
  sort

df2 = df %>%
    filter(assay == 'cMono', ID %in% genes) 

hcl = df2 %>%
  select(ID, t, coef) %>%
  pivot_wider(names_from = ID, values_from=t) %>%
  column_to_rownames(var='coef') %>%
  t %>%
  dist(.) %>%
  hclust

df2 = df2 %>%
      mutate(ID = factor(ID, hcl$labels[hcl$order]))

df2 %>%
    ggplot(aes(ID, coef, fill=t)) +
    geom_tile() +
    scale_fill_gradient2(low="blue", mid="white", high="red") +
    theme_classic() +
    coord_equal() +
    facet_wrap(~assay, ncol=4) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
    ggtitle(gs1) +
    geom_text(data = df2 %>% filter(adj.P.Val < 0.05), label='*', hjust=0.5, vjust=0.8) 

# Heatmap of Gene set results
df_gsa = res.msigdb %>%
          as_tibble %>%
          filter(Geneset %in% gs.array) %>%
          filter(grepl("Mono|cDC", assay)) %>%
          mutate(coef = factor(coef, coef_array)) %>%
          mutate(tstat = delta/se)         

zmax = max(abs(df_gsa$tstat))             

fig = df_gsa %>%
      ggplot(aes(assay , coef, fill=tstat)) +
        geom_tile() +
        scale_fill_gradient2(low="blue", mid="white", high="red", limits=c(-zmax, zmax)) +
        theme_classic() +
        coord_equal()  +
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
        geom_text(data = df_gsa %>% filter(FDR < 0.05), label='*', hjust=0.5, vjust=0.8) +
        facet_wrap(~Geneset, nrow=1)
fig
```


# Figures
```{r figures, fig.height=20, fig.width=12}
gs.array = c("GO0006119: oxidative phosphorylation",
"GO0070469: respirasome",
"GO0031012: extracellular matrix",
"GO0000502: proteasome complex",
"GO0038024: cargo receptor activity",
"GO0035456: response to interferon-beta",
"GO0034340: response to type I interferon")

gs1 = "GO0034340: response to type I interferon"

genes = sort(geneIds(go.gs[[gs1]]))
# Summarize differential expression for each coef and assay
df = coef_array %>%
  map_df(function(x) 
    topTable(fit, coef = x, number=Inf) %>%
    as_tibble %>%
    mutate(coef = factor(x, coef_array))) %>%
  mutate(assay = factor(assay, assayNames(pb))) %>%
  filter(grepl("Mono", assay), ID %in% genes) %>%
  droplevels 

# heatmap of significant genes
genes = df %>% 
  filter(assay == 'cMono') %>%
  group_by(ID) %>%
  summarize(FDR = min(adj.P.Val)) %>%
  filter(FDR < 0.01) %>%
  pull(ID) %>%
  sort

df2 = df %>%
    filter(assay == 'cMono', ID %in% genes) 

hcl = df2 %>%
  select(ID, t, coef) %>%
  pivot_wider(names_from = ID, values_from=t) %>%
  column_to_rownames(var='coef') %>%
  t %>%
  dist(.) %>%
  hclust

df2 = df2 %>%
      mutate(ID = factor(ID, hcl$labels[hcl$order]))

df2 %>%
    ggplot(aes(ID, coef, fill=t)) +
    geom_tile() +
    scale_fill_gradient2(low="blue", mid="white", high="red") +
    theme_classic() +
    coord_equal() +
    facet_wrap(~assay, ncol=4) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
    ggtitle(gs1) +
    geom_text(data = df2 %>% filter(adj.P.Val < 0.05), label='*', hjust=0.5, vjust=0.8) 

# Heatmap of Gene set results
df_gsa = res.gsa %>%
          as_tibble %>%
          filter(Geneset %in% gs.array) %>%
          filter(grepl("Mono|cDC", assay)) %>%
          mutate(coef = factor(coef, coef_array)) %>%
          mutate(tstat = delta/se)         

zmax = max(abs(df_gsa$tstat))          

fig = df_gsa %>%
      ggplot(aes(coef, Geneset, fill=tstat)) +
        geom_tile() +
        scale_fill_gradient2(low="blue", mid="white", high="red", limits=c(-zmax, zmax)) +
        theme_classic() +
        coord_equal()  +
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
        geom_text(data = df_gsa %>% filter(FDR < 0.05), label='*', hjust=0.5, vjust=0.8) +
        facet_wrap(~assay)
fig
```



```{r figures2, fig.height=8, fig.width=12}
library(paletteer)

# cluster, GEX_region, minor_subset, major_subset
incl = grep("cDC|Mono", sce$GEX_region)
figA = plotProjection(sce[,incl], "X_umap", "minor_subset", pointsize=2) + scale_color_manual(values = adjustcolor(paletteer_d("ggsci::default_nejm"), offset=c(.3, .3, .3, 0)))

figB = plotVoom(res.proc[["cMono"]])

incl = res.vp$assay == "cMono"
col = c('#3A488AFF', '#DABD61FF', '#BE3428FF', "grey85")
figC = plotVarPart( sortCols(res.vp[incl,]), col=col ) 

# statistics
sortCols(res.vp[incl,]) %>%
  as_tibble %>%
  summarize(median = median(Source), 
    N10 = sum(Source> .25), 
    length = length(Source))


figD = plotVolcano(fit[["cMono"]], coef="SourceCOVID_MILD")

df.vp = res.vp %>%
  as_tibble %>%
  filter(assay == "cMono") %>%
  arrange(Source, Age, sex, Residuals) %>%
  column_to_rownames(var='gene') %>%
  select(-assay)

genes = c("NFKBIA", "PIM3", "CLU", "PRDX2", "XIST", "CETN2")
genes = c(genes, 'NLRC5', 'GIGYF2', 'FADD', 'PTPN2', 'IFI27')
genes = genes[rev(order(match(genes, rownames(df.vp))))]
figE = plotPercentBars(sortCols(df.vp), genes = genes, assays="cMono", col=col) + theme(aspect.ratio=1, legend.position="none")

plot_grid(figA, figB, figC, figD, figE, ncol=4, labels=LETTERS)
```

```{r heatmap}

gs.array = c("GO0006119: oxidative phosphorylation",
"GO0070469: respirasome",
"GO0031012: extracellular matrix",
"GO0000502: proteasome complex",
"GO0038024: cargo receptor activity",
"GO0035456: response to interferon-beta",
"GO0034340: response to type I interferon")

gs1 = gs.array[6]
genes = sort(geneIds(go.gs[[gs1]]))

# Summarize differential expression for each coef and assay
df = coef_array %>%
  map_df(function(x) 
    topTable(fit, coef = x, number=Inf) %>%
    as_tibble %>%
    mutate(coef = factor(x, coef_array))) %>%
  mutate(assay = factor(assay, assayNames(pb))) %>%
  filter(grepl("Mono", assay), ID %in% genes) %>%
  droplevels 
  
fig = df %>% 
  ggplot(aes(coef, ID, fill=t)) +
    geom_tile() +
    scale_fill_gradient2(low="blue", mid="white", high="red") +
    theme_classic() +
    coord_equal() +
    facet_wrap(~assay, ncol=4) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
    ggtitle(gs1) +
    geom_text(data = df %>% filter(adj.P.Val < 0.05), label='*', hjust=0.5, vjust=0.8)



# heatmap of significant genes
genes = df %>% 
  filter(assay == 'cMono') %>%
  group_by(ID) %>%
  summarize(FDR = min(adj.P.Val)) %>%
  filter(FDR < 0.05) %>%
  pull(ID) %>%
  sort

df2 = df %>%
    filter(assay == 'cMono', ID %in% genes) 

hcl = df2 %>%
  select(ID, t, coef) %>%
  pivot_wider(names_from = ID, values_from=t) %>%
  column_to_rownames(var='coef') %>%
  t %>%
  dist(.) %>%
  hclust

df2 = df2 %>%
      mutate(ID = factor(ID, hcl$labels[hcl$order]))

fig = df2 %>%
    ggplot(aes(ID, coef, fill=t)) +
    geom_tile() +
    scale_fill_gradient2(low="blue", mid="white", high="red") +
    theme_classic() +
    coord_equal() +
    facet_wrap(~assay, ncol=4) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
    ggtitle(gs1) +
    geom_text(data = df2 %>% filter(adj.P.Val < 0.05), label='*', hjust=0.5, vjust=0.8) 


df_gsa = res.gsa %>%
          as_tibble %>%
          filter(grepl("nterferon", Geneset)) %>%
          filter(grepl("Mono|cDC", assay))

fig = df_gsa %>%
      ggplot(aes(coef, Geneset, fill=delta/se)) +
        geom_tile() +
        scale_fill_gradient2(low="blue", mid="white", high="red") +
        theme_classic() +
        coord_equal()  +
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
        geom_text(data = df_gsa %>% filter(FDR < 0.05), label='*', hjust=0.5, vjust=0.8) +
        facet_wrap(~assay, nrow=1)
fig

# Summarize differential expression for each coef and assay
df = coef_array %>%
  map_df(function(x) 
    topTable(fit, coef = x, number=Inf) %>%
    as_tibble %>%
    mutate(coef = factor(x, coef_array))) %>%
  mutate(assay = factor(assay, assayNames(pb))) %>%
  filter(grepl("Mono|cDC", assay), ID %in% genes) %>%
  droplevels 

figList = lapply(genes, function(gene){
  df %>%
  filter(ID == gene) %>%
  mutate(se = logFC / t) %>%
  ggplot(aes(coef, logFC,  color = -log10(pmax(1e-4, adj.P.Val)))) +
    geom_point() +
    theme_classic() +
    geom_errorbar(aes(ymin = logFC - 1.96 * se, ymax = logFC + 1.96 * se), width = 0) +
    facet_wrap(~assay) +
    theme(aspect.ratio=1) +
    coord_flip() +
    scale_color_gradient(name = bquote(-log[10] ~ FDR), low = "grey", high = "red", limits = c(0, 4)) +
    ggtitle(gene) +
    ylab(bquote(log[2] ~ fold ~ change)) +
    geom_hline(yintercept = 0, linetype = "dashed") 
})

figList


# heatmap
###########
genes = c("NLRC5","PTPN2", "FADD")

# Summarize differential expression for each coef and assay
df = coef_array %>%
  map_df(function(x) 
    topTable(fit, coef = x, number=Inf) %>%
    as_tibble %>%
    mutate(coef = factor(x, coef_array))) %>%
  mutate(assay = factor(assay, assayNames(pb))) %>%
  filter(grepl("Mono|cDC", assay)) %>%
  droplevels 


df2 = df %>%
  filter(ID %in% genes, assay == "cMono") %>%
  mutate(se = logFC / t)

zmax = 1.05*max(abs(df2$logFC))
fig1 = df2 %>%  
  ggplot(aes(coef, ID, fill=logFC)) +
    geom_tile() +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))  +
    scale_fill_gradient2("logFC", low = "blue", mid = "white", high = "red", limits = c(-zmax, zmax), na.value = "grey70") +
    geom_text(data = df2 %>% filter(adj.P.Val < 0.05), label='*', hjust=0.5, vjust=0.8) +
    coord_equal()

genes = c("NFKBIA", "CLU", "IFI27", "PIM3")

df2 = df %>%
  filter(ID %in% genes, assay == "cMono") %>%
  mutate(se = logFC / t)

zmax = 1.05*max(abs(df2$logFC))
fig2 = df2 %>%  
  ggplot(aes(coef, ID, fill=logFC)) +
    geom_tile() +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
    coord_equal() +
    scale_fill_gradient2("logFC", low = "blue", mid = "white", high = "red", limits = c(-zmax, zmax), na.value = "grey70") +
    geom_text(data = df2 %>% filter(adj.P.Val < 0.05), label='*', hjust=0.5, vjust=0.8) +
    coord_equal()

plot_grid(fig1 + theme(axis.text.x =  element_blank()), fig2, ncol=1, align="h", axis="lr")
```


```{r figures3, fig.height=8, fig.width=12, eval=FALSE, echo=FALSE}
library(paletteer)

# cluster, GEX_region, minor_subset, major_subset
incl = grep("cDC|Mono", sce$GEX_region)
figA = plotProjection(sce[,incl], "X_umap", "minor_subset", pointsize=2) + scale_color_manual(values = paletteer_d("ggsci::default_nejm"))

figB = plotVoom(res.proc[["cMono"]])

incl = res.vp$assay == "cMono"
col = c('#3A488AFF', '#DABD61FF', '#BE3428FF', "grey85")
figC = plotVarPart( sortCols(res.vp[incl,]), col=col ) 

figD = plotVolcano(fit[["cMono"]], coef="SourceCOVID_MILD")

df.vp = res.vp %>%
  as_tibble %>%
  filter(assay == "cMono") %>%
  column_to_rownames(var='gene') %>%
  select(-assay)

# df.vp %>%
#   arrange(-Age) %>%
#   head

# df.vp %>%
#   arrange(-Residuals) %>%
#   head

# df.vp %>%
#   arrange(-Source) %>%
#   head(10)


genes = c("NFKBIA", "PIM3", "CLU", "PRDX2", "XIST", "CETN2")
genes = c(genes, 'NLRC5', 'GIGYF2', 'AZI2', 'FADD', 'PTPN2', 'IFI27')
figE = plotPercentBars(sortCols(df.vp), genes = genes, assays="cMono", col=col) + theme(aspect.ratio=1, legend.position="none")


fig = plot_grid(figA, figB, figC, figD, figE, ncol=4, labels=LETTERS)
ggsave(file="~/www/test.pdf", fig, height=8, width=12)

# plotGeneHeatmap
CT = c('ncMono', 'cMono', 'cMono.cyc', 'cMono.PLT', 'cDC' )
fig = plotGeneHeatmap(fit, coef = "SourceCOVID_MILD", genes = genes, assays=CT) + ggtitle("SourceCOVID_MILD")
ggsave(file="~/www/test.pdf", fig, height=8, width=12)


# plotForest
fig = plotForest(fit, coef="SourceCOVID_MILD", gene="NFKBIA", assays=CT) +
  theme(aspect.ratio=1)

# Summarize differential expression for each coef and assay
df = coef_array %>%
  map_df(function(x) 
    topTable(fit, coef = x, number=Inf) %>%
    as_tibble %>%
    mutate(coef = factor(x, coef_array))) %>%
  mutate(assay = factor(assay, assayNames(pb))) %>%
  filter(grepl("Mono|cDC", assay), ID %in% genes) %>%
  droplevels 
  
zmax = 2.5
fig = df %>%
  filter(ID == "NFKBIA") %>%
  ggplot(aes(coef, assay, fill=logFC)) +
    geom_tile() +
    theme_classic() +
    coord_equal() +
    scale_fill_gradient2("z-statistic", low = "blue", mid = "white", high = "red", limits = c(-zmax, zmax), na.value = "grey70") 

ggsave(file="~/www/test.pdf", fig, height=8, width=12)


figList = lapply(genes, function(gene){
  df %>%
  filter(ID == gene) %>%
  mutate(se = logFC / t) %>%
  ggplot(aes(coef, logFC,  color = -log10(pmax(1e-4, adj.P.Val)))) +
    geom_point() +
    theme_classic() +
    geom_errorbar(aes(ymin = logFC - 1.96 * se, ymax = logFC + 1.96 * se), width = 0) +
    facet_wrap(~assay) +
    theme(aspect.ratio=1) +
    coord_flip() +
    scale_color_gradient(name = bquote(-log[10] ~ FDR), low = "grey", high = "red", limits = c(0, 4)) +
    ggtitle(gene) +
    ylab(bquote(log[2] ~ fold ~ change)) +
    geom_hline(yintercept = 0, linetype = "dashed") 
})

pdf("~/www/test.pdf", height=7, width=10)
figList
dev.off()

# Forest
genes = c("NFKBIA","PIM3","IFI27", "CLU","NLRC5","PTPN2", "FADD")

fig = df %>%
  filter(ID %in% genes, assay == "cMono") %>%
  mutate(se = logFC / t) %>%
  ggplot(aes(coef, logFC,  color = -log10(pmax(1e-4, adj.P.Val)))) +
    geom_point() +
    theme_classic() +
    geom_errorbar(aes(ymin = logFC - 1.96 * se, ymax = logFC + 1.96 * se), width = 0) +
    facet_wrap(~ ID, scales="free_x") +
    theme(aspect.ratio=1) +
    coord_flip() +
    scale_color_gradient(name = bquote(-log[10] ~ FDR), low = "grey", high = "red", limits = c(0, 4)) +
    ylab(bquote(log[2] ~ fold ~ change)) +
    geom_hline(yintercept = 0, linetype = "dashed") 

ggsave(file="~/www/test.pdf", fig, height=8, width=12)

# heatmap
###########
genes = c("NLRC5","PTPN2", "FADD")

# Summarize differential expression for each coef and assay
df = coef_array %>%
  map_df(function(x) 
    topTable(fit, coef = x, number=Inf) %>%
    as_tibble %>%
    mutate(coef = factor(x, coef_array))) %>%
  mutate(assay = factor(assay, assayNames(pb))) %>%
  filter(grepl("Mono|cDC", assay)) %>%
  droplevels 

df2 = df %>%
  filter(ID %in% genes, assay == "cMono") %>%
  mutate(se = logFC / t)

zmax = 1.05*max(abs(df2$logFC))
fig1 = df2 %>%  
  ggplot(aes(coef, ID, fill=logFC)) +
    geom_tile() +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))  +
    scale_fill_gradient2("logFC", low = "blue", mid = "white", high = "red", limits = c(-zmax, zmax), na.value = "grey70") +
    geom_text(data = df2 %>% filter(adj.P.Val < 0.05), label='*', hjust=0.5, vjust=0.8) +
    coord_equal()

genes = c("NFKBIA", "CLU", "IFI27", "PIM3")

df2 = df %>%
  filter(ID %in% genes, assay == "cMono") %>%
  mutate(se = logFC / t)

zmax = 1.05*max(abs(df2$logFC))
fig2 = df2 %>%  
  ggplot(aes(coef, ID, fill=logFC)) +
    geom_tile() +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
    coord_equal() +
    scale_fill_gradient2("logFC", low = "blue", mid = "white", high = "red", limits = c(-zmax, zmax), na.value = "grey70") +
    geom_text(data = df2 %>% filter(adj.P.Val < 0.05), label='*', hjust=0.5, vjust=0.8) +
    coord_equal()

fig = plot_grid(fig1 + theme(axis.text.x =  element_blank()), fig2, ncol=1, align="h", axis="lr")


ggsave(file="~/www/test.pdf", fig, height=8, width=12)








# BOXPLOT
CT = c('ncMono', 'cMono', 'cMono.cyc', 'cMono.PLT', 'cDC' )
df = lapply(CT, function(x){
  df = extractData(res.proc, assay = x, genes=genes[1:4])
  df$assay = x
  df
  })
df = bind_rows(df) %>%
  mutate(Source = factor(Source, c("HV", "COVID_HCW_MILD", "COVID_MILD", "COVID_SEV", "COVID_CRIT", "Sepsis")))

fig = ggplot(df, aes(assay, NFKBIA, fill=Source)) +
  geom_boxplot() +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5)) + 
  ggtitle('NFKBIA')


ggsave(file="~/www/test.pdf", fig, height=8, width=12)


```



<!---
NFKBIA/PIM3/IFI27 at the same level for all COVID levels
CLU/NLRC5/PTPN2 increases with severity and Sepsis
GIGYF2 returns to normal with severity
FADD is only sepsis
--->