Correlates of inverse IQ and methylation

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:

limma_blood_ilanguage.csv
limma_guthrie_ilanguage.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")
  library("globaltest")
  library("ebGSEA")
  data("dualmap850kEID")
})

source("meth_functions.R")

Load data

bl_mvals <- readRDS(file="bl_mvals.rds")
gu_mvals <- readRDS(file="gu_mvals.rds")

bl_design <- readRDS(file="bl_design_ilanguage.rds")
gu_design <- readRDS(file="gu_design_ilanguage.rds")

bl_lim <- read.csv("limma_blood_ilanguage.csv")
gu_lim <- read.csv("limma_guthrie_ilanguage.csv")

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.00

Gene level enrichment for blood

Effect size is estimated with the median t-statistic for each gene. Significance is estimated using the fry test. Results are merged for the user.

PMEA step 196.063 sec elapsed on 8 cores
Fry step 293.569 sec elapsed on 8 cores
Fry single core 1205.519 sec

tic()
res_bl <- main(mval=bl_mvals,design=bl_design,sets,cores=detectCores()/2)
res_bl[[1]] <- subset(res_bl[[1]],nprobes!=0) # remove genes with no probes obvs
res_bl_df <- res_bl[[1]] # enrichment results
res_bl_top <- res_bl[[2]] # limma results
toc() ## 502-522 sec elapsed on 8 cores (8.3-8.7 mins)

## 503.424 sec elapsed

head(res_bl_df,30) %>%
  kbl(caption="Top DM genes in blood identified with GMEA") %>%
  kable_paper("hover", full_width = F)

Top DM genes in blood identified with GMEA
	nprobes	median	PValue	FDR
SNORD6	7	-1.9882610	2.10e-06	0.0106485
MIXL1	16	-1.6090148	2.10e-06	0.0106485
AP3S2	33	-0.9182507	2.60e-06	0.0106485
UCK1	13	-1.3850986	2.70e-06	0.0106485
ATL1	48	-0.2831590	3.40e-06	0.0106485
SYCN	12	-2.0614308	5.60e-06	0.0106485
PITX3	31	-1.2894677	6.30e-06	0.0106485
U2SURP	8	-2.1584287	6.90e-06	0.0106485
TMEFF1	14	-1.3394303	6.90e-06	0.0106485
ABO	18	-2.3115247	7.10e-06	0.0106485
FAM162B	16	-1.2203038	7.10e-06	0.0106485
DCLK2	62	-1.1264628	7.80e-06	0.0106485
ACBD3	27	-1.1025352	7.80e-06	0.0106485
SQLE	24	-1.5889093	8.20e-06	0.0106485
EMCN	16	-2.7046781	8.50e-06	0.0106485
MRPL21	26	-0.6758746	9.30e-06	0.0106485
SMIM8	7	-2.6510976	9.90e-06	0.0106485
KPNA1	28	-1.2206313	1.06e-05	0.0106485
SELO	34	-0.9730533	1.09e-05	0.0106485
MARCH8	43	-1.3617938	1.14e-05	0.0106485
MAP2K1	38	-1.0096544	1.15e-05	0.0106485
MRPL55	26	-0.3675478	1.16e-05	0.0106485
SPECC1L	39	-1.4993315	1.16e-05	0.0106485
COL1A2	27	-2.1262718	1.21e-05	0.0106485
VTA1	26	-0.1880468	1.24e-05	0.0106485
CARD17	5	-0.9824821	1.24e-05	0.0106485
COX19	36	-1.4354743	1.31e-05	0.0106485
EYA4	90	-1.0553683	1.33e-05	0.0106485
AK6	15	-1.4297739	1.34e-05	0.0106485
LOC103352541	2	-2.2128620	1.36e-05	0.0106485

par(mfrow=c(2,1))
hist(res_bl_top$t,xlab="probe t-stat",main="blood at assessment") ; grid() ; abline(v=0,lty=1,lwd=3)
hist(res_bl_df$median,xlab="gene median t-stat",main="blood at assessment") ; grid() ; abline(v=0,lty=1,lwd=3)

pdf("hist_bl_ilanguage.pdf")
par(mfrow=c(2,1))
hist(res_bl_top$t,xlab="probe t-stat",main="blood at assessment") ; grid() ; abline(v=0,lty=1,lwd=3)
plot(1)
hist(res_bl_df$median,xlab="gene median t-stat",main="blood at assessment") ; grid() ; abline(v=0,lty=1,lwd=3)
plot(1)
dev.off()

## png 
##   2

par(mfrow=c(1,1))

Now some visualisations.

volcano_plot(res_bl_df)

probe_bias(res_bl_df)

gmea_boxplot(res_bl,sets)

Gene level enrichment for guthrie cards

tic()
res_gu <- main(mval=gu_mvals,design=gu_design,sets,cores=detectCores()/2)
res_gu[[1]] <- subset(res_gu[[1]],nprobes!=0) # remove genes with no probes obvs
res_gu_df <- res_gu[[1]] # enrichment results
res_gu_top <- res_gu[[2]] # limma results
toc() ## 502-522 sec elapsed on 8 cores (8.3-8.7 mins)

## 465.308 sec elapsed

head(res_gu_df,30) %>%
  kbl(caption="Top DM genes in Guthrie card identified with GMEA") %>%
  kable_paper("hover", full_width = F)

Top DM genes in Guthrie card identified with GMEA
	nprobes	median	PValue	FDR
C19orf12	22	-0.4985087	0.0000021	0.0568644
IRAIN	1	-6.3640135	0.0000240	0.3282428
TREML3	2	-2.5712583	0.0000551	0.5026194
LOC401324	3	-2.6940139	0.0001016	0.6953240
FOXD2-AS1	1	-4.3553728	0.0003364	0.9484434
LINC01128	3	-2.3829515	0.0005174	0.9484434
IGFL1	9	-1.4207310	0.0005684	0.9484434
TAS1R1	30	-0.5131941	0.0005828	0.9484434
CDC45L	20	0.8369956	0.0006354	0.9484434
CST2	2	-2.2099663	0.0006831	0.9484434
LOC440600	8	-0.9689782	0.0006925	0.9484434
ASPG	41	-0.6379497	0.0007419	0.9484434
ARHGDIG	15	-0.7880962	0.0009776	0.9484434
C1orf25	19	0.3381455	0.0009878	0.9484434
LGALS3BP	19	-0.6718457	0.0010060	0.9484434
WDR38	11	-0.8231076	0.0010258	0.9484434
MYO15A	66	-0.6493167	0.0010537	0.9484434
LSM14B	17	-0.7444866	0.0010955	0.9484434
RTP2	15	-1.1318082	0.0011933	0.9484434
MIR936	5	-1.2529337	0.0012362	0.9484434
HAR1A	14	-0.9846030	0.0012574	0.9484434
SIGLEC11	10	-0.7773793	0.0012784	0.9484434
MIR6716	2	-2.1719422	0.0013003	0.9484434
GPR144	14	-0.7707788	0.0015848	0.9484434
SLC22A8	20	-1.2120095	0.0016961	0.9484434
USE1	19	-0.4198341	0.0017451	0.9484434
LOC101928766	8	-1.6265892	0.0018903	0.9484434
RXRG	29	-0.7129962	0.0021397	0.9484434
LINC01276	11	-1.2166790	0.0021467	0.9484434
C14orf132	50	-0.6914275	0.0022193	0.9484434

par(mfrow=c(2,1))
hist(res_gu_top$t,xlab="probe t-stat",main="neonatal Guthrie card") ; grid() ; abline(v=0,lty=1,lwd=3)
hist(res_gu_df$median,xlab="gene median t-stat",main="neonatal Guthrie card") ; grid() ; abline(v=0,lty=1,lwd=3)

pdf("hist_gu_ilanguage.pdf")
par(mfrow=c(2,1))
hist(res_gu_top$t,xlab="probe t-stat",main="neonatal Guthrie card") ; grid() ; abline(v=0,lty=1,lwd=3)
plot(1)
hist(res_gu_df$median,xlab="gene median t-stat",main="neonatal Guthrie card") ; grid() ; abline(v=0,lty=1,lwd=3)
plot(1)
dev.off()

## png 
##   2

par(mfrow=c(1,1))

Now some visualisations.

volcano_plot(res_gu_df)

probe_bias(res_gu_df)

gmea_boxplot(res_gu,sets)

Comparison of blood and guthrie card

Now compare blood and guthrie card

res_bl_df$bl <- res_bl_df$median
res_gu_df$gu <- res_gu_df$median
m1 <- merge(res_bl_df,res_gu_df,by=0)
m2 <- m1[,c("Row.names","bl","gu")]
rownames(m2) <- m2$Row.names
m2$Row.names=NULL

plot(m2$bl,m2$gu,xlab="bl",ylab="gu",
  col = rgb(red = 0.6, green = 0.6, blue = 0.6, alpha = 0.5) ,
  pch=19,cex=0.9)
grid()
abline(v=0,lty=1,lwd=3,col="black") ; abline(h=0,lty=1,lwd=3,col="black")
mtext("aggregated t-statistic for each gene (median)")

#plot(rank(m1$bl),rank(m1$gu),xlab="bl rank",ylab="gu rank",cex=0.6,pch=19) ; grid()
rm2 <- apply(m2,2,rank)
mydiff <- apply(m2,2,function(x) { length(which(x<0)) } )
rm3 <- rm2
rm3[,1] <- rm3[,1] - mydiff[1]
rm3[,2] <- rm3[,2] - mydiff[2]
rnk <- rm3
    palette <- colorRampPalette(c("white", "yellow", "orange", "red",
        "darkred", "black"))
    xmin = min(rnk[, 1])
    xmax = max(rnk[, 1])
    ymin = min(rnk[, 2])
    ymax = max(rnk[, 2])
    k <- MASS::kde2d(rnk[, 1], rnk[, 2])
    X_AXIS = paste("Rank in contrast", colnames(rnk)[1])
    Y_AXIS = paste("Rank in contrast", colnames(rnk)[2])
    filled.contour(k, xlim = c(xmin, xmax), ylim = c(ymin, ymax), 
        color.palette = palette, plot.title = {
            abline(v = 0, h = 0, lty = 2, lwd = 2, col = "blue")
            title(main = "Rank-rank plot of all genes", xlab = X_AXIS, 
            ylab = Y_AXIS)
        })

pdf("genedens_ilanguage.pdf")
par(mfrow=c(2,2))
plot(1)
plot(1)

plot(m2$bl,m2$gu,xlab="bl",ylab="gu",
  col = rgb(red = 0.6, green = 0.6, blue = 0.6, alpha = 0.5) ,
  pch=19,cex=0.9)
grid()
abline(v=0, h=0, lty=2, lwd=2, col="blue")
mtext("aggregated t-statistic for each gene (median)")

    filled.contour(k, xlim = c(xmin, xmax), ylim = c(ymin, ymax), 
        color.palette = palette, plot.title = {
            abline(v = 0, h = 0, lty = 2, lwd = 2, col = "blue")
            title(main = "Rank-rank plot of all genes", xlab = X_AXIS, 
            ylab = Y_AXIS)
        })

dev.off()

## png 
##   2

Now we need a table of genes with top scores.

m3 <- m2
m3$med <- apply(rnk,1,median)
mg <- merge(res_gu_df,res_bl_df,by=0)
mg$gu = mg$bl = NULL
mg$med <- m3$med
mg2 <- data.frame(mg,rnk)
mg2$med <- apply(rnk,1,median)
mg2 <- mg2[order(mg2$med),]
mg2 <-subset(mg2,nprobes.x>=5)

head(mg2,20) %>%
  kbl(caption="Genes with coordinated hypomethylation in neonatal Guthrie card and blood at assessment") %>%
  kable_paper("hover", full_width = F)

Genes with coordinated hypomethylation in neonatal Guthrie card and blood at assessment
	Row.names	nprobes.x	median.x	PValue.x	FDR.x	nprobes.y	median.y	PValue.y	FDR.y	med	bl	gu
