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

Introduction

Here we are analysing Guthrie card methylation in 9 twin pairs.

baseDir=getwd()
dataDir=paste(baseDir,"/ASD_EPIC_DATA",sep="")

suppressPackageStartupMessages({
  library("missMethyl")
  library("limma")
  library("minfi")
  library("IlluminaHumanMethylation450kanno.ilmn12.hg19")
  library("IlluminaHumanMethylationEPICanno.ilm10b2.hg19")
  library("ruv")
  library("RColorBrewer")
  library("matrixStats")
  library("gplots")
  library("WGCNA")
  library("FlowSorted.Blood.450k")
  library("reshape2")
  library("ggplot2")
  library("missMethyl")
  library("DMRcate")
  library("FlowSorted.Blood.EPIC")
  library("mitch")
  library("kableExtra")
})

source("meth_functions.R")

Load data

Load the annotation data and the Epic methylation data.

This analysis is to be conducted on Ubuntu with R4.

ann = getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b2.hg19)

ann_sub = ann[,c("chr","pos","strand","Name","Islands_Name",
    "Relation_to_Island","UCSC_RefGene_Name","UCSC_RefGene_Group")]

targets_gen = read.metharray.sheet(dataDir, pattern = "ASD_guthrie_summarysheet.csv")
## [read.metharray.sheet] Found the following CSV files:
## [1] "/home/mdz/projects/asd_meth/ASD_EPIC_DATA/ASD_guthrie_summarysheet.csv"
#targets$ID = paste(targets$Sample_Group,targets_gen$Sample_Name,sep=".")
rgSet = read.metharray.exp(targets = targets_gen)
sampleNames(rgSet) = targets_gen$SampleID

Testing MZ status

snpBetas = getSnpBeta(rgSet)
d = dist(t(snpBetas))
hr = hclust(d, method = "complete", members=NULL)
plot(hr)

Quality control

detP = detectionP(rgSet)
qcReport(rgSet, sampNames = targets_gen$Sample_Name,
  pdf="qc-report_ASD_guthrie_nov27.pdf")
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE

## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## png 
##   2
cols=brewer.pal(4,"Set1")

barplot(apply(detP,2,mean),
  col=as.numeric(factor(targets_gen$Sample_Name)),
  las=2,cex.names= 0.8, cex.axis=0.75,
  main="Mean detection p-values of probe signals",
  ylab="Mean detection p-value")

barplot(apply(detP,2,mean),
  col=as.numeric(factor(targets_gen$Sample_Name)),
  las=2,cex.names= 0.8, cex.axis=0.75,ylim=c(0,0.010),
  main="Mean detection p-values of probe signals",
  ylab="Mean detection p-value")

Cell type composition analysis

Using old method.

rgSet$Slide <- as.numeric(rgSet$Slide)
rgSet$Sex <- as.character(rgSet$Sex)
#rgSet$Sample_Name <- as.character(rgSet$Sample_Name)
#sampleNames(rgSet) = targets_gen$SampleID

cellCounts_new <- estimateCellCounts(rgSet, compositeCellType = "Blood", 
  processMethod = "auto", probeSelect = "auto", 
  cellTypes = c("CD8T","CD4T", "NK","Bcell","Mono","Gran"),
  referencePlatform = c("IlluminaHumanMethylation450k"),
  returnAll = FALSE, meanPlot = TRUE)
## [estimateCellCounts] Combining user data with reference (flow sorted) data.
## Warning in DataFrame(sampleNames = c(colnames(rgSet),
## colnames(referenceRGset)), : 'stringsAsFactors' is ignored
## [estimateCellCounts] Processing user and reference data together.
## [preprocessQuantile] Mapping to genome.
## [preprocessQuantile] Fixing outliers.
## [preprocessQuantile] Quantile normalizing.
## [estimateCellCounts] Picking probes for composition estimation.
## [estimateCellCounts] Estimating composition.

#plot cell type composition by sample group
a = cellCounts_new[targets_gen$Diagnosis == "0",]
b = cellCounts_new[targets_gen$Diagnosis == "1",]
c = cellCounts_new[targets_gen$Diagnosis == "2",]
age.pal <- brewer.pal(8,"Set1")

cellCounts_long <- melt(cellCounts_new, id = "celltype")

cellCounts_long$Var1 <- targets_gen[match(cellCounts_long$Var1,  targets_gen$SampleID),"Diagnosis"]

colnames(cellCounts_long) <- c("diagnosis","celltype","value")

ggplot(cellCounts_long, aes(x = factor(celltype), y = value, fill = factor(diagnosis))) + 
  geom_boxplot() 

wx <- lapply(1:ncol(cellCounts_new) , function(i) {
  wilcox.test(a[,i],c[,i], paired=FALSE)
} )
## Warning in wilcox.test.default(a[, i], c[, i], paired = FALSE): cannot compute
## exact p-value with ties
names(wx) <- colnames(cellCounts_new)
wx
## $CD8T
## 
##  Wilcoxon rank sum exact test
## 
## data:  a[, i] and c[, i]
## W = 17, p-value = 0.2674
## alternative hypothesis: true location shift is not equal to 0
## 
## 
## $CD4T
## 
##  Wilcoxon rank sum exact test
## 
## data:  a[, i] and c[, i]
## W = 37, p-value = 0.3196
## alternative hypothesis: true location shift is not equal to 0
## 
## 
## $NK
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  a[, i] and c[, i]
## W = 22, p-value = 0.491
## alternative hypothesis: true location shift is not equal to 0
## 
## 
## $Bcell
## 
##  Wilcoxon rank sum exact test
## 
## data:  a[, i] and c[, i]
## W = 36, p-value = 0.3773
## alternative hypothesis: true location shift is not equal to 0
## 
## 
## $Mono
## 
##  Wilcoxon rank sum exact test
## 
## data:  a[, i] and c[, i]
## W = 23, p-value = 0.6612
## alternative hypothesis: true location shift is not equal to 0
## 
## 
## $Gran
## 
##  Wilcoxon rank sum exact test
## 
## data:  a[, i] and c[, i]
## W = 26, p-value = 0.913
## alternative hypothesis: true location shift is not equal to 0

Try estimatecellcounts2.

cells <- estimateCellCounts2(rgSet, referencePlatform= "IlluminaHumanMethylationEPIC",
  returnAll = TRUE)
## snapshotDate(): 2022-10-31
## snapshotDate(): 2022-10-31
## see ?FlowSorted.Blood.EPIC and browseVignettes('FlowSorted.Blood.EPIC') for documentation
## loading from cache
## [estimateCellCounts2] Combining user data with reference (flow sorted) data.
## Warning in asMethod(object): NAs introduced by coercion
## [estimateCellCounts2] Processing user and reference data together.
## [estimateCellCounts2] Using IDOL L-DMR probes for composition estimation.
## [estimateCellCounts2] Estimating proportion composition (prop), if you provide cellcounts those will be provided as counts in the composition estimation.
mset <- cells$normalizedData

cellCounts_new <- cells[[1]]
#plot cell type composition by sample group
a = cellCounts_new[targets_gen$Diagnosis == "0",]
b = cellCounts_new[targets_gen$Diagnosis == "1",]
c = cellCounts_new[targets_gen$Diagnosis == "2",]
age.pal <- brewer.pal(8,"Set1")

cellCounts_long <- melt(cellCounts_new, id = "celltype")

cellCounts_long$Var1 <- targets_gen[match(cellCounts_long$Var1,  targets_gen$SampleID),"Diagnosis"]

colnames(cellCounts_long) <- c("diagnosis","celltype","value")

ggplot(cellCounts_long, aes(x = factor(celltype), y = value, fill = factor(diagnosis))) +
  geom_boxplot()

pdf("cells_guthrie.pdf")

par(mar=c(3,8,2,2))

ggplot(cellCounts_long, aes(x = factor(celltype), y = value, fill = factor(diagnosis))) +
  geom_boxplot()

