DGE_Epithelial.Rmd

---
title: 'DGE Analysis Notebook - Epithelial_response: HighInducerVsLowInducer'
output:
  html_document: 
    toc: yes
    fig_width: 8
    fig_height: 8
    fig_caption: yes
    number_sections: yes
    toc_depth: 4
    df_print: tibble
---

```{r, message = FALSE, warning = FALSE}
suppressPackageStartupMessages(library("dplyr"))
suppressPackageStartupMessages(library("DESeq2"))
suppressPackageStartupMessages(library("pheatmap"))
suppressPackageStartupMessages(library("PoiClaClu"))
suppressPackageStartupMessages(library("RColorBrewer"))
suppressPackageStartupMessages(library('tidyverse'))
suppressPackageStartupMessages(library("PoiClaClu"))
suppressPackageStartupMessages(library("vsn"))
suppressPackageStartupMessages(library('EnhancedVolcano'))
suppressPackageStartupMessages(library('gplots'))
suppressPackageStartupMessages(library('org.Mm.eg.db'))
suppressPackageStartupMessages(library('stringr'))
suppressPackageStartupMessages(library("genefilter"))
suppressPackageStartupMessages(library("dplyr"))
suppressPackageStartupMessages(library("ggplot2"))
suppressPackageStartupMessages(library("glmpca"))
suppressPackageStartupMessages(library('org.Mm.eg.db'))
suppressPackageStartupMessages(library("AnnotationDbi"))
suppressPackageStartupMessages(library("apeglm"))
suppressPackageStartupMessages(library("ComplexHeatmap"))
suppressPackageStartupMessages(library("clusterProfiler"))
suppressPackageStartupMessages(library('ggrepel'))
suppressPackageStartupMessages(library('corrplot'))
suppressPackageStartupMessages(library("GO.db"))
suppressPackageStartupMessages(library('edgeR'))
suppressPackageStartupMessages(library('GOstats'))
suppressPackageStartupMessages(library('pathview'))
suppressPackageStartupMessages(library("gage"))
suppressPackageStartupMessages(library("gageData"))
suppressPackageStartupMessages(library('GOSemSim'))
suppressPackageStartupMessages(library('DOSE'))
suppressPackageStartupMessages(library('enrichplot'))
suppressPackageStartupMessages(library('ggnewscale'))
suppressPackageStartupMessages(library('glue'))
suppressPackageStartupMessages(library("ggupset"))
#suppressPackageStartupMessages(library("SPIA"))
suppressPackageStartupMessages(library("stats"))
suppressPackageStartupMessages(library("FactoMineR"))
suppressPackageStartupMessages(library("factoextra"))
suppressPackageStartupMessages(library("pcaExplorer"))
```

# Differential Gene Expression Analysis

## Creating metadata for the DGE Analysis

DESeq2 needs sample information (metadata) for performing DGE analysis.
Let's create the sample information

```{r}
# Read the csv file and change the column name. the samples.csv is a list of sample names, ie, the names of bam files.
sample_ID <- read.csv("/home/keshavprasadgubbi/Documents/AlinaRnaSeq/samples.csv")
head(sample_ID)
```

```{r}
condition <- c("Infected", "Infected", "Infected", "Infected", "Infected", "Infected",
               "Infected", "Infected", "Infected", "Infected", "Infected", "Infected",
               "Infected", "Infected", "control", "control")
coldata <- data.frame(sample_ID, condition)
colnames(coldata) <- c('Sample_Name','condition') # change name of one of the columns
# The metadata can be found in a df called coldata!
head(coldata)
```

### Tidying up the names for plots later!

#### First from coldata

```{r}
# tidying up the names od samples in both columns that list of samples
#coldata$Samples <- str_remove_all(coldata$Sampls, pattern = "run6_trimmed_|_.bam|_S\\d\\d|_S\\d")
coldata$Sample_Name <- str_remove_all(coldata$Sample_Name,
                                  pattern = "run6_trimmed_|_.bam|_S\\d\\d|_S\\d" )
coldata$condition <- as.factor(coldata$condition)
```

Changing the names of samples (as per Alina)

```{r}
coldata[coldata == '476_R1'] <- 'T'
coldata[coldata == '754_R1'] <- 'S54'
coldata[coldata == '755_R1'] <- 'S55'
coldata[coldata == '757_R1'] <- 'L57'
coldata[coldata == '758_R1'] <- 'A58'
coldata[coldata == '760_R1'] <- 'L60'
coldata[coldata == '761_R1'] <- 'S61'
coldata[coldata == '762_R1'] <- 'A62'
coldata[coldata == '763_R1'] <- 'L63'
coldata[coldata == '764_R1'] <- 'A64'
coldata[coldata == '765_R1'] <- 'S65'
coldata[coldata == '766_R1'] <- 'L66'
coldata[coldata == '768_R1'] <- 'A68'
coldata[coldata == '769_R1'] <- 'L69'
coldata[coldata == 'Ctrl1_R1'] <- 'C1'
coldata[coldata == 'Ctrl2_R2'] <- 'C2'
# convert column1 with sample names to row.names of coldata
rownames(coldata) <- coldata$Sample_Name
coldata
```

### Adding the groupings by Alina for further Metadata Information

```{r}
coldata$Epithelial_response <- c("LowInducer", "LowInducer", "HighInducer",
                                  "HighInducer", "LowInducer", "LowInducer",
                                  "HighInducer", "HighInducer", "LowInducer",
                                  "HighInducer", "LowInducer", "HighInducer",
                                  "LowInducer", "LowInducer", 'NR', 'NR')
coldata$clinical_outcome <- c('symptomatic', 'symptomatic', 'symptomatic',
                              'Lethal', 'asymptomatic', 'Lethal', 'symptomatic',
                              'asymptomatic', 'Lethal', 'symptomatic', 'symptomatic',
                              'Lethal', 'asymptomatic', 'Lethal', 'NR', 'NR')
coldata$microcolonies <- c('Low', 'Low', 'Low', 'High', 'Low', 'Low',
                           'High', 'High', 'Low', 'High', 'Low', 'High', 'Low',
                           'Low', 'NR', 'NR')
coldata$ER_microcolonies <- c("LI_LM", "LI_LM", "HI_LM", "HI_HM", "LI_LM", "LI_LM",
                              "HI_HM", "HI_HM", "LI_LM", "HI_HM", "LI_LM", "HI_HM",
                              "LI_LM", "LI_LM", 'NR', 'NR')
coldata$phylogenomic_lineage <- c("EPEC1", "EPEC10", "EPEC9", "EPEC9", "NC", "EPEC5",
                                  "EPEC8", "NC", "EPEC7", "NC", "EPEC2", "EPEC9",
                                  "EPEC2", "EPEC2", 'NR', 'NR')
coldata$phylogroup <- c("B2", "A", "B2", "B2", "B1", "A", "B2", "B2", "B1", "B2", "B1",
                        "B2", "B1", "B2", 'NR', 'NR')
coldata$Intimin_Type <- c("alpha", "ND", "lambda", "lambda", "epsilon", "epsilon",
                          "mu", "lambda", "beta", "kappa", "beta", "alpha", "beta",
                          "beta", 'NR', 'NR')
```