KERA	KERA	5	-1.442164	0.2406691	0.9484434	5	-4.198001	0.0115255	0.0359133	-20171.5	-20388	-19955
MIR1251	MIR1251	6	-1.566666	0.1250110	0.9484434	6	-3.965967	0.0026541	0.0180924	-20153.0	-20263	-20043
HISPPD1	HISPPD1	5	-1.393458	0.3159473	0.9484434	5	-4.233575	0.0032922	0.0195897	-20151.0	-20397	-19905
LINC00383	LINC00383	6	-1.272004	0.2857798	0.9484434	6	-4.470898	0.0024793	0.0176489	-20118.0	-20478	-19758
LOC102467212	LOC102467212	5	-1.315068	0.2336061	0.9484434	5	-4.256651	0.0076799	0.0286468	-20111.5	-20411	-19812
YSK4	YSK4	5	-1.452411	0.0741538	0.9484434	5	-3.915358	0.0011743	0.0144625	-20100.5	-20240	-19961
MIR490	MIR490	5	-1.676489	0.2033219	0.9484434	5	-3.739026	0.0015671	0.0153564	-20095.0	-20094	-20096
BZW1L1	BZW1L1	5	-1.575017	0.2071137	0.9484434	5	-3.713624	0.0031207	0.0192492	-20059.0	-20070	-20048
OR52J3	OR52J3	5	-1.617323	0.1109103	0.9484434	5	-3.690170	0.0001926	0.0116137	-20057.5	-20043	-20072
HOXA10-AS	HOXA10-AS	11	-1.371958	0.0418893	0.9484434	11	-3.835491	0.0009276	0.0136843	-20031.0	-20182	-19880
KLHL41	KLHL41	6	-1.401598	0.1821345	0.9484434	6	-3.749549	0.0004109	0.0119116	-20011.0	-20110	-19912
LOC101927020	LOC101927020	5	-1.262644	0.3879680	0.9484434	5	-3.999246	0.0121889	0.0371423	-20010.0	-20286	-19734
LOC101926944	LOC101926944	5	-1.391290	0.2190700	0.9484434	5	-3.761719	0.0012037	0.0144625	-20008.0	-20116	-19900
SMIM23	SMIM23	5	-1.190943	0.3617565	0.9484434	5	-4.261368	0.0074248	0.0281989	-19994.0	-20413	-19575
IFI44	IFI44	6	-1.263797	0.2734447	0.9484434	6	-3.860024	0.0022181	0.0171120	-19967.5	-20197	-19738
CLEC1B	CLEC1B	5	-1.153377	0.4763260	0.9484434	5	-4.284786	0.0005164	0.0122195	-19950.5	-20427	-19474
TCHHL1	TCHHL1	5	-1.170389	0.2641570	0.9484434	6	-4.156266	0.0017691	0.0158957	-19946.5	-20368	-19525
MUC7	MUC7	11	-1.382363	0.1264465	0.9484434	11	-3.638143	0.0017684	0.0158957	-19934.0	-19978	-19890
TAS2R8	TAS2R8	5	-1.192002	0.4351908	0.9484434	6	-3.947778	0.0079505	0.0291009	-19915.5	-20253	-19578
ADH4	ADH4	5	-1.116192	0.2308001	0.9484434	6	-4.284673	0.0025514	0.0177930	-19894.5	-20426	-19363

mg2 <- mg2[order(-mg2$med),]

head(mg2,20) %>%
  kbl(caption="Genes with coordinated hypermethylation in neonatal Guthrie card and blood at assessment") %>%
  kable_paper("hover", full_width = F)

Genes with coordinated hypermethylation in neonatal Guthrie card and blood at assessment
	Row.names	nprobes.x	median.x	PValue.x	FDR.x	nprobes.y	median.y	PValue.y	FDR.y	med	bl	gu
LOC401397	LOC401397	5	0.9605015	0.2802885	0.9484434	5	2.771171	0.0044878	0.0222873	5192.0	5063.0	5321.0
BPNT1	BPNT1	15	0.9351762	0.1482424	0.9484434	15	2.630442	0.0007931	0.0132026	5179.0	5045.0	5313.0
POLR2K	POLR2K	12	0.7551125	0.2707521	0.9484434	12	3.513846	0.0010772	0.0140939	5154.0	5108.0	5200.0
LOC100506860	LOC100506860	6	0.9864499	0.3962992	0.9484434	7	2.171558	0.0440996	0.0875212	5152.0	4970.0	5334.0
RGS1	RGS1	5	0.8663298	0.4016069	0.9484434	6	2.043556	0.0486597	0.0938187	5111.5	4943.0	5280.0
LOC101928489	LOC101928489	7	0.9339362	0.0428327	0.9484434	9	1.853009	0.0355676	0.0749805	5105.0	4898.0	5312.0
FLJ35776	FLJ35776	9	0.7464442	0.1083689	0.9484434	9	2.114034	0.0033159	0.0196169	5073.0	4957.0	5189.0
GZMA	GZMA	7	0.8124345	0.0053753	0.9484434	9	1.552742	0.0108655	0.0346636	5009.5	4776.0	5243.0
CLEC2B	CLEC2B	15	0.5951499	0.3019104	0.9484434	20	2.325400	0.0375173	0.0779569	5008.5	5009.0	5008.0
HINT2	HINT2	14	0.5952225	0.1812992	0.9484434	14	2.236879	0.0019749	0.0164264	4998.0	4987.0	5009.0
C6orf150	C6orf150	9	0.7032948	0.0975029	0.9484434	10	1.604680	0.0489260	0.0941830	4974.5	4801.0	5148.0
OR4X1	OR4X1	5	0.5201176	0.3356452	0.9484434	5	3.078235	0.0000400	0.0106485	4967.0	5087.0	4847.0
HIST1H2AD	HIST1H2AD	8	0.9449629	0.1488417	0.9484434	8	1.282174	0.0844084	0.1427182	4963.0	4609.5	5316.5
HIST1H2BF	HIST1H2BF	8	0.9449629	0.1488417	0.9484434	8	1.282174	0.0844084	0.1427182	4963.0	4609.5	5316.5
EVI2A	EVI2A	10	0.6023904	0.0058496	0.9484434	10	1.710295	0.0308407	0.0677994	4933.0	4844.0	5022.0
ARL5A	ARL5A	15	0.7472006	0.0536447	0.9484434	15	1.340718	0.1643645	0.2424228	4926.0	4660.0	5192.0
LRRIQ1	LRRIQ1	18	0.5026435	0.4079856	0.9484434	19	2.666761	0.0021149	0.0167612	4924.5	5050.0	4799.0
TMEM194A	TMEM194A	13	0.5931801	0.2026753	0.9484434	14	1.633018	0.0097268	0.0325834	4909.0	4813.0	5005.0
CSTF2T	CSTF2T	16	0.8742612	0.1752447	0.9484434	16	1.197256	0.0645877	0.1163832	4906.0	4526.0	5286.0
DLGAP1-AS2	DLGAP1-AS2	8	0.8743515	0.0941520	0.9484434	8	1.189417	0.0416177	0.0840774	4903.5	4520.0	5287.0

Gene set level enrichment

genesets <- gmt_import("ReactomePathways.gmt")

cores = 8

Mitch enrichment analysis

Here I’m using the median t-statistic for downstream enrichment. Pathways are from REACTOME and analysis using mitch.

res_bl_df$bl <- res_bl_df$median
res_gu_df$gu <- res_gu_df$median

bl <-  res_bl_df[,"bl",drop=FALSE]
mitch_bl <- mitch_calc(bl,genesets,priority="effect")
head( mitch_bl$enrichment_result,30) %>%
  kbl(caption="BLOOD: Top effect size pathways found with mitch") %>%
  kable_paper("hover", full_width = F)

sig <- subset(mitch_bl$enrichment_result,`p.adjustANOVA`<0.05)
head(sig[order(-abs(sig$s.dist)),],30) %>%
  kbl(caption="BLOOD: Top effect size pathways found with mitch after 1% FDR filtering") %>%
  kable_paper("hover", full_width = F)

gu <-  res_gu_df[,"gu",drop=FALSE]
mitch_gu <- mitch_calc(gu,genesets,priority="effect")
head( mitch_gu$enrichment_result,30) %>%
  kbl(caption="GUTHRIE: Top effect size pathways found with mitch") %>%
  kable_paper("hover", full_width = F)

sig <- subset(mitch_gu$enrichment_result,`p.adjustANOVA`<0.05)
head(sig[order(-abs(sig$s.dist)),],30) %>%
  kbl(caption="GUTHRIE: Top effect size pathways found with mitch after 1% FDR filtering") %>%
  kable_paper("hover", full_width = F)

Mitch joint analysis of blood and guthrie card. As using the median value is apparently the most sensitive, we’ll go with that.

