RELIEF: overview analysis

Introduction

Here we will generate some graphs:

• To show the number of increased /decreased genes / probes in each contrast

• To show the number of increased /decreased reactome genesets in each contrast

• Other ideas include Euler/Venn and UpSet diagrams

Packages

suppressPackageStartupMessages({
    library("gplots")
    library("mitch")
    library("kableExtra")
    library("eulerr")
    library("UpSetR")    
    library(ggplot2)
})

Data import

There are 3 contrasts we are focusing on:

• post-op samples of patients with low CRP compared to high CRP

• baseline samples of patients with low CRP compared to high CRP

• compare pre-op samples to post-op samples

rna_podcrp <- read.table("pod_crp_rna.tsv",header=TRUE,sep="\t")

rna_t0crp <- read.table("t0_crp_rna.tsv",header=TRUE,sep="\t")

rna_t0_v_pod <- read.table("t0_v_pod_rna.tsv",header=TRUE,sep="\t")

meth_podcrp <- read.table("dm11.tsv",header=TRUE,sep="\t")

meth_t0crp <- read.table("dm12.tsv",header=TRUE,sep="\t")

meth_t0_v_pod <- read.table("dm10.tsv",header=TRUE,sep="\t")

Number of genes and probes

Here the methylation data is represented by probes, not summarised to gene names.

rna_podcrp_up <- rownames(subset(rna_podcrp, padj < 0.05 & log2FoldChange > 0 ))
rna_podcrp_dn <- rownames(subset(rna_podcrp, padj < 0.05 & log2FoldChange < 0 ))

rna_t0crp_up <- rownames(subset(rna_t0crp, padj < 0.05 & log2FoldChange > 0 ))
rna_t0crp_dn <- rownames(subset(rna_t0crp, padj < 0.05 & log2FoldChange <0 ))

rna_t0_v_pod_up <- rownames(subset(rna_t0_v_pod, padj < 0.05 & log2FoldChange > 0 ))
rna_t0_v_pod_dn <- rownames(subset(rna_t0_v_pod, padj < 0.05 & log2FoldChange < 0 ))

meth_podcrp_up <- rownames(subset(meth_podcrp, adj.P.Val < 0.05 & logFC > 0 ))
meth_podcrp_dn <- rownames(subset(meth_podcrp, adj.P.Val < 0.05 & logFC < 0 ))

meth_t0crp_up <- rownames(subset(meth_t0crp, adj.P.Val < 0.05 & logFC > 0 ))
meth_t0crp_dn <- rownames(subset(meth_t0crp, adj.P.Val < 0.05 & logFC < 0 ))

meth_t0_v_pod_up <- rownames(subset(meth_t0_v_pod, adj.P.Val < 0.05 & logFC > 0 ))
meth_t0_v_pod_dn <- rownames(subset(meth_t0_v_pod, adj.P.Val < 0.05 & logFC < 0 ))


xl <- list("RNA POD CRP up"=rna_podcrp_up,
  "RNA POD CRP down"=rna_podcrp_dn,
  "RNA t0 CRP up"=rna_t0crp_up,
  "RNA t0 CRP down"=rna_t0crp_dn,
  "RNA t0 Vs POD up"=rna_t0_v_pod_up,
  "RNA t0 Vs POD down"=rna_t0_v_pod_dn,
  "Meth POD CRP up"=meth_podcrp_up,
  "Meth POD CRP down"=meth_podcrp_dn,
  "Meth t0 CRP up"=meth_t0crp_up,
  "Meth t0 CRP down"=meth_t0crp_dn,
  "Meth t0 Vs POD up"=meth_t0_v_pod_up,
  "Meth t0 Vs POD down"=meth_t0_v_pod_dn)

xlr <- sapply(xl,length)

xlr_up <- xlr[grep("up",names(xlr))]

xlr_dn <- xlr[grep("down",names(xlr))]

MAX=max(xlr_up)*1.2
MIN=-max(xlr_dn)*1.2

par(mar=c(5,5,5,5))
names(xlr_up) <- gsub(" up","",names(xlr_up))
barplot(xlr_up,ylim=c(MIN,MAX),cex.names=0.8,ylab="no. genes/probes up- and down-regulated",col="lightblue")
text(x = (1:length(xlr_up)*1.2)-0.5  , y = xlr_up , label = xlr_up, pos = 3, cex = 1, col = "black")
barplot(-xlr_dn,ylim=c(MIN,MAX),cex.names=0.5, col="pink",add=T, xaxt='n')
text(x = (1:length(xlr_up)*1.2)-0.5  , y = -xlr_dn-5000 , label = xlr_dn, pos = 3, cex = 1, col = "black")

