2.2_networkAnalysis.Rmd

---
title: "RNA-seq Mouse Placenta"
author: "Ha T. H. Vu"
output: html_document
---
  
```{r setup, include=FALSE}
options(max.print = "75")
knitr::opts_chunk$set(
  echo = TRUE,
  collapse = TRUE,
  comment = "#>",
  fig.path = "Files/",
  fig.width = 8,
  prompt = FALSE,
  tidy = FALSE,
  message = FALSE,
  warning = TRUE
)
knitr::opts_knit$set(width = 75)
```

Documentation for the analysis of mouse placental RNA-seq data at e7.5, e8.5 and e9.5. <br>
All analyses in this step were done in R 4.1.1 unless stated otherwise.

## 1. STRINGdb network analysis:
To obtain networks from STRING, use https://string-db.org/, then upload the lists of interested genes. After obtaining the largest connected component, use GLay algorithm in Cytoscape to subcluster the networks and obtain centralities. <br>
After getting each subnetwork, we look at their functional analysis to determine the relevance of the networks with package `clusterProfiler`. The following function is to run the GO analysis. All GO analysis was carried our in R 4.0.2 (see specific package versions in Methods).
```{r, eval=F}
library("clusterProfiler")
library("org.Mm.eg.db")

go <- function(subnetworkGenes) {
  GO <- enrichGO(gene = subnetworkGenes, OrgDb=org.Mm.eg.db, ont = "BP", pvalueCutoff = 1, qvalueCutoff = 1, maxGSSize=1000, readable = T, keyType = "ENSEMBL")
  res <- as.data.frame(GO@result)
  res$Rank <- seq(1, nrow(res))
  res$GeneRatio <- sapply(strsplit(res$GeneRatio, "/"), function(x) as.numeric(x[1])/as.numeric(x[2]))
  res$BgRatio <- sapply(strsplit(res$BgRatio, "/"), function(x) as.numeric(x[1])/as.numeric(x[2]))
  res$Fold <- as.numeric(res$GeneRatio)/as.numeric(res$BgRatio)
  #res <- res[res$qvalue <= 0.05 & res$Fold >= 2 & res$Count >= 5,]
  
  return(res)
}
```
We now run GO analysis for every subnetwork.
```{r, eval=F}
for (time in c("e7.5", "e8.5", "e9.5")) {
  if (time == "e7.5") {
    for (i in 1) {
      geneList <- read.table(paste0("Files/STRING/", time, "/largestComponent/", time, "_", i, "/", time, "_", i, "_TSS.bed"), header = F, sep = "\t")
      goTerms <- go(geneList$V4)
      #save(goTerms, file = paste0("Files/STRING/", time, "/largestComponent/", time, "_", i, "/", time, "_", i, "BP_filtered.rda") 
    }
  } else if (time == "e8.5") {
    for (i in 1:3) {
      geneList <- read.table(paste0("Files/STRING/", time, "/largestComponent/", time, "_", i, "/", time, "_", i, "_TSS.bed"), header = F, sep = "\t")
      goTerms <- go(geneList$V4)
      #save(goTerms, file = paste0("Files/STRING/", time, "/largestComponent/", time, "_", i, "/", time, "_", i, "BP_filtered.rda") 
    }
  } else {
    for (i in 1:5) {
      geneList <- read.table(paste0("Files/STRING/", time, "/largestComponent/", time, "_", i, "/", time, "_", i, "_TSS.bed"), header = F, sep = "\t")
      goTerms <- go(geneList$V4)
      #save(goTerms, file = paste0("Files/STRING/", time, "/largestComponent/", time, "_", i, "/", time, "_", i, "BP_filtered.rda") 
    }
  }
}
```
We started out with one subnetwork from the largest connected component of e7.5, three of e8.5 and fave of e9.5. After looking at the GO terms, we can narrow down to focus on some particular networks (`e7.5_1_STRING`, `e8.5_1_STRING`, `e9.5_1_STRING`, `e9.5_2_STRING`, `e9.5_3_STRING`, and `e9.5_4_STRING`). We then get the hub genes,  defined to be nodes that have degree, betweenness and closeness centralities in 90th quantile in their subnetworks. The following code is to obtain hub genes:
```{r}
library("tidyverse")
t2g <- read.table("Files/t2g.txt", header = T, sep = "\t")

#function to analyze STRING networks
hubSFunc <- function(network) {
  network <- inner_join(network, nameMap[,c("queryItem", "preferredName")], by = c("name" = "preferredName"))
  network <- inner_join(network, t2g[,2:3], by = c("queryItem" = "ens_gene"))
  network <- distinct(network)
  network <- subset(network, network$queryItem %in% coordinates$V4)
  #network analysis
  topDegree <- subset(network, network$Degree > quantile(network$Degree, 0.9))
  topCloseness <- subset(network, network$ClosenessCentrality > quantile(network$ClosenessCentrality, 0.9)) 
  topBetweenness <- subset(network, network$BetweennessCentrality > quantile(network$BetweennessCentrality, 0.9)) 

  all <- intersect(intersect(topDegree$ext_gene, topCloseness$ext_gene), topBetweenness$ext_gene)
  return(all)
}
```
Now we get hub genes for each subnetwork of interest.
```{r}
#load networks
for (time in c("e7.5", "e8.5", "e9.5")) {
  coordinates <- read.table(paste0("Files/", time, "specific_TSS.bed"), header = F, sep = "\t")
  nameMap <- read.table(paste0("Files/STRING/", time, "/", time, "_string_mapping.tsv"), sep = "\t", quote = "", header = F)
  colnames(nameMap) <- c("queryIndex", "queryItem",	"stringId",	"preferredName",	"annotation")
  if (time == "e7.5") {
    for (i in 1) {
      network <- read.table(paste0("Files/STRING/", time, "/largestComponent/", time, "_", i, "/", time, "_", i, "_nodeTable.csv"), sep = ",", header = T)
      all <- hubSFunc(network)
      print(paste0("Number of genes in the ", time, "_", i, "_STRING network is ", length(unique(network$name))))
      print(paste0("Hub genes of the ", time, "_", i, "_STRING network are:"))
      print(all)
      #save(all, file=paste0("Files/STRING/", time, "/largestComponent/", time, "_", i, "/", time, "_", i, "_hubGenes.rda"))
    }
  } else if (time == "e8.5") {
    for (i in 1) {
      network <- read.table(paste0("Files/STRING/", time, "/largestComponent/", time, "_", i, "/", time, "_", i, "_nodeTable.csv"), sep = ",", header = T)
      all <- hubSFunc(network)
      print(paste0("Number of genes in the ", time, "_", i, "_STRING network is ", length(unique(network$name))))
      print(paste0("Hub genes of the ", time, "_", i, "_STRING network are:"))
      print(all)
      #save(all, file=paste0("Files/STRING/", time, "/largestComponent/", time, "_", i, "/", time, "_", i, "_hubGenes.rda"))
    }
  } else {
    for (i in 1:4) {
      network <- read.table(paste0("Files/STRING/", time, "/largestComponent/", time, "_", i, "/", time, "_", i, "_nodeTable.csv"), sep = ",", header = T)
      all <- hubSFunc(network)
      print(paste0("Number of genes in the ", time, "_", i, "_STRING network is ", length(unique(network$name))))
      print(paste0("Hub genes of the ", time, "_", i, "_STRING network are:"))
      print(all)
      #save(all, file=paste0("Files/STRING/", time, "/largestComponent/", time, "_", i, "/", time, "_", i, "_hubGenes.rda"))
    }
  }
}
```