#### then fix Countsmatrix:

NOTE:

1.  From the manuals the countsData must be a numeric matrix
2.  It is IMPORTANT to keep the names of the genes in the rownames

```{r}
# Readin  countsmatrix
#countsmatrix <-as.matrix(read.csv("~/R/Rtuts/Data/Alina_EPEC_project/counts.csv"))
countsmatrix <- read.csv("/home/keshavprasadgubbi/Documents/AlinaRnaSeq/newcounts.csv")
#countsmatrix <- as.data.frame(countsmatrix)
```

```{r}
## Removal of Gender Genes from ENSEMBL ID itself
countsmatrix <- countsmatrix %>% filter(countsmatrix$X != "ENSMUSG00000086503",
                                  countsmatrix$X != "ENSMUSG00000097571",
                                  countsmatrix$X != "ENSMUSG00000086370",
                                  countsmatrix$X != "ENSMUSG00000031329")
nrow(countsmatrix)
#countsmatrix <- as.matrix(countsmatrix)
```

```{r}
#tidying up these names again
#colnames(countsmatrix) <- str_remove_all(colnames(countsmatrix), pattern = "run6_trimmed_|_.bam|_S\\d\\d|_S\\d")
rownames(countsmatrix) <- countsmatrix[,1] #converting first column of gene names into rownames, to be used for sanity check later
# It is IMPORTANT to keep the names of the genes in the rownames
countsmatrix <- subset(countsmatrix, select = - X)#dropping the X column
# the elements from Sample_Name from coldata must the the colnames of countsmatrix
colnames(countsmatrix) <- coldata$Sample_Name
# Display the column names
colnames(countsmatrix)
```

## Annotating and Exporting ENSEMBL ID into Gene Symbols

Adding genes annotated from ENSEMBL ID to Gene symbols and ENTREZ Id to countsmatrix table. Will be keeping the symbols and entrez columsn to be added later into results table as it is for later use
```{r}
cm_row <- rownames(countsmatrix)
head(cm_row)
# Mapping the ENSEMBL ID to Symbol and ENTREZ ID
symbols <- mapIds(
  org.Mm.eg.db,
  keys = cm_row,
  column = c('SYMBOL'),
  keytype = 'ENSEMBL',
  multiVals = "first"
)

```
```{r}
symbols <- symbols[!is.na(symbols)]
symbols <- symbols[match(rownames(countsmatrix), names(symbols))]
head(symbols, 25)

# Creating a new column called genename and putting in the symbols and entrez columns into count matrix
countsmatrix$genename <- symbols

# Removing all rows with NA values for genenames, so that those rows are filtered out.
countsmatrix <- unique(countsmatrix[rowSums(is.na(countsmatrix)) == 0, ]) # Apply rowSums & is.na
nrow(countsmatrix)
# Moving the ENSEMBL ID from rownames into separate column for itself.
countsmatrix <- tibble::rownames_to_column(countsmatrix, "E_ID")
# Removing the duplicated genes so that then these genes can be made into rownames for countsmatrix
countsmatrix <- distinct(countsmatrix[!duplicated(countsmatrix$genename), ])
```

```{r}
# Now make the ganename column into rownames of count matrix
rownames(countsmatrix) <- countsmatrix[,"genename"]
# Keeping this version of countsmatrix for later use
cm_table <- countsmatrix
# dropping the column E_ID, genenames so that only numeric values are present in it as an input of DESEq Object.
countsmatrix <- subset(countsmatrix, select = -c(genename, E_ID))#

# Changing countsmatrix into Matrix of numeric values so that only numeric values are present in it as an input of DESEq Object.
countsmatrix <- as.matrix(countsmatrix)
class(countsmatrix) <- "numeric"
```


# Calculating CPM Values

```{r}
# as DGEList
dge_er <- DGEList(counts = countsmatrix)

dim(dge_er)
colnames(dge_er)
#dge_er$samples

## calculate norm. factors
nr <- calcNormFactors(dge_er)

## get normalized counts
cpmvalues <- cpm(nr)
cpmvalues_d <- cpm.default(nr)
```

# Differential Gene Expression analysis using DESeq2

Now, construct DESeqDataSet for DGE analysis.

But before that, a sanity check : It is essential to have the name of
the columns in the count matrix in the same order as that in name of the
samples (rownames in coldata).

```{r}
all(rownames(coldata) %in% colnames(countsmatrix))
ncol(countsmatrix) == nrow(coldata)
dim(countsmatrix)
```

## Creating the DESeq Data set Object

```{r}
dds_epithelial <- DESeqDataSetFromMatrix(countData = countsmatrix,
                              colData = coldata, 
                              design = ~ Epithelial_response)
nrow(dds_epithelial)
```

```{r}
# Function to save generic plots
saveplot <- function(plot,name ){
  # Function to save the plots
  ggsave(filename = 
           glue('/home/keshavprasadgubbi/Documents/AlinaRnaSeq/ER/{name}.png'),
       plot = plot,
       dpi = 300,
       width = 10,
       height = 10,
       units = "in")
}
```

## Exploratory Data Analysis and Visualization

### Pre-filtering the dataset

Our count matrix with our DESeqDataSet contains many rows with only
zeros, and additionally many rows with only a few fragments total. In
order to reduce the size of the object, and to increase the speed of our
functions, we can remove the rows that have no or nearly no information
about the amount of gene expression.

Applying the most minimal filtering rule: removing rows of the
DESeqDataSet that have no counts, or only a single count across all
samples. Additional weighting/filtering to improve power is applied at a
later step in the workflow.

```{r}
keep <- rowSums(counts(dds_epithelial)) > 1
dds_epithelial <- dds_epithelial[keep,]
nrow(dds_epithelial)
```

### The variance stabilizing transformation

## Applying VST transformation

```{r}
vsd <- vst(dds_epithelial, blind = FALSE)
#head(assay(vsd), 3)
colData(vsd)
vsd_coldata <- colData(vsd)
```

```{r}
dds_epithelial <- estimateSizeFactors(dds_epithelial)
dds_epithelial
```

## Sample Distances

useful first step in an RNA-seq analysis is often to assess overall
similarity between samples:

1.  Which samples are similar to each other, which are different?
2.  Does this fit to the expectation from the experiment's design?

### Euclidean Distance between samples

dist to calculate the Euclidean distance between samples - useful for
ONLY normalized data. To ensure we have a roughly equal contribution
from all genes, we use it on the VST data.

```{r}
sampleDists <- dist(t(assay(vsd)))
sampleDistMatrix <- as.matrix( sampleDists )
rownames(sampleDistMatrix) <- vsd$Sample_Name
colnames(sampleDistMatrix) <- vsd$Sample_Name

colors <- colorRampPalette( rev(brewer.pal(9, "RdYlBu")) )(255)
(EuclideanDistanceHeatmap <- pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         main = "Sample-to-Sample Euclidean Distance of BL6 - LowInducer_vs_HighInducer",
         col = colors))

```

```{r}
#saveplot(EuclideanDistanceHeatmap, "EuclideanDistanceHeatmap")
```

### Poisson Distance between Samples

