set_analysis.Rmd

---
title: "Characterising newly made gene sets"
author: "The DEE2 Gene Signatures Group"
date: "`r Sys.Date()`"
output:
  html_document
---

Source: https://github.com/markziemann/dee2_gene_signatures

## Background

In this analysis we will be characterising the new gene sets.

## Libraries

```{r, libs}

suppressPackageStartupMessages({
  library("getDEE2")
  library("mitch")
  library("triwise")
  library("dplyr")
  library("gplots")
  library("reshape2")
  library("network")
  library("eulerr")
})

```

## Reference gene list

This is the background - all genes in the "universe" according to Ensembl version 90, which is the same annotation set used by DEE2. 
You can see that there are a large number of non-protein coding genes.

```{r,genelist}

# get universe of gene names and biotypes
if (! file.exists("hs.gtf.gz") ) {
  download.file("ftp://ftp.ensembl.org/pub/release-90/gtf/homo_sapiens/Homo_sapiens.GRCh38.90.gtf.gz",destfile= "hs.gtf.gz")
}

g<-read.table("hs.gtf.gz",sep="\t")
g <- g[grep("gene",g$V3),9] 
universe <- sapply(strsplit(g," "),"[[",6)
universe <- gsub(";","",universe)

biotypes <- sapply(strsplit(g," "),"[[",10)
biotypes <- gsub("_"," ",biotypes)
biotypes <- gsub(";","",biotypes)

biotypes_df <- data.frame(universe,biotypes)

universe <- unique(universe)

mytable <- table(biotypes)
mytable <- mytable[order(mytable)]
mytable <- mytable[which(mytable>100)]

par(mar=c(3,13,1,1)); barplot(mytable,horiz=TRUE,las=2,cex.names = 0.6,cex.axis = 0.6, 
  main="number of genes in each biotype class",xlim = c(0,20000)) ;grid()

mytable <- mytable/sum(mytable)*100

par(mar=c(3,13,1,1)); barplot(mytable,horiz=TRUE,las=2,cex.names = 0.6,cex.axis = 0.6, 
  main="proportion (%) of genes in each biotype class" ,xlim=c(0,40)) ;grid()

# Take note how many transcripts are protein coding versus non-protein coding
paste("protein coding:",length(which(biotypes=="protein coding")))
paste("non-protein coding:",length(which(biotypes!="protein coding")))

```

## Curated gene sets

### Reactome

Here I download the current Reactome gene set library.
I would like to see how many Ensembl genes have some sort of annotated function.
Breaking it down into protein coding and non-protein coding, we can see that non-protein coding genes are severely underrepresented in Reactome.
Reactome has 2400 sets and a total of 11193 genes as of 29/Sep/2020.

```{r,reactome}

#o
download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip", destfile="ReactomePathways.gmt.zip")
#unzip("ReactomePathways.gmt.zip")
reactome <- gmt_import("ReactomePathways.gmt")
reactome_genes <- unique(unname(unlist(reactome)))
length(reactome)
length(reactome_genes)
v1 <- list("Reactome"=reactome_genes, "Ensembl universe"=universe)
library("eulerr")
plot(euler(v1),quantities = TRUE)
length(reactome_genes)/length(universe)*100

prot <- biotypes_df[which(biotypes_df$biotypes=="protein coding"),1] 
prot <- unique(prot)
nprot <- biotypes_df[which(biotypes_df$biotypes!="protein coding"),1]
nprot <- unique(nprot)
intersect(prot,nprot)
nprot <-   setdiff(nprot,prot)
v1 <- list("Reactome"=reactome_genes, "protein coding"=prot, "non-protein coding"=nprot)
plot(euler(v1),quantities = TRUE)

```

### Gene Ontology

Now let's do the same with GO sets downloaded from Ensembl biomart.
GO has 12264 sets and 8412 genes as of 26/Sep/2020.

```{r,go}

go <- read.table("biomart_2020-09-26.txt.gz",sep="\t", fill = TRUE)
colnames(go) <- go[1,]
go <- go[2:nrow(go),]
go <- go[which(go$`GO term accession`!=""),]
go_genes <- unique(go$`Gene name`)
length(unique(go$`GO term accession`))
length(go_genes)
v1 <- list("GO"=go_genes, "protein coding"=prot, "non-protein coding"=nprot)
plot(euler(v1),quantities = TRUE)
length(unique(go$`GO term accession`))
head(go)

```

### Human Phenotype Ontology

Now same for Human phenotype ontology.
Downloaded from https://hpo.jax.org/app/download/annotation.
There are 2958 sets and 4366 genes as of 27/Sep/2020.