## 2. GENIE3 networks
First, we use `GENIE3` to infer networks. At each timepoint, as inputs for GENIE3, timepoint-specific transcripts with average TPM at the timepoint >= 5 were aggregated to obtain gene counts with the R package `tximport` (see script `2.1_getAbundance.Rmd` for details how to get the transcript expression as inputs for gene abundance).
```{r, eval=F}
library("tximport")
library("GENIE3")
library("dplyr")
setwd("/work/LAS/geetu-lab/hhvu/project1_2/rna-seq/PlacentaRNA-seq")

#building transcript names
t2g <- read.table("Files/t2g.txt", header = T, sep = "\t")

for (time in c("e7.5", "e8.5", "e9.5")) {
  sample <- list.files(path="Files/GENIE3/kallisto_GENIE3", pattern = time)
  files <- file.path("Files/GENIE3/kallisto_GENIE3/", sample)
  kallisto.tsv <- tximport(files, type = "kallisto", tx2gene = t2g, ignoreAfterBar = TRUE, countsFromAbundance="lengthScaledTPM")
  exprMatr <- as.matrix(kallisto.tsv$counts)
  exprMatr <- t(scale(t(exprMatr)))
  regulators <- read.table(paste0("Files/", time, "specific_TFensGenes.txt"), header = F)
  regulators <- regulators[regulators$V1 %in% rownames(exprMatr), "V1"]
  set.seed(123) # For reproducibility of results
  weightMat <- GENIE3(exprMatr, regulators = as.vector(unique(regulators)))
  linkList <- getLinkList(weightMat)
  linkList2 <- getLinkList(weightMat, threshold=quantile(linkList$weight, 0.9))
  print(time)
  print(linkList2[1:3,])
  #write.table(linkList2, paste0("Files/GENIE3/", time, "/", time, "_0.9.txt"), row.names = F, quote = F, sep = "\t")
}
```
Similar to STRING network analysis, we carry out GO analysis first to evaluate the biological relevant of the networks. We can use the same function `go()`.
```{r, eval=F}
for (time in c("e7.5", "e8.5", "e9.5")) {
  for (i in 1:3) {
      geneList <- read.table(paste0("Files/GENIE3/", time, "/", time, "_", i, "/", time, "_", i, "_TSS.bed"), header = F, sep = "\t")
      goTerms <- go(geneList$V4)
      #save(goTerms, file = paste0("Files/GENIE3/", time, "/", time, "_", i, "/", time, "_", i, "BP_filtered.rda") 
  }
}
```
We identified `e7.5_2_GENIE3`, `e8.5_2_GENIE3`, `e9.5_1_GENIE3`, `e9.5_2_GENIE3` and `e9.5_3_GENIE3` as ones to focus on. Get hub genes for each subnetwork of interest.
```{r}
hubGFunc <- function(network) {
  network <- inner_join(network, t2g[,2:3], by = c("name" = "ens_gene"))
network <- distinct(network)
  topDegree <- subset(network, network$Degree > quantile(network$Degree, 0.9))
  topCloseness <- subset(network, network$ClosenessCentrality > quantile(network$ClosenessCentrality, 0.9)) 
  topBetweenness <- subset(network, network$BetweennessCentrality > quantile(network$BetweennessCentrality, 0.9)) 

  all <- intersect(intersect(topDegree$ext_gene, topCloseness$ext_gene), topBetweenness$ext_gene)
  
  return(all)
}
```

