Kfoury_CancerCell_2021.Rmd

---
title: "Human prostate cancer bone metastases have an actionable immunosuppressive microenvironment (GSE143791)"
subtitle: '[Kfoury, et al. Cancer Cell (2021)](https://doi.org/10.1016/j.ccell.2021.09.005)'
author: "Developed by [Gabriel Hoffman](http://gabrielhoffman.github.io/)"
date: "Run on `r Sys.time()`"
documentclass: article
output: 
  html_document:
    toc: true
    smart: false
vignette: >
  %\VignetteIndexEntry{Human prostate cancer bone metastases have an actionable immunosuppressive microenvironment}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
  %\usepackage[utf8]{inputenc}
---

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

# Load packages
```{r load.packages, cache=FALSE}
library(data.table)
library(SingleCellExperiment)
library(zellkonverter)
library(stringr)
library(GSEABase)
library(dreamlet)
library(scater)
library(zenith)
library(knitr)
library(kableExtra)
library(scattermore)
library(cowplot)
library(ggplot2)
library(qvalue)
library(tidyverse)
library(RColorBrewer)
library(BiocParallel)
library(DelayedArray)
```

# Load data
```{r load.data}
# read H5AD file
file = "/sc/arion/projects/CommonMind/hoffman/scRNAseq_data/Kfoury_CancerCell_2021/Kfoury_CancerCell_2021.h5ad"
sce = readH5AD(file, use_hdf5=TRUE)   

sce$cells = factor(sce$cells, sort(levels(sce$cells)))

# create pseudobulk 
pb <- aggregateToPseudoBulk(sce,
    assay = "counts",     
    cluster_id = "cells",  
    sample_id = "ID")
```

# UMAP embedding
```{r embeding}
colData(sce) %>%
  data.frame() %>%
  ggplot(aes(Coord1, Coord2, color=cells)) + 
  geom_scattermore() +
  theme_classic() + 
  theme(aspect.ratio=1, axis.text=element_blank()) +
  scale_color_discrete(name="Cell type")
```


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

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

df = list()
df$`# Samples` = ncol(pb)
df$`# Subjects` = nlevels(sce$subject.id)
df$`# Cells` = ncol(sce)
df = unlist(df)
df = data.frame(Propert=names(df), count = df)

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

## Process pseudobulk data to estimate precision weights
```{r voom}
types = c("Benign", "Distal", "Involved", "Tumor")
pb$Status = factor(pb$Status, types)

# Normalize and apply voom/voomWithDreamWeights
form = ~ (1|subject.id) + (1|Status) + (1|subject.status) 
res.proc = processAssays( pb, 
                          form,  
                          min.cells = 10,
                          min.counts = 10,
                          BPPARAM = SnowParam(6))
```

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

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


# Variance partitioning analysis
```{r vp}
# run variance partitioning analysis
form = ~ (1|subject.id) + (1|Status) + (1|subject.status)  
vp.lst = fitVarPart( res.proc, form, BPPARAM = SnowParam(6, progressbar=TRUE))
```

```{r vp.plot, fig.height=20, fig.width=12}
# Summarize variance fractions genome-wide for each cell type
plotVarPart(sortCols(vp.lst), label.angle=60, ncol=4)   
```

# dreamlet analysis
```{r dreamlet}
# Differential expression analysis for Status within each cell type,
# evaluated on the voom normalized data 
form = ~ 0 + Status + (1|subject.id) 

# Specify contrasts to compare regression coefficients
# For example, 
# Tumor_Involved = 'StatusTumor - StatusInvolved' tests if 
# expression in Tumor is different that in Involved samples
contrasts = c(Tumor_Involved  = 'StatusTumor - StatusInvolved',
              Tumor_Distal    = 'StatusTumor - StatusDistal',
              Involved_Distal = 'StatusInvolved - StatusDistal')

# dreamlet analysis
res.dl = dreamlet(res.proc, 
                  form, 
                  contrasts=contrasts, 
                  BPPARAM = SnowParam(6, progressbar=TRUE))
```


# Volcano plots for each contrast
### Tumor_Involved  = 'StatusTumor - StatusInvolved'
```{r volcano.Tumor_Involved, fig.height=16, fig.width=10}
plotVolcano( res.dl, coef = 'Tumor_Involved', ncol=4)
```

### Tumor_Distal = 'StatusTumor - StatusDistal'
```{r volcano.Tumor_Distal, fig.height=16, fig.width=10}
plotVolcano( res.dl, coef = 'Tumor_Distal', ncol=4)
```

### Involved_Distal = 'StatusInvolved - StatusDistal'
```{r volcano.Involved_Distal, fig.height=16, fig.width=10}
plotVolcano( res.dl, coef = 'Involved_Distal', ncol=4)
```


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