```{r,hpo}

hpo <- read.table("phenotype_to_genes.txt",sep="\t")
hpo_genes <- unique(hpo$V4)
head(hpo)
length(unique(hpo$V1))
length(hpo_genes)
v1 <- list("HPO"=hpo_genes, "protein coding"=prot, "non-protein coding"=nprot)
plot(euler(v1),quantities = TRUE)

```

### MSigDB

Now for MSigDB version 7.2 accessed 27/Sep/2020.
There are 31120 sets and 40044 genes.

```{r,msigdb}
msigdb <- gmt_import("msigdb.v7.2.symbols.gmt")
msigdb_genes <- unique(unname(unlist(msigdb)))
length(msigdb)
length(msigdb_genes)
v1 <- list("MSigDB"=msigdb_genes,"protein coding"=prot, "non-protein coding"=nprot)
plot(euler(v1),quantities = TRUE)

```

## Defining the overrepresentation analysis function

Here I'm using the triwise::testEnrichment function which uses a Fisher test under the hood.

Reference
van de Laar L, Saelens W, De Prijck S, Martens L, Scott CL, Van Isterdael G, Hoffmann E, Beyaert R, Saeys Y, Lambrecht BN, Guilliams M. Yolk Sac Macrophages, Fetal Liver, and Adult Monocytes Can Colonize an Empty Niche and Develop into Functional Tissue-Resident Macrophages. Immunity. 2016 Apr 19;44(4):755-68. doi: 10.1016/j.immuni.2016.02.017. Epub 2016 Mar 15. PMID: 26992565.

```{r,orafunc}

ora <- function(gs,genesets,universe){ 
  res <- triwise::testEnrichment(gs, genesets, universe, minknown = 3, mindiffexp = 2, maxknown = 5000)
  if( nrow(res)>0 ) {
    res <- res[order(res$pval),]
    res <- head(res,5)
    hits <- log(res[which(res$qval<0.01),2])
    names(hits) <- res[which(res$qval<0.01),4]
    return(hits)
  } else {
    return(0)
  }
}

```

Here is a function that converts list of vectors to a network diagram.
There are 2x edges than nodes.
Only the edges with highest similarity are retained, as per jaccard.
The size of the gene set is proportional to the node size (sqrt).

```{r,gsnet}

gs2net <- function(gset){
  mydf <- bind_rows(lapply(gset, as.data.frame.list))
  rownames(mydf) <- names(gset)
  j <- apply(mydf,1,function(x) {
    apply(mydf,1,function(y) {
      length(intersect(x,y) ) / length(union(x,y))
    })
  })
  j[lower.tri(j)] <- NA
  j[lower.tri(j,diag=TRUE)] <- 0
  jl <- melt(j)
  jl <- jl[which(jl$Var1 != jl$Var2),]
  jl <- jl[which(jl$value != 1),]
  jl <- jl[order(-jl$value),]
  jl <- head(jl,length(gset)*2)
  jl$edgeSize = with(jl, jl$value/sum(jl$value))
  lengths <- unlist(lapply(gset,length))
  lengths <- sqrt(lengths/sum(lengths)*100)
  jl$vertexsize <- lengths[match(as.character(jl$Var1),names(lengths))]
  jlnet = network(jl[,1:2])
  plot(jlnet, displaylabels = TRUE, label.col = "steelblue",
       edge.lwd = c(jl$edgeSize) * 100,
       arrowhead.cex = 0,
       label.cex = 1, vertex.border = "white",vertex.cex = jl$vertexsize,
       vertex.col = "blue", edge.col = rgb(0, 0, 0, alpha = 0.5))
}

```


## Epilepsy

Let's analyse the new epilepsy gene sets.
There are only 12 sets.
We begin with a histogram of the set size.

```{r,epilepsy1}

epi <- gmt_import("epilepsy_genesymbols.gmt")
epi <- diab[which(!names(epi)=='')]
epi <- epi[which(!duplicated(names(epi)))]
length(epi)
setnames <- names(epi)
names(epi) <-paste("GS",1:length(epi))
numgenes <- unlist(lapply(epi,length))
hist(numgenes,breaks = 15, xlab="number of genes per set",main="Epilepsy gene sets")

```

Next we can look at the similarity between the gene sets with a network diagram.
GS8, GS12 and GS10 have some overlap, as do GS7 with GS11.

Next we can see that the gene sets have diverse representation of biotypes. 
There is a high proportion of protein coding genes but also lincRNA.