```{r}
poisd <- PoissonDistance(t(counts(dds_epithelial))) # raw counts or unnormalised data
samplePoisDistMatrix <- as.matrix( poisd$dd )
rownames(samplePoisDistMatrix) <- dds_epithelial$Sample_Name
colnames(samplePoisDistMatrix) <- dds_epithelial$Sample_Name

colors <- colorRampPalette( rev(brewer.pal(9, "RdYlBu")) )(255)
(poisson_dist_plot <- pheatmap(samplePoisDistMatrix,
         clustering_distance_rows = poisd$dd,
         clustering_distance_cols = poisd$dd,
         main = "Sample-to-Sample Poisson Distance of BL6 - LowInducer_vs_HighInducer",
         col = colors))
         
```

```{r}
#saveplot(poisson_dist_plot, "poisson_dist_plot")
```

# PCA Plot
```{r}
### Functions for Plot aethetics and saving PCA Plots
color_values <- c("red", "red","red","red", "black","black","red","red", "red",
                  "red", "red","red","red", "red", "red" ,"blue")

# the basic set of common aesthetic settings for PCA plots, 
theme.my.own <- list(theme_bw() ,
                      geom_point(size = 3),
                      coord_fixed() ,
                      scale_y_continuous(breaks = seq(-20, 20, 5), 
                                         sec.axis = sec_axis(~. *1, 
                                                             labels = NULL,
                                                             breaks = NULL)) ,
                      scale_x_continuous(breaks = seq(-20, 20, 5), 
                                         sec.axis = sec_axis(~. *1, 
                                                             labels = NULL,
                                                             breaks = NULL)) ,
                      theme_classic() ,
                      geom_hline(yintercept = 0, color = "gray", size = 1) ,
                      geom_vline(xintercept = 0, color = "gray", size = 1) ,
                      theme(text = element_text(size = 15),
                            axis.text = element_text(size = 15),
                            legend.position = "right",
                            aspect.ratio = 1) ,
                      #geom_text(size = 4, hjust = 0, vjust = 0)
                      geom_text_repel(size = 5,min.segment.length = 0.5)
  )

```

## Calculating all PCA Values

```{r}
plotPCA_local = function(object,
                         intgroup = "condition",
                         ntop = 500,
                         returnData = TRUE,
                         nPC = 4)
{
  # calculate the variance for each gene
  rv <- rowVars(assay(object))
  ntop <- 500
  
  # select the ntop genes by variance
  select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
  
  # perform a PCA on the data in assay(x) for the selected genes
  pca <- prcomp(t(assay(object)[select, ]))
  #summary(pca)
  
  # the contribution to the total variance for each component
  percentVar <- pca$sdev ^ 2 / sum(pca$sdev ^ 2)
  
  if (!all(intgroup %in% names(colData(object)))) {
    stop("the argument 'intgroup' should specify columns of colData(dds)")
  }
  
  intgroup.df <-
    as.data.frame(colData(object)[, intgroup, drop = FALSE])
  
  # add the intgroup factors together to create a new grouping factor
  group <- if (length(intgroup) > 1) {
    factor(apply(intgroup.df, 1, paste, collapse = ":"))
  } else {
    colData(object)[[intgroup]]
  }
  
  # assembly the data for the plot
  d <- cbind(pca$x[, seq_len(min(nPC, ncol(pca$x))), drop = FALSE],
             data.frame(group = group, intgroup.df, name = colnames(object)))
  
  if (returnData) {
    attr(d, "percentVar") <- percentVar[1:nPC]
    # l <- list(pca,d)
    # return(l)
     return(d)
  }
}
  
```

## PCA Plot with VST Data

### Function for calculating percentvar

```{r}
percentvar_calculation <- function(pcaData_variable){
  # function to calculate percentvar for different variables
  percentvar_variable <- round(100 * attr(pcaData_variable, "percentVar"), digits = 3 )
  return(percentvar_variable)
}

savingFunction <- function(plotname, metadatacolumn){
  # Function to save the PCA plots
  ggsave(filename = 
           glue('/home/keshavprasadgubbi/Documents/AlinaRnaSeq/ER/PCAplot_HighVsLowInducer_{metadatacolumn}.png'),
       plot = plotname,
       dpi = 300,
       width = 10,
       height = 10,
       units = "in")
}
```

```{r}
pcaData_epithelial <- plotPCA_local(vsd, intgroup = c("Epithelial_response","Sample_Name"),returnData = T)
pcaData_epithelial
percentVar_epithelial <- percentvar_calculation(pcaData_epithelial)
```

```{r }
percentVar_epithelial
```
```{r }
(PCAplot_vst <- ggplot(pcaData_epithelial,
                      aes(x = PC1,
                          y = PC2,
                          color = Sample_Name,
                          label = Sample_Name)) +
        xlab(paste0("PC1: ", percentVar_epithelial[1], "% variance")) +
        ylab(paste0("PC2: ", percentVar_epithelial[2], "% variance")) +
        ggtitle("PCA Plot - High Vs Low Inducer") +
        scale_colour_manual(values = color_values) +
        theme.my.own )

savingFunction(PCAplot_vst, "Epithelial_response")
```
## PCA Plot for PC2 vs PC3

```{r}

(PCAplot_vst23 <- ggplot(pcaData_epithelial,
                      aes(x = PC2,
                          y = PC3,
                          color = Sample_Name,
                          label = Sample_Name)) +
        xlab(paste0("PC2: ", percentVar_epithelial[2], "% variance")) +
        ylab(paste0("PC3: ", percentVar_epithelial[3], "% variance")) +
        ggtitle("PCA Plot - High Vs Low Inducer") +
        scale_colour_manual(values = color_values) +
        theme.my.own )

ggsave(filename = '/home/keshavprasadgubbi/Documents/AlinaRnaSeq/ER/PCAplot23_HighVsLowInducer_epithelial.png',
       plot = PCAplot_vst23,
       dpi = 300,
       width = 10,
       height = 10,
       units = "in")

```


## PCA Plot for PC3 vs PC4

```{r}

(PCAplot_vst34 <- ggplot(pcaData_epithelial,
                      aes(x = PC3,
                          y = PC4,
                          color = Sample_Name,
                          label = Sample_Name)) +
        xlab(paste0("PC3: ", percentVar_epithelial[3], "% variance")) +
        ylab(paste0("PC4: ", percentVar_epithelial[4], "% variance")) +
        ggtitle("PCA Plot - High Vs Low Inducer") +
        scale_colour_manual(values = color_values) +
        theme.my.own )

ggsave(filename = '/home/keshavprasadgubbi/Documents/AlinaRnaSeq/ER/PCAplot34_HighVsLowInducer_epithelial.png',
       plot = PCAplot_vst34,
       dpi = 300,
       width = 10,
       height = 10,
       units = "in")

```


## PCA Plot for different groupings of metadata

### PCA for Microcolonies

```{r}
PCAdata_epithelial_MC <- plotPCA(vsd, intgroup = c("Sample_Name", "microcolonies"), returnData = TRUE)
percentVar_epithelial_MC <- percentvar_calculation(PCAdata_epithelial_MC)

(PCAplot_MC <- ggplot(PCAdata_epithelial_MC,
                     aes(x = PC1, y = PC2,
                         color = microcolonies,
                         label= Sample_Name)) +
        xlab(paste0("PC1: ", percentVar_epithelial_MC[1], "% variance")) +
        ylab(paste0("PC2: ", percentVar_epithelial_MC[2], "% variance")) +
        ggtitle("PCA Plot - High Vs Low Inducer - Microcolonies") +
        theme.my.own)

savingFunction(PCAplot_MC, "microcolonies")
```

