SARS_micro_array_analysis.Rmd

---
title: "Cellular Immuneregulation upon SARS-Cov2 Infection"
author: "Florian Heigwer"
date: "`r Sys.Date()`"
abstract: >
 Lorem Ipsum
bibliography: references.bib
output: BiocStyle::pdf_document
vignette: >
  %\VignetteIndexEntry{Vignette Title}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

# About

This document contains code to analyse a micro-array experiment performed by Chrisptopher Neufeldt and collegues at the University Clinic Heidelberg "Bartenschlager Group".

The data comprises multiple Affymetrix chips performed by the DKFZ? core facility.

Purpose is to find genes differentially expressed upon infection of epithelial lung cancer cell lines with SARS-Cov2 virus isolates.

Cell lines:

  A549, cells ACE2 overexpression transduced
  Calu-3, cells wild type
  
Virus:
  
  Sars-Cov2, MOI ?
  
Data was produced in replicates sampling the treatment after 4, 8, 12 and 24 hours.

# Dependencies

We load a number of packages whose functions are needed throughout the analysis


```{r message=F, warning=F, results='hide'}
library(limma)
library(ggrepel)
#library(lumi)
library(patchwork)
library(ggrastr)
library(pheatmap)
library(fgsea)
library(tidyverse)
#library(affy)
#library(affydata)
#library(affycoretools)
library(gplots)
library(RColorBrewer)
library(eulerr)

```

#Lets define some themes

This we do as a quality of life step to give all figures that we produce a common look and feel.

```{r theme, include=FALSE}

theme_b110 <- function(){
  theme_classic() +
  theme(
    axis.text = element_text(size = 10), 
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 12,hjust = 0.5,face = "bold"),
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 10),
    legend.position = "bottom"
    )
}

sgi_blue    = '#5087C8'
sgi_yellow1 = '#F2EE35'
sgi_yellow2 = '#FED98E'
b110_grey   = '#808080'
b110_grey_light   = '#909090'
b110_transparent_black = alpha('#000000',0.5)
google_red = '#dd4b39'
google_green = '#0F9D58'
google_yellow = '#F4B400'
google_blue = '#4285F4'

```

#Read in data

Data is used as provided by the DKFZ core-facility. Data is assumed to be quantil normalized per experiment not already log-transformed and summarized by averaging on gene level. Gene - Level View is presented.

Excel files were collected and saved as tab-delimited files to manually check integrity and column names.

Of note: In initial data Calu2-24h-mock.mean was wrongly labeled as Calu1-12h-infected.mean. I corrected this manually based on missing value and doubled naming plus its position in the data table

```{r}

calu_data <- read_delim("transcriptome/all.qnorm.gv.calu3.txt",delim = "\t")

a549_data <- read_delim("transcriptome/all.qnorm.gv.a549.txt",delim = "\t")

joined_data <- calu_data %>% left_join(a549_data,by = "ProbeID")

```

prepare data for GEO submission

```{r}

normed_data <- joined_data %>% 
  select(ProbeID,ends_with(c(".mean",".sd",".p",".nbeads"))) %>%
  filter(grepl("^T",ProbeID))

colnames(normed_data) <- gsub("Calu(\\d)(.+)","Calu3-\\1\\2",gsub("\\.","-",colnames(normed_data)))

normed_data %<>% 
  gather(dat,value,-ProbeID) %>%
  separate(dat,c("cell_line","replicate","time","treatment","type"),remove = FALSE) %>%
  mutate(time =factor(time,levels = c("4h","8h","12h","24h"))) %>%
  mutate(treatment  =factor(treatment ,levels = c("mock","infected"))) %>% 
  arrange(cell_line,replicate,time,treatment) %>%
  mutate(dat = factor(dat,levels = unique(dat))) %>%
  select(ProbeID,dat,value) %>%
  spread(dat,value) %>%
  rename(ID_REF=ProbeID)

suspects <- colnames(normed_data)

suspects = gsub(".*-sd","sd",gsub(".*-p","p-value",gsub(".*-nbeads","nbeads",gsub("(.*)-mean","\\1",suspects))))

colnames(normed_data) <- suspects

write_csv(normed_data,"transcriptome/GEO/raw_CEL_files/all.qnorm.gv.corrected.csv.gz")
write_delim(normed_data,"transcriptome/GEO/raw_CEL_files/all.qnorm.gv.corrected.txt",delim = "\t")

raw_data <- read_delim("transcriptome/GEO/raw_CEL_files/affy_chips01_all.raw.txt.geo.GPL23159.txt",delim = "\t") %>% 
  select(-stop) %>%
  filter(grepl("^T",ID_REF))

mappings <- read_delim("transcriptome/GEO/raw_CEL_files/sample description.csv",delim = ";") 

colnames(raw_data) <- gsub(".*(CN\\w{2}).* (\\w+)$","\\1_\\2",  colnames(raw_data))

raw_data %<>%
  gather(dat,value,-ID_REF) %>%
  separate(dat,c("sampleID","type"),remove = FALSE) %>%
  left_join(mappings)  %>%
  mutate(time =factor(time,levels = c("4h","8h","12h","24h"))) %>%
  mutate(treatment  =factor(treatment ,levels = c("mock","infected"))) %>% 
  arrange(cell_line,replicate,time,treatment) %>%
  unite("clname",cell_line,replicate,time,treatment,type,sep="-") %>%
  mutate(clname = factor(clname,levels = unique(clname))) %>%
  select(ID_REF,clname,value) %>%
  spread(clname,value)

suspects <- colnames(raw_data)

suspects = gsub(".*-pvalue","pvalue",gsub(".*-nbeads","nbeads",gsub("(.*)-mean","\\1",suspects)))

colnames(raw_data) <- suspects

write_csv(raw_data,"transcriptome/GEO/raw_CEL_files/all.raw.corrected.csv.gz")
write_delim(raw_data,"transcriptome/GEO/raw_CEL_files/all.raw.corrected.txt",delim = "\t")

```


