RNA_SEQ_PIPELINE_2024.Rmd

---
title: "RNASeq - Pipeline"
author: "Sara Nicholson"
date: "2024-04-18"
output: html_document
---

```{r}
setwd("/storage1/fs1/leyao.wang/Active/saran")
```

# Quantification
```{r}
library(Rsubread)

path = "/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/RSEQC/new/RSV"   # RSV
path = "/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/RSEQC/new/LACT"  # LACT
path = "/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/RSEQC/new/MOCK"  # MOCK
path = "/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/RSEQC/new/COINF" # COINF

path = "/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/counts/"

# List demultiplexed sample fastq files
samples <- list.files(path)  
bam <- sort(list.files(path, pattern = "*csv$", full.names = TRUE))
bam

count = 1
# Run Feature Counts
for (i in bam) {
    counts <- featureCounts(files = i, isPairedEnd=TRUE, GTF.featureType="exon", GTF.attrType="gene_id", annot.ext= "/storage1/fs1/leyao.wang/Active/Users/run/rna.attempts/genomeDir/gencode.v42.primary_assembly.annotation.gff3", isGTFAnnotationFile=TRUE, countMultiMappingReads = FALSE, countReadPairs=T)
  get.sample.name <- function(fname) (basename(fname))
remove.ext <- function(fname) (tools::file_path_sans_ext(fname, compression=TRUE))
sample.names <- unname(sapply(bam, get.sample.name))
sample.names <- unname(sapply(sample.names, remove.ext))
sample.names <- substr(sample.names, 1, 7 )
  write.csv(counts$counts,  paste('/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/counts', sample.names[count], 'csv', sep = '.'))
  count = count + 1
}

#### edit to replace script below
path = "/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/counts"
files <- list.files(path, full.names = T)
features <- read.csv(files[1], header=T, sep=",")[,1]   # gene names
df    <- do.call(cbind,lapply(files,function(fn)read.csv(fn,header=T, sep=",")[,2]))
df <- cbind(features,df)
colnames(df) <- c("features","Co1", "Co2", "Co3", "Co4", "Lac1", "Lac2", "Lac3", "Lac4", "Mock1", "Mock2", "Mock3", "Mock4", "RSV1", "RSV2", "RSV3", "RSV4")
as.data.frame(df)

write.csv(df, "/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/Count_Table.csv")
```

Combine all Quantified CSV Files into one chart, You can do this within R or in excel/libreOffice Calc as I have done; Samples should be columns and genes/features as rows. 

# EdgeR #

### Initializations
```{r}
library(edgeR)
library(dplyr)
library(tibble)
```

**Read Counts File**
```{r}
# Read in your count table
x <- read.csv("/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/Count_Table.csv", row.names="features")
head(x)
x <- x %>% dplyr::select(-X)
### Subset to desired samples for Different Comparisons
x <- x %>% dplyr::select("Mock1", "Mock2", "Mock3", "Mock4", "RSV1", "RSV2", "RSV3", "RSV4")
x <- x %>% dplyr::select("Mock1", "Mock2", "Mock3", "Mock4", "Co1", "Co2", "Co3", "Co4")
x <- x %>% dplyr::select("Mock1", "Mock2", "Mock3", "Mock4", "Lac1", "Lac2", "Lac3", "Lac4")
x <- x %>% dplyr::select("RSV1", "RSV2", "RSV3", "RSV4", "Co1", "Co2", "Co3", "Co4")

library(biomaRt)


biomart_fx <- function(x){
   # Create vector of ensembl gene IDs
ENSG <- as.vector(row.names(x))

# define biomart object
mart <- useMart(biomart = "ensembl", dataset = "hsapiens_gene_ensembl")
# query biomart - retrieve corresponding gene names to ENSG
getBM(attributes = c("ensembl_gene_id_version", "external_gene_name", "ensembl_gene_id"),
               filters = "ensembl_gene_id_version", values = ENSG,
               mart = mart)
}

results <- biomart_fx(x)

# Write conversion chart to CSV
write.csv(results, "/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/ENSG_symbol_conversions.csv")
```

**Create DGElist object for EDgeR**
```{r}
# Create Factored Groups from Column Names 
# Must be in same order as column names on count table
DataGroups <- c("co", "co", "co", "co", "lac", "lac", "lac",  "lac","mock", "mock", "mock", "mock","rsv", "rsv", "rsv", "rsv" )
DataGroups <- c(rep("mock", 4), rep("rsv", 4)) 


# Create DGEList for EdgeR
d <- DGEList(counts=x,group=factor(DataGroups))
d$samples
```

### Filter Out Lowly Expressed Genes
Use the cpm function & remove features that are not expressed over 1 CPM in at least 2 samples.
```{r}
# Check Prior Dimensions
dim(d)

# Function to remove features not expressed over 1 CPM in at least 3 samples
keep <- rowSums(cpm(d)>10) >=3 
d <- d[keep,]

# Check Dimensions after Filtering
dim(d)
```