### PCA for Clinical Outcome

```{r}
PCAdata_epithelial_CO <- plotPCA(vsd, 
                                 intgroup = c("Sample_Name", "clinical_outcome"), 
                                 returnData = TRUE)
percentVar_epithelial_CO <- percentvar_calculation(PCAdata_epithelial_CO)

(PCAplot_CO <- ggplot(PCAdata_epithelial_CO,
                     aes(x = PC1, y = PC2,
                         color = clinical_outcome,
                         label = Sample_Name)) +
        xlab(paste0("PC1: ", percentVar_epithelial_CO[1], "% variance")) +
        ylab(paste0("PC2: ", percentVar_epithelial_CO[2], "% variance")) +
        ggtitle("PCA Plot - High Vs Low Inducer - Clinical Outcome") +
        theme.my.own)

savingFunction(PCAplot_CO, "clinical_outcome")
```

### PCA for ER_Microcolonies

```{r}
PCAdata_epithelial_ERMC <- plotPCA(vsd, intgroup = c("Sample_Name", "ER_microcolonies"), returnData = TRUE)
percentVar_epithelial_ERMC <- percentvar_calculation(PCAdata_epithelial_ERMC)

(PCAplot_ERMC <- ggplot(PCAdata_epithelial_ERMC,
                       aes(x = PC1, y = PC2,
                           color = ER_microcolonies,
                           label = Sample_Name)) +
        xlab(paste0("PC1: ", percentVar_epithelial_ERMC[1], "% variance")) +
        ylab(paste0("PC2: ", percentVar_epithelial_ERMC[2], "% variance")) +
        ggtitle("PCA Plot - High Vs Low Inducer - ER_Microcolonies") +
        theme.my.own)

savingFunction(PCAplot_ERMC, "ER_microcolonies")
```

### PCA for Phylogenomic Lineage

```{r}
PCAdata_epithelial_PL <- plotPCA(vsd, intgroup = c("Sample_Name", 
                                                   "phylogenomic_lineage"), 
                                 returnData = TRUE)
percentVar_epithelial_PL <- percentvar_calculation(PCAdata_epithelial_PL)

(PCAplot_PL <- ggplot(PCAdata_epithelial_PL,
                       aes(x = PC1, y = PC2,
                           color = phylogenomic_lineage,
                           label = Sample_Name)) +
        xlab(paste0("PC1: ", percentVar_epithelial_PL[1], "% variance")) +
        ylab(paste0("PC2: ", percentVar_epithelial_PL[2], "% variance")) +
        ggtitle("PCA Plot - High Vs Low Inducer - Phylogenomic Lineage") +
        theme.my.own)

savingFunction(PCAplot_PL, "phylogenomic_lineage")
```

### PCA for Phylogroup

```{r}
PCAdata_epithelial_PG <- plotPCA(vsd, intgroup = c("Sample_Name", "phylogroup"), returnData = TRUE)
percentVar_epithelial_PG <- percentvar_calculation(PCAdata_epithelial_PG)

(PCAplot_PG <- ggplot(PCAdata_epithelial_PG,
                     aes(x = PC1, y = PC2,
                         color = phylogroup,
                         label= Sample_Name)) +
        xlab(paste0("PC1: ", percentVar_epithelial_PG[1], "% variance")) +
        ylab(paste0("PC2: ", percentVar_epithelial_PG[2], "% variance")) +
        ggtitle("PCA Plot - High Vs Low Inducer - Phylogroup") +
        theme.my.own)

savingFunction(PCAplot_PG, "phylogroup")
```

### PCA for Intimin Group

```{r}
PCAdata_epithelial_IT <- plotPCA(vsd, intgroup = c("Sample_Name", "Intimin_Type"), returnData = TRUE)
head(PCAdata_epithelial_IT)
percentVar_epithelial_IT <- round(100 * attr(PCAdata_epithelial_IT, "percentVar"), digits = 4)

(PCAplot_IT <- ggplot(PCAdata_epithelial_IT,
                     aes(x = PC1, y = PC2,
                         color = Intimin_Type,
                         label = Sample_Name)) +
        xlab(paste0("PC1: ", percentVar_epithelial_IT[1], "% variance")) +
        ylab(paste0("PC2: ", percentVar_epithelial_IT[2], "% variance")) +
        ggtitle("PCA Plot - High Vs Low Inducer - Intimin type") +
        theme.my.own)

savingFunction(PCAplot_IT, "Intimin_Type")
```

```{r}
# calculate the variance for top 500 gene
rv <- rowVars(assay(vsd))
ntop <- 500
# select the ntop genes by variance
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
df1 <- t(assay(vsd)[select,])
```
```{r}
res.pca <- PCA(df1,  graph = FALSE, scale.unit = FALSE)
summary.PCA(res.pca)
```
```{r}
# Visualize eigenvalues/variances
fviz_screeplot(res.pca, addlabels = TRUE)
```
```{r}
library("factoextra")
eig.val <- get_eigenvalue(res.pca)
eig.val
```
```{r}
var <- get_pca_var(res.pca)
#fviz_pca_var(res.pca, repel = TRUE)
```

```{r fig.width=10}
library("corrplot")
corrplot(var$cos2, is.corr=T)
```

## Genes + PCA Biplots

```{r}
fviz_pca_biplot(res.pca, repel = TRUE,
                gradient.cols = c("pink", "blue", "yellow", "green", "red", "black"))
```

```{r fig.height=10, fig.width=10}
heat.colors <- brewer.pal(6, "RdYlBu")
fviz_pca_var(res.pca, col.var = "contrib", repel = TRUE,
             gradient.cols = c("Gray", "blue", "yellow","orange", "green", "red", "black"),
             )
```

```{r}
# Contributions of variables to PC2
fviz_contrib(res.pca, choice = "var", axes = 2, top = 25)
```
```{r}
# Contributions of variables to PC1
fviz_contrib(res.pca, choice = "var", axes = 1, top = 25)
```

## Hierarchical Clustering

### applying rlog Transformation

```{r}
rld <- rlog(dds_epithelial, blind = FALSE)
head(assay(rld), 3)
```

```{r}
### Extract the rlog matrix from the object
rld_mat <- assay(rld) #assay() is function from the "SummarizedExperiment" package that was loaded when you loaded DESeq2

### Compute pairwise correlation values
rld_cor <- cor(rld_mat)    ## cor() is a base R function
head(rld_cor)   ## check the output of cor(), make note of the rownames and colnames
```

```{r}
### Plot heatmap
heat.colors <- brewer.pal(6, "RdYlBu")
(Hclust_plot <- pheatmap(rld_cor, 
                        color = heat.colors,
                        main = 'Heirarchical Clustering of Samples - ER - Correlation Matrix'
                        #filename = '/home/keshavprasadgubbi/Documents/AlinaRnaSeq/ER/Hclust_plot.tiff'
                       ))
#Hclust_plot
```

