L2fc.Rmd

---
title: ' L2FC table'
output:
  html_document: 
    toc: yes
    fig_width: 8
    fig_height: 8
    fig_caption: yes
    number_sections: yes
    toc_depth: 4
    df_print: tibble
---
<!-- # Log2FoldChnage Table  -->
<!-- # Author : Keshava Prasad Gubbi -->


```{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/countmatrices/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/countmatrices/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"

# Creating a list of samplenmes that can be provided
samplelist <- as.list(coldata$Sample_Name)

```

# Design a generic fucntion that does this job for everything else

```{r}
run_individualDESeq <- function(namelist, cd, cm) {
  #cd is coldata, cm is countsmatrix
  
  # Create an empty list for the results
  resultslist <- list()
  
  for (i in 1:14) {
    
    name = namelist[[i]]
    print(glue("Working with Strain:{name}"))
    
    cm_input <- cm[, c(i, 15:16)]
    coldata_input <- cd[c(i, 15:16), ]
    #print(head(cm_input, 3))
    #print(head(coldata_input, 3))
    
    rownames(coldata_input)
    colnames(cm_input)
    
    all(rownames(coldata_input) %in% colnames(cm_input))
    ncol(cm_input) == nrow(coldata_input)
    
    dds <- DESeqDataSetFromMatrix(countData = cm_input, colData = coldata_input, design = ~ condition)
    # nrow(dds)
    dds <- DESeq(dds)
    dds1 <- DESeq(dds, minReplicatesForReplace = Inf)
    
    res <- results(dds1, cooksCutoff = FALSE, independentFiltering = FALSE,
                   contrast = c("condition", "Infected", "control"))
    # head(res, 30)
    
    res_df <- as.data.frame(res) # convert the results table to a df
    # nrow(res_df)

    res_df <- tibble::rownames_to_column(res_df, "symbol")
    res_df <- res_df %>% filter(!is.na(log2FoldChange))
    # nrow(res_df)
    colnames(res_df)[which(names(res_df) == "log2FoldChange")] <- glue("log2FoldChange_{name}")
    res_df <- as.data.frame(res_df)
    print(glue("completed {name} strain DESeq calculations and produced result tables!!"))
  
    # Appending each new results table into a resultslist
    resultslist[[length(resultslist) + 1]] <- res_df
  }
  return(resultslist)
}
```


```{r}
# Calling the function to run DESeq for each Strain
resulttable <- run_individualDESeq(namelist = samplelist, cd = coldata, cm = countsmatrix)

# Mapping -> converting the list of objects to list of dataframes
list_df <- Map(as.data.frame, resulttable)

#merge all data frames in list based on symbol column
table <- list_df %>% reduce(full_join, by = 'symbol')
```

```{r}
l2fctable <- table %>% select(
  "symbol" ,
  "log2FoldChange_T",
  "log2FoldChange_S54",
  "log2FoldChange_S55",
  "log2FoldChange_L57",
  "log2FoldChange_L57",
  "log2FoldChange_A58",
  "log2FoldChange_L60" ,
  "log2FoldChange_S61",
  "log2FoldChange_A62",
  "log2FoldChange_L63",
  "log2FoldChange_A64",
  "log2FoldChange_S65",
  "log2FoldChange_L66",
  "log2FoldChange_A68",
  "log2FoldChange_L69"
)
```

```{r}
l2fctable
```

```{r}
write.csv(l2fctable, "/home/keshavprasadgubbi/Documents/AlinaRnaSeq/l2fctable_InfectedvsControl.csv")
```


# Calculating CPM Values


```{r}
run_individualDESeq2 <- function(namelist, cd, cm) {
  #cd is coldata, cm is countsmatrix
  
  # Create an empty list for the results
  cpmlist <- list()
  
  for (i in 1:14) {
    
    name = namelist[[i]]
    print(glue("Working with Strain:{name}"))
    
    cm_input <- cm[, i]
    coldata_input <- cd[i, ]
    #print(head(cm_input, 3))
    #print(head(coldata_input, 3))
    
    # CPM Values
    # as DGEList
    dge_er <- DGEList(counts = cm_input)
  
    dim(dge_er)
    colnames(dge_er)
    #dge_er$samples
  
    ## calculate norm. factors
    nr <- calcNormFactors(dge_er)
  
    ## get normalized counts
    cpm_df <- as.data.frame(cpm(nr))
    
    colnames(cpm_df) <- glue("CPM_{name}")
    
    cpm_df <- tibble::rownames_to_column(cpm_df,"symbol")
    
    # Appending each new result table into cpmlist
    cpmlist[[length(cpmlist) + 1]] <- cpm_df
  }
    return(cpmlist)
}
```


```{r}
cpmT <- run_individualDESeq2(namelist = samplelist, cd = coldata, cm = countsmatrix)
```
```{r}
# Mapping -> converting the list of objects to list of dataframes
list_df2 <- Map(as.data.frame, cpmT)

#merge all data frames in list based on symbol column
cpmtable <- list_df2 %>% reduce(full_join, by = 'symbol')
nrow(cpmtable)

rownames(cpmtable) <- cpmtable[,1]
cpmtable <- subset(cpmtable, select = -symbol)
keep <- rowSums(cpmtable) > 0
cpmtable <- cpmtable[keep,]
nrow(cpmtable)

```

```{r}
write.csv(cpmtable, "/home/keshavprasadgubbi/Documents/AlinaRnaSeq/CPMtable_BL6.csv")
```