```{r}
library("dplyr")
t2g <- read.table("Files/t2g.txt", header = T, sep = "\t")
for (time in c("e7.5", "e8.5", "e9.5")) {
  coordinates <- read.table(paste0("Files/", time, "specific_TSS.bed"), header = F, sep = "\t")
  if (time == "e7.5") {
    for (i in 2) {
      network <- read.table(paste0("Files/GENIE3/", time, "/", time, "_", i, "/", time, "_", i, "_nodeTable.csv"), sep = ",", header = T)
      all <- hubGFunc(network)
      print(paste0("Number of genes in the ", time, "_", i, "_GENIE3 network is ", length(unique(network$name))))
      print(paste0("Hub genes of the ", time, "_", i, "_GENIE3 network are:"))
      print(all)
      #save(all, file=paste0("Files/GENIE3/", time, "/", time, "_", i, "/", time, "_", i, "_hubGenes.rda")
    }
  } else if (time == "e8.5") {
    for (i in 2) {
      network <- read.table(paste0("Files/GENIE3/", time, "/", time, "_", i, "/", time, "_", i, "_nodeTable.csv"), sep = ",", header = T)
      all <- hubGFunc(network)
      print(paste0("Number of genes in the ", time, "_", i, "_GENIE3 network is ", length(unique(network$name))))
      print(paste0("Hub genes of the ", time, "_", i, "_GENIE3 network are:"))
      print(all)
      #save(all, file=paste0("Files/GENIE3/", time, "/", time, "_", i, "/", time, "_", i, "_hubGenes.rda")
    }
  } else {
    for (i in 1:3) {
      network <- read.table(paste0("Files/GENIE3/", time, "/", time, "_", i, "/", time, "_", i, "_nodeTable.csv"), sep = ",", header = T)
      all <- hubGFunc(network)
      print(paste0("Number of genes in the ", time, "_", i, "_GENIE3 network is ", length(unique(network$name))))
      print(paste0("Hub genes of the ", time, "_", i, "_GENIE3 network are:"))
      print(all)
      #save(all, file=paste0("Files/GENIE3/", time, "/", time, "_", i, "/", time, "_", i, "_hubGenes.rda")
    }
  }
}
```

