DGE_Alina.Rmd

---
title: "Differential Gene Expression Analysis Notebook - Alina's Data Analysis Project"
output:
  pdf_document:
    toc: yes
  html_document:
    toc: yes
    number_sections: yes
---

                                             
The DGE analysis will be performed using the raw integer read counts for control and infected conditions. The goal here is to identify the differentially expressed genes under infected condition.

```{r message=F, warning=F, paged.print=TRUE}
library("DESeq2")
library("PoiClaClu")
library("pheatmap")
library("RColorBrewer")
library('tidyverse')
library("PoiClaClu")
library("vsn")
library('EnhancedVolcano')
```

```{r}
cm1 <-as.matrix(read.csv(
      "C:/Users/kesha/Documents/R/Rtuts/Data/WG__Discussion_Alina's_EPEC_project/fc/counts.csv",
      row.names = "X"
    ))
dim(cm1)
```


```{r}
head(cm1)
```


DESeq2 needs sample information (metadata) for performing DGE analysis. Let’s create the sample information (you can also import sample information if you have it in a file).

```{r}
sample_ID <- list(read.csv("C:/Users/kesha/Documents/R/Rtuts/Data/WG__Discussion_Alina's_EPEC_project/fc/samples.csv", row.names = "X"))
sample_ID
```


```{r}
condition = c("Infected","Infected","Infected","Infected","Infected","Infected","Infected","Infected","Infected","Infected","Infected","Infected","Infected","Infected","Infected","Infected","Infected","Infected","Infected","control","control")
```


```{r}
coldata <- data.frame(sample_ID, condition)
colnames(coldata) <- c('sample','condition') # change name of one of the columns
coldata
```

```{r}
library(stringr)
coldata$sample = str_remove_all(coldata$sample, pattern = '.markdup.sorted') # tidying up the
coldata
```


```{r}
coldata$condition <- as.factor(coldata$condition)
coldata <- data.frame(coldata, row.names = 1) # convert column to row.names
coldata 
```

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}
colnames(cm1) = str_remove_all(colnames(cm1), pattern = '.markdup.sorted') #tidying up these names again
colnames(cm1)
```


```{r}
all(rownames(coldata) %in% colnames(cm1))
```
Now, construct DESeqDataSet for DGE analysis:


```{r}
dds <- DESeqDataSetFromMatrix(countData = cm1, 
                              colData = coldata, 
                              design = ~ condition)
```

```{r}
dds
```

# Normalization


## Count Normalization via DESq2

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

Can check the normalization factors using:
```{r}
sizeFactors(dds)
```

```{r}
normalized_counts <- counts(dds, normalized = TRUE)
```

Saving the normalized counts table.
```{r}
write.table(normalized_counts, 
          file="C:/Users/kesha/Documents/R/Rtuts/Data/WG__Discussion_Alina's_EPEC_project/results/normalized_counts.txt",
          sep = "\t",
          quote = FALSE,
          col.names = NA)
```

ASK Q1!


# Quality Control


Assess overall similarity between samples and understand the following:

1. Which samples are similar to each other and which are different?
2. Does this fit to the expectation from experiment's design?
3. What are the major sources of variation in dataset?

## Count Data Tansformations - Extracting Transformed Values


### Transform normalised counts using rlog transformation for Visualization (only!)

```{r}
rld <- rlog(dds, blind = TRUE) 
# blind flag turned on to keep the transformation unbiased to sample condition information
vsd <- vst(dds, blind = TRUE)

```


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

### Effects of transformations on the variance


```{r}
# this gives log2(n + 1)
ntd <- normTransform(dds)

meanSdPlot(assay(ntd))
```

```{r}
meanSdPlot(assay(vsd))
```
```{r}
meanSdPlot(assay(rld))
```
## Data Quality assessment via sample clustering and visualization


### Heatmap of Count Matrix

```{r}

select <- order(rowMeans(counts(dds,normalized=TRUE)),
                decreasing=TRUE)[1:21]
df <- as.data.frame(colData(dds))
```


```{r}