dev.off()
## png 
##   2
wx <- lapply(1:ncol(cellCounts_new) , function(i) {
  wilcox.test(a[,i],c[,i], paired=FALSE)
} )
names(wx) <- colnames(cellCounts_new)
wx
## $CD8T
## 
##  Wilcoxon rank sum exact test
## 
## data:  a[, i] and c[, i]
## W = 32, p-value = 0.6612
## alternative hypothesis: true location shift is not equal to 0
## 
## 
## $CD4T
## 
##  Wilcoxon rank sum exact test
## 
## data:  a[, i] and c[, i]
## W = 36, p-value = 0.3773
## alternative hypothesis: true location shift is not equal to 0
## 
## 
## $NK
## 
##  Wilcoxon rank sum exact test
## 
## data:  a[, i] and c[, i]
## W = 19, p-value = 0.3773
## alternative hypothesis: true location shift is not equal to 0
## 
## 
## $Bcell
## 
##  Wilcoxon rank sum exact test
## 
## data:  a[, i] and c[, i]
## W = 39, p-value = 0.2212
## alternative hypothesis: true location shift is not equal to 0
## 
## 
## $Mono
## 
##  Wilcoxon rank sum exact test
## 
## data:  a[, i] and c[, i]
## W = 20, p-value = 0.4409
## alternative hypothesis: true location shift is not equal to 0
## 
## 
## $Neu
## 
##  Wilcoxon rank sum exact test
## 
## data:  a[, i] and c[, i]
## W = 29, p-value = 0.913
## alternative hypothesis: true location shift is not equal to 0

Filtering

None of the p-values are less than 0.05.

Now filter for high detection p-value and overlap with SNP.

Need to get a copy of the Xreact probes from Namitha. In the mean time I have downloaded from github - link below.

detP <- detectionP(rgSet)
keep <- rowSums(detP < 0.01) == ncol(rgSet)
mset <- mset[keep,]
gmset <- mapToGenome(mset)

#remove SNPs
gmset_flt = dropLociWithSnps(gmset, snps = c("CpG", "SBE"))

#Removing cross-reactive probes
XURL="https://raw.githubusercontent.com/sirselim/illumina450k_filtering/master/EPIC/13059_2016_1066_MOESM1_ESM.csv"
Xreact <- read.csv(XURL)

#Xreact = read.csv(file="/group/canc2/puumba/Data/InfiniumData/NamithaData/Rprojects/Autism/Analysis_Sept11/EPIC_850k_crossreactiveProbes.csv", stringsAsFactors=FALSE)
#Xreact = read.csv(file="~/48639-non-specific-probes-Illumina450k.csv", stringsAsFactors=FALSE)
noXreact <-  !(featureNames(gmset) %in% Xreact$X)

gmset <- gmset[noXreact,]

#Removing probes on X and Y chromosomes
autosomes <- !(featureNames(gmset) %in% ann$Name[ann$chr %in% c("chrX","chrY")])
gmset_flt <- gmset[autosomes,]

#Relative log expression (RLE plot)
mvals = getM(gmset_flt)
medSq = apply(mvals, 1, median)
YSq = mvals - medSq

boxplot(YSq,outline=FALSE,ylim=c(-1.5,1.5), 
  ylab="Relative Log Methylation Value", 
  cols=as.character(factor(targets_gen$Social_interaction_on_ADOS,)) )

MDS plots generation after filtering

pal = brewer.pal(8, "Dark2")
mds1Sq = plotMDS(mvals, top=1000, gene.selection="common",dim.plot=c(1,2))

mds2Sq = plotMDS(mvals, top=1000, gene.selection="common",dim.plot=c(1,3))

mds3Sq = plotMDS(mvals, top=1000, gene.selection="common",dim.plot=c(2,3))

mds4Sq = plotMDS(mvals, top=1000, gene.selection="common",dim.plot=c(3,4))

plotMDS(mds1Sq, xlab="Dimension 1", ylab="Dimension 2",
  col=pal[as.factor(targets_gen$Diagnosis)],
  dim=c(1,2), labels=targets_gen$Sample_Name)
legend("bottomright",bg="white",col=pal,cex=.7,pch=1,legend=0:1)

plotMDS(mds2Sq, xlab="Dimension 1", ylab="Dimension 3",
  col=pal[as.factor(targets_gen$Diagnosis)],dim=c(1,3),
  labels=targets_gen$Sample_Name)
legend("bottomright",bg="white",col=pal,cex=.7,pch=1,legend=0:1)

plotMDS(mds3Sq, xlab="Dimension 2", ylab="Dimension 3",
  col=pal[as.factor(targets_gen$Diagnosis)],dim=c(2,3),
  labels=targets_gen$Sample_Name)
legend("bottomright",bg="white",col=pal,cex=.7,pch=1,legend=0:1)

plotMDS(mds4Sq, xlab="Dimension 3", ylab="Dimension 4",
  col=pal[as.factor(targets_gen$Diagnosis)],dim=c(3,4),
  labels=targets_gen$Sample_Name)
legend("bottomright",bg="white",col=pal,cex=.7,pch=1,legend=0:1)

Principal Component Analysis (PCA)

fit <- prcomp(t(mvals),center = TRUE, scale = TRUE,retx=TRUE)
loadings = fit$x
plot(fit,type="lines",col="blue")

saveRDS(fit,"scree_guthrie_fircorrect.Rds")

fitnc <- readRDS("scree_guthrie_fitnc.Rds")
plot(fitnc,type="lines",col="blue")

