---
title: "Set-level analysis of gene expression data with Gene ontology"
output:
  html_document:
    df_print: paged
---

Gene Ontology is the world’s largest knowledgebase of the functions of genes. In this notebook we will demonstrate its utilization for set-level analysis of gene expression data. The main goals will be:

* Access and understand gene-ontology annotation data, use MSigDB or Bioconductor packages.
* Learn and compare two ways of set-level gene expression analysis:
   + over-representation analysis (ORA),
   + gene-set enrichment analysis (GSEA).
* Visualize and interpret the reached results.

```{r prepare, echo=FALSE, message=FALSE, warning=FALSE, paged.print=FALSE}

# autism data, 
# see https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0049475
# GEO: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE18123

# DEG tutorials
# see https://rpubs.com/jrgonzalezISGlobal/enrichment
# see https://bioinformatics-core-shared-training.github.io/cruk-summer-school-2018/RNASeq2018/html/06_Gene_set_testing.nb.html

# deal with ExpressionSets
library(Biobase)
# differential expression for microarray data, 
# use limma instead of DESeq2 suitable for count data (such as RNA-seq)
library(limma)
# annotation library for the given microarray: Affymetrix Human Genome U133 Plus 2.0 Array
library(hgu133plus2.db)
# boxplots
library(ggplot2)
# enrichment analysis
library(clusterProfiler)
# fast gene-set enrichment analysis
library(fgsea)
# visualize cluterProfiler's results
library(enrichplot)
# pathway based data integration and visualization
library(pathview)
# deframe function
library(tibble) 
# classification with elastic nets
library(glmnet)
```

## The input data

[Kong et al.](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0049475 "Kong, Sek Won, et al. “Characteristics and predictive value of blood transcriptome signature in males with autism spectrum disorders.” PloS one 7.12 (2012): e49475.") presents the results of a large blood transcriptome study that aims to evaluate the utility of gene expression profiling in the diagnosis of Autism Spectrum Disorders (ASD). The full gene expression dataset is available through Gene Expression Omnibus. In this notebook, we will work with a RData expression file that contains 54,613 features (detection probes) and 99 samples (individuals). 66 individuals suffer with an ASD, 33 individuals represent the control healthy group. This file corresponds to the train dataset referenced in the paper.

There are four classes in the expression set: CONTROL, AUTISM, PDD-NOS, ASPERGER'S DISORDER. We will merge the three non-control groups into a single category referred to as Autism spectrum disorder (ASD).

```{r load, echo=FALSE, message=FALSE, warning=FALSE, paged.print=FALSE}
load("GSE18123.Rdata")
colnames(phenoData(gse18123))
phenoData(gse18123)$group
head(exprs(gse18123))

# merge three autism disorders?
merge_disorders=T

if (merge_disorders){
  experiments_of_interest <- c("ASD", "CONTROL")
  levels(phenoData(gse18123)$group)<-c("CONTROL","AUTISM","AUTISM","AUTISM")
  levels(phenoData(gse18123)$group)<-c("CONTROL","ASD")
} else experiments_of_interest <- c("AUTISM", "CONTROL")

columns_of_interest <- which(phenoData(gse18123)$group %in% experiments_of_interest)
grouping <- droplevels(phenoData(gse18123)$group[columns_of_interest])
counts_of_interest <- exprs(gse18123)[,columns_of_interest]
```

## Differential gene expression

The first step is to quantify differential expression for the individual probes. We deal with microarray data in this task, we will use limma library. Essentially, linear regression is applied.

```{r diffexp}

design <- model.matrix(~ grouping)

# deals with log2 scaled data 
fit <- limma::lmFit(log2(counts_of_interest), design)
result <- limma::eBayes(fit)
topTable(result)

# plot the top results, their ordering not kept
top12<-rownames(topTable(result,number=12))
df<-data.frame(probe=rep(top12,length(columns_of_interest)),class=rep(grouping,each=12),data=as.vector(exprs(gse18123)[rownames(exprs(gse18123)) %in% top12,columns_of_interest]))
ggplot(df, aes(x=probe, y=data, fill=class)) + 
  geom_boxplot() +
  facet_wrap(~probe, scale="free")
```

## Over-representation analysis

Over-representation analysis determines whether genes from a pre-defined set (in here, those beloging to a specific GO term) are present more than would be expected (over-represented) in a subset of your data (in here, a set of differentially expressed genes). In fact, it is a repetitive application of hypergeometric test.