```{r}
### Edge R recommended Filtering
# keep.exprs <- filterByExpr(d, group=d$samples$group)
# d <- d[keep.exprs,, keep.lib.sizes=FALSE]
# dim(d)
```

*17586 genes remain in analysis*

### Normalize for Library Size: TMM Normalization.
*'If a small proportion of highly expressed genes consume a substantial proportion of the total library size for a particular sample, this will cause the remaining genes to be under-sampled for that sample.....The calcNormFactors function normalizes the library sizes by finding a set of scaling factors for the library sizes that minimizes the log-fold changes between the samples for most genes.' -EdgeR*
```{r}
# set method to TMM or CPM normalization will be implemented
d <- calcNormFactors(d, method = "TMM")
# view Normalization Factors
d$samples
dim(d)

# use cpm() function to get normalized counts - applies normalization factors so it is TMM when in combination with above function
dtmm <- as.data.frame(cpm(d))
head(dtmm)

write.csv(dtmm,"/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/ANALYSIS/dtmm_ALLSAMP.csv")
```

### PCoA Plot + Relative Log Expression
```{r}
library(EDASeq)
x <- as.factor(DataGroups)
set <- newSeqExpressionSet(as.matrix(dtmm),
                           phenoData = data.frame(x, row.names=colnames(d)))
set

library(RColorBrewer)
colors <- brewer.pal(3, "Set2")
plotRLE(set, outline=FALSE, ylim=c(-4, 4), col=colors[x])

plotPCA(set, col=colors[x], cex=1.2)
```

# Limma #
### Initializations
```{r}
library(limma)
library(Glimma)
library(dplyr)
library("AnnotationDbi")
library("org.Hs.eg.db")
library(biomaRt)
```

**Unsupervised Clustering of Samples**
```{r}
# layout(matrix(1:2, nrow = 1))
# MDS/PCoA Plots
layout(matrix(1:2, nrow = 1))
plotMDS(d, col = as.numeric(d$samples$group))
plotMDS(d, col = as.numeric(d$samples$group), gene.selection = "common", cex = 1, pch = 19)

####################### Another PCA plot with Clusters surrounded
library(factoextra)
library(ggplot2)
library(ggforce)

dtmm1 <- t(as.matrix(d$counts))
dtmm1 <- as.data.frame(dtmm1)
dtmm1$samp <- c("mock", "mock", "mock", "mock", "lac", "lac", "lac", "lac")
pca <- prcomp(dtmm1[,1:293],  
                   scale = F)
 
summary(pca)

fviz_pca_ind(pca, 
             habillage=dtmm1$samp, label = 'none', repel = T) +  ggforce::geom_mark_ellipse(aes(fill = Groups,
                        color = Groups)) +
  theme(legend.position = 'bottom') +
  coord_equal()


```



```{r}
# Create Design Matrix (Model by Group : RSV vs. MOCK)
design <- model.matrix(~d$samples$group)
rownames(design) <- colnames(d)
design

#### Edit if needed
colnames(design) <- c("intercept", "Co-Inf")
design[1:4,2] <- 0
design[5:8,2] <- 1
design

####### Contrast Matrices when comparing all samples at once
colnames(design) <- c("co","lac", "mock","rsv")
contr.matrix <- makeContrasts(
   RSVvsMock = rsv - mock,
   RSVvsLac = rsv - lac,
   LacvsMock = lac - mock,
   COvsMock = co - mock,
   COvsLac = co - lac,
   COvsRSV = co - rsv,
   levels = colnames(design))
contr.matrix
```

### Differential Expression- Limma
```{r}
# Calculate Weighted Likelihoods, prepare to be linearly modeled
v <- voomWithQualityWeights(d, design, plot = TRUE)
v

# Write VoomWithQualityWeights output to CSV for future input into WGCNA
wgcna_v <- as.data.frame(v$E)
write.csv(wgcna_v, "/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/ANALYSIS/wgcna_RSVMOCK.csv")

# linear model
# fit linear model
vfit <- lmFit(v, design)
# vfit <- contrasts.fit(vfit, contrasts=contr.matrix)

# Empirical Bayes Transform
vfit <- eBayes(vfit)
plotSA(vfit)
# Variance no longer dependent on mean

# Check # of up + down regulated genes
dt <- decideTests(vfit)
summary(dt)

# TOP Differential expressed GENES
topTable(vfit, coef=2, sort.by = "P")
top <- topTable(vfit,coef=2,number=Inf,sort.by="P", adjust.method = "fdr")

conversions <- read.csv("/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/ENSG_symbol_conversions.csv")
top_gene_table <- left_join(tibble::rownames_to_column(top), results, by=c("rowname" = "ensembl_gene_id_version"))

top_gene_table <- top_gene_table[which(top_gene_table$adj.P.Val < 0.05),]
top_up <- top_gene_table[which(top_gene_table$adj.P.Val < 0.05 & top_gene_table$logFC > 0),]
top_down <- top_gene_table[which(top_gene_table$adj.P.Val < 0.05 & top_gene_table$logFC < 0),]
# write.csv(top_gene_table, "/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/ANALYSIS/RSVvMOCK_top_gene_table.csv")
```
### Differential Expression - EdgeR
```{r}
dge <- estimateDisp(d, design, robust = T)

fit <- glmQLFit(dge, design)

fit <- glmQLFTest(fit)

res_edgeR=as.data.frame(topTags(fit, n=Inf))
head(res_edgeR)
res_edgeR$ensembl_gene_id_version <- row.names(res_edgeR)
head(res_edgeR)

conversions <- read.csv("/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/ENSG_symbol_conversions.csv")
res <- left_join(res_edgeR, conversions, by=c("ensembl_gene_id_version"))
# Differential Expressed Gene Table
res_edgeR <- res %>% dplyr::select(external_gene_name, logCPM, PValue, logFC, FDR, ensembl_gene_id_version)
```