m <- merge(res_bl_df,res_gu_df,by=0)
rownames(m) <- m$Row.names
m2 <- m[,c("bl","gu")]
m2[is.na(m2)] <- 0
res <- mitch_calc(m2, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be 
##             statistically significant.

mitch_report(res,"blgu_mitch.html",overwrite=TRUE)

## Note: overwriting existing report

## Dataset saved as " /tmp/RtmpVDIxly/blgu_mitch.rds ".

## 
## 
## processing file: mitch.Rmd

## 
  |                                               
  |                                         |   0%
  |                                               
  |.                                        |   3%                             
  |                                               
  |..                                       |   6% (checklibraries)            
  |                                               
  |....                                     |   9%                             
  |                                               
  |.....                                    |  12% (peek)                      
  |                                               
  |......                                   |  15%                             
  |                                               
  |.......                                  |  18% (metrics)                   
  |                                               
  |........                                 |  21%                             
  |                                               
  |..........                               |  24% (scatterplot)               
  |                                               
  |...........                              |  26%                             
  |                                               
  |............                             |  29% (contourplot)               
  |                                               
  |.............                            |  32%                             
  |                                               
  |..............                           |  35% (input_geneset_metrics1)    
  |                                               
  |................                         |  38%                             
  |                                               
  |.................                        |  41% (input_geneset_metrics2)    
  |                                               
  |..................                       |  44%                             
  |                                               
  |...................                      |  47% (input_geneset_metrics3)    
  |                                               
  |....................                     |  50%                             
  |                                               
  |......................                   |  53% (echart1d)                  
  |                                               
  |.......................                  |  56% (echart2d)                  
  |                                               
  |........................                 |  59%                             
  |                                               
  |.........................                |  62% (heatmap)                   
  |                                               
  |...........................              |  65%                             
  |                                               
  |............................             |  68% (effectsize)                
  |                                               
  |.............................            |  71%                             
  |                                               
  |..............................           |  74% (results_table)             
  |                                               
  |...............................          |  76%                             
  |                                               
  |.................................        |  79% (results_table_complete)    
  |                                               
  |..................................       |  82%                             
  |                                               
  |...................................      |  85% (detailed_geneset_reports1d)
  |                                               
  |....................................     |  88%                             
  |                                               
  |.....................................    |  91% (detailed_geneset_reports2d)
  |                                               
  |.......................................  |  94%                             
  |                                               
  |........................................ |  97% (session_info)              
  |                                               
  |.........................................| 100%

## output file: /home/mdz/projects/asd_meth/mitch.knit.md

## /usr/bin/pandoc +RTS -K512m -RTS /home/mdz/projects/asd_meth/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpVDIxly/mitch_report.html --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/latex-div.lua --self-contained --variable bs3=TRUE --section-divs --template /usr/local/lib/R/site-library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --mathjax --variable 'mathjax-url=https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' --include-in-header /tmp/RtmpVDIxly/rmarkdown-str3a6a6e6ffec546.html

## 
## Output created: /tmp/RtmpVDIxly/mitch_report.html

## [1] TRUE

mitch_plots(res=res,outfile="blgu_mitch_ilanguage.pdf")

## png 
##   2

sig <- subset(res$enrichment_result,p.adjustMANOVA<0.01)

head(sig[order(-abs(sig$s.dist)),],30) %>%
  kbl(caption="Top pathways found with mitch (effect size after 1% FDR filtering)") %>%
  kable_paper("hover", full_width = F)

Top pathways found with mitch (effect size after 1% FDR filtering)
	set	setSize	pMANOVA	s.bl	s.gu	p.bl	p.gu	s.dist	SD	p.adjustMANOVA
1365	TICAM1-dependent activation of IRF3/IRF7	12	0.0003971	0.5293469	0.5523005	0.0014958	0.0009223	0.7650125	0.0162307	0.0022323
1163	SARS-CoV-1 modulates host translation machinery	34	0.0000000	0.6209671	0.4416586	0.0000000	0.0000083	0.7620121	0.1267902	0.0000000
910	Pexophagy	11	0.0010629	0.5277749	0.5290481	0.0024349	0.0023766	0.7472873	0.0009003	0.0054414
212	Class C/3 (Metabotropic glutamate/pheromone receptors)	39	0.0000000	-0.5027237	-0.5115133	0.0000001	0.0000000	0.7172008	0.0062152	0.0000000
903	Peptide chain elongation	84	0.0000000	0.5807640	0.4018308	0.0000000	0.0000000	0.7062257	0.1265249	0.0000000
441	Folding of actin by CCT/TriC	10	0.0009512	0.6782541	0.1754221	0.0002036	0.3367526	0.7005723	0.3555559	0.0049513
397	Eukaryotic Translation Termination	87	0.0000000	0.5775933	0.3875647	0.0000000	0.0000000	0.6955720	0.1343705	0.0000000
395	Eukaryotic Translation Elongation	88	0.0000000	0.5624193	0.4009282	0.0000000	0.0000000	0.6906945	0.1141914	0.0000000
882	PINK1-PRKN Mediated Mitophagy	21	0.0000165	0.4998656	0.4692390	0.0000730	0.0001967	0.6856026	0.0216563	0.0001210
1496	Viral mRNA Translation	84	0.0000000	0.5609079	0.3860822	0.0000000	0.0000000	0.6809384	0.1236204	0.0000000
1145	Response of EIF2AK4 (GCN2) to amino acid deficiency	95	0.0000000	0.5396613	0.3982391	0.0000000	0.0000000	0.6706927	0.1000006	0.0000000
1208	Selenocysteine synthesis	87	0.0000000	0.5553780	0.3738088	0.0000000	0.0000000	0.6694607	0.1283888	0.0000000
829	Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC)	89	0.0000000	0.5471612	0.3740007	0.0000000	0.0000000	0.6627684	0.1224430	0.0000000
516	Glucuronidation	25	0.0000080	-0.3986159	-0.5062894	0.0005592	0.0000117	0.6443784	0.0761367	0.0000610
1169	SARS-CoV-2 modulates host translation machinery	46	0.0000000	0.5085026	0.3949018	0.0000000	0.0000036	0.6438341	0.0803278	0.0000000
1187	SRP-dependent cotranslational protein targeting to membrane	105	0.0000000	0.5176559	0.3721199	0.0000000	0.0000000	0.6375271	0.1029095	0.0000000
449	Formation of a pool of free 40S subunits	94	0.0000000	0.5166662	0.3543593	0.0000000	0.0000000	0.6265097	0.1147683	0.0000000
532	Golgi Cisternae Pericentriolar Stack Reorganization	14	0.0020052	0.3514790	0.5110581	0.0227726	0.0009286	0.6202563	0.1128394	0.0093373
1172	SCF(Skp2)-mediated degradation of p27/p21	58	0.0000000	0.4597786	0.4143119	0.0000000	0.0000000	0.6189109	0.0321499	0.0000000
493	GTP hydrolysis and joining of the 60S ribosomal subunit	104	0.0000000	0.5110067	0.3306852	0.0000000	0.0000000	0.6086711	0.1275066	0.0000000
565	Hh mutants are degraded by ERAD	54	0.0000000	0.4411985	0.4143596	0.0000000	0.0000001	0.6052686	0.0189780	0.0000000
1363	TICAM1, RIP1-mediated IKK complex recruitment	18	0.0006614	0.4617783	0.3856331	0.0006929	0.0046127	0.6016245	0.0538427	0.0035626
492	GSK3B and BTRC:CUL1-mediated-degradation of NFE2L2	50	0.0000000	0.4424686	0.4070849	0.0000001	0.0000006	0.6012459	0.0250201	0.0000000
828	Nonsense Mediated Decay (NMD) enhanced by the Exon Junction Complex (EJC)	108	0.0000000	0.4851059	0.3421526	0.0000000	0.0000000	0.5936296	0.1010833	0.0000000
830	Nonsense-Mediated Decay (NMD)	108	0.0000000	0.4851059	0.3421526	0.0000000	0.0000000	0.5936296	0.1010833	0.0000000
13	APC/C:Cdc20 mediated degradation of Cyclin B	24	0.0000600	0.3361974	0.4881691	0.0043516	0.0000346	0.5927375	0.1074603	0.0004016
1404	The role of GTSE1 in G2/M progression after G2 checkpoint	58	0.0000000	0.4333958	0.4034549	0.0000000	0.0000001	0.5921214	0.0211714	0.0000000
663	L13a-mediated translational silencing of Ceruloplasmin expression	103	0.0000000	0.4943754	0.3237760	0.0000000	0.0000000	0.5909635	0.1206320	0.0000000
165	Cap-dependent Translation Initiation	111	0.0000000	0.5010715	0.3121452	0.0000000	0.0000000	0.5903451	0.1335911	0.0000000
396	Eukaryotic Translation Initiation	111	0.0000000	0.5010715	0.3121452	0.0000000	0.0000000	0.5903451	0.1335911	0.0000000

Camera enrichment analysis

Camera is designed to control for correlation bias. Using the median value gives better results than the directional significance or the combined metric. Blood dataset. Need to filter by FDR then sort by effect size.

stat <- res_bl_df$median
names(stat) <- rownames(res_bl_df)
stat[is.na(stat)] <- 0

tic()
cres <- cameraPR(statistic=stat, index=genesets, use.ranks = FALSE, inter.gene.cor=0.01, sort = TRUE)
cres <- subset(cres,NGenes>4)
cres$FDR <- p.adjust(cres$PValue,method="fdr")
mymedians <- sapply(1:nrow(cres),function(i) {
  myset <- genesets[[which(names(genesets) == rownames(cres)[i])]]
  mystats <- stat[names(stat) %in% myset]
  median(mystats)
})
cres$median <- mymedians
toc() # 1.0 sec

## 3.596 sec elapsed

sig <- subset(cres,FDR<0.05)

head(sig[order(-abs(sig$median)),],20) %>%
  kbl(caption = "BLOOD: Top effect size pathways found with CAMERA after 5% FDR filtering") %>%
  kable_paper("hover", full_width = F)

BLOOD: Top effect size pathways found with CAMERA after 5% FDR filtering
	NGenes	Direction	PValue	FDR	median
Post-transcriptional silencing by small RNAs	7	Down	0.0038263	0.0486340	-2.5747291
Class C/3 (Metabotropic glutamate/pheromone receptors)	39	Down	0.0000000	0.0000042	-1.7242924
Expression and translocation of olfactory receptors	359	Down	0.0000012	0.0000912	-1.6808688
Olfactory Signaling Pathway	366	Down	0.0000019	0.0001341	-1.6117618
Sensory perception of sweet, bitter, and umami (glutamate) taste	41	Down	0.0000233	0.0012524	-1.4092232
Sensory perception of taste	47	Down	0.0003067	0.0075968	-1.2863848
Sensory Perception	570	Down	0.0000571	0.0025066	-1.2359325
Aspirin ADME	44	Down	0.0008114	0.0158350	-1.1545689
Glucuronidation	25	Down	0.0001057	0.0039143	-1.1497817
Regulation of TP53 Activity through Phosphorylation	88	Up	0.0021269	0.0318544	-0.4823377
Cell Cycle Checkpoints	248	Up	0.0024968	0.0359985	-0.4111076
SARS-CoV-1 Infection	134	Up	0.0023168	0.0339097	-0.4089573
G2/M Checkpoints	129	Up	0.0026179	0.0371897	-0.3977334
Mitotic Anaphase	221	Up	0.0011797	0.0198185	-0.3968072
Mitotic Metaphase and Anaphase	222	Up	0.0010299	0.0179257	-0.3948076
Antigen processing-Cross presentation	100	Up	0.0031954	0.0428711	-0.3908201
SARS-CoV-1 modulates host translation machinery	34	Up	0.0000013	0.0000987	0.3873579
TNFR2 non-canonical NF-kB pathway	94	Up	0.0022917	0.0337978	-0.3859942
UCH proteinases	81	Up	0.0010896	0.0186286	-0.3839069
Separation of Sister Chromatids	161	Up	0.0010889	0.0186286	-0.3812494

nrow(subset(cres,FDR<0.05 & Direction=="Up"))

## [1] 145

nrow(subset(cres,FDR<0.05 & Direction=="Down"))

## [1] 9

cres_bl <- cres

Now drilling down to the top set.

myset <- rownames(head(sig[order(-abs(sig$median)),],1))
myset_members <- genesets[[which(names(genesets) == myset)]]
mystats <- stat[which(names(stat) %in% myset_members)]
myl <- list("allgenes"=stat,myset=mystats)
boxplot(myl,cex=0)

library(vioplot)

## Loading required package: sm

## Package 'sm', version 2.2-5.7: type help(sm) for summary information

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

vioplot(myl,main=myset)

mygenes <- head(mystats[order(mystats)],5)
mygenes

##      AGO1      AGO4      AGO3    TNRC6A    TNRC6C 
## -3.285122 -2.984469 -2.577816 -2.574729 -1.810748

myprobes <- unique(unlist(sets[which(names(sets) %in% names(mygenes))]))
myprobedat <- res_bl_top[which(rownames(res_bl_top) %in% myprobes),]

myprobedat %>%
  kbl(caption = "Probe data for top 5 genes") %>%
  kable_paper("hover", full_width = F)

Probe data for top 5 genes
	logFC	AveExpr	t	P.Value	adj.P.Val	B
