Patients_GSE157256.Rmd

---
title: "GSE157256"
author: "Christina Schmidt"
date: "12 April 2021"
output:
  html_document:
    toc: true
    toc_depth: 3
    toc_float: true
editor_options: 
  chunk_output_type: console
---

```{r, message=FALSE, warning=FALSE}
library(tidyverse) # used for data manipulation
library(rmarkdown) # used for paged_table function
library(DESeq2)
```

# Introduction - Dataset
The dataset used here was produced by the Linhean lab and published as part of the [paper](https://pubmed.ncbi.nlm.nih.gov/33402335/) of Crooks. et. al. 2021. The RNA-seq data of HLRCC patients tumours/adjacent tissue (5 primary tumours and paired adjacent tissue, and 16 metastatic tumours) were downloaded from the GEO database, under accession [GSE157256](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE157256).\

```{r}
sampleinfo <-read.csv("InputData_Patients_GSE157256/SampleInfo_GSE157256.csv", header=TRUE)
```

# Quality control: Raw count data
## The data
For this the raw counts "GSE157256_RSEM_counts.txt" have been downloaded.\
Filtering the genes: I will keep all genes where the total number of reads across all samples is > than 5.
```{r, message=FALSE, warning=FALSE}
df_raw <- read.delim("InputData_Patients_GSE157256/GSE157256_RSEM_counts.txt", sep = "")%>%
  separate(col = 1, into = c("ENSEMBL", "Gene"), sep = "\\|")

#Filtering the genes: keep all genes where the rowSum >=  5
#Firstformat the table:
countdata <- df_raw[,-c(1:2)]%>%
               as.matrix()

# Filtering the data:
dim(countdata)
keep <- rowSums(countdata) >= 5
countdata <- countdata[keep,]
dim(countdata)
```

## Library sizes bar plot
First, we plot how many reads each sample has. Whilst normalisation can account for imbalance in coverage across the samples, extreme differences may be indicative of underlying problems in the samples.
```{r, message=FALSE, warning=FALSE}
librarySizes <- colSums(countdata)

barplot(librarySizes, 
        names=names(librarySizes), 
        las=2, 
        main="Barplot of library sizes")
abline(h=20e6, lty=2)
```

## Count distribution boxplots
Count data is not normally distributed, so if we want to examine the distributions of the raw counts it is helpful to transform the data on to a log scale. Typically we use a `log2` transformation, however, because the data is count data and will contain many `0`s we need to add a count of 1 to every 
value in order to prevent attempting `log2(0)` from creating errors.\
The boxplots show that overall the density distributions of raw log-intensities are a bit variable. If a sample is really far above or below the sample median (horizontal line) this could cause problems downstream, so we need to keep this in mind
```{r, message=FALSE, warning=FALSE}
# Get log2 counts per million
logcounts <- log2(countdata + 1)

#statusCol <- as.numeric(factor(sampleinfo$TissueType)) + 1
# Check distributions of samples using boxplots
boxplot(logcounts, 
        xlab="", 
        ylab="Log2(Counts)",
        las=2
        #col=statusCol
        )
# Let's add a blue horizontal line that corresponds to the median logCPM
abline(h=median(as.matrix(logcounts)), col="black")
```

## PCA
To run the PCA we should first normalise our data for library size and transform to a log scale using DESeq2 command `rlog` that performs a log2 scale transformation in a way that compensates for differences between samples for genes with low read  count and also normalizes between samples for library size.\
*[Information](http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#count-data-transformations)  about `rlog`, it's alternative `vst` and the comparison between the two.\
\
Since we are not interested in the metastatic samples, we will remove them from this analysis.\
Here we can see that the primary tumour and adjacent tissue samples cluster nicely apart on PC1.\
```{r ,message=FALSE, warning=FALSE}
library(ggfortify)

RawCounts <- read.delim("InputData_Patients_GSE157256/GSE157256_RSEM_counts.txt", sep = "")
RawCounts <- as.data.frame(t(RawCounts[,-1]))%>%
  rownames_to_column("Linehan_SampleID_R")
RawCounts <- merge(sampleinfo,RawCounts, by=c(1))

#Extract the Data for the downstream comparisons
#Remove metastatic samples
RawCounts_TvWT <- RawCounts%>%
  subset(TissueType!="Metastasis")
row.names(RawCounts_TvWT) <- RawCounts_TvWT$GSM_SampleID
sampleinfo_TvWT <-RawCounts_TvWT[,c(1:7)]
RawCounts_TvWT <- as.data.frame(t(RawCounts_TvWT[,-c(1:7)]))

# data normalization
rlogcounts <- t(rlog(as.matrix(RawCounts_TvWT)))
rlogcounts <- cbind(rownames(rlogcounts), data.frame(rlogcounts, row.names=NULL))
rlogcounts <- merge(x=sampleinfo,y=rlogcounts, by.x=3, by.y=1)
rownames(rlogcounts) <- rlogcounts$Linehan_SampleID_R

# run PCA
PCA<- autoplot (prcomp(rlogcounts[,-c(1:7)]),
         data= rlogcounts,
         colour = "TissueType", #colour = row including the sample information to colour code
         label=T,
         label.size=3,
         label.repel = TRUE,
         #loadings=T, #draws Eigenvectors
         #loadings.label = TRUE,
         #loadings.label.vjust = 1.2,
         #loadings.label.size=2,
         #loadings.colour="grey10",
         #loadings.label.colour="grey10",
         color = "black",
         fill = "TissueType",
         alpha = 0.3,
         shape = 21,
         size = 8
         ) +
    ggtitle("rlog normalised countdata", subtitle = "PCA color = Tissue" )
plot(PCA)
```

# Differential expression analysis
## Primary_Tumor versus Adjacent_Normal
For this we also use the raw counts "GSE157256_RSEM_counts.txt".
DESeq2:\
- Filtering the genes: keep all genes where the rowSum >=  5\
- We performed the DESeq analysis setting the significance level for independent filtering to 0.05.
```{r,message=FALSE, warning=FALSE}
#The Rawcounts data:
RawCounts <- read.delim("InputData_Patients_GSE157256/GSE157256_RSEM_counts.txt", sep = "")
RawCounts <- as.data.frame(t(RawCounts[,-1]))%>%
  rownames_to_column("Linehan_SampleID_R")
RawCounts <- merge(sampleinfo,RawCounts, by=c(1))

#Extract the Data for the downstream comparisons
#Tumour versus Adjacent 
RawCounts_TvWT <- RawCounts%>%
  subset(TissueType!="Metastasis")
row.names(RawCounts_TvWT) <- RawCounts_TvWT$GSM_SampleID
sampleinfo_TvWT <-RawCounts_TvWT[,c(1:7)]
RawCounts_TvWT <- as.data.frame(t(RawCounts_TvWT[,-c(1:7)]))

#Filtering the genes: keep all genes where the rowSum >=  5
#First we want to format the table:
countdata <- RawCounts_TvWT%>%
               as.matrix()

# Filtering the data:
dim(countdata)
keep <- rowSums(countdata) >= 5
countdata <- countdata[keep,]
dim(countdata)

# first lets check that our rows and columns match
all(sampleinfo_TvWT$GSM_SampleID == colnames(countdata))

# Safe the raw count matrix:
write.csv(countdata, "OutputData_Patients_GSE157256/Filtered_RawCountMatrix_GSE157256.csv")

# Create a DESeqDataSet object with the raw data:
# Convert the counts to DESeqDataSet Object. It is an object used by `DESeq2` to store count data -->It has a number of slots for storing count data, sample information, the model design for the differential expression analysis, and various other parameters about the data
ddsObj.raw <- DESeqDataSetFromMatrix(countData = countdata,
                                 colData = sampleinfo_TvWT,
                                 design = ~ TissueType)
# DESeq2:
# 1. Calculate the "median ratio" normalisation size factors:
ddsObj <- estimateSizeFactors(ddsObj.raw)
sizeFactors(ddsObj)

#2.Estimate dispersion (Gene wise dispersion)
ddsObj <- estimateDispersions(ddsObj)
#3.Apply Negative Binomial GLM fitting and calculate Wald statistics
ddsObj <- nbinomWaldTest(ddsObj)

# Generate a results table
#--> We can generate a table of differential expression results from the DDS object using the `results` function of DESeq2. res <- results(ddsObj, alpha=0.05) # default significance level for independent filtering is `0.1, we use 0.05!
res_TvWT <- results(ddsObj, alpha=0.05, contrast=c("TissueType","Primary_Tumor","Adjacent_Normal"))
resdata_TvWT<- as.data.frame(res_TvWT)

TvWT <- as.data.frame(res_TvWT)%>%
            rownames_to_column("ID")
TvWT <- separate(data = TvWT, col = ID, into = c("ENSEMBL", "Gene"), sep = "\\|")

#Remove duplicated genes: Function 
RemoveDublons <-function(MyData, OutputFileName){
  MyData <- MyData[complete.cases(MyData),]
  print("Number of genes in CS_DESeq_TvWT.csv dataset:")
  print(nrow(MyData))
  doublons <- as.character(MyData[duplicated(MyData$Gene),"Gene"])
  # Keep the entry with the greatest Log2FC:
  MyData$absLogFC <- abs(MyData$log2FoldChange)
  MyData <- MyData[ order(MyData$absLogFC), ]
  MyData_Select <- MyData[!duplicated(MyData$Gene),]
  print("Number of genes after removing duplicated gene names from dataset:")
  print(nrow(MyData_Select))
  #Safe:
  write.csv(MyData_Select, paste("",OutputFileName), row.names= TRUE)
  OutputFileName <-  MyData_Select
}

TvWT_ND <- RemoveDublons(MyData= TvWT,
              OutputFileName="DESeq_TvWT_Patients_GSE157256.csv")
```

# Gene set enrichment analysis (GSEA)
The analysis below are done using the DESeq2 result of the comparison primary tumour versus adjacent normal.\
\
For this we use gene set collections of biological pathways or molecular network information about a biological system, such as the Molecular signatures database (MsigDB) to find sets of altered signalling pathways or processes.\
\
GSEA was run using the t-value of the DESeq2 results for ranking.\
First I downloaded the signatures of interest from the [MsigDB](http://software.broadinstitute.org/gsea/msigdb) including "KEGG", "Reactome", "Biocarta", "Hallmarks" and "NFE2L2" (=NRF2) as well as "EMT" (= manually curated EMT-signature based on [Taube et al.](https://pubmed.ncbi.nlm.nih.gov/20713713/).).\
In order to be able to compare the results of the different signatures of interest, we combined the signatures and run them together. 

```{r, message=FALSE, warning=FALSE}
df_GSEA <- read.csv(" DESeq_TvWT_Patients_GSE157256.csv")

#Prepare data matrix for GSEA:
Column_Gene <- as.character(df_GSEA$Gene)
Column_tval <- as.numeric(df_GSEA$stat)
MyData_Extracted <- data.frame(cbind(Column_Gene, Column_tval), stringsAsFactors = F)
MyData_Extracted$Column_tval <- as.numeric(MyData_Extracted$Column_tval)

t_val <- as.numeric(MyData_Extracted$Column_tval)#Make the data into a vector
names(t_val) <- MyData_Extracted$Column_Gene

#Load the Pathways
library(fgsea)
library(GSEABase)
gmt_to_csv <- function(gmtfile, fast = T){
  if(fast)
  {
    genesets = GSEABase::getGmt(con = gmtfile)
    genesets = unlist(genesets)
    
    gene_to_term =plyr::ldply(genesets,function(geneset){
      temp <- geneIds(geneset)
      temp2 <- setName(geneset)
      temp3 <- as.data.frame(cbind(temp,rep(temp2,length(temp))))
      
    },.progress = plyr::progress_text())
    names(gene_to_term) <- c("gene","term")
    return(gene_to_term[complete.cases(gene_to_term),])
  }
  else
  {
    genesets = getGmt(con = gmtfile)
    genesets = unlist(genesets)
    
    gene_to_term <- data.frame(NA,NA)
    names(gene_to_term) <- c("gene","term")
    for (geneset in genesets)
    {
      temp <- geneIds(geneset)
      temp2 <- setName(geneset)
      temp3 <- as.data.frame(cbind(temp,rep(temp2,length(temp))))
      names(temp3) <- c("gene","term")
      gene_to_term <- rbind(gene_to_term,temp3)
    }
    
    return(gene_to_term[complete.cases(gene_to_term),])
  }
}

KEGG_Human <- gmt_to_csv("Input_MSigDB_Signatures/c2.cp.kegg.v6.2.symbols.gmt")
Reactome_Human <- gmt_to_csv("Input_MSigDB_Signatures/c2.cp.reactome.v6.2.symbols.gmt")
Biocarta_Human <- gmt_to_csv("Input_MSigDB_Signatures/c2.cp.biocarta.v6.2.symbols.gmt")
Hallmarks_Human <- gmt_to_csv("Input_MSigDB_Signatures/h.all.v6.2.symbols.gmt")
NRF2_Human <- gmt_to_csv("Input_MSigDB_Signatures/NFE2L2.gmt")
EMT_Human <-read.csv("Input_MSigDB_Signatures/EMT-signature.csv")

#Run the GSEA analysis
pathways <- rbind(KEGG_Human, Reactome_Human, Biocarta_Human, Hallmarks_Human,NRF2_Human,EMT_Human)
pathway_list <- list()
for(pathway in unique(pathways$term))
{pathway_list[[pathway]] <- as.character(pathways[pathways$term == pathway, 1])}

gsea_result1 <- fgsea(pathways = pathway_list, stats = t_val, nperm = 1000)
write_csv(gsea_result1[,c(-8)], "GSEA_result_KEGG-Hallmark-Reactome-Biocarta-NRF2-EMT_Patients_GSE157256.csv")

```

## GSEA Volcano Plots
We represented the results in a volcano plot. Important to note, for GSEA results a p.adj of < 25% (0.25) is considered significant.\
```{r,message=FALSE, warning=FALSE}
#Load the Data:
GSEA_Patients <-read.csv("OutputData_Patients_GSE157256/GSEA_result_KEGG-Hallmark-Reactome-Biocarta-NRF2-EMT_Patients_GSE157256.csv")
Patients <- read.csv("OutputData_Patients_GSE157256/DESeq_TvWT_Patients_GSE157256.csv")

#Establish the VolcanoPlot functions
library(ggrepel)
library(EnhancedVolcano)
VolcanoPlot_RNAseq_Patients <- function(Signature_Merge_RNA, Signature, OutputPlotName, NES, p.adj, Title){
  VolcanoPlot <-EnhancedVolcano (Signature_Merge_RNA,
                lab = Signature_Merge_RNA$symbol,#Metabolite name
                selectLab = c("CDC20", "CDC45", "MYC", "KPNA2", "CDKN2A", "CDKN1A"),
                x = "log2FoldChange",#Log2FC
                y = "padj",#p-value or q-value
                xlab = bquote(~Log[2]~ "FC"),
                ylab = bquote(~-Log[10]~p.adj),#(~-Log[10]~adjusted~italic(P))
                pCutoff = 0.05,
                FCcutoff = 0.5,#Cut off Log2FC, automatically 2
                pointSize = 4,
                labSize = 3,
                labFace = 'bold',
                boxedLabels = FALSE,
                titleLabSize = 12,
                subtitleLabSize = 8,
                col=c("darkseagreen3", "mediumaquamarine", "mediumaquamarine", "cyan4"),
                colAlpha = 0.5,
                title=(file=paste("RNAseq:", Title, " " ,sep="")),
                subtitle = (file=paste("", OutputPlotName, " (NES=", NES, "p.adj.=", p.adj, ")",sep="")),
                caption = paste0("total = ", nrow(Signature_Merge_RNA), " genes of ", nrow(Signature), " genes in pathway"),
                xlim = c((((Reduce(min,Signature_Merge_RNA$log2FoldChange))))-0.3,(((Reduce(max,Signature_Merge_RNA$log2FoldChange))))+0.3),
                ylim = c(0,((-log10(Reduce(min,Signature_Merge_RNA$padj))))+0.2),
                drawConnectors = TRUE,
                widthConnectors = 0.5,
                colConnectors = "black",
                cutoffLineType = "dashed",
                cutoffLineCol = "black",
                cutoffLineWidth = 0.5,
                #legendLabels=c('No changes',"-0.5< Log2FC <0.5","-0.5< Log2FC <0.5", 'p.adj<0.05 & -0.5< Log2FC <0.5"'),
                legendPosition = 'right',
                legendLabSize = -1,
                legendIconSize = -1
                )
  ggsave(file=paste("Figures_Patients_GSE157256/VolcanoPlot_GSEA", OutputPlotName,Title, "SelectedLabels.pdf", sep="_"), plot=VolcanoPlot, width=10, height=8)
  #ggsave(file=paste("VolcanoPlot_GSEA", OutputPlotName,Title, "SelectedLabels.png", sep="_"), plot=VolcanoPlot, width=10, height=8)
  plot(VolcanoPlot)
}

VolcanoPlot_GSEA_Patients <- function(GSEA_Results, OutputfileName){
  Volcano1 <- separate(GSEA_Results,"pathway", into = c("signature", "rest"), sep= "_", remove=FALSE)%>% 
    mutate(colour = case_when(signature =="KEGG" ~ 'blue',
                              signature =="BIOCARTA" ~ 'gold4',
                              signature =="HALLMARK" ~ 'deeppink4',
                              signature =="REACTOME" ~ 'seagreen4',
                              signature =="NFE2L2.V2" ~ 'pink',
                              signature =="EMTdown" ~ 'red',
                              signature =="EMTup" ~ 'red',
                                  TRUE ~ 'Not_Detected'))
  keyvals <- ifelse(
    Volcano1$colour == "blue", "blue",
    ifelse(Volcano1$colour == "gold4", "gold4",
    ifelse(Volcano1$colour == "deeppink4", "deeppink4",
    ifelse(Volcano1$colour == "seagreen4", "seagreen4",
    ifelse(Volcano1$colour == "pink", "pink",
    ifelse(Volcano1$colour == "red", "red",
           "black"))))))
  keyvals[is.na(keyvals)] <- 'black'
  names(keyvals)[keyvals == 'blue'] <- "KEGG"
  names(keyvals)[keyvals == 'gold4'] <- "BIOCARTA"
  names(keyvals)[keyvals == 'deeppink4'] <- "HALLMARK"
  names(keyvals)[keyvals == 'seagreen4'] <- "REACTOME"
  names(keyvals)[keyvals == 'pink'] <- "NRF2"
  names(keyvals)[keyvals == 'red'] <- "EMT"
    names(keyvals)[keyvals == 'black'] <- 'X'
  VolcanoPlot <- EnhancedVolcano (Volcano1,
                lab = Volcano1$pathway,#Metabolite name
                selectLab = c("HALLMARK_E2F_TARGETS", "HALLMARK_G2M_CHECKPOINT", "REACTOME_DNA_REPLICATION", "REACTOME_CELL_CYCLE_CHECKPOINTS", "HALLMARK_MYC_TARGETS_V1"),
                x = "NES",#Log2FC
                y = "padj",#p-value or q-value
                xlab = "NES",
                ylab = bquote(~-Log[10]~p.adj),#(~-Log[10]~adjusted~italic(P))
                pCutoff = 0.25,
                FCcutoff = 0.5,#Cut off Log2FC, automatically 2
                pointSize = 3,
                labSize = 1,
                labFace = 'bold',
                boxedLabels = FALSE,
                colCustom = keyvals,
                titleLabSize = 16,
                col=c("black", "grey", "grey", "purple"),#if you want to change colors
                colAlpha = 0.5,
                title=paste(OutputfileName),
                subtitle = bquote(italic("GSEA")),
                caption = paste0("total = ", (nrow(Volcano1)), " Pathways"),
                #xlim = c(-3.5,2.5),
                #ylim = c(0,1.5),
                ylim = c(0,((-log10(Reduce(min,GSEA_Results$padj))))+0.05),
                drawConnectors = TRUE,
                widthConnectors = 0.5,
                colConnectors = "black",
                cutoffLineType = "dashed",
                cutoffLineCol = "black",
                cutoffLineWidth = 0.5,
                legendLabels=c('No changes',"-0.5< NES <0.5", 'p.adj<0.25 & -0.5< NES <0.5"'),
                legendPosition = 'right',
                legendLabSize = 12,
                legendIconSize = 5.0
                )
  ggsave(file=paste("Figures_Patients_GSE157256/VolcanoPlot_GSEA_KEGG-Hallmark-Reactome-Biocarta-EMT-NRF2", OutputfileName, "SelectedLabels.pdf", sep="_"), plot=VolcanoPlot, width=10, height=8)
  #ggsave(file=paste("Figures_Patients_GSE157256/VolcanoPlot_GSEA_KEGG-Hallmark-Reactome-Biocarta-EMT-NRF2", OutputfileName, "SelectedLabels.png", sep="_"), plot=VolcanoPlot, width=10, height=8)
  plot(VolcanoPlot)
}

```
\
Make an overview plot of the GSEA results:
```{r, message=FALSE, warning=FALSE}
VolcanoPlot_GSEA_Patients(GSEA_Results= GSEA_Patients,
                 OutputfileName="Patients_TvN")
```
\
- Signatures to show: HALLMARK_E2F_TARGETS; HALLMARK_G2M_CHECKPOINT; REACTOME_DNA_REPLICATION; REACTOME_CELL_CYCLE_CHECKPOINTS; HALLMARK_MYC_TARGETS_V1.\
- From those signatures, highlight the following genes: Cdc20, Cdc45, Myc, Kpna2, Cdkn2a, Cdkn1a.
```{r, message=FALSE, warning=FALSE}
#-----------------------------------------------------------------------
#HALLMARK_E2F_TARGETS
Signature <- subset(Hallmarks_Human, term == "HALLMARK_E2F_TARGETS")
Signature_Merge_RNA <- merge(x=Signature,y=Patients, by.x="gene", by.y="Gene", all.x=TRUE)%>%
  na.omit()%>%
  rename("gene"="symbol")

VolcanoPlot_RNAseq_Patients(Signature_Merge_RNA=Signature_Merge_RNA, 
                   Signature=Signature, 
                   OutputPlotName="HALLMARK_E2F_TARGETS", 
                   NES="3.43", 
                   p.adj="0.02",
                   Title="HLRCC Patients (Tumour vs. Normal")

#-----------------------------------------------------------------------
#HALLMARK_MYC_TARGETS_V1
Signature <- subset(Hallmarks_Human, term == "HALLMARK_MYC_TARGETS_V1")
Signature_Merge_RNA <- merge(x=Signature,y=Patients, by.x="gene", by.y="Gene", all.x=TRUE)%>%
  na.omit()%>%
  rename("gene"="symbol")

VolcanoPlot_RNAseq_Patients(Signature_Merge_RNA=Signature_Merge_RNA, 
                   Signature=Signature, 
                   OutputPlotName="HALLMARK_MYC_TARGETS_V1", 
                   NES="2.76", 
                   p.adj="0.02",
                   Title="HLRCC Patients (Tumour vs. Normal")

```

# Information about packge used and versions 
```{r}
sessionInfo()
```