**Explore DE Genes**
```{r}
# Get Gene Symbols to match Ensembl Gene IDs
ensembl <- as.data.frame(rownames(d))
colnames(ensembl) <- "ensembl_gene_id_version"
conversions <- read.csv("/storage1/fs1/leyao.wang/Active/RNA_SEQ_DATA_2024/ENSG_symbol_conversions.csv")
gene_names <- left_join(ensembl, conversions, by=c("ensembl_gene_id_version"))

# Add Gene Names to DGElist Object
d$genes <- data.frame(ENSEMBL=rownames(d), SYMBOL=gene_names$external_gene_name) 

# Create Interactive MD plot
glMDPlot(vfit, counts=dtmm, groups= d$samples$group, side.main = "SYMBOL", anno = d$genes, launch=FALSE, status = dt)

# View Differential Expressed Gene Counts at significance = 5%
summary(decideTests(vfit))

# Volcano Plot
volcanoplot(vfit, coef=2, highlight = 392, names = d$genes$SYMBOL)

# MD Plot
plotMD(vfit, status = d$samples$group, hl.pch = 16, hl.cex = 0.7, bg.cex =0.7, bg.pch = 21) 
```

**HEATMAP*
```{r}
i <- which(d$genes$ENSEMBL %in% top_up$rowname)
i <- which(d$genes$ENSEMBL %in% top_down$rowname)
library(gplots)
library(heatmap3)

dtmm <- as.matrix(dtmm)
heatmap.2(dtmm[i,], scale="row", labRow=d$genes$SYMBOL[i], labCol= colnames(dtmm), trace="none", density.info="none", dendrogram="column", cexRow = 0.8, cexCol = 1, Rowv = FALSE, srtCol = 45, col =redgreen(100))
```

**Average Expression DTMM - show expression across all samples of the 102 genes found to be significant in prior RNA-Seq Experiment**
```{r}
dtmm102 <- read.csv("/storage1/fs1/leyao.wang/Active/saran/RNA_Sara/finalRNAdata/EdgeR_Limma/COvMock102.csv")
dtmm102

dtmm_table <- dtmm %>% filter(row.names(dtmm) %in% dtmm102$rowname)
dtmm_table$Co <- apply(dtmm_table[,1:4], 1, mean)
dtmm_table$Lact <- apply(dtmm_table[,5:8], 1, mean)
dtmm_table$Mock <- apply(dtmm_table[,9:12], 1, mean)
dtmm_table$RSV <- apply(dtmm_table[,13:16], 1, mean)
dtmm_table <- dtmm_table %>% dplyr::select(Co, Mock, Lact, RSV)

DTMM <- left_join(tibble::rownames_to_column(dtmm_table), results, by=c("rowname" = "ensembl_gene_id_version"))

DTMM <- DTMM %>% dplyr::select(Mock, RSV, Lact, Co,external_gene_name)
DTMM[which(DTMM$external_gene_name == "" | is.na(DTMM$external_gene_name)), "external_gene_name"] <- DTMM[which(DTMM$external_gene_name == "" | is.na(DTMM$external_gene_name)), "rowname"] 
row.names(DTMM) <- DTMM$external_gene_name
DTMM <- DTMM %>% dplyr::select(Mock, RSV, Lact, Co)
DTMM <- DTMM %>% dplyr::select("Co1","Co2", "Co3", "Co4", "Lac1",  "Lac2", "Lac3", "Lac4",  "Mock1",  "Mock2",  "Mock3", "Mock4", "RSV1", "RSV2", "RSV3","RSV4" )

par(mar = c(5,5,8,10))
hm <- heatmap.2(as.matrix(DTMM), scale="row", labRow=row.names(DTMM), labCol=  colnames(DTMM), trace="none", density.info="none", dendrogram="column", cexRow = 0.5, cexCol = 1.5, Rowv = FALSE, srtCol = 0, col =bluered(100))
# dev.off()

hm

# Write Heatmap to PDF
pdf("AveExpr_Heatmap.pdf", width = 13,height = 0.4*102)
hm <- heatmap.2(as.matrix(DTMM), scale="row", labRow=row.names(DTMM), labCol=  colnames(DTMM), trace="none", density.info="none", dendrogram="column", cexRow = 0.5, cexCol = 1.5, Rowv = FALSE, srtCol = 0, col =bluered(100))
dev.off()
```