```{r}
#saveplot(Hclust_plot, "Hclust_plot")
```

## DGE Results

### Running the differential expression pipeline

```{r}
dds1_epithelial <- DESeq(dds_epithelial)
#str(dds1)
```

### Building the results table

```{r}
res_epithelial <- results(dds1_epithelial,
                          contrast = c("Epithelial_response",
                                        "HighInducer", 
                                        "LowInducer"))
head(res_epithelial, 30)
```

```{r}
summary(res_epithelial)
```

A Problem: there a number of rows with padj that has NA.

Reaons:

1.  If within a row, all samples have zero counts, the baseMean column
    will be zero, and the log2 fold change estimates, p value and
    adjusted p value will all be set to NA. ------\>\>\>\>\>\> checked
    and not true!

2.  If a row contains a sample with an extreme count outlier then the p
    value and adjusted p value will be set to NA. These outlier counts
    are detected by Cook's distance. Customization of this outlier
    filtering and description of functionality for replacement of
    outlier counts and refitting is described below. ------\>\>\>\>\>\>
    checked and not true!

3.  If a row is filtered by automatic independent filtering, for having
    a low mean normalized count, then only the adjusted p value will be
    set to NA. Description and customization of independent filtering is
    described below. ------\>\>\>\>\>\> could be true!

=====\>\>\>\>\>\> Most likely caused by automatic independent filtering
due to the presence of low mean normalized counts.

Solution: Obtain unfiltered DESeq2 results! (counts without automatic
independent filtering and outlier removal). Then, only the genes with
ALL counts set to zero will have NA for pvalues.

```{r}
dds2_epithelial <- DESeq(dds_epithelial, minReplicatesForReplace = Inf)
res2_epithelial <- results(
  dds2_epithelial,
  cooksCutoff = FALSE,
  independentFiltering = FALSE,
  contrast = c("Epithelial_response",
               "HighInducer",
               "LowInducer")
)
head(res2_epithelial, 30)
```

### Results

```{r}
res2df_epithelial <- as.data.frame(res2_epithelial) # convert the results table to a df
head(res2df_epithelial)
```

## MA Plot

```{r}
plotMA_res2_epithelial <- plotMA(res2_epithelial, ylim = c(-2, 2))
```

### Histogram of p-values

```{r}
hist(res2_epithelial$pvalue, breaks = 50, col = "grey50", border = "blue")
```

Further Filtering: baseMean \> 1

```{r}
hist(res2_epithelial$pvalue[res2_epithelial$baseMean > 10],
     breaks = 50, col = "grey50", border = "blue")
```
```{r}
head(res2df_epithelial)
```
Now I have only Gene names and no longer have ENSEMBLID. So do the following to obtain Gene symbols and ENTREZ ID.
1. make  new column of gene names from rownames and keep the rownames as well as it is.
2. map the gene symbols into Entrez Id and make it into a separate column in itself.

```{r}
res2df_epithelial <- tibble::rownames_to_column(res2df_epithelial, "symbol")

gn <- res2df_epithelial$symbol

# Mapping the Symbol to ENTREZ ID

entrez <- mapIds(
  org.Mm.eg.db,
  keys = gn,
  column = "ENTREZID",
  keytype = "SYMBOL",
  multiVals = "first"
)

ensembl_id <- mapIds(
  org.Mm.eg.db,
  keys = gn,
  column = "ENSEMBL",
  keytype = "SYMBOL",
  multiVals = "first"
)
res2df_epithelial$entrez <- entrez
res2df_epithelial$ensemblID <- ensembl_id
```

Omit NA values from symbol and respective rows!

```{r}
res3df_epithelial <- res2df_epithelial %>% filter(!is.na(symbol) & !is.na(entrez))
nrow(res3df_epithelial)
```

## Saving the Results

```{r}
resOrdered_epithelial <- res3df_epithelial[order(res3df_epithelial$pvalue),]
head(resOrdered_epithelial)

write.csv(as.data.frame(resOrdered_epithelial),
          file = "/home/keshavprasadgubbi/Documents/AlinaRnaSeq/ER/results_DGE_ER.csv")
```

## Heatmap of count matrix

To explore a count matrix, it is often instructive to look at it as a heatmap.

```{r}
select <- order(rowMeans(counts(dds2_epithelial,normalized = FALSE)),decreasing = TRUE)[1:20]
df <- as.data.frame(colData(dds2_epithelial)[,c("Epithelial_response","Sample_Name")])
Heatmap(assay(vsd)[select,], 
         cluster_columns = TRUE
        )
```

## Effect of Transformations on Variance

-   These set of plots depict the standard deviation of transformed data
    (across samples), against mean.

### Based on Shifted Log Transformation

```{r}
ntd <- normTransform(dds2_epithelial)
meanSdPlot(assay(ntd))
```

### Based on Rlog Transformation

```{r}
meanSdPlot(assay(rld))
```

### Based on Variance Stabilizing Transformation

```{r}
meanSdPlot(assay(vsd))
```

## Dispersion Plots

-   Its a useful diagnostic to plot the dispersion estimates.

```{r}
plotDispEsts(dds2_epithelial)
```

# Volcano Plots

## Volcano Plots based on Enhanced Volcano

```{r fig.height= 10, fig.width= 10}
p1 <- EnhancedVolcano(res3df_epithelial,
    lab = res3df_epithelial$symbol,
    x = 'log2FoldChange',
    y = 'pvalue',
    xlab = bquote(~Log[2]~'FoldChange'),
    pCutoff = 0.05,
    FCcutoff = 1.0,
    title = 'DE genes: Log2FoldChange Vs -Log10 Pvalue',
    subtitle = bquote(~Log[2]~'|FoldChange| = 1, pvalue < 0.05') ,
    pointSize = 2.0,
    labSize = 5.0,
    boxedLabels = FALSE,
    gridlines.major = FALSE,
    gridlines.minor = FALSE,
    colAlpha = 0.5,
    legendPosition = 'bottom',
    legendLabSize = 12,
    legendIconSize = 4.0,
    drawConnectors = T,
    widthConnectors = 0.75,
    max.overlaps = 4
    )
(volcano1 <- p1 + scale_y_continuous(limits = c(0, 8),
                                    breaks = seq(0, 8, 1),
                                    sec.axis = sec_axis(~. *1,
                                                        labels = NULL,
                                                        breaks = NULL)) +
 
                  scale_x_continuous(limits = c(-4, 5), 
                                     breaks = seq(-4, 5, 1),
                                    sec.axis = sec_axis(~. *1,
                                                        labels = NULL,
                                                        breaks = NULL)) +
                      xlab(expression(DownRegulated %<->% UpRegulated))
  )

ggsave(filename = '/home/keshavprasadgubbi/Documents/AlinaRnaSeq/ER/Volcano1_L2fcVsPvalue.png',
       plot = volcano1,
       dpi = 300,
       width = 5,
       height = 5,
       units = "in")
```
## volcano Plot of Log2FC vs Padj

