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

Introduction

Here we are analysing buccal swab methylation in 6 twin pairs.

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

suppressPackageStartupMessages({
  library("missMethyl")
  library("limma")
  library("minfi")
  library("IlluminaHumanMethylation450kanno.ilmn12.hg19")
  library("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
  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")
  library("vioplot")
  library("RhpcBLASctl")
})

source("meth_functions.R")

RhpcBLASctl::blas_set_num_threads(1)

Load data

Load the annotation data and the Epic methylation data.

This analysis is to be conducted on Ubuntu with R4.

ann = getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.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 = "ILMLEPIC-16276_SampleSheet.csv")
## [read.metharray.sheet] Found the following CSV files:
## [1] "/asd_meth/cp_asd_data/ILMLEPIC-16276_SampleSheet.csv"
targets_gen <- targets_gen[grep("ASD",targets_gen$Sample_Name),]

# Remove 204375410103_R07C01 because it is corrupted
targets_gen <- targets_gen[grep("204375410103_R07C01",targets_gen$Basename, invert=TRUE),]

# Remove 204375410103_R08C01 beacuse its twin is missing
targets_gen <- targets_gen[grep("204375410103_R08C01",targets_gen$Basename, invert=TRUE),]

#targets$ID = paste(targets$Sample_Group,targets_gen$Sample_Name,sep=".")
rgSet = read.metharray.exp(targets = targets_gen)
sampleNames(rgSet) = targets_gen$SampleID

# patient data
ss <- read.table("ASD_EPIC_DATA/ASD_guthrie_blood_combined_summarysheet.csv",sep=",",header=TRUE)
ss <- ss[grep("G",ss$Sample_Name),]
ss$Sample_Name <- gsub("\\.","_",ss$Sample_Name)
ss$Sample_Name <- gsub("_G","",ss$Sample_Name)
ss <- ss[ss$Sample_Name %in% targets_gen$Sample_Name,]

ss$Sample_Name
##  [1] "ASD002_1" "ASD002_2" "ASD005_1" "ASD005_2" "ASD008_1" "ASD008_2"
##  [7] "ASD009_1" "ASD009_2" "ASD013_1" "ASD013_2"
targets_gen$Sample_Name
##  [1] "ASD002_1" "ASD002_2" "ASD005_1" "ASD005_2" "ASD008_1" "ASD008_2"
##  [7] "ASD009_1" "ASD009_2" "ASD013_1" "ASD013_2"
targets_gen <- ss

colnames(rgSet) <- targets_gen$Sample_Name

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_buccal_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")

Preprocessing

mset = preprocessRaw(rgSet)

Data exploration

Using Multi-dimensional scaling (MDS) plots before filtering.

mdsPlot(mset, sampGroups = targets_gen$Sample_Name,
  sampNames=targets_gen$Social_interaction_on_ADOS,legendPos="bottom")

mdsPlot(mset, sampGroups = targets_gen$Sex,
  sampNames=targets_gen$SampleID,legendPos="bottom")

Cell type composition analysis

Buccal swabs contain epithelial and leucocytes. Get the abundance of neutrophils as a proxy marker.

#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)

cells <- estimateCellCounts2(rgSet,
  referencePlatform= "IlluminaHumanMethylationEPIC",
  returnAll = TRUE)
## 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.
cells$prop
##            CD8T   CD4T     NK  Bcell   Mono    Neu
## ASD002_1 0.0000 0.4942 0.0000 0.0970 0.1416 0.2656
## ASD002_2 0.0000 0.4554 0.0000 0.0926 0.1213 0.3328
## ASD005_1 0.0000 0.5309 0.0000 0.0936 0.1592 0.1992
## ASD005_2 0.0000 0.4832 0.0000 0.0943 0.1365 0.2879
## ASD008_1 0.0000 0.4181 0.0230 0.0926 0.1233 0.1846
## ASD008_2 0.0844 0.2596 0.1168 0.1576 0.0022 0.1451
## ASD009_1 0.0000 0.4870 0.0000 0.1026 0.1357 0.2792
## ASD009_2 0.0000 0.3078 0.0000 0.0840 0.0757 0.5568
## ASD013_1 0.0000 0.5408 0.0000 0.0977 0.1555 0.1886
## ASD013_2 0.0000 0.4860 0.0000 0.1022 0.1257 0.2889

Filtering

Remove samples with more than 5% of probes with detection p-values greater than 0.01.

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_samples <- apply(detP,2, function(x) { length(which(x>0.01))/length(x) } ) < 0.05
detP <- detP[,keep_samples]
mset <- mset[,keep_samples]
targets_gen <- targets_gen[keep_samples,]

mset <- mset[rowMeans(detP) < 0.01,]
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,]

#getBeta
beta = getM(gmset_flt)
saveRDS(beta,"buccal_beta.rds")
df <- data.frame(t(t(colMeans(beta))))
colnames(df) = "gwam_bu"
write.table(df,file="buccal_gwam.tsv")

#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)

mvals <- mvals[apply(mvals,1,sd)!=0,]
mvals <- mvals[which(!is.na(rowSums(mvals))),]

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

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_buccal.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
design_tw <- model.matrix(~ Twin_Group + ADOS )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))
##        (Intercept) Twin_Group5 Twin_Group9 Twin_Group13   ADOS
## Down        332718         982        1025          892      0
## NotSig      169973      799215      799261       799480 801260
## Up          298569        1063         974          888      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] 37749
output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output, file="limma_buccal_ADOS.csv",row.names=FALSE)
output <- subset(output,P.Value<1e-4)
head(output,30) %>% kbl() %>% kable_paper("hover", full_width = F)
Name chr pos strand Islands_Name Relation_to_Island UCSC_RefGene_Name UCSC_RefGene_Group logFC AveExpr t P.Value adj.P.Val B
cg05143078 chr15 74988056
chr15:74988166-74988669 N_Shore EDC3;EDC3;EDC3 5’UTR;5’UTR;5’UTR 0.2755362 -3.788836 12.160630 8.40e-06 0.9047411 3.096494
cg15365028 chr2 236442817
chr2:236443815-236444025 N_Shore AGAP1;AGAP1 Body;Body -0.1870254 5.516570 -10.051787 2.82e-05 0.9047411 2.411940
cg10535017 chr17 61699755
chr17:61699007-61700728 Island MAP3K3;MAP3K3 TSS200;TSS200 0.1709695 -3.960242 9.504179 4.01e-05 0.9047411 2.190087
cg25208142 chr10 134611553
chr10:134611552-134613137 Island -0.1775579 1.344564 -9.409393 4.27e-05 0.9047411 2.149478
cg05616303 chr18 65321470
OpenSea LOC643542 Body -0.1770356 1.602827 -9.011285 5.60e-05 0.9047411 1.971274
cg06818937 chr19 50887645
chr19:50887181-50887933 Island POLD1;POLD1;POLD1;POLD1;POLD1 1stExon;1stExon;5’UTR;5’UTR;Body 0.1834566 -3.789042 8.749652 6.73e-05 0.9047411 1.847084
cg05075067 chr2 28113413
chr2:28113060-28113764 Island BRE;BRE;BRE;RBKS;BRE;LOC100302650;BRE TSS200;TSS200;TSS200;TSS200;TSS200;Body;TSS200 0.1785022 -5.457695 8.635742 7.30e-05 0.9047411 1.791165
cg01157070 chr16 56228511
chr16:56224731-56224980 S_Shelf DKFZP434H168;GNAO1;GNAO1 TSS200;Body;Body 0.1532936 -4.143216 8.601444 7.48e-05 0.9047411 1.774102
cg11534677 chr7 32932328
chr7:32930730-32931956 S_Shore KBTBD2 TSS1500 0.1794302 -5.468210 8.502515 8.03e-05 0.9047411 1.724290
cg23781719 chr6 43398988
chr6:43395188-43395894 S_Shelf ABCC10 TSS1500 -0.1495713 1.121903 -8.449884 8.35e-05 0.9047411 1.697425
cg14466709 chr9 95087590
chr9:95087949-95088606 N_Shore NOL8;NOL8;CENPP;NOL8 Body;5’UTR;TSS200;1stExon 0.2193702 -4.378564 8.428112 8.48e-05 0.9047411 1.686236
cg20554049 chr5 96064099
OpenSea CAST;CAST;CAST;CAST;CAST Body;Body;Body;Body;Body -0.1457575 1.808843 -8.314543 9.22e-05 0.9047411 1.627153
cg02514062 chr4 106629846
chr4:106629575-106630214 Island INTS12;INTS12;GSTCD;INTS12;INTS12 5’UTR;1stExon;TSS200;5’UTR;1stExon 0.3480929 -4.319272 8.226062 9.85e-05 0.9047411 1.580270 
saveRDS(design_tw, "buccal_design_ados.rds")
saveRDS(mvals, "buccal_mvals.rds")