names(contrasts) %>%
  map_df(function(x) 
    topTable(res.dl, 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 = names(contrasts) %>%
  map_df(function(x) 
    topTable(res.dl, coef = x, number=Inf) %>%
    as_tibble %>%
    mutate(coef = x)) %>%
  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=1, legend.position="none") +
    scale_x_continuous(limits=c(0,ymax), expand=c(0,0)) +
    xlab("# genes expressed") +
    ylab("Cell type") +
    facet_wrap(~ coef)

ymax = 1.05*max(df$nDE)
fig2 = ggplot(df, aes(nDE, assay, fill=assay)) + 
    geom_bar(stat="identity") + 
    theme_classic() +
    theme(aspect.ratio=1, 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)

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

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


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

# Run zenith gene set analysis on result of dreamlet
res_zenith.Tumor_Involved = zenith_gsa(res.dl, coef = 'Tumor_Involved', go.gs)
res_zenith.Tumor_Distal = zenith_gsa(res.dl, coef = 'Tumor_Distal', go.gs)
res_zenith.Involved_Distal = zenith_gsa(res.dl, coef = 'Involved_Distal', go.gs)
```



### Tumor_Involved
```{r zenith.plot_Tumor_Involved, fig.height=18, fig.width=12}
plotZenithResults(res_zenith.Tumor_Involved, 5, 3)
```

### Tumor_Distal
```{r zenith.plot_Tumor_Distal, fig.height=18, fig.width=12}
plotZenithResults(res_zenith.Tumor_Distal, 5, 3)
```

### Involved_Distal
```{r zenith.plot_Involved_Distal, fig.height=22, fig.width=13}
plotZenithResults(res_zenith.Involved_Distal, 5, 3)
```



Tumor vs Involved
Monocytes show increase MHC class Involved_Distal
  and interferon signalling
Increased translation in tumor inflammatory monocytes (TIMs) 
Erythroid show increased antigen presentation
B-cell progenentors show icnrease positive regulation of cell growth

CTL exhaustion, loss of mitochrondrial translation


# Combine figures into panels
```{r combine, fig.height=5, fig.width=12}
CT = "Mono3"

fig.embed = colData(sce) %>%
  data.frame() %>%
  ggplot(aes(Coord1, Coord2, color=ifelse(cells==CT, CT, ''))) + 
  geom_scattermore() +
  theme_classic() + 
  scale_color_manual(name='', values=c("grey", "orange")) +
  theme(aspect.ratio=1, 
    axis.text=element_blank(),
    axis.ticks=element_blank(), 
    legend.position="bottom") +
  xlab("Coord 1") +
  ylab("Coord 2")

thm = theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), aspect.ratio=1)
col = c(brewer.pal(n=3, "Set1"), "grey85")

fig.voom = plotVoom( res.proc[[CT]]) + theme_classic() + thm
fig.vp = plotVarPart( sortCols(vp.lst), label.angle=20, assays=CT, col=col) + thm 

genes = c('CD52', 'IDH3G', "CCNL1", 'NAMPT', 'CHMP1B', 'IFNAR1', 'VIM', 'UFC1', 'NOTCH2')

fig.percent = plotPercentBars( sortCols(vp.lst), assays = CT, genes = genes, col=col) + theme(aspect.ratio=1, legend.position="bottom")

fig.vol = plotVolcano( res.dl[[CT]], coef = 'Tumor_Involved') + ggtitle("Tumor vs Involved")
  
data = extractData(res.proc, CT)

fig.strat0 = plotStratify( CCNL1 ~ Status, data, legend=FALSE, x.labels=TRUE, sort=FALSE) + 
  theme(aspect.ratio=2, 
    plot.title = element_text(hjust = 0.5),
    axis.text.x=element_text(hjust=1, vjust=1, angle=30)) +
  ggtitle('CCNL1') +
  ylab(bquote(log[2]~CPM)) +
  xlab('')

fig.strat1 = plotStratify( VIM ~ Status, data, legend=FALSE, x.labels=TRUE, sort=FALSE) + 
  theme(aspect.ratio=2, 
    plot.title = element_text(hjust = 0.5),
    axis.text.x=element_text(hjust=1, vjust=1, angle=30)) +
  ggtitle('VIM') +
  ylab(bquote(log[2]~CPM)) +
  xlab('')
 
fig.strat2 = plotStratify( CD52 ~ subject.id, data, legend=FALSE, x.labels=TRUE) + 
  theme(aspect.ratio=1, 
    plot.title = element_text(hjust = 0.5),
    axis.text.x=element_text(hjust=1, vjust=1, angle=70)) +
  ggtitle('CD52') +
  ylab(bquote(log[2]~CPM)) +
  xlab('')