```{r fig.height= 10, fig.width= 10}
p2 <- EnhancedVolcano(res3df_epithelial,
    lab = res3df_epithelial$symbol,
    x = 'log2FoldChange',
    y = 'padj',
    xlab = bquote(~Log[2]~'FoldChange'),
    ylab = bquote(~Log[10]~'Padj'),
    pCutoff = 0.05,
    FCcutoff = 1.0,
    title = 'DE genes: Log2FoldChange Vs -Log10 Padj',
    subtitle = bquote(~Log[2]~'|FoldChange| = 1, pvalue < 0.05') ,
    pointSize = 2.0,
    labSize = 5.0,
    boxedLabels = FALSE,
    gridlines.major = FALSE,
    gridlines.minor = FALSE,
    colAlpha = 0.5,
    legendPosition = 'bottom',
    legendLabSize = 12,
    legendIconSize = 4.0,
    drawConnectors = T,
    widthConnectors = 0.75,
    max.overlaps = 4
    )

(volcano2 <- p2 + scale_y_continuous(limits = c(0, 8),
                                    breaks = seq(0, 8, 1),
                                    sec.axis = sec_axis(~. *1,
                                                        labels = NULL,
                                                        breaks = NULL)) +
                  scale_x_continuous(limits = c(-4, 5), 
                                     breaks = seq(-4, 5, 1),
                                    sec.axis = sec_axis(~. *1,
                                                        labels = NULL,
                                                        breaks = NULL))+
                      xlab(expression(DownRegulated %<->% UpRegulated))
  )

ggsave(filename = '/home/keshavprasadgubbi/Documents/AlinaRnaSeq/ER/Volcano2_L2fcVsPadj.png',
       plot = volcano2,
       dpi = 300,
       width = 10,
       height = 10,
       units = "in")
```

## Significant Differentially Expressed Genes

```{r}
res3df_epithelial$diffexpressed <- "NS"

# if log2Foldchange > 1.0 and pvalue < 0.05, set as "UP"
res3df_epithelial$diffexpressed[res3df_epithelial$log2FoldChange > 1.0 & res3df_epithelial$pvalue < 0.05] <- "UP"

# if log2Foldchange < -1.0 and pvalue < 0.05, set as "DOWN"
res3df_epithelial$diffexpressed[res3df_epithelial$log2FoldChange < -1.0 & res3df_epithelial$pvalue < 0.05] <- "DOWN"

# Create a new column "delabel" to de, that will contain the name of genes differentially expressed (NA in case they are not)
res3df_epithelial$delabel <- NA
res3df_epithelial$delabel[res3df_epithelial$diffexpressed != "NS"] <- res3df_epithelial$symbol[res3df_epithelial$diffexpressed != "NS"]
```

### Arrive at relevant genes by imposing thresholds.

```{r}
sigsdf_epithelial <- res3df_epithelial[(abs(res3df_epithelial$log2FoldChange) > 1) & (res3df_epithelial$pvalue < 0.05),]

nrow(sigsdf_epithelial)
```
### Gender Genes Removed from Table

```{r}
sigsdf_epithelial <- filter(sigsdf_epithelial,
                          symbol != "Xist",
                          symbol != "Jpx", 
                          symbol != "Ftx", 
                          symbol != "Tsx",
                          symbol != "Cnbp2" )

nrow(sigsdf_epithelial)
#head(sigsdf_epithelial)
```

Therefore, gender genes werent so much in play in terms of significance!

```{r}

write.csv(sigsdf_epithelial ,"/home/keshavprasadgubbi/Documents/AlinaRnaSeq/ER/SignificantGenes_ER.csv")
```

### Number of Genes from different strains that are contributing to UP/DOWN regulation.

```{r}
sigsdf_epithelial_UP <- sigsdf_epithelial[(sigsdf_epithelial$log2FoldChange) > 1, ] #UP Regulation Table
sigsdf_epithelial_DOWN <- sigsdf_epithelial[(sigsdf_epithelial$log2FoldChange) < -1, ] #DOWN Regulation Table

nrow(sigsdf_epithelial_UP)
nrow(sigsdf_epithelial_DOWN)
```

### Determine Top20 UP Genes and Top20 DOWN Genes

```{r}
UpGene_epithelial <- sigsdf_epithelial_UP[order(-sigsdf_epithelial_UP$log2FoldChange,decreasing = TRUE), ]$symbol
DownGene_epithelial <- sigsdf_epithelial_DOWN[order(-sigsdf_epithelial_DOWN$log2FoldChange,decreasing = TRUE), ]$symbol
DE_Genes_table_epithelial <- as.data.frame(cbind(UpGene_epithelial,DownGene_epithelial)) # sorted based on highest Log2FC value
```

```{r}
write.csv(DE_Genes_table_epithelial ,"/home/keshavprasadgubbi/Documents/AlinaRnaSeq/ER/DEGenesTable_ER.csv")
head(DE_Genes_table_epithelial, 25)
```

```{r}
head(sigsdf_epithelial)
```

```{r}
# Determining the significant Genes based on Log2FC and pvalue thresholds
sigs2df_epithelial <- res2df_epithelial[(abs(res2df_epithelial$log2FoldChange) > 1) &  (res2df_epithelial$pvalue < 0.05),]
```

## Z-score based Gene Heatmaps

### with Whole table (all genes together!)

```{r}
#mat <- counts(dds2_epithelial, normalized = TRUE)[rownames(sigsdf_epithelial),]
mat <- counts(dds2_epithelial, normalized = TRUE)[(sigsdf_epithelial$symbol),]
mat.zs <- t(apply(mat, 1, scale)) # Calculating the zscore for each row
colnames(mat.zs) <- coldata$Sample_Name # need to provide correct sample names for each of the columns
head(mat.zs)
```

```{r}
mat2 <- counts(dds2_epithelial, normalized = TRUE)[(sigsdf_epithelial$symbol),]
mat2.zs <- t(apply(mat2, 1, scale)) # Calculating the zscore for each row
colnames(mat2.zs) <- coldata$Sample_Name # need to provide correct sample names for each of the columns
head(mat2.zs)
#pheatmap(mat2.zs, cluster_cols = TRUE, cluster_rows = FALSE, show_rownames = FALSE)
```
```{r}
newHP <- Heatmap(mat2.zs, 
        cluster_columns = TRUE, 
        cluster_rows = TRUE, 
        column_labels = colnames(mat2.zs),
        name = 'DE Genes Heatmap',
        show_row_names = FALSE, 
        use_raster = TRUE, 
        raster_quality = 5,
        #width = unit(8, "in"), 
        #height = unit(8, "in")
        row_labels = sigs2df_epithelial[rownames(mat2.zs),]$symbol
        )
newHP
```
```{r fig.height= 100, fig.width= 10}
LongHeatMap_Allgenes <- Heatmap(mat2.zs, 
        cluster_columns = TRUE, 
        cluster_rows = TRUE, 
        column_labels = colnames(mat2.zs),
        row_labels = sigs2df_epithelial$symbol,
        name = 'All DE Genes',
        show_row_names = TRUE, 
        use_raster = TRUE, 
        raster_quality = 5,
        )
LongHeatMap_Allgenes
#saveplot(LongHeatMap_Allgenes, "LongHeatMap_Allgenes")
```

Need to filter these results to accommodate better the heat maps and also volcano plots

### with tighter constraints (all genes together!)