## 3. Hub gene analyses:
Hub genes that are known or implicated. Known = annotated as related to placenta/trophoblast/trophectoderm/chorion in GO biological process and MGI phenotype terms. Implicated = annotated as related to embryo in GO biological process and MGI phenotype terms.
```{r}
mart <- biomaRt::useMart("ensembl", host="http://sep2019.archive.ensembl.org", dataset = "mmusculus_gene_ensembl")
mm10 <- biomaRt::getBM(attributes = c("ensembl_gene_id", "external_gene_name", "mgi_id", "gene_biotype"), mart = mart)
placenta <- read.table("/work/LAS/geetu-lab/hhvu/project1_2/rna-seq/hubGeneCuration/MGI_GO_placentaChorioTBtrophec_genes.txt", header = F)
placenta <- dplyr::inner_join(placenta, mm10, by = c("V1" = "mgi_id"))
embryo <- read.table("/work/LAS/geetu-lab/hhvu/project1_2/rna-seq/hubGeneCuration/MGI_GO_embryo_genes.txt", header = F)
embryo <- dplyr::inner_join(embryo, mm10, by = c("V1" = "mgi_id"))
annot <- unique(c(placenta$external_gene_name, embryo$external_gene_name))
```
E7.5_1_STRING:
```{r}
load("Files/STRING/e7.5/largestComponent/e7.5_1/e7.5_1_hubGenes.rda")

length(intersect(all, placenta$external_gene_name))
length(setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name))


intersect(all, placenta$external_gene_name)
setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name)

intersect(all, placenta$external_gene_name)[order(intersect(all, placenta$external_gene_name))]
setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name)[order(setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name))]
setdiff(all, annot)[order(setdiff(all, annot))]
all[order(all)]
```
E7.5_2_GENIE3:
```{r}
load("Files/GENIE3/e7.5/e7.5_2/e7.5_2_hubGenes.rda")

length(intersect(all, placenta$external_gene_name))
length(setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name))

intersect(all, placenta$external_gene_name)
setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name)

intersect(all, placenta$external_gene_name)[order(intersect(all, placenta$external_gene_name))]
setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name)[order(setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name))]
setdiff(all, annot)[order(setdiff(all, annot))]
all[order(all)]
```
E8.5_1_STRING:
```{r}
load("Files/STRING/e8.5/largestComponent/e8.5_1/e8.5_1_hubGenes.rda")

length(intersect(all, placenta$external_gene_name))
length(setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name))

intersect(all, placenta$external_gene_name)
setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name)

intersect(all, placenta$external_gene_name)[order(intersect(all, placenta$external_gene_name))]
setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name)[order(setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name))]
setdiff(all, annot)[order(setdiff(all, annot))]
all[order(all)]
```
E8.5_2_GENIE3:
```{r}
load("Files/GENIE3/e8.5/e8.5_2/e8.5_2_hubGenes.rda")

length(intersect(all, placenta$external_gene_name))
length(setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name))

intersect(all, placenta$external_gene_name)
setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name)

intersect(all, placenta$external_gene_name)[order(intersect(all, placenta$external_gene_name))]
setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name)[order(setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name))]
setdiff(all, annot)[order(setdiff(all, annot))]
all[order(all)]
```
E9.5 networks:
```{r}
load("Files/STRING/e9.5/largestComponent/e9.5_1/e9.5_1_hubGenes.rda")
temp <- all
load("Files/STRING/e9.5/largestComponent/e9.5_2/e9.5_2_hubGenes.rda")
temp <- c(all, temp)
load("Files/STRING/e9.5/largestComponent/e9.5_3/e9.5_3_hubGenes.rda")
temp <- c(all, temp)
load("Files/STRING/e9.5/largestComponent/e9.5_4/e9.5_4_hubGenes.rda")
temp <- c(all, temp)
load("Files/GENIE3/e9.5/e9.5_1/e9.5_1_hubGenes.rda")
temp <- c(all, temp)
load("Files/GENIE3/e9.5/e9.5_2/e9.5_2_hubGenes.rda")
temp <- c(all, temp)
load("Files/GENIE3/e9.5/e9.5_3/e9.5_3_hubGenes.rda")
temp <- c(all, temp)
all <- temp

length(intersect(all, placenta$external_gene_name))
length(setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name))

intersect(all, placenta$external_gene_name)
setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name)

intersect(all, placenta$external_gene_name)[order(intersect(all, placenta$external_gene_name))]
setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name)[order(setdiff(intersect(all, embryo$external_gene_name), placenta$external_gene_name))]
setdiff(all, annot)[order(setdiff(all, annot))]
all[order(all)]
```