## Separate annotation and data

We assume the data is indeed already normalized and non-log transformed. So we do the log-transform of the expression data here.

```{r}

expression_data <- joined_data %>% select(ProbeID,ends_with("mean")) %>% mutate_if(is.numeric,log) %>% column_to_rownames("ProbeID") 

colnames(expression_data) <- gsub("Calu(\\d)","Calu3-\\1",colnames(expression_data)) 

bead_count <- joined_data %>% select(ends_with("nbeads"),ProbeID) %>% column_to_rownames("ProbeID")

annotation_data <- joined_data %>% select(Symbol=Symbol...89,Name=...90.x,ProbeID) 

```

##Overview density plots

So we vizualize a common density plot.

```{r}

full <- expression_data %>%
  rownames_to_column("ProbeID") %>%
  gather(sample,value,-ProbeID) %>%
  ggplot(aes(value,col=sample,fill=NULL)) +
    geom_density() +
    theme_classic()

full

ggsave(plot = full,filename = "plots/all_distributions.pdf")
```

Indeed the data seems sufficiently normalized as all dsitributions looks very similar and overlap to the extend expected from normalized data.

# PCA analysis

By PCA we can see if treatments or time points are driving the variation. Thus one should always see replicated close together, samples lined up by increasing time, and treatments separated in different corners.

```{r}

pc <- prcomp(expression_data %>% drop_na() )

pc <- pc$rotation  %>% as.data.frame() %>% rownames_to_column("sample") %>% as_tibble %>% separate(sample,c("cell","repl","time","treatment"))

p1 <- pc %>%
  ggplot(aes(x=PC1,y=PC2)) +
    geom_point(aes(col=time,shape=cell),size=4) +
    theme_b110() +
    ggtitle("by time and cell")

p2 <- pc %>%
  ggplot(aes(x=PC1,y=PC2)) +
    geom_point(aes(col=treatment,shape=cell),size=4) +
    theme_b110()+
    ggtitle("by treatment and cell")

pc_calu <- prcomp(expression_data %>% drop_na() %>% select(starts_with("Calu")) )

pc_calu <- pc_calu$rotation  %>% as.data.frame() %>% rownames_to_column("sample") %>% as_tibble %>% separate(sample,c("cell","repl","time","treatment"))

p3 <- pc_calu %>%
  ggplot(aes(x=PC1,y=PC2)) +
    geom_point(aes(col=time,shape=treatment),size=4) +
    theme_b110() +
    ggtitle("Calu - 3")


pc_a549<- prcomp(expression_data %>% drop_na() %>% select(starts_with("A549")) )

pc_a549 <- pc_a549$rotation  %>% as.data.frame() %>% rownames_to_column("sample") %>% as_tibble %>% separate(sample,c("cell","repl","time","treatment"))

p4 <- pc_a549 %>%
  ggplot(aes(x=PC1,y=PC2)) +
    geom_point(aes(col=time,shape=treatment),size=4) +
    theme_b110() +
    ggtitle("A549")


p1 + p2 + p3 + p4

ggsave("plots/PCA_analysis.pdf")
```

As the cell line marks the biggest difference between all samples we will analyse them separately. 

Likewise we can observe that treatment and time nicely separate replicates and sample in the Calu-3 samples but not in the A549 samples. There is a separation by time in PC3 of A549, but not by treatment.

#Differential gene expression analysis

##Let's call DE genes for Calu -3 only

Next we can run limma's moderated t-test to find differentially expressed genes between virus and mock treated samples.

For simplicity of the gene expression model we fused the 4 h and 8 h timepoints to "early" and the remaining time points to "late".

```{r, results='hide', warning=F, message=F}

expression_calu <- expression_data %>% select(starts_with("Calu"))
## vector specifying treatment per sample
treatment <- factor(ifelse(grepl('mock', colnames(expression_calu)), 'mock', 'infected'))

#or simply define the first two as early and the last two as late
time_agg <- factor(
  ifelse(grepl('24h', colnames(expression_calu)), 'late',
         ifelse(grepl('8h', colnames(expression_calu)), 'early', 
                ifelse(grepl('12h', colnames(expression_calu)), 
                'late', 'early'))),levels = c("early","late"))

## model matrix for DGE analysis
mm <- model.matrix(~treatment + time_agg) #+ cell:treatment 

colnames(mm)=gsub(':','_',colnames(mm))


## first we need to define the contrasts
contr <- makeContrasts(
  treat=treatmentmock,
  time=time_agglate,
  levels=mm)

## perform DGE analysis
fit_affy <- lmFit(expression_calu, mm)
fit_affy <- eBayes(fit_affy)

fit_affy2 <- contrasts.fit(fit_affy, contr)
fit_affy2 <- eBayes(fit_affy2)

# results Trametinib rescue experiments
calu_model <- map(1:2, function(i){
  topTable(fit_affy2, n=Inf, coef=i) %>% 
    rownames_to_column('ProbeID') %>% tbl_df %>% 
    left_join(distinct(annotation_data)) %>%
    return()
})
names(calu_model) <-colnames(contr)


```

