Autism Spectrum Disorder Methylation analysis - comparison of Guthrie card and fresh blood
Mark Ziemann
2023-01-19
Source: https://github.com/markziemann/asd_meth
Introduction
Here I will be running a comparison of differential methylation data from guthrie and fresh blood samples.
Here are the files that I’m using:
ASD_blood_top_dmps_onADOS_withintw_Nov27_limma.csv
ASD_guthrie_top_dmps_onADOS_withintw_Nov27_limma.csv
suppressPackageStartupMessages({
library("parallel")
library("mitch")
library("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
source("https://raw.githubusercontent.com/markziemann/gmea/main/meth_functions.R")
library("data.table")
library("kableExtra")
library("eulerr")
library("RIdeogram")
library("GenomicRanges")
library("tictoc")
})Load data
bl_design <- readRDS(file="bl_design.rds")
bl_mvals <- readRDS(file="bl_mvals.rds")
gu_design <- readRDS(file="gu_design.rds")
gu_mvals <- readRDS(file="gu_mvals.rds")Probe sets
For each gene, extract out the probes.
anno <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
myann <- data.frame(anno[,c("UCSC_RefGene_Name","Regulatory_Feature_Group","Islands_Name","Relation_to_Island")])
promoters <- grep("Prom",myann$Regulatory_Feature_Group)
promoters <- myann[promoters,]
gp <- myann[,"UCSC_RefGene_Name",drop=FALSE]
gp2 <- strsplit(gp$UCSC_RefGene_Name,";")
names(gp2) <- rownames(gp)
sets <- split(rep(names(gp2), lengths(gp2)), unlist(gp2))
summary(unlist(lapply(sets,length)))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 9.00 24.00 49.68 55.00 6778.00Load Limma and RUV data for manhattan plot
First I will generate plots for limmma analysis.
ASD_blood_top_dmps_onADOS_Nov27_withintw_limma_top1k.csv
ASD_guthrie_top_dmps_onADOS_Nov27_withintw_limma_top1k.csv
# blood at assessment
top <- read.csv("ASD_blood_top_dmps_onADOS_withintw_Nov27_limma.csv")
nrow(top)## [1] 802647top <- subset(top,P.Value<1e-2)
nrow(top)## [1] 3573-log10(min(top$P.Value))## [1] 4.641884top$chr <- as.integer(gsub("chr","",top$chr))
top$snp <- paste(top$Name,top$UCSC_RefGene_Name)
top$snp <- sapply(strsplit(top$snp,";"),"[[",1)
up <- subset(top,logFC>0)
dn <- subset(top,logFC<0)
par(mfrow=c(2,1))
par(mar=c(4,4,3,5))
manhattan(x=up,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 ,main="Blood at assessment limma hypermethylated")
manhattan(x=dn,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 , main="Blood at assessment limma hypomethylated")
pdf("manhat_limma_bl.pdf")
par(mfrow=c(2,1))
par(mar=c(4,4,3,5))
manhattan(x=up,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 ,main="Blood at assessment limma hypermethylated")
manhattan(x=dn,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 , main="Blood at assessment limma hypomethylated")
dev.off()## png
## 2# neonatal guthrie card
top <- read.csv("ASD_guthrie_top_dmps_onADOS_withintw_Nov27_limma.csv")
nrow(top)## [1] 790658top <- subset(top,P.Value<1e-2)
nrow(top)## [1] 3052-log10(min(top$P.Value))## [1] 5.396928top$chr <- as.integer(gsub("chr","",top$chr))
top$snp <- paste(top$Name,top$UCSC_RefGene_Name)
top$snp <- sapply(strsplit(top$snp,";"),"[[",1)
up <- subset(top,logFC>0)
dn <- subset(top,logFC<0)
par(mfrow=c(2,1))
par(mar=c(4,4,3,5))
manhattan(x=up,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 ,main="Neonatal Guthrie card limma hypermethylated",
annotateTop=FALSE)
manhattan(x=dn,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 , main="Neonatal Guthrie card limma hypomethylated",
annotateTop=FALSE)
pdf("manhat_limma_gu.pdf")
par(mfrow=c(2,1))
par(mar=c(4,4,3,5))
manhattan(x=up,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 ,main="Neonatal Guthrie card limma hypermethylated",
annotateTop=FALSE)
manhattan(x=dn,chr="chr",bp="pos",p="P.Value",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-04 , main="Neonatal Guthrie card limma hypomethylated",
annotateTop=FALSE)
dev.off()## png
## 2Now I will do the same for the RUV analysis.
ASD_blood_topruv80pc_1kdmps_onADOS_withintw_Nov29.csv
ASD_guthrie_topruv80%_1kdmps_onADOS_withintw_Nov29.csv
# blood at assessment
top <- read.csv("ASD_blood_topruv80pc_dmps_onADOS_withintw_Nov29.csv")
nrow(top)## [1] 802647top <- subset(top,F.p<1e-2)
nrow(top)## [1] 22779-log10(min(top$F.p))## [1] 6.028082top$chr <- as.integer(gsub("chr","",top$chr))
top$snp <- paste(top$Name,top$UCSC_RefGene_Name)
top$snp <- sapply(strsplit(top$snp,";"),"[[",1)
up <- subset(top,b_X1>0)
dn <- subset(top,b_X1<0)
par(mfrow=c(2,1))
par(mar=c(4,4,3,5))
manhattan(x=up,chr="chr",bp="pos",p="F.p",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-05 ,main="Blood at assessment RUV hypermethylated",
annotateTop=FALSE)
manhattan(x=dn,chr="chr",bp="pos",p="F.p",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-05 ,main="Blood at assessment RUV hypomethylated",
annotateTop=FALSE)
pdf("manhat_ruv_bl.pdf")
par(mfrow=c(2,1))
par(mar=c(4,4,3,5))
manhattan(x=up,chr="chr",bp="pos",p="F.p",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-05 ,main="Blood at assessment RUV hypermethylated",
annotateTop=FALSE)
manhattan(x=dn,chr="chr",bp="pos",p="F.p",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 1e-05 ,main="Blood at assessment RUV hypomethylated",
annotateTop=FALSE)
dev.off()## png
## 2# Neonatal guthrie card
top <- read.csv("ASD_guthrie_topruv80pc_dmps_onADOS_withintw_Nov29.csv")
nrow(top)## [1] 790658top <- subset(top,F.p<1e-2)
nrow(top)## [1] 9277-log10(min(top$F.p))## [1] 5.397192top$chr <- as.integer(gsub("chr","",top$chr))
top$snp <- paste(top$Name,top$UCSC_RefGene_Name)
top$snp <- sapply(strsplit(top$snp,";"),"[[",1)
up <- subset(top,b_X1>0)
dn <- subset(top,b_X1<0)
par(mfrow=c(2,1))
par(mar=c(4,4,3,5))
manhattan(x=up,chr="chr",bp="pos",p="F.p",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 5e-05 ,main="Neonatal Guthrie card RUV hypermethylated",
annotateTop=FALSE)
manhattan(x=dn,chr="chr",bp="pos",p="F.p",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 5e-05 ,main="Neonatal Guthrie card RUV hypomethylated",
annotateTop=FALSE)
pdf("manhat_ruv_gu.pdf")
par(mfrow=c(2,1))
par(mar=c(4,4,3,5))
manhattan(x=up,chr="chr",bp="pos",p="F.p",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 5e-05 ,main="Neonatal Guthrie card RUV hypermethylated",
annotateTop=FALSE)
manhattan(x=dn,chr="chr",bp="pos",p="F.p",snp="snp",suggestiveline = -log10(1e-04),
ylim=c(2,6), annotatePval= 5e-05 ,main="Neonatal Guthrie card RUV hypomethylated",
annotateTop=FALSE)
dev.off()## png
## 2Compartment enrichment
Looking for enrichments in different genomic compartments.
compartment_enrichment <- function(dma) {
up <- subset(dma,b_X1>0 & F.p<1e-4)
dn <- subset(dma,b_X1<0 & F.p<1e-4)
all <- table(unique(dma)$Regulatory_Feature_Group)
up <- table(unique(up)$Regulatory_Feature_Group)
dn <- table(unique(dn)$Regulatory_Feature_Group)
xx=NULL
xx <- merge(as.data.frame(all, row.names = 1),as.data.frame(up,row.names = 1),by=0, all = TRUE)
rownames(xx) <- xx[,1]
rownames(xx)[1] <- "Intergenic"
xx[,1] = NULL
colnames(xx) <- c("all","up")
xx[is.na(xx)] <- 0
head(xx)
x=xx$up
m=xx$all
n=sum(xx$all)-xx$all
k=sum(xx$up)
xl <- apply(xx,1,function(x) {
mat <- matrix(c(x[2],x[1]-x[2], sum(xx$up)-x[2], sum(xx$all) - sum(xx$up) -x [1] + x[2] ),2,2)
mat
fisher.test(mat)
})
xx$OR <- unname(unlist(lapply(X=xl, FUN = function(x) {x$estimate})))
xx$fisherPval <- unname(unlist(lapply(X=xl, FUN = function(x) {x$p.value})))
xx$lowerCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[1]]})))
xx$upperCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[2]]})))
up_comp <- xx
xx=NULL
xx <- merge(as.data.frame(all, row.names = 1),as.data.frame(dn,row.names = 1),by=0, all = TRUE)
rownames(xx) <- xx[,1]
rownames(xx)[1] <- "Intergenic"
xx[,1] = NULL
colnames(xx) <- c("all","dn")
xx[is.na(xx)] <- 0
x=xx$dn
m=xx$all
n=sum(xx$all)-xx$all
k=sum(xx$dn)
xl <- apply(xx,1,function(x) {
mat <- matrix(c(x[2],x[1]-x[2], sum(xx$dn)-x[2], sum(xx$all) - sum(xx$dn) -x [1] + x[2] ),2,2)
mat
fisher.test(mat)
})
xx$OR <- unname(unlist(lapply(X=xl, FUN = function(x) {x$estimate})))
xx$fisherPval <- unname(unlist(lapply(X=xl, FUN = function(x) {x$p.value})))
xx$lowerCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[1]]})))
xx$upperCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[2]]})))
dn_comp <- xx
list("up_comp"=up_comp,"dn_comp"=dn_comp)
}
make_forest_plots <- function(comp) {
comp_data <-
structure(list(
"mean" = comp$up_comp$OR ,
"lower" = comp$up_comp$lowerCI ,
"upper" = comp$up_comp$upperCI ,
.Names = c("mean", "lower", "upper"),
row.names = c(NA, -11L),
class = "data.frame"))
comp_data <- as.data.frame(comp_data[1:3],row.names = rownames(comp$up_comp) )
forestplot(comp_data,title = "hypermethylated",
labeltext = as.list(rownames(comp_data)),
mean=mean,lower=lower,upper=upper)
comp_data <-
structure(list(
"mean" = comp$dn_comp$OR ,
"lower" = comp$dn_comp$lowerCI ,
"upper" = comp$dn_comp$upperCI ,
.Names = c("mean", "lower", "upper"),
row.names = c(NA, -11L),
class = "data.frame"))
comp_data <- as.data.frame(comp_data[1:3],row.names = rownames(comp$dn_comp) )
forestplot(comp_data,title = "hypomethylated",
labeltext = as.list(rownames(comp_data)),
mean=mean,lower=lower,upper=upper)
}
par(mfrow=c(2,1))
# guthrie
dma1 <- read.csv("ASD_guthrie_topruv80pc_dmps_onADOS_withintw_Nov29.csv")
comp <- compartment_enrichment(dma1)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp## $up_comp
## all up OR fisherPval lowerCI
## Intergenic 586499 20 0.7735433 0.5264496 0.3366793
## Promoter_Associated 103401 3 0.7669001 1.0000000 0.1485921
## Unclassified 39902 2 1.3937071 0.6566854 0.1606349
## Unclassified_Cell_type_specific 47673 4 2.4937161 0.0944339 0.6307097
## upperCI
## Intergenic 1.929357
## Promoter_Associated 2.502477
## Unclassified 5.546076
## Unclassified_Cell_type_specific 7.222354
##
## $dn_comp
## all dn OR fisherPval lowerCI upperCI
## Intergenic 586499 9 0.7832138 0.75163247 0.2185802 3.480538
## Promoter_Associated 103401 4 2.9540820 0.07904552 0.6647220 10.584600make_forest_plots(comp)
# blood
dma2 <- read.csv("ASD_blood_topruv80pc_dmps_onADOS_withintw_Nov29.csv")
comp <- compartment_enrichment(dma2)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp## $up_comp
## all up OR fisherPval
## Intergenic 597228 185 0.7313138 0.018019994
## Gene_Associated_Cell_type_specific 2722 1 1.0844344 0.603132121
## NonGene_Associated 1469 2 4.0439384 0.089376202
## Promoter_Associated 103706 38 1.0945386 0.587623365
## Promoter_Associated_Cell_type_specific 6922 1 0.4241172 0.735543881
## Unclassified 40194 24 1.8361662 0.007532495
## Unclassified_Cell_type_specific 48172 21 1.3105087 0.248151382
## lowerCI upperCI
## Intergenic 0.56382230 0.9548131
## Gene_Associated_Cell_type_specific 0.02734596 6.0986789
## NonGene_Associated 0.48682443 14.7803141
## Promoter_Associated 0.75546916 1.5470624
## Promoter_Associated_Cell_type_specific 0.01069972 2.3837941
## Unclassified 1.15431372 2.7966714
## Unclassified_Cell_type_specific 0.79704067 2.0476139
##
## $dn_comp
## all dn OR fisherPval
## Intergenic 597228 130 1.4424506 0.07061926
## Promoter_Associated 103706 13 0.5919452 0.07679495
## Promoter_Associated_Cell_type_specific 6922 1 0.7184342 1.00000000
## Unclassified 40194 6 0.7342520 0.58746392
## Unclassified_Cell_type_specific 48172 11 1.1485873 0.61742215
## lowerCI upperCI
## Intergenic 0.9687511 2.209144
## Promoter_Associated 0.3078952 1.043872
## Promoter_Associated_Cell_type_specific 0.0180750 4.062341
## Unclassified 0.2653658 1.634975
## Unclassified_Cell_type_specific 0.5612410 2.114927make_forest_plots(comp)
Compartments top 1000
compartment_enrichment2 <- function(dma) {
all <- table(unique(dma)$Regulatory_Feature_Group)
dma <- head(dma,1000)
up <- subset(dma,b_X1>0 )
dn <- subset(dma,b_X1<0 )
up <- table(unique(up)$Regulatory_Feature_Group)
dn <- table(unique(dn)$Regulatory_Feature_Group)
xx=NULL
xx <- merge(as.data.frame(all, row.names = 1),as.data.frame(up,row.names = 1),by=0, all = TRUE)
rownames(xx) <- xx[,1]
rownames(xx)[1] <- "Intergenic"
xx[,1] = NULL
colnames(xx) <- c("all","up")
xx[is.na(xx)] <- 0
head(xx)
x=xx$up
m=xx$all
n=sum(xx$all)-xx$all
k=sum(xx$up)
xl <- apply(xx,1,function(x) {
mat <- matrix(c(x[2],x[1]-x[2], sum(xx$up)-x[2], sum(xx$all) - sum(xx$up) -x [1] + x[2] ),2,2)
mat
fisher.test(mat)
})
xx$OR <- unname(unlist(lapply(X=xl, FUN = function(x) {x$estimate})))
xx$fisherPval <- unname(unlist(lapply(X=xl, FUN = function(x) {x$p.value})))
xx$lowerCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[1]]})))
xx$upperCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[2]]})))
up_comp <- xx
xx=NULL
xx <- merge(as.data.frame(all, row.names = 1),as.data.frame(dn,row.names = 1),by=0, all = TRUE)
rownames(xx) <- xx[,1]
rownames(xx)[1] <- "Intergenic"
xx[,1] = NULL
colnames(xx) <- c("all","dn")
xx[is.na(xx)] <- 0
x=xx$dn
m=xx$all
n=sum(xx$all)-xx$all
k=sum(xx$dn)
xl <- apply(xx,1,function(x) {
mat <- matrix(c(x[2],x[1]-x[2], sum(xx$dn)-x[2], sum(xx$all) - sum(xx$dn) -x [1] + x[2] ),2,2)
mat
fisher.test(mat)
})
xx$OR <- unname(unlist(lapply(X=xl, FUN = function(x) {x$estimate})))
xx$fisherPval <- unname(unlist(lapply(X=xl, FUN = function(x) {x$p.value})))
xx$lowerCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[1]]})))
xx$upperCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[2]]})))
dn_comp <- xx
list("up_comp"=up_comp,"dn_comp"=dn_comp)
}
par(mfrow=c(2,1))
# guthrie
dma1 <- read.csv("ASD_guthrie_topruv80pc_dmps_onADOS_withintw_Nov29.csv")
comp <- compartment_enrichment2(dma1)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp## $up_comp
## all up OR fisherPval
## Intergenic 586499 506 1.2065930 4.876302e-02
## Gene_Associated 1950 2 1.2447575 6.777329e-01
## Gene_Associated_Cell_type_specific 2678 2 0.9052893 1.000000e+00
## Promoter_Associated 103401 46 0.5043021 1.231374e-06
## Promoter_Associated_Cell_type_specific 6837 5 0.8858826 1.000000e+00
## Unclassified 39902 35 1.0673950 6.550857e-01
## Unclassified_Cell_type_specific 47673 56 1.4648926 8.291708e-03
## lowerCI upperCI
## Intergenic 1.0017218 1.4607631
## Gene_Associated 0.1503354 4.5202421
## Gene_Associated_Cell_type_specific 0.1093455 3.2866902
## Promoter_Associated 0.3651993 0.6811623
## Promoter_Associated_Cell_type_specific 0.2865813 2.0790122
## Unclassified 0.7364585 1.5010562
## Unclassified_Cell_type_specific 1.0931475 1.9289815
##
## $dn_comp
## all dn OR fisherPval
## Intergenic 586499 238 0.7530572 0.016713246
## Gene_Associated 1950 2 2.3392145 0.212214450
## Gene_Associated_Cell_type_specific 2678 3 2.5603634 0.115588582
## NonGene_Associated 1463 2 3.1209575 0.136427860
## Promoter_Associated 103401 37 0.7906693 0.202712147
## Promoter_Associated_Cell_type_specific 6837 5 1.6716751 0.234181666
## Unclassified 39902 21 1.2084103 0.389861123
## Unclassified_Cell_type_specific 47673 40 2.0248912 0.000103447
## lowerCI upperCI
## Intergenic 0.5982076 0.9528608
## Gene_Associated 0.2819776 8.5254474
## Gene_Associated_Cell_type_specific 0.5251045 7.5537539
## NonGene_Associated 0.3761146 11.3796560
## Promoter_Associated 0.5463753 1.1141385
## Promoter_Associated_Cell_type_specific 0.5391825 3.9406503
## Unclassified 0.7378654 1.8787372
## Unclassified_Cell_type_specific 1.4184788 2.8216572make_forest_plots(comp)
# blood
dma2 <- read.csv("ASD_blood_topruv80pc_dmps_onADOS_withintw_Nov29.csv")
comp <- compartment_enrichment2(dma2)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp## $up_comp
## all up OR fisherPval
## Intergenic 597228 392 0.7247122 0.0004074401
## Gene_Associated_Cell_type_specific 2722 2 1.0204165 0.7239169330
## NonGene_Associated 1469 3 2.8492119 0.0910107247
## Promoter_Associated 103706 82 1.1142721 0.3524437133
## Promoter_Associated_Cell_type_specific 6922 3 0.5995968 0.5007512117
## Unclassified 40194 46 1.6409242 0.0021670963
## Unclassified_Cell_type_specific 48172 50 1.4836366 0.0107321764
## lowerCI upperCI
## Intergenic 0.6070772 0.8676483
## Gene_Associated_Cell_type_specific 0.1232130 3.7065876
## NonGene_Associated 0.5852007 8.3858998
## Promoter_Associated 0.8710170 1.4101953
## Promoter_Associated_Cell_type_specific 0.1232574 1.7618248
## Unclassified 1.1865714 2.2206049
## Unclassified_Cell_type_specific 1.0869761 1.9854700
##
## $dn_comp
## all dn OR fisherPval
## Intergenic 597228 347 1.5916450 1.763790e-04
## Gene_Associated 1976 1 0.9624463 1.000000e+00
## Gene_Associated_Cell_type_specific 2722 1 0.6979273 1.000000e+00
## Promoter_Associated 103706 20 0.3351772 2.495065e-08
## Promoter_Associated_Cell_type_specific 6922 5 1.3786143 4.214949e-01
## Unclassified 40194 13 0.6028010 7.328577e-02
## Unclassified_Cell_type_specific 48172 35 1.4167569 5.127794e-02
## lowerCI upperCI
## Intergenic 1.23685872 2.0715126
## Gene_Associated 0.02429785 5.3972125
## Gene_Associated_Cell_type_specific 0.01762221 3.9128087
## Promoter_Associated 0.20253516 0.5246626
## Promoter_Associated_Cell_type_specific 0.44516373 3.2443629
## Unclassified 0.31819539 1.0425638
## Unclassified_Cell_type_specific 0.97234103 2.0051041make_forest_plots(comp)
Looking for enrichments in proximity to CGI contexts
cgi_enrichment <- function(dma) {
dma$Relation_to_Island <- gsub("N_","",dma$Relation_to_Island)
dma$Relation_to_Island <- gsub("S_","",dma$Relation_to_Island)
all <- table(unique(dma)$Relation_to_Island)
up <- subset(dma,b_X1>0 & F.p < 1e-4)
dn <- subset(dma,b_X1<0 & F.p < 1e-4)
up <- table(unique(up)$Relation_to_Island)
dn <- table(unique(dn)$Relation_to_Island)
xx=NULL
xx <- merge(as.data.frame(all, row.names = 1),as.data.frame(up,row.names = 1),by=0, all = TRUE)
rownames(xx) <- xx[,1]
xx[,1] = NULL
colnames(xx) <- c("all","up")
xx[is.na(xx)] <- 0
head(xx)
x=xx$up
m=xx$all
n=sum(xx$all)-xx$all
k=sum(xx$up)
xl <- apply(xx,1,function(x) {
mat <- matrix(c(x[2],x[1]-x[2], sum(xx$up)-x[2], sum(xx$all) - sum(xx$up) -x [1] + x[2] ),2,2)
mat
fisher.test(mat)
})
xx$OR <- unname(unlist(lapply(X=xl, FUN = function(x) {x$estimate})))
xx$fisherPval <- unname(unlist(lapply(X=xl, FUN = function(x) {x$p.value})))
xx$lowerCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[1]]})))
xx$upperCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[2]]})))
up_comp <- xx
xx=NULL
xx <- merge(as.data.frame(all, row.names = 1),as.data.frame(dn,row.names = 1),by=0, all = TRUE)
rownames(xx) <- xx[,1]
xx[,1] = NULL
colnames(xx) <- c("all","dn")
xx[is.na(xx)] <- 0
x=xx$dn
m=xx$all
n=sum(xx$all)-xx$all
k=sum(xx$dn)
xl <- apply(xx,1,function(x) {
mat <- matrix(c(x[2],x[1]-x[2], sum(xx$dn)-x[2], sum(xx$all) - sum(xx$dn) -x [1] + x[2] ),2,2)
mat
fisher.test(mat)
})
xx$OR <- unname(unlist(lapply(X=xl, FUN = function(x) {x$estimate})))
xx$fisherPval <- unname(unlist(lapply(X=xl, FUN = function(x) {x$p.value})))
xx$lowerCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[1]]})))
xx$upperCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[2]]})))
dn_comp <- xx
list("up_comp"=up_comp,"dn_comp"=dn_comp)
}
par(mfrow=c(2,1))
# guthrie
comp <- cgi_enrichment(dma1)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp## $up_comp
## all up OR fisherPval lowerCI upperCI
## Island 150182 3 0.4920734 0.3425155 0.09532123 1.605679
## OpenSea 442427 17 1.1150168 0.8528643 0.50188448 2.558864
## Shelf 55491 1 0.4731480 0.7196312 0.01157372 2.862179
## Shore 142558 8 1.7319296 0.2220033 0.66339359 4.070420
##
## $dn_comp
## all dn OR fisherPval lowerCI upperCI
## OpenSea 442427 4 0.3498128 0.09224566 0.07869422 1.253381
## Shelf 55491 4 5.8876860 0.01036871 1.32500314 21.098177
## Shore 142558 5 2.8413884 0.06846775 0.73136878 9.851305make_forest_plots(comp)
# blood
comp <- cgi_enrichment(dma2)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp## $up_comp
## all up OR fisherPval lowerCI upperCI
## Island 150186 35 0.6414918 0.01263722 0.4363857 0.9175731
## OpenSea 452360 161 1.1231731 0.35951610 0.8764598 1.4436709
## Shelf 56174 11 0.5599464 0.05658912 0.2760940 1.0187081
## Shore 143927 65 1.4373430 0.01396625 1.0704631 1.9079306
##
## $dn_comp
## all dn OR fisherPval lowerCI upperCI
## Island 150186 20 0.6161739 0.04275833 0.3652212 0.988168
## OpenSea 452360 99 1.2365200 0.20382121 0.8912209 1.727134
## Shelf 56174 17 1.5689439 0.08708702 0.8893292 2.601255
## Shore 143927 25 0.8412901 0.47295813 0.5259715 1.295663make_forest_plots(comp)
CGI enrichments in top 1000 probes
cgi_enrichment2 <- function(dma) {
dma$Relation_to_Island <- gsub("N_","",dma$Relation_to_Island)
dma$Relation_to_Island <- gsub("S_","",dma$Relation_to_Island)
all <- table(unique(dma)$Relation_to_Island)
dma <- head(dma,1000)
up <- subset(dma,b_X1>0 )
dn <- subset(dma,b_X1<0 )
up <- table(unique(up)$Relation_to_Island)
dn <- table(unique(dn)$Relation_to_Island)
xx=NULL
xx <- merge(as.data.frame(all, row.names = 1),as.data.frame(up,row.names = 1),by=0, all = TRUE)
rownames(xx) <- xx[,1]
xx[,1] = NULL
colnames(xx) <- c("all","up")
xx[is.na(xx)] <- 0
head(xx)
x=xx$up
m=xx$all
n=sum(xx$all)-xx$all
k=sum(xx$up)
xl <- apply(xx,1,function(x) {
mat <- matrix(c(x[2],x[1]-x[2], sum(xx$up)-x[2], sum(xx$all) - sum(xx$up) -x [1] + x[2] ),2,2)
mat
fisher.test(mat)
})
xx$OR <- unname(unlist(lapply(X=xl, FUN = function(x) {x$estimate})))
xx$fisherPval <- unname(unlist(lapply(X=xl, FUN = function(x) {x$p.value})))
xx$lowerCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[1]]})))
xx$upperCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[2]]})))
up_comp <- xx
xx=NULL
xx <- merge(as.data.frame(all, row.names = 1),as.data.frame(dn,row.names = 1),by=0, all = TRUE)
rownames(xx) <- xx[,1]
xx[,1] = NULL
colnames(xx) <- c("all","dn")
xx[is.na(xx)] <- 0
x=xx$dn
m=xx$all
n=sum(xx$all)-xx$all
k=sum(xx$dn)
xl <- apply(xx,1,function(x) {
mat <- matrix(c(x[2],x[1]-x[2], sum(xx$dn)-x[2], sum(xx$all) - sum(xx$dn) -x [1] + x[2] ),2,2)
mat
fisher.test(mat)
})
xx$OR <- unname(unlist(lapply(X=xl, FUN = function(x) {x$estimate})))
xx$fisherPval <- unname(unlist(lapply(X=xl, FUN = function(x) {x$p.value})))
xx$lowerCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[1]]})))
xx$upperCI <- unname(unlist(lapply(X=xl, FUN = function(x) {x$conf.int[[2]]})))
dn_comp <- xx
list("up_comp"=up_comp,"dn_comp"=dn_comp)
}
par(mfrow=c(2,1))
# guthrie
comp <- cgi_enrichment2(dma1)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp## $up_comp
## all up OR fisherPval lowerCI upperCI
## Island 150182 79 0.5877619 3.117728e-06 0.4584982 0.7447366
## OpenSea 442427 372 1.0457548 5.808564e-01 0.8930576 1.2256760
## Shelf 55491 42 0.9121162 6.450546e-01 0.6505930 1.2479057
## Shore 142558 159 1.4667425 5.370682e-05 1.2187456 1.7574979
##
## $dn_comp
## all dn OR fisherPval lowerCI upperCI
## Island 150182 63 0.9426708 0.73255118 0.7057912 1.242379
## OpenSea 442427 178 0.8240650 0.07470790 0.6640092 1.022907
## Shelf 55491 30 1.2499824 0.24707566 0.8294635 1.820371
## Shore 142558 77 1.2918868 0.05071192 0.9895641 1.669752make_forest_plots(comp)
# blood
comp <- cgi_enrichment2(dma2)
comp <- lapply(comp,function(x) {subset(x,OR!=0)} )
comp## $up_comp
## all up OR fisherPval lowerCI upperCI
## Island 150186 86 0.7592429 0.018757791 0.5965829 0.956475
## OpenSea 452360 332 1.0451040 0.614836821 0.8835046 1.237595
## Shelf 56174 31 0.7529520 0.141334722 0.5065357 1.081569
## Shore 143927 129 1.3152090 0.007767201 1.0727433 1.603292
##
## $dn_comp
## all dn OR fisherPval lowerCI upperCI
## Island 150186 55 0.6509374 0.0021787704 0.4810497 0.8660852
## OpenSea 452360 272 1.4044035 0.0008239567 1.1463658 1.7258799
## Shelf 56174 36 1.2395016 0.2142960665 0.8551875 1.7467519
## Shore 143927 59 0.7437787 0.0359749180 0.5549769 0.9816042make_forest_plots(comp)
remove(dma1)
remove(dma2)GMEA functions
This is an enrichment technique.
# get the median, mean and 1-sample t-test result for each gene
pmea <- function(mval,design,sets,cores=2) {
fit <- lmFit(mval, design)
fit <- eBayes(fit)
top <- topTable(fit,coef=ncol(design),num=Inf, sort.by = "P")
l <- mclapply(seq(1,length(sets)), function(i) {
g <- names(sets[i])
tstats <- top[rownames(top) %in% sets[[i]],"t"]
myn <- length(tstats)
mymean <- mean(tstats)
mymedian <- median(tstats)
if ( length(tstats) < 2 ) {
pval=1
} else {
wtselfcont <- t.test(tstats)
pval=wtselfcont$p.value
}
c("gene"=g,"nprobes"=myn,"mean"=mymean,"median"=mymedian,
"P.Value"=pval)
} , mc.cores=cores)
df <- do.call(rbind, l)
rownames(df) <- df[,1]
df <- df[,-1]
tmp <- apply(df,2,as.numeric)
rownames(tmp) <- rownames(df)
df <- as.data.frame(tmp)
df$sig <- -log10(df[,4])
df <- df[order(-df$sig),]
df$FDR <- p.adjust(df$P.Value)
out <- list("df"=df,"toptable"=top)
return(out)
}
# pmea_res <- pmea(mval=mval,design=design,sets=head(sets,20),cores=detectCores()/2)
# Run the fry test for each gene. This is a more conservative test.
run_fry <- function(mval,design,sets,cores=2) {
split_sets <- split(sets, ceiling(seq_along(sets)/200))
fry_l <- mclapply(split_sets,function(l) {
fry(y=mval, index = l, design = design,
contrast = ncol(design) )
} , mc.cores=cores )
fry_res <- do.call(rbind,fry_l)
rownames(fry_res) <- sub("\\.","@",rownames(fry_res))
rownames(fry_res) <- sapply(strsplit(rownames(fry_res),"@"),"[[",2)
fry_res[is.na(fry_res$PValue),"PValue"] <- 1
fry_res <- fry_res[order(fry_res$PValue),]
fry_res$FDR <- p.adjust(fry_res$PValue,method="fdr")
return(fry_res)
}
#fry_res <- run_fry(mval=mval,design=design,sets=sets,cores=cores)
# main function to perform 1-sample t-test and fry test and merge the results.
main <- function(mval,design,sets,cores=2){
pmea_res <- pmea(mval=mval,design=design,sets=sets,cores=cores)
pmea_df <- pmea_res[[1]]
limma_df <- pmea_res[[2]]
fry_res <- run_fry(mval=mval,design=design,sets=sets,cores=cores)
m <- merge(pmea_df,fry_res,by=0)
rownames(m) <- m$Row.names
m$Row.names = NULL
m <- m[,c("nprobes","median","PValue","FDR.y")]
colnames(m) <- c("nprobes","median","PValue","FDR")
m <- m[order(m$PValue),]
out <- list("res"=m,"limma_df"=limma_df)
return(out)
}
#res <- main(mval,design,sets,cores=detectCores()/2)
probe_bias <- function(res) {
res$sig <- -log10(res$PValue)
sig <- subset(res,FDR < 0.05)
plot(res$nprobes,res$sig,log="x",
pch=19,cex=0.6,
xlab="no. probes",ylab="-log10(p-value)")
points(sig$nprobes,sig$sig,col="red",pch=19,cex=0.62)
SIG = nrow(sig)
TOT = nrow(res)
HEADER <- paste(TOT, "genes in total.", SIG, "with FDR<0.05.")
mtext(HEADER)
}
volcano_plot <- function(res) {
res$sig <- -log10(res$PValue)
sig <- subset(res,FDR < 0.05)
plot(res$median,res$sig,
pch=19,cex=0.6,
xlab="median t-statistic",ylab="-log10(p-value)")
points(sig$median,sig$sig,col="red",pch=19,cex=0.62)
SIG = nrow(sig)
UP = nrow(subset(sig,median>0))
DN = nrow(subset(sig,median<0))
TOT = nrow(res)
HEADER <- paste(TOT, "genes in total.", SIG, "with FDR<0.05;",DN,"down,",UP,"up")
mtext(HEADER)
}
gmea_boxplot <- function(res,sets,n=50) {
df <- res[[1]]
limma_df <- res[[2]]
# smallest pval
par(mfrow=c(1,2))
gs <- head(rownames(df),n)
mysets <- sets[names(sets) %in% gs]
tstats <- lapply(mysets, function(set) {
limma_df[rownames(limma_df) %in% set,"t"]
})
tstats <- tstats[order(unlist(lapply(tstats,median)))]
boxplot(tstats,horizontal=TRUE,las=1,
main="smallest p-val",cex.axis=0.6,
xlab="t-statistic")
grid()
# biggest effect size (median)
sig <- subset(df,FDR < 0.05)
gs <- head(rownames(sig[order(-abs(sig$median)),]),n)
if ( length(gs) >2 ) {
tstats <- lapply(gs, function(g) {
df[which(df$genes==g),"tvals"]
})
names(tstats) <- gs
tstats <- tstats[order(unlist(lapply(tstats,median)))]
boxplot(tstats,horizontal=TRUE,las=1,
main="biggest effect size(median)",cex.axis=0.6,
xlab="t-statistic")
grid()
} else {
plot(1)
mtext("too few significant genes found")
}
par(mfrow=c(1,1))
}Session information
For reproducibility
sessionInfo()## R version 4.2.2 Patched (2022-11-10 r83330)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.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] grid stats4 parallel stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] tictoc_1.1
## [2] RIdeogram_0.2.2
## [3] kableExtra_1.3.4
## [4] data.table_1.14.6
## [5] ENmix_1.34.0
## [6] doParallel_1.0.17
## [7] qqman_0.1.8
## [8] RCircos_1.2.2
## [9] beeswarm_0.4.0
## [10] forestplot_3.1.1
## [11] abind_1.4-5
## [12] checkmate_2.1.0
## [13] reshape2_1.4.4
## [14] gplots_3.1.3
## [15] eulerr_7.0.0
## [16] GEOquery_2.66.0
## [17] RColorBrewer_1.1-3
## [18] IlluminaHumanMethylation450kmanifest_0.4.0
## [19] topconfects_1.14.0
## [20] DMRcatedata_2.16.0
## [21] ExperimentHub_2.6.0
## [22] AnnotationHub_3.6.0
## [23] BiocFileCache_2.6.0
## [24] dbplyr_2.3.0
## [25] DMRcate_2.12.0
## [26] limma_3.54.0
## [27] missMethyl_1.32.0
## [28] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1
## [29] R.utils_2.12.2
## [30] R.oo_1.25.0
## [31] R.methodsS3_1.8.2
## [32] plyr_1.8.8
## [33] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [34] minfi_1.44.0
## [35] bumphunter_1.40.0
## [36] locfit_1.5-9.7
## [37] iterators_1.0.14
## [38] foreach_1.5.2
## [39] Biostrings_2.66.0
## [40] XVector_0.38.0
## [41] SummarizedExperiment_1.28.0
## [42] Biobase_2.58.0
## [43] MatrixGenerics_1.10.0
## [44] matrixStats_0.63.0
## [45] GenomicRanges_1.50.2
## [46] GenomeInfoDb_1.34.6
## [47] IRanges_2.32.0
## [48] S4Vectors_0.36.1
## [49] BiocGenerics_0.44.0
## [50] mitch_1.10.0
##
## loaded via a namespace (and not attached):
## [1] rappdirs_0.3.3 rtracklayer_1.58.0
## [3] GGally_2.1.2 tidyr_1.2.1
## [5] ggplot2_3.4.0 bit64_4.0.5
## [7] knitr_1.41 DelayedArray_0.24.0
## [9] rpart_4.1.19 KEGGREST_1.38.0
## [11] RCurl_1.98-1.9 AnnotationFilter_1.22.0
## [13] generics_0.1.3 GenomicFeatures_1.50.3
## [15] preprocessCore_1.60.1 RSQLite_2.2.20
## [17] bit_4.0.5 tzdb_0.3.0
## [19] webshot_0.5.4 xml2_1.3.3
## [21] httpuv_1.6.8 assertthat_0.2.1
## [23] xfun_0.36 hms_1.1.2
## [25] jquerylib_0.1.4 evaluate_0.20
## [27] promises_1.2.0.1 fansi_1.0.3
## [29] restfulr_0.0.15 scrime_1.3.5
## [31] progress_1.2.2 caTools_1.18.2
## [33] readxl_1.4.1 DBI_1.1.3
## [35] geneplotter_1.76.0 htmlwidgets_1.6.1
## [37] reshape_0.8.9 purrr_1.0.1
## [39] ellipsis_0.3.2 dplyr_1.0.10
## [41] backports_1.4.1 permute_0.9-7
## [43] calibrate_1.7.7 grImport2_0.2-0
## [45] annotate_1.76.0 biomaRt_2.54.0
## [47] deldir_1.0-6 sparseMatrixStats_1.10.0
## [49] vctrs_0.5.1 ensembldb_2.22.0
## [51] withr_2.5.0 cachem_1.0.6
## [53] Gviz_1.42.0 BSgenome_1.66.2
## [55] GenomicAlignments_1.34.0 prettyunits_1.1.1
## [57] mclust_6.0.0 svglite_2.1.1
## [59] cluster_2.1.4 RPMM_1.25
## [61] lazyeval_0.2.2 crayon_1.5.2
## [63] genefilter_1.80.2 edgeR_3.40.1
## [65] pkgconfig_2.0.3 nlme_3.1-161
## [67] ProtGenerics_1.30.0 nnet_7.3-18
## [69] rlang_1.0.6 lifecycle_1.0.3
## [71] filelock_1.0.2 dichromat_2.0-0.1
## [73] rsvg_2.4.0 cellranger_1.1.0
## [75] rngtools_1.5.2 base64_2.0.1
## [77] Matrix_1.5-3 Rhdf5lib_1.20.0
## [79] base64enc_0.1-3 viridisLite_0.4.1
## [81] png_0.1-8 rjson_0.2.21
## [83] bitops_1.0-7 KernSmooth_2.23-20
## [85] rhdf5filters_1.10.0 blob_1.2.3
## [87] DelayedMatrixStats_1.20.0 doRNG_1.8.6
## [89] stringr_1.5.0 nor1mix_1.3-0
## [91] readr_2.1.3 jpeg_0.1-10
## [93] scales_1.2.1 memoise_2.0.1
## [95] magrittr_2.0.3 zlibbioc_1.44.0
## [97] compiler_4.2.2 BiocIO_1.8.0
## [99] illuminaio_0.40.0 Rsamtools_2.14.0
## [101] cli_3.6.0 DSS_2.46.0
## [103] htmlTable_2.4.1 Formula_1.2-4
## [105] MASS_7.3-58.1 tidyselect_1.2.0
## [107] stringi_1.7.12 highr_0.10
## [109] yaml_2.3.6 askpass_1.1
## [111] latticeExtra_0.6-30 sass_0.4.4
## [113] VariantAnnotation_1.44.0 tools_4.2.2
## [115] rstudioapi_0.14 foreign_0.8-84
## [117] bsseq_1.34.0 gridExtra_2.3
## [119] digest_0.6.31 BiocManager_1.30.19
## [121] shiny_1.7.4 quadprog_1.5-8
## [123] Rcpp_1.0.9 siggenes_1.72.0
## [125] BiocVersion_3.16.0 later_1.3.0
## [127] org.Hs.eg.db_3.16.0 httr_1.4.4
## [129] AnnotationDbi_1.60.0 biovizBase_1.46.0
## [131] colorspace_2.0-3 rvest_1.0.3
## [133] XML_3.99-0.13 splines_4.2.2
## [135] statmod_1.5.0 multtest_2.54.0
## [137] systemfonts_1.0.4 xtable_1.8-4
## [139] jsonlite_1.8.4 dynamicTreeCut_1.63-1
## [141] R6_2.5.1 echarts4r_0.4.4
## [143] Hmisc_4.7-2 pillar_1.8.1
## [145] htmltools_0.5.4 mime_0.12
## [147] glue_1.6.2 fastmap_1.1.0
## [149] BiocParallel_1.32.5 interactiveDisplayBase_1.36.0
## [151] beanplot_1.3.1 codetools_0.2-18
## [153] utf8_1.2.2 lattice_0.20-45
## [155] bslib_0.4.2 tibble_3.1.8
## [157] curl_5.0.0 gtools_3.9.4
## [159] openssl_2.0.5 interp_1.1-3
## [161] survival_3.5-0 rmarkdown_2.19
## [163] munsell_0.5.0 rhdf5_2.42.0
## [165] GenomeInfoDbData_1.2.9 HDF5Array_1.26.0
## [167] impute_1.72.2 gtable_0.3.1