### Volcano Plot 
```{r}
# Results from Limma/EdgeR
data <- top
data <- left_join(tibble::rownames_to_column(top), results, by=c("rowname" = "ensembl_gene_id_version"))

diffexpr <- which(data$adj.P.Val <= 0.05)
data$FDR_Significance <- "Not Significant"

# Remove cases with NA Values
de <- data[complete.cases(data), ]
# if log2Foldchange > 1 and pvalue < 0.05, set as "UP" 
de$FDR_Significance[de$logFC > 0.5 & de$adj.P.Val < 0.05] <- "Significantly Up-Regulated"
# if log2Foldchange < -1 and pvalue < 0.05, set as "DOWN"
de$FDR_Significance[de$logFC < -0.5 & de$adj.P.Val < 0.05] <- "Significantly Down-Regulated"
# if log2Foldchange < -1 and pvalue < 0.05, set as "DOWN"
# Labels for Significant
de$delabel <- NA
de$delabel[de$FDR_Significance != "Not Significant"] <- de$external_gene_name[de$FDR_Significance != "Not Significant"]


library(ggrepel)
options(ggrepel.max.overlaps = Inf)

hm <- ggplot(data=de, aes(x=logFC, y=-log10(adj.P.Val), col=FDR_Significance, label=delabel)) +
        geom_point() + 
        theme_minimal() +
        geom_text_repel(size=3.5) +
        scale_color_manual(values=c("black", "blue",  "red", "green",  "purple")) +
        geom_vline(xintercept=c(-0.5, 0.5), col="blue", linetype=2) +
        geom_hline(yintercept=-log10(0.05), col="blue", linetype=2) + 
        labs(y = "-log10(Adjusted.Pvalue)", x = "Log2 Fold Change") + 
        theme(axis.title.x = element_text(size=16, face="bold"), axis.title.y = element_text( size=16, face="bold"), legend.text = element_text(size=), legend.title = element_text(size=12)) + xlim(-5, 12) #+ ylim(0,7)
hm
```


## Boxplots
```{r}
# Open pdf file
pdf(file= "./boxplots.pdf")
  
library(ggplot2)
d$norm <- dtmm

top_genes <- which(d$genes$ENSEMBL %in% top_gene_table$rowname)

# create a 2X2 grid
par( mfrow= c(11,2))
layout(matrix(1:4, nrow = 2))
par(mar = c(5,3,2,2))

for (i in top_genes) {
boxplot(unlist(d$norm[i,]) ~ d$samples$group, ylab = "Expression (TMM normalized)", xlab = d$genes$SYMBOL[i])
stripchart(unlist(d$norm[i,])~ d$samples$group,              # Data
           method = "jitter", # Random noise
           pch = 19,          # Pch symbols
           col = 4,           # Color of the symbol
           vertical = TRUE,   # Vertical mode
           add = TRUE)        # Add it over
}
```

# GSEA #
### initializations
```{r, message=FALSE}
library(dplyr)
library(tibble)
library(ggplot2)
library(pathview)
library(gage)
library(gageData)
library(annotate)
library(clusterProfiler)
library(DOSE)
library(enrichplot)
```

## Read in Contrasts
Read in the Results from performing Differential Expression Analysis with EdgeR/Limma. 
```{r}
res <- top_gene_table
res[,"ensembl_gene_id_version"] <- res[,"rowname"]

# remove unnecessary columns
res <- res %>% dplyr::select(ensembl_gene_id_version, logFC, AveExpr, t, P.Value, adj.P.Val, B)
head(res)
```

**Remove Version**
Must get ENSG IDs without version for AnnotationDbi to retrieve entrez IDs for GAGE.
```{r}
# Join tables & re-format
table <- left_join(res, results, by=c("ensembl_gene_id_version"))
res <- table %>% dplyr::select(ensembl_gene_id, everything())
head(res)
```

**Use AnnotationDbi to retrieve Entrez IDs and other annotations (optional).**
```{r}
library("AnnotationDbi")
library("org.Hs.eg.db")
columns(org.Hs.eg.db)

# Add Annotations
res$entrez = mapIds(org.Hs.eg.db,
                     keys=res$ensembl_gene_id, 
                     column="ENTREZID",
                     keytype="ENSEMBL",
                     multiVals="first")
res$symbol = mapIds(org.Hs.eg.db,
                     keys=res$ensembl_gene_id, 
                     column="SYMBOL",
                     keytype="ENSEMBL",
                     multiVals="first")
res$name =   mapIds(org.Hs.eg.db,
                     keys=res$ensembl_gene_id, 
                     column="GENENAME",
                     keytype="ENSEMBL",
                     multiVals="first")
res$go =   mapIds(org.Hs.eg.db,
                     keys=res$ensembl_gene_id, 
                     column="GO",
                     keytype="ENSEMBL",
                     multiVals="first")
```