##  CALU Lets produce overview plots

Time component

```{r}

this_data <- calu_model$time %>%
  mutate(significant = if_else(adj.P.Val<=0.1,"significant","non-significant"))

a <- this_data %>%
  ggplot(aes(x=AveExpr,y=logFC,col=significant)) +
    geom_point(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_smooth(data = subset(this_data,significant=="non-significant"),col="black") +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue)) + 
    ggtitle("Time differential MA-plot")

b <- this_data %>%
  ggplot(aes(x=logFC,y=-log10(P.Value),col=significant))  +
    geom_point(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_vline(xintercept = 0,lty=2,col=b110_grey_light) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue))+ 
    ggtitle("Time differential Volcano-plot")

a+b

ggsave("plots/calu_time_MA_volcano.pdf")
```

Time component rastered

```{r}

this_data <- calu_model$time %>%
  mutate(significant = if_else(adj.P.Val<=0.1,"significant","non-significant"))

a <- this_data %>%
  ggplot(aes(x=AveExpr,y=logFC,col=significant)) +
    geom_point_rast(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_smooth(data = subset(this_data,significant=="non-significant"),col="black") +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue)) + 
    ggtitle("Time differential MA-plot")

b <- this_data %>%
  ggplot(aes(x=logFC,y=-log10(P.Value),col=significant))  +
    geom_point_rast(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_vline(xintercept = 0,lty=2,col=b110_grey_light) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue))+ 
    ggtitle("Time differential Volcano-plot")

a+b

ggsave("plots/calu_time_MA_volcano_rastered.pdf")
```

Treatment component

High logFC means upregulated

```{r}

  this_data <- calu_model$treat %>%
  mutate(significant = if_else(adj.P.Val<=0.1,"significant","non-significant"))

a <- this_data %>%
  ggplot(aes(x=AveExpr,y=-logFC,col=significant)) +
    geom_point(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_smooth(data = subset(this_data,significant=="non-significant"),col="black") +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue)) + 
    ggtitle("Infection differential MA-plot")

b <- this_data %>%
  ggplot(aes(x=-logFC,y=-log10(P.Value),col=significant))  +
    geom_point(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_vline(xintercept = 0,lty=2,col=b110_grey_light) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue))+ 
    ggtitle("Infection differential Volcano-plot")

a+b


ggsave("plots/calu_treat_MA_volcano.pdf")
```


Treatment component rastered

High logFC means upregulated

```{r}

  this_data <- calu_model$treat %>%
  mutate(significant = if_else(adj.P.Val<=0.1,"significant","non-significant"))

a <- this_data %>%
  ggplot(aes(x=AveExpr,y=-logFC,col=significant)) +
    geom_point_rast(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_smooth(data = subset(this_data,significant=="non-significant"),col="black") +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue)) + 
    ggtitle("Infection differential MA-plot")

b <- this_data %>%
  ggplot(aes(x=-logFC,y=-log10(P.Value),col=significant))  +
    geom_point_rast(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_vline(xintercept = 0,lty=2,col=b110_grey_light) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue))+ 
    ggtitle("Infection differential Volcano-plot")

a+b


ggsave("plots/calu_treat_MA_volcano_rast.pdf")
```

This code can be used to pull out and vizualize some of the results

```{r}

candidate <- "TC0200016651.hg.1"

a<-expression_calu %>%
  rownames_to_column("ProbeID") %>%
  gather(sample,value,-ProbeID) %>%
  filter(ProbeID==candidate) %>%
  separate(sample,c("cell","rep","time","treat"),sep = "-") %>%
  #mutate(time=factor(time,levels = c("4h","8h","12h","24h"))) %>%
  mutate(time=if_else(time %in% c("4h","8h"),1,2)) %>%
  ggplot(aes(x=time,y=value,col=treat,grp=treat)) +
    geom_point() +
    geom_smooth(method = "lm",se = F) +
    theme_classic() +
    ggtitle("Calu3")

b<-expression_calu %>%
  rownames_to_column("ProbeID") %>%
  gather(sample,value,-ProbeID) %>%
  filter(ProbeID==candidate) %>%
  separate(sample,c("cell","rep","time","treat"),sep = "-") %>%
  mutate(time=factor(time,levels = c("4h","8h","12h","24h"))) %>%
   # mutate(time=if_else(time %in% c("4h","8h"),"early","late"),treat=factor(treat,levels = c("mock.mean","infected.mean"))) %>%
  ggplot(aes(x=treat,y=value,col=time,grp=treat)) +
    geom_boxplot() +
    geom_point(position = position_dodge(0.9))+
    theme_b110() +
    ggtitle("Calu3")

a + b


```


#Then we vizualise the data as heatmap and the top differential genes als as a fucntion of time

# Time Heatmap

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

# here we threshold the data simply in upregulated and down regulated according to foldchange and sort by FDR
# for enrichment of functional annotation clustering it make further sense to threshold by fold change and FDR

res <- calu_model$treat %>%
  mutate(Symbol = gsub(pattern = "\\s",replacement = "",x = Symbol)) %>% 
  na.omit() %>% 
  distinct() %>% 
  group_by(Symbol) %>% 
  summarize(logFC=-mean(logFC)) %>%
  ungroup() %>%
  arrange(-(abs(logFC))) %>%
  head(n=100) %>%
  left_join(annotation_data %>% mutate(Symbol = gsub(pattern = "\\s",replacement = "",x = Symbol))) %>%
  left_join(expression_calu %>% 
  rownames_to_column(var = "ProbeID")) %>%
  gather(sample,value,-Symbol,-logFC,-Name,-ProbeID) %>% 
  separate(sample,c("cell","rep","time","treat"),sep = "-") %>%
  mutate(time=factor(time,levels = c("4h","8h","12h","24h")))