Test if the number of known/implicated genes in the networks is significant: <br>
E7.5 genes in networks
```{r}
t2g <- read.table("Files/t2g.txt", header = T)

allGenes <- read.table("Files/e7.5specific_ensGenes.txt", header = F)
allGenes <- unique(t2g[t2g$ens_gene %in% allGenes$V1, "ext_gene"])

load("Files/STRING/e7.5/largestComponent/e7.5_1/e7.5_1_hubGenes.rda")
temp <- all
load("Files/GENIE3/e7.5/e7.5_2/e7.5_2_hubGenes.rda")
all <- c(all, temp)

set.seed(1234)
t <- 0
t2 <- c()
for (i in 1:10000) {
  s <- sample(allGenes, 17)
  t2 <- c(t2, length(intersect(s, annot)))
  if (length(intersect(s, annot)) > length(intersect(all, annot))) {
    t <- t+1
  }
}
t/10000
t2<-as.data.frame(t2)

ggplot(t2, aes(x = t2)) + geom_histogram(binwidth = 1, color="darkgreen", fill="#b8e186") +
    geom_vline(aes(xintercept = length(intersect(all, annot))), color="black", linetype="dashed", size=1) +
  ggtitle("e7.5 hub genes")
```

E8.5 genes in networks
```{r}
allGenes <- read.table("Files/e8.5specific_ensGenes.txt", header = F)
allGenes <- unique(t2g[t2g$ens_gene %in% allGenes$V1, "ext_gene"])

load("Files/STRING/e8.5/largestComponent/e8.5_1/e8.5_1_hubGenes.rda")
temp <- all
load("Files/GENIE3/e8.5/e8.5_2/e8.5_2_hubGenes.rda")
all <- c(all, temp)

set.seed(1234)
t <- 0
t2 <- c()
for (i in 1:10000) {
  s <- sample(allGenes, 28)
  #print(length(intersect(s, placenta$external_gene_name)))
  t2 <- c(t2, length(intersect(s, annot)))
  if (length(intersect(s, annot)) >= length(intersect(all, annot))) {
    t <- t+1
  }
}
t/10000
t2<-as.data.frame(t2)

ggplot(t2, aes(x = t2)) + geom_histogram(binwidth = 1, color="darkblue", fill="lightblue") +
    geom_vline(aes(xintercept = length(intersect(all, annot))), color="black", linetype="dashed", size=1) +
  ggtitle("e8.5 hub genes")
```