# GAGE
```{r}
# Prepare Input for GAGE - all log fold change values and entrez IDs
foldchanges = res$logFC
names(foldchanges) = res$entrez

# GO - Gene Ontology Database
data(go.sets.hs)
data(go.subs.hs)

# GO - Biological Process
gobpsets = go.sets.hs[go.subs.hs$BP]
gobpres = gage(foldchanges, gsets=gobpsets, same.dir=TRUE)

# GO - Cellular Component
goCCsets = go.sets.hs[go.subs.hs$CC]
goCCres = gage(foldchanges, gsets=goCCsets, same.dir=TRUE)

# GO - Molecular Function
goMFsets = go.sets.hs[go.subs.hs$MF]
goMFres = gage(foldchanges, gsets=goMFsets, same.dir=TRUE)

## MSigDB - Molecular Signatures Database
msigdb <- readList("/storage1/fs1/leyao.wang/Active/saran/RNA_Sara/h.all.v2022.1.Hs.entrez.gmt") 
msigres = gage(foldchanges, gsets=msigdb, same.dir=TRUE)
```

# GO Pathway Analysis Pathways - Biological Process
```{r}
### Greater
neg.pval <-  -log(gobpres$greater[,4])
neg.p.sort <- sort(neg.pval, decreasing =  T)
head(neg.p.sort, 16)
par(mar=c(5,25,4,1)+.1)
barplot(neg.p.sort[1:10], col="red",horiz=TRUE, cex.names=0.1, las =2, main = "GO Biological Processes", xlab = "-log(adj.p.val)")

### Less
neg.pval <-  -log(gobpres$less[,4])
neg.p.sort <- sort(neg.pval, decreasing =  T)
head(neg.p.sort)
par(mar=c(5,18,4,1)+.1)
barplot(rev(neg.p.sort[1:10]), col="blue",horiz=TRUE, cex.names=0.8, las =2, main = "GO Biological Processes", xlab = "-log(adj.p.val)")

### Combined
df1 <- as.data.frame(gobpres$greater) %>% filter(q.val <= 0.05)
df2 <- as.data.frame(gobpres$less) %>% filter(q.val <= 0.05)
df <- rbind(df1,df2)
par(mar=c(5,22,4,1)+.1)
barplot(rev(-log(df$q.val[1:26])), names.arg = rownames(df), col="red", horiz=TRUE, cex.names=0.9, las =2, main = "GO Biological Processes", xlab = "-log(q.val)")

par(mar=c(5,22,4,1)+.1)
barplot(df$stat.mean[1:26], names.arg = rownames(df), col="red", horiz=TRUE, cex.names=0.9, las =2, main = "GO Biological Processes", xlab = "stat.mean")
```

# MSigDB Pathway Analysis Pathways
```{r}
### Greater
length(which(msigres$greater[,4] < 0.05))
neg.pval <-  -log(msigres$greater[,4])
neg.p.sort <- sort(neg.pval, decreasing =  T)
head(neg.p.sort)
par(mar=c(5,20,4,1)+.1)
barplot(neg.p.sort[1:2], col="red", horiz=TRUE, cex.names=0.8, las =2, main = "MSigDB Biological Processes", xlab = "-log(adj.p.val)")

### Less
length(which(msigres$less[,4] < 0.05))
neg.pval <-  -log(msigres$greater[,4])
neg.p.sort <- sort(neg.pval, decreasing =  T)
head(neg.p.sort)
par(mar=c(5,20,4,1)+.1)
barplot(neg.p.sort[1:2], col="blue", horiz=TRUE, cex.names=0.8, las =2, main = "MSigDB Biological Processes", xlab = "-log(adj.p.val)")

###  Combined
df1 <- as.data.frame(msigres$greater) %>% filter(q.val <= 0.05)
df2 <- as.data.frame(msigres$less) %>% filter(q.val <= 0.05)
df <- rbind(df1,df2)
par(mar=c(5,22,4,1)+.1)
barplot(df$stat.mean[1:7], names.arg = rownames(df), col="red", horiz=TRUE, cex.names=0.9, las =2, main = "MSigDB Biological Processes", xlab = "stat.mean")

par(mar=c(5,18,4,1)+.1)
barplot(-log(df$q.val[1:10]), names.arg = rownames(df), col="red", horiz=TRUE, cex.names=0.9, las =2, main = "MSigDB Biological Processes", xlab = "-log(q.val)")
```

# WGCNA #####################################################################################################################################

```{r, echo=FALSE, message=FALSE} 
# Initialization
library(ggplot2)
library(dplyr)
library(WGCNA)
options(stringsAsFactors = FALSE)
enableWGCNAThreads()
allowWGCNAThreads()
```