```{r ora}
# identify differentially expressed genes first, adjusted p.value must be smaller than 0.05, 412 genes selected (around 0.7%)
DEGs<-topTable(result,number=Inf,p.value=0.05)

# get their annotation
dbq<-AnnotationDbi::select(hgu133plus2.db, keys=rownames(gse18123), keytype="PROBEID",columns = c("SYMBOL","ENTREZID"),multiVals="list")
annot<-dbq[match(unique(dbq$PROBEID), dbq$PROBEID),]
DEGs_ENTREZID<-annot$ENTREZID[annot$PROBEID %in% rownames(DEGs)]
DEGs_ENTREZID<-DEGs_ENTREZID[!is.na(DEGs_ENTREZID)]

# use clusterProfiler with built-in gene lists
ans.go <- enrichGO(gene = DEGs_ENTREZID, ont = "BP",
                   OrgDb ="org.Hs.eg.db",
                   universe = annot$ENTREZID,
                   readable=TRUE,
                   pvalueCutoff = 1, qvalueCutoff=1)

tab.go <- as.data.frame(ans.go)
tab.go<- subset(tab.go, Count>5)
tab.go[1:5, 1:6]

# show multiple scores at once, adjusted p-value as well as count, sort by adjusted p-value
# count is the number of DE genes that belong to a given gene-set
barplot(ans.go, showCategory=10)

# visualizes gene sets as a network, mutually overlapping gene sets tend to cluster together, making it easier for interpretation.
#emapplot(ans.go)
emapplot(pairwise_termsim(ans.go))

ans.kegg <- enrichKEGG(gene = DEGs_ENTREZID,
                       organism = 'hsa',
                       universe = annot$ENTREZID,
                       pvalueCutoff = 0.05, qvalueCutoff=1)
tab.kegg <- as.data.frame(ans.kegg)
tab.kegg<- subset(tab.kegg, Count>5)
tab.kegg[1:5, 1:6]

# show multiple scores at once, GeneRatio, adjusted p-value and count
# GeneRatio=count/setSize, where 'count' is the number of DE genes that belong to a given gene-set, while 'setSize' is the total number of genes in the gene-set
# pushes up larger gene sets, such as Endocytosis (contains 249 genes)
dotplot(ans.kegg, showCategory=12) + ggtitle("KEGG")

# visualize particular pathway
logFC <- DEGs$logFC
names(logFC) <- DEGs_ENTREZID
pathview(gene.data = logFC, 
         pathway.id = "hsa04140", 
         species = "hsa", 
         limit = list(gene=2, cpd=1))

# deal with external gene lists
c3.tf <- read.gmt("c3.tft.v7.5.1.entrez.gmt")
# by default, abnalyze only the sets with minimum size 10 and maximum size 500
ans.tf <- enricher(DEGs_ENTREZID, TERM2GENE=c3.tf,universe = annot$ENTREZID)
#ans.tf <- enricher(DEGs_ENTREZID, TERM2GENE=c3.tf,universe = annot$ENTREZID, pvalueCutoff = 1, qvalueCutoff = 1) # full result to compare with fgsea
tab.tf <- as.data.frame(ans.tf)
tab.tf.large<- subset(tab.tf, Count>5)
tab.tf.large[1:5,1:5]

# visualizes gene sets as a network, mutually overlapping gene sets tend to cluster together, making it easier for interpretation.
#emapplot(ans.tf)
emapplot(pairwise_termsim(ans.tf))

# run Fisher test for the NFMUE1_Q6 list
# inspect the list first
head(c3.tf)
NFMUE1_Q6<-c3.tf$gene[c3.tf$term=="NFMUE1_Q6"] # older call: NFMUE1_Q6<-c3.tf$gene[c3.tf$ont=="NFMUE1_Q6"]
annot[annot$ENTREZID %in% NFMUE1_Q6,] # 723 probes
# construct the contingency table and run the hypergeometric test
f_NFMUE1_Q6<-factor(unique(annot$ENTREZID) %in% NFMUE1_Q6)
f_DEGs_ENTREZID<-factor(unique(annot$ENTREZID) %in% DEGs_ENTREZID)
table(f_DEGs_ENTREZID,f_NFMUE1_Q6)
fisher.test(f_DEGs_ENTREZID,f_NFMUE1_Q6)
```


# Gene-set enrichment analysis

