analysis.rmd

---
title: "MicroarrayDataAnalysis"
author: "B243232"
date: "2024-04-06"
output: html_document
---

# The code is from basicanalysisscript.rmd provided in class with modifications made in certain parts

```{r setup, include=FALSE}
knitr::opts_chunk$set(eval = FALSE)
```

## Install all packages
```{r}
packages <- c("affy","simpleaffy","clusterProfiler","GEOquery","Biobase","mouse4302.db","enrichplot")

for (package in packages){
 BiocManager::install(package) 
}
```


## Loading dependencies
```{r message=FALSE, warning=FALSE}
library(affy)
library(simpleaffy)
library(limma)
#library(clusterProfiler)

library(annotate)
library(GEOquery)
library(Biobase)
library(mouse4302.db)

#library(enrichplot)
library(ggplot2)
library(scatterplot3d)
```

## Downloading data and setting up the enviroment
```{r message=FALSE,eval=FALSE}
# Taking input from user
repeat {
  geo.number <- readline(prompt="Please enter a valid GEO accession number (e.g., GSE12345): ")
  if (grepl("^GSE\\d+$", geo.number)) {
    break
  } else {
    cat("Invalid input. Please enter a GEO accession number in the correct format (e.g., GSE12345).\n")
  }
}
# Extracting the file contents
system("mkdir DownloadedFiles")
getGEOSuppFiles(geo.number, baseDir = "./DownloadedFiles")
tar.filename <- paste0("",geo.number,"_RAW.tar")
tar.command <- paste0("tar -xvf ./DownloadedFiles/",geo.number,"/",tar.filename," -C ./DownloadedFiles/")
system(tar.command)
Sys.sleep(3)
system("gunzip ./DownloadedFiles/*.gz ")
system("cp ./DownloadedFiles/* .")
system("rm -r ./DownloadedFiles")
system("mkdir Plots")
system("mkdir OutputTables")
```



> Need to provide code for creation of targets file

## Starting the Analysis

```{r message=FALSE,warning=FALSE}
raw.data <- ReadAffy() # Reading in .CEL files from the directory
adf.table <- read.table("GSE13074_targets.txt", header = TRUE, row.names = 1, as.is = TRUE) # Reading in targets file
annotated.table <- AnnotatedDataFrame(adf.table)
```


### Pre-Normalization Quality Assessment


```{r message=FALSE,warning=FALSE}
png("./PreNormDensity.png", width = 8, height = 8, units = 'in', res = 500)
#tiff("./Plots/PreNormDensity.tiff", width = 8, height = 8, units = 'in', res = 500, compression = 'lzw')
hist(raw.data,main="Pre-Normalization Density Plot")
```

> Plotting Pre Normalized Box Plot

```{r message=FALSE,warning=FALSE}
colors <- colors <- c(rep("#c20433",3),rep("#041ac2",3)) 
names <- c("DARCR1","DARCR2","DARCR3","LDRR1","LDRR2","LDRR3")

png("./plots/PreNormBox.png", width = 8, height = 8, units = 'in', res = 500)
#tiff("./Plots/PreNormBox.tiff", width = 8, height = 8, units = 'in', res = 500, compression = 'lzw')
boxplot(raw.data, col = colors, las=1,main="Pre-Normalization Box plot",names=names, xlab="Sample",ylab="Log Intensity")
```


> Plotting Pre Normalized MA plot

```{r message=FALSE,warning=FALSE}
png("./plots/PreNormMA.png", width = 12, height = 12, units = 'in', res = 500)
#tiff("./Plots/PreNormMA.tiff", width = 12, height = 12, units = 'in', res = 500, compression = 'lzw')
mva.pairs(pm(raw.data),main="Pre-Normalization MVA plot")
```


> Creating RNA degradation plot and QC statistic plot using simpleaffy

```{r}
sampleNames(raw.data) <- c("DARCR1","DARCR2","DARCR3","LDRR1","LDRR2","LDRR3")
deg <- AffyRNAdeg(raw.data)
#png("./Plots/RNAdeg.png", width = 8, height = 8, units = 'in', res = 500)
png("./plots/RNAdeg.png", width = 10, height = 10, units = 'in', res = 500)
plotAffyRNAdeg(deg, cols = c(2,2,2,3,3,3))
```

```{r}
qc.data <- qc(raw.data) #Simple Affy function
```