cg04525998	-0.6695065	2.2051331	-6.2619018	0.0000290	0.0161108	2.7574772
cg02324079	-0.5424994	2.9446621	-4.8607052	0.0003098	0.0162658	0.5539728
cg15863419	-0.4489136	2.2747176	-4.7957116	0.0003480	0.0162658	0.4444189
cg23873145	-0.9637875	2.4859603	-4.6491424	0.0004534	0.0164874	0.1952065
cg11739924	-0.3919682	4.3163159	-4.4720815	0.0006263	0.0169870	-0.1096477
cg07093496	-0.6183371	3.2706702	-4.4717062	0.0006268	0.0169870	-0.1102981
cg24008619	-0.8388309	3.0434555	-4.3467804	0.0007890	0.0175717	-0.3277312
cg05416591	-0.8681966	3.0121867	-4.2825918	0.0008886	0.0179032	-0.4401537
cg13143099	-0.4843671	3.8957234	-4.2317610	0.0009767	0.0181631	-0.5295020
cg27078729	-0.2996983	3.4065867	-4.2106833	0.0010158	0.0182968	-0.5666318
cg19699973	-0.3565108	3.6719598	-4.1759472	0.0010838	0.0185338	-0.6279221
cg02183851	-0.2759769	3.9710753	-4.1156463	0.0012132	0.0189932	-0.7346068
cg26340657	-0.6167166	3.1259482	-4.1133370	0.0012185	0.0190039	-0.7386993
cg02393980	-0.4691323	3.3251791	-4.0822594	0.0012916	0.0192731	-0.7938250
cg15744368	-0.9274538	2.7917669	-4.0733450	0.0013134	0.0193587	-0.8096541
cg06476845	-0.4287926	4.3039100	-4.0724997	0.0013155	0.0193650	-0.8111555
cg09024808	-0.4040012	4.0126904	-4.0665396	0.0013303	0.0194201	-0.8217430
cg21512216	-0.5290094	3.0663717	-4.0142147	0.0014679	0.0199110	-0.9148300
cg00496251	-0.3922419	4.0874038	-3.9981303	0.0015131	0.0200893	-0.9434920
cg05404890	-0.3830171	4.2705700	-3.9767853	0.0015752	0.0203343	-0.9815616
cg27186436	-0.4491451	3.3707054	-3.9709918	0.0015925	0.0204015	-0.9919009
cg11158641	-0.6055694	3.3392500	-3.8248604	0.0021001	0.0222368	-1.2535351
cg18752890	-0.5793243	3.4698101	-3.8193919	0.0021220	0.0223074	-1.2633549
cg22463563	-0.4456983	4.1707013	-3.7959537	0.0022187	0.0226423	-1.3054637
cg16785092	-0.7078907	3.2401719	-3.7863951	0.0022594	0.0227857	-1.3226462
cg05353134	-0.7673800	3.5450705	-3.7816542	0.0022799	0.0228544	-1.3311706
cg06282219	-0.3860057	4.3852101	-3.7602166	0.0023748	0.0232053	-1.3697323
cg16124834	-0.7577977	3.3265298	-3.7596420	0.0023774	0.0232144	-1.3707663
cg03156546	-0.6219926	0.7655621	-3.7194939	0.0025664	0.0239217	-1.4430540
cg07787614	-0.3017845	-0.0859365	-3.7138235	0.0025943	0.0240227	-1.4532705
cg03664025	-0.5650885	3.6750238	-3.6484367	0.0029393	0.0253324	-1.5711891
cg01303655	-0.6364708	3.2961779	-3.5728679	0.0033965	0.0270207	-1.7076851
cg18748888	-0.4477384	4.6275557	-3.5658700	0.0034424	0.0271833	-1.7203346
cg25644447	0.4545088	-3.8647746	3.5296147	0.0036901	0.0280771	-1.7858920
cg00018443	-0.4273985	4.1250086	-3.5268472	0.0037097	0.0281479	-1.7908976
cg09816107	-0.4363462	3.8036341	-3.5167630	0.0037822	0.0284102	-1.8091382
cg04215186	-0.3572516	4.3731017	-3.5142105	0.0038007	0.0284795	-1.8137556
cg00152838	0.4055649	-3.2941953	3.4099238	0.0046435	0.0315083	-2.0024928
cg04256065	-0.6784260	4.9860043	-3.3844855	0.0048763	0.0323123	-2.0485437
cg25362361	-0.3329820	4.3100507	-3.3811627	0.0049076	0.0324129	-2.0545590
cg09889256	-0.4377801	4.1830191	-3.3334738	0.0053792	0.0340459	-2.1408852
cg01329260	-0.2454907	-1.7412320	-3.2851222	0.0059039	0.0358583	-2.2283871
cg14781700	-0.6158522	3.5650395	-3.2793078	0.0059704	0.0360832	-2.2389069
cg12489727	-0.5166535	3.7335510	-3.2619872	0.0061729	0.0367487	-2.2702399
cg20035319	-0.3875523	3.5427599	-3.2256520	0.0066204	0.0382595	-2.3359462
cg02391633	-0.3269045	4.1732022	-3.1884739	0.0071119	0.0398784	-2.4031344
cg07530460	-0.1872637	4.2125512	-3.1593226	0.0075228	0.0412098	-2.4557799
cg21498511	-0.8152181	1.3841073	-3.0362037	0.0095359	0.0474423	-2.6776452
cg17525025	0.2920549	-4.5574238	3.0027304	0.0101705	0.0493206	-2.7377956
cg07157902	-0.2515190	3.0134901	-2.9430522	0.0114076	0.0528780	-2.8448091
cg18045878	-0.4974457	1.5292289	-2.9395362	0.0114849	0.0531025	-2.8511040
cg02785870	-0.3151234	4.4977700	-2.9103544	0.0121474	0.0549721	-2.9033040
cg26200971	-0.2366234	4.2669840	-2.8293820	0.0141897	0.0606864	-3.0476736
cg14619832	-0.3001402	4.2596628	-2.7925788	0.0152263	0.0634480	-3.1130348
cg07768777	-0.2482032	4.1245728	-2.7565339	0.0163136	0.0662926	-3.1768758
cg15473007	-0.3304667	0.8130780	-2.7538733	0.0163968	0.0665004	-3.1815809
cg05873271	-0.1945508	-1.2457826	-2.7488219	0.0165559	0.0669027	-3.1905114
cg09589586	-0.2520902	4.2032768	-2.7124034	0.0177487	0.0699907	-3.2547849
cg04268643	-0.2089903	-0.1359135	-2.7044724	0.0180194	0.0706668	-3.2687554
cg18833170	-0.1660208	4.2577083	-2.6988210	0.0182148	0.0711449	-3.2787043
cg11696355	-0.2483818	4.0212366	-2.6896295	0.0185369	0.0719489	-3.2948744
cg01560476	0.2351854	1.0426569	2.6703051	0.0192325	0.0736429	-3.3288265
cg14908168	-0.2549218	3.5559920	-2.5747291	0.0230635	0.0828902	-3.4957936
cg07710734	0.1922728	3.4265711	2.5415460	0.0245594	0.0863448	-3.5533599
cg00423030	-0.1890065	0.0994404	-2.4420327	0.0296283	0.0977030	-3.7245943
cg09029192	-0.2064013	1.5186864	-2.3939405	0.0324248	0.1036938	-3.8065252
cg02958960	-0.1746421	3.7728018	-2.3538315	0.0349485	0.1089666	-3.8744090
cg12700283	-0.1624570	0.1119807	-2.3106033	0.0378776	0.1149705	-3.9470906
cg27407707	0.2213543	-4.1243624	2.1937469	0.0470045	0.1328068	-4.1408463
cg16791170	-0.1374864	4.0862128	-2.1194274	0.0538454	0.1453663	-4.2618059
cg09146982	-0.1554287	0.5138050	-2.0813172	0.0577035	0.1522741	-4.3230879
cg00672192	0.1414764	0.9711645	2.0638595	0.0595552	0.1555467	-4.3509837
cg08844988	0.2835730	-3.8678881	2.0229884	0.0641077	0.1634583	-4.4158429
cg03660034	-0.1416407	3.0410668	-1.9405154	0.0742847	0.1803717	-4.5447116
cg08898446	-0.1300190	4.0241655	-1.9384527	0.0745572	0.1808302	-4.5478988
cg24207009	-0.2024391	-0.5670785	-1.9300394	0.0756782	0.1826237	-4.5608800
cg10134076	-0.1776110	4.3350578	-1.8585948	0.0858334	0.1987661	-4.6698721
cg01880908	-0.1873937	5.5418080	-1.8234163	0.0912731	0.2071395	-4.7226924
cg06420129	-0.1271357	3.3301889	-1.8204965	0.0917382	0.2078366	-4.7270506
cg12069304	-0.1350360	3.7279495	-1.8120408	0.0930972	0.2099609	-4.7396489
cg05478790	-0.1262252	2.3018531	-1.8107477	0.0933067	0.2102819	-4.7415726
cg04762133	0.1483586	3.0203085	1.7103734	0.1109036	0.2361423	-4.8883912
cg03760856	-0.1524466	4.3306718	-1.6679057	0.1191923	0.2478494	-4.9489617
cg21571594	0.1047939	3.7478583	1.6176545	0.1296966	0.2624011	-5.0193810
cg00150128	-0.1020941	1.5623377	-1.5974342	0.1341450	0.2684354	-5.0473218
cg27579805	0.1021259	1.9909200	1.5918128	0.1354049	0.2700968	-5.0550487
cg12963342	-0.1287899	4.1860140	-1.5458737	0.1460887	0.2842685	-5.1175120
cg05890693	-0.0942396	-1.1002304	-1.5016208	0.1570541	0.2984775	-5.1765073
cg16792948	0.1025535	-5.1924596	1.4770478	0.1634390	0.3065416	-5.2087545
cg14744741	-0.1383403	4.2764331	-1.4763657	0.1636193	0.3067736	-5.2096444
cg04498014	-0.1168476	0.1044817	-1.3789110	0.1911516	0.3405091	-5.3337477
cg26740249	-0.0860303	0.6715219	-1.3292028	0.2066037	0.3586301	-5.3946540
cg08446633	-0.0955655	4.1630504	-1.3173185	0.2104447	0.3630025	-5.4089675
cg25012186	-0.1487844	3.3153320	-1.3165556	0.2106932	0.3632768	-5.4098830
cg01053405	-0.1281078	3.4620805	-1.3021381	0.2154347	0.3687872	-5.4271095
cg00472027	-0.1156890	4.0733751	-1.2872626	0.2204165	0.3744855	-5.4447317
cg08500967	-0.1227000	-4.4052158	-1.2141356	0.2462599	0.4035178	-5.5290800
cg05689121	0.1147050	-4.6870701	1.1494938	0.2710307	0.4302352	-5.6003878
cg15831500	-0.0873823	4.0140877	-1.0843509	0.2978809	0.4583332	-5.6690409
cg13800273	-0.0973113	3.9304125	-1.0528284	0.3115676	0.4722478	-5.7010726
cg18878381	0.0648999	3.8046251	1.0100753	0.3308633	0.4914766	-5.7432471
cg08076643	-0.0738249	0.1244664	-0.9334574	0.3675729	0.5269123	-5.8150834
cg06037603	0.1248534	3.2685268	0.8800858	0.3947676	0.5526268	-5.8622076
cg08083403	0.0687916	-4.8491599	0.8594170	0.4056566	0.5627234	-5.8797993
cg23972793	-0.0692379	3.7567842	-0.8369091	0.4177408	0.5738424	-5.8985334
cg12958612	0.0515338	-5.0038565	0.8104116	0.4322677	0.5868372	-5.9200173
cg01155522	0.0538019	-4.6516656	0.7979305	0.4392225	0.5930934	-5.9299213
cg12027254	-0.0814866	1.8287541	-0.7495008	0.4668827	0.6173693	-5.9670300
cg27503577	-0.0549950	4.0942280	-0.6327310	0.5378640	0.6777133	-6.0476783
ch.16.759907F	0.0520217	-4.1980441	0.6038856	0.5562959	0.6926339	-6.0656377
cg04508233	0.0356501	-6.3243602	0.5942408	0.5625350	0.6975714	-6.0714664
cg04047827	-0.0366427	-4.9183537	-0.4938947	0.6296079	0.7505222	-6.1268010
cg12072376	0.0503900	3.6848109	0.4580591	0.6544613	0.7692745	-6.1441859
cg12229081	-0.0405821	-3.5553608	-0.4331645	0.6719850	0.7824518	-6.1555175
cg13280512	-0.0508799	1.9516830	-0.4103848	0.6881969	0.7945286	-6.1653477
cg17172331	-0.0218578	0.3146128	-0.3865893	0.7053045	0.8069989	-6.1750638
cg00175441	-0.0327317	3.9293636	-0.3692814	0.7178543	0.8159564	-6.1817750
cg09657673	-0.0297928	3.2897807	-0.3552005	0.7281274	0.8233196	-6.1870129
cg15661240	-0.0287192	3.3347891	-0.3524342	0.7301522	0.8247530	-6.1880185
cg08073142	-0.0246797	-6.5432840	-0.3482911	0.7331887	0.8269337	-6.1895102
cg17969276	-0.0404229	-4.5513599	-0.3211742	0.7531758	0.8410786	-6.1988454
cg22094953	-0.0190553	4.0095283	-0.2871631	0.7785068	0.8587384	-6.2095012
cg15489593	-0.0362374	-3.1739684	-0.2798278	0.7840060	0.8625615	-6.2116453
cg20857275	-0.0219502	0.1503920	-0.2416258	0.8128340	0.8821357	-6.2219245
cg14230917	-0.0118790	-4.8081088	-0.1918441	0.8508240	0.9074510	-6.2330779
cg00382632	-0.0147249	5.3193775	-0.1788270	0.8608263	0.9139885	-6.2355744
cg06223856	-0.0304672	0.8331048	-0.1716157	0.8663783	0.9175856	-6.2368824
cg09646181	0.0076717	3.5037929	0.1162083	0.9092608	0.9448257	-6.2451422
cg01436550	0.0130110	3.5053408	0.1144580	0.9106211	0.9456974	-6.2453514
cg10943443	0.0044324	3.6927332	0.0577953	0.9547893	0.9729798	-6.2504120
cg05736847	-0.0020957	3.6288991	-0.0270197	0.9788538	0.9874880	-6.2517660

Guthrie dataset.

stat <- res_gu_df$median
names(stat) <- rownames(res_gu_df)
stat[is.na(stat)] <- 0

tic()
cres <- cameraPR(statistic=stat, index=genesets, use.ranks = FALSE, inter.gene.cor=0.01, sort = TRUE)
cres <- subset(cres,NGenes>4)
cres$FDR <- p.adjust(cres$PValue,method="fdr")
mymedians <- sapply(1:nrow(cres),function(i) {
  myset <- genesets[[which(names(genesets) == rownames(cres)[i])]]
  mystats <- stat[names(stat) %in% myset]
  median(mystats)
})
cres$median <- mymedians
toc() # 1.0 sec

## 2.863 sec elapsed

sig <- subset(cres,FDR<0.05)

head(sig[order(-abs(sig$median)),],20) %>%
  kbl(caption = "GUTHRIE: Top effect size pathways found with CAMERA after 5% FDR filtering") %>%
  kable_paper("hover", full_width = F)

GUTHRIE: Top effect size pathways found with CAMERA after 5% FDR filtering
	NGenes	Direction	PValue	FDR	median