```{r}
sigs1df_epithelial <- res2df_epithelial[(res2df_epithelial$baseMean > 500) & (abs(res2df_epithelial$log2FoldChange) > 1) & (res2df_epithelial$pvalue < 0.05),]
mat1 <- counts(dds2_epithelial, normalized = TRUE)[(sigs1df_epithelial$symbol),]
mat1.zs <- t(apply(mat1, 1, scale)) # Calculating the zscore for each row
colnames(mat1.zs) <- coldata$Sample_Name # need to provide correct sample names for each of the columns
head(mat1.zs)
```

```{r fig.height= 10, fig.width= 10}
Heatmap(mat1.zs,
        cluster_columns = TRUE,
        cluster_rows = TRUE,
        column_labels = colnames(mat1.zs),
        name = 'DE Genes',
        row_labels = sigs1df_epithelial$symbol)
```

# GO Terms with clusterProfiler

## GO Terms for UP Regulated Genes

### GO over-representation analysis for UP Regulated Genes
```{r}
UPgene_ENS_ID <- (sigsdf_epithelial_UP$ensemblID)

GO_UPRegResults_epithelial <- enrichGO(gene = UPgene_ENS_ID,
                                       OrgDb = "org.Mm.eg.db",
                                       keyType = "ENSEMBL",
                                       ont = "BP",
                                       pAdjustMethod = "BH",
                                       pvalueCutoff  = 0.01,
                                       qvalueCutoff  = 0.05,
                                       readable      = TRUE)
```

```{r}
GO_UpRegdf_epithelial <- as.data.frame(GO_UPRegResults_epithelial)
head(GO_UpRegdf_epithelial)
```

```{r}
write.csv(GO_UpRegdf_epithelial ,"/home/keshavprasadgubbi/Documents/AlinaRnaSeq/ER/GO_UpReg_results_Epithelial.csv")
```

```{r fig.height=10, fig.width=10}
GO_UPReg_Barplot_epithelial <- plot(barplot(GO_UPRegResults_epithelial, showCategory = 22))
saveplot(GO_UPReg_Barplot_epithelial, "GO_UPReg_Barplot_epithelial")
```

```{r fig.height=10, fig.width=10}
GO_UPReg_Dotplot_epithelial <- plot(dotplot(GO_UPRegResults_epithelial, showCategory = 23))

saveplot(GO_UPReg_Dotplot_epithelial, "GO_UPReg_Dotplot_epithelial")
```

```{r fig.height=10, fig.width=12}
GO_UPReg_Cnetplot_epithelial <- plot(cnetplot(GO_UPRegResults_epithelial, showCategory = 13))

saveplot(GO_UPReg_Cnetplot_epithelial, "GO_UPReg_Cnetplot_epithelial")
```

### Upset Plot

```{r fig.height= 8, fig.width=12}
GO_UPReg_Upsetplot <- plot( upsetplot(GO_UPRegResults_epithelial))
saveplot(GO_UPReg_Upsetplot, "GO_UPReg_Upsetplot")
```

### Heatplot

The heatplot is similar to cnetplot, while displaying the relationships
as a heatmap. The gene-concept network may become too complicated if
user want to show a large number significant terms. The heatplot can
simplify the result and more easy to identify expression patterns.

```{r fig.height= 12, fig.width=15}
GO_UPReg_Heatplot <- plot( heatplot(GO_UPRegResults_epithelial))
saveplot(GO_UPReg_Heatplot, "GO_UPReg_Heatplot")
```

### Tree Plot of Enriched Terms

```{r fig.height= 10, fig.width=15}
edox2 <- pairwise_termsim(GO_UPRegResults_epithelial)
GO_UPReg_enrichtreeplot <- plot(treeplot(edox2))
saveplot(GO_UPReg_enrichtreeplot, "GO_UPReg_enrichtreeplot")
```

## GO Terms for Down Regulated Genes

### GO over-representation analysis for DOWN Regulated Genes

```{r}
DOWNgene_ENS_ID <- (sigsdf_epithelial_DOWN$ensemblID)

GO_DOWNRegResults_epithelial <- enrichGO(gene = DOWNgene_ENS_ID,
                       OrgDb = "org.Mm.eg.db",
                       keyType = "ENSEMBL",
                       ont = "BP",
                       pAdjustMethod = "BH",
                       pvalueCutoff  = 0.01,
                       qvalueCutoff  = 0.05,
                       readable      = TRUE
                       )
GO_DOWNRegdf_epithelial <- as.data.frame(GO_DOWNRegResults_epithelial)
head(GO_DOWNRegdf_epithelial)
```

```{r}
write.csv(GO_DOWNRegdf_epithelial,
          "/home/keshavprasadgubbi/Documents/AlinaRnaSeq/ER/GO_DOWNReg_results_Epithelial.csv")
```

```{r}
GO_DOWNReg_Barplot_epithelial <- plot(barplot(GO_DOWNRegResults_epithelial))
saveplot(GO_DOWNReg_Barplot_epithelial, "GO_DOWNReg_Barplot_epithelial")
```

```{r}
GO_DOWNReg_Dotplot_epithelial <- plot(dotplot(GO_DOWNRegResults_epithelial))

saveplot(GO_DOWNReg_Dotplot_epithelial, "GO_DOWNReg_Dotplot_epithelial")
```

```{r }
GO_DOWNReg_Cnetplot_epithelial <- plot(cnetplot(GO_DOWNRegResults_epithelial, 
                                                showCategory = 13, 
                                                circular = FALSE, 
                                                colorEdge = TRUE
                                                ))

saveplot(GO_DOWNReg_Cnetplot_epithelial, "GO_DOWNReg_Cnetplot_epithelial")
```

### Upset Plot

The upsetplot is an alternative to cnetplot for visualizing the complex
association between genes and gene sets. It emphasizes the gene
overlapping among different gene sets.

```{r}
GO_DOWNReg_Upsetplot <- plot( upsetplot(GO_DOWNRegResults_epithelial))
saveplot(GO_DOWNReg_Upsetplot, "GO_DOWNReg_Upsetplot")
```

### Heat Plot

```{r}
GO_DOWNReg_Heatplot <- plot(heatplot(GO_DOWNRegResults_epithelial))
saveplot(GO_DOWNReg_Heatplot, "GO_DOWNReg_Heatplot")
```

### Tree Plot of Enriched Terms

```{r}
edox1 <- pairwise_termsim(GO_DOWNRegResults_epithelial)
#GO_DOWNReg_enrichtreeplot <- plot(treeplot(edox1))
#saveplot(GO_DOWNReg_enrichtreeplot, "GO_DOWNReg_enrichtreeplot")
```

## Pathway Analysis

### KEGG Pathways

1.  The gageData package has pre-compiled databases mapping genes to
    KEGG pathways and GO terms for common organisms.
2.  kegg.sets.mm is a named list. Each element is a character vector of
    member gene Entrez IDs for a single KEGG pathway. sigmet.idx.hs is
    an index of numbers of sinaling and metabolic pathways in
    kegg.set.gs.
3.  In other words, KEGG pathway include other types of pathway
    definitions, like "Global Map" and "Human Diseases", which may be
    undesirable in pathway analysis. Therefore,
    kegg.sets.mm[sigmet.idx.mm] gives the "cleaner" gene sets of
    sinaling and metabolic pathways only.