nGenes = nrow(mvals)
nSamples = ncol(mvals)
datTraits = targets_gen[,7:15]
moduleTraitCor = cor(loadings[,1:6], datTraits, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
par(cex=0.75, mar = c(6, 8.5, 3, 3))
textMatrix = paste(signif(moduleTraitCor, 2), "\n(", 
  signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)

labeledHeatmap(Matrix = t(moduleTraitCor), 
  xLabels = colnames(loadings)[1:ncol(t(moduleTraitCor))],
  yLabels = names(datTraits), colorLabels = FALSE, colors = blueWhiteRed(6),
  textMatrix = t(textMatrix), setStdMargins = FALSE, cex.text = 0.5,
  cex.lab.y = 0.6, zlim = c(-1,1), 
  main = paste("PCA-trait relationships: Top principal components"))

pdf("pca_guthrie.pdf")

par(mar=c(3,8,2,2))

labeledHeatmap(Matrix = t(moduleTraitCor), 
  xLabels = colnames(loadings)[1:ncol(t(moduleTraitCor))],
  yLabels = names(datTraits), colorLabels = FALSE, colors = blueWhiteRed(6),
  textMatrix = t(textMatrix), setStdMargins = FALSE, cex.text = 0.5,
  cex.lab.y = 0.6, zlim = c(-1,1), 
  main = paste("PCA-trait relationships: Top principal components"))

dev.off()
## png 
##   2

Regression analysis - within twin pair on ADOS

ADOS is the outcome variable (continuous)

Twin_Group <- factor(targets_gen$Family_ID)
ADOS <- targets_gen$Social_interaction_on_ADOS
names(ADOS) <- targets_gen$SampleID
motor <- targets_gen$Motor_skills
diagnosis <- targets_gen$Diagnosis
SRS<- targets_gen$SRS_social_scores
iiq <- 1/targets_gen$IQ
# fix infinite
iiq[iiq==Inf] <- 0.025
iiq <- as.vector(scale(iiq))
ilanguage <- as.vector(scale(1/targets_gen$language))

#ADOS:0,4.009335e-06
design_tw <- model.matrix(~ Twin_Group + ADOS )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))
##        (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down        245993        2344        2566        2439        1909        3492
## NotSig       59114      786047      786358      786292      786913      785103
## Up          485551        2267        1734        1927        1836        2063
##        Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13   ADOS
## Down         48678        2099        2939         1854         2557      0
## NotSig      728154      786567      785798       786943       786409 790658
## Up           13826        1992        1921         1861         1692      0
top <- topTable(fit2_tw,coef="ADOS",num=Inf, sort.by = "P")
nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)
## [1] 21151
output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_guthrie_ADOS.csv",row.names=FALSE)
saveRDS(design_tw, "gu_design_ados.rds")
saveRDS(mvals, "gu_mvals.rds")
head(output,30) %>% kbl(caption="ADOS") %>% kable_paper("hover", full_width = F)
ADOS
Name chr pos strand Islands_Name Relation_to_Island UCSC_RefGene_Name UCSC_RefGene_Group logFC AveExpr t P.Value adj.P.Val B
cg12565057 chr7 100482617
chr7:100484425-100484651 N_Shore SRRT;SRRT;SRRT;SRRT Body;Body;Body;Body 0.0826621 3.0460140 7.387285 0.0000040 0.9999951 4.581264
cg00546248 chr1 2412743
chr1:2411125-2411404 S_Shore PLCH2 Body 0.1034780 1.8014843 6.865500 0.0000089 0.9999951 3.850448
cg15740041 chr16 87781336
OpenSea KLHDC4;KLHDC4;KLHDC4 Body;Body;Body 0.1317611 1.6008461 6.823044 0.0000095 0.9999951 3.789115
cg05116656 chr5 94417892
OpenSea MCTP1;MCTP1;MCTP1 TSS1500;TSS1500;Body 0.0915714 -3.3750586 6.719336 0.0000112 0.9999951 3.638100
cg11632273 chr22 20102580
chr22:20104374-20105729 N_Shore TRMT2A;TRMT2A;TRMT2A;RANBP1;MIR6816 Body;Body;Body;TSS1500;TSS1500 -0.1332566 2.2681145 -6.338476 0.0000208 0.9999951 3.068838
cg01656620 chr12 124788698
chr12:124788559-124788884 Island FAM101A 5’UTR -0.1330939 3.0511986 -6.203748 0.0000259 0.9999951 2.861912
cg10881749 chr1 234792390
OpenSea 0.0804887 0.9589997 6.114550 0.0000301 0.9999951 2.723318
cg13699963 chr17 71548861
OpenSea SDK2 Body -0.0663828 2.9801617 -5.750106 0.0000558 0.9999951 2.143861
cg22656126 chr17 1637206
chr17:1637053-1637492 Island WDR81;WDR81;WDR81;WDR81 Body;Body;Body;Body 0.1395930 5.7130309 5.698989 0.0000609 0.9999951 2.060902
cg07162704 chr6 91113687
OpenSea 0.1009010 0.5704454 5.695702 0.0000613 0.9999951 2.055553
cg13277464 chr2 43393040
OpenSea 0.1211560 2.3025832 5.606707 0.0000715 0.9999951 1.910100
cg18731055 chr13 49108131
chr13:49106411-49107463 S_Shore RCBTB2 TSS1500 0.1016068 0.0882235 5.584670 0.0000743 0.9999951 1.873891
cg03188919 chr10 31372663
OpenSea -0.1030369 2.0662356 -5.526631 0.0000823 0.9999951 1.778167
cg01128603 chr11 65820391
chr11:65819367-65820224 S_Shore SF3B2 Body 0.0946490 -3.2440188 5.486222 0.0000883 0.9999951 1.711213
cg03177876 chr4 71701931
chr4:71704664-71705736 N_Shelf GRSF1;GRSF1 5’UTR;Body 0.0982965 -1.2874739 5.481370 0.0000891 0.9999951 1.703158
cg05791256 chr9 96077420
chr9:96080113-96080330 N_Shelf WNK2;WNK2 Body;Body -0.0749658 3.6607111 -5.477688 0.0000896 0.9999951 1.697042
cg25723331 chr11 44613980
OpenSea CD82;CD82 5’UTR;5’UTR 0.1123597 0.8048882 5.413638 0.0001004 0.9999951 1.590320
cg01016829 chr16 46682692
OpenSea 0.0860670 2.2697378 5.369826 0.0001084 0.9999951 1.516959
cg23157290 chr7 138602695
OpenSea KIAA1549;KIAA1549 Body;Body 0.1583347 3.2345901 5.356147 0.0001111 0.9999951 1.493994
cg25652742 chr2 240152334
chr2:240153476-240153792 N_Shore HDAC4 Body 0.1179861 -1.9251649 5.354083 0.0001115 0.9999951 1.490528
cg23061257 chr1 54355272
chr1:54355354-54355656 N_Shore YIPF1 5’UTR -0.0648082 -3.3673834 -5.351267 0.0001121 0.9999951 1.485795
cg08008144 chr8 1983127
OpenSea -0.0927890 2.5786632 -5.319054 0.0001187 0.9999951 1.431580
cg23404427 chr19 1653402
OpenSea TCF3 TSS1500 -0.0903479 -5.0294197 -5.300923 0.0001226 0.9999951 1.400998
cg10619752 chr2 233863457
OpenSea NGEF 5’UTR -0.0820301 3.0526878 -5.290482 0.0001249 0.9999951 1.383363
cg08225137 chr3 54915807
OpenSea CACNA2D3 Body -0.0561511 2.6468160 -5.227687 0.0001397 0.9999951 1.276965
cg15856662 chr22 38380501
chr22:38379093-38379964 S_Shore SOX10;SOX10;POLR2F;POLR2F 5’UTR;1stExon;Body;Body -0.0734815 2.4748310 -5.180008 0.0001522 0.9999951 1.195788
cg11355135 chr2 108602860
chr2:108602824-108603467 Island SLC5A7 TSS200 0.1104571 -3.8226413 5.113981 0.0001714 0.9999951 1.082825
cg07873041 chr1 91301840
chr1:91300979-91301891 Island 0.0679299 -5.4645582 5.090814 0.0001788 0.9999951 1.043041
cg15165643 chr8 145518600
chr8:145514719-145515959 S_Shelf HSF1 Body -0.0814095 3.0943937 -5.080871 0.0001820 0.9999951 1.025942
cg00183468 chr3 142926178
OpenSea 0.0609273 1.9139298 5.080236 0.0001822 0.9999951 1.024849

Other models:

  • motor skills impairment

  • diagnosis

  • SRS

  • inverse IQ

  • inverse langguage

Regression analysis - within twin pair on motor skills impairment