Expression and translocation of olfactory receptors	356	Down	0.0000000	0.0000000	-0.7322374
Olfactory Signaling Pathway	363	Down	0.0000000	0.0000000	-0.7313853
Aspirin ADME	44	Down	0.0000156	0.0041099	-0.6863595
Class C/3 (Metabotropic glutamate/pheromone receptors)	39	Down	0.0000021	0.0010364	-0.6851519
Glucuronidation	25	Down	0.0002288	0.0094041	-0.6840526
Keratinization	213	Down	0.0000102	0.0039535	-0.6022933
Sensory Perception	567	Down	0.0000000	0.0000020	-0.5774976
Sensory perception of sweet, bitter, and umami (glutamate) taste	41	Down	0.0001855	0.0088423	-0.5262750
Sensory perception of taste	47	Down	0.0001922	0.0088423	-0.5262750
Activation of Matrix Metalloproteinases	31	Down	0.0021035	0.0315035	-0.4875238
Gene expression (Transcription)	1420	Up	0.0036337	0.0452926	-0.1807918
G2/M Transition	175	Up	0.0026165	0.0368989	-0.1755054
Mitotic G2-G2/M phases	177	Up	0.0026393	0.0369501	-0.1755054
Cell Cycle	601	Up	0.0018695	0.0296036	-0.1677739
Transcriptional Regulation by TP53	346	Up	0.0028613	0.0389300	-0.1674637
Late Phase of HIV Life Cycle	126	Up	0.0020618	0.0313315	-0.1641372
Cell Cycle, Mitotic	476	Up	0.0014178	0.0242407	-0.1632072
M Phase	335	Up	0.0022860	0.0339727	-0.1623486
Signaling by ROBO receptors	203	Up	0.0012825	0.0225257	-0.1515047
Cellular responses to stimuli	696	Up	0.0009957	0.0204917	-0.1496043

cres_gu <- cres

Compare CAMERA data

Compare median set score for both contrasts. In the graph below, the larger points represent larger gene sets. Red fill means significant in blood contrast. Blue ring means significant in guthrie contrast.

m4 <- merge(cres_bl,cres_gu,by=0)
sigx <- subset(m4,FDR.x<0.05)
sigy <- subset(m4,FDR.y<0.05)

plot(m4$median.x,m4$median.y,col="gray",pch=19,cex=log10(m4$NGenes.x),
  xlab="blood",ylab="Guthrie",main="median t-stat")
grid()
points(sigx$median.x,sigx$median.y,col="red",pch=19,cex=log10(sigx$NGenes.x))
points(sigy$median.x,sigy$median.y,col="blue",pch=1,cex=log10(sigy$NGenes.x))
abline(v=0,h=0,lty=2,lwd=3)

par(mfrow=c(1,1))
pdf("camera_pw_ilanguage.pdf")

plot(m4$median.x,m4$median.y,col="gray",pch=19,cex=log10(m4$NGenes.x),
  xlab="blood",ylab="Guthrie",main="median t-stat")
grid()
points(sigx$median.x,sigx$median.y,col="red",pch=19,cex=log10(sigx$NGenes.x))
points(sigy$median.x,sigy$median.y,col="blue",pch=1,cex=log10(sigy$NGenes.x))
abline(v=0,h=0,lty=2,lwd=3)

dev.off()

## png 
##   2

subset(m4,FDR.x<0.05 & FDR.y<0.05 ) %>%
  kbl(caption = "Pathways significant in blood (x) and guthrie cards (y) at 5%FDR") %>%
  kable_paper("hover", full_width = F)

Pathways significant in blood (x) and guthrie cards (y) at 5%FDR
	Row.names	NGenes.x	Direction.x	PValue.x	FDR.x	median.x	NGenes.y	Direction.y	PValue.y	FDR.y	median.y