# produce a time seriesed heatmap
genes_to_order <-res %>%
  group_by(rep,time,ProbeID) %>% mutate(value=value-mean(value[treat=="mock.mean"],na.rm=T)) %>% 
  group_by(rep,treat,Symbol,time) %>% summarise(value=mean(x = value,na.rm=T)) %>%
  ungroup() %>% 
  filter(treat=="infected.mean") %>% 
  group_by(Symbol,time) %>% summarise(value=mean(x = value,na.rm=T)) %>% 
  spread(time,value) %>% column_to_rownames("Symbol") 

gene_order <- rownames(genes_to_order )[genes_to_order %>% dist() %>% hclust() %>% .$order] %>% gsub(pattern = "ENSG.+_(.+)$",replacement = "\\1")

p2 <- res %>% 
  group_by(rep,time,ProbeID) %>% mutate(value=value-mean(value[treat=="mock.mean"],na.rm=T)) %>% 
  group_by(rep,treat,Symbol,time) %>% summarise(value=mean(x = value,na.rm=T)) %>%
  ungroup() %>% 
  filter(treat=="infected.mean",Symbol %in% rownames(genes_to_order)) %>% 
  drop_na() %>%
    mutate(Symbol=factor(Symbol,levels = gene_order )) %>% 
    ggplot(aes(x=time,y=Symbol,fill=value)) + 
    geom_tile() + 
    scale_fill_gradient2(low = "blue",mid = "white",high = "red",midpoint = 0) +
    theme_b110()

ggsave(p2,filename = "plots/gene_expression_over_time_heatmap_reltoMockpertime_top100.pdf",height = 16)

#plot the top 12 up regulated

res <- calu_model$treat %>%
  mutate(Symbol = gsub(pattern = "\\s",replacement = "",x = Symbol)) %>% 
  na.omit() %>% 
  distinct() %>% 
  group_by(Symbol) %>% 
  summarize(t=-mean(t)) %>%
  ungroup() %>%
  #arrange(-(t)) %>%
  #head(.,n=12) %>%
  filter(grepl("NUAK2",Symbol)) %>%
  left_join(annotation_data %>% mutate(Symbol = gsub(pattern = "\\s",replacement = "",x = Symbol))) %>%
  left_join(expression_calu %>% 
  rownames_to_column(var = "ProbeID")) %>%
  gather(sample,value,-Symbol,-t,-Name,-ProbeID) %>% 
  separate(sample,c("cell","rep","time","treat"),sep = "-") %>%
  #mutate(treatment=factor(treatment,levels = c("Cov2","Mock"))) %>%
  mutate(time=
  ifelse(grepl('24h',time), 24,
         ifelse(grepl('8h', time), 8, 
                ifelse(grepl('12h', time),12, 4))))#%>%
  #mutate(time=factor(time,levels = c("4h","8h","12h","24h")))

p1 <- res %>% 
  ggplot(aes(x=time,y = value,col=treat))+
    geom_point() +
    geom_smooth() +
    facet_wrap(~Symbol,scales = "free") +
    theme_b110() +
    scale_color_manual(values = c("red",b110_grey))


res <- calu_model$treat %>%
  mutate(Symbol = gsub(pattern = "\\s",replacement = "",x = Symbol)) %>% 
  na.omit() %>% 
  distinct() %>% 
  group_by(Symbol) %>% 
  summarize(t=-mean(t)) %>%
  ungroup() %>%
  arrange((t)) %>%
  head(.,n=12) %>%
  left_join(annotation_data %>% mutate(Symbol = gsub(pattern = "\\s",replacement = "",x = Symbol))) %>%
  left_join(expression_calu %>% 
  rownames_to_column(var = "ProbeID")) %>%
  gather(sample,value,-Symbol,-t,-Name,-ProbeID) %>% 
  separate(sample,c("cell","rep","time","treat"),sep = "-") %>%
  #mutate(treatment=factor(treatment,levels = c("Cov2","Mock"))) %>%
  mutate(time=
  ifelse(grepl('24h',time), 24,
         ifelse(grepl('8h', time), 8, 
                ifelse(grepl('12h', time),12, 4))))#%>%
  #mutate(time=factor(time,levels = c("4h","8h","12h","24h")))

p2 <- res %>% 
  ggplot(aes(x=time,y = value,col=treat))+
    geom_point() +
    geom_smooth() +
    facet_wrap(~Symbol,scales = "free") +
    theme_b110() +
    scale_color_manual(values = c("red",b110_grey))

ggsave(p1,filename = "plots/top_upregulated_overtime_calu3.pdf")

ggsave(p2,filename = "plots/top_downregulated_overtime_calu3.pdf")

```


##Let's call DE genes for A549 -3 only

Next we can run limma's moderated t-test to find differentially expressed genes between Trametinib and DMSO treated samples.

```{r, results='hide', warning=F, message=F}

expression_calu <- expression_data %>% select(starts_with("A549"))
## vector specifying treatment per sample
treatment <- factor(ifelse(grepl('mock', colnames(expression_calu)), 'mock', 'infected'))

#or simply define the first two as early and the last two as late
time_agg <- factor(
  ifelse(grepl('24h', colnames(expression_calu)), 'late',
         ifelse(grepl('8h', colnames(expression_calu)), 'early', 
                ifelse(grepl('12h', colnames(expression_calu)), 
                'late', 'early'))),levels = c("early","late"))