#motor:820,1.053055e-12
design_tw <- model.matrix(~ Twin_Group + motor )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))
##        (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down        222835        3100        3073        1578        1451        2457
## NotSig      109032      784572      785915      788025      788150      786947
## Up          458791        2986        1670        1055        1057        1254
##        Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13  motor
## Down        152974        2956        3138         1373         1855    433
## NotSig      607875      784741      785529       788188       787768 789838
## Up           29809        2961        1991         1097         1035    387
top <- topTable(fit2_tw,coef="motor",num=Inf, sort.by = "P")
nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)
## [1] 53567
output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_guthrie_motor.csv",row.names=FALSE)
head(output,30) %>% kbl(caption="motor") %>% kable_paper("hover", full_width = F)
motor
Name chr pos strand Islands_Name Relation_to_Island UCSC_RefGene_Name UCSC_RefGene_Group logFC AveExpr t P.Value adj.P.Val B
cg22735331 chr11 2369390
OpenSea -4.487554 2.7456613 -24.58356 0 0.0000008 11.331647
cg01593834 chr5 145718401
chr5:145718289-145720095 Island POU4F3 TSS200 -3.047554 -1.4836163 -19.69127 0 0.0000080 10.529614
cg08840551 chr5 56718349
OpenSea -2.630660 -0.3473951 -18.30316 0 0.0000140 10.208855
cg09002832 chr9 4300181
chr9:4297817-4300182 Island GLIS3 TSS200 2.576547 -4.7977912 17.43241 0 0.0000198 9.978451
cg01059824 chr15 69706539
chr15:69706516-69707352 Island KIF23;KIF23;KIF23 TSS200;TSS200;TSS200 2.125246 -4.7737263 16.53201 0 0.0000281 9.712490
cg05098590 chr10 111767379
chr10:111767087-111768355 Island ADD3;ADD3;ADD3 5’UTR;TSS1500;TSS1500 2.378836 -5.5155179 16.45865 0 0.0000281 9.689455
cg23199907 chr13 33305966
OpenSea PDS5B Body -3.696135 2.7389798 -16.07868 0 0.0000312 9.566639
cg25259284 chr1 52521812
chr1:52521780-52522296 Island TXNDC12;TXNDC12;BTF3L4;BTF3L4;BTF3L4;TXNDC12;TXNDC12 1stExon;5’UTR;TSS200;TSS200;TSS200;Body;Body 2.467314 -4.1955426 15.94548 0 0.0000312 9.522146
cg09804439 chr17 40540457
chr17:40539837-40540775 Island STAT3;STAT3;STAT3;STAT3;STAT3 5’UTR;1stExon;1stExon;5’UTR;TSS200 2.103803 -5.6192796 15.82500 0 0.0000312 9.481235
cg04055615 chr6 31697723
chr6:31695894-31698245 Island DDAH2 5’UTR 2.215811 -4.9498235 15.15245 0 0.0000449 9.240656
cg07241733 chr8 68255213
chr8:68255129-68256274 Island ARFGEF1 Body 2.256471 -3.1368420 15.14960 0 0.0000449 9.239591
cg19571293 chr19 1028488
chr19:1028319-1028641 Island CNN2;CNN2 Body;Body 2.420525 -1.9882125 15.03487 0 0.0000455 9.196366
cg15592945 chr1 226736711
chr1:226736354-226737412 Island C1orf95 1stExon 1.738758 -5.7015620 13.92959 0 0.0000939 8.744156
cg17121412 chr15 22833400
chr15:22833282-22833911 Island TUBGCP5;TUBGCP5;TUBGCP5;TUBGCP5 1stExon;5’UTR;5’UTR;1stExon -1.746412 -4.7364030 -13.85471 0 0.0000939 8.711002
cg07762054 chr11 75273187
chr11:75272950-75273616 Island SERPINH1;SERPINH1 5’UTR;1stExon 1.789843 -5.0227036 13.78260 0 0.0000939 8.678756
cg14883700 chr2 173096328
chr2:173099482-173100098 N_Shelf -1.964926 3.6945378 -13.77481 0 0.0000939 8.675253
cg07220356 chr6 33263442
chr6:33266302-33267582 N_Shelf RGL2;RGL2 Body;Body -1.844915 3.5902570 -13.75131 0 0.0000939 8.664668
cg22157239 chr1 18972096
chr1:18971729-18972097 Island PAX7;PAX7;PAX7 Body;Body;Body 2.008880 -4.8748956 13.71440 0 0.0000939 8.647966
cg23422659 chr17 44929215
chr17:44928287-44929690 Island WNT9B Body -1.862788 -5.2099662 -13.70637 0 0.0000939 8.644320
cg04114460 chr18 29239388
OpenSea B4GALT6 Body -2.312635 2.8380218 -13.65741 0 0.0000939 8.622013
cg09669099 chr16 69419557
chr16:69419316-69420086 Island TERF2 1stExon 1.608697 -5.1869754 13.44182 0 0.0001095 8.521973
cg01957330 chr7 26904470
chr7:26903869-26904753 Island SKAP2 TSS200 1.999471 -4.4333864 13.35145 0 0.0001140 8.479147
cg04426862 chr14 37131008
chr14:37131181-37132785 N_Shore PAX9 5’UTR 1.655081 -4.4928784 13.20356 0 0.0001256 8.407898
cg26763618 chr16 14726469
chr16:14726516-14727366 N_Shore BFAR TSS200 1.907022 -4.4007952 13.04177 0 0.0001408 8.328258
cg26516701 chr12 53845683
chr12:53845668-53846606 Island PCBP2;PCBP2;PCBP2;PCBP2;PCBP2;PCBP2;PCBP2 TSS1500;TSS1500;TSS1500;TSS1500;TSS1500;TSS1500;TSS1500 1.502112 -4.5229391 12.94985 0 0.0001479 8.282206
cg13523875 chr7 117823728
chr7:117823984-117824485 N_Shore LSM8 TSS1500 1.508651 -5.5317694 12.88476 0 0.0001516 8.249235
cg21434530 chr4 4870486
chr4:4868440-4869173 S_Shore 1.814807 -4.8330389 12.78776 0 0.0001607 8.199539
cg20221020 chr11 67209257
OpenSea CORO1B;CORO1B Body;Body -1.974616 4.9962484 -12.63877 0 0.0001797 8.121881
cg02005873 chr6 10385839
chr6:10383525-10384114 S_Shore 1.833420 -1.2450526 12.55614 0 0.0001885 8.078111
cg20729301 chr2 233368115
chr2:233367780-233368577 Island 1.782130 -6.2407714 12.46508 0 0.0001998 8.029276
saveRDS(design_tw, "gu_design_motor.rds")

Regression analysis - within twin pair on diagnosis

#diagnosis:0,3.080605e-05
design_tw <- model.matrix(~ Twin_Group + diagnosis )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))
##        (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down        246805        2371        2834        2476        1954        4329
## NotSig       54218      786027      785932      786232      786824      783943
## Up          489635        2260        1892        1950        1880        2386
##        Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13 diagnosis
## Down         45182        2607        2861         1993         2546         0
## NotSig      733235      785493      785921       786617       786426    790658
## Up           12241        2558        1876         2048         1686         0
top <- topTable(fit2_tw,coef="diagnosis",num=Inf, sort.by = "P")
nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)
## [1] 16978
output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_guthrie_diagnosis.csv",row.names=FALSE)
head(output,30) %>% kbl(caption="diagnosis") %>% kable_paper("hover", full_width = F)
diagnosis
Name chr pos strand Islands_Name Relation_to_Island UCSC_RefGene_Name UCSC_RefGene_Group logFC AveExpr t P.Value adj.P.Val B
cg14507658 chr21 46408142
chr21:46407925-46408143 Island 0.5671208 3.3489041 6.094954 0.0000308 0.9905792 -1.590186
cg05281810 chr18 77542644
chr18:77542920-77543323 N_Shore -0.3448119 3.1050953 -6.089657 0.0000311 0.9905792 -1.592175
cg11355135 chr2 108602860
chr2:108602824-108603467 Island SLC5A7 TSS200 0.6328982 -3.8226413 5.874010 0.0000448 0.9905792 -1.675536
cg08008144 chr8 1983127
OpenSea -0.5206780 2.5786632 -5.716288 0.0000587 0.9905792 -1.739575
cg02344606 chr7 112062297
OpenSea IFRD1 TSS1500 0.4806032 -2.0898626 5.667407 0.0000639 0.9905792 -1.759966
cg18731055 chr13 49108131
chr13:49106411-49107463 S_Shore RCBTB2 TSS1500 0.5571933 0.0882235 5.572878 0.0000753 0.9905792 -1.800145
cg00633277 chr12 7250385
OpenSea C1RL;C1RL;C1RL Body;Body;Body -0.4834382 0.7795169 -5.473622 0.0000896 0.9905792 -1.843412
cg15740041 chr16 87781336
OpenSea KLHDC4;KLHDC4;KLHDC4 Body;Body;Body 0.6813155 1.6008461 5.469788 0.0000902 0.9905792 -1.845106
cg24267485 chr18 74766743
chr18:74770351-74770590 N_Shelf MBP;MBP Body;Body 0.5571721 -2.8069857 5.386618 0.0001045 0.9905792 -1.882264
cg14529256 chr1 245843607
OpenSea KIF26B Body 0.5967061 2.3027143 5.337372 0.0001141 0.9905792 -1.904644
cg25272125 chr1 31538096
chr1:31538412-31538938 N_Shore PUM1;PUM1 5’UTR;5’UTR 0.4303036 -5.3122575 5.324817 0.0001166 0.9905792 -1.910395
cg16239628 chr13 46626330
chr13:46626454-46627019 N_Shore ZC3H13 5’UTR -0.5717242 -4.6148948 -5.229213 0.0001384 0.9905792 -1.954801
cg11712265 chr19 7968106
chr19:7968325-7969421 N_Shore MAP2K7;MAP2K7;MAP2K7 TSS1500;TSS1500;TSS1500 0.3435839 -3.7499297 5.220197 0.0001406 0.9905792 -1.959045
cg27149150 chr1 7449851
OpenSea CAMTA1 Body -0.4871695 -0.3731033 -5.217092 0.0001414 0.9905792 -1.960509
cg23404427 chr19 1653402
OpenSea TCF3 TSS1500 -0.4921666 -5.0294197 -5.184932 0.0001499 0.9905792 -1.975739
cg04900486 chr1 151320117
chr1:151319326-151319545 S_Shore RFX5;RFX5 TSS1500;TSS1500 0.4810737 -4.3365457 5.163718 0.0001557 0.9905792 -1.985854
cg08277291 chr6 32037426
OpenSea TNXB Body -0.4216451 3.0373254 -5.158534 0.0001572 0.9905792 -1.988333
cg17904575 chr14 102290370
OpenSea PPP2R5C;PPP2R5C;PPP2R5C;PPP2R5C;PPP2R5C Body;Body;Body;Body;Body -0.3901626 1.0410297 -5.134116 0.0001642 0.9905792 -2.000058
cg05116656 chr5 94417892
OpenSea MCTP1;MCTP1;MCTP1 TSS1500;TSS1500;Body 0.4638695 -3.3750586 5.081264 0.0001807 0.9905792 -2.025685
cg08471713 chr17 41738893
OpenSea MEOX1;MEOX1;MEOX1 1stExon;1stExon;5’UTR -0.5408241 -2.6803562 -5.078389 0.0001817 0.9905792 -2.027089
cg01497262 chr2 43635730
OpenSea THADA;THADA;THADA Body;Body;Body 0.5232107 0.3214082 5.007192 0.0002068 0.9905792 -2.062179
cg06595448 chr7 50727315
OpenSea GRB10;GRB10;GRB10;GRB10 Body;Body;Body;Body 0.3863605 2.2935124 4.959104 0.0002258 0.9905792 -2.086235
cg26684767 chr2 231000057
OpenSea 0.4316435 -1.9244246 4.932767 0.0002370 0.9905792 -2.099533
cg25481455 chr6 22147551
OpenSea NBAT1;CASC15 TSS200;Body 0.4057296 -4.6757321 4.928824 0.0002387 0.9905792 -2.101531
cg25022405 chr22 26137668
OpenSea MYO18B TSS1500 0.3427957 3.6944563 4.919226 0.0002429 0.9905792 -2.106403
cg06753153 chr14 65186982
OpenSea PLEKHG3 5’UTR -0.3732891 3.1661977 -4.914021 0.0002453 0.9905792 -2.109050
cg11250279 chr13 24463401
chr13:24463119-24463677 Island C1QTNF9B-AS1;MIPEP;C1QTNF9B-AS1 5’UTR;1stExon;TSS200 0.3386238 -4.5577692 4.896740 0.0002532 0.9905792 -2.117864
cg21175642 chr3 48701770
chr3:48698335-48701667 S_Shore CELSR3 TSS1500 0.5201876 -3.2165700 4.887731 0.0002574 0.9905792 -2.122473
cg20492597 chr5 161274081
OpenSea GABRA1;GABRA1;GABRA1;GABRA1 TSS1500;TSS1500;TSS200;TSS1500 0.5847567 -2.6606980 4.885151 0.0002586 0.9905792 -2.123795
cg00010078 chr2 109967172
OpenSea SH3RF3 Body 0.4001143 -0.0677274 4.882349 0.0002600 0.9905792 -2.125231
saveRDS(design_tw, "gu_design_diagnosis.rds")