pheatmap(assay(ntd)[select,], cluster_rows=TRUE, show_rownames=TRUE,
         cluster_cols=FALSE, annotation_col=df)
```


```{r}
pheatmap(assay(vsd)[select,], cluster_rows=TRUE, show_rownames=TRUE,
         cluster_cols=FALSE, annotation_col=df)
```

### Heatmap of sample-to-sample distances using the variance stabilizing transformed values.


```{r}
sampleDists <- dist(t(assay(vsd)))

sampleDistMatrix <- as.matrix( sampleDists )
rownames(sampleDistMatrix) <- names(vsd$sizeFactor)
colnames(sampleDistMatrix) <- NULL
```


```{r}

colors <- colorRampPalette( rev(brewer.pal(9, "Spectral")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         col= colors
         )
```
### Plot PCA of Samples

```{r}
plotPCA(vsd, intgroup="condition")
```


### Hierarchical Clustering

```{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
```


```{r}
### 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
pheatmap(rld_cor, border_color=NA, 
         fontsize = 8, 
         fontsize_row = 10, 
         height=50)
```

```{r}
heat.colors <- brewer.pal(6, "Spectral")
pheatmap(rld_cor, color = heat.colors, border_color=NA, fontsize = 7, 
			fontsize_row = 10, height=50)
```


# DGE Analysis

Pre-filter the genes which have low counts. Here, I will remove the genes which have < 10 reads (this can vary based on research goal) in total across all the samples. Pre-filtering helps to remove genes that have very few mapped reads, reduces memory, and increases the speed of the DESeq2 analysis.

```{r}
dds <- dds[rowSums(counts(dds)) >= 1,]
```

Now, select the reference level for condition comparisons. The reference level can set using ref parameter. The comparisons of other conditions will be compared against this reference i.e, the log2 fold changes will be calculated based on ref value (infected/control) . If this parameter is not set, comparisons will be based on alphabetical order of the levels.

```{r}
dds$condition <- relevel(dds$condition, ref = "control")
```

Perform Differntial Gene Expression Analysis:

```{r}
dds1 <- DESeq(dds)
```

```{r}
resultsNames(dds1)
```

### Result 

```{r}
res <- results(dds1, independentFiltering=FALSE, name = "condition_Infected_vs_control")
res
```

Order gene expression table by adjusted p value (Benjamini-Hochberg FDR method)

```{r}
res_adj <- res[order(res$padj),]  
res_adj
```


### Note on NA:  some genes with p value set to NA. 

1. This is due to all samples have zero counts for a gene ----->>> baseMean is 0 and all other values are NA.
2. There is a row with extreme outlier count for a gene (determined by cook's distance) ----->>> pvalue set to NA.

Cook's distance : DESeq function calculates, for every gene and for every sample, a diagnostic test for outliers called Cook’s distance. Cook’s distance is a measure of how much a single sample is influencing the fitted coefficients for a gene, and a large value of Cook’s distance is intended to indicate an outlier count. The Cook’s distances are stored as a matrix available in assays(dds)[["cooks"]].

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.



```{r}
write.csv(as.data.frame(res[order(res$padj),] ), 
          file="C:/Users/kesha/Documents/R/Rtuts/Data/WG__Discussion_Alina's_EPEC_project/results/condition_infected_vs_control_dge.csv")
```

```{r}
summary(results(dds1, alpha=0.05))
```
```{r}
r<- results(dds1, alpha=0.05, lfcThreshold =1)
r
```


```{r}
summary(results(dds1, alpha=0.05, lfcThreshold =1))
```

```{r}
ir <- results(dds1, alpha=0.05, lfcThreshold =1)
ir
```


## Results - with padj and LFC thresholds

```{r}
res.05 <- results(dds1, alpha = 0.05)
table(res.05$padj < 0.05)
```

```{r}
resLFC1 <- results(dds1, lfcThreshold=1)
table(resLFC1$padj < 0.05)
```

```{r}
new_res <- results(dds1, alpha = 0.05, lfcThreshold=1)
new_res
```


### Order gene expression table by adjusted p value (Benjamini-Hochberg FDR method)
```{r}
#
new_res_adj <- new_res[order(new_res$padj),]  
new_res_adj

```
saving these results table :

```{r}
write.csv(as.data.frame(new_res_adj), 
          file="C:/Users/kesha/Documents/R/Rtuts/Data/WG__Discussion_Alina's_EPEC_project/results/condition_infected_vs_control_dge_withthresholdsimposed.csv")
```



# Exploratory Data analysis

```{r}
plotMA(new_res)
```



```{r}
plotDispEsts(dds1)
```

### Calculating sample distances is to use the Poisson Distance 

This measure of dissimilarity between counts also takes the inherent variance structure of counts into consideration when calculating the distances between samples. The PoissonDistance function takes the original count matrix (not normalized) with samples as rows instead of columns.

```{r}

poisd <- PoissonDistance(t(counts(dds1)))
print(poisd)
```


```{r}
samplePoisDistMatrix <- as.matrix( poisd$dd )
dim(samplePoisDistMatrix)
#rownames(samplePoisDistMatrix) <- paste( dds1$dex, dds1$cell, sep=" - " )
#colnames(samplePoisDistMatrix) <- NULL

```

```{r}

pheatmap(samplePoisDistMatrix,
         clustering_distance_rows = poisd$dd,
         clustering_distance_cols = poisd$dd,
         col = colors)
```



## Log fold change shrinkage for visualization and ranking

```{r}
resLFC <- lfcShrink(dds1, coef="condition_Infected_vs_control", type="apeglm")
resLFC
```

```{r}
plotMA(resLFC)
```



# Plotting Results

## Counts Plot

A quick way to visualize the counts for a particular gene is to use the plotCounts function that takes as arguments the DESeqDataSet, a gene name, and the group over which to plot the counts

```{r}
topGene <- rownames(res)[which.max(res_adj$padj)]
topGene
```

```{r}
plotCounts(dds1, gene = topGene, intgroup="condition")
```
### Normalized counts for a single gene over treatment group
```{r}
library("ggbeeswarm")
geneCounts <- plotCounts(dds1, gene = topGene, 
                         intgroup="condition",
                         returnData = TRUE)
geneCounts
```
```{r}
ggplot(geneCounts, aes(x = condition, y = count, color = condition)) +
  scale_y_log10() +  geom_beeswarm(cex = 3)
```


# MDS Plot

# MA Plot

## Gene Clustering

Sample distance heatmap made previously, the dendrogram at the side shows us a hierarchical clustering of the samples. Such a clustering can also be performed for the genes.For demonstration, let us select the 20 genes with the highest variance across samples.

```{r}
library("genefilter")
topVarGenes <- head(order(rowVars(assay(vsd)), decreasing = TRUE), 20)
topVarGenes
```
The heatmap becomes more interesting if we do not look at absolute expression strength but rather at the amount by which each gene deviates in a specific sample from the gene’s average across all samples. Hence, we center each genes’ values across samples, and plot a heatmap (figure below). We provide a data.frame that instructs the pheatmap function how to label the columns.
```{r}
mat  <- assay(vsd)[ topVarGenes, ]
mat  <- mat - rowMeans(mat)
print(mat)
anno <- as.data.frame(colData(vsd)[, "condition"])
```
```{r}
all(rownames(anno) %in% colnames(assay(vsd)))
# this needs to match in order for the genen clustering to take place

```

### Volcano Plots

```{r}
EnhancedVolcano(new_res,
    lab = rownames(new_res),
    x = 'log2FoldChange',
    y = 'pvalue',
    
    ylim = c(0, 40),
    pCutoff = 0.05,
    FCcutoff = 1.0,
    pointSize = 3.0,
    labSize = 3.0,
    labCol = 'black',
    labFace = 'bold',
    boxedLabels = TRUE,
    colAlpha = 1,

    legendLabSize = 10,
    legendIconSize = 4.0,
    drawConnectors = TRUE,
    widthConnectors = 0.75)
```