```{r}
# Read in log2 CPM values from Limma VoomWithQualityWeights
v <- read.csv("/storage1/fs1/leyao.wang/Active/saran/RNA2/WGCNA/wgcna_COMOCK.csv", header = T)
wcgna_v <- v %>% dplyr::select(Mock1,Mock2,Mock3,Co1,Co2,Co3) #,Lac1,Lac2,Lac3,RSV1,RSV2,RSV3
# Set Rownames to ENSG version IDs
rownames(wcgna_v) <- v$X
# mockrsv <- wcgna_v %>% select(Mock1,Mock2,Mock3,RSV1,RSV2,RSV3)
# lacmock <- wcgna_v %>% select(Mock1,Mock2,Mock3,Lac1,Lac2,Lac3)
# comock <- wcgna_v %>% select(Mock1,Mock2,Mock3,Co1,Co2,Co3)
# Transpose table
wcgna_v <- t(wcgna_v)
# mockrsv <- t(mockrsv)
# lacmock <- t(lacmock)
# comock <- t(comock)
```

```{r}
# Read in MetaData to specify Mock & RSV samples
setwd("/storage1/fs1/leyao.wang/Active/saran/RNA2/WGCNA/")
metaData <- read.csv('coldata.csv', header = TRUE, sep = ",")
metaData
datTraits <- metaData %>% dplyr::select(condition)
row.names(datTraits) <- metaData$sample
datTraits

datTraits <- metaData[c(1,2,3,4,5,6),] %>% dplyr::select(condition)
row.names(datTraits) <- metaData$sample[c(1,2,3,4,5,6)]
```

### cluster samples to detect outliers
```{r}
sampleTree = hclust(dist(wcgna_v), method = "average");

datTraits$condition <- as.numeric(as.factor(datTraits$condition))
# Color by Condition
traitColors = numbers2colors(datTraits, signed = FALSE)
# Plot the sample dendrogram and the colors underneath.
plotDendroAndColors(sampleTree, traitColors,
                    groupLabels = names(datTraits), 
                    main = "Sample dendrogram and trait heatmap")
```

```{r}
# Choose a set of soft-thresholding powers
powers = c(c(1:10), seq(from = 12, to=30, by=2))
# Call the network topology analysis function
sft = pickSoftThreshold(wcgna_v, powerVector = powers, networkType = "unsigned", RsquaredCut = 0.8, verbose = 5)

# Plot the results:
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit, unsigned R^2",type="n",
     main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.8,col="red")
```

```{r}
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers,col="red")
```

### Run WCGNA
```{r}
# This function performs automatic network construction and module detection on large expression datasets in a block-wise manner.
library(WGCNA)
bwnet <- blockwiseModules(wcgna_v,
  maxBlockSize = ncol(wcgna_v), # What size chunks (how many genes) the calculations should be run in
  networkType = "unsigned",
  TOMType = "unsigned", # topological overlap matrix
  power = 9, # soft threshold for network construction - 9 suggested for unsigned network
  corType = "pearson",
  numericLabels = F, # Let's use numbers instead of colors for module labels
  randomSeed = 1234, # there's some randomness associated with this calculation
  minModuleSize = 30,
  mergeCutHeight = 0.25, # threshold for merging modules
  # saveTOMs = T,
  # saveTOMFileBase = "wgcna_313_blockwise"
)

# Plot the dendrogram and the module colors underneath for block 1(all genes)
module_eigengenes <- bwnet$MEs
module_colors <- bwnet$colors
table(module_colors)
plotDendroAndColors(bwnet$dendrograms[[1]], module_colors[bwnet$blockGenes[[1]]], "Module Colors",
                     dendroLabels = FALSE, hang = 0.03,
                     addGuide = TRUE, guideHang = 0.05,
                     main = "Gene dendrogram and module colors")

# Print out a preview
head(module_eigengenes)
MEs0 = orderMEs(module_eigengenes) #Reorder given (eigen-)vectors such that similar ones (as measured by correlation) are next to each other.
```

```{r}
plotEigengeneNetworks(module_eigengenes, "Eigengene adjacency heatmap",
                      marDendro = c(3,3,2,4),
                      marHeatmap = c(3,4,2,2), plotDendrograms = T,
                      xLabelsAngle = 90)

```

```{r}
metaData <- metaData[c(1,2,3,10,11,12),]

metaData$condition <- as.factor(metaData$condition)
design <- model.matrix(~0+metaData$condition)
colnames(design) <- levels(metaData$condition) #get the group

nGenes = ncol(wcgna_v);
nSamples = nrow(wcgna_v);
moduleTraitCor = cor(module_eigengenes, design, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

sizeGrWindow(12,8)
# Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(1, 25, 4, 4));

# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = colnames(design),
yLabels = names(module_eigengenes),
ySymbols = names(module_eigengenes),
colorLabels = FALSE,
colors = blueWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.4,
zlim = c(-1,1),
main = paste("Module-trait relationships"))
```

```{r}
library(dplyr)
# Boxplot
row.names(metaData) <- metaData$sample
metadata <- metaData %>% dplyr::select(condition)
if(T){
mes_group <- merge(MEs0,metadata,by="row.names") 
library(gplots)
library(ggpubr)
library(grid)
library(gridExtra)
draw_ggboxplot <- function(data,Module="Module",group="group"){
  ggboxplot(data,x=group, y=Module,
            ylab = paste0(Module),
            xlab = group,
            fill = group,
            palette = "jco",
            add="jitter",
            legend = "") +stat_compare_means(method = "t.test")
}

colorNames <- names(MEs0)
pdf("Module-trait-relationshipCOMOCK_boxplot.pdf", width = 7.5,height = 1.6*ncol(MEs0))
p <- lapply(colorNames,function(x) {
  draw_ggboxplot(mes_group, Module = x, group = "condition")
})
do.call(grid.arrange,c(p,ncol=2)) 
dev.off()
}
```