Regression analysis - within twin pair on SRS

#SRS:0,1.895198e-07
design_tw <- model.matrix(~ Twin_Group + SRS )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))
##        (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down        216305        2079        3241        1609        1985       54923
## NotSig      130496      786758      785600      787930      786968      727301
## Up          443857        1821        1817        1119        1705        8434
##        Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13    SRS
## Down        203038        2981        3165         1345         1041      0
## NotSig      553849      784700      785433       788210       789043 790658
## Up           33771        2977        2060         1103          574      0
top <- topTable(fit2_tw,coef="SRS",num=Inf, sort.by = "P")
head(top)
##                  logFC   AveExpr         t      P.Value adj.P.Val        B
## cg19720702 -0.03841315  3.071925 -9.525594 1.895198e-07 0.1498453 7.225728
## cg23884297  0.02680846  2.572387  7.478187 3.217331e-06 0.4586088 4.265063
## cg08281515 -0.04465879 -3.213167 -7.424466 3.488795e-06 0.4586088 4.180428
## cg15865511  0.02908008  2.876906  7.373374 3.769491e-06 0.4586088 4.099582
## cg03767531 -0.03515160 -4.926828 -6.952022 7.228624e-06 0.4586088 3.419673
## cg13534438 -0.02551836  2.702726 -6.877745 8.127541e-06 0.4586088 3.297352
nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)
## [1] 85217
output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_guthrie_SRS.csv",row.names=FALSE)
head(output,30) %>% kbl(caption="SRS score") %>% kable_paper("hover", full_width = F)
SRS score
Name chr pos strand Islands_Name Relation_to_Island UCSC_RefGene_Name UCSC_RefGene_Group logFC AveExpr t P.Value adj.P.Val B
cg19720702 chr1 53856531
OpenSea -0.0384131 3.0719249 -9.525594 2.00e-07 0.1498453 7.225728
cg23884297 chr20 30124811
OpenSea HM13;HM13;HM13;HM13 Body;Body;Body;Body 0.0268085 2.5723870 7.478187 3.20e-06 0.4586088 4.265063
cg08281515 chr11 111956885
chr11:111957352-111957681 N_Shore TIMM8B;SDHD;TIMM8B Body;TSS1500;Body -0.0446588 -3.2131669 -7.424466 3.50e-06 0.4586088 4.180428
cg15865511 chr2 43252424
OpenSea 0.0290801 2.8769063 7.373374 3.80e-06 0.4586088 4.099582
cg03767531 chr2 182545572
chr2:182547873-182549177 N_Shelf NEUROD1 TSS200 -0.0351516 -4.9268280 -6.952022 7.20e-06 0.4586088 3.419673
cg13534438 chr3 197229854
chr3:197225631-197226150 S_Shelf -0.0255184 2.7027261 -6.877745 8.10e-06 0.4586088 3.297352
cg08748969 chr12 69327779
chr12:69327021-69327532 S_Shore CPM;CPM;CPM TSS1500;TSS1500;5’UTR 0.0253835 -1.2268838 6.696612 1.09e-05 0.4586088 2.995932
cg19868495 chr13 61995249
OpenSea 0.0325808 0.8882802 6.667498 1.14e-05 0.4586088 2.947068
cg13206794 chr2 242010703
OpenSea SNED1 Body -0.0254861 2.0212841 -6.654190 1.16e-05 0.4586088 2.924695
cg03013237 chr16 1802436
OpenSea MAPK8IP3;MAPK8IP3 Body;Body -0.0339334 2.3972206 -6.523127 1.44e-05 0.4586088 2.703060
cg15616400 chr1 16062894
OpenSea SLC25A34;SLC25A34 5’UTR;1stExon -0.0347249 2.2495978 -6.514356 1.46e-05 0.4586088 2.688143
cg06874119 chr5 1892431
OpenSea CTD-2194D22.4 Body 0.0247034 -1.0871226 6.487037 1.52e-05 0.4586088 2.641619
cg17916960 chr15 79447300
OpenSea 0.0207799 -2.1616055 6.463602 1.58e-05 0.4586088 2.601628
cg01201782 chr16 1652449
OpenSea IFT140 Body -0.0237944 3.6792856 -6.460084 1.59e-05 0.4586088 2.595618
cg17497176 chr7 150763829
OpenSea SLC4A2 Body -0.0322948 2.6257672 -6.381261 1.81e-05 0.4586088 2.460518
cg21401292 chr6 134570669
OpenSea SGK1 Body 0.0263213 1.8937236 6.278702 2.15e-05 0.4586088 2.283470
cg22593006 chr6 170603682
chr6:170604411-170606479 N_Shore FAM120B;FAM120B TSS1500;Body 0.0302294 -0.7067929 6.253230 2.24e-05 0.4586088 2.239277
cg13385114 chr1 9656493
OpenSea TMEM201;TMEM201 Body;Body -0.0318501 2.3856163 -6.245681 2.27e-05 0.4586088 2.226162
cg15697902 chr19 14676213
chr19:14676201-14676542 Island TECR Body 0.0383106 3.4063584 6.206997 2.42e-05 0.4586088 2.158836
cg16134349 chr10 61900552
OpenSea ANK3;ANK3;ANK3;ANK3;ANK3 1stExon;5’UTR;Body;Body;Body -0.0194575 2.8615804 -6.178204 2.54e-05 0.4586088 2.108592
cg09451572 chr22 46760120
chr22:46759899-46760121 Island CELSR1 Body -0.0263675 2.3694795 -6.129153 2.76e-05 0.4586088 2.022739
cg14266530 chr13 52511468
OpenSea ATP7B;ATP7B;ATP7B Body;Body;Body -0.0285055 3.7038678 -6.128550 2.76e-05 0.4586088 2.021682
cg11920999 chr22 37417225
chr22:37414319-37416111 S_Shore MPST;MPST;TST;MPST 5’UTR;5’UTR;TSS1500;Body 0.0254667 3.2107679 6.122789 2.79e-05 0.4586088 2.011577
cg13709529 chr6 54001895
OpenSea C6orf142 Body 0.0351276 2.9963215 6.106824 2.87e-05 0.4586088 1.983550
cg01038640 chr5 157097845
chr5:157098263-157099041 N_Shore C5orf52 TSS1500 -0.0344627 -3.1684180 -6.100173 2.90e-05 0.4586088 1.971864
cg23657393 chr20 42790697
chr20:42788216-42789168 S_Shore JPH2 Body -0.0280709 2.5652689 -6.097772 2.91e-05 0.4586088 1.967644
cg09487134 chr10 35195634
OpenSea -0.0194752 0.2801082 -6.079980 3.00e-05 0.4586088 1.936345
cg00647512 chr8 42669499
OpenSea -0.0351548 2.2233507 -6.071817 3.04e-05 0.4586088 1.921973
cg18764121 chr15 25017456
chr15:25018174-25018533 N_Shore -0.0261292 0.1184601 -6.052347 3.14e-05 0.4586088 1.887653
cg17594278 chr12 112192763
OpenSea ACAD10;ACAD10 Body;Body -0.0281927 2.5565854 -6.043254 3.19e-05 0.4586088 1.871607
saveRDS(design_tw, "gu_design_srs.rds")