## model matrix for DGE analysis
mm <- model.matrix(~treatment + time_agg) #+ cell:treatment 

colnames(mm)=gsub(':','_',colnames(mm))


## first we need to define the contrasts
contr <- makeContrasts(
  treat=treatmentmock,
  time=time_agglate,
  levels=mm)

## perform DGE analysis
fit_affy <- lmFit(expression_calu, mm)
fit_affy <- eBayes(fit_affy)

fit_affy2 <- contrasts.fit(fit_affy, contr)
fit_affy2 <- eBayes(fit_affy2)

# results Trametinib rescue experiments
a549_model <- map(1:2, function(i){
  topTable(fit_affy2, n=Inf, coef=i) %>% 
    rownames_to_column('ProbeID') %>% tbl_df %>% 
    left_join(distinct(annotation_data)) %>%
    return()
})
names(a549_model) <-colnames(contr)


```

## A549 Lets produce overview plots

Time component

```{r}

this_data <- a549_model$time %>%
  mutate(significant = if_else(adj.P.Val<=0.2,"significant","non-significant"))

a <- this_data %>%
  ggplot(aes(x=AveExpr,y=logFC,col=significant)) +
    geom_point(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_smooth(data = subset(this_data,significant=="non-significant"),col="black") +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue)) + 
    ggtitle("Time differential MA-plot")

b <- this_data %>%
  ggplot(aes(x=logFC,y=-log10(P.Value),col=significant))  +
    geom_point(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_vline(xintercept = 0,lty=2,col=b110_grey_light) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue))+ 
    ggtitle("Time differential Volcano-plot")

a+b

ggsave("plots/a549_time_MA_volcano.pdf")
```

Time component rastered

```{r}

this_data <- a549_model$time %>%
  mutate(significant = if_else(adj.P.Val<=0.2,"significant","non-significant"))

a <- this_data %>%
  ggplot(aes(x=AveExpr,y=logFC,col=significant)) +
    geom_point_rast(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_smooth(data = subset(this_data,significant=="non-significant"),col="black") +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue)) + 
    ggtitle("Time differential MA-plot")

b <- this_data %>%
  ggplot(aes(x=logFC,y=-log10(P.Value),col=significant))  +
    geom_point_rast(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_vline(xintercept = 0,lty=2,col=b110_grey_light) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue))+ 
    ggtitle("Time differential Volcano-plot")

a+b

ggsave("plots/a549_time_MA_volcano_rast.pdf")
```

Treatment component

High logFC means upregulated

```{r}

  this_data <- a549_model$treat %>%
  mutate(significant = if_else(adj.P.Val<=0.2,"significant","non-significant"))

a <- this_data %>%
  ggplot(aes(x=AveExpr,y=-logFC,col=significant)) +
    geom_point(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_smooth(data = subset(this_data,significant=="non-significant"),col="black") +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue)) + 
    ggtitle("Infection differential MA-plot")

b <- this_data %>%
  ggplot(aes(x=-logFC,y=-log10(P.Value),col=significant))  +
    geom_point(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_vline(xintercept = 0,lty=2,col=b110_grey_light) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue))+ 
    ggtitle("Infection differential Volcano-plot")

a+b

ggsave("plots/a549_treat_MA_volcano.pdf")
```

Treatment component rastered

High logFC means upregulated

```{r}

  this_data <- a549_model$treat %>%
  mutate(significant = if_else(adj.P.Val<=0.2,"significant","non-significant"))

a <- this_data %>%
  ggplot(aes(x=AveExpr,y=-logFC,col=significant)) +
    geom_point_rast(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_smooth(data = subset(this_data,significant=="non-significant"),col="black") +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue)) + 
    ggtitle("Infection differential MA-plot")

b <- this_data %>%
  ggplot(aes(x=-logFC,y=-log10(P.Value),col=significant))  +
    geom_point_rast(data = subset(this_data,significant=="non-significant"),alpha=0.3) +
    geom_point(data = subset(this_data,significant=="significant")) +
    geom_vline(xintercept = 0,lty=2,col=b110_grey_light) +
    geom_label_repel(data = this_data %>% filter(significant=="significant") %>% arrange(-abs(logFC)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue))+ 
    ggtitle("Infection differential Volcano-plot")

a+b

ggsave("plots/a549_treat_MA_volcano_rast.pdf")
```

#Interpretation

Gene set annotations are downloaded from MSigDB
https://www.gsea-msigdb.org/gsea/msigdb/collections.jsp#H
29.05.2020: 14:42 CET
C5 : GO Bioprocess

C7 : Immunology

h.all : Hallmark gene sets

Citation for hallmark collection:

Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015 Dec 23;1(6):417-425. 
 

Signature enrichment

```{r}

# Load the pathways into a named list
pathways.hallmark <- gmtPathways("data/h.all.v7.1.symbols.gmt")

pathways.immune <- gmtPathways("data/c7.all.v7.1.symbols.gmt")

pathways.go <- gmtPathways("data/c5.bp.v7.1.symbols.gmt")



calu_treat_res <- calu_model$treat %>%
  select(stat=t,Symbol) %>%
  mutate(Symbol = gsub(pattern = "\\s",replacement = "",x = Symbol)) %>% 
  na.omit() %>% 
  distinct() %>% 
  group_by(Symbol) %>% 
  summarize(stat=-mean(stat)) %>%
  ungroup() %>%
  deframe()

#GSEA of hallmark sets