```{r,epilepsy2,fig.width=8,fig.height=6}

gs2net(epi)

# classes
mytable <- lapply(epi,function(x) {
  table(biotypes_df[which(x %in% biotypes_df$universe),2])
})
mydf <- bind_rows(lapply(mytable, as.data.frame.list))
rownames(mydf) <- names(epi)
mydf_n <- mydf/rowSums(mydf)
mydf <- mydf[,head(order(-colSums(mydf_n)),8)]
colfunc <- colorRampPalette(c("white", "yellow","orange","red","darkred"))(n=25)
heatmap.2(t(as.matrix(mydf)),col=colfunc,trace="none",margin=c(5,10),cexCol = 0.9, cexRow = 0.9,
  scale="col",dendrogram = "none",main="Epilepsy gene sets by biotype",na.color = "white")

```

The overall biotype has a high representaiton of protein coding genes similar to Reactome.
We also observed a higher representation of lincRNAs in these gene sets.
Pseudogenes were underrepresented.

```{r,epilepsy3}

mydf_n[is.na(mydf_n)] <- 0
mytable <-colSums(mydf_n)
mytable <- mytable/sum(mytable)*100
mytable<-mytable[order(mytable)]
par(mar=c(3,13,1,1)); barplot(mytable,horiz=TRUE,las=2,cex.names = 0.7,cex.axis = 0.8, 
  main="Overall biotype class representation (%)" ) ;grid()

```

Next we wanted to know the Reactome functional categories that are enriched in each of the gene sets.
Using a Fisher text, we performed an enrichment analysis with Reactome genes using Ensembl as the background (universe).
We are aware that it may the improper background and the best background would be all genes that were detected in the original experiment.
As expected, some of these sets are related to neurological functions.

```{r,epilepsy4,fig.width=10,fig.height=8}

# ORA
res <- lapply(epi, function(x) { ora(gs=x,genesets=reactome, universe=universe)  }  )
res <- res[which(lapply(res,length)>0)]
myres <- bind_rows(lapply(res, as.data.frame.list))
myres[is.na(myres)] <- 0
rownames(myres) <- names(res)
colfunc <- colorRampPalette(c("white", "lightblue","blue", "darkblue"))(n=25)
heatmap.2(t(as.matrix(myres)),col=colfunc,trace="none",margin=c(10,25),cexRow = 0.6,
  dendrogram = "none",main="Epilepsy gene sets by Reactome enrichment")

```


## Diabetes

Now lets analyse the new diabetes gene sets.
There are 129 gene sets.

```{r,diabetes1}

diab <- gmt_import("diabetes_genesymbols.gmt")
diab <- diab[which(!names(diab)=='')]
diab <- diab[which(!duplicated(names(diab)))]
length(diab)
setnames <- names(diab)
names(diab) <-paste("GS",1:length(diab))
numgenes <- unlist(lapply(diab,length))
hist(numgenes,breaks = 15, xlab="number of genes per set",main="Diabetes gene sets")

```

The network diagram indicates clustering into a few groups. 
A large group with GS39 at the centre (right side).
A smaller cluster around GS121 (left side).
A smaller cluster around GS120 (top).

The biotype heatmap indicates representation of protein coding, lincRNAs and processed pseudgenes.
In some sets, non-coding genes outnumber protein coding genes (left).
In most sets, protein coding genes are more abundant (right).

```{r,diabetes2,fig.width=8,fig.height=6}

gs2net(diab)

# classes
mytable <- lapply(diab,function(x) {
  table(biotypes_df[which(x %in% biotypes_df$universe),2])
})
mydf <- bind_rows(lapply(mytable, as.data.frame.list))
rownames(mydf) <- names(diab)
mydf_n <- mydf/rowSums(mydf)
mydf <- mydf[,head(order(-colSums(mydf_n)),8)]
colfunc <- colorRampPalette(c("white", "yellow","orange","red","darkred"))(n=25)
heatmap.2(t(as.matrix(mydf)),col=colfunc,trace="none",margin=c(5,10),cexCol = 0.3, 
  cexRow = 0.6 ,scale="col",dendrogram = "none",main="Diabetes gene sets by biotype",na.color = "white")

```

In these gene sets, psudogenes are underrepresented while protein coding, lincRNA and antisense RNA genes are overrepresented.

```{r,diabetes3}

mydf_n[is.na(mydf_n)] <- 0
mytable <-colSums(mydf_n)
mytable <- mytable/sum(mytable)*100
mytable<-mytable[order(mytable)]
par(mar=c(3,13,1,1)); barplot(mytable,horiz=TRUE,las=2,cex.names = 0.7,cex.axis = 0.8, 
  main="Overall biotype class representation (%)" ) ;grid()

```

