MetaStudy.Rmd

---
title: "MetaStudy"
author: "Moritz Madern"
date: '2023-01-14'
output:
  pdf_document: default
  html_document: default
---


This script compiles differential expression data from files of different experiments (each provided in a standardized format & containing gene-wise FC vs p-values) and performs meta-analysis on specified genes. 

```{r read in packages, echo=FALSE, warning=FALSE, message=FALSE}

## read in packages
library(ggplot2)
library(tidyverse)
library(RColorBrewer)
library(rlist)
library(poolr)
library(superheat)
library(gplots)

```


<br><br><br>

## Specify parameters

```{r specify parameters}

## specify file path to datasets to read in (all datasets should be located in the same folder)
filepath_datasets = "./Datasets"


## specify file path to table that defines the conditions within those datasets of interest
filepath_tableConditions = "Condition_datasets_11012023.csv"


## specify genes to be analyzed (note that they will be concverted to capital letters to make pattern searching case-invariant)
vector_geneNames = c("Cdc45", "eIF2A","Rpi7","Gapdh","Ncor1","Stat1", "Gata3", "Stat4")


## specify subset of those genes (from vector_geneNames) to be illustrated in a Heatmap. 
vector_geneNames_heatmap = c("Cdc45","Ncor1","Stat1", "Stat4", "Gata3")


## specify heatmap color scheme. Suppurted colorschemes are "viridis" and "redblue"
heatmap_colorscheme = "redblue"


```



<br><br><br>

## Read in datasets

Read in conditions table:

```{r read in tableConditions and specify some variables, echo=FALSE}

## read in the tabe, then print it
tableConditions = read.delim(filepath_tableConditions, header=TRUE, sep=",")
tableConditions$dataset <- as.character(tableConditions$dataset)
print(tableConditions)


## add combined column (unique identifier)
tableConditions$uniqueID = paste0("Dataset", tableConditions$dataset, "_", tableConditions$condition.number)


## get unique dataset number vector
vector_uniqueDatasetNumbers = unique(tableConditions$dataset) %>% as.character()


## give unique color to each dataset
colpal = colorRampPalette( c(brewer.pal(n=8, "Set1"), brewer.pal(n=8, "Set2")))
DatasetCols = setNames(colpal(length(vector_uniqueDatasetNumbers)), nm=vector_uniqueDatasetNumbers)

```

Read in all datasets as specified by filepath_datasets:

```{r read in each dataset and store in a list, warning=FALSE, echo=FALSE, message=FALSE}

## read in all datasets
files = list.files(filepath_datasets, full.names = TRUE)
m = files %>% length()
list_datasets = vector(mode="list")
for (i in files){
  print(i)
  list_datasets[[i]] <- read.table(file = i, header = TRUE, sep = "\t")
}         


## extract dataset numbers from filenames. Rename list based on dataset numbers
DatasetNumbers = sub(files, pattern=" .*", replacement = "") %>% sub(., pattern="[.]/Datasets/Dataset", replacement = "") 
names(list_datasets) = DatasetNumbers


## extract dataset descriptions, and name them by DatasetNumbers. Print dataframe that shows 1:1 relation
DatasetDescription = sub(x=files, pattern=".* ", replacement = "") %>% sub(., pattern=".csv", replacement = "")
names(DatasetDescription) = DatasetNumbers


```


```{r extract relevant data for plotting, warning=FALSE, echo=FALSE, message=FALSE}

## initiate list
list_of_dfs = list()

## go over each pair of dataset and condition in tableConditions 
for (i in 1:nrow(tableConditions)){
  tableConditions_i <- tableConditions[i,]
  dataset_i <- tableConditions_i$dataset
  condition_i <- tableConditions_i$condition.number
  dataset_condition_i <- tableConditions_i$uniqueID
  
  # extract dataset i
  df_dataset <- list_datasets[[dataset_i]]
  
  # select relevant column indices (if multiple match, take first one!)
  ind_condition.foldChange <- grepl(names(df_dataset), pattern=paste0(condition_i,"[.]foldChange")) %>% which()
  ind_condition.foldChange <- ind_condition.foldChange[1]
  ind_condition.pValue <- grepl(names(df_dataset), pattern=paste0(condition_i,"[.]pValue")) %>% which()
  ind_condition.pValue <- ind_condition.pValue[1]
  
  # extract relevant columns and store in dataframe
  df_i <- data.frame(DatasetNumber_Condition = dataset_condition_i,
                     GeneName =  df_dataset$Gene.Name,
                     foldChange = df_dataset[,ind_condition.foldChange],
                     pValue = df_dataset[,ind_condition.pValue])
  list_of_dfs[[i]] <- df_i
}


## merge to a single dataframe
df_all <- do.call(what=rbind, args=list_of_dfs)


## Set gene names to capital letters, then complete missing information (for unique combinations of GeneName and DatasetNumber_Condition)
df_all$GeneName <- toupper(df_all$GeneName)
df_all_complete <- complete(data=df_all, DatasetNumber_Condition, GeneName)

```