```{r warning=FALSE}
#png("./Plots/QCStats.png", width = 8, height = 8, units = 'in', res = 500)
png("./plots/QCStats.png", width = 10, height = 10, units = 'in', res = 500)
plot(qc.data)
```

> Calcute the average background, scaling factor and percent calls for cross checking the QC statistics plot

```{r}
# average background
avbg(qc.data)
```

```{r}
# scaling factor
sfs(qc.data)
```

```{r}
# percent calls
percent.present(qc.data)
```


```{r}
# Calculating the percent calls
ratios(qc.data)[,1:4]
```



### Normalization

```{r message=FALSE,warning=FALSE}
norm.expression <- rma(raw.data, normalize = TRUE, background = TRUE, verbose = TRUE) # Performing RMA normalization on the expression set object (and converting the affy batch object to an expression set object)
norm.matrix <- exprs(norm.expression) # Extracting the normalized expression matrix values
```

It is worth nothing that *norm.expression* is an expression set object and *norm.matrix* is the expression matrix values.


### Post Normalization Quality Assessment

> Post Normalization Box Plot

```{r message=FALSE,warning=FALSE}
png("./plots/PostNormBox.png", width = 8, height = 8, units = 'in', res = 500)
#tiff("./Plots/PostNormBox.tiff", width = 8, height = 8, units = 'in', res = 500, compression = 'lzw')
boxplot(norm.matrix, col = colors, las=1,main="Post-Normalization Box plot",names=names, xlab="Sample",ylab="Log Intensity")
```

> Post Normalization MA Plot

```{r message=FALSE,warning=FALSE}
png("./plots/PostNormMA.png", width = 8, height = 8, units = 'in', res = 500)
mva.pairs(norm.matrix,main="Post-Normalization MVA plot")
```


### Exploratory Data Analysis - Clustering and Dimensionality Reduction

#### Annotation and Heirarchical Clustering

```{r message=FALSE,warning=FALSE}
colnames(norm.matrix) <- rownames(pData(annotated.table))
hc<-hclust(as.dist(1-cor(norm.matrix, method="pearson")), method="average")
png("./plots/Dendrogram.png", width = 8, height = 8, units = 'in', res = 500)
plot(hc)
```

#### Principal Component Analysis

```{r message=FALSE,warning=FALSE}
pca <- prcomp(t(norm.matrix), scale=T)
png("./plots/PCAscatter.png", width = 8, height = 8, units = 'in', res = 500)
s3d<-scatterplot3d(pca$x[,1:3], pch=19, color=rainbow(1))
s3d.coords <- s3d$xyz.convert(pca$x[,1:3])
text(s3d.coords$x, s3d.coords$y, labels = colnames(norm.matrix),pos
= 3,offset = 0.5)
```

#### Updated way of visualizing the PCA results
```{r fig.height=8, fig.width=8}
# Perform PCA
pca <- prcomp(t(norm.matrix), scale = TRUE)

# Calculate variance explained
variance_explained <- (pca$sdev^2) / sum(pca$sdev^2) * 100
variance_explained_df <- data.frame(
    PC = paste0("PC", 1:length(variance_explained)),
    Variance = variance_explained
)

# Create a data frame with PCA scores and sample names
pca_scores <- as.data.frame(pca$x[, 1:2])  # Extract PC1 and PC2
colnames(pca_scores) <- c("PC1", "PC2")
pca_scores$Sample <- colnames(norm.matrix)  # Add sample names

# 2D PCA Scatter Plot
pca_plot <- ggplot(pca_scores, aes(x = PC1, y = PC2, label = Sample)) +
    geom_point(size = 3, color = "blue") +
    geom_text(vjust = -1, size = 3) +
    labs(
        title = "PCA of the Samples",
        x = paste0("PC1 (", round(variance_explained[1], 1), "% Variance)"),
        y = paste0("PC2 (", round(variance_explained[2], 1), "% Variance)")
    ) +
    theme_classic()
```

```{r}
ggsave(filename = "./plots/PCA_Plot.png",plot=pca_plot,height = 6, width = 8, units = "in", dpi = 300)
```



### Differential Gene Expression Analysis

#### Build the Fold Change Table

```{r message=FALSE,warning=FALSE}
exprs.vals <- exprs(norm.expression)
exprs.vals10 <- 2^exprs.vals

# Extracting sample names and probe ids
sample.names <- sampleNames(norm.expression) 
probe.ids <- featureNames(norm.expression) 

entrez.ids <- mapIds(mouse4302.db, keys=probe.ids, column="ENTREZID", keytype="PROBEID", multiVals="first") # Mapping the probeIDs to EntrezIds
```