Over-representation analysis may suffer from a couple of limitations. First, after correcting for multiple hypotheses testing, no individual gene may meet the threshold for statistical significance. The relevant biological differences could be modest relative to the noise inherent to the expression measurement. Second, the selection of a set of DEGs leads to information loss as the quantitative information about gene expression is binarized. Again, subtle differential expression patterns spread accross entire groups of mutually linked genes could be omitted. GSEA ([Subramanian et al.](https://www.pnas.org/doi/10.1073/pnas.0506580102 " Subramanian, Aravind, et al. “Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles.” Proceedings of the National Academy of Sciences 102.43 (2005): 15545-15550.")) avoids these bottlenecks. 

```{r gsea}
probe.tstats<-topTable(result,number=Inf)$t
names(probe.tstats)<-annot$ENTREZID[match(rownames(topTable(result,number=Inf)),annot$PROBEID)]
probe.tstats<-probe.tstats[!is.na(names(probe.tstats))]
gene.tstats<-deframe(aggregate(probe.tstats,list(names(probe.tstats)),function(x) x[which(abs(x)==max(abs(x)))]))
gene.tstats<-deframe(aggregate(probe.tstats,list(names(probe.tstats)),mean))
c3.tf.lists <- gmtPathways("c3.tft.v7.5.1.entrez.gmt")
fgseaRes <- fgsea(pathways=c3.tf.lists, stats=gene.tstats, minSize=10, maxSize=500,nperm=100000)
View(fgseaRes[order(pval)])

# low p-values found again, however low ability to distinguish between gene sets
fgseaRes[fgseaRes$pathway %in% tab.tf[1:5,]$ID,]

# compare p-val ranking
cor.test(fgseaRes$pval,tab.tf$pvalue[match(fgseaRes$pathway,tab.tf$ID)],method="spearman")

convgo<-function(go_name){
  terms<-strsplit(go_name," ")[[1]]
  return(toupper(paste(c("GOBP",terms),collapse="_")))
}

c5.go.bp.lists <- gmtPathways("c5.go.bp.v7.5.1.entrez.gmt")
fgseaRes <- fgsea(pathways=c5.go.bp.lists, stats=gene.tstats, minSize=10, maxSize=500,nperm=100000)
View(fgseaRes[order(pval)])
# low p-values found again, however low ability to distinguish between gene sets
fgseaRes[fgseaRes$pathway %in% sapply(tab.go[1:5,]$Description,convgo)]
# compare p-val ranking
cor.test(fgseaRes$pval,tab.go$pvalue[match(fgseaRes$pathway,sapply(tab.go$Description,convgo))],method="spearman")
```

## Your task for this tutorial
  
See the chunk in the Rmd script attached below. Inspect it, run it and explain whether:
  
* the peripheral blood expression profiling for ASD detection/classification may be feasible, 
* the outlined gene set aggregation helps in ASD detection, 
* is the outlined evaluation of gene set aggregation unbiased, 
* the proposed GO BP terms match the terms proposed earlier in terms of ORA and GSEA and whether they should match them in theory. 

Eventually, propose your own way of gene-set-based feature extraction for ASD detection. 

```{r classify}

# classify ASD and control samples
# start with probes

# learn an elastic network, obtain its cross-validated performance
cvfit <- cv.glmnet(t(counts_of_interest),grouping,family="binomial",type.measure="auc",nfolds=5)
plot(cvfit,main="GLMNET: ASD classification with probes")
# get selected probes, i.e., the probes used for classification
tmp_coeffs <- coef(cvfit, s = "lambda.min")
lasso_features<-data.frame(name = tmp_coeffs@Dimnames[[1]][tmp_coeffs@i+1], coefficient = tmp_coeffs@x)[-1,]
annot[annot$PROBEID %in% lasso_features$name,]

# use SVD to identify the first principal component in each GS
svdGS<-function(TrainingDataSet){
  decomp<-svd(TrainingDataSet) 
  u<-decomp$u  
  inv_u<-tryCatch(leftinverse(u),error=function(err){        
    inv_u<-t(u)*NA   
  })
  dv<-diag(decomp$d)%*%t(decomp$v)
  return(dv[1,])
}

# precompute aggregate expressions for GO terms
c5.go.bp.lengths<-sapply(c5.go.bp.lists,length)
gs_counts<-t(sapply(c5.go.bp.lists[c5.go.bp.lengths>10 & c5.go.bp.lengths<100],function(x) svdGS(counts_of_interest[annot$ENTREZID %in% x,])))
# learn an elastic network to classify with them
cvfit <- cv.glmnet(t(gs_counts),grouping,family="binomial",type.measure="auc",nfolds=5)
plot(cvfit,main="GLMNET: ASD classification with GO terms")
tmp_coeffs <- coef(cvfit, s = "lambda.1se")
lasso_features<-data.frame(name = tmp_coeffs@Dimnames[[1]][tmp_coeffs@i+1], coefficient = tmp_coeffs@x)[-1,]
lasso_features
tab.go[sapply(tab.go$Description,convgo) %in% as.character(lasso_features$name),] # no big overlap with ORA
fgseaRes[fgseaRes$pathway %in% as.character(lasso_features$name)] # neither with GSEA
```