<br><br><br>

## Plot data for each specified gene

Plots show foldChange vs experiment, points are colored by statistical significance. Note that experiments in which the respective gene was not quantified are not shown. Furthermore, for each gene, the individual p-values were combined to an overall p-value via Fisher's method (https://en.wikipedia.org/wiki/Fisher%27s_method), rounded to 6 digits. This method assumes that the individual p-values to be pooled are independent but test for the same Null hypothesis. If these assumptions are not met, ignore the combined p-value.

```{r plot fc vs p-value, echo = FALSE, message = FALSE, warning = FALSE, fig.align="center", fig.path='figures/', dev=c('png', 'pdf')}

## initiate a list to be filled with foldChanges (including missing values) to be rearranged in wide-format. Could be used to plot heatmaps. Currently not implemented in this script
list_wide <- vector(mode="list", length = length(vector_geneNames))
names(list_wide) = toupper(vector_geneNames)


## initiate plot list
plot_list = list()


## go over each gene and plot fc vs p-value. Make sure the search is case-invariant by converting all letters to capital
vector_geneNames <- toupper(vector_geneNames)
for (g in vector_geneNames){
  
  # select data for gene g
  df_g <- df_all_complete %>% as.data.frame() %>% filter(GeneName == g)
  rownames(df_g) <- df_g$DatasetNumber_Condition
  
  # catch exception where gene is not found in any dataset
  if (nrow(df_g) == 0){
    paste0("Gene ", g, "could not be found in any dataset")
    next()
  }
  
  # reorder, then save fold changes in list. Next drop missing values
  df_g <- df_g[tableConditions$uniqueID,]
  list_wide[[g]] <- matrix(df_g$foldChange, ncol=1, dimnames = list(df_g$DatasetNumber_Condition, g))
  bool_keep <- !is.na(df_g$foldChange)
  df_g <- df_g[bool_keep,]
  
  # convert DatasetNumber_Condition to factor
  conditions_all = tableConditions$uniqueID
  factorlevels_condition_g <- conditions_all[conditions_all %in% df_g$DatasetNumber_Condition]
  df_g$DatasetNumber_Condition <- factor(df_g$DatasetNumber_Condition, levels = rev(factorlevels_condition_g))
  
  ## create significance vector
  df_g$significance <- cut(df_g$pValue, breaks=c(0, 0.001, 0.01, 0.05, 1.01), right=FALSE)
  levels(df_g$significance) <- c("p-value < 0.001", "p-value < 0.01", "p-value < 0.05", "non-significant")
  
  # calculate combined p-value via Fisher's method (assuming independent p-values that test the same hypothesis)
  p_combined_res <- poolr::fisher(p=df_g$pValue)
  p_combined <- round(p_combined_res$p, digits=6)
  
  # plot data for gene g
  xlim <- df_g$foldChange %>% abs() %>% max()
  gg_singleGene <- ggplot(df_g) +
                      geom_point(aes(x=foldChange, y=DatasetNumber_Condition, col=significance, pval=pValue), size=6) +
                      scale_color_manual(values=c(rgb(red=255, green=50, blue=50, alpha=255, maxColorValue = 255),
                                                  rgb(red=245, green=130, blue=0, alpha=255, maxColorValue = 255),
                                                  rgb(red=250, green=200, blue=0, alpha=255, maxColorValue = 255),
                                                  rgb(red=200, green=200, blue=200, alpha=255, maxColorValue = 255)),
                                         drop=FALSE)+
                      xlim(-xlim,xlim) +
                      theme_bw() + ylab("Experiment") + xlab("Fold change [log2]") +
                      geom_vline(xintercept=0, linetype=2, color="black") +
                      ggtitle(label = g, subtitle = paste("combined p-value:", p_combined ))

                      
  
  
  gg_singleGene
  plot_list <- list.append(plot_list, gg_singleGene)
}


```