Here is the reactome enrichment plot for the diabetes gene sets. 
The text is too small to make out from this chart so I will need to look at making a PDF version.
One thing we can see from the chart is that not many of the enriched Reactomes are common apart from the horizontal strip near the bottom.

```{r,diabetes4,fig.width=10,fig.height=8}

# reactome enrichment
res <- lapply(diab, function(x) { ora(gs=x,genesets=reactome, universe=universe)  }  )
res <- res[which(lapply(res,length)>0)]
myres <- bind_rows(lapply(res, as.data.frame.list))
myres[is.na(myres)] <- 0
rownames(myres) <- names(res)
colfunc <- colorRampPalette(c("white", "lightblue","blue", "darkblue"))(n=25)
heatmap.2(t(as.matrix(myres)),col=colfunc,trace="none",margin=c(5,10),cexRow = 0.2, 
  cexCol = 0.3 ,dendrogram = "none",main="Diabetes gene sets by Reactome enrichment")

```


## Heart disease

Let's analyse the new heart disease gene sets.
There are 50 gene sets.

```{r,heartdisease1}

hd <- gmt_import("heartdisease_genesymbols.gmt")
hd <- diab[which(!names(hd)=='')]
hd <- hd[which(!duplicated(names(hd)))]
length(hd)
setnames <- names(hd)
names(hd) <-paste("GS",1:length(hd))
numgenes <- unlist(lapply(hd,length))
hist(numgenes,breaks = 15, xlab="number of genes per set",main="Heart disease gene sets")

```

The network chart indicatesmost sets are linked around GS50. 
There are several very small groups on the right of  the chart.

The biotype heatmap shows overall genes sets are composed motly of protein coding genes and a small number of sets have a majority of lincRNA and pseudo genes.

```{r,heartdisease2,fig.width=8,fig.height=6}

gs2net(hd)

# classes
mytable <- lapply(hd,function(x) {
  table(biotypes_df[which(x %in% biotypes_df$universe),2])
})
mydf <- bind_rows(lapply(mytable, as.data.frame.list))
rownames(mydf) <- names(hd)
mydf_n <- mydf/rowSums(mydf)
mydf <- mydf[,head(order(-colSums(mydf_n)),8)]
colfunc <- colorRampPalette(c("white", "yellow","orange","red","darkred"))(n=25)
heatmap.2(t(as.matrix(mydf)),col=colfunc,trace="none",margin=c(5,10),cexCol = 0.2,
  cexRow = 0.8,scale="col",dendrogram = "none",main="Heart disease gene sets by biotype",na.color = "white")

```

Similar to the other sets, pseudogenes are underrepresented and protein coding, lincRNA, and antisense genes are overrepresented.

```{r,heartdisease3}

mydf_n[is.na(mydf_n)] <- 0
mytable <-colSums(mydf_n)
mytable <- mytable/sum(mytable)*100
mytable<-mytable[order(mytable)]
par(mar=c(3,13,1,1)); barplot(mytable,horiz=TRUE,las=2,cex.names = 0.7,cex.axis = 0.8, 
  main="Overall biotype class representation (%)" ) ;grid()

```

The reactome heatmap shows a cluster gene sets that are enriched for ECM, collagen and IGF signaling (lower).
There is a strip showing metabolism, immune signaling and signal transduction are common.
Apart from that, the enriched gene sets are quite distinct.

```{r,heartdisease4,fig.width=10,fig.height=8}

# reactome enrichment
res <- lapply(hd, function(x) { ora(gs=x,genesets=reactome, universe=universe)  }  )
res <- res[which(lapply(res,length)>0)]
myres <- bind_rows(lapply(res, as.data.frame.list))
myres[is.na(myres)] <- 0
rownames(myres) <- names(res)
colfunc <- colorRampPalette(c("white", "lightblue","blue", "darkblue"))(n=25)
heatmap.2(t(as.matrix(myres)),col=colfunc,trace="none",margin=c(5,15),cexRow = 0.4, cexCol = 0.5, 
  dendrogram = "none",main="Heart disease gene sets by Reactome enrichment")

```

## SARS, MERS and SARS-CoV-2

Let's analyse the new virus infection gene sets.
There are 66 sets.

```{r,sars1}

sars <- gmt_import("sarsmers_genesymbols.gmt")
sars <- diab[which(!names(sars)=='')]
sars <- sars[which(!duplicated(names(sars)))]
length(sars)
setnames <- names(sars)
names(sars) <-paste("GS",1:length(sars))
numgenes <- unlist(lapply(sars,length))
hist(numgenes,breaks = 15, xlab="number of genes per set",main="SARS, MERS and SARS-CoV-2 gene sets")

```