30	Activation of APC/C and APC/C:Cdc20 mediated degradation of mitotic proteins	75	Up	0.0000881	0.0034723	-0.1020616	75	Up	0.0001223	0.0074514	-0.0782185
44	Activation of NF-kappaB in B cells	64	Up	0.0004799	0.0106577	-0.1080817	64	Up	0.0017341	0.0279196	-0.0977134
101	Antigen processing-Cross presentation	100	Up	0.0031954	0.0428711	-0.3908201	100	Up	0.0034851	0.0442973	-0.1142653
107	APC:Cdc20 mediated degradation of cell cycle proteins prior to satisfation of the cell cycle checkpoint	72	Up	0.0001479	0.0044227	-0.0823142	72	Up	0.0001475	0.0074986	-0.0708212
108	APC/C-mediated degradation of cell cycle proteins	86	Up	0.0001331	0.0043159	-0.1616983	86	Up	0.0001404	0.0074514	-0.0808179
110	APC/C:Cdc20 mediated degradation of mitotic proteins	74	Up	0.0000951	0.0036735	-0.0823142	74	Up	0.0001039	0.0074514	-0.0708212
111	APC/C:Cdc20 mediated degradation of Securin	66	Up	0.0001520	0.0044227	-0.0614369	66	Up	0.0001267	0.0074514	-0.0595988
112	APC/C:Cdh1 mediated degradation of Cdc20 and other APC/C:Cdh1 targeted proteins in late mitosis/early G1	72	Up	0.0001534	0.0044227	-0.0823142	72	Up	0.0001427	0.0074514	-0.0755749
128	Aspirin ADME	44	Down	0.0008114	0.0158350	-1.1545689	44	Down	0.0000156	0.0041099	-0.6863595
134	Assembly of the pre-replicative complex	82	Up	0.0021045	0.0317653	-0.1585581	82	Up	0.0002203	0.0092515	-0.0708212
136	Asymmetric localization of PCP proteins	62	Up	0.0009039	0.0166326	-0.1194611	62	Up	0.0011428	0.0216971	-0.0916172
144	AUF1 (hnRNP D0) binds and destabilizes mRNA	53	Up	0.0000843	0.0033926	0.0178378	53	Up	0.0009970	0.0204917	-0.0782185
146	Autodegradation of Cdh1 by Cdh1:APC/C	62	Up	0.0002957	0.0075176	-0.0479101	62	Up	0.0000710	0.0071176	-0.0489751
147	Autodegradation of the E3 ubiquitin ligase COP1	50	Up	0.0001255	0.0042858	-0.0479101	50	Up	0.0011953	0.0219794	-0.0804587
198	Cap-dependent Translation Initiation	111	Up	0.0000000	0.0000057	0.0919582	111	Up	0.0005492	0.0141466	-0.1381658
217	Cdc20:Phospho-APC/C mediated degradation of Cyclin A	71	Up	0.0001364	0.0043159	-0.0625668	71	Up	0.0001025	0.0074514	-0.0634240
220	CDK-mediated phosphorylation and removal of Cdc6	71	Up	0.0003032	0.0075968	-0.0982138	71	Up	0.0001230	0.0074514	-0.0557737
222	Cell Cycle Checkpoints	248	Up	0.0024968	0.0359985	-0.4111076	248	Up	0.0023082	0.0340409	-0.1446161
234	Cellular response to hypoxia	71	Up	0.0009140	0.0166593	-0.3179755	71	Up	0.0016747	0.0274196	-0.1024464
235	Cellular response to starvation	147	Up	0.0000225	0.0012524	-0.1263260	147	Up	0.0000637	0.0071176	-0.0722396
262	Class C/3 (Metabotropic glutamate/pheromone receptors)	39	Down	0.0000000	0.0000042	-1.7242924	39	Down	0.0000021	0.0010364	-0.6851519
324	Cyclin A:Cdk2-associated events at S phase entry	83	Up	0.0009961	0.0174955	-0.1757481	83	Up	0.0006873	0.0166322	-0.0869686
327	Cyclin E associated events during G1/S transition	81	Up	0.0007965	0.0157017	-0.1517451	81	Up	0.0011028	0.0215220	-0.0872553
347	Dectin-1 mediated noncanonical NF-kB signaling	60	Up	0.0006392	0.0130729	-0.2213814	60	Up	0.0009905	0.0204917	-0.0827369
356	Defective CFTR causes cystic fibrosis	59	Up	0.0000360	0.0016941	-0.0355133	59	Up	0.0006899	0.0166322	-0.0634240
387	Degradation of AXIN	53	Up	0.0009834	0.0174301	-0.1757481	53	Up	0.0024930	0.0356778	-0.0872553
390	Degradation of DVL	55	Up	0.0003387	0.0079805	-0.1278626	55	Up	0.0012447	0.0220623	-0.0782185
391	Degradation of GLI1 by the proteasome	58	Up	0.0013652	0.0223518	-0.0826639	58	Up	0.0036074	0.0452610	-0.0916172
392	Degradation of GLI2 by the proteasome	58	Up	0.0009325	0.0168375	-0.1392545	58	Up	0.0034710	0.0442973	-0.0916172
445	DNA Replication	125	Up	0.0007450	0.0149935	-0.2608607	125	Up	0.0003276	0.0115052	-0.0848001
447	DNA Replication Pre-Initiation	97	Up	0.0011251	0.0190678	-0.1757481	97	Up	0.0003255	0.0115052	-0.0848001
457	Downstream signaling events of B Cell Receptor (BCR)	78	Up	0.0013183	0.0217685	-0.3369284	78	Up	0.0034416	0.0442973	-0.1346760
499	ER-Phagosome pathway	87	Up	0.0029766	0.0407860	-0.3070056	87	Up	0.0014617	0.0247724	-0.0994477
518	Eukaryotic Translation Elongation	88	Up	0.0000000	0.0000041	0.2406206	88	Up	0.0000436	0.0064415	-0.0173977
519	Eukaryotic Translation Initiation	111	Up	0.0000000	0.0000057	0.0919582	111	Up	0.0005492	0.0141466	-0.1381658
520	Eukaryotic Translation Termination	87	Up	0.0000000	0.0000037	0.2431649	87	Up	0.0000776	0.0071423	-0.0508882
522	Expression and translocation of olfactory receptors	359	Down	0.0000012	0.0000912	-1.6808688	356	Down	0.0000000	0.0000000	-0.7322374
533	FBXL7 down-regulates AURKA during mitotic entry and in early mitosis	53	Up	0.0001758	0.0048528	-0.0355133	53	Up	0.0008785	0.0196655	-0.0782185
538	FCERI mediated NF-kB activation	74	Up	0.0015663	0.0254289	-0.3217050	74	Up	0.0023921	0.0344889	-0.1111544
566	Formation of a pool of free 40S subunits	94	Up	0.0000001	0.0000070	0.1483608	94	Up	0.0001921	0.0088423	-0.0714718
613	G1/S DNA Damage Checkpoints	66	Up	0.0001722	0.0048204	-0.1815829	66	Up	0.0003251	0.0115052	-0.0619587
614	G1/S Transition	126	Up	0.0013134	0.0217685	-0.2423762	126	Up	0.0008100	0.0186302	-0.1045016
617	G2/M Checkpoints	129	Up	0.0026179	0.0371897	-0.3977334	129	Up	0.0018715	0.0296036	-0.1106350
644	GLI3 is processed to GLI3R by the proteasome	58	Up	0.0008913	0.0165580	-0.1392545	58	Up	0.0029269	0.0395443	-0.0916172
652	Glucuronidation	25	Down	0.0001057	0.0039143	-1.1497817	25	Down	0.0002288	0.0094041	-0.6840526
682	GSK3B and BTRC:CUL1-mediated-degradation of NFE2L2	50	Up	0.0001496	0.0044227	0.0191124	50	Up	0.0004063	0.0120775	-0.0595988
683	GTP hydrolysis and joining of the 60S ribosomal subunit	104	Up	0.0000000	0.0000057	0.1140907	104	Up	0.0003390	0.0115052	-0.1131203
699	Hedgehog ligand biogenesis	63	Up	0.0005840	0.0125370	-0.2670147	63	Up	0.0008856	0.0196655	-0.0657709
705	Hh mutants abrogate ligand secretion	57	Up	0.0001264	0.0042858	-0.0625668	57	Up	0.0003508	0.0115052	-0.0580928
706	Hh mutants are degraded by ERAD	54	Up	0.0000731	0.0030689	-0.0479101	54	Up	0.0002967	0.0114663	-0.0569332
712	HIV Infection	216	Up	0.0005068	0.0111276	-0.3633244	216	Up	0.0003425	0.0115052	-0.1213172
713	HIV Life Cycle	139	Up	0.0019147	0.0295942	-0.3748624	139	Up	0.0011455	0.0216971	-0.1274081
719	Host Interactions of HIV factors	122	Up	0.0020441	0.0313436	-0.3441441	122	Up	0.0004930	0.0136077	-0.0969390
747	Influenza Infection	148	Up	0.0000003	0.0000274	-0.0020937	148	Up	0.0001405	0.0074514	-0.1024548
748	Influenza Viral RNA Transcription and Replication	129	Up	0.0000004	0.0000331	0.0306721	129	Up	0.0000956	0.0074514	-0.0741485
838	L13a-mediated translational silencing of Ceruloplasmin expression	103	Up	0.0000001	0.0000114	0.0919582	103	Up	0.0004346	0.0125316	-0.1047502
843	Late Phase of HIV Life Cycle	126	Up	0.0036443	0.0469385	-0.3633244	126	Up	0.0020618	0.0313315	-0.1641372
874	Major pathway of rRNA processing in the nucleolus and cytosol	171	Up	0.0000000	0.0000042	0.0306721	171	Up	0.0000443	0.0064415	-0.1027438
922	Metabolism of polyamines	55	Up	0.0006319	0.0130729	-0.1757481	55	Up	0.0036078	0.0452610	-0.0872553
925	Metabolism of RNA	639	Up	0.0000019	0.0001341	-0.2900232	639	Up	0.0000378	0.0064415	-0.1324440
952	Mitochondrial translation	93	Up	0.0000161	0.0010028	-0.1164349	93	Up	0.0009846	0.0204917	-0.1081294
953	Mitochondrial translation elongation	87	Up	0.0000299	0.0015181	-0.1677323	87	Up	0.0020679	0.0313315	-0.1273132
954	Mitochondrial translation initiation	87	Up	0.0000210	0.0012308	-0.1677323	87	Up	0.0009283	0.0202055	-0.1081294
955	Mitochondrial translation termination	87	Up	0.0000168	0.0010141	-0.1309188	87	Up	0.0009608	0.0204917	-0.0981363
959	Mitotic Anaphase	221	Up	0.0011797	0.0198185	-0.3968072	221	Up	0.0016972	0.0275551	-0.1289578
960	Mitotic G1 phase and G1/S transition	144	Up	0.0016354	0.0256875	-0.3062859	144	Up	0.0007754	0.0180499	-0.1297066
962	Mitotic Metaphase and Anaphase	222	Up	0.0010299	0.0179257	-0.3948076	222	Up	0.0016255	0.0268412	-0.1281830
976	mRNA Splicing	184	Up	0.0000124	0.0007957	-0.2965636	184	Up	0.0000145	0.0041099	-0.0971896
977	mRNA Splicing - Major Pathway	174	Up	0.0000309	0.0015310	-0.3214684	174	Up	0.0000557	0.0069965	-0.1003278
978	mRNA Splicing - Minor Pathway	51	Up	0.0000606	0.0026023	-0.2915931	51	Up	0.0001529	0.0075725	-0.0854447
1017	Negative regulation of NOTCH4 signaling	52	Up	0.0003219	0.0077737	-0.0088378	52	Up	0.0004141	0.0121215	-0.0595988
1042	NIK–>noncanonical NF-kB signaling	57	Up	0.0003659	0.0085171	-0.1248621	57	Up	0.0010206	0.0205390	-0.0872553
1048	Nonsense Mediated Decay (NMD) enhanced by the Exon Junction Complex (EJC)	108	Up	0.0000001	0.0000141	0.0623098	108	Up	0.0002162	0.0092515	-0.0731940
1049	Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC)	89	Up	0.0000000	0.0000046	0.2056701	89	Up	0.0001150	0.0074514	-0.0688872
1050	Nonsense-Mediated Decay (NMD)	108	Up	0.0000001	0.0000141	0.0623098	108	Up	0.0002162	0.0092515	-0.0731940
1051	NoRC negatively regulates rRNA expression	43	Up	0.0038680	0.0486491	-0.3517863	43	Up	0.0028043	0.0386987	-0.1102945
1096	Olfactory Signaling Pathway	366	Down	0.0000019	0.0001341	-1.6117618	363	Down	0.0000000	0.0000000	-0.7313853
1101	Orc1 removal from chromatin	69	Up	0.0029715	0.0407860	-0.1757481	69	Up	0.0012216	0.0220579	-0.0959791
1112	Oxygen-dependent proline hydroxylation of Hypoxia-inducible Factor Alpha	62	Up	0.0007923	0.0157017	-0.2869486	62	Up	0.0019000	0.0296036	-0.0977134
1116	p53-Dependent G1 DNA Damage Response	64	Up	0.0001385	0.0043159	-0.1191575	64	Up	0.0003514	0.0115052	-0.0619587
1117	p53-Dependent G1/S DNA damage checkpoint	64	Up	0.0001385	0.0043159	-0.1191575	64	Up	0.0003514	0.0115052	-0.0619587
1118	p53-Independent DNA Damage Response	50	Up	0.0002215	0.0058628	-0.0614369	50	Up	0.0003912	0.0119448	-0.0595988
1119	p53-Independent G1/S DNA damage checkpoint	50	Up	0.0002215	0.0058628	-0.0614369	50	Up	0.0003912	0.0119448	-0.0595988
1137	Peptide chain elongation	84	Up	0.0000000	0.0000037	0.2446386	84	Up	0.0000467	0.0064415	-0.0173977
1178	PINK1-PRKN Mediated Mitophagy	21	Up	0.0016209	0.0256690	-0.0886773	21	Up	0.0029589	0.0396982	0.1153080
1230	Processing of Capped Intron-Containing Pre-mRNA	232	Up	0.0000271	0.0014141	-0.2976054	232	Up	0.0000170	0.0041099	-0.0982329
1296	Regulation of activated PAK-2p34 by proteasome mediated degradation	48	Up	0.0002736	0.0070483	-0.0479101	48	Up	0.0010484	0.0206856	-0.0708212
1297	Regulation of APC/C activators between G1/S and early anaphase	79	Up	0.0002527	0.0065967	-0.1476485	79	Up	0.0001157	0.0074514	-0.0782185
1298	Regulation of Apoptosis	51	Up	0.0009552	0.0170874	-0.0625668	51	Up	0.0026132	0.0368989	-0.0855276
1305	Regulation of expression of SLITs and ROBOs	159	Up	0.0000002	0.0000209	-0.0331799	159	Up	0.0000813	0.0071435	-0.0782185
1326	Regulation of mitotic cell cycle	86	Up	0.0001331	0.0043159	-0.1616983	86	Up	0.0001404	0.0074514	-0.0808179
1327	Regulation of mRNA stability by proteins that bind AU-rich elements	87	Up	0.0005483	0.0119021	-0.3498033	87	Up	0.0030772	0.0407199	-0.1258677
1333	Regulation of ornithine decarboxylase (ODC)	49	Up	0.0006154	0.0129229	-0.0625668	49	Up	0.0020758	0.0313315	-0.0797944
1343	Regulation of RUNX3 expression and activity	53	Up	0.0004381	0.0099566	-0.0603069	53	Up	0.0009308	0.0202055	-0.0634240
1370	Respiratory electron transport	90	Up	0.0006428	0.0130729	-0.1796005	90	Up	0.0003957	0.0119448	-0.0879771
1371	Respiratory electron transport, ATP synthesis by chemiosmotic coupling, and heat production by uncoupling proteins.	95	Up	0.0008812	0.0165291	-0.2761915	95	Up	0.0004838	0.0136077	-0.1032458
1373	Response of EIF2AK4 (GCN2) to amino acid deficiency	95	Up	0.0000000	0.0000043	0.1695889	95	Up	0.0000277	0.0053602	-0.0416488
1455	rRNA processing	186	Up	0.0000001	0.0000070	-0.0269377	186	Up	0.0000579	0.0069965	-0.1037470
1457	rRNA processing in the nucleus and cytosol	180	Up	0.0000000	0.0000048	0.0167065	180	Up	0.0000260	0.0053602	-0.0864008
1481	S Phase	157	Up	0.0032244	0.0429618	-0.3028758	157	Up	0.0007521	0.0177209	-0.1224841
1484	SARS-CoV-1 Infection	134	Up	0.0023168	0.0339097	-0.4089573	134	Up	0.0005379	0.0141466	-0.1150736
1485	SARS-CoV-1 modulates host translation machinery	34	Up	0.0000013	0.0000987	0.3873579	34	Up	0.0010493	0.0206856	-0.0033424
1487	SARS-CoV-1-host interactions	90	Up	0.0000798	0.0032797	-0.3008345	90	Up	0.0003943	0.0119448	-0.1003646
1491	SARS-CoV-2 modulates host translation machinery	46	Up	0.0000107	0.0007139	0.1064032	46	Up	0.0011256	0.0216971	-0.0448705
1498	SCF-beta-TrCP mediated degradation of Emi1	53	Up	0.0003208	0.0077737	-0.0603069	53	Up	0.0006602	0.0163518	-0.0634240
1499	SCF(Skp2)-mediated degradation of p27/p21	58	Up	0.0000335	0.0016160	0.0191124	58	Up	0.0002163	0.0092515	-0.0523222
1502	Selenoamino acid metabolism	103	Up	0.0000002	0.0000182	0.0679048	103	Up	0.0003186	0.0115052	-0.0802736
1503	Selenocysteine synthesis	87	Up	0.0000000	0.0000042	0.2056701	87	Up	0.0001209	0.0074514	-0.0707040
1510	Senescence-Associated Secretory Phenotype (SASP)	53	Up	0.0035504	0.0463470	-0.2311687	53	Up	0.0005953	0.0149371	-0.0524834
1512	Sensory Perception	570	Down	0.0000571	0.0025066	-1.2359325	567	Down	0.0000000	0.0000020	-0.5774976
1514	Sensory perception of sweet, bitter, and umami (glutamate) taste	41	Down	0.0000233	0.0012524	-1.4092232	41	Down	0.0001855	0.0088423	-0.5262750
1515	Sensory perception of taste	47	Down	0.0003067	0.0075968	-1.2863848	47	Down	0.0001922	0.0088423	-0.5262750
1519	Separation of Sister Chromatids	161	Up	0.0010889	0.0186286	-0.3812494	161	Up	0.0036781	0.0454093	-0.1409840
1603	Signaling by NOTCH4	79	Up	0.0008262	0.0159619	-0.3068824	79	Up	0.0012396	0.0220623	-0.1272730
1623	Signaling by ROBO receptors	203	Up	0.0000233	0.0012524	-0.2307918	203	Up	0.0012825	0.0225257	-0.1515047
1659	SRP-dependent cotranslational protein targeting to membrane	105	Up	0.0000000	0.0000057	0.0681744	105	Up	0.0000737	0.0071176	-0.0707040
1660	Stabilization of p53	55	Up	0.0001201	0.0042199	-0.0603069	55	Up	0.0005578	0.0141788	-0.0634240
1688	Switching of origins to a post-replicative state	90	Up	0.0022643	0.0336507	-0.1726080	90	Up	0.0001347	0.0074514	-0.0755749
1702	Synthesis of DNA	118	Up	0.0020686	0.0314691	-0.2814605	118	Up	0.0002522	0.0101496	-0.0819130
1758	The role of GTSE1 in G2/M progression after G2 checkpoint	58	Up	0.0001135	0.0040608	-0.1191575	58	Up	0.0002791	0.0110029	-0.0595988
1774	TNFR2 non-canonical NF-kB pathway	94	Up	0.0022917	0.0337978	-0.3859942	94	Up	0.0016043	0.0267202	-0.0898906
1831	Translation	262	Up	0.0000001	0.0000088	-0.1376113	262	Up	0.0000681	0.0071176	-0.1208969
1875	tRNA processing in the nucleus	55	Up	0.0033577	0.0444315	-0.1619091	55	Up	0.0011989	0.0219794	-0.0981367
1884	Ubiquitin Mediated Degradation of Phosphorylated Cdc25A	50	Up	0.0002215	0.0058628	-0.0614369	50	Up	0.0003912	0.0119448	-0.0595988
1885	Ubiquitin-dependent degradation of Cyclin D	50	Up	0.0001615	0.0045883	-0.0479101	50	Up	0.0015553	0.0261296	-0.0827369
1903	Vif-mediated degradation of APOBEC3G	50	Up	0.0004036	0.0092836	-0.0479101	50	Up	0.0012059	0.0219794	-0.0708212
1905	Viral mRNA Translation	84	Up	0.0000000	0.0000042	0.2406206	84	Up	0.0000850	0.0071436	-0.0552680
1918	Vpu mediated degradation of CD4	50	Up	0.0003367	0.0079805	-0.0479101	50	Up	0.0006973	0.0166322	-0.0595988