Regression analysis - within twin pair on inverse IQ

#iiq:0,1.326967e-05
design_tw <- model.matrix(~ Twin_Group + iiq )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))
##        (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down        247642        2515        3381        2813        1956       14775
## NotSig       47504      786029      785507      785809      786712      771811
## Up          495512        2114        1770        2036        1990        4072
##        Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13    iiq
## Down        182628        2926        2792         2019         3029      0
## NotSig      575220      784703      786190       786679       785893 790658
## Up           32810        3029        1676         1960         1736      0
top <- topTable(fit2_tw,coef="iiq",num=Inf, sort.by = "P")
head(top)
##                 logFC   AveExpr         t      P.Value adj.P.Val           B
## cg13470919  0.4700907 -3.017373  6.591035 1.326967e-05 0.3977294 -0.05845352
## cg09714307 -0.6303321  3.218485 -6.556696 1.402703e-05 0.3977294 -0.07677677
## cg08210568  0.4937255 -4.368887  6.459055 1.643964e-05 0.3977294 -0.12970089
## cg15934248  0.4688290 -4.799787  6.399469 1.812315e-05 0.3977294 -0.16260397
## cg03922146  0.4240616  2.534101  6.378132 1.876925e-05 0.3977294 -0.17449911
## cg26869412 -0.5128173  5.366690 -6.229826 2.398524e-05 0.3977294 -0.25885486
nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)
## [1] 99359
output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_guthrie_iIQ.csv",row.names=FALSE)
head(output,30) %>% kbl(caption="inverse IQ") %>% kable_paper("hover", full_width = F)
inverse IQ
Name chr pos strand Islands_Name Relation_to_Island UCSC_RefGene_Name UCSC_RefGene_Group logFC AveExpr t P.Value adj.P.Val B
cg13470919 chr12 108955320
chr12:108956188-108956759 N_Shore ISCU;ISCU;SART3 TSS1500;TSS1500;TSS200 0.4700907 -3.0173732 6.591035 1.33e-05 0.3977294 -0.0584535
cg09714307 chr5 10483440
OpenSea -0.6303321 3.2184848 -6.556697 1.40e-05 0.3977294 -0.0767768
cg08210568 chr3 15141216
chr3:15140035-15140890 S_Shore ZFYVE20 TSS1500 0.4937255 -4.3688871 6.459054 1.64e-05 0.3977294 -0.1297009
cg15934248 chr2 86332749
chr2:86332659-86333650 Island POLR1A;PTCD3 Body;TSS1500 0.4688290 -4.7997874 6.399469 1.81e-05 0.3977294 -0.1626040
cg03922146 chr6 30656499
chr6:30654392-30654934 S_Shore KIAA1949;NRM;KIAA1949 TSS1500;Body;TSS1500 0.4240616 2.5341006 6.378132 1.88e-05 0.3977294 -0.1744991
cg26869412 chr10 3023801
OpenSea -0.5128173 5.3666898 -6.229826 2.40e-05 0.3977294 -0.2588549
cg07241312 chr16 87051295
OpenSea 0.5246257 -1.2528635 6.078011 3.09e-05 0.3977294 -0.3483084
cg05845592 chr16 28634766
chr16:28634440-28635031 Island SULT1A1;SULT1A1 5’UTR;1stExon 0.6166146 -5.6410886 5.938857 3.92e-05 0.3977294 -0.4331377
cg15957908 chr6 158923478
OpenSea TULP4;TULP4 Body;Body -0.3848645 3.6984208 -5.887559 4.27e-05 0.3977294 -0.4651084
cg16336665 chr17 73028495
chr17:73030676-73031160 N_Shelf KCTD2 TSS200 0.7384996 -2.2320450 5.869774 4.41e-05 0.3977294 -0.4762820
cg07943548 chr15 86304357
chr15:86301949-86302213 S_Shelf KLHL25 3’UTR -0.4935287 4.2351942 -5.824888 4.76e-05 0.3977294 -0.5046867
cg13709580 chr5 10668678
OpenSea -0.3628510 2.2360219 -5.820586 4.79e-05 0.3977294 -0.5074243
cg23905944 chr20 62571879
chr20:62571738-62572556 Island UCKL1 Body -0.5878267 3.4256349 -5.807980 4.90e-05 0.3977294 -0.5154632
cg09375488 chr7 95115026
OpenSea ASB4;ASB4 TSS1500;TSS1500 -0.4619154 0.2983563 -5.806454 4.91e-05 0.3977294 -0.5164376
cg08848154 chr10 6278954
OpenSea -0.5018118 5.9889411 -5.805757 4.92e-05 0.3977294 -0.5168831
cg21419752 chr3 139514052
OpenSea -0.4099476 1.0313638 -5.766283 5.26e-05 0.3977294 -0.5422187
cg26362488 chr12 14924051
chr12:14922877-14923974 S_Shore HIST4H4;HIST4H4 1stExon;5’UTR 0.8214279 -4.6320444 5.765836 5.27e-05 0.3977294 -0.5425070
cg07690326 chr8 145019032
chr8:145018815-145019214 Island PLEC1;PLEC1;PLEC1;PLEC1;PLEC1;PLEC1 Body;Body;Body;TSS1500;Body;TSS200 -0.5396723 5.3537324 -5.753420 5.38e-05 0.3977294 -0.5505243
cg25849262 chr3 19189070
chr3:19189688-19190100 N_Shore KCNH8 TSS1500 0.4565090 -4.5611537 5.750801 5.41e-05 0.3977294 -0.5522186
cg14570121 chr19 45874153
chr19:45873339-45874033 S_Shore ERCC2;ERCC2 TSS1500;TSS1500 0.3972347 -4.7541665 5.725689 5.65e-05 0.3977294 -0.5685141
cg23644589 chr1 228404256
chr1:228404148-228404397 Island OBSCN;OBSCN Body;Body -0.6047790 3.8404119 -5.712356 5.78e-05 0.3977294 -0.5772043
cg27325462 chr11 129903118
chr11:129902811-129903119 Island -0.6266766 4.3295504 -5.698560 5.92e-05 0.3977294 -0.5862246
cg07019009 chr15 102265803
chr15:102264097-102264750 S_Shore TARSL2 TSS1500 -0.5618290 1.0241536 -5.697717 5.93e-05 0.3977294 -0.5867766
cg03873049 chr5 149625194
OpenSea CAMK2A;CAMK2A Body;Body -0.3784796 3.3104549 -5.692352 5.98e-05 0.3977294 -0.5902926
cg04465078 chr16 15766942
OpenSea NDE1;NDE1 Body;Body 0.9780949 -2.3163770 5.691853 5.99e-05 0.3977294 -0.5906199
cg02025435 chr16 89986317
chr16:89985913-89986356 Island MC1R 1stExon -0.5241210 5.4439103 -5.686550 6.04e-05 0.3977294 -0.5941000
cg19686637 chr15 85427743
OpenSea SLC28A1;SLC28A1 TSS200;TSS200 -0.3449407 3.2605496 -5.685110 6.06e-05 0.3977294 -0.5950455
cg19166616 chr20 62259876
chr20:62257589-62259476 S_Shore GMEB2 TSS1500 -0.4381218 -3.4166360 -5.684392 6.07e-05 0.3977294 -0.5955177
cg21645826 chr19 19313533
chr19:19314137-19314416 N_Shore NR2C2AP Body -0.4032165 -2.4565009 -5.673476 6.18e-05 0.3977294 -0.6026977
cg21353154 chr13 99381820
OpenSea SLC15A1 Body -0.4783244 -0.3243186 -5.654775 6.39e-05 0.3977294 -0.6150408
saveRDS(design_tw, "gu_design_iiq.rds")