The network diagram shows a large cluster of sets around GS21 (left side).
A smaller cluster exists around GS24 (right side) and GSGS29 (bottom).
The biotype chart shows most gene sets are dominated with protein coding genes but about 1/3 have dominance of lincRNA and pseudogenes.


```{r,sars2,fig.width=8,fig.height=6}

gs2net(sars)

# classes
mytable <- lapply(sars,function(x) {
  table(biotypes_df[which(x %in% biotypes_df$universe),2])
})
mydf <- bind_rows(lapply(mytable, as.data.frame.list))
rownames(mydf) <- names(sars)
mydf_n <- mydf/rowSums(mydf)
mydf <- mydf[,head(order(-colSums(mydf_n)),8)]
colfunc <- colorRampPalette(c("white", "yellow","orange","red","darkred"))(n=25)
heatmap.2(t(as.matrix(mydf)),col=colfunc,trace="none",margin=c(5,10),cexCol = 0.4, cexRow = 0.8,
  scale="col",dendrogram = "none",main="SARS, MERS and SARS-CoV2 gene sets by biotype",na.color = "white")

```

The overall biotype representation is similar to the other sets with overrepresentation of lincRNA and antisense RNA genes.

```{r,sars3}

mydf_n[is.na(mydf_n)] <- 0
mytable <-colSums(mydf_n)
mytable <- mytable/sum(mytable)*100
mytable<-mytable[order(mytable)]
par(mar=c(3,13,1,1)); barplot(mytable,horiz=TRUE,las=2,cex.names = 0.7,cex.axis = 0.8, 
  main="Overall biotype class representation (%)" ) ;grid()

```

There is a cluster of immune related enrichment in several gene sets as expected (lower left).
Moreover there are some clustering of integrins (lower left), and cell cycle (middle).
Post translational modification and signal transduction were also common.

```{r,sars4,fig.width=10,fig.height=8}

# Reactome enrichment
res <- lapply(sars, function(x) { ora(gs=x,genesets=reactome, universe=universe)  }  )
res <- res[which(lapply(res,length)>0)]
myres <- bind_rows(lapply(res, as.data.frame.list))
myres[is.na(myres)] <- 0
rownames(myres) <- names(res)
colfunc <- colorRampPalette(c("white", "lightblue","blue", "darkblue"))(n=25)
heatmap.2(t(as.matrix(myres)),col=colfunc,trace="none",margin=c(5,10),cexRow = 0.3, cexCol = 0.4,
  dendrogram = "none",main="SARS, MERS and SARS-CoV-2 gene sets genes by Reactome enrichment")

```

## Overall analysis

There are 257 gene sets in total.
The median set size was smallest for epilepsy and largest for diabetes.
In terms of number of genes in each library, epilepsy has 2175, diabetes has 11179, heart disease has 8727 and virus has 6827.
There were 731 genes that were common to all three libraries.
Altogether these new sets include 15522 genes, of these 12824 are protein coding and 2698 are non-protein coding.


```{r,overall1}

names(epi) <- paste("Epilepsy",names(epi) )
names(diab) <- paste("Diabetes",names(diab) )
names(hd) <- paste("Heart disease",names(hd) )
names(sars) <- paste("Virus related",names(sars) )
all <- c(epi,diab,hd,sars)
length(all)

barplot( c(length(epi) , length(diab) , length(hd) , length(sars) ) , 
         names.arg = c("Epilepsy","Diabetes","Heart Disease","Virus") ,
         ylab= "number of gene sets")

epil <- sapply(epi,length)
diabl <- sapply(diab,length)
hdl <- sapply(hd,length)
sarsl <- sapply(sars,length)

boxplot(epil,diabl,hdl,sarsl,names=c("Epilepsy","Diabetes","Heart Disease","Virus"))

epig <- unique(unname(unlist(epi)))
diabg <- unique(unname(unlist(diab)))
hdg <- unique(unname(unlist(hd)))
sarsg <- unique(unname(unlist(sars)))
length(epig)
length(diabg)
length(hdg)
length(sarsg)

v1 <- list("Epilepsy"=epig, "Diabetes"=diabg, "Heart diabetes"=hdg, "Virus related" = sarsg)
plot(euler(v1),quantities = TRUE)

allg <- unique(unname(unlist(all)))
length(allg)
v1 <- list("DEE2 sets"=allg,"protein coding"=prot, "non-protein coding"=nprot)
plot(euler(v1),quantities = TRUE)

```


## Session information

```{r sessioninfo}
sessionInfo()
```