## Find Genes in Module Red, Brown, Yellow
```{r}
gene_module_key <- tibble::enframe(bwnet$colors, name = "gene", value = "module") %>%
  # Let's add the `ME` part so its more clear what these numbers are and it matches elsewhere
  dplyr::mutate(module = paste0("ME", module))

meS_genes <- gene_module_key %>%
  dplyr::filter(module == "MEpurple")

# meB_genes <- gene_module_key %>%
#   dplyr::filter(module == "MEbrown")
# 
# meR_genes <- gene_module_key %>%
#   dplyr::filter(module == "MEred")
# head(meR_genes)
```

```{r}
# load biomart package
library(biomaRt)
library(dplyr)

allgenes <- rbind(meS_genes)
# Create vector of ensembl gene IDs
ENSGv <- as.vector(allgenes$gene)

# define biomart object
mart <- useMart(biomart = "ensembl", dataset = "hsapiens_gene_ensembl")
# listAttributes(mart)

# query biomart - retrieve corresponding gene names to ENSGv
results <- getBM(attributes = c("ensembl_gene_id", "ensembl_gene_id_version"),
               filters = "ensembl_gene_id_version", values = ENSGv,
               mart = mart)
head(results)

# Join tables - SALMON
table <- left_join(meS_genes, results, by=c("gene" = "ensembl_gene_id_version"))
meS_genes <- table %>% dplyr::select(ensembl_gene_id, everything())
head(meS_genes)
```

```{r}
library("AnnotationDbi")
library("org.Hs.eg.db")
columns(org.Hs.eg.db)

# Add Annotations - SALMON
meS_genes$symbol = mapIds(org.Hs.eg.db,
                     keys=meS_genes$ensembl_gene_id,
                     column="SYMBOL",
                     keytype="ENSEMBL",
                     multiVals="first")
meS_genes$entrez = mapIds(org.Hs.eg.db,
                     keys=meS_genes$ensembl_gene_id,
                     column="ENTREZID",
                     keytype="ENSEMBL",
                     multiVals="first")
meS_genes$name =   mapIds(org.Hs.eg.db,
                     keys=meS_genes$ensembl_gene_id,
                     column="GENENAME",
                     keytype="ENSEMBL",
                     multiVals="first")
dfPRP <- as.data.frame(apply(meS_genes,2,as.character))
write.csv(dfPRP, "/storage1/fs1/leyao.wang/Active/saran/RNA2/WGCNA/PurpleCOMOCK_ModuleGenes.csv")
```

```{r}
make_module_heatmap <- function(module_name,
                                expression_mat = wcgna_v,
                                metadata_df = metaData,
                                gene_module_key_df = gene_module_key,
                                module_eigengenes_df = module_eigengenes) {
  # Create a summary heatmap of a given module.
  #
  # Args:
  # module_name: a character indicating what module should be plotted, e.g. "ME19"
  # expression_mat: The full gene expression matrix. Default is `normalized_counts`.
  # metadata_df: a data frame with refinebio_accession_code and time_point
  #              as columns. Default is `metadata`.
  # gene_module_key: a data.frame indicating what genes are a part of what modules. Default is `gene_module_key`.
  # module_eigengenes: a sample x eigengene data.frame with samples as row names. Default is `module_eigengenes`.
  #
  # Returns:
  # A heatmap of expression matrix for a module's genes, with a barplot of the
  # eigengene expression for that module.

  # Set up the module eigengene with its refinebio_accession_code
  module_eigengene <- module_eigengenes_df %>%
    dplyr::select(all_of(module_name)) %>%
    tibble::rownames_to_column("sample")

  # Set up column annotation from metadata
  col_annot_df <- metadata_df %>%
    # Only select the treatment and sample ID columns
    dplyr::select(sample, condition) %>%
    # Add on the eigengene expression by joining with sample IDs
    dplyr::inner_join(module_eigengene, by = "sample") %>%
    # Arrange by patient and time point
    dplyr::arrange(condition,sample) %>%
    # Store sample
    tibble::column_to_rownames("sample")

  # Create the ComplexHeatmap column annotation object
  col_annot <- ComplexHeatmap::HeatmapAnnotation(
    # Supply treatment labels
    condition = col_annot_df$condition,
    # Add annotation barplot
    module_eigengene = ComplexHeatmap::anno_barplot(dplyr::select(col_annot_df, module_name)),
    # Pick colors for each experimental group in time_point
    col = list(condition = c("mock" = "#f1a340", "co" = "#998ec3"))
  )

  # Get a vector of the Ensembl gene IDs that correspond to this module
  module_genes <- gene_module_key_df %>%
    dplyr::filter(module == module_name) %>%
    dplyr::pull(gene)

  # Set up the gene expression data frame
  mod_mat <- expression_mat %>%
    t() %>%
    as.data.frame() %>%
    # Only keep genes from this module
    dplyr::filter(rownames(.) %in% module_genes) %>%
    # Order the samples to match col_annot_df
    dplyr::select(rownames(col_annot_df)) %>%
    # Data needs to be a matrix
    as.matrix()

  # Normalize the gene expression values
  mod_mat <- mod_mat %>%
    # Scale can work on matrices, but it does it by column so we will need to
    # transpose first
    t() %>%
    scale() %>%
    # And now we need to transpose back
    t()

  # Create a color function based on standardized scale
  color_func <- circlize::colorRamp2(
    c(-2, 0, 2),
    c("#67a9cf", "#f7f7f7", "#ef8a62")
  )

  # Plot on a heatmap
  heatmap <- ComplexHeatmap::Heatmap(mod_mat,
    name = module_name,
    # Supply color function
    col = color_func,
    # Supply column annotation
    bottom_annotation = col_annot,
    # We don't want to cluster samples
    cluster_columns = FALSE,
    # We don't need to show sample or gene labels
    show_row_names = FALSE,
    show_column_names = FALSE
  )

  # Return heatmap
  return(heatmap)
}

mod_S_heatmap <- make_module_heatmap(module_name = "MEpurple")
mod_S_heatmap
```