```{r message=FALSE,warning=FALSE}
# Calculating the group means
DARCR.mean <- apply(exprs.vals10[,c("DARCR1","DARCR2","DARCR3")],1,mean)
LDRR.mean <- apply(exprs.vals10[,c("LDRR1","LDRR2","LDRR3")],1,mean)

# Calculating the fold change
LDRR.DARCR.FC <- LDRR.mean/DARCR.mean # Fold change of LDRR wrt DARCR

# Creating and Saving the fold change table
fc.summary.table <- cbind(DARCR.mean,LDRR.mean,LDRR.DARCR.FC)
write.table(fc.summary.table, file = "./FoldChangeTable.txt",quote=FALSE, sep="\t")
```

> Annotating the expression set object *(norm.expression)* using feature data

```{r message=FALSE,warning=FALSE}
sampleNames(norm.expression) <- c("DARCR1","DARCR2","DARCR3","LDRR1","LDRR2","LDRR3")

ID <- featureNames(norm.expression)
Symbol <- getSYMBOL(ID, "mouse4302.db")
Name <- as.character(lookUp(ID, "mouse4302.db", "GENENAME"))

tmp <- data.frame(ID=ID, EntrezID=entrez.ids, Symbol=Symbol, Name=Name, stringsAsFactors=F) # Creating the metadata data frame
tmp[tmp=="NA"] <- NA

fData(norm.expression) <- tmp # Assigning the feature data to the expressionset object

```

#### Statistical analysis using limma

```{r message=FALSE,warning=FALSE}
# Creating the design matrix
design.matrix <- model.matrix(~-1+factor(c(1,1,1,2,2,2)))
colnames(design.matrix) <- c("DARCR","LDRR")

# Creating the contrast matrix
contrast.matrix <- makeContrasts(LDRR-DARCR,
                                levels=design.matrix)

## Fitting the linear model
fit <- lmFit(norm.expression, design.matrix)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
```

> Using topTable() for differential gene expression analysis

```{r message=FALSE,warning=FALSE}
toptable.results <- topTable(fit2,coef=1, adjust="fdr", number=nrow(norm.expression))
toptable.results <- toptable.results[!is.na(toptable.results$Symbol) & toptable.results$Symbol != "", ]
write.table(toptable.results,"./topTableResults.txt",sep="\t")
```

> Using topTreat() for differential gene expression analysis

```{r}
fit2.treat <- treat(fit2)
```

```{r message=FALSE,warning=FALSE}
toptreat.results <- topTreat(fit2.treat, coef=1, adjust.method="fdr", number=nrow(norm.expression))
toptreat.results <- toptreat.results[!is.na(toptreat.results$Symbol) & toptreat.results$Symbol != "", ]
write.table(toptreat.results,"./topTreatResults.txt",sep="\t")
```


> Creating volcano plot using ggplot2

# The code to create this volcano plot was inspired by https://biostatsquid.com/volcano-plots-r-tutorial/#step1

```{r message=FALSE,warning=FALSE}
# The code to create this volcano plot was inspired by https://biostatsquid.com/volcano-plots-r-tutorial/#step1
toptable.results$neg.log10.pva.adj <- -log10(toptable.results$adj.P.Val) # Using adjusted p value rather than p value

volcano.plot <- ggplot(toptable.results, aes(x = logFC, y = neg.log10.pva.adj)) +
  geom_point(aes(color = adj.P.Val < 0.05 & abs(logFC) > 1), show.legend = FALSE) +
  scale_color_manual(values = c("black", "red")) +
  labs(title = "", x = "Log Fold Change", y = "-Log10 Adjusted P-value") +
  geom_vline(xintercept = c(-1, 1), color = "green", linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), color = "blue", linetype = "dashed") +
  theme_minimal()

volcano.plot <- volcano.plot +
  coord_cartesian(xlim = c(-6, 6), ylim = c(0, 6)) +
  scale_x_continuous(labels = scales::comma) +
  scale_y_continuous(labels = scales::comma)

top.genes <- toptable.results[order(toptable.results$adj.P.Val, -abs(toptable.results$logFC)), ]
top.genes <- head(top.genes, 5)

volcano.plot <- volcano.plot +
  geom_text(data = top.genes, aes(label = Symbol), vjust = 0, hjust = -0.2, color = "blue", position = position_nudge(x = 0.05), check_overlap = FALSE, angle = 45, size = 2.8) + theme_classic()

ggsave("./plots/VolcanoDEG.png", plot = volcano.plot, width = 8, height = 8, unit="in",dpi = 300)

```