MAX=max(xlr_up)*1.4
MIN=-max(xlr_dn)*1.4

par(mar=c(5,15,5,5))
barplot(xlr_up,xlim=c(MIN,MAX),cex.names=1,main="no. genes/probes up- and down-regulated",col="lightblue",horiz=TRUE,las=1)
text(y = (1:length(xlr_up)*1.2)-0.75  , x = xlr_up + 5000 , label = xlr_up, pos = 3, cex = 1, col = "black")
barplot(-xlr_dn,ylim=c(MIN,MAX),cex.names=0.5, col="pink",add=T, yaxt='n',horiz=TRUE,las=1)
text(y = (1:length(xlr_up)*1.2)-0.75  , x = -xlr_dn-5000 , label = xlr_dn, pos = 3, cex = 1, col = "black")

Number of genes

Here the methylation data is summarised to gene names. RNA-seq data is also summarised to gene names.

rna_podcrp_up <- rownames(subset(rna_podcrp, padj < 0.05 & log2FoldChange > 0 ))
if (length(rna_podcrp_up)>0) { rna_podcrp_up <- unique(sapply(strsplit(rna_podcrp_up," "),"[[",2)) }

rna_podcrp_dn <- rownames(subset(rna_podcrp, padj < 0.05 & log2FoldChange < 0 ))
if (length(rna_podcrp_dn)>0) { rna_podcrp_dn <- unique(sapply(strsplit(rna_podcrp_dn," "),"[[",2)) }

rna_t0crp_up <- rownames(subset(rna_t0crp, padj < 0.05 & log2FoldChange > 0 ))
if (length(rna_t0crp_up)>0) { rna_t0crp_up <- unique(sapply(strsplit(rna_t0crp_up," "),"[[",2)) }

rna_t0crp_dn <- rownames(subset(rna_t0crp, padj < 0.05 & log2FoldChange <0 ))
if (length(rna_t0crp_dn)>0) { rna_t0crp_dn <- unique(sapply(strsplit(rna_t0crp_dn," "),"[[",2)) }

rna_t0_v_pod_up <- rownames(subset(rna_t0_v_pod, padj < 0.05 & log2FoldChange > 0 ))
if (length(rna_t0_v_pod_up)>0) { rna_t0_v_pod_up <- unique(sapply(strsplit(rna_t0_v_pod_up," "),"[[",2)) }

rna_t0_v_pod_dn <- rownames(subset(rna_t0_v_pod, padj < 0.05 & log2FoldChange < 0 ))
if (length(rna_t0_v_pod_dn)>0) { rna_t0_v_pod_dn <- unique(sapply(strsplit(rna_t0_v_pod_dn," "),"[[",2)) }

g <- subset(meth_podcrp, adj.P.Val < 0.05 & logFC > 0 )$UCSC_RefGene_Name
g<-g[which(g!="")]
g<-sapply(strsplit(g,";"),"[[",1)
meth_podcrp_up <- unique(g)

g <- subset(meth_podcrp, adj.P.Val < 0.05 & logFC < 0 )$UCSC_RefGene_Name
g<-g[which(g!="")]
g<-sapply(strsplit(g,";"),"[[",1)
meth_podcrp_dn <- unique(g)

g <- meth_podcrp$UCSC_RefGene_Name
g<-g[which(g!="")]
g<-sapply(strsplit(g,";"),"[[",1)
meth_podcrp_all <- unique(g)

g <- subset(meth_t0crp, adj.P.Val < 0.05 & logFC > 0 )$UCSC_RefGene_Name
g<-g[which(g!="")]
g<-sapply(strsplit(g,";"),"[[",1)
meth_t0crp_up <- unique(g)

g <- subset(meth_t0crp, adj.P.Val < 0.05 & logFC < 0 )$UCSC_RefGene_Name
g<-g[which(g!="")]
g<-sapply(strsplit(g,";"),"[[",1)
meth_t0crp_dn <- unique(g)

g <- subset(meth_t0_v_pod, adj.P.Val < 0.05 & logFC > 0 )$UCSC_RefGene_Name
g<-g[which(g!="")]
g<-sapply(strsplit(g,";"),"[[",1)
meth_t0_v_pod_up <- unique(g)

g <- subset(meth_t0_v_pod, adj.P.Val < 0.05 & logFC < 0 )$UCSC_RefGene_Name
g<-g[which(g!="")]
g<-sapply(strsplit(g,";"),"[[",1)
meth_t0_v_pod_dn <- unique(g)

g <- meth_t0_v_pod$UCSC_RefGene_Name
g<-g[which(g!="")]
g<-sapply(strsplit(g,";"),"[[",1)
meth_t0_v_pod_all <- unique(g)