###### ClusterProfiler Enrichment Analysis
```{r}
library(clusterProfiler)
library(dplyr)
library(DOSE)
de <- res$entrez[1:102]
de <- dfPRP$entrez

ego <- enrichGO(de, OrgDb = org.Hs.eg.db, ont = "BP")

library(enrichplot)
barplot(ego, showCategory=10) 

mutate(ego, qscore = -log(p.adjust, base=10)) %>% 
    barplot(x="qscore")

library(DOSE)
ego <- enrichGO(gene          = de,
                OrgDb         = org.Hs.eg.db,
                ont           = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.05,
        readable      = TRUE)
dotplot(ego, showCategory=15)

########################################################################################
top_genes <- read.csv("/storage1/fs1/leyao.wang/Active/saran/RNA2/Limma-EDgeR/top_DEGs_COvMock_cpm10.csv")

# Join tables
table <- left_join(dfPRP, top_genes, by=c("gene" = "rowname"))
prp_genes <- table %>% dplyr::select(entrez, everything())
head(prp_genes)
  
  prp_up <- prp_genes[prp_genes$logFC > 0, ]
  prp_up <- prp_up[order(prp_up$logFC, decreasing = T),]
  
  prp_down <- prp_genes[prp_genes$logFC < 0, ]
  prp_down <- prp_down[order(prp_down$logFC, decreasing = T),]

# Join tables
table <- left_join(dfYEL, top_genes, by=c("gene" = "rowname"))
yel_genes <- table %>% dplyr::select(entrez, everything())
head(yel_genes)
  
  yel_up <- yel_genes[yel_genes$logFC > 0, ]
  yel_up <- yel_up[order(yel_up$logFC, decreasing = T),]
  
  yel_down <- yel_genes[yel_genes$logFC < 0, ]
  yel_down <- yel_down[order(yel_down$logFC, decreasing = T),]

# Join tables
table <- left_join(dfBR, top_genes, by=c("gene" = "rowname"))
brn_genes <- table %>% dplyr::select(entrez, everything())
head(brn_genes)
  
  brn_up <- brn_genes[brn_genes$logFC > 0, ]
  brn_up <- brn_up[order(brn_up$logFC, decreasing = T),]
  
  brn_down <- brn_genes[brn_genes$logFC < 0, ]
  brn_down <- brn_down[order(brn_down$logFC, decreasing = T),]
 
df_up <- rbind(brn_up, prp_up, yel_up)
###run go analysis
  formula_res <- compareCluster(
    entrez ~ module,
    data = prp_down, 
    fun = "enrichGO",
    OrgDb = org.Hs.eg.db,
    ont = "BP",  
    pAdjustMethod = "BH",
    pvalueCutoff = 0.05,
    qvalueCutoff = 0.05,
    readable = T
  )
  
  lineage1_ego <- clusterProfiler::simplify( 
    formula_res,
    cutoff= 0.05,
    by="p.adjust",
    select_fun=min
  )
  
dotplot(lineage1_ego, showCategory=15)


###########################
ego_up <- enrichGO(gene     = prp_up$entrez, # change to prp_up or yel_up
              OrgDb         = org.Hs.eg.db,
              ont           = "BP",
              pAdjustMethod = "BH",
              pvalueCutoff  = 0.05,
              qvalueCutoff  = 0.05,
              readable      = TRUE)
dotplot(ego_up, showCategory=15)

lineage1_ego <- clusterProfiler::simplify( 
  ego_up,
  cutoff= 0.05,
  by="p.adjust",
  select_fun=min
)

dotplot(lineage1_ego, showCategory=15)

```