E9.5 genes in networks
```{r}
allGenes <- read.table("Files/e9.5specific_ensGenes.txt", header = F)
allGenes <- unique(t2g[t2g$ens_gene %in% allGenes$V1, "ext_gene"])

load("Files/STRING/e9.5/largestComponent/e9.5_1/e9.5_1_hubGenes.rda")
temp <- all
load("Files/STRING/e9.5/largestComponent/e9.5_2/e9.5_2_hubGenes.rda")
temp <- c(all, temp)
load("Files/STRING/e9.5/largestComponent/e9.5_3/e9.5_3_hubGenes.rda")
temp <- c(all, temp)
load("Files/STRING/e9.5/largestComponent/e9.5_4/e9.5_4_hubGenes.rda")
temp <- c(all, temp)
load("Files/GENIE3/e9.5/e9.5_1/e9.5_1_hubGenes.rda")
temp <- c(all, temp)
load("Files/GENIE3/e9.5/e9.5_2/e9.5_2_hubGenes.rda")
temp <- c(all, temp)
load("Files/GENIE3/e9.5/e9.5_3/e9.5_3_hubGenes.rda")
temp <- c(all, temp)
all <- temp

set.seed(1234)
t <- 0
t2 <- c()
for (i in 1:10000) {
  s <- sample(allGenes, 127)
  #print(length(intersect(s, placenta$external_gene_name)))
  t2 <- c(t2, length(intersect(s, annot)))
  if (length(intersect(s, annot)) >= length(intersect(all, annot))) {
    t <- t+1
  }
}
t/10000
t2<-as.data.frame(t2)

ggplot(t2, aes(x = t2)) + geom_histogram(binwidth = 1, alpha=.2, fill="#FF6666", color="#841F27") +
    geom_vline(aes(xintercept = length(intersect(all, annot))), color="black", linetype="dashed", size=1) +
  ggtitle("e9.5 hub genes")
```

Hub genes in one method are also nodes in another method: <br>
E7.5
```{r}
e7.5_1_STRING <- read.table("Files/STRING/e7.5/largestComponent/e7.5_1/e7.5_1_nodeTable.csv", header = T, sep = ",")
e7.5_2_GENIE3 <- read.table("Files/GENIE3/e7.5/e7.5_2/e7.5_2_nodeTable.csv", header = T, sep = ",")
e7.5_2_GENIE3 <- inner_join(e7.5_2_GENIE3, mm10[,1:2], by = c("name" = "ensembl_gene_id"))
c <- intersect(e7.5_1_STRING$shared.name, e7.5_2_GENIE3$external_gene_name)
load("Files/STRING/e7.5/largestComponent/e7.5_1/e7.5_1_hubGenes.rda")
intersect(all, c)
load("Files/GENIE3/e7.5/e7.5_2/e7.5_2_hubGenes.rda")
intersect(all, c)
```
E8.5
```{r}
e8.5_1_STRING <- read.table("Files/STRING/e8.5/largestComponent/e8.5_1/e8.5_1_nodeTable.csv", header = T, sep = ",")
e8.5_2_GENIE3 <- read.table("Files/GENIE3/e8.5/e8.5_2/e8.5_2_nodeTable.csv", header = T, sep = ",")
e8.5_2_GENIE3 <- inner_join(e8.5_2_GENIE3, mm10[,1:2], by = c("name" = "ensembl_gene_id"))
c <- intersect(e8.5_1_STRING$shared.name, e8.5_2_GENIE3$external_gene_name)

load("Files/STRING/e8.5/largestComponent/e8.5_1/e8.5_1_hubGenes.rda")
intersect(all, c)

load("Files/GENIE3/e8.5/e8.5_2/e8.5_2_hubGenes.rda")
intersect(all, c)
```