```{r plot the gene-wise information, echo = FALSE, message = FALSE, warning = FALSE, fig.align="center", fig.path='figures/', dev=c('png', 'pdf')}

## plot the plots outside a loop
plot_list

```

Note: Upon compiling the script to pdf, each plot is saved as separate pdf and png file in the newly generated folder called "figures" in your working directory.



<br><br><br>

## Plot Heatmap


For selected genes, a heatmap is plotted that depicts fold changes. Missing values (i.e. NAs) are colored as grey.

```{r write Heatmap function, echo = FALSE, message = FALSE, warning = FALSE}


### Heatmap plot function ###
heatmap_plot <- function(m, groups, legend_colors, samplenames, dendrogram="column", labrow="", bool_rowv=TRUE, bool_colv = TRUE, plot_path=NULL, title=""){
  
  # load packages
  library(gplots)
  
  # create groups
  names(colors) <- levels(groups)
  
  # replaces NAs with 0
  m[is.na(m)] <- 0
  colnames(m) <- samplenames
  
  # should rows be reordered
  if (bool_rowv){
    rowv <- as.dendrogram(hclust(dist(m)))
  } else {
    rowv <- FALSE
  }
  
  # should columns be reordered
  if (is.logical(bool_colv)){
    if (bool_colv){
      colv <- as.dendrogram(hclust(dist(t(m))))
    } else {
      colv <- FALSE
    }
  } else {
    colv <- bool_colv
  }
  
  # specify colors
  if(is.null(legend_colors)){
    sidecolors <- rep("white", times=ncol(m))
  } else{
    sidecolors <- legend_colors[groups]
  }
  
  # create color palette
  heatmap_pal <- colorRampPalette(rev(brewer.pal(11, "RdBu")))
  
  # plot heatmap for centered log2 intensities if centered
  min_m <- min(m, na.rm=TRUE)
  max_m <- max(m, na.rm=TRUE)
  heatmap.2(m,         
            Rowv = rowv,
            Colv=colv,
            labRow=labrow, margins=c(8,8), ColSideColors = sidecolors, trace="none",col=heatmap_pal(50),
            breaks = seq(from=-2,to=2, length.out=51), 
            symkey = F,
            dendrogram=dendrogram, main=title)
  
  # save plot if specified 
  if(!is.null(plot_path)){
    pdf(file=plot_path, width = 6, height = 6)
    heatmap.2(m,         
              Rowv = rowv,
              Colv=colv,
              labRow=labrow, margins=c(8,8), ColSideColors = sidecolors, trace="none",col=heatmap_pal(50),
              breaks = seq(from=-2,to=2, length.out=51), 
              symkey = F,
              dendrogram=dendrogram,
              main=title)
    dev.off()
  }
}

```


```{r plot heatmap, fig.align="center", echo = FALSE, message = FALSE, warning = FALSE, fig.path='figures/', dev=c('png', 'pdf'),fig.width=6}

## from list_wide, select those genes to be visualized in heatmap as specified in the parameter section
list_wide_heatmap <- list_wide[toupper(vector_geneNames_heatmap)]


## condense to matrix
m_heatmap <- do.call(what =cbind, args = list_wide_heatmap)


## remove rows with only NAs
bool_kick <- apply(is.na(m_heatmap), FUN=all, MARGIN = 1)
m_heatmap <- m_heatmap[!bool_kick,]


## prepare heatmap colors
heatmap_pal_RdBU<- colorRampPalette(rev(brewer.pal(11, "RdBu")))
if (heatmap_colorscheme == "viridis"){
  heatmap_pal <- colorRampPalette(viridis::viridis(21))
}
if (heatmap_colorscheme == "redblue"){
  heatmap_pal <- colorRampPalette(rev(brewer.pal(11, "RdBu")))
}


## plot heatmap
heatmap.2(m_heatmap, margins=c(8,8), trace="none",col=heatmap_pal(30), na.color = "grey",
            breaks = seq(from=-2,to=2, length.out=31),
            symkey = F, scale="none",cexCol = 1, cexRow = 0.8,
            dendrogram="none", Rowv=NA, Colv=NA, xlab="Gene Name", density.info = "none")



```

```{r}

```