plot_grid(fig.strat0, fig.strat1, fig.strat2, 
  nrow=1, ncol=3, labels=LETTERS[1:10], align="hv", axis="t")
```


```{r combine2, fig.height=7, fig.width=12}
gs = c('GO0022624', 'GO1905369', 'GO0007186', 'GO0050920', 'GO0042613', 'GO0044183', 'GO1903561', 'GO0023026', 'GO0042026')

res_zenith = res_zenith.Tumor_Involved

idx = grep(paste0(gs, collapse="|"), res_zenith$Geneset)

res_zenith$assay = factor(res_zenith$assay)
 
fig.zenith = plotZenithResults(res_zenith[idx,], 5, 3, transpose=FALSE) + theme(axis.text.x=element_text(size=7), legend.position="right",legend.key.size = unit(.4, 'cm'))
  
lst = list(fig.embed, fig.voom, fig.vp, fig.percent, fig.strat0, fig.strat1, fig.vol) 
plot_grid(plotlist=lst, nrow=2, ncol=4, labels=LETTERS[1:10], align="hv", axis="t") 
```

# results counts
```{r counts}
# variance partition results
vp.lst %>% 
  as_tibble %>%
  filter(assay == "Mono3") %>%
  summarize(count = sum(Status > 0.1), ngenes = length(Status))

# count DE genes
names(contrasts) %>%
  map_df(function(x) 
    topTable(res.dl, coef = x, number=Inf) %>%
    as_tibble %>%
    mutate(coef = x)) %>%
    filter(assay == "Mono3") %>%
    group_by(coef) %>%
    summarize(count = sum(adj.P.Val < 0.05), ngenes = length(adj.P.Val))
```





```{r combine3, fig.height=5, fig.width=12}
res_zenith = res_zenith.Tumor_Involved 
idx = grep(paste0(gs, collapse="|"), res_zenith$Geneset)

res_zenith$assay = factor(res_zenith$assay)
fig.zenith1 = plotZenithResults(res_zenith[idx,], 5, 3, sortByGeneset=FALSE) + 
  theme(axis.text.x=element_text(size=7, angle=90),
    legend.position="right",
    plot.title = element_text(hjust = 0.5),
    legend.key.size = unit(.4, 'cm'))+
  ggtitle('Tumor_Involved')
  
res_zenith = res_zenith.Tumor_Distal
idx = grep(paste0(gs, collapse="|"), res_zenith$Geneset)

res_zenith$assay = factor(res_zenith$assay)
fig.zenith2 = plotZenithResults(res_zenith[idx,], 5, 3, sortByGeneset=FALSE) + 
  theme(axis.text.x=element_text(size=7, angle=90),
    legend.position="right",
    plot.title = element_text(hjust = 0.5),
    legend.key.size = unit(.4, 'cm')) +
  ggtitle('Tumor_Distal')

plot_grid(fig.zenith1, fig.zenith2, ncol=1)
```


### Highlight specific genesets

```{r plotGeneHeatmap, fig.height=6, fig.width=8}
gs1 = "GO0044183: protein folding chaperone"
plotGeneHeatmap( res.dl, coef="Tumor_Involved", genes=sort(geneIds(go.gs[[gs1]])), transpose=TRUE) + 
  ggtitle(gs1) + 
  theme(legend.position = "bottom",
    axis.text.x=element_text(size=10, angle=60),
    axis.text.y=element_text(size=10)) + 
  xlab('') + ylab('') 
```

```{r plotGeneHeatmap2, fig.height=3, fig.width=8}
gs1 = "GO0044183: protein folding chaperone"
plotGeneHeatmap( res.dl, coef="Tumor_Involved", genes=sort(geneIds(go.gs[[gs1]])), assay=CT, transpose=TRUE) + 
  ggtitle(gs1) + 
  theme(legend.position = "bottom",
    axis.text.x=element_text(size=10, angle=45),
    axis.text.y=element_text(size=10)) + 
  xlab('') + ylab('') 
```






# Session info
<details>
```{r sessionInfo}
sessionInfo()
```
</details>

```{r exit123, cache=FALSE, echo=FALSE}
knitr::knit_exit()
```







# Cell composition
```{r composition.plot}
plotCellComposition(pb[,order(pb$Status)])     
```

```{r fracs, fig.height=12}
frac = apply(cellCounts(pb), 1, function(x) x / sum(x))

df = reshape2::melt(frac)
df = merge(df, colData(pb), by.x='Var2', by.y="row.names")

ggplot(as.data.frame(df), aes(Var1, value, fill=Status)) + geom_boxplot() + theme_classic() + theme(aspect.ratio=3) + ylab("Fraction") + xlab('') + coord_flip()
```