> Creating volcano plot using limma

```{r message=FALSE,warning=FALSE}
# Volcano Plot using limma
png("./plots/VolcanoDEGLim.png", width = 10, height = 10, units = 'in', res = 500)
volcanoplot(fit2, coef=1, names=toptable.results$genes)
abline(h=-log10(0.05), col="blue", lty=2)
abline(v=c(-1, 1), col="red", lty=2)

```

### Functional Enrichment Analysis

#### Pathway Enrichment Analysis using Camera
```{r message=FALSE, warning=FALSE}
# Loading the Gene Set
gene.set.rds <- readRDS("./Mm.h.all.v7.1.entrez.rds") # Class Provided Gene Set

res <- select(mouse4302.db, keys = rownames(norm.expression), columns = c("ENTREZID", "ENSEMBL","SYMBOL"), keytype="PROBEID")
idx <- match(rownames(norm.expression), res$PROBEID)

fData(norm.expression) <- res[idx, ]
eset_t<-norm.expression[is.na(fData(norm.expression)$ENTREZID)==0,]

H.indices <- ids2indices(gene.set.rds,fData(eset_t)$ENTREZID)
camera.enrichment <-camera(eset_t,index=H.indices,design=design.matrix,contrast=contrast.matrix[,1],adjust.method = "BH")

write.table(camera.enrichment, "./output_files/CameraEnrichment.txt", sep="\t", quote=FALSE)
write.csv(camera.enrichment[order(camera.enrichment$FDR, decreasing = FALSE),], "./output_files/SortedCameraEnrichment.csv")
```


### Gene Ontology Enrichment (Not compatible with R 4.0.5)

```{r message=FALSE,warning=FALSE}
# Selecting only the most significant genes
gene.en.id <- toptable.results$EntrezID[toptable.results$adj.P.Val < 0.05]

# GO enrichment for biological function
go.enrich.result <- enrichGO(gene = gene.en.id, OrgDb = org.Mm.eg.db, keyType = "ENTREZID", ont = "BP", pAdjustMethod = "fdr", qvalueCutoff = 0.05, minGSSize = 10, maxGSSize = 100)
```

```{r message=FALSE,warning=FALSE}
# Saving the plot
tiff("./plots/GOenrichment.tiff", width = 11, height = 11, units = 'in', res = 500, compression = 'lzw')
barplot(go.enrich.result,showCategory=20)

# Saving the results for future use
write.table(go.enrich.result, "./OutputTables/GOEnrichment.txt", sep="\t", quote=FALSE)
```


#### Top 10 significantly expressed genes by topTable()

```{r}
top.10.genes.by.limma <- toptable.results[order(toptable.results$adj.P.Val, decreasing = FALSE),]
top.10.genes.by.limma <- head(top.10.genes.by.limma, 10)
top.10.genes.by.limma.list <- top.10.genes.by.limma$Symbol
write.csv(top.10.genes.by.limma, "./Top10DEGLM.csv", row.names = FALSE)
```

```{r}
top.10.genes.by.fc <- toptable.results[order(toptable.results$logFC, decreasing = TRUE),]
top.10.genes.by.fc <- head(top.10.genes.by.fc, 10)
top.10.genes.by.fc.list <- top.10.genes.by.fc$Symbol
write.csv(top.10.genes.by.fc, "./Top10DEGFC.csv", row.names = FALSE)
```

#### Top 10 significantly expressed genes by topTreat()

```{r}
top.10.genes.by.limma.t <- toptreat.results[order(toptreat.results$adj.P.Val, decreasing = FALSE),]
top.10.genes.by.limma.t <- head(top.10.genes.by.limma.t, 10)
top.10.genes.by.limma.list.t <- top.10.genes.by.limma.t$Symbol
write.csv(top.10.genes.by.limma.list.t, "./Top10DEGLMTr.csv", row.names = FALSE)
```

```{r}
top.10.genes.by.fc.t <- toptreat.results[order(toptreat.results$logFC, decreasing = TRUE),]
top.10.genes.by.fc.t <- head(top.10.genes.by.fc, 10)
top.10.genes.by.fc.list.t <- top.10.genes.by.fc$Symbol
write.csv(top.10.genes.by.fc.t, "./Top10DEGFCTr.csv", row.names = FALSE)
```