xl <- list("RNA POD CRP up"=rna_podcrp_up,
  "RNA POD CRP down"=rna_podcrp_dn,
  "RNA t0 CRP up"=rna_t0crp_up,
  "RNA t0 CRP down"=rna_t0crp_dn,
  "RNA t0 Vs POD up"=rna_t0_v_pod_up,
  "RNA t0 Vs POD down"=rna_t0_v_pod_dn,
  "Meth POD CRP up"=meth_podcrp_up,
  "Meth POD CRP down"=meth_podcrp_dn,
  "Meth t0 CRP up"=meth_t0crp_up,
  "Meth t0 CRP down"=meth_t0crp_dn,
  "Meth t0 Vs POD up"=meth_t0_v_pod_up,
  "Meth t0 Vs POD down"=meth_t0_v_pod_dn)

xlr <- sapply(xl,length)

xlr_up <- xlr[grep("up",names(xlr))]

xlr_dn <- xlr[grep("down",names(xlr))]

MAX=max(xlr_up)*1.2
MIN=-max(xlr_dn)*1.2

par(mar=c(5,5,5,5))
names(xlr_up) <- gsub(" up","",names(xlr_up))
barplot(xlr_up,ylim=c(MIN,MAX),cex.names=0.8,ylab="no. genes/probes up- and down-regulated",col="lightblue")
text(x = (1:length(xlr_up)*1.2)-0.5  , y = xlr_up , label = xlr_up, pos = 3, cex = 1, col = "black")
barplot(-xlr_dn,ylim=c(MIN,MAX),cex.names=0.5, col="pink",add=T, xaxt='n')
text(x = (1:length(xlr_up)*1.2)-0.5  , y = -xlr_dn-1000 , label = xlr_dn, pos = 3, cex = 1, col = "black")

MAX=max(xlr_up)*1.4
MIN=-max(xlr_dn)*1.4

par(mar=c(5,15,5,5))
barplot(xlr_up,xlim=c(MIN,MAX),cex.names=1,main="no. genes/probes up- and down-regulated",col="lightblue",horiz=TRUE,las=1)
text(y = (1:length(xlr_up)*1.2)-0.75  , x = xlr_up + 1000 , label = xlr_up, pos = 3, cex = 1, col = "black")
barplot(-xlr_dn,ylim=c(MIN,MAX),cex.names=0.5, col="pink",add=T, yaxt='n',horiz=TRUE,las=1)
text(y = (1:length(xlr_up)*1.2)-0.75  , x = -xlr_dn-1000 , label = xlr_dn, pos = 3, cex = 1, col = "black")

Venn diagrams

Now to identify the gene overlaps (same omics different contrasts).

v1 <- list("RNA POD CRP up"=rna_podcrp_up,
  "RNA POD CRP down"=rna_podcrp_dn,
  "RNA t0 CRP up"=rna_t0crp_up,
  "RNA t0 CRP down"=rna_t0crp_dn,
  "RNA t0 Vs POD up"=rna_t0_v_pod_up,
  "RNA t0 Vs POD down"=rna_t0_v_pod_dn)

plot(euler(v1),quantities = TRUE)

## Warning in colSums(id & !empty) == 0 | merged_sets: longer object length is not
## a multiple of shorter object length

v1 <- list( "Meth POD CRP up"=meth_podcrp_up,
  "Meth POD CRP down"=meth_podcrp_dn,
  "Meth t0 CRP up"=meth_t0crp_up,
  "Meth t0 CRP down"=meth_t0crp_dn,
  "Meth t0 Vs POD up"=meth_t0_v_pod_up,
  "Meth t0 Vs POD down"=meth_t0_v_pod_dn)

plot(euler(v1),quantities = TRUE)

## Warning in colSums(id & !empty) == 0 | merged_sets: longer object length is not
## a multiple of shorter object length

Now to identify the overlaps (same contrast different omics).

v1 <- list("RNA POD CRP up"=rna_podcrp_up,
  "RNA POD CRP down"=rna_podcrp_dn,
  "Meth POD CRP up"=meth_podcrp_up,
  "Meth POD CRP down"=meth_podcrp_dn)

plot(euler(v1),quantities = TRUE)

v1 <-  list( "RNA t0 CRP up"=rna_t0crp_up,
  "RNA t0 CRP down"=rna_t0crp_dn,
  "Meth t0 CRP up"=meth_t0crp_up,
  "Meth t0 CRP down"=meth_t0crp_dn)

plot(euler(v1),quantities = TRUE)