Regression analysis - within twin pair on inverse language skills

#ilanguage:0,9.149182e-08
design_tw <- model.matrix(~ Twin_Group + ilanguage )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))
##        (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down        247342        1826        2269        2268        1972        5818
## NotSig       52060      787197      786886      786664      786825      782068
## Up          491256        1635        1503        1726        1861        2772
##        Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13 ilanguage
## Down        138499        2740        2456         1792         2752         0
## NotSig      624750      785140      786604       787138       786170    790658
## Up           27409        2778        1598         1728         1736         0
top <- topTable(fit2_tw,coef="ilanguage",num=Inf, sort.by = "P")
head(top)
##                 logFC   AveExpr          t      P.Value  adj.P.Val          B
## cg22165507 -0.7964738 -4.457186 -10.209898 9.149182e-08 0.07233874 -0.7254298
## cg13435526 -0.5169839  3.227245  -7.450780 3.624389e-06 0.99022496 -1.2323101
## cg25679269  0.9219619 -4.314010   6.921194 8.134347e-06 0.99022496 -1.3773828
## cg22533406 -0.6174444  4.882065  -6.864434 8.890001e-06 0.99022496 -1.3941709
## cg09476130 -0.4752016 -3.483859  -6.629915 1.289124e-05 0.99022496 -1.4662838
## cg18698884  0.7123973 -3.397628   6.556641 1.450054e-05 0.99022496 -1.4897535
nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)
## [1] 27881
output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_guthrie_ilanguage.csv",row.names=FALSE)
head(output,30) %>% kbl(caption="inverse language") %>% kable_paper("hover", full_width = F)
inverse language
Name chr pos strand Islands_Name Relation_to_Island UCSC_RefGene_Name UCSC_RefGene_Group logFC AveExpr t P.Value adj.P.Val B
cg22165507 chr1 165601007
chr1:165599563-165600574 S_Shore MGST3 5’UTR -0.7964738 -4.457186 -10.209898 1.00e-07 0.0723387 -0.7254298
cg13435526 chr8 22472050
OpenSea CCAR2;CCAR2 Body;Body -0.5169839 3.227245 -7.450780 3.60e-06 0.9902250 -1.2323101
cg25679269 chr11 118747550
OpenSea 0.9219619 -4.314011 6.921194 8.10e-06 0.9902250 -1.3773828
cg22533406 chr19 49106849
chr19:49106882-49107178 N_Shore FAM83E Body -0.6174444 4.882065 -6.864434 8.90e-06 0.9902250 -1.3941709
cg09476130 chr1 159870086
chr1:159869901-159870143 Island CCDC19 TSS200 -0.4752016 -3.483859 -6.629915 1.29e-05 0.9902250 -1.4662838
cg18698884 chr11 16632724
chr11:16632508-16632725 Island 0.7123973 -3.397628 6.556641 1.45e-05 0.9902250 -1.4897535
cg10021749 chr7 2557358
chr7:2558371-2559967 N_Shore LFNG;LFNG Body;TSS200 -0.5742813 2.799210 -6.470328 1.67e-05 0.9902250 -1.5179938
cg12843467 chr15 99190301
chr15:99190446-99194559 N_Shore IRAIN;IGF1R;IGF1R Body;TSS1500;TSS1500 -0.4936995 -4.971481 -6.364014 1.98e-05 0.9902250 -1.5536832
cg13938969 chr15 29409443
chr15:29407695-29408229 S_Shore APBA2;APBA2 3’UTR;3’UTR -0.5900858 2.762155 -6.346409 2.04e-05 0.9902250 -1.5596911
cg23404435 chr7 1588396
chr7:1590387-1591020 N_Shore TMEM184A Body -0.5840304 2.412321 -6.306655 2.18e-05 0.9902250 -1.5733623
cg13724550 chr5 142001033
OpenSea FGF1;FGF1;FGF1;FGF1;FGF1;FGF1;FGF1 5’UTR;TSS200;5’UTR;5’UTR;5’UTR;Body;Body -0.5416578 3.493387 -6.289624 2.24e-05 0.9902250 -1.5792638
cg20781140 chr2 175352216
chr2:175351364-175352154 S_Shore GPR155;GPR155 TSS1500;TSS1500 0.5757550 -4.026189 6.268960 2.32e-05 0.9902250 -1.5864597
cg19107578 chr5 493262
chr5:494803-495306 N_Shore SLC9A3 Body 0.5244209 3.194001 6.174053 2.71e-05 0.9902250 -1.6200247
cg03547220 chr19 7069348
chr19:7069396-7069921 N_Shore ZNF557;ZNF557;ZNF557 TSS200;TSS200;TSS200 0.5614201 -2.007877 6.150281 2.82e-05 0.9902250 -1.6285656
cg01369908 chr11 34375707
chr11:34378229-34379993 N_Shelf ABTB2 Body -0.4638439 3.189187 -6.143036 2.86e-05 0.9902250 -1.6311793
cg01231620 chr14 23298884
chr14:23298984-23299313 N_Shore MRPL52;MRPL52;MRPL52;MRPL52;MRPL52;MRPL52 TSS1500;TSS1500;TSS1500;TSS1500;TSS1500;TSS1500 0.5300644 -4.348508 6.073912 3.21e-05 0.9902250 -1.6563723
cg11913416 chr1 13909012
chr1:13910137-13910868 N_Shore PDPN;PDPN TSS1500;TSS1500 0.4995019 -0.696540 6.009515 3.57e-05 0.9902250 -1.6802627
cg02882857 chr5 149381215
chr5:149379964-149380802 S_Shore HMGXB3;TIGD6;TIGD6 5’UTR;TSS1500;TSS1500 0.6426011 -4.189717 5.997624 3.65e-05 0.9902250 -1.6847187
cg09219421 chr19 13779879
OpenSea -0.4059665 2.954756 -5.922746 4.14e-05 0.9902250 -1.7131055
cg27594834 chr19 11708137
chr19:11708277-11708683 N_Shore ZNF627 TSS200 -0.5832143 -4.512216 -5.851430 4.67e-05 0.9902250 -1.7406710
cg09459188 chr3 158520176
chr3:158519509-158520329 Island MFSD1;MFSD1 Body;Body 0.8657676 -4.775907 5.837325 4.79e-05 0.9902250 -1.7461847
cg03081615 chr4 8007156
OpenSea ABLIM2;ABLIM2;ABLIM2;ABLIM2;ABLIM2;ABLIM2 Body;Body;Body;Body;Body;Body -0.4584577 4.502980 -5.804816 5.06e-05 0.9902250 -1.7589717
cg27583671 chr6 151670624
OpenSea AKAP12;AKAP12 Body;Body -0.5434487 2.662483 -5.714514 5.91e-05 0.9902250 -1.7950719
cg17506328 chr11 46868210
chr11:46867194-46868077 S_Shore CKAP5;CKAP5 TSS1500;TSS1500 -0.3991032 -3.443528 -5.703255 6.03e-05 0.9902250 -1.7996333
cg27522733 chr1 193448591
OpenSea 0.5184504 -2.431826 5.680984 6.26e-05 0.9902250 -1.8086967
cg09733297 chr16 54155568
chr16:54155953-54156180 N_Shore -0.5840849 3.388928 -5.658058 6.52e-05 0.9902250 -1.8180817
cg17427349 chr15 38988004
OpenSea C15orf53 TSS1500 -0.4379188 3.008786 -5.618339 6.98e-05 0.9902250 -1.8344761
cg19572637 chr2 66660175
chr2:66660452-66660794 N_Shore 0.3265461 -4.483144 5.594407 7.28e-05 0.9902250 -1.8444370
cg00637015 chr1 31468126
OpenSea PUM1;PUM1 Body;Body 0.6517028 1.224266 5.593356 7.29e-05 0.9902250 -1.8448759
cg04263436 chr6 33245467
chr6:33244677-33245554 Island B3GALT4 1stExon 0.4051744 -5.021011 5.568245 7.62e-05 0.9902250 -1.8553980
saveRDS(design_tw, "gu_design_ilanguage.rds")

