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.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")
  library("vioplot")
})

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 = "ILMLEPIC-16276_SampleSheet.csv")
## [read.metharray.sheet] Found the following CSV files:
## [1] "/home/mdz/projects/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)
## snapshotDate(): 2023-04-24
## snapshotDate(): 2023-04-24
## 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 = "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(): 2023-04-24
## 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   6515597   6515624    28      *       2     2.821586e-04 0.967935
## 2      chr1  11174802  11174850    49      *       2     1.897069e-06 0.967935
## 3      chr1  16010724  16010732     9      *       2     6.987128e-04 0.967935
## 4      chr1 150739309 150739533   225      *       2     3.928301e-04 0.967935
## 5      chr1 152087267 152087315    49      *       2     3.535447e-04 0.967935
## 6      chr1 167035992 167036150   159      *       2     5.829688e-06 0.967935
## 7      chr1 177669385 177669462    78      *       2     6.729763e-04 0.967935
## 8      chr1 183560913 183560923    11      *       2     2.461174e-04 0.967935
## 9      chr1 218695426 218695459    34      *       2     4.032992e-04 0.967935
## 10     chr1 222227333 222227490   158      *       2     5.134055e-04 0.967935
## 11     chr2   7909341   7909485   145      *       2     8.740165e-05 0.967935
## 12     chr2  17941268  17941593   326      *       2     1.629182e-07 0.967935
## 13     chr2  27342798  27342877    80      *       2     4.941410e-04 0.967935
## 14     chr2  45879684  45880618   935      *       2     2.380826e-04 0.967935
## 15     chr2  73616832  73617110   279      *       2     2.416732e-04 0.967935
## 16     chr2 100106147 100106285   139      *       2     2.346173e-05 0.967935
## 17     chr2 101868356 101868516   161      *       2     8.728933e-04 0.967935
## 18     chr2 136742286 136742428   143      *       2     1.785300e-05 0.967935
## 19     chr2 153363436 153364045   610      *       2     1.496503e-06 0.967935
## 20     chr2 182548750 182548786    37      *       2     5.880878e-04 0.967935
##        HMFDR   Fisher      maxdiff      meandiff overlapping.genes
## 1  0.9047411 0.982442 -0.005000101 -0.0046987065              ESPN
## 2  0.9047411 0.982442 -0.011220365 -0.0089553767              MTOR
## 3  0.9047411 0.982442  0.001854973  0.0017242523              <NA>
## 4  0.9047411 0.982442 -0.016094921 -0.0106117989              <NA>
## 5  0.9047411 0.982442 -0.010921456 -0.0105541593              <NA>
## 6  0.9047411 0.982442  0.003724200  0.0036397926             GPA33
## 7  0.9047411 0.982442 -0.013327799 -0.0038219775              <NA>
## 8  0.9047411 0.982442 -0.009560729 -0.0003340177              SMG7
## 9  0.9047411 0.982442 -0.012709739 -0.0093789342          C1orf143
## 10 0.9047411 0.982442 -0.017413191 -0.0119660084     RP11-400N13.3
## 11 0.9047411 0.982442 -0.020877842 -0.0071641853              <NA>
## 12 0.9047411 0.982442 -0.017672562 -0.0124820782        GEN1, SMC6
## 13 0.9047411 0.982442 -0.006619968 -0.0065196253              <NA>
## 14 0.9047411 0.982442 -0.012151324 -0.0049426212             PRKCE
## 15 0.9047411 0.982442 -0.013602305 -0.0089560660             ALMS1
## 16 0.9047411 0.982442  0.001872926  0.0004131779              REV1
## 17 0.9047411 0.982442  0.011164324  0.0091928148            TBC1D8
## 18 0.9047411 0.982442  0.004480379  0.0031129951              DARS
## 19 0.9047411 0.982442 -0.018027597 -0.0132216403             FMNL2
## 20 0.9047411 0.982442  0.004084872  0.0027096631        AC013733.3
write.csv(results.ranges, file = "dmrcoutput_buccal.csv", row.names = TRUE)

Session information

sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 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       
## 
## time zone: Australia/Melbourne
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] DMRcatedata_2.18.0                                 
##  [2] IlluminaHumanMethylationEPICmanifest_0.3.0         
##  [3] vioplot_0.4.0                                      
##  [4] zoo_1.8-12                                         
##  [5] sm_2.2-5.7.1                                       
##  [6] kableExtra_1.3.4                                   
##  [7] mitch_1.12.0                                       
##  [8] FlowSorted.Blood.EPIC_2.4.2                        
##  [9] ExperimentHub_2.8.1                                
## [10] AnnotationHub_3.8.0                                
## [11] BiocFileCache_2.8.0                                
## [12] dbplyr_2.3.3                                       
## [13] DMRcate_2.14.1                                     
## [14] ggplot2_3.4.3                                      
## [15] reshape2_1.4.4                                     
## [16] FlowSorted.Blood.450k_1.38.0                       
## [17] WGCNA_1.72-1                                       
## [18] fastcluster_1.2.3                                  
## [19] dynamicTreeCut_1.63-1                              
## [20] gplots_3.1.3                                       
## [21] RColorBrewer_1.1-3                                 
## [22] ruv_0.9.7.1                                        
## [23] IlluminaHumanMethylationEPICanno.ilm10b2.hg19_0.6.0
## [24] limma_3.56.2                                       
## [25] missMethyl_1.34.0                                  
## [26] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [27] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1 
## [28] minfi_1.46.0                                       
## [29] bumphunter_1.42.0                                  
## [30] locfit_1.5-9.8                                     
## [31] iterators_1.0.14                                   
## [32] foreach_1.5.2                                      
## [33] Biostrings_2.68.1                                  
## [34] XVector_0.40.0                                     
## [35] SummarizedExperiment_1.30.2                        
## [36] Biobase_2.60.0                                     
## [37] MatrixGenerics_1.12.3                              
## [38] matrixStats_1.0.0                                  
## [39] GenomicRanges_1.52.0                               
## [40] GenomeInfoDb_1.36.3                                
## [41] IRanges_2.34.1                                     
## [42] S4Vectors_0.38.1                                   
## [43] BiocGenerics_0.46.0                                
## 
## loaded via a namespace (and not attached):
##   [1] DSS_2.48.0                    ProtGenerics_1.32.0          
##   [3] bitops_1.0-7                  webshot_0.5.5                
##   [5] httr_1.4.7                    doParallel_1.0.17            
##   [7] tools_4.3.1                   doRNG_1.8.6                  
##   [9] backports_1.4.1               utf8_1.2.3                   
##  [11] R6_2.5.1                      HDF5Array_1.28.1             
##  [13] lazyeval_0.2.2                Gviz_1.44.1                  
##  [15] rhdf5filters_1.12.1           permute_0.9-7                
##  [17] withr_2.5.0                   GGally_2.1.2                 
##  [19] prettyunits_1.1.1             gridExtra_2.3                
##  [21] base64_2.0.1                  preprocessCore_1.62.1        
##  [23] cli_3.6.1                     sass_0.4.7                   
##  [25] readr_2.1.4                   genefilter_1.82.1            
##  [27] askpass_1.2.0                 systemfonts_1.0.4            
##  [29] Rsamtools_2.16.0              foreign_0.8-85               
##  [31] svglite_2.1.1                 siggenes_1.74.0              
##  [33] illuminaio_0.42.0             R.utils_2.12.2               
##  [35] dichromat_2.0-0.1             scrime_1.3.5                 
##  [37] BSgenome_1.68.0               readxl_1.4.3                 
##  [39] rstudioapi_0.15.0             impute_1.74.1                
##  [41] RSQLite_2.3.1                 generics_0.1.3               
##  [43] BiocIO_1.10.0                 gtools_3.9.4                 
##  [45] dplyr_1.1.3                   GO.db_3.17.0                 
##  [47] Matrix_1.6-1                  interp_1.1-4                 
##  [49] fansi_1.0.4                   abind_1.4-5                  
##  [51] R.methodsS3_1.8.2             lifecycle_1.0.3              
##  [53] edgeR_3.42.4                  yaml_2.3.7                   
##  [55] rhdf5_2.44.0                  grid_4.3.1                   
##  [57] blob_1.2.4                    promises_1.2.1               
##  [59] crayon_1.5.2                  lattice_0.21-8               
##  [61] echarts4r_0.4.5               GenomicFeatures_1.52.2       
##  [63] annotate_1.78.0               KEGGREST_1.40.0              
##  [65] pillar_1.9.0                  knitr_1.44                   
##  [67] beanplot_1.3.1                rjson_0.2.21                 
##  [69] codetools_0.2-19              glue_1.6.2                   
##  [71] data.table_1.14.8             vctrs_0.6.3                  
##  [73] png_0.1-8                     cellranger_1.1.0             
##  [75] gtable_0.3.4                  cachem_1.0.8                 
##  [77] xfun_0.40                     mime_0.12                    
##  [79] S4Arrays_1.0.6                survival_3.5-7               
##  [81] statmod_1.5.0                 ellipsis_0.3.2               
##  [83] interactiveDisplayBase_1.38.0 nlme_3.1-163                 
##  [85] bit64_4.0.5                   bsseq_1.36.0                 
##  [87] progress_1.2.2                filelock_1.0.2               
##  [89] bslib_0.5.1                   nor1mix_1.3-0                
##  [91] KernSmooth_2.23-22            rpart_4.1.19                 
##  [93] colorspace_2.1-0              DBI_1.1.3                    
##  [95] Hmisc_5.1-1                   nnet_7.3-19                  
##  [97] tidyselect_1.2.0              bit_4.0.5                    
##  [99] compiler_4.3.1                curl_5.0.2                   
## [101] rvest_1.0.3                   htmlTable_2.4.1              
## [103] xml2_1.3.5                    DelayedArray_0.26.7          
## [105] rtracklayer_1.60.1            checkmate_2.2.0              
## [107] scales_1.2.1                  caTools_1.18.2               
## [109] quadprog_1.5-8                rappdirs_0.3.3               
## [111] stringr_1.5.0                 digest_0.6.33                
## [113] rmarkdown_2.24                GEOquery_2.68.0              
## [115] htmltools_0.5.6               pkgconfig_2.0.3              
## [117] jpeg_0.1-10                   base64enc_0.1-3              
## [119] sparseMatrixStats_1.12.2      highr_0.10                   
## [121] fastmap_1.1.1                 ensembldb_2.24.0             
## [123] rlang_1.1.1                   htmlwidgets_1.6.2            
## [125] shiny_1.7.5                   DelayedMatrixStats_1.22.6    
## [127] jquerylib_0.1.4               jsonlite_1.8.7               
## [129] BiocParallel_1.34.2           mclust_6.0.0                 
## [131] R.oo_1.25.0                   VariantAnnotation_1.46.0     
## [133] RCurl_1.98-1.12               magrittr_2.0.3               
## [135] Formula_1.2-5                 GenomeInfoDbData_1.2.10      
## [137] Rhdf5lib_1.22.1               munsell_0.5.0                
## [139] Rcpp_1.0.11                   stringi_1.7.12               
## [141] zlibbioc_1.46.0               MASS_7.3-60                  
## [143] plyr_1.8.8                    org.Hs.eg.db_3.17.0          
## [145] deldir_1.0-9                  splines_4.3.1                
## [147] multtest_2.56.0               hms_1.1.3                    
## [149] rngtools_1.5.2                biomaRt_2.56.1               
## [151] BiocVersion_3.17.1            XML_3.99-0.14                
## [153] evaluate_0.21                 latticeExtra_0.6-30          
## [155] biovizBase_1.48.0             BiocManager_1.30.22          
## [157] httpuv_1.6.11                 tzdb_0.4.0                   
## [159] tidyr_1.3.0                   openssl_2.1.0                
## [161] purrr_1.0.2                   reshape_0.8.9                
## [163] xtable_1.8-4                  restfulr_0.0.15              
## [165] AnnotationFilter_1.24.0       later_1.3.1                  
## [167] viridisLite_0.4.2             tibble_3.2.1                 
## [169] beeswarm_0.4.0                memoise_2.0.1                
## [171] AnnotationDbi_1.62.2          GenomicAlignments_1.36.0     
## [173] cluster_2.1.4