v1 <- list("RNA t0 Vs POD up"=rna_t0_v_pod_up,
  "RNA t0 Vs POD down"=rna_t0_v_pod_dn,
  "Meth t0 Vs POD up"=meth_t0_v_pod_up,
  "Meth t0 Vs POD down"=meth_t0_v_pod_dn)

plot(euler(v1),quantities = TRUE)

Overrepresentaton test

Using the fisher test based on this tutorial https://seqqc.wordpress.com/2019/07/25/how-to-use-phyper-in-r/

Firstly for baseline vs post-op.

• over representation meth_t0_v_pod_dn and rna_t0_v_pod_up OR=2.1 p=2.2e-16

• under representation meth_t0_v_pod_up and rna_t0_v_pod_up OR=1.3 p=1

• over representation meth_t0_v_pod_up and rna_t0_v_pod_dn OR=1.2 p=2.5e-6

• under representation meth_t0_v_pod_dn and rna_t0_v_pod_dn OR=0.81 p-6.6e-11

# over representation meth_t0_v_pod_dn and rna_t0_v_pod_up
s1 = meth_t0_v_pod_dn
s2 = rna_t0_v_pod_up
tot = meth_t0_v_pod_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1,s2)) ,
  length(s2) - length(intersect(s1,s2)),
  length(tot) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='greater')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(tot) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is greater than 1
## 95 percent confidence interval:
##  1.987566      Inf
## sample estimates:
## odds ratio 
##   2.101154

# under representation meth_t0_v_pod_up and rna_t0_v_pod_up
s1 = meth_t0_v_pod_up
s2 = rna_t0_v_pod_up
tot = meth_t0_v_pod_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1 , s2)) ,
  length(s2) - length(intersect(s1 , s2)),
  length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='less')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value = 1
## alternative hypothesis: true odds ratio is less than 1
## 95 percent confidence interval:
##  0.000000 1.381368
## sample estimates:
## odds ratio 
##   1.292886

# over representation meth_t0_v_pod_up and rna_t0_v_pod_dn
s1 = meth_t0_v_pod_up
s2 = rna_t0_v_pod_dn
tot = meth_t0_v_pod_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1,s2)) ,
  length(s2) - length(intersect(s1,s2)),
  length(tot) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='greater')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(tot) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value = 2.491e-06
## alternative hypothesis: true odds ratio is greater than 1
## 95 percent confidence interval:
##  1.122921      Inf
## sample estimates:
## odds ratio 
##    1.19877

# under representation meth_t0_v_pod_dn and rna_t0_v_pod_dn
s1 = meth_t0_v_pod_dn
s2 = rna_t0_v_pod_dn
tot = meth_t0_v_pod_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1 , s2)) ,
  length(s2) - length(intersect(s1 , s2)),
  length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='less')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value = 6.601e-11
## alternative hypothesis: true odds ratio is less than 1
## 95 percent confidence interval:
##  0.0000000 0.8535117
## sample estimates:
## odds ratio 
##  0.8074922

Next for low vs high CRP.

• over representation meth_podcrp_dn and rna_podcrp_up OR=2.8 p=1.1e-8

• under representation meth_podcrp_up and rna_podcrp_up OR=0.79 p=0.23

• over representation meth_podcrp_up and rna_podcrp_dn OR=5.3 p=2.2e-16

• under representation meth_podcrp_dn and rna_podcrp_dn OR=1.2 p=0.80

# over representation meth_podcrp_dn and rna_podcrp_up
s1 = meth_podcrp_dn
s2 = rna_podcrp_up
tot = meth_podcrp_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1,s2)) ,
  length(s2) - length(intersect(s1,s2)),
  length(tot) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='greater')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(tot) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value = 1.09e-08
## alternative hypothesis: true odds ratio is greater than 1
## 95 percent confidence interval:
##  2.12314     Inf
## sample estimates:
## odds ratio 
##   2.823862

# under representation meth_podcrp_up and rna_podcrp_up
s1 = meth_podcrp_up
s2 = rna_podcrp_up
tot = meth_podcrp_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1 , s2)) ,
  length(s2) - length(intersect(s1 , s2)),
  length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='less')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value = 0.2327
## alternative hypothesis: true odds ratio is less than 1
## 95 percent confidence interval:
##  0.000000 1.248886
## sample estimates:
## odds ratio 
##  0.7892876

# over representation meth_podcrp_up and rna_podcrp_dn
s1 = meth_podcrp_up
s2 = rna_podcrp_dn
tot = meth_podcrp_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1,s2)) ,
  length(s2) - length(intersect(s1,s2)),
  length(tot) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='greater')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(tot) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is greater than 1