Converting to beta from M values

bDat = ilogit2(mvals)
bDat_new = getBeta(gmset_flt)
#View(bDat)
write.csv(bDat,file="ASD_guthrie_beta_onADOS_withintw_Nov27.csv",row.names=TRUE)

Plotting effect sizes of top DMPs

res <- read.csv("limma_guthrie_ADOS.csv")
TOPPROBES <- head(res$Name,9)
par(mfrow=c(3,3))
sapply(TOPPROBES,effect_plot)

## $cg12565057
## NULL
## 
## $cg00546248
## NULL
## 
## $cg15740041
## NULL
## 
## $cg05116656
## NULL
## 
## $cg11632273
## NULL
## 
## $cg01656620
## NULL
## 
## $cg10881749
## NULL
## 
## $cg13699963
## NULL
## 
## $cg22656126
## NULL
pdf("effect_guthrie.pdf")
par(mfrow=c(3,3))
sapply(TOPPROBES,effect_plot)
## $cg12565057
## NULL
## 
## $cg00546248
## NULL
## 
## $cg15740041
## NULL
## 
## $cg05116656
## NULL
## 
## $cg11632273
## NULL
## 
## $cg01656620
## NULL
## 
## $cg10881749
## NULL
## 
## $cg13699963
## NULL
## 
## $cg22656126
## NULL
dev.off()
## png 
##   2
par(mfrow=c(1,1))

Gene ontology of top DMPs LIMMA data

Gene Ontology analysis (gometh): top 1000 probes (limma data).

#res <- read.csv("ASD_guthrie_top_dmps_onADOS_withintw_Nov27_limma.csv", header = TRUE)
sigCpGs_1k = res$Name[1:1000]

sigCpGs_1k = as.character(sigCpGs_1k)
all = res$Name
length(all)

# kegg
gometh_kegg <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "KEGG", prior.prob=TRUE)
gometh_kegg <- subset(gometh_kegg,DE>2)
gometh_kegg$FDR <- p.adjust(gometh_kegg$P.DE)
gometh_kegg <- gometh_kegg[order(gometh_kegg$P.DE),]
head( gometh_kegg , 20)

# GO terms
gometh_go <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "GO" , prior.prob=TRUE)
gometh_go <- subset(gometh_go,DE>2)
gometh_go$FDR <- p.adjust(gometh_go$P.DE)
gometh_go <- gometh_go[order(gometh_go$P.DE),]
head( gometh_go , 20)

Now specifically analyse hypermethylated probes

res2 <- subset(res,logFC>0)
sigCpGs_1k = res2$Name[1:1000]
sigCpGs_1k = as.character(sigCpGs_1k)

# kegg
gometh_kegg <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "KEGG", prior.prob=TRUE)
gometh_kegg <- subset(gometh_kegg,DE>2)
gometh_kegg$FDR <- p.adjust(gometh_kegg$P.DE)
gometh_kegg <- gometh_kegg[order(gometh_kegg$P.DE),]
head( gometh_kegg , 20)

# GO terms
gometh_go <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "GO" , prior.prob=TRUE)
gometh_go <- subset(gometh_go,DE>2)
gometh_go$FDR <- p.adjust(gometh_go$P.DE)
gometh_go <- gometh_go[order(gometh_go$P.DE),]
head( gometh_go , 20)

Now specifically analyse hypomethylated probes

res2 <- subset(res,logFC<0)
sigCpGs_1k = res2$Name[1:1000]
sigCpGs_1k = as.character(sigCpGs_1k)

# kegg
gometh_kegg <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "KEGG", prior.prob=TRUE)
gometh_kegg <- subset(gometh_kegg,DE>2)
gometh_kegg$FDR <- p.adjust(gometh_kegg$P.DE)
gometh_kegg <- gometh_kegg[order(gometh_kegg$P.DE),]
head( gometh_kegg , 20)

# GO terms
gometh_go <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "GO" , prior.prob=TRUE)
gometh_go <- subset(gometh_go,DE>2)
gometh_go$FDR <- p.adjust(gometh_go$P.DE)
gometh_go <- gometh_go[order(gometh_go$P.DE),]
head( gometh_go , 20)

DMRCate - differentially methylated region analysis

#design matrix in regression 
design_tw <- model.matrix(~ADOS+Twin_Group)
design_tw_dmrc <- model.matrix(~ADOS)
#myannotation <- cpg.annotate("array", mDat, analysis.type="differential", design=design_tw, coef=2, fdr=1)
myannotation <- cpg.annotate("array", mvals, what = "M",
  arraytype = "EPIC", analysis.type="differential", design_tw, coef=2)
## Your contrast returned no individually significant probes. Try increasing the fdr. Alternatively, set pcutoff manually in dmrcate() to return DMRs, but be warned there is an increased risk of Type I errors.
dmrcoutput <- dmrcate(myannotation, lambda=1000, C=2, pcutoff=0.001)
## Fitting chr1...
## Fitting chr2...
## Fitting chr3...
## Fitting chr4...
## Fitting chr5...
## Fitting chr6...
## Fitting chr7...
## Fitting chr8...
## Fitting chr9...
## Fitting chr10...
## Fitting chr11...
## Fitting chr12...
## Fitting chr13...
## Fitting chr14...
## Fitting chr15...
## Fitting chr16...
## Fitting chr17...
## Fitting chr18...
## Fitting chr19...
## Fitting chr20...
## Fitting chr21...
## Fitting chr22...
## Demarcating regions...
## Done!
results.ranges <- data.frame(extractRanges(dmrcoutput, genome = "hg19") )
## snapshotDate(): 2022-10-31
## see ?DMRcatedata and browseVignettes('DMRcatedata') for documentation
## loading from cache
head(results.ranges, 20)
##    seqnames     start       end width strand no.cpgs min_smoothed_fdr Stouffer
## 1      chr1 234792390 234792752   363      *       2     8.755567e-07        1
## 2      chr2  46525649  46525678    30      *       2     9.543493e-04        1
## 3      chr6  91113506  91113687   182      *       2     8.564702e-05        1
## 4      chr9  96077410  96077420    11      *       2     9.668722e-04        1
## 5     chr11  44613980  44614255   276      *       2     4.084688e-05        1
## 6     chr22  20102216  20102580   365      *       3     1.418446e-06        1
## 7     chr22  31318147  31318444   298      *       7     9.447998e-04        1
## 8     chr12  75784541  75785295   755      *      11     8.243423e-09        1
## 9      chr2  33359198  33359688   491      *      10     9.097454e-04        1
## 10     chr7 100481990 100482960   971      *       5     3.223939e-08        1
## 11    chr16  87780848  87781501   654      *       5     5.201435e-07        1
## 12    chr11  65820391  65820516   126      *       4     4.325928e-04        1
## 13    chr12 124788627 124789036   410      *       4     2.317340e-05        1
##        HMFDR Fisher      maxdiff     meandiff     overlapping.genes
## 1  0.9999951      1  0.011843398  0.007349988           RP4-781K5.4
## 2  0.9999951      1 -0.003715122 -0.002901168                 EPAS1
## 3  0.9999951      1  0.015641676  0.011832067                  <NA>
## 4  0.9999951      1 -0.004274340 -0.001878717                  WNK2
## 5  0.9999951      1  0.017079330  0.015007099                  CD82
## 6  0.9999951      1 -0.012596545 -0.003140646                TRMT2A
## 7  0.9999951      1  0.009559386  0.003625967             MORC2-AS1
## 8  0.9999951      1  0.006045794  0.001967226       GLIPR1L2, CAPS2
## 9  0.9999951      1  0.015386270  0.008549527                 LTBP1
## 10 0.9999951      1  0.005702107  0.002285903                  SRRT
## 11 0.9999951      1  0.015904041  0.006726073 KLHDC4, RP11-278A23.2
## 12 0.9999951      1  0.006046969 -0.001162702                 SF3B2
## 13 0.9999951      1 -0.009740326 -0.002426397               FAM101A
write.csv(results.ranges, file = "dmrcoutput_guthrie.csv", row.names = TRUE)

Session information

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