E9.5
```{r}
e9.5_1_STRING <- read.table("Files/STRING/e9.5/largestComponent/e9.5_1/e9.5_1_nodeTable.csv", header = T, sep = ",")
e9.5_2_STRING <- read.table("Files/STRING/e9.5/largestComponent/e9.5_2/e9.5_2_nodeTable.csv", header = T, sep = ",")
e9.5_3_STRING <- read.table("Files/STRING/e9.5/largestComponent/e9.5_3/e9.5_3_nodeTable.csv", header = T, sep = ",")
e9.5_4_STRING <- read.table("Files/STRING/e9.5/largestComponent/e9.5_4/e9.5_4_nodeTable.csv", header = T, sep = ",")

e9.5_1_GENIE3 <- read.table("Files/GENIE3/e9.5/e9.5_1/e9.5_1_nodeTable.csv", header = T, sep = ",")
e9.5_1_GENIE3 <- inner_join(e9.5_1_GENIE3, mm10[,1:2], by = c("name" = "ensembl_gene_id"))
e9.5_2_GENIE3 <- read.table("Files/GENIE3/e9.5/e9.5_2/e9.5_2_nodeTable.csv", header = T, sep = ",")
e9.5_2_GENIE3 <- inner_join(e9.5_2_GENIE3, mm10[,1:2], by = c("name" = "ensembl_gene_id"))
e9.5_3_GENIE3 <- read.table("Files/GENIE3/e9.5/e9.5_3/e9.5_3_nodeTable.csv", header = T, sep = ",")
e9.5_3_GENIE3 <- inner_join(e9.5_3_GENIE3, mm10[,1:2], by = c("name" = "ensembl_gene_id"))

c <- intersect(c(e9.5_1_STRING$shared.name, e9.5_2_STRING$shared.name, e9.5_3_STRING$shared.name, e9.5_4_STRING$shared.name), c(e9.5_1_GENIE3$external_gene_name, e9.5_2_GENIE3$external_gene_name, e9.5_3_GENIE3$external_gene_name))
length(c)

load("Files/STRING/e9.5/largestComponent/e9.5_1/e9.5_1_hubGenes.rda")
intersect(all, c)
load("Files/STRING/e9.5/largestComponent/e9.5_2/e9.5_2_hubGenes.rda")
intersect(all, c)
load("Files/STRING/e9.5/largestComponent/e9.5_3/e9.5_3_hubGenes.rda")
intersect(all, c)
load("Files/STRING/e9.5/largestComponent/e9.5_4/e9.5_4_hubGenes.rda")
intersect(all, c)
load("Files/GENIE3/e9.5/e9.5_1/e9.5_1_hubGenes.rda")
intersect(all, c)
load("Files/GENIE3/e9.5/e9.5_2/e9.5_2_hubGenes.rda")
intersect(all, c)
load("Files/GENIE3/e9.5/e9.5_3/e9.5_3_hubGenes.rda")
intersect(all, c)
```

## 4. Randomization tests:
To test if the networks are significantly enriched for the GO terms of interest, compared to random gene sets of the same size. The following script is to run the randomization test. It will take a long time to run 10000 times, hence not run here.
```{r, eval=F}
#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)

library("doParallel")

go <- function(subnetworkGenes) {
  GO <- enrichGO(gene = subnetworkGenes, OrgDb=org.Mm.eg.db, ont = "BP", qvalueCutoff = 1, pvalueCutoff = 1, maxGSSize=1000, readable = T, keyType = "ENSEMBL")
  goSimplify <- as.data.frame(GO)

  return(goSimplify)
}

randFunc <- function(pop, size, go_term, i) {
  library("clusterProfiler", suppressMessages())
  library("org.Mm.eg.db", suppressMessages())

  gene_set <- sample(pop, size)
  res <- go(gene_set)
  if (go_term %in% res$Description) {
    p <- res[res$Description == go_term, "qvalue"]
  } else {
    p <- 1
  }

  print(i)
  return(p)
}

pop <- read.table(args[1], header = T)
sub <- read.table(args[2], sep = ",", header = T)
no_cores <- 20
cl <- makeCluster(no_cores)
registerDoParallel(cl)
result <- foreach(i=1:10000) %dopar% randFunc(unique(pop$ens_gene), length(sub$name), args[4], i)
stopCluster(cl)

result2 <- unlist(result)
write.table(result2, args[3], row.names = F, sep = "\t", quote = F)
```
Example command:
```{BASH, eval=F}
module load gcc/7.3.0-xegsmw4
module load r/4.0.2-py3-icvulwq
module load gsl/2.5-fpqcpxf
module load udunits/2.2.24-yldmp4h
module load gdal/2.4.4-nw2drgf
module load geos/3.8.1-2m7gav4

Rscript --vanilla randomizationGO.R "/work/LAS/geetu-lab/hhvu/project1_2/rna-seq/PlacentaRNA-seq/Files/e7.5Group.txt" "/work/LAS/geetu-lab/hhvu/project1_2/rna-seq/PlacentaRNA-seq/Files/STRING/e7.5/largestComponent/e7.5_1/e7.5_1_nodeTable.csv" "/work/LAS/geetu-lab/hhvu/project1_2/rna-seq/randomizationGO/e7.5_1_STRING_inflammatory response_qval.txt" "inflammatory response"
```
```{r}
sessionInfo()
```