## 95 percent confidence interval:
##  4.094566      Inf
## sample estimates:
## odds ratio 
##   5.280571

# under representation meth_podcrp_dn and rna_podcrp_dn
s1 = meth_podcrp_dn
s2 = rna_podcrp_dn
tot = meth_podcrp_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1 , s2)) ,
  length(s2) - length(intersect(s1 , s2)),
  length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='less')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value = 0.8023
## alternative hypothesis: true odds ratio is less than 1
## 95 percent confidence interval:
##  0.000000 1.876394
## sample estimates:
## odds ratio 
##   1.200895

Number of gene sets

Starting with the RNA analysis.

#download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip", destfile="ReactomePathways.gmt.zip")
#unzip("ReactomePathways.gmt.zip")
genesets <- gmt_import("ReactomePathways.gmt")

mitch_bubbleplot <- function(res,n) {
top <- head(res$enrichment_result,n)
top <- top[order(top$s.dist),]
top$set <- substr(top$set,start=1,stop=60)
top$set <- factor(top$set, levels = top$set[order(top$s.dist)])
ggplot(top, aes(s.dist, set, size = setSize)) + geom_point(aes(colour=-log10(top$p.adjustANOVA)))
}

gt <- as.data.frame(rownames(rna_podcrp))
gt$genename <- sapply(strsplit(gt[,1]," "),"[[",2)
y <-mitch_import(x=rna_podcrp,DEtype="deseq2",geneTable=gt)

## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.

## Note: Mean no. genes in input = 21664

## Note: no. genes in output = 21631

## Note: estimated proportion of input genes in output = 0.998

yres <- mitch_calc(y,genesets=genesets)

## Note: When prioritising by significance (ie: small
##             p-values), large effect sizes might be missed.

mitch_bubbleplot(yres,30)

rna_podcrp_up <-  subset(yres$enrichment_result,p.adjustANOVA < 0.05 & s.dist > 0)$set
rna_podcrp_dn <-  subset(yres$enrichment_result,p.adjustANOVA < 0.05 & s.dist < 0)$set

gt <- as.data.frame(rownames(rna_t0crp))
gt$genename <- sapply(strsplit(gt[,1]," "),"[[",2)
y <-mitch_import(x=rna_t0crp,DEtype="deseq2",geneTable=gt)

## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.

## Note: Mean no. genes in input = 22158

## Note: no. genes in output = 22124

## Note: estimated proportion of input genes in output = 0.998

yres <- mitch_calc(y,genesets=genesets)

## Note: When prioritising by significance (ie: small
##             p-values), large effect sizes might be missed.

mitch_bubbleplot(yres,30)

rna_t0crp_up <-  subset(yres$enrichment_result,p.adjustANOVA < 0.05 & s.dist > 0)$set
rna_t0crp_dn <-  subset(yres$enrichment_result,p.adjustANOVA < 0.05 & s.dist < 0)$set

gt <- as.data.frame(rownames(rna_t0_v_pod))
gt$genename <- sapply(strsplit(gt[,1]," "),"[[",2)
y <-mitch_import(x=rna_t0_v_pod,DEtype="deseq2",geneTable=gt)

## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.

## Note: Mean no. genes in input = 21940

## Note: no. genes in output = 21907

## Note: estimated proportion of input genes in output = 0.998

yres <- mitch_calc(y,genesets=genesets)

## Note: When prioritising by significance (ie: small
##             p-values), large effect sizes might be missed.

mitch_bubbleplot(yres,30)

rna_t0_v_pod_up <-  subset(yres$enrichment_result,p.adjustANOVA < 0.05 & s.dist > 0)$set
rna_t0_v_pod_dn <-  subset(yres$enrichment_result,p.adjustANOVA < 0.05 & s.dist < 0)$set

run_mitch_rank <-function(dma){
  dmap <- dma[grep("Promoter_Associated",dma$Regulatory_Feature_Group),]
  dmap[which(dmap$UCSC_RefGene_Name==""),2] <- "NA"
  dmap$genename <- sapply(strsplit(dmap$UCSC_RefGene_Name,";"),"[[",1)
  dmap2 <- dmap[,c("genename","t")]
  rank <- aggregate(. ~ genename,dmap2,mean)
  rownames(rank) <- rank$genename
  rank$genename=NULL
  return(rank)
}

run_mitch_1d <- function(dma,name) {
  library("mitch")
  rank <- run_mitch_rank(dma)
  res <- mitch_calc(x = rank,genesets = genesets, priority = "significance",resrows=20)
  return(res)
}


yres <- run_mitch_1d(dma= meth_podcrp , name="meth_podcrp")