fgsea_hallmark <- fgsea(pathways=pathways.hallmark, stats=calu_treat_res, nperm=1000)
fgsea_hallmarktdy <- fgsea_hallmark %>%
  as_tibble() %>%
  arrange(desc(NES))

fgsea_hallmarktdy %>% 
  arrange(padj) %>%
  head(.,40) %>%
  ggplot(aes(reorder(pathway, NES), NES)) +
  geom_col(aes(fill=padj<0.05)) +
  coord_flip() +
  labs(x="Pathway", y="Normalized Enrichment Score",
       title="Hallmark pathways NES from GSEA") + 
  theme_b110()+
  scale_fill_manual(values = c(b110_grey,google_blue))

ggsave("plots/enriched_signatures.pdf")

#GSEA of GO-term sets

fgsea_go <- fgsea(pathways=pathways.go, stats=calu_treat_res, nperm=1000)
fgsea_gotdy <- fgsea_go %>%
  as_tibble() %>%
  arrange(desc(NES))

fgsea_gotdy %>% 
  arrange(padj) %>%
  head(.,40) %>%
  ggplot(aes(reorder(pathway, NES), NES)) +
  geom_col(aes(fill=padj<0.05)) +
  coord_flip() +
  labs(x="Pathway", y="Normalized Enrichment Score",
       title="Hallmark pathways NES from GSEA") + 
  theme_minimal()

#GSEA of imuunity experimental sets

fgsea_immune <- fgsea(pathways=pathways.immune, stats=calu_treat_res, nperm=1000)
fgsea_immunetdy <- fgsea_immune %>%
  as_tibble() %>%
  arrange(desc(NES))

fgsea_immunetdy %>% 
  arrange(padj) %>%
  head(.,40) %>%
  ggplot(aes(reorder(pathway, NES), NES)) +
  geom_col(aes(fill=padj<0.05)) +
  coord_flip() +
  labs(x="Pathway", y="Normalized Enrichment Score",
       title="Hallmark pathways NES from GSEA") + 
  theme_minimal()

```


### GSEA analysis

We perform gene set enrichment analysis using the Broad Institute's [GSEA](http://software.broadinstitute.org/gsea/index.jsp) [@pmid17644558]. An R version of the algorithm is implemented in the `fgsea` algorithm [@Sergushichev060012], which we use for analysis.

We want to visualize the results as a barcode plot. `fgsea` already implements a nice barcode plot, which we cusotomize a bit to adapt it according to our expectations.

Code for the abrcode plot was inspired by [@pmid31097693] . c Benedikt Rauscher

```{r, results='hide', warning=F, message=F}
custom_barcode_plot <- function(df, sig){
  ## named vector of gene-level stats
  stat_vector <- setNames(df$t, df$Symbol)
  ## genes in signature
  sig_genes <- pathways.hallmark[[sig]]
  
  ## generate barcode plot
  bc_plot <- plotEnrichment(sig_genes, stat_vector)
  
  ## remove unwanted layers
  bc_plot$layers <- list()
  
  ## add barcode at the bottom
  lowest_pos <- min(bc_plot$data[,2])
  dash_length <- abs(reduce(range(bc_plot$data[,2]), `-`)*0.1)
  middle <- which.min(abs(sort(df$t, decreasing=T)))
  
  bc_plot_custom <- bc_plot + geom_segment(aes(x=x, xend=x), y=lowest_pos,
                           yend=lowest_pos-dash_length) + 
    geom_line(colour='#4daf4a') + 
    geom_hline(yintercept=lowest_pos, colour='#cccccc') + 
    geom_hline(yintercept=0, colour='#cccccc') + xlab('') +
    theme_classic() +
    geom_tile(data=tibble(rank=1:length(stat_vector), 
                          y=lowest_pos-(1.25*dash_length)), 
              aes(x=rank, y=y, fill=rank),
                  width=1,
                  height=0.5*dash_length) +
    scale_fill_gradient2(low ='#b2182b', high='#2166ac', 
                         mid='#f7f7f7', midpoint = middle) + 
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    theme(panel.grid=element_blank(), 
          axis.text.x=element_blank(),
          axis.ticks.x = element_blank(),
          legend.position = 'none') + 
    ggtitle(paste(sig, 'signature')) +
    ylab('Enrichment score')
  
  return(bc_plot_custom)
}
```

barcode plots

```{r}

this_data <- calu_model$treat %>%
  select(stat=t,Symbol) %>%
  mutate(Symbol = gsub(pattern = "\\s",replacement = "",x = Symbol)) %>% 
  na.omit() %>% 
  distinct() %>% 
  group_by(Symbol) %>% 
  summarize(t=-mean(stat)) %>%
  ungroup()

down_fgsea <- fgsea_hallmark %>% arrange(NES)

## generate plots
bc_plots_down <- map(1:4, function(j){#nrow()
  bcp <- custom_barcode_plot(this_data, down_fgsea$pathway[j]) + 
    annotate('text', x=Inf , y=Inf, hjust=1, vjust=1, 
             label=paste('NES =', round(down_fgsea$NES[j], 2), 
                         '\nFDR =', round(down_fgsea$padj[j], 3)))
  
  return(bcp)
})

up_fgsea <- fgsea_hallmark %>% arrange(-NES)

## generate plots
bc_plots_up <- map(c(1,2,5,8), function(j){#nrow()
  bcp <- custom_barcode_plot(this_data, up_fgsea$pathway[j]) + 
    annotate('text', x=Inf , y=Inf, hjust=1, vjust=1, 
             label=paste('NES =', round(up_fgsea$NES[j], 2), 
                         '\nFDR =', round(up_fgsea$padj[j], 3)))
  
  return(bcp)
})