The gage() function requires a named vector of fold changes, where the
names of the values are the Entrez gene IDs.

clusterProfiler requires the genes to be identified using Entrez IDs for
all genes in our results dataset. We also need to remove the NA values
and duplicates (due to gene ID conversion) prior to the analysis

```{r}
## Remove any NA values (reduces the data by quite a bit) and Remove any Entrez duplicates
sigsdf_epithelial_unique <- unique(dplyr::filter(sigsdf_epithelial, entrez != "NA"))
```

extract and name the fold changes

```{r}
foldchanges <-  sigsdf_epithelial$log2FoldChange
names(foldchanges) <-  sigsdf_epithelial$entrez

## Sort fold changes in decreasing order
foldchanges <- sort(foldchanges, decreasing = TRUE)
head(foldchanges)
```

```{r}
data(kegg.sets.mm)
data(sigmet.idx.mm)
kegg.sets.mm <-  kegg.sets.mm[sigmet.idx.mm]
```

Here, I am using same.dir = TRUE, which will give us separate lists for
pathways that are upregulated versus pathways that are downregulated.
Let's look at the first few results from each.

```{r}
# Get the results
keggres = gage(exprs = foldchanges, 
               gsets = kegg.sets.mm, 
               same.dir = TRUE)
```

```{r}
# Look at both up (greater), down (less), and statatistics.
lapply(keggres, head)

```

Now, process the results to pull out the top 20 up-regulated pathways,
then further process that just to get the IDs. I can use these KEGG
pathway IDs downstream for plotting.

```{r}
# Get the pathways

keggrespathways = data.frame(id = rownames(keggres$greater),
                             keggres$greater) %>% 
                  tibble::as_tibble() %>% 
                  filter(row_number() <= 10) %>% 
                  .$id %>% 
                  as.character()
keggrespathways
```

```{r}
# Get the IDs.
keggresids = substr(keggrespathways, start = 1, stop = 8)
keggresids
```

The pathview() function in the pathview package makes the plots. The
function is written below to loop through and draw plots for the top 10
pathways we created above.

```{r}
devtools::install_github("javadnoorb/pathview")
library("pathview")
# Define plotting function for applying later
plot_pathway = function(pid) pathview(gene.data = foldchanges, 
                                      pathway.id = pid, 
                                      species = "mouse", 
                                      new.signature = FALSE)


# plot multiple pathways (plots saved to disk and returns a throwaway list object)
tmp = sapply(keggresids, function(pid) pathview(gene.data = foldchanges, 
                                                pathway.id = pid, 
                                                species = "mouse"))
```
## GSEA - Gene Set Enrichment Analysis

```{r}
head(foldchanges)
```

Perform the GSEA using KEGG gene sets:

```{r}
# library("stats")
# ### GSEA using gene sets from KEGG pathways
# gseaKEGG <- gseKEGG(geneList = foldchanges, # ordered named vector of fold changes (Entrez IDs are the associated names)
# organism = "mmu", # mmu is Mus Musculus
# nPerm = 1000, # default number permutations
# minGSSize = 20, # minimum gene set size (# genes in set) - change to test more sets or recover sets with fewer # genes
# pvalueCutoff = 0.05, # padj cutoff value
# verbose = T)
# 
# ## Extract the GSEA results
# gseaKEGG_results <- gseaKEGG@result
```

How many pathways are enriched? View the enriched pathways:

## Write GSEA results to file

```{r}
# head(gseaKEGG_results)
# 
# write.csv(gseaKEGG_results, "/home/keshavprasadgubbi/Documents/AlinaRnaSeq/ER/gseaOE_kegg.csv")
```

Explore the GSEA plot of enrichment of one of the pathways in the ranked
list:

```{r fig.width=10, fig.height=6}
## Plot the GSEA plot for a single enriched pathway, `mmu04621`
#gseaplot(gseaKEGG, geneSetID = 'mmu04621')
```

Use the Pathview R package to integrate the KEGG pathway data from
clusterProfiler into pathway images:

```{r}
# detach("package:dplyr", unload = TRUE) # first unload dplyr to avoid conflicts
# 
# ## Output images for a single significant KEGG pathway
# pathview(gene.data = foldchanges,
#               pathway.id = "mmu04621",
#               species = "mmu",
#               limit = list(gene = 2, # value gives the max/min limit for foldchanges
#               cpd = 1))
```

# Gene Ontology - Hyper Geometric Test

```{r}
selectUPGenes_entrezid = unique(sigsdf_epithelial_UP['entrez'] )
selectDOWNGenes_entrezid = unique(sigsdf_epithelial_DOWN['entrez'] )
UniverseGenes = unique(sigsdf_epithelial$entrez)
cutOff = 0.01 # setting the cutoff at 1%
```

```{r}
upParams = new("GOHyperGParams",
               geneIds = selectUPGenes_entrezid,
               universeGeneIds = UniverseGenes,
               annotation = "org.Mm.eg.db",
               ontology = "BP",
               pvalueCutoff = cutOff,
               testDirection = "over")
downParams = new("GOHyperGParams",
               geneIds = selectDOWNGenes_entrezid,
               universeGeneIds = UniverseGenes,
               annotation = "org.Mm.eg.db",
               ontology = "BP",
               pvalueCutoff = cutOff,
               testDirection = "over")
```

```{r}
upBP = hyperGTest(upParams)
upBP_df <- as.data.frame(summary(upBP))
plot(upBP_df)
```

```{r}
htmlReport(upBP, file = "upBP.html")
```

```{r}
downBP = hyperGTest(downParams)
downBP_df <- as.data.frame(summary(downBP))
plot(downBP_df)
```
```{r}
htmlReport(downBP, file = "downBP.html")
```

# Signaling Pathway Impact Analysis

```{r}
foldchanges_vector <- as.vector(foldchanges)
```

```{r}
selectGenes_entrezid <- as.character(UniverseGenes)
```


```{r}
#spia(de = foldchanges_vector, all = selectGenes_entrezid, organism = "mmu")
```

The analysis with SPIA requires the conversion of gene-gene networks in a suitable format.
This conversion is performed by the function prepareSPIA that must be executed before the
analysis command runSPIA. The SPIA data will be saved in the current working directory;
every time you change it, you should also re-run prepareSPIA. Edges not included in SPIA
have been coerced into the admitted SPIA types. Metabolite mediated interactions (edges of
propagation through metabolites) annotated in graphite with the "indirect" type are mapped
into the SPIA edge type "indirect effect". For more details see the SPIA package [5, 6, 4].

```{r}

# library(SPIA)
# data(colorectalcancer)
# library(hgu133plus2.db)

```

```{r}
# top$ENTREZ <- mapIds(hgu133plus2.db, top$ID, "ENTREZID", "PROBEID", multiVals = "first")
# top <- top[!is.na(top$ENTREZ) & !duplicated(top$ENTREZ), ]
# top$ENTREZ <- paste("ENTREZID", top$ENTREZ, sep = ":")
# tg1 <- top[top$adj.P.Val < 0.05, ]
# DE_Colorectal = tg1$logFC
# names(DE_Colorectal) <- tg1$ENTREZ
# ALL_Colorectal <- top$ENTREZ

```

```{r}
sessionInfo()
```