## Note: When prioritising by significance (ie: small
##             p-values), large effect sizes might be missed.

mitch_bubbleplot(yres,30)

meth_podcrp_up <-  subset(yres$enrichment_result,p.adjustANOVA < 0.05 & s.dist > 0)$set
meth_podcrp_dn <-  subset(yres$enrichment_result,p.adjustANOVA < 0.05 & s.dist < 0)$set

yres <- run_mitch_1d(dma= meth_t0crp , name="meth_t0crp")

## Note: When prioritising by significance (ie: small
##             p-values), large effect sizes might be missed.

mitch_bubbleplot(yres,30)

meth_t0crp_up <-  subset(yres$enrichment_result,p.adjustANOVA < 0.05 & s.dist > 0)$set
meth_t0crp_dn <-  subset(yres$enrichment_result,p.adjustANOVA < 0.05 & s.dist < 0)$set

yres <- run_mitch_1d(dma= meth_t0_v_pod , name="meth_t0_v_pod")

## Note: When prioritising by significance (ie: small
##             p-values), large effect sizes might be missed.

mitch_bubbleplot(yres,30)

meth_t0_v_pod_up <-  subset(yres$enrichment_result,p.adjustANOVA < 0.05 & s.dist > 0)$set
meth_t0_v_pod_dn <-  subset(yres$enrichment_result,p.adjustANOVA < 0.05 & s.dist < 0)$set

Bar plots on gene sets

xl <- list("RNA POD CRP up"=rna_podcrp_up,
  "RNA POD CRP down"=rna_podcrp_dn,
  "RNA t0 CRP up"=rna_t0crp_up,
  "RNA t0 CRP down"=rna_t0crp_dn,
  "RNA t0 Vs POD up"=rna_t0_v_pod_up,
  "RNA t0 Vs POD down"=rna_t0_v_pod_dn,
  "Meth POD CRP up"=meth_podcrp_up,
  "Meth POD CRP down"=meth_podcrp_dn,
  "Meth t0 CRP up"=meth_t0crp_up,
  "Meth t0 CRP down"=meth_t0crp_dn,
  "Meth t0 Vs POD up"=meth_t0_v_pod_up,
  "Meth t0 Vs POD down"=meth_t0_v_pod_dn)

xlr <- sapply(xl,length)

xlr_up <- xlr[grep("up",names(xlr))]

xlr_dn <- xlr[grep("down",names(xlr))]

MAX=max(xlr_up)*1.2
MIN=-max(xlr_dn)*1.2

par(mar=c(5,5,5,5))
names(xlr_up) <- gsub(" up","",names(xlr_up))
barplot(xlr_up,ylim=c(MIN,MAX),cex.names=0.7,ylab="no. gene sets up- and down-regulated",col="lightblue")
text(x = (1:length(xlr_up)*1.2)-0.5  , y = xlr_up , label = xlr_up, pos = 3, cex = 1, col = "black")
barplot(-xlr_dn,ylim=c(MIN,MAX),cex.names=0.5, col="pink",add=T, xaxt='n')
text(x = (1:length(xlr_up)*1.2)-0.5  , y = -xlr_dn-50 , label = xlr_dn, pos = 3, cex = 1, col = "black")

MAX=max(xlr_up)*1.4
MIN=-max(xlr_dn)*1.4

par(mar=c(5,15,5,5))
barplot(xlr_up,xlim=c(MIN,MAX),cex.names=1,main="no. gene sets up- and down-regulated",col="lightblue",horiz=TRUE,las=1)
text(y = (1:length(xlr_up)*1.2)-0.75  , x = xlr_up + 50 , label = xlr_up, pos = 3, cex = 1, col = "black")
barplot(-xlr_dn,ylim=c(MIN,MAX),cex.names=0.5, col="pink",add=T, yaxt='n',horiz=TRUE,las=1)
text(y = (1:length(xlr_up)*1.2)-0.75  , x = -xlr_dn-50 , label = xlr_dn, pos = 3, cex = 1, col = "black")

Venn diagrams on gene sets

Now to identify the gene set overlaps (same omics different contrasts).

v1 <- list("RNA POD CRP up"=rna_podcrp_up,
  "RNA POD CRP down"=rna_podcrp_dn,
  "RNA t0 CRP up"=rna_t0crp_up,
  "RNA t0 CRP down"=rna_t0crp_dn,
  "RNA t0 Vs POD up"=rna_t0_v_pod_up,
  "RNA t0 Vs POD down"=rna_t0_v_pod_dn)

plot(euler(v1),quantities = TRUE)