make scatterplots of the top sets

XMAX=max(res$ranked_profile[,1])
XMIN=min(res$ranked_profile[,1])
YMAX=max(res$ranked_profile[,2])
YMIN=min(res$ranked_profile[,2])
plot(res$detailed_sets[[1]], pch=19,cex=2,col="lightgray",
  xlim=c(XMIN,XMAX),ylim=c(YMIN,YMAX),
  xlab="blood",ylab="Guthrie")
text(res$detailed_sets[[1]],labels=rownames(res$detailed_sets[[1]]))
mtext(names(res$detailed_sets)[[1]])
grid() ; abline(v=0,h=0,lty=2,lwd=3)

pdf("c3_metabotropic_ilanguage.pdf")

plot(res$detailed_sets[[1]], pch=19,cex=2,col="lightgray",
  xlim=c(XMIN,XMAX),ylim=c(YMIN,YMAX), 
  xlab="blood",ylab="Guthrie")
text(res$detailed_sets[[1]],labels=rownames(res$detailed_sets[[1]]))
mtext(names(res$detailed_sets)[[1]])
grid() ; abline(v=0,h=0,lty=2,lwd=3)

dev.off()

## png 
##   2

Distinguishing ASD from birth

Here I will use MDS plot to cluster samples.

First I will use all probes, then I will focus on a subset of candidate probes.

plot(cmdscale(dist(t(bl_mvals))), xlab="Coordinate 1", ylab="Coordinate 2",
  type = "n",main="Blood")
mtext("Numbers indicate ilanguage score")
text(cmdscale(dist(t(bl_mvals))), labels=bl_design[,ncol(bl_design)], )

plot(cmdscale(dist(t(gu_mvals))), xlab="Coordinate 1", ylab="Coordinate 2", 
  type = "n",main="Guthrie card")
mtext("Numbers indicate ilanguage score")
text(cmdscale(dist(t(gu_mvals))), labels=gu_design[,ncol(gu_design)], )

There is no apparent clustering by ilanguage score, indicating a subtle effect.

Now I will examine the probes associated with the pathway of interest.

myprobes

##   [1] "cg01329260"    "cg06476845"    "cg05873271"    "cg12958612"   
##   [5] "cg14230917"    "cg16791170"    "cg16785092"    "cg15863419"   
##   [9] "cg17969276"    "cg21498511"    "cg08076643"    "cg06282219"   
##  [13] "cg25644447"    "cg05416591"    "cg00018443"    "cg04762133"   
##  [17] "cg18752890"    "cg06223856"    "cg04525998"    "cg14781700"   
##  [21] "cg13143099"    "cg18878381"    "cg11696355"    "cg12963342"   
##  [25] "cg21512216"    "cg08844988"    "cg16792948"    "cg00382632"   
##  [29] "cg04508233"    "cg08073142"    "cg01155522"    "cg09024808"   
##  [33] "cg27078729"    "cg04215186"    "cg09889256"    "cg06037603"   
##  [37] "cg14908168"    "cg02785870"    "cg23873145"    "cg07689148"   
##  [41] "cg06420129"    "cg12069304"    "cg15473007"    "cg16124834"   
##  [45] "cg12229081"    "cg01436550"    "cg04047827"    "cg02183851"   
##  [49] "cg25012186"    "cg01053405"    "cg00152838"    "cg05689121"   
##  [53] "cg02324079"    "cg03156546"    "cg01560476"    "cg02393980"   
##  [57] "cg27407707"    "cg22463563"    "cg08083403"    "cg05404890"   
##  [61] "cg09657673"    "cg09589586"    "cg27186436"    "cg07768777"   
##  [65] "cg10134076"    "cg07157902"    "cg03664025"    "cg24527964"   
##  [69] "cg09816107"    "cg11739924"    "cg14744741"    "cg18045878"   
##  [73] "cg02391633"    "cg07710734"    "cg14619832"    "cg18833170"   
##  [77] "cg17525025"    "cg05353134"    "cg08500967"    "ch.16.759907F"
##  [81] "cg01880908"    "cg04256065"    "cg01303655"    "cg12700283"   
##  [85] "cg20035319"    "cg00150128"    "cg18748888"    "cg12489727"   
##  [89] "cg15661240"    "cg12027254"    "cg11158641"    "cg03660034"   
##  [93] "cg05736847"    "cg17172331"    "cg05478790"    "cg26340657"   
##  [97] "cg19699973"    "cg09646181"    "cg13800273"    "cg08898446"   
## [101] "cg12072376"    "cg08446633"    "cg24207009"    "cg07787614"   
## [105] "cg02958960"    "cg20857275"    "cg23972793"    "cg27503577"   
## [109] "cg25362361"    "cg09029192"    "cg03760856"    "cg15744368"   
## [113] "cg07093496"    "cg09146982"    "cg24008619"    "cg13280512"   
## [117] "cg00423030"    "cg15489593"    "cg00175441"    "cg00496251"   
## [121] "cg00672192"    "cg00472027"    "cg26200971"    "cg27579805"   
## [125] "cg10943443"    "cg21571594"    "cg22094953"    "cg15831500"   
## [129] "cg05890693"    "cg04268643"    "cg04498014"    "cg07530460"   
## [133] "cg26740249"

myprobedat <- bl_mvals[which(rownames(bl_mvals) %in% myprobes),]
plot(cmdscale(dist(t(myprobedat))), xlab="Coordinate 1", ylab="Coordinate 2",
  type = "n",main="Blood")
text(cmdscale(dist(t(myprobedat))), labels=bl_design[,ncol(bl_design)], )

myprobedat <- gu_mvals[which(rownames(gu_mvals) %in% myprobes),]
plot(cmdscale(dist(t(myprobedat))), xlab="Coordinate 1", ylab="Coordinate 2",
  type = "n",main="Guthrie card")
text(cmdscale(dist(t(myprobedat))), labels=gu_design[,ncol(gu_design)], )

Again, there is no apparent clustering by ilanguage score. This indicates that the top probes cannot be used diagnostically.

myprobes2 <- head(bl_lim$Name,50)
myprobedat2 <- bl_mvals[which(rownames(bl_mvals) %in% myprobes2),]
plot(cmdscale(dist(t(myprobedat2))), xlab="Coordinate 1", ylab="Coordinate 2",
  type = "n",main="Blood")
text(cmdscale(dist(t(myprobedat2))), labels=bl_design[,ncol(bl_design)], )

myprobes2 <- head(gu_lim$Name,50)
myprobedat2 <- gu_mvals[which(rownames(gu_mvals) %in% myprobes2),]
plot(cmdscale(dist(t(myprobedat2))), xlab="Coordinate 1", ylab="Coordinate 2",
  type = "n", main="Guthrie")
text(cmdscale(dist(t(myprobedat2))), labels=gu_design[,ncol(gu_design)], )

MDS analysis of top probes identified without correction for covariates.

bl_design2 <- bl_design[,c(1,10)]
fit <- lmFit(bl_mvals, bl_design2)
fit <- eBayes(fit)
summary(decideTests(fit))

##        (Intercept) ilanguage
## Down        256222         0
## NotSig       16895    802647
## Up          529530         0

top <- topTable(fit,coef="ilanguage",num=Inf, sort.by = "P")
head(top)

##                 logFC    AveExpr         t      P.Value adj.P.Val            B
## cg17692763 -0.6457889  0.8069554 -5.915963 9.206791e-06 0.7379969  0.226393453
## cg04367545  0.4769302  0.8087474  5.737223 1.363368e-05 0.7379969  0.072659929
## cg03553358 -0.5139536 -2.3186818 -5.708938 1.451255e-05 0.7379969  0.047878940
## cg25692524 -0.3920523 -3.6852899 -5.648790 1.657936e-05 0.7379969 -0.005229242
## cg03146706 -0.4352132  1.7198109 -5.596459 1.862176e-05 0.7379969 -0.051892199
## cg12959820 -0.4433110 -0.3527632 -5.523459 2.190894e-05 0.7379969 -0.117697054

myprobes3 <- rownames(top)[1:20]

myprobedat3 <- bl_mvals[which(rownames(bl_mvals) %in% myprobes3),]
plot(cmdscale(dist(t(myprobedat3))), xlab="Coordinate 1", ylab="Coordinate 2",
  type = "n",main="Blood")
text(cmdscale(dist(t(myprobedat3))), labels=bl_design[,ncol(bl_design)], )

gu_design2 <- gu_design[,c(1,12)]
fit <- lmFit(gu_mvals, gu_design2)
fit <- eBayes(fit)
summary(decideTests(fit))

##        (Intercept) ilanguage
## Down        262383         0
## NotSig       22814    790658
## Up          505461         0

top <- topTable(fit,coef="ilanguage",num=Inf, sort.by = "P")
head(top)

##                 logFC    AveExpr         t      P.Value adj.P.Val          B
## cg02239677  0.4371694 -1.0143109  6.599736 8.543598e-07 0.6755064  1.0165260
## cg08971667 -0.4722478 -2.4549604 -6.308589 1.712085e-06 0.6768370  0.7667963
## cg20771596 -0.4323488  0.1339515 -5.773075 6.300529e-06 0.9999836  0.2747645
## cg04986899  0.3344568  2.0528811  5.090863 3.439218e-05 0.9999836 -0.4129538
## cg11100281 -0.4588872  1.1950259 -5.009866 4.215605e-05 0.9999836 -0.4989519
## cg09098195  0.5807053  1.9778800  4.991100 4.419382e-05 0.9999836 -0.5190026

myprobes3 <- rownames(top)[1:20]

myprobedat3 <- gu_mvals[which(rownames(gu_mvals) %in% myprobes3),]
plot(cmdscale(dist(t(myprobedat3))), xlab="Coordinate 1", ylab="Coordinate 2", 
  type = "n",main="Guthrie")
text(cmdscale(dist(t(myprobedat3))), labels=gu_design[,ncol(gu_design)], )

Again it was unsuccessful.

Try ebGSEA

ebGSEA uses globaltest to identify DMGs from infinium data. Some points to note:

only 1 sided test - does not distinguish between up and down methylated genes.
cannot use complex experiment designs.

Therefore I modify how the ebGSEA code uses globaltest.

item mapEIDto850k.lv : A list mapping Entrez Gene ID to 850k probes

#bl_gt <- doGT(bl_design,bl_mvals,array="850k",ncores=16)
data.m <- bl_mvals

subsets <- mclapply(mapEIDto850k.lv,intersect,rownames(data.m),mc.cores = 16)
nrep.v <- unlist(lapply(subsets,length))
selG.idx <- which(nrep.v>0)

gt.o <- gt(response=bl_design[,ncol(bl_design)], alternative=t(data.m), model="linear",
  directional = TRUE, standardize = FALSE, permutations = 0,
  subsets=subsets[selG.idx],trace=F)

gt_res <- result(gt.o)
gt_res <- gt_res[order(gt_res$`p-value`),]

Gene set level enrichment - ilanguage

genesets <- gmt_import("ReactomePathways.gmt")
cores = 8

Aggregate mval to genes and examine pathway enrichment in blood and guthries.

tic()
bl_gsmea <- gsmea(mval=bl_mvals,design=bl_design,probesets=sets,genesets=genesets,cores=8)
toc()

## 198.092 sec elapsed

bl_gsmea$res <- subset(bl_gsmea$res,nprobes>=5)
bl_gsmea$res$FDR <- p.adjust(bl_gsmea$res$PValue,method="fdr")
head(bl_gsmea$res,20)  %>%  kbl(caption="Top DM pathways in blood identified with GSMEA") %>%
  kable_paper("hover", full_width = F)

Top DM pathways in blood identified with GSMEA
	nprobes	median	PValue	FDR