## plot to canvas
reduce(c(bc_plots_up), `+`) + plot_layout(ncol=2)

ggsave("plots/signature_enrichment_upregulated.pdf")

## plot to canvas
reduce(c(bc_plots_down), `+`) + plot_layout(ncol=2)

ggsave("plots/signature_enrichment_downregulated.pdf")


```

#write the top tables

```{r}

  write_delim(calu_model$treat,"results/calu_top_table_treatment.txt",delim = "\t")

write_delim(calu_model$time,"results/calu_top_table_time.txt",delim = "\t")

write_delim(a549_model$treat,"results/a549_top_table_treatment.txt",delim = "\t")

write_delim(calu_model$time,"results/a549_top_table_time.txt",delim = "\t")

```


#Cell line comparison

##Heatmaps

We see a lot of overlap in treatment responses of the two cell lines, thus we take the top 300 differential genes of Calu3 and plot their expression data in a heatmap for each cell line after filtering them for relevance in the top 7 and bottom 7 most enriched hallmark gene signatures.

```{r}

top_genes_symbol <- calu_model$treat %>% arrange(-abs(logFC)) %>% head(.,500) %>% pull(Symbol) %>% gsub(x = .,pattern = "\\s",replacement="")

relevant_genes <-unique(c(pathways.hallmark[[up_fgsea$pathway[1]]],
  pathways.hallmark[[up_fgsea$pathway[2]]],
  pathways.hallmark[[up_fgsea$pathway[3]]],
  pathways.hallmark[[up_fgsea$pathway[4]]],
  pathways.hallmark[[up_fgsea$pathway[5]]],
  pathways.hallmark[[up_fgsea$pathway[6]]],
  pathways.hallmark[[up_fgsea$pathway[7]]]),
  c(pathways.hallmark[[down_fgsea$pathway[1]]],
  pathways.hallmark[[down_fgsea$pathway[2]]],
  pathways.hallmark[[down_fgsea$pathway[3]]],
  pathways.hallmark[[down_fgsea$pathway[4]]],
  pathways.hallmark[[down_fgsea$pathway[5]]],
  pathways.hallmark[[down_fgsea$pathway[6]]],
  pathways.hallmark[[down_fgsea$pathway[7]]]))

focus_genes <- intersect(top_genes_symbol,relevant_genes)

focus_probes <- annotation_data %>% mutate(Symbol=gsub(x = Symbol,pattern = "\\s",replacement="")) %>% filter(Symbol %in% focus_genes) %>% group_by(Symbol) %>% do(head(.,1)) %>% distinct() %>% pull(ProbeID) %>% unique()

top_data_calu <- expression_data[focus_probes,]  %>% select(starts_with("Calu")) #%>% select(contains("24h")) 

top_data_calu %<>% rownames_to_column("ProbeID") %>% left_join(annotation_data)  %>% select(-Name,-ProbeID) %>% column_to_rownames("Symbol")

treatment <- factor(ifelse(grepl('mock', colnames(top_data_calu)), 'mock', 'infected'))

#we can either specify each time point as level 
time <- 
  ifelse(grepl('24h', colnames(top_data_calu)), 24,
         ifelse(grepl('8h', colnames(top_data_calu)), 8, 
                ifelse(grepl('12h', colnames(top_data_calu)), 
                12, 4)))

annotation_col = data.frame(
                    treatment = treatment, 
                    time = time
                )
rownames(annotation_col) = colnames(top_data_calu)

dev.off()
pdf("plots/Calu-3-heatmap-focused.pdf")
pheatmap(top_data_calu,main = "Calu -3 differential expressed genes",annotation_col = annotation_col,color = brewer.pal(9, "Reds"))#,annotation_col = annotation_col
dev.off()
```

We see a lot of overlap in treatment responses of the two cell lines, thus we take the top 300 differential genes of Calu3 and plot their expression data in a heatmap

```{r}

top_genes <- calu_model$treat %>% arrange(-abs(logFC)) %>% head(.,100) %>% pull(ProbeID)

top_data_calu <- expression_data[top_genes,]  %>% select(starts_with("Calu")) #%>% select(contains("24h")) 

top_data_calu %<>% rownames_to_column("ProbeID") %>% left_join(annotation_data)  %>% select(-Name,-ProbeID) %>% column_to_rownames("Symbol")

treatment <- factor(ifelse(grepl('mock', colnames(top_data_calu)), 'mock', 'infected'))

#we can either specify each time point as level 
time <- 
  ifelse(grepl('24h', colnames(top_data_calu)), 24,
         ifelse(grepl('8h', colnames(top_data_calu)), 8, 
                ifelse(grepl('12h', colnames(top_data_calu)), 
                12, 4)))

annotation_col = data.frame(
                    treatment = treatment, 
                    time = time
                )
rownames(annotation_col) = colnames(top_data_calu)

dev.off()
pdf("plots/Calu-3-heatmap-top_100_logFC.pdf",height = 16)
pheatmap(top_data_calu,main = "Calu -3 differential expressed genes",annotation_col = annotation_col,color = brewer.pal(9, "Reds"))#,annotation_col = annotation_col
dev.off()

```


Heatmap for A549

```{r}

top_data_a549<- expression_data[focus_probes,]  %>% select(starts_with("A549")) #%>% select(contains("24h")) 