v1 <- list( "Meth POD CRP up"=meth_podcrp_up,
  "Meth POD CRP down"=meth_podcrp_dn,
  "Meth t0 CRP up"=meth_t0crp_up,
  "Meth t0 CRP down"=meth_t0crp_dn,
  "Meth t0 Vs POD up"=meth_t0_v_pod_up,
  "Meth t0 Vs POD down"=meth_t0_v_pod_dn)

plot(euler(v1),quantities = TRUE)

Now to identify the gene set overlaps (same contrast different omics).

v1 <- list("RNA POD CRP up"=rna_podcrp_up,
  "RNA POD CRP down"=rna_podcrp_dn,
  "Meth POD CRP up"=meth_podcrp_up,
  "Meth POD CRP down"=meth_podcrp_dn)

plot(euler(v1),quantities = TRUE)

## Warning in colSums(id & !empty) == 0 | merged_sets: longer object length is not
## a multiple of shorter object length

v1 <-  list( "RNA t0 CRP up"=rna_t0crp_up,
  "RNA t0 CRP down"=rna_t0crp_dn,
  "Meth t0 CRP up"=meth_t0crp_up,
  "Meth t0 CRP down"=meth_t0crp_dn)

plot(euler(v1),quantities = TRUE)

v1 <- list("RNA t0 Vs POD up"=rna_t0_v_pod_up,
  "RNA t0 Vs POD down"=rna_t0_v_pod_dn,
  "Meth t0 Vs POD up"=meth_t0_v_pod_up,
  "Meth t0 Vs POD down"=meth_t0_v_pod_dn)

plot(euler(v1),quantities = TRUE)

Overrepresentaton test on gene sets

Using the fisher test based on this tutorial https://seqqc.wordpress.com/2019/07/25/how-to-use-phyper-in-r/

Firstly for baseline vs post-op.

• over representation meth_t0_v_pod_dn and rna_t0_v_pod_up OR=2.1 p=2.2e-16

• under representation meth_t0_v_pod_up and rna_t0_v_pod_up OR=1.3 p=1

• over representation meth_t0_v_pod_up and rna_t0_v_pod_dn OR=1.2 p=2.5e-6

• under representation meth_t0_v_pod_dn and rna_t0_v_pod_dn OR=0.81 p-6.6e-11

# over representation meth_t0_v_pod_dn and rna_t0_v_pod_up
s1 = meth_t0_v_pod_dn
s2 = rna_t0_v_pod_up
tot = meth_t0_v_pod_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1,s2)) ,
  length(s2) - length(intersect(s1,s2)),
  length(tot) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='greater')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(tot) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value = 1
## alternative hypothesis: true odds ratio is greater than 1
## 95 percent confidence interval:
##    0 Inf
## sample estimates:
## odds ratio 
##          0

# under representation meth_t0_v_pod_up and rna_t0_v_pod_up
s1 = meth_t0_v_pod_up
s2 = rna_t0_v_pod_up
tot = meth_t0_v_pod_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1 , s2)) ,
  length(s2) - length(intersect(s1 , s2)),
  length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='less')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value = 0.9999
## alternative hypothesis: true odds ratio is less than 1
## 95 percent confidence interval:
##   0.00000 17.90319
## sample estimates:
## odds ratio 
##   7.744478

# over representation meth_t0_v_pod_up and rna_t0_v_pod_dn
s1 = meth_t0_v_pod_up
s2 = rna_t0_v_pod_dn
tot = meth_t0_v_pod_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1,s2)) ,
  length(s2) - length(intersect(s1,s2)),
  length(tot) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='greater')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(tot) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is greater than 1
## 95 percent confidence interval:
##  327.1081      Inf
## sample estimates:
## odds ratio 
##   626.0282

# under representation meth_t0_v_pod_dn and rna_t0_v_pod_dn
s1 = meth_t0_v_pod_dn
s2 = rna_t0_v_pod_dn
tot = meth_t0_v_pod_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1 , s2)) ,
  length(s2) - length(intersect(s1 , s2)),
  length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='less')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value = 0.9931
## alternative hypothesis: true odds ratio is less than 1
## 95 percent confidence interval:
##     0.000 1026.969
## sample estimates:
## odds ratio 
##          0

Next for low vs high CRP.

• over representation meth_podcrp_dn and rna_podcrp_up OR=2.8 p=1.1e-8

• under representation meth_podcrp_up and rna_podcrp_up OR=0.79 p=0.23

• over representation meth_podcrp_up and rna_podcrp_dn OR=5.3 p=2.2e-16

• under representation meth_podcrp_dn and rna_podcrp_dn OR=1.2 p=0.80