VLDL assembly	5	-1.8955934	0.0000008	0.0016351
Vitamin D (calciferol) metabolism	12	-2.4124560	0.0000050	0.0047951
Methionine salvage pathway	6	-3.1134366	0.0000096	0.0052252
RHOH GTPase cycle	36	-1.1361085	0.0000129	0.0052252
Regulation of TLR by endogenous ligand	17	-1.6911724	0.0000135	0.0052252
Assembly of the ORC complex at the origin of replication	10	-3.9842186	0.0000189	0.0060902
Defective HLCS causes multiple carboxylase deficiency	7	-2.7445180	0.0000411	0.0063163
Dual Incision in GG-NER	39	-0.7927163	0.0000427	0.0063163
Inhibition of replication initiation of damaged DNA by RB1/E2F1	11	-1.7628462	0.0000586	0.0063163
Competing endogenous RNAs (ceRNAs) regulate PTEN translation	16	-3.0762735	0.0000736	0.0063163
Release of Hh-Np from the secreting cell	8	-1.8096843	0.0000774	0.0063163
HDL clearance	5	-2.0842172	0.0000846	0.0063163
Formation of Incision Complex in GG-NER	40	-1.4406897	0.0000853	0.0063163
Linoleic acid (LA) metabolism	7	-1.1580701	0.0000853	0.0063163
Synthesis of very long-chain fatty acyl-CoAs	23	-2.4423450	0.0000865	0.0063163
Role of LAT2/NTAL/LAB on calcium mobilization	13	-2.0112126	0.0000873	0.0063163
Defects in biotin (Btn) metabolism	8	-2.5911369	0.0000937	0.0063163
RHOD GTPase cycle	48	-1.9086189	0.0000949	0.0063163
Muscarinic acetylcholine receptors	5	-2.0532666	0.0001001	0.0063163
Fatty acyl-CoA biosynthesis	36	-2.3192477	0.0001027	0.0063163

tic()
gu_gsmea <- gsmea(mval=gu_mvals,design=gu_design,probesets=sets,genesets=genesets,cores=8)
toc()

## 194.893 sec elapsed

gu_gsmea$res <- subset(gu_gsmea$res,nprobes>=5)
gu_gsmea$res$FDR <- p.adjust(gu_gsmea$res$PValue,method="fdr")
head(gu_gsmea$res,20)  %>%  kbl(caption="Top DM pathways in Guthrie card identified with GSMEA") %>%
  kable_paper("hover", full_width = F)

Top DM pathways in Guthrie card identified with GSMEA
	nprobes	median	PValue	FDR
Recycling of eIF2:GDP	7	1.0530744	0.0025353	0.8838114
Sensory perception of salty taste	6	-1.9930074	0.0027402	0.8838114
Diseases associated with surfactant metabolism	9	-1.2577513	0.0058503	0.8838114
KSRP (KHSRP) binds and destabilizes mRNA	17	-0.6957056	0.0131495	0.8838114
Chylomicron remodeling	10	-0.8226941	0.0139064	0.8838114
Creatine metabolism	9	-0.8008186	0.0140569	0.8838114
Defective CSF2RA causes SMDP4	7	-1.2577513	0.0210804	0.8838114
Defective CSF2RB causes SMDP5	7	-1.2577513	0.0210804	0.8838114
Tandem pore domain potassium channels	12	-1.1411864	0.0218081	0.8838114
Inactivation of APC/C via direct inhibition of the APC/C complex	21	0.4044948	0.0235526	0.8838114
Inhibition of the proteolytic activity of APC/C required for the onset of anaphase by mitotic spindle checkpoint components	21	0.4044948	0.0235526	0.8838114
Cation-coupled Chloride cotransporters	7	-1.1202600	0.0341080	0.8838114
Na+/Cl- dependent neurotransmitter transporters	18	-0.9272093	0.0352764	0.8838114
Activated NTRK3 signals through PI3K	6	-1.0512534	0.0392557	0.8838114
APC-Cdc20 mediated degradation of Nek2A	26	0.3854902	0.0399509	0.8838114
Assembly of active LPL and LIPC lipase complexes	17	-0.9537194	0.0418350	0.8838114
Phase 4 - resting membrane potential	19	-0.9553124	0.0431821	0.8838114
Signaling by NTRK3 (TRKC)	17	-0.6776604	0.0546364	0.8838114
Chylomicron assembly	10	-0.5346036	0.0579233	0.8838114
AKT phosphorylates targets in the nucleus	9	-1.1352247	0.0604976	0.8838114

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.2 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] vioplot_0.4.0                                      
##  [2] zoo_1.8-11                                         
##  [3] sm_2.2-5.7.1                                       
##  [4] pkgload_1.3.2                                      
##  [5] GGally_2.1.2                                       
##  [6] ggplot2_3.4.1                                      
##  [7] gtools_3.9.4                                       
##  [8] tibble_3.2.0                                       
##  [9] dplyr_1.1.0                                        
## [10] echarts4r_0.4.4                                    
## [11] ebGSEA_0.1.0                                       
## [12] globaltest_5.52.0                                  
## [13] survival_3.5-5                                     
## [14] tictoc_1.1                                         
## [15] RIdeogram_0.2.2                                    
## [16] kableExtra_1.3.4                                   
## [17] data.table_1.14.8                                  
## [18] ENmix_1.34.0                                       
## [19] doParallel_1.0.17                                  
## [20] qqman_0.1.8                                        
## [21] RCircos_1.2.2                                      
## [22] beeswarm_0.4.0                                     
## [23] forestplot_3.1.1                                   
## [24] abind_1.4-5                                        
## [25] checkmate_2.1.0                                    
## [26] reshape2_1.4.4                                     
## [27] gplots_3.1.3                                       
## [28] eulerr_7.0.0                                       
## [29] GEOquery_2.66.0                                    
## [30] RColorBrewer_1.1-3                                 
## [31] IlluminaHumanMethylation450kmanifest_0.4.0         
## [32] topconfects_1.14.0                                 
## [33] DMRcatedata_2.16.0                                 
## [34] ExperimentHub_2.6.0                                
## [35] AnnotationHub_3.6.0                                
## [36] BiocFileCache_2.6.0                                
## [37] dbplyr_2.3.1                                       
## [38] DMRcate_2.12.0                                     
## [39] limma_3.54.0                                       
## [40] missMethyl_1.32.0                                  
## [41] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1 
## [42] R.utils_2.12.2                                     
## [43] R.oo_1.25.0                                        
## [44] R.methodsS3_1.8.2                                  
## [45] plyr_1.8.8                                         
## [46] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [47] minfi_1.44.0                                       
## [48] bumphunter_1.40.0                                  
## [49] locfit_1.5-9.7                                     
## [50] iterators_1.0.14                                   
## [51] foreach_1.5.2                                      
## [52] Biostrings_2.66.0                                  
## [53] XVector_0.38.0                                     
## [54] SummarizedExperiment_1.28.0                        
## [55] Biobase_2.58.0                                     
## [56] MatrixGenerics_1.10.0                              
## [57] matrixStats_0.63.0                                 
## [58] GenomicRanges_1.50.2                               
## [59] GenomeInfoDb_1.34.6                                
## [60] IRanges_2.32.0                                     
## [61] S4Vectors_0.36.1                                   
## [62] BiocGenerics_0.44.0                                
## [63] mitch_1.10.0                                       
## 
## loaded via a namespace (and not attached):
##   [1] rappdirs_0.3.3                rtracklayer_1.58.0           
##   [3] tidyr_1.3.0                   bit64_4.0.5                  
##   [5] knitr_1.42                    DelayedArray_0.24.0          
##   [7] rpart_4.1.19                  KEGGREST_1.38.0              
##   [9] RCurl_1.98-1.10               AnnotationFilter_1.22.0      
##  [11] generics_0.1.3                GenomicFeatures_1.50.3       
##  [13] preprocessCore_1.60.2         RSQLite_2.3.0                
##  [15] bit_4.0.5                     tzdb_0.3.0                   
##  [17] webshot_0.5.4                 xml2_1.3.3                   
##  [19] httpuv_1.6.9                  assertthat_0.2.1             
##  [21] xfun_0.37                     hms_1.1.2                    
##  [23] jquerylib_0.1.4               evaluate_0.20                
##  [25] promises_1.2.0.1              fansi_1.0.4                  
##  [27] restfulr_0.0.15               scrime_1.3.5                 
##  [29] progress_1.2.2                caTools_1.18.2               
##  [31] readxl_1.4.2                  DBI_1.1.3                    
##  [33] geneplotter_1.76.0            htmlwidgets_1.6.2            
##  [35] reshape_0.8.9                 purrr_1.0.1                  
##  [37] ellipsis_0.3.2                backports_1.4.1              
##  [39] permute_0.9-7                 calibrate_1.7.7              
##  [41] grImport2_0.2-0               annotate_1.76.0              
##  [43] biomaRt_2.54.0                deldir_1.0-6                 
##  [45] sparseMatrixStats_1.10.0      vctrs_0.6.0                  
##  [47] ensembldb_2.22.0              withr_2.5.0                  
##  [49] cachem_1.0.7                  Gviz_1.42.0                  
##  [51] BSgenome_1.66.2               GenomicAlignments_1.34.0     
##  [53] prettyunits_1.1.1             mclust_6.0.0                 
##  [55] svglite_2.1.1                 cluster_2.1.4                
##  [57] RPMM_1.25                     lazyeval_0.2.2               
##  [59] crayon_1.5.2                  genefilter_1.80.3            
##  [61] labeling_0.4.2                edgeR_3.40.2                 
##  [63] pkgconfig_2.0.3               nlme_3.1-162                 
##  [65] ProtGenerics_1.30.0           nnet_7.3-18                  
##  [67] rlang_1.1.0                   lifecycle_1.0.3              
##  [69] filelock_1.0.2                dichromat_2.0-0.1            
##  [71] rsvg_2.4.0                    cellranger_1.1.0             
##  [73] rngtools_1.5.2                base64_2.0.1                 
##  [75] Matrix_1.5-3                  Rhdf5lib_1.20.0              
##  [77] base64enc_0.1-3               viridisLite_0.4.1            
##  [79] png_0.1-8                     rjson_0.2.21                 
##  [81] bitops_1.0-7                  KernSmooth_2.23-20           
##  [83] rhdf5filters_1.10.0           blob_1.2.4                   
##  [85] DelayedMatrixStats_1.20.0     doRNG_1.8.6                  
##  [87] stringr_1.5.0                 nor1mix_1.3-0                
##  [89] readr_2.1.4                   jpeg_0.1-10                  
##  [91] scales_1.2.1                  memoise_2.0.1                
##  [93] magrittr_2.0.3                zlibbioc_1.44.0              
##  [95] compiler_4.2.2                BiocIO_1.8.0                 
##  [97] illuminaio_0.40.0             Rsamtools_2.14.0             
##  [99] cli_3.6.0                     DSS_2.46.0                   
## [101] htmlTable_2.4.1               Formula_1.2-5                
## [103] MASS_7.3-58.3                 tidyselect_1.2.0             
## [105] stringi_1.7.12                highr_0.10                   
## [107] yaml_2.3.7                    askpass_1.1                  
## [109] latticeExtra_0.6-30           sass_0.4.5                   
## [111] VariantAnnotation_1.44.0      tools_4.2.2                  
## [113] rstudioapi_0.14               foreign_0.8-84               
## [115] bsseq_1.34.0                  gridExtra_2.3                
## [117] farver_2.1.1                  digest_0.6.31                
## [119] BiocManager_1.30.20           shiny_1.7.4                  
## [121] quadprog_1.5-8                Rcpp_1.0.10                  
## [123] siggenes_1.72.0               BiocVersion_3.16.0           
## [125] later_1.3.0                   org.Hs.eg.db_3.16.0          
## [127] httr_1.4.5                    AnnotationDbi_1.60.0         
## [129] biovizBase_1.46.0             colorspace_2.1-0             
## [131] rvest_1.0.3                   XML_3.99-0.13                
## [133] splines_4.2.2                 statmod_1.5.0                
## [135] kpmt_0.1.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                      Hmisc_5.0-1                  
## [143] pillar_1.8.1                  htmltools_0.5.4              
## [145] mime_0.12                     glue_1.6.2                   
## [147] fastmap_1.1.1                 BiocParallel_1.32.5          
## [149] interactiveDisplayBase_1.36.0 beanplot_1.3.1               
## [151] codetools_0.2-19              utf8_1.2.3                   
## [153] lattice_0.20-45               bslib_0.4.2                  
## [155] curl_5.0.0                    openssl_2.0.6                
## [157] interp_1.1-3                  rmarkdown_2.20               
## [159] munsell_0.5.0                 rhdf5_2.42.0                 
## [161] GenomeInfoDbData_1.2.9        HDF5Array_1.26.0             
## [163] impute_1.72.3                 gtable_0.3.2

Correlates of inverse IQ and methylation - enrichment analysis

Mark Ziemann

2023-03-20