Other models:

  • motor skills impairment

  • diagnosis

  • SRS

  • inverse IQ

  • inverse langguage

Regression analysis - within twin pair on motor skills

#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_Group5 Twin_Group9 Twin_Group13  motor
## Down        325839         386         389          401      0
## NotSig      170771      800420      800531       800494 801260
## Up          304650         454         340          365      0
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] 7186
output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_buccal_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
cg22703975 chr11 108765859
OpenSea DDX10 Body 2.490644 2.1414771 12.214028 0.0000050 1 -3.958841
cg14858551 chr20 18269244
chr20:18268544-18269405 Island ZNF133;ZNF133;ZNF133;ZNF133 1stExon;5’UTR;5’UTR;1stExon -2.163516 -4.3573914 -11.655568 0.0000068 1 -3.962402
cg01015483 chr4 186311121
OpenSea -1.996494 1.7240814 -11.360672 0.0000081 1 -3.964482
cg00870694 chr7 106775657
OpenSea PRKAR2B Body 1.869406 1.5565260 9.671466 0.0000240 1 -3.979909
cg01852341 chr18 46065250
chr18:46065061-46065263 Island CTIF;CTIF TSS200;TSS200 2.040689 -6.1585631 9.640040 0.0000245 1 -3.980267
cg24207248 chr17 78965708
chr17:78965416-78966484 Island CHMP6;CHMP6 1stExon;5’UTR 2.120551 -5.6823886 9.287217 0.0000314 1 -3.984512
cg10462093 chr12 125348448
chr12:125347784-125348688 Island SCARB1;SCARB1;SCARB1;SCARB1 1stExon;5’UTR;1stExon;5’UTR 1.542419 -3.9940376 9.123933 0.0000353 1 -3.986626
cg05293994 chr7 122526843
chr7:122525748-122526959 Island CADPS2;CADPS2;CADPS2 TSS200;TSS200;TSS200 1.840815 -5.6423722 9.062853 0.0000369 1 -3.987443
cg08984202 chr6 30421649
chr6:30418844-30419630 S_Shelf 2.901236 -4.6391852 8.957668 0.0000399 1 -3.988884
cg13846098 chr9 139720972
OpenSea RABL6;RABL6;RABL6 Body;Body;Body -1.634802 6.5903509 -8.858770 0.0000429 1 -3.990281
cg16860786 chr8 141542833
chr8:141545521-141545749 N_Shelf AGO2;AGO2 Body;Body -1.679968 1.5440340 -8.787735 0.0000452 1 -3.991310
cg13898527 chr11 82905656
chr11:82904478-82905660 Island ANKRD42;ANKRD42;ANKRD42;ANKRD42;ANKRD42;ANKRD42;ANKRD42;ANKRD42;ANKRD42;ANKRD42;ANKRD42;ANKRD42;ANKRD42 5’UTR;5’UTR;5’UTR;5’UTR;5’UTR;5’UTR;1stExon;1stExon;1stExon;1stExon;1stExon;1stExon;Body -1.640271 -5.1816482 -8.758695 0.0000462 1 -3.991737
cg01933251 chr12 16035274
chr12:16034931-16035836 Island STRAP TSS200 1.895662 -5.1276725 8.631727 0.0000509 1 -3.993648
cg23690488 chr2 114299609
chr2:114299035-114300177 Island -1.630336 -4.8256301 -8.618148 0.0000514 1 -3.993856
cg18347921 chr6 29975366
chr6:29974220-29975369 Island HLA-J;NCRNA00171 Body;Body -2.469043 -4.1263403 -8.588885 0.0000525 1 -3.994309
cg22118604 chr17 37353833
chr17:37353205-37353877 Island CACNB1;CACNB1;CACNB1;CACNB1;CACNB1;CACNB1 1stExon;5’UTR;5’UTR;1stExon;1stExon;5’UTR 1.966332 -3.8842252 8.577241 0.0000530 1 -3.994490
cg10938586 chr10 97849795
chr10:97849568-97850095 Island -1.617124 -4.2229001 -8.430095 0.0000594 1 -3.996838
cg15352224 chr3 10157326
chr3:10157129-10157612 Island C3orf10 TSS200 1.515095 -4.9683485 8.366541 0.0000624 1 -3.997886
cg18686665 chr2 629121
chr2:628909-629190 Island -1.619918 2.8460845 -8.338448 0.0000638 1 -3.998355
cg20552688 chr14 105956371
chr14:105956175-105958197 Island C14orf80;C14orf80;C14orf80 TSS1500;TSS1500;TSS200 -1.970191 -5.1523055 -8.299465 0.0000657 1 -3.999014
cg14949948 chr9 96215894
chr9:96213340-96214953 S_Shore FAM120AOS;FAM120A;FAM120A;FAM120A;FAM120A TSS200;Body;Body;Body;Body 1.925057 -4.4737027 8.294947 0.0000660 1 -3.999091
cg00685135 chr9 36136348
chr9:36136188-36137229 Island GLIPR2;GLIPR2;GLIPR2;GLIPR2;GLIPR2;GLIPR2;GLIPR2;GLIPR2;GLIPR2;GLIPR2;GLIPR2 TSS1500;TSS1500;TSS200;TSS200;TSS200;TSS200;TSS200;TSS200;TSS200;TSS200;TSS200 1.593536 -2.3932661 8.195157 0.0000714 1 -4.000816
cg13950210 chr1 45452319
chr1:45452037-45452444 Island EIF2B3;EIF2B3;EIF2B3;EIF2B3 1stExon;5’UTR;1stExon;5’UTR -1.392866 -5.8662807 -8.176482 0.0000725 1 -4.001144
cg17584296 chr10 75007026
chr10:75006629-75007184 Island DNAJC9 TSS200 2.186909 -5.9114782 8.092926 0.0000775 1 -4.002640
cg09293560 chr7 150068240
chr7:150068756-150070051 N_Shore REPIN1;REPIN1;REPIN1;REPIN1 TSS200;5’UTR;Body;5’UTR -1.680800 -0.4769369 -7.842289 0.0000950 1 -4.007373
cg01305596 chr16 16228299
OpenSea ABCC1;ABCC1;ABCC1;ABCC1;ABCC1 Body;Body;Body;Body;Body -2.240010 6.5922180 -7.809762 0.0000975 1 -4.008016
cg04440717 chr10 133938266
chr10:133937695-133938295 Island JAKMIP3 Body -1.845305 5.9518297 -7.740463 0.0001033 1 -4.009408
cg02267838 chr4 26862852
chr4:26862031-26863238 Island STIM2;STIM2;STIM2 1stExon;1stExon;1stExon 1.819857 -5.9190003 7.709248 0.0001060 1 -4.010046
cg23719290 chr17 56912132
OpenSea PPM1E;PPM1E Body;Body -2.323096 1.3137953 -7.697173 0.0001071 1 -4.010294
cg22851378 chr10 133747932
chr10:133751399-133751634 N_Shelf PPP2R2D;PPP2R2D;PPP2R2D TSS200;TSS200;TSS200 -2.172887 4.2620734 -7.580277 0.0001182 1 -4.012751 
saveRDS(design_tw, "buccal_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_Group5 Twin_Group9 Twin_Group13 diagnosis
## Down        344954         968         947         1017         0
## NotSig      132868      799225      799404       799249    801260
## Up          323438        1067         909          994         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] 43305
output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_buccal_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
cg03068039 chr8 146219928
OpenSea ZNF252;TMED10P Body;TSS1500 -0.9402151 1.5348067 -9.426553 0.0000378 0.7958854 -3.028003
cg05069807 chr4 6945702
OpenSea TBC1D14;TBC1D14 Body;Body -0.8338688 4.3534744 -8.939754 0.0000530 0.7958854 -3.046314
cg00958781 chr5 140018929
chr5:140019034-140019502 N_Shore TMCO6 TSS200 0.8662803 -5.6843175 8.587218 0.0000683 0.7958854 -3.061280
cg05536227 chr10 68771427
OpenSea LRRTM3;LOC101928961;LRRTM3;CTNNA3;CTNNA3 TSS1500;Body;Body;Body;Body -0.7915378 1.0721115 -8.086689 0.0000994 0.7958854 -3.085441
cg14303187 chr9 91722423
OpenSea SHC3 Body -0.7877157 1.7725928 -8.028744 0.0001040 0.7958854 -3.088485
cg05391486 chr1 201123919
chr1:201123245-201123746 S_Shore TMEM9 TSS1500 0.9579550 -4.5443969 8.024512 0.0001043 0.7958854 -3.088710
cg14466709 chr9 95087590
chr9:95087949-95088606 N_Shore NOL8;NOL8;CENPP;NOL8 Body;5’UTR;TSS200;1stExon 1.1634882 -4.3785636 7.897945 0.0001152 0.7958854 -3.095565
cg15365028 chr2 236442817
chr2:236443815-236444025 N_Shore AGAP1;AGAP1 Body;Body -0.9775296 5.5165699 -7.756175 0.0001289 0.7958854 -3.103576
cg20169823 chr12 116354837
chr12:116354787-116355187 Island 0.7301161 -4.5226081 7.547239 0.0001525 0.7958854 -3.116065
cg15265204 chr5 76326036
chr5:76326120-76326871 N_Shore AGGF1 TSS200 0.7669029 -4.3824796 7.462371 0.0001635 0.7958854 -3.121384
cg25208142 chr10 134611553
chr10:134611552-134613137 Island -0.9282856 1.3445644 -7.457083 0.0001642 0.7958854 -3.121721
cg24149295 chr17 34058557
chr17:34058549-34058778 Island RASL10B TSS200 0.6971515 -4.8031902 7.446120 0.0001657 0.7958854 -3.122420
cg23781719 chr6 43398988
chr6:43395188-43395894 S_Shelf ABCC10 TSS1500 -0.7885140 1.1219031 -7.445303 0.0001658 0.7958854 -3.122472
cg20554049 chr5 96064099
OpenSea CAST;CAST;CAST;CAST;CAST Body;Body;Body;Body;Body -0.7691136 1.8088429 -7.410573 0.0001706 0.7958854 -3.124704
cg05616303 chr18 65321470
OpenSea LOC643542 Body -0.9271867 1.6028268 -7.363175 0.0001775 0.7958854 -3.127792
cg24514137 chr16 72042695
chr16:72042347-72042888 Island DHODH Body 0.7540640 -4.4345296 7.359750 0.0001780 0.7958854 -3.128017
cg21233275 chr2 136634095
chr2:136633363-136634229 Island MCM6 TSS200 0.7593088 -5.2925296 7.303621 0.0001865 0.7958854 -3.131741
cg01216369 chr3 172859022
OpenSea SPATA16;SPATA16 1stExon;5’UTR -0.7423681 4.2262335 -7.302414 0.0001867 0.7958854 -3.131822
cg25185057 chr22 20073996
chr22:20073363-20073931 S_Shore DGCR8 Body -0.9151545 3.7513658 -7.279766 0.0001903 0.7958854 -3.133345
cg16641915 chr11 76778247
chr11:76777569-76778453 Island CAPN5 5’UTR 0.7327447 -4.9600350 7.253054 0.0001946 0.7958854 -3.135156
cg15520174 chr3 57234383
OpenSea HESX1 TSS200 -0.7008273 1.4913223 -7.221470 0.0001999 0.7958854 -3.137319
cg05143078 chr15 74988056
chr15:74988166-74988669 N_Shore EDC3;EDC3;EDC3 5’UTR;5’UTR;5’UTR 1.4138775 -3.7888359 7.186865 0.0002058 0.7958854 -3.139714
cg22037518 chr7 142491168
chr7:142494563-142495248 N_Shelf 0.7855042 -3.2927716 7.176212 0.0002077 0.7958854 -3.140457
cg22339221 chr12 49063107
OpenSea KANSL2 Body -0.7328616 1.4331794 -7.166741 0.0002094 0.7958854 -3.141120
cg13618906 chr4 109684291
chr4:109683729-109684362 Island AGXT2L1;AGXT2L1;AGXT2L1;AGXT2L1;AGXT2L1 TSS200;TSS200;TSS200;TSS200;TSS200 -0.6505212 -4.5837045 -7.133985 0.0002153 0.7958854 -3.143428
cg07939214 chr12 111872727
OpenSea SH2B3;SH2B3 1stExon;Body 0.7826665 -3.4072007 7.075174 0.0002264 0.7958854 -3.147638
cg19825037 chr11 30263038
OpenSea -0.7426653 1.1861244 -7.072070 0.0002270 0.7958854 -3.147862
cg02264416 chr11 104215127
OpenSea -0.6726719 0.7011332 -7.030571 0.0002353 0.7958854 -3.150886
cg22563390 chr19 55973338
chr19:55972697-55973442 Island ISOC2;ISOC2;ISOC2 TSS1500;TSS1500;TSS1500 0.9757795 -4.0096962 6.975951 0.0002467 0.7958854 -3.154932
cg07323885 chr3 137779756
OpenSea -0.6785586 1.0912169 -6.961496 0.0002498 0.7958854 -3.156016 
saveRDS(design_tw, "buccal_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_Group5 Twin_Group9 Twin_Group13    SRS
## Down        300937         693          69           95      0
## NotSig      299748      799875      801143       801080 801260
## Up          200575         692          48           85      0
top <- topTable(fit2_tw,coef="SRS",num=Inf, sort.by = "P")
head(top)
##                  logFC    AveExpr         t      P.Value adj.P.Val         B
## cg00893875  0.06855150 -2.4625349 10.790009 1.394412e-05 0.9400125 3.8550478
## cg09454438  0.04874571 -4.5061172  8.455495 6.783051e-05 0.9400125 2.1839588
## cg11613015  0.04392531 -4.4956212  7.075659 2.078179e-04 0.9400125 0.9794038
## cg14083017  0.04632029 -4.9332621  7.033319 2.156620e-04 0.9400125 0.9393432
## cg18302315 -0.03831760  0.8019458 -6.797351 2.659883e-04 0.9400125 0.7123816
## cg03041989  0.04241424 -5.8316602  6.726795 2.835133e-04 0.9400125 0.6432793
nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)
## [1] 20203
output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_buccal_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
cg00893875 chr10 91597593
chr10:91596974-91597792 Island 0.0685515 -2.4625349 10.790009 0.0000139 0.9400125 3.8550478
cg09454438 chr10 38383228
chr10:38383200-38383962 Island ZNF37A;ZNF37A TSS200;TSS200 0.0487457 -4.5061172 8.455495 0.0000678 0.9400125 2.1839588
cg11613015 chr5 112073406
OpenSea APC;APC;APC TSS200;TSS200;5’UTR 0.0439253 -4.4956212 7.075659 0.0002078 0.9400125 0.9794038
cg14083017 chr16 66914211
chr16:66914167-66914782 Island PDP2 TSS200 0.0463203 -4.9332621 7.033319 0.0002157 0.9400125 0.9393432
cg18302315 chr1 75046010
OpenSea ERICH3-AS1;ERICH3 Body;Body -0.0383176 0.8019458 -6.797351 0.0002660 0.9400125 0.7123816
cg03041989 chr16 30022548
chr16:30022651-30023007 N_Shore DOC2A TSS200 0.0424142 -5.8316602 6.726795 0.0002835 0.9400125 0.6432793
cg11860760 chr5 179448895
OpenSea RNF130 Body -0.0396291 3.7857502 -6.668426 0.0002990 0.9400125 0.5856734
cg07126263 chr18 57364862
chr18:57363714-57365001 Island CCBE1 TSS1500 0.0391097 -5.4896650 6.558316 0.0003309 0.9400125 0.4759057
cg12848864 chr11 108093218
chr11:108093211-108093969 Island ATM;NPAT TSS1500;Body 0.0401666 -4.4454585 6.537421 0.0003373 0.9400125 0.4549116
cg16039034 chr5 102148039
OpenSea -0.0394080 1.3829630 -6.495197 0.0003509 0.9400125 0.4123273
cg08529295 chr1 245132782
chr1:245133028-245134711 N_Shore EFCAB2;EFCAB2;EFCAB2;EFCAB2;EFCAB2 TSS1500;TSS1500;TSS1500;TSS1500;TSS1500 0.0360587 -2.0990154 6.442066 0.0003687 0.9400125 0.3584359
cg06360347 chr1 168023375
chr1:168025508-168025722 N_Shelf DCAF6;DCAF6 Body;Body -0.0371685 1.4357195 -6.344286 0.0004044 0.9400125 0.2583523
cg18122786 chr8 117778686
chr8:117778487-117779146 Island UTP23 TSS200 0.0374002 -6.1212670 6.278368 0.0004306 0.9400125 0.1902124
cg23018873 chr11 94501467
chr11:94501366-94502696 Island AMOTL1 TSS200 0.0413040 -5.4088566 6.257302 0.0004393 0.9400125 0.1683218
cg07861448 chr8 72757787
chr8:72755783-72756667 S_Shore MSC TSS1500 0.0368915 -3.4099017 6.199354 0.0004645 0.9400125 0.1078177
cg23302962 chr3 60112912
OpenSea FHIT;FHIT Body;Body -0.0369701 1.0869682 -6.192707 0.0004675 0.9400125 0.1008501
cg12474080 chr11 60929017
chr11:60928560-60929262 Island VPS37C TSS200 0.0353511 -3.4844951 6.183529 0.0004717 0.9400125 0.0912204
cg26234786 chr3 99575786
OpenSea FILIP1L;FILIP1L;FILIP1L;MIR548G;CMSS1;FILIP1L;FILIP1L 5’UTR;5’UTR;5’UTR;Body;Body;Body;Body -0.0381561 1.7421040 -6.142250 0.0004910 0.9400125 0.0477778
cg00691529 chr17 61904925
chr17:61904726-61905184 Island FTSJ3;FTSJ3;PSMC5;PSMC5 5’UTR;1stExon;TSS200;Body 0.0367419 -5.1801180 6.131896 0.0004959 0.9400125 0.0368475
cg06113789 chr6 43253097
chr6:43252734-43253476 Island TTBK1 Body 0.0361879 -5.1934779 6.119667 0.0005019 0.9400125 0.0239191
cg06538774 chr22 24584247
chr22:24584178-24584384 Island SUSD2 Body 0.0393077 6.4169419 6.081224 0.0005211 0.9400125 -0.0168465
cg01659632 chr6 144262297
OpenSea PLAGL1;PLAGL1;PLAGL1;PLAGL1;PLAGL1;PLAGL1;PLAGL1;PLAGL1 3’UTR;3’UTR;3’UTR;3’UTR;3’UTR;3’UTR;3’UTR;3’UTR -0.0418459 1.9423631 -6.067526 0.0005281 0.9400125 -0.0314179
cg05065230 chr3 49395807
chr3:49394855-49395942 Island GPX1;GPX1 TSS200;TSS200 -0.0417374 -4.5441581 -6.066474 0.0005286 0.9400125 -0.0325372
cg20738928 chr2 99871988
OpenSea LYG2 TSS1500 -0.0380426 1.5453416 -6.032863 0.0005464 0.9400125 -0.0683994
cg10390736 chr1 40420692
chr1:40420487-40421366 Island MFSD2A;MFSD2A TSS200;TSS200 0.0384756 -4.6046350 6.032020 0.0005468 0.9400125 -0.0693006
cg08628777 chr7 158050156
OpenSea PTPRN2;PTPRN2;PTPRN2 Body;Body;Body -0.0410000 5.0853243 -5.932697 0.0006034 0.9400125 -0.1761469
cg26474110 chr5 96143884
chr5:96143102-96143885 Island ERAP1;ERAP1;ERAP1;ERAP1;ERAP1 1stExon;1stExon;5’UTR;5’UTR;5’UTR 0.0476642 -4.2126886 5.919062 0.0006116 0.9400125 -0.1909159
cg22900057 chr10 134481296
chr10:134477298-134477515 S_Shelf INPP5A Body -0.0327685 1.6150439 -5.909406 0.0006176 0.9400125 -0.2013890
cg13455700 chr11 122753504
chr11:122753386-122753951 Island C11orf63;C11orf63;C11orf63;C11orf63 1stExon;5’UTR;5’UTR;1stExon 0.0340082 -5.3274470 5.900140 0.0006233 0.9400125 -0.2114518
cg00183355 chr21 34100294
chr21:34099716-34100796 Island SYNJ1;SYNJ1;SYNJ1;SYNJ1 TSS200;1stExon;TSS200;1stExon 0.0360371 -5.1366017 5.897747 0.0006248 0.9400125 -0.2140522 
saveRDS(design_tw, "buccal_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_Group5 Twin_Group9 Twin_Group13    iiq
## Down        350438         759        1207         1054      0
## NotSig      103618      799700      798902       799203 801260
## Up          347204         801        1151         1003      0
top <- topTable(fit2_tw,coef="iiq",num=Inf, sort.by = "P")
head(top)
##                logFC   AveExpr         t      P.Value adj.P.Val        B
## cg24734172  1.438352 -5.373293  14.70365 3.908877e-06 0.3949487 1.947570
## cg01692340  1.567375 -5.942862  12.99056 8.380591e-06 0.3949487 1.743991
## cg26579311  1.511178 -6.097256  12.69117 9.668200e-06 0.3949487 1.701574
## cg15640339 -1.796122  6.003061 -12.63942 9.913214e-06 0.3949487 1.694001
## cg15039313  1.055687 -5.246787  11.56420 1.705723e-05 0.3949487 1.518779
## cg01762070  1.023393 -3.591630  11.45904 1.803170e-05 0.3949487 1.499621
nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)
## [1] 84921
output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_buccal_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
cg24734172 chr22 37172326
chr22:37171766-37172669 Island IFT27;IFT27;IFT27 TSS200;TSS200;TSS200 1.4383516 -5.3732928 14.703649 3.90e-06 0.3949487 1.947570
cg01692340 chr7 100318319
chr7:100318106-100318684 Island EPO TSS200 1.5673746 -5.9428619 12.990558 8.40e-06 0.3949487 1.743991
cg26579311 chr8 119124103
chr8:119123974-119124432 Island EXT1 TSS200 1.5111779 -6.0972559 12.691175 9.70e-06 0.3949487 1.701574
cg15640339 chr8 2017508
OpenSea MYOM2 Body -1.7961220 6.0030614 -12.639419 9.90e-06 0.3949487 1.694001
cg15039313 chr9 132565418
chr9:132565316-132566033 Island TOR1B TSS200 1.0556873 -5.2467874 11.564204 1.71e-05 0.3949487 1.518779
cg01762070 chr14 21101007
chr14:21100838-21101043 Island 1.0233934 -3.5916299 11.459044 1.80e-05 0.3949487 1.499621
cg11025641 chr12 104324133
chr12:104323017-104324811 Island TTC41P;MIR3652;HSP90B1;HSP90B1 TSS200;TSS200;5’UTR;1stExon 1.0557674 -4.3220337 11.246323 2.02e-05 0.3949487 1.459640
cg27446054 chr8 146203972
OpenSea ZNF252P Body -1.0161387 1.3552784 -10.732521 2.68e-05 0.3949487 1.355817
cg26235537 chr10 129924700
chr10:129923813-129924962 Island MKI67;MKI67 TSS1500;TSS1500 0.9474233 -3.4469476 10.719992 2.70e-05 0.3949487 1.353150
cg16652621 chr1 27614186
OpenSea WDTC1;WDTC1 Body;Body -1.1110180 4.4124362 -10.536645 3.00e-05 0.3949487 1.313325
cg24679395 chr11 45687275
chr11:45686160-45687495 Island CHST1 TSS200 1.1244319 -6.2870159 10.356229 3.33e-05 0.3949487 1.272656
cg12352347 chr13 46785668
OpenSea 1.0764194 -3.3964482 10.136828 3.78e-05 0.3949487 1.221122
cg03157605 chr6 30685409
chr6:30684836-30685503 Island MDC1;MDC1 1stExon;5’UTR 1.0852189 -4.0161674 10.077081 3.92e-05 0.3949487 1.206678
cg05701691 chr19 18682687
chr19:18682374-18683034 Island UBA52;UBA52;UBA52;UBA52 1stExon;1stExon;5’UTR;5’UTR 1.0050193 -5.6639318 10.060190 3.96e-05 0.3949487 1.202562
cg14898116 chr10 17272764
chr10:17270430-17272617 S_Shore VIM Body 1.0318489 -3.0539520 10.039893 4.01e-05 0.3949487 1.197597
cg02301929 chr22 51066512
chr22:51066107-51067051 Island ARSA;ARSA;ARSA;ARSA;ARSA;ARSA;ARSA;ARSA;ARSA;ARSA 1stExon;5’UTR;5’UTR;1stExon;5’UTR;1stExon;5’UTR;1stExon;1stExon;5’UTR 0.9454864 -5.1998520 9.967915 4.19e-05 0.3949487 1.179818
cg18144544 chr7 6197890
chr7:6198495-6198731 N_Shore USP42 Body -0.9331595 4.1019710 -9.871836 4.44e-05 0.3949487 1.155666
cg00505936 chr16 69139418
chr16:69139628-69141900 N_Shore HAS3 TSS1500 0.9580267 -3.4977151 9.840330 4.52e-05 0.3949487 1.147639
cg12128740 chr1 62207888
chr1:62207827-62208818 Island INADL TSS1500 0.8817804 -3.6754241 9.784752 4.68e-05 0.3949487 1.133350
cg26651784 chr16 67062934
chr16:67062613-67064012 Island CBFB;CBFB TSS200;TSS200 0.9311122 -5.5471283 9.741141 4.81e-05 0.3949487 1.122019
cg01842890 chr6 30292264
chr6:30294169-30295071 N_Shore HCG18;HCG18 Body;Body -1.0729160 2.6099062 -9.727309 4.85e-05 0.3949487 1.118404
cg17234198 chr20 47894957
chr20:47894197-47895199 Island C20orf199;SNORD12C;C20orf199;ZNFX1;C20orf199 Body;TSS1500;TSS1500;TSS1500;TSS200 1.1085744 -5.1530924 9.631370 5.14e-05 0.3949487 1.093034
cg10176510 chr6 97337802
OpenSea NDUFAF4 3’UTR -0.9143135 1.4237125 -9.527062 5.49e-05 0.3949487 1.064858
cg02498477 chr18 44471411
OpenSea PIAS2;PIAS2 Body;Body -0.8727067 0.7771203 -9.511021 5.54e-05 0.3949487 1.060470
cg15148879 chr11 47788577
chr11:47788200-47789022 Island FNBP4 Body -0.9603499 -5.3239285 -9.506272 5.56e-05 0.3949487 1.059168
cg13760252 chr19 5455500
chr19:5455481-5456262 Island ZNRF4;ZNRF4 5’UTR;1stExon -0.8639407 4.3884380 -9.403198 5.93e-05 0.3949487 1.030577
cg02605007 chr19 10541304
chr19:10542662-10543120 N_Shore PDE4A;PDE4A Body;TSS1500 1.6812926 -3.6873959 9.386516 6.00e-05 0.3949487 1.025890
cg07058528 chr17 6918698
chr17:6917569-6918539 S_Shore C17orf49;C17orf49;C17orf49 Body;Body;Body 0.8741340 -4.7528918 9.346481 6.15e-05 0.3949487 1.014574
cg05375741 chr17 40897145
chr17:40896658-40897234 Island EZH1 TSS200 1.0259873 -5.2611217 9.343781 6.16e-05 0.3949487 1.013807
cg21422208 chr3 53164667
chr3:53164206-53164698 Island RFT1 TSS200 0.8693483 -4.6939422 9.335322 6.20e-05 0.3949487 1.011402 
saveRDS(design_tw, "buccal_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_Group5 Twin_Group9 Twin_Group13 ilanguage
## Down        341316         522         966          627         0
## NotSig      124052      800218      799396       800085    801260
## Up          335892         520         898          548         0
top <- topTable(fit2_tw,coef="ilanguage",num=Inf, sort.by = "P")
head(top)
##                logFC   AveExpr         t      P.Value adj.P.Val         B
## cg14982506 -1.534375  2.453888 -12.46237 5.144798e-06 0.9920423 -3.292601
## cg26062048  1.511903 -5.160393  11.91519 6.946143e-06 0.9920423 -3.299917
## cg12654743  1.351536 -4.562291  11.20733 1.044224e-05 0.9920423 -3.310881
## cg13944275 -1.283127 -5.041450 -10.66305 1.452019e-05 0.9920423 -3.320698
## cg05351078  1.243993 -5.553076  10.61673 1.494370e-05 0.9920423 -3.321597
## cg03068874 -1.260082  1.653308 -10.16936 1.984429e-05 0.9920423 -3.330859
nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)
## [1] 18026
output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_buccal_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
cg14982506 chr10 24832948
OpenSea KIAA1217;KIAA1217;KIAA1217 Body;Body;Body -1.5343752 2.453888 -12.462372 0.0000051 0.9920423 -3.292601
cg26062048 chr5 153418412
chr5:153418183-153418669 Island MFAP3;MFAP3;FAM114A2;MFAP3 TSS200;TSS200;5’UTR;TSS200 1.5119026 -5.160393 11.915188 0.0000069 0.9920423 -3.299917
cg12654743 chr17 8023707
chr17:8023960-8024239 N_Shore ALOXE3 TSS1500 1.3515362 -4.562291 11.207331 0.0000104 0.9920423 -3.310881
cg13944275 chr9 140024112
chr9:140023747-140025010 Island -1.2831266 -5.041450 -10.663051 0.0000145 0.9920423 -3.320698
cg05351078 chr10 45455109
chr10:45455239-45455474 N_Shore RASSF4 TSS200 1.2439930 -5.553076 10.616735 0.0000149 0.9920423 -3.321597
cg03068874 chr7 158152212
OpenSea PTPRN2;PTPRN2;PTPRN2 Body;Body;Body -1.2600819 1.653308 -10.169364 0.0000198 0.9920423 -3.330859
cg20007021 chr14 24780404
chr14:24779874-24780932 Island LTB4R2;LTB4R2;LTB4R2;LTB4R;CIDEB;CIDEB Body;1stExon;Body;TSS1500;1stExon;5’UTR -1.4387543 -2.815103 -10.071784 0.0000211 0.9920423 -3.333028
cg06149604 chr1 113933334
chr1:113932580-113934070 Island MAGI3;MAGI3 TSS200;TSS200 1.3941705 -5.706267 9.991134 0.0000223 0.9920423 -3.334863
cg05845592 chr16 28634766
chr16:28634440-28635031 Island SULT1A1;SULT1A1 5’UTR;1stExon -1.2684225 -4.900838 -9.731058 0.0000265 0.9920423 -3.341062
cg24583179 chr5 180229965
chr5:180229375-180230147 Island MGAT1;MGAT1;MGAT1;MGAT1;MGAT1;MGAT1 1stExon;5’UTR;5’UTR;5’UTR;5’UTR;5’UTR 1.1928096 -4.782711 9.209242 0.0000379 0.9920423 -3.354924
cg07198656 chr10 48952399
chr10:48952215-48952880 Island GLUD1P7;BMS1P5 TSS200;Body -1.8564093 -6.386467 -9.131953 0.0000401 0.9920423 -3.357157
cg06240216 chr6 129963901
OpenSea ARHGAP18 Body -1.6454080 1.801600 -9.034550 0.0000430 0.9920423 -3.360042
cg17752088 chr5 78810367
chr5:78809347-78810468 Island HOMER1 TSS1500 1.4223703 -6.066901 8.912833 0.0000469 0.9920423 -3.363763
cg10573158 chr6 30689496
chr6:30688624-30689530 Island TUBB Body 1.0718563 -5.179513 8.805741 0.0000507 0.9920423 -3.367149
cg10136354 chr2 73151457
chr2:73151200-73152060 Island EMX1 Body 1.0334153 -4.314413 8.452436 0.0000660 0.9920423 -3.379115
cg02687055 chr2 24232865
chr2:24232680-24233260 Island MFSD2B TSS200 1.1593473 -6.197381 8.435392 0.0000668 0.9920423 -3.379726
cg17147121 chr19 52674949
chr19:52674327-52674855 S_Shore ZNF836 TSS200 0.9938850 -3.180525 8.272116 0.0000757 0.9920423 -3.385734
cg02991063 chr10 120938274
chr10:120938028-120938418 Island PRDX3;PRDX3 1stExon;1stExon 1.1109398 -5.612401 8.271481 0.0000758 0.9920423 -3.385758
cg19145477 chr17 73348462
OpenSea GRB2;GRB2 Body;Body 0.9682111 -4.329982 8.232074 0.0000781 0.9920423 -3.387254
cg01983185 chr19 8428990
chr19:8428965-8429492 Island ANGPTL4;ANGPTL4;ANGPTL4 TSS200;TSS200;TSS200 1.0404459 -4.243620 8.092829 0.0000871 0.9920423 -3.392690
cg17732216 chr6 167178280
OpenSea RPS6KA2 Body -1.0868877 3.323186 -8.062800 0.0000892 0.9920423 -3.393893
cg09499699 chr13 80146405
OpenSea -1.3181366 2.253812 -8.040193 0.0000908 0.9920423 -3.394807
cg24204556 chr22 22222030
chr22:22221055-22222367 Island MAPK1;MAPK1 TSS200;TSS200 1.2279562 -5.233175 8.040000 0.0000908 0.9920423 -3.394815
cg00158122 chr10 101290029
chr10:101290025-101290338 Island 0.9875524 -1.614181 8.020809 0.0000922 0.9920423 -3.395595
cg03342084 chr11 43902343
chr11:43902254-43902528 Island ALKBH3 TSS200 1.3664376 -5.255032 7.959508 0.0000968 0.9920423 -3.398121
cg06924657 chr7 150381630
OpenSea GIMAP2 TSS1500 -1.1076854 1.570041 -7.950493 0.0000975 0.9920423 -3.398497
cg21846876 chr2 211341373
chr2:211341181-211341536 Island LANCL1;LANCL1;LANCL1;CPS1;LANCL1;LANCL1 1stExon;5’UTR;5’UTR;TSS1500;1stExon;5’UTR 0.9942232 -4.872929 7.841211 0.0001065 0.9920423 -3.403136
cg15081882 chr6 31588692
chr6:31587779-31589024 Island BAT2 5’UTR 1.0444103 -5.953226 7.806689 0.0001095 0.9920423 -3.404636
cg19450240 chr16 74734251
chr16:74734152-74734572 Island MLKL;MLKL 5’UTR;5’UTR -0.9404235 -4.974311 -7.775852 0.0001123 0.9920423 -3.405989
cg14490783 chr17 39184887
OpenSea KRTAP1-5 TSS1500 -0.9283875 2.577990 -7.775659 0.0001123 0.9920423 -3.405998 
saveRDS(design_tw, "buccal_design_ilanguage.rds")

Converting to beta from M values

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

Genome wide average methylation

vioplot(bDat)

#gwam <- colMeans(bDat)
gwam <- apply(bDat,2,median)

message("GWAM assocation with ADOS")
## GWAM assocation with ADOS
mylm <- lm(gwam~ADOS)
summary(mylm)
## 
## Call:
## lm(formula = gwam ~ ADOS)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0105251 -0.0033729 -0.0005064  0.0019995  0.0180868 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.5101155  0.0074943  68.067 6.76e-10 ***
## ADOS        0.0004096  0.0007353   0.557    0.598    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.009442 on 6 degrees of freedom
## Multiple R-squared:  0.04918,    Adjusted R-squared:  -0.1093 
## F-statistic: 0.3104 on 1 and 6 DF,  p-value: 0.5976
message("GWAM assocation with SRS")
## GWAM assocation with SRS
mylm <- lm(gwam~SRS)
summary(mylm)
## 
## Call:
## lm(formula = gwam ~ SRS)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.007020 -0.005317 -0.003939  0.003673  0.015638 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.5066154  0.0065454   77.40 3.13e-10 ***
## SRS         0.0001012  0.0000810    1.25    0.258    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.008625 on 6 degrees of freedom
## Multiple R-squared:  0.2065, Adjusted R-squared:  0.0743 
## F-statistic: 1.562 on 1 and 6 DF,  p-value: 0.2579
message("GWAM assocation with motor impairment")
## GWAM assocation with motor impairment
mylm <- lm(gwam~motor)
summary(mylm)
## 
## Call:
## lm(formula = gwam ~ motor)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0145370 -0.0028985  0.0005501  0.0032915  0.0157135 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.5145370  0.0043063 119.485 2.32e-11 ***
## motor       -0.0009113  0.0035161  -0.259    0.804    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.009629 on 6 degrees of freedom
## Multiple R-squared:  0.01107,    Adjusted R-squared:  -0.1538 
## F-statistic: 0.06717 on 1 and 6 DF,  p-value: 0.8042
message("GWAM assocation with diagnosis")
## GWAM assocation with diagnosis
mylm <- lm(gwam~diagnosis)
summary(mylm)
## 
## Call:
## lm(formula = gwam ~ diagnosis)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.013759 -0.003534 -0.000589  0.003880  0.016492 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.138e-01  5.229e-03  98.250 7.49e-11 ***
## diagnosis   9.496e-05  3.953e-03   0.024    0.982    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.009682 on 6 degrees of freedom
## Multiple R-squared:  9.617e-05,  Adjusted R-squared:  -0.1666 
## F-statistic: 0.0005771 on 1 and 6 DF,  p-value: 0.9816
message("GWAM assocation with iIQ")
## GWAM assocation with iIQ
mylm <- lm(gwam~iiq)
summary(mylm)
## 
## Call:
## lm(formula = gwam ~ iiq)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.012389 -0.002407 -0.000714  0.001730  0.017685 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.513854   0.003270  157.16 4.48e-12 ***
## iiq         0.002657   0.003495    0.76    0.476    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.009248 on 6 degrees of freedom
## Multiple R-squared:  0.08785,    Adjusted R-squared:  -0.06418 
## F-statistic: 0.5779 on 1 and 6 DF,  p-value: 0.476
plot(gwam~iiq)
abline(mylm)

message("GWAM assocation with ilanguage")
## GWAM assocation with ilanguage
mylm <- lm(gwam~ilanguage)
summary(mylm)
## 
## Call:
## lm(formula = gwam ~ ilanguage)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0094937 -0.0042126 -0.0003427  0.0020640  0.0152664 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.513854   0.002816 182.468 1.83e-12 ***
## ilanguage   0.005097   0.003011   1.693    0.141    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.007965 on 6 degrees of freedom
## Multiple R-squared:  0.3233, Adjusted R-squared:  0.2105 
## F-statistic: 2.867 on 1 and 6 DF,  p-value: 0.1414
plot(gwam~ilanguage)
abline(mylm)

probes_body <- rownames(ann_sub[grep("Body",ann_sub$UCSC_RefGene_Group),])
bDat_body <- bDat[which(rownames(bDat) %in% probes_body),]
gwam <- apply(bDat_body,2,median)

probes_body <- rownames(ann_sub[grep("Island",ann_sub$Relation_to_Island),])
bDat_body <- bDat[which(rownames(bDat) %in% probes_body),]
gwam <- apply(bDat_body,2,median)

Plotting effect sizes of top DMPs

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

## $cg05143078
## NULL
## 
## $cg15365028
## NULL
## 
## $cg10535017
## NULL
## 
## $cg25208142
## NULL
## 
## $cg05616303
## NULL
## 
## $cg06818937
## NULL
## 
## $cg05075067
## NULL
## 
## $cg01157070
## NULL
## 
## $cg11534677
## NULL
pdf("effect_buccal.pdf")
par(mfrow=c(3,3))
sapply(TOPPROBES,effect_plot)
## $cg05143078
## NULL
## 
## $cg15365028
## NULL
## 
## $cg10535017
## NULL
## 
## $cg25208142
## NULL
## 
## $cg05616303
## NULL
## 
## $cg06818937
## NULL
## 
## $cg05075067
## NULL
## 
## $cg01157070
## NULL
## 
## $cg11534677
## 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).

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 = "EPICv1", analysis.type="differential", design=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.
## Warning in S4Vectors:::normarg_mcols(mcols, Class, ans_len): You supplied metadata columns of length 801260 to set on an object of
##   length 802647. However please note that the latter is not a multiple of
##   the former.
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") )
## see ?DMRcatedata and browseVignettes('DMRcatedata') for documentation
## loading from cache
head(results.ranges, 20)
##    seqnames     start       end width strand no.cpgs min_smoothed_fdr
## 1     chr21  46238957  46239952   996      *       3     8.894480e-39
## 2     chr10  45469396  45470562  1167      *       5     3.494039e-20
## 3      chr8 144810854 144811385   532      *       3     1.593597e-14
## 4     chr11  64135217  64136263  1047      *      13     7.898096e-14
## 5      chr1  85405081  85405313   233      *       3     1.219110e-13
## 6     chr12  69198743  69199037   295      *       2     1.906223e-13
## 7      chr4    818664    819151   488      *       2     4.003904e-13
## 8     chr11  61716821  61717780   960      *      10     1.055967e-12
## 9      chr4  40565477  40566405   929      *       3     1.055967e-12
## 10     chr3 115717990 115718524   535      *       2     1.055967e-12
## 11     chr2  10223562  10224752  1191      *       7     2.599837e-12
## 12    chr11 105892558 105893693  1136      *      16     3.225752e-11
## 13     chr6  32030160  32030554   395      *       2     1.414488e-10
## 14    chr19   1904554   1905739  1186      *      10     1.649791e-10
## 15     chr6  30847733  30848178   446      *       2     1.986771e-10
## 16    chr12   4140243   4141193   951      *       6     2.267371e-10
## 17     chr4 110880215 110880666   452      *       2     3.533034e-10
## 18    chr17  43210053  43211248  1196      *       7     4.058657e-10
## 19    chr11  77532707  77532772    66      *       3     7.154523e-10
## 20     chr3 178253375 178254281   907      *       7     9.325203e-10
##       Stouffer       HMFDR      Fisher      maxdiff      meandiff
## 1  0.004703397 0.004953677 0.004441480  0.010101048  0.0045015049
## 2  0.117076477 0.004953677 0.004441480  0.006543836  0.0017272009
## 3  0.017269376 0.004953677 0.004441480  0.002678706  0.0001939040
## 4  0.004703397 0.004953677 0.004441480 -0.021390093 -0.0017306087
## 5  0.028163912 0.004953677 0.004441480 -0.021740837 -0.0061419923
## 6  0.083276948 0.004953677 0.004441480  0.002641622  0.0012332835
## 7  0.004703397 0.004953677 0.004441480 -0.004999324 -0.0008278370
## 8  0.235928140 0.004953677 0.021889323  0.005625006  0.0004115625
## 9  0.008014137 0.004953677 0.004441480 -0.015619994 -0.0038180681
## 10 0.092247660 0.004953677 0.004441480  0.015513750  0.0076934189
## 11 0.202359561 0.004953677 0.014249428 -0.008242199 -0.0013574566
## 12 0.004703397 0.007475687 0.004441480 -0.010388034 -0.0003756070
## 13 0.004703397 0.004953677 0.004441480 -0.005786158 -0.0018945286
## 14 0.565453433 0.004953677 0.038088204 -0.022620670 -0.0031239417
## 15 0.004703397 0.004953677 0.004441480 -0.024473249 -0.0160562177
## 16 0.092693699 0.004953677 0.017458525 -0.020586862 -0.0033703102
## 17 0.004703397 0.004953677 0.004441480 -0.021003147 -0.0096557712
## 18 0.133658256 0.004953677 0.007456698 -0.007912741  0.0002020913
## 19 0.130078847 0.004953677 0.005962454 -0.005628352 -0.0006570470
## 20 0.076892078 0.004953677 0.008415333 -0.009714773 -0.0010701446
##       overlapping.genes
## 1                  <NA>
## 2      RASSF4, C10orf10
## 3                FAM83H
## 4      RPS6KA4, MIR1237
## 5                MCOLN2
## 6          RP11-611O2.2
## 7                 CPLX1
## 8                 BEST1
## 9                 RBM47
## 10                LSAMP
## 11           AC104794.4
## 12              MSANTD4
## 13                 TNXB
## 14        SCAMP4, ADAT3
## 15                 DDR1
## 16                 <NA>
## 17                  EGF
## 18         ACBD4, PLCD3
## 19                AAMDC
## 20 KCNMB2, RP11-385J1.2
write.csv(results.ranges, file = "dmrcoutput_buccal.csv", row.names = TRUE)

Session information

sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 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       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] DMRcatedata_2.22.0                                 
##  [2] IlluminaHumanMethylation450kmanifest_0.4.0         
##  [3] IlluminaHumanMethylationEPICmanifest_0.3.0         
##  [4] RhpcBLASctl_0.23-42                                
##  [5] vioplot_0.5.0                                      
##  [6] zoo_1.8-12                                         
##  [7] sm_2.2-6.0                                         
##  [8] kableExtra_1.4.0                                   
##  [9] mitch_1.16.0                                       
## [10] FlowSorted.Blood.EPIC_2.8.0                        
## [11] ExperimentHub_2.12.0                               
## [12] AnnotationHub_3.12.0                               
## [13] BiocFileCache_2.12.0                               
## [14] dbplyr_2.5.0                                       
## [15] DMRcate_3.0.8                                      
## [16] ggplot2_3.5.1                                      
## [17] reshape2_1.4.4                                     
## [18] FlowSorted.Blood.450k_1.42.0                       
## [19] WGCNA_1.72-5                                       
## [20] fastcluster_1.2.6                                  
## [21] dynamicTreeCut_1.63-1                              
## [22] gplots_3.1.3.1                                     
## [23] RColorBrewer_1.1-3                                 
## [24] limma_3.60.4                                       
## [25] missMethyl_1.38.0                                  
## [26] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [27] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1 
## [28] minfi_1.50.0                                       
## [29] bumphunter_1.46.0                                  
## [30] locfit_1.5-9.10                                    
## [31] iterators_1.0.14                                   
## [32] foreach_1.5.2                                      
## [33] Biostrings_2.72.1                                  
## [34] XVector_0.44.0                                     
## [35] SummarizedExperiment_1.34.0                        
## [36] Biobase_2.64.0                                     
## [37] MatrixGenerics_1.16.0                              
## [38] matrixStats_1.3.0                                  
## [39] GenomicRanges_1.56.1                               
## [40] GenomeInfoDb_1.40.1                                
## [41] IRanges_2.38.1                                     
## [42] S4Vectors_0.42.1                                   
## [43] BiocGenerics_0.50.0                                
## 
## loaded via a namespace (and not attached):
##   [1] ProtGenerics_1.36.0       bitops_1.0-8             
##   [3] httr_1.4.7                doParallel_1.0.17        
##   [5] tools_4.4.1               doRNG_1.8.6              
##   [7] backports_1.5.0           utf8_1.2.4               
##   [9] R6_2.5.1                  HDF5Array_1.32.1         
##  [11] lazyeval_0.2.2            Gviz_1.48.0              
##  [13] rhdf5filters_1.16.0       permute_0.9-7            
##  [15] withr_3.0.0               GGally_2.2.1             
##  [17] prettyunits_1.2.0         gridExtra_2.3            
##  [19] base64_2.0.1              preprocessCore_1.61.0    
##  [21] cli_3.6.3                 labeling_0.4.3           
##  [23] sass_0.4.9                readr_2.1.5              
##  [25] genefilter_1.86.0         askpass_1.2.0            
##  [27] systemfonts_1.1.0         Rsamtools_2.20.0         
##  [29] foreign_0.8-86            svglite_2.1.3            
##  [31] siggenes_1.78.0           illuminaio_0.46.0        
##  [33] R.utils_2.12.3            dichromat_2.0-0.1        
##  [35] scrime_1.3.5              BSgenome_1.72.0          
##  [37] readxl_1.4.3              rstudioapi_0.16.0        
##  [39] impute_1.78.0             RSQLite_2.3.7            
##  [41] generics_0.1.3            BiocIO_1.14.0            
##  [43] gtools_3.9.5              dplyr_1.1.4              
##  [45] GO.db_3.19.1              Matrix_1.7-0             
##  [47] interp_1.1-6              fansi_1.0.6              
##  [49] abind_1.4-5               R.methodsS3_1.8.2        
##  [51] lifecycle_1.0.4           edgeR_4.2.1              
##  [53] yaml_2.3.8                rhdf5_2.48.0             
##  [55] SparseArray_1.4.8         grid_4.4.1               
##  [57] blob_1.2.4                promises_1.3.0           
##  [59] crayon_1.5.3              lattice_0.22-6           
##  [61] echarts4r_0.4.5           GenomicFeatures_1.56.0   
##  [63] annotate_1.82.0           KEGGREST_1.44.1          
##  [65] pillar_1.9.0              knitr_1.47               
##  [67] beanplot_1.3.1            rjson_0.2.22             
##  [69] codetools_0.2-20          glue_1.7.0               
##  [71] data.table_1.15.4         vctrs_0.6.5              
##  [73] png_0.1-8                 cellranger_1.1.0         
##  [75] gtable_0.3.5              cachem_1.1.0             
##  [77] xfun_0.45                 mime_0.12                
##  [79] S4Arrays_1.4.1            survival_3.7-0           
##  [81] statmod_1.5.0             nlme_3.1-165             
##  [83] bit64_4.0.5               bsseq_1.40.0             
##  [85] progress_1.2.3            filelock_1.0.3           
##  [87] bslib_0.7.0               nor1mix_1.3-3            
##  [89] KernSmooth_2.23-24        rpart_4.1.23             
##  [91] colorspace_2.1-1          DBI_1.2.3                
##  [93] Hmisc_5.1-3               nnet_7.3-19              
##  [95] tidyselect_1.2.1          bit_4.0.5                
##  [97] compiler_4.4.1            curl_5.2.1               
##  [99] httr2_1.0.1               htmlTable_2.4.3          
## [101] BiasedUrn_2.0.12          xml2_1.3.6               
## [103] DelayedArray_0.30.1       rtracklayer_1.64.0       
## [105] checkmate_2.3.2           scales_1.3.0             
## [107] caTools_1.18.2            quadprog_1.5-8           
## [109] rappdirs_0.3.3            stringr_1.5.1            
## [111] digest_0.6.36             rmarkdown_2.27           
## [113] GEOquery_2.72.0           htmltools_0.5.8.1        
## [115] pkgconfig_2.0.3           jpeg_0.1-10              
## [117] base64enc_0.1-3           sparseMatrixStats_1.16.0 
## [119] highr_0.11                fastmap_1.2.0            
## [121] ensembldb_2.28.1          rlang_1.1.4              
## [123] htmlwidgets_1.6.4         UCSC.utils_1.0.0         
## [125] shiny_1.8.1.1             DelayedMatrixStats_1.26.0
## [127] farver_2.1.2              jquerylib_0.1.4          
## [129] jsonlite_1.8.8            BiocParallel_1.38.0      
## [131] mclust_6.1.1              R.oo_1.26.0              
## [133] VariantAnnotation_1.50.0  RCurl_1.98-1.16          
## [135] magrittr_2.0.3            Formula_1.2-5            
## [137] GenomeInfoDbData_1.2.12   Rhdf5lib_1.26.0          
## [139] munsell_0.5.1             Rcpp_1.0.12              
## [141] stringi_1.8.4             zlibbioc_1.50.0          
## [143] MASS_7.3-61               plyr_1.8.9               
## [145] org.Hs.eg.db_3.19.1       ggstats_0.6.0            
## [147] deldir_2.0-4              splines_4.4.1            
## [149] multtest_2.60.0           hms_1.1.3                
## [151] rngtools_1.5.2            biomaRt_2.60.1           
## [153] BiocVersion_3.19.1        XML_3.99-0.17            
## [155] evaluate_0.24.0           latticeExtra_0.6-30      
## [157] biovizBase_1.52.0         BiocManager_1.30.23      
## [159] httpuv_1.6.15             tzdb_0.4.0               
## [161] tidyr_1.3.1               openssl_2.2.0            
## [163] purrr_1.0.2               reshape_0.8.9            
## [165] xtable_1.8-4              restfulr_0.0.15          
## [167] AnnotationFilter_1.28.0   later_1.3.2              
## [169] viridisLite_0.4.2         tibble_3.2.1             
## [171] beeswarm_0.4.0            memoise_2.0.1            
## [173] AnnotationDbi_1.66.0      GenomicAlignments_1.40.0 
## [175] cluster_2.1.6