# over representation meth_podcrp_dn and rna_podcrp_up
s1 = meth_podcrp_dn
s2 = rna_podcrp_up
tot = meth_podcrp_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1,s2)) ,
  length(s2) - length(intersect(s1,s2)),
  length(tot) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='greater')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(tot) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value = 1
## alternative hypothesis: true odds ratio is greater than 1
## 95 percent confidence interval:
##    0 Inf
## sample estimates:
## odds ratio 
##          0

# under representation meth_podcrp_up and rna_podcrp_up
s1 = meth_podcrp_up
s2 = rna_podcrp_up
tot = meth_podcrp_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1 , s2)) ,
  length(s2) - length(intersect(s1 , s2)),
  length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='less')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value = 0.9986
## alternative hypothesis: true odds ratio is less than 1
## 95 percent confidence interval:
##    0.0000 107.3482
## sample estimates:
## odds ratio 
##    19.0969

# over representation meth_podcrp_up and rna_podcrp_dn
s1 = meth_podcrp_up
s2 = rna_podcrp_dn
tot = meth_podcrp_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1,s2)) ,
  length(s2) - length(intersect(s1,s2)),
  length(tot) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='greater')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(tot) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value = 8.766e-15
## alternative hypothesis: true odds ratio is greater than 1
## 95 percent confidence interval:
##  77.13462      Inf
## sample estimates:
## odds ratio 
##    225.731

# under representation meth_podcrp_dn and rna_podcrp_dn
s1 = meth_podcrp_dn
s2 = rna_podcrp_dn
tot = meth_podcrp_all

fisher.test(matrix(c(
  length(intersect(s1 , s2)),
  length(s1) - length(intersect(s1 , s2)) ,
  length(s2) - length(intersect(s1 , s2)),
  length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1,s2))
), 2, 2), alternative='less')

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(length(intersect(s1, s2)), length(s1) - length(intersect(s1, s2)), length(s2) - length(intersect(s1, s2)), length(meth_t0_v_pod_all) - length(s1) - length(s2) + length(intersect(s1, s2))), 2, 2)
## p-value = 1
## alternative hypothesis: true odds ratio is less than 1
## 95 percent confidence interval:
##    0 Inf
## sample estimates:
## odds ratio 
##          0

Session Info

Session information.

sessionInfo()

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggplot2_3.3.3    UpSetR_1.4.0     eulerr_6.1.0     kableExtra_1.3.4
## [5] mitch_1.4.0      gplots_3.1.1    
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.4.2         sass_0.4.0         jsonlite_1.7.2     viridisLite_0.4.0 
##  [5] gtools_3.8.2       bslib_0.2.5.1      shiny_1.6.0        assertthat_0.2.1  
##  [9] highr_0.9          yaml_2.2.1         pillar_1.6.1       glue_1.4.2        
## [13] digest_0.6.27      RColorBrewer_1.1-2 promises_1.2.0.1   polyclip_1.10-0   
## [17] rvest_1.0.0        colorspace_2.0-1   htmltools_0.5.1.1  httpuv_1.6.1      
## [21] plyr_1.8.6         pkgconfig_2.0.3    purrr_0.3.4        xtable_1.8-4      
## [25] scales_1.1.1       webshot_0.5.2      svglite_2.0.0      later_1.2.0       
## [29] tibble_3.1.2       echarts4r_0.4.0    generics_0.1.0     farver_2.1.0      
## [33] ellipsis_0.3.2     withr_2.4.2        magrittr_2.0.1     crayon_1.4.1      
## [37] mime_0.10          evaluate_0.14      GGally_2.1.1       fansi_0.5.0       
## [41] MASS_7.3-54        xml2_1.3.2         beeswarm_0.3.1     tools_4.1.0       
## [45] lifecycle_1.0.0    stringr_1.4.0      munsell_0.5.0      compiler_4.1.0    
## [49] jquerylib_0.1.4    caTools_1.18.2     systemfonts_1.0.2  rlang_0.4.11      
## [53] grid_4.1.0         rstudioapi_0.13    htmlwidgets_1.5.3  labeling_0.4.2    
## [57] bitops_1.0-7       rmarkdown_2.8      gtable_0.3.0       DBI_1.1.1         
## [61] reshape_0.8.8      reshape2_1.4.4     R6_2.5.0           gridExtra_2.3     
## [65] knitr_1.33         dplyr_1.0.6        fastmap_1.1.0      utf8_1.2.1        
## [69] KernSmooth_2.23-20 polylabelr_0.2.0   stringi_1.6.2      parallel_4.1.0    
## [73] Rcpp_1.0.6         vctrs_0.3.8        tidyselect_1.1.1   xfun_0.23