top_data_a549 %<>% rownames_to_column("ProbeID") %>% left_join(annotation_data)  %>% select(-Name,-ProbeID) %>% column_to_rownames("Symbol")

treatment <- factor(ifelse(grepl('mock', colnames(top_data_a549)), 'mock', 'infected'))

#we can either specify each time point as level 
time <- 
  ifelse(grepl('24h', colnames(top_data_calu)), 24,
         ifelse(grepl('8h', colnames(top_data_calu)), 8, 
                ifelse(grepl('12h', colnames(top_data_a549)), 
                12, 4)))

annotation_col = data.frame(
                    treatment = treatment, 
                    time = time
                )
rownames(annotation_col) = colnames(top_data_a549)

dev.off()
pdf("plots/A549-heatmap-focused.pdf")
pheatmap(top_data_a549,main = "A549 differential expressed genes",annotation_col = annotation_col)#,annotation_col = annotation_col
dev.off()
```

##Correlation Scatterplot

```{r}


this_data <- calu_model$treat %>%
  select(ProbeID,FC_calu = logFC, Expr_calu = AveExpr,padj_calu = adj.P.Val) %>%
  mutate(FC_calu = -FC_calu) %>%
  left_join(a549_model$treat) %>%
  mutate(FC_a549 = -logFC,Expr_a549 = -AveExpr) %>%
  mutate(significant = if_else(padj_calu <= 0.1,"significant","non-significant"))


this_data %>%
  ggplot(aes(x = FC_calu,y = FC_a549,col = significant)) +
    geom_abline(slope = 1) +
    geom_point(data = subset(this_data,significant == "non-significant"),alpha=0.3) +
    geom_point(data = subset(this_data,significant == "significant")) +
    geom_vline(xintercept = 0,lty=2,col=b110_grey_light) +
    geom_label_repel(data = this_data %>% filter(significant == "significant") %>% arrange(-abs(FC_calu)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue))+ 
    ggtitle("Infection comparing Cell lines")

ggsave("plots/cell-line-comparison.pdf")


this_data %>%
  ggplot(aes(x = FC_calu,y = FC_a549,col = significant)) +
    geom_abline(slope = 1) +
    geom_point_rast(data = subset(this_data,significant == "non-significant"),alpha=0.3) +
    geom_point(data = subset(this_data,significant == "significant")) +
    geom_vline(xintercept = 0,lty=2,col=b110_grey_light) +
    geom_label_repel(data = this_data %>% filter(significant == "significant") %>% arrange(-abs(FC_calu)) %>% head(.,15),aes(label=Symbol),col=b110_grey) +
    theme_b110() +
    scale_color_manual(values = c(b110_grey_light,google_blue))+ 
    ggtitle("Infection comparing Cell lines")

ggsave("plots/cell-line-comparison_rast.pdf")

```

##Venn Diagrams

```{r}


calu_infection_upregulated <- calu_model$treat %>%
  filter(logFC < 0,adj.P.Val < 0.1) %>%
  pull(Symbol) %>% 
  gsub(pattern = "\\s",replacement = "",x = .)

calu_infection_downregulated <- calu_model$treat %>%
  filter(logFC > 0,adj.P.Val < 0.1) %>%
  pull(Symbol) %>% 
  gsub(pattern = "\\s",replacement = "",x = .)

a549_infection_upregulated <- a549_model$treat %>%
  filter(logFC < 0,adj.P.Val<0.2) %>%
  pull(Symbol) %>% 
  gsub(pattern = "\\s",replacement = "",x = .)

a549_infection_downregulated <- a549_model$treat %>%
  filter(logFC > 0,adj.P.Val < 0.2) %>%
  pull(Symbol) %>% 
  gsub(pattern = "\\s",replacement = "",x = .)

pdf("plots/Venn_diagrams_collection.pdf")

euler(list(calu_up=unique(calu_infection_upregulated),a549_up=unique(a549_infection_upregulated))) %>% plot(main="up-regulated genes")

euler(list(calu_down=unique(calu_infection_downregulated),a549_down=unique(a549_infection_downregulated))) %>% plot(main="down-regulated genes")

venn(list(calu_up=unique(calu_infection_upregulated),a549_up=unique(a549_infection_upregulated))) %>%
  plot(main="up-regulated genes")

venn(list(calu_down=unique(calu_infection_downregulated),a549_down=unique(a549_infection_downregulated))) %>% plot(main="down-regulated genes")

dev.off()

```

##Get consitently regulated genes

Up

```{r}
base::intersect(unique(calu_infection_upregulated),unique(a549_infection_upregulated))
```

Down

```{r}
base::intersect(unique(calu_infection_downregulated),unique(a549_infection_downregulated))
```

# Data integration

```{r}

a <- calu_model$treat %>%
  mutate(Symbol = gsub(pattern = "\\s",replacement = "",x = Symbol)) %>% 
  na.omit() %>% 
  distinct() %>% 
  group_by(Symbol) %>% 
  summarize(Calu3_microarray=-mean(logFC)) %>%
  ungroup() 

b <- a549_model$treat %>%
  mutate(Symbol = gsub(pattern = "\\s",replacement = "",x = Symbol)) %>% 
  na.omit() %>% 
  distinct() %>% 
  group_by(Symbol) %>% 
  summarize(A549_microarray=-mean(logFC)) %>%
  ungroup() 

a %>% full_join(b) %>%
  write_delim("results/foldchanges_microarray.txt",delim = "\t")

```

# Session info

```{r}
writeLines(capture.output(sessionInfo()), "results/SessionInfo.txt")
```

# Bibliography