Autism Spectrum Disorder Methylation analysis - Guthrie card analysis
Namitha Mohandras, Mark Ziemann
2024-10-08
Source: https://github.com/markziemann/asd_meth
Introduction
Here we are analysing Guthrie card methylation in 9 twin pairs.
baseDir=getwd()
dataDir=paste(baseDir,"/ASD_EPIC_DATA",sep="")
suppressPackageStartupMessages({
library("missMethyl")
library("limma")
library("minfi")
library("IlluminaHumanMethylation450kanno.ilmn12.hg19")
library("IlluminaHumanMethylationEPICanno.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 = "ASD_guthrie_summarysheet.csv")## [read.metharray.sheet] Found the following CSV files:## [1] "/asd_meth/ASD_EPIC_DATA/ASD_guthrie_summarysheet.csv"#targets$ID = paste(targets$Sample_Group,targets_gen$Sample_Name,sep=".")
rgSet = read.metharray.exp(targets = targets_gen)
sampleNames(rgSet) = targets_gen$SampleIDTesting MZ status
snpBetas = getSnpBeta(rgSet)
d = dist(t(snpBetas))
hr = hclust(d, method = "complete", members=NULL)
plot(hr)
Quality control
detP = detectionP(rgSet)
qcReport(rgSet, sampNames = targets_gen$Sample_Name,
pdf="qc-report_ASD_guthrie_nov27.pdf")## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE
## Warning in type.convert.default(X[[i]], ...): 'as.is' should be specified by
## the caller; using TRUE## png
## 2cols=brewer.pal(4,"Set1")
barplot(apply(detP,2,mean),
col=as.numeric(factor(targets_gen$Sample_Name)),
las=2,cex.names= 0.8, cex.axis=0.75,
main="Mean detection p-values of probe signals",
ylab="Mean detection p-value")
barplot(apply(detP,2,mean),
col=as.numeric(factor(targets_gen$Sample_Name)),
las=2,cex.names= 0.8, cex.axis=0.75,ylim=c(0,0.010),
main="Mean detection p-values of probe signals",
ylab="Mean detection p-value")
Cell type composition analysis
Using old method.
rgSet$Slide <- as.numeric(rgSet$Slide)
rgSet$Sex <- as.character(rgSet$Sex)
#rgSet$Sample_Name <- as.character(rgSet$Sample_Name)
#sampleNames(rgSet) = targets_gen$SampleID
cellCounts_new <- estimateCellCounts(rgSet, compositeCellType = "Blood",
processMethod = "auto", probeSelect = "auto",
cellTypes = c("CD8T","CD4T", "NK","Bcell","Mono","Gran"),
referencePlatform = c("IlluminaHumanMethylation450k"),
returnAll = FALSE, meanPlot = TRUE)## [estimateCellCounts] Combining user data with reference (flow sorted) data.## Warning in DataFrame(sampleNames = c(colnames(rgSet),
## colnames(referenceRGset)), : 'stringsAsFactors' is ignored## [estimateCellCounts] Processing user and reference data together.## [preprocessQuantile] Mapping to genome.## [preprocessQuantile] Fixing outliers.## [preprocessQuantile] Quantile normalizing.## [estimateCellCounts] Picking probes for composition estimation.## [estimateCellCounts] Estimating composition.
#plot cell type composition by sample group
a = cellCounts_new[targets_gen$Diagnosis == "0",]
b = cellCounts_new[targets_gen$Diagnosis == "1",]
c = cellCounts_new[targets_gen$Diagnosis == "2",]
age.pal <- brewer.pal(8,"Set1")
cellCounts_long <- melt(cellCounts_new, id = "celltype")
cellCounts_long$Var1 <- targets_gen[match(cellCounts_long$Var1, targets_gen$SampleID),"Diagnosis"]
colnames(cellCounts_long) <- c("diagnosis","celltype","value")
ggplot(cellCounts_long, aes(x = factor(celltype), y = value, fill = factor(diagnosis))) +
geom_boxplot() 
wx <- lapply(1:ncol(cellCounts_new) , function(i) {
wilcox.test(a[,i],c[,i], paired=FALSE)
} )## Warning in wilcox.test.default(a[, i], c[, i], paired = FALSE): cannot compute
## exact p-value with tiesnames(wx) <- colnames(cellCounts_new)
wx## $CD8T
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 17, p-value = 0.2674
## alternative hypothesis: true location shift is not equal to 0
##
##
## $CD4T
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 37, p-value = 0.3196
## alternative hypothesis: true location shift is not equal to 0
##
##
## $NK
##
## Wilcoxon rank sum test with continuity correction
##
## data: a[, i] and c[, i]
## W = 35.5, p-value = 0.3289
## alternative hypothesis: true location shift is not equal to 0
##
##
## $Bcell
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 36, p-value = 0.3773
## alternative hypothesis: true location shift is not equal to 0
##
##
## $Mono
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 23, p-value = 0.6612
## alternative hypothesis: true location shift is not equal to 0
##
##
## $Gran
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 26, p-value = 0.913
## alternative hypothesis: true location shift is not equal to 0Try estimatecellcounts2.
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.mset <- cells$normalizedData
cellCounts_new <- cells[[1]]
#plot cell type composition by sample group
a = cellCounts_new[targets_gen$Diagnosis == "0",]
b = cellCounts_new[targets_gen$Diagnosis == "1",]
c = cellCounts_new[targets_gen$Diagnosis == "2",]
age.pal <- brewer.pal(8,"Set1")
cellCounts_long <- melt(cellCounts_new, id = "celltype")
cellCounts_long$Var1 <- targets_gen[match(cellCounts_long$Var1, targets_gen$SampleID),"Diagnosis"]
colnames(cellCounts_long) <- c("diagnosis","celltype","value")
ggplot(cellCounts_long, aes(x = factor(celltype), y = value, fill = factor(diagnosis))) +
geom_boxplot()
pdf("cells_guthrie.pdf")
par(mar=c(3,8,2,2))
ggplot(cellCounts_long, aes(x = factor(celltype), y = value, fill = factor(diagnosis))) +
geom_boxplot()
dev.off()## png
## 2wx <- lapply(1:ncol(cellCounts_new) , function(i) {
wilcox.test(a[,i],c[,i], paired=FALSE)
} )
names(wx) <- colnames(cellCounts_new)
wx## $CD8T
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 34, p-value = 0.5096
## alternative hypothesis: true location shift is not equal to 0
##
##
## $CD4T
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 36, p-value = 0.3773
## alternative hypothesis: true location shift is not equal to 0
##
##
## $NK
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 19, p-value = 0.3773
## alternative hypothesis: true location shift is not equal to 0
##
##
## $Bcell
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 40, p-value = 0.1804
## alternative hypothesis: true location shift is not equal to 0
##
##
## $Mono
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 20, p-value = 0.4409
## alternative hypothesis: true location shift is not equal to 0
##
##
## $Neu
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 29, p-value = 0.913
## alternative hypothesis: true location shift is not equal to 0Filtering
None of the p-values are less than 0.05.
Now filter for high detection p-value and overlap with SNP.
Need to get a copy of the Xreact probes from Namitha. In the mean time I have downloaded from github - link below.
detP <- detectionP(rgSet)
keep <- rowSums(detP < 0.01) == ncol(rgSet)
mset <- mset[keep,]
gmset <- mapToGenome(mset)
#remove SNPs
gmset_flt = dropLociWithSnps(gmset, snps = c("CpG", "SBE"))
#Removing cross-reactive probes
XURL="https://raw.githubusercontent.com/sirselim/illumina450k_filtering/master/EPIC/13059_2016_1066_MOESM1_ESM.csv"
Xreact <- read.csv(XURL)
#Xreact = read.csv(file="/group/canc2/puumba/Data/InfiniumData/NamithaData/Rprojects/Autism/Analysis_Sept11/EPIC_850k_crossreactiveProbes.csv", stringsAsFactors=FALSE)
#Xreact = read.csv(file="~/48639-non-specific-probes-Illumina450k.csv", stringsAsFactors=FALSE)
noXreact <- !(featureNames(gmset) %in% Xreact$X)
gmset <- gmset[noXreact,]
#Removing probes on X and Y chromosomes
autosomes <- !(featureNames(gmset) %in% ann$Name[ann$chr %in% c("chrX","chrY")])
gmset_flt <- gmset[autosomes,]
#getBeta
beta = getM(gmset_flt)
saveRDS(beta,"gu_beta.rds")
df <- data.frame(t(t(colMeans(beta))))
colnames(df) = "gwam_gu"
write.table(df,file="gu_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)
fit <- prcomp(t(mvals),center = TRUE, scale = TRUE,retx=TRUE)
loadings = fit$x
plot(fit,type="lines",col="blue")
saveRDS(fit,"scree_guthrie_fircorrect.Rds")
#fitnc <- readRDS("scree_guthrie_fitnc.Rds")
plot(fit,type="lines",col="blue")
nGenes = nrow(mvals)
nSamples = ncol(mvals)
datTraits = targets_gen[,7:15]
moduleTraitCor = cor(loadings[,1:6], datTraits, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
par(cex=0.75, mar = c(6, 8.5, 3, 3))
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
labeledHeatmap(Matrix = t(moduleTraitCor),
xLabels = colnames(loadings)[1:ncol(t(moduleTraitCor))],
yLabels = names(datTraits), colorLabels = FALSE, colors = blueWhiteRed(6),
textMatrix = t(textMatrix), setStdMargins = FALSE, cex.text = 0.5,
cex.lab.y = 0.6, zlim = c(-1,1),
main = paste("PCA-trait relationships: Top principal components"))
pdf("pca_guthrie.pdf")
par(mar=c(3,8,2,2))
labeledHeatmap(Matrix = t(moduleTraitCor),
xLabels = colnames(loadings)[1:ncol(t(moduleTraitCor))],
yLabels = names(datTraits), colorLabels = FALSE, colors = blueWhiteRed(6),
textMatrix = t(textMatrix), setStdMargins = FALSE, cex.text = 0.5,
cex.lab.y = 0.6, zlim = c(-1,1),
main = paste("PCA-trait relationships: Top principal components"))
dev.off()## png
## 2Regression analysis - within twin pair on ADOS
ADOS is the outcome variable (continuous)
Twin_Group <- factor(targets_gen$Family_ID)
ADOS <- targets_gen$Social_interaction_on_ADOS
names(ADOS) <- targets_gen$SampleID
motor <- targets_gen$Motor_skills
diagnosis <- targets_gen$Diagnosis
SRS<- targets_gen$SRS_social_scores
iiq <- 1/targets_gen$IQ
# fix infinite
iiq[iiq==Inf] <- 0.025
iiq <- as.vector(scale(iiq))
ilanguage <- as.vector(scale(1/targets_gen$language))
#ADOS:0,4.009335e-06
design_tw <- model.matrix(~ Twin_Group + ADOS )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))## (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down 246029 2316 2580 2465 1947 3518
## NotSig 59324 786067 786324 786237 786853 785072
## Up 485305 2275 1754 1956 1858 2068
## Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13 ADOS
## Down 48180 2078 2902 1867 2607 0
## NotSig 728758 786577 785834 786912 786340 790658
## Up 13720 2003 1922 1879 1711 0top <- topTable(fit2_tw,coef="ADOS",num=Inf, sort.by = "P")
nrow(subset(top,adj.P.Val < 0.05))## [1] 0nrow(subset(top,P.Value < 1e-4))## [1] 14output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_guthrie_ADOS.csv",row.names=FALSE)
saveRDS(design_tw, "gu_design_ados.rds")
saveRDS(mvals, "gu_mvals.rds")
head(output,30) %>% kbl(caption="ADOS") %>% kable_paper("hover", full_width = F)| Name | chr | pos | strand | Islands_Name | Relation_to_Island | UCSC_RefGene_Name | UCSC_RefGene_Group | logFC | AveExpr | t | P.Value | adj.P.Val | B |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| cg00546248 | chr1 | 2412743 |
|
chr1:2411125-2411404 | S_Shore | PLCH2 | Body | 0.1030486 | 1.7980001 | 6.898679 | 0.0000090 | 0.9999949 | 3.8484926 |
| cg15740041 | chr16 | 87781336 |
|
OpenSea | KLHDC4;KLHDC4;KLHDC4 | Body;Body;Body | 0.1314927 | 1.5969029 | 6.830272 | 0.0000100 | 0.9999949 | 3.7505238 | |
| cg05116656 | chr5 | 94417892 |
|
OpenSea | MCTP1;MCTP1;MCTP1 | TSS1500;TSS1500;Body | 0.0915052 | -3.3779970 | 6.785304 | 0.0000108 | 0.9999949 | 3.6857228 | |
| cg11632273 | chr22 | 20102580 |
|
chr22:20104374-20105729 | N_Shore | TRMT2A;TRMT2A;TRMT2A;RANBP1;MIR6816 | Body;Body;Body;TSS1500;TSS1500 | -0.1327794 | 2.2652795 | -6.340401 | 0.0000219 | 0.9999949 | 3.0273497 |
| cg01656620 | chr12 | 124788698 |
|
chr12:124788559-124788884 | Island | FAM101A | 5’UTR | -0.1328796 | 3.0466588 | -6.207538 | 0.0000272 | 0.9999949 | 2.8246260 |
| cg10881749 | chr1 | 234792390 |
|
OpenSea | 0.0804076 | 0.9538435 | 6.180843 | 0.0000284 | 0.9999949 | 2.7835535 | |||
| cg13699963 | chr17 | 71548861 |
|
OpenSea | SDK2 | Body | -0.0664415 | 2.9744674 | -5.859973 | 0.0000486 | 0.9999949 | 2.2809795 | |
| cg07162704 | chr6 | 91113687 |
|
OpenSea | 0.1007837 | 0.5651694 | 5.696819 | 0.0000641 | 0.9999949 | 2.0191661 | |||
| cg13277464 | chr2 | 43393040 |
|
OpenSea | 0.1210637 | 2.2975670 | 5.616519 | 0.0000736 | 0.9999949 | 1.8887704 | |||
| cg05791256 | chr9 | 96077420 |
|
chr9:96080113-96080330 | N_Shelf | WNK2;WNK2 | Body;Body | -0.0750945 | 3.6558579 | -5.535122 | 0.0000847 | 0.9999949 | 1.7555692 |
| cg03188919 | chr10 | 31372663 |
|
OpenSea | -0.1027692 | 2.0636053 | -5.534005 | 0.0000849 | 0.9999949 | 1.7537327 | |||
| cg01128603 | chr11 | 65820391 |
|
chr11:65819367-65820224 | S_Shore | SF3B2 | Body | 0.0946669 | -3.2456959 | 5.519339 | 0.0000871 | 0.9999949 | 1.7296217 |
| cg23157290 | chr7 | 138602695 |
|
OpenSea | KIAA1549;KIAA1549 | Body;Body | 0.1523237 | 3.6464280 | 5.504452 | 0.0000893 | 0.9999949 | 1.7051115 | |
| cg03177876 | chr4 | 71701931 |
|
chr4:71704664-71705736 | N_Shelf | GRSF1;GRSF1 | 5’UTR;Body | 0.0982485 | -1.2925574 | 5.485795 | 0.0000923 | 0.9999949 | 1.6743466 |
| cg23061257 | chr1 | 54355272 |
|
chr1:54355354-54355656 | N_Shore | YIPF1 | 5’UTR | -0.0648132 | -3.3661385 | -5.434827 | 0.0001009 | 0.9999949 | 1.5900303 |
| cg25723331 | chr11 | 44613980 |
|
OpenSea | CD82;CD82 | 5’UTR;5’UTR | 0.1123393 | 0.7995294 | 5.404970 | 0.0001063 | 0.9999949 | 1.5404530 | |
| cg01016829 | chr16 | 46682692 |
|
OpenSea | 0.0860348 | 2.2645018 | 5.388257 | 0.0001095 | 0.9999949 | 1.5126430 | |||
| cg10619752 | chr2 | 233863457 |
|
OpenSea | NGEF | 5’UTR | -0.0821679 | 3.0471681 | -5.339104 | 0.0001194 | 0.9999949 | 1.4306045 | |
| cg08008144 | chr8 | 1983127 |
|
OpenSea | -0.0927043 | 2.5733842 | -5.335420 | 0.0001201 | 0.9999949 | 1.4244401 | |||
| cg25652742 | chr2 | 240152334 |
|
chr2:240153476-240153792 | N_Shore | HDAC4 | Body | 0.1178984 | -1.9304540 | 5.330430 | 0.0001212 | 0.9999949 | 1.4160888 |
| cg23404427 | chr19 | 1653402 |
|
OpenSea | TCF3 | TSS1500 | -0.0903214 | -5.0274142 | -5.322169 | 0.0001230 | 0.9999949 | 1.4022536 | |
| cg08225137 | chr3 | 54915807 |
|
OpenSea | CACNA2D3 | Body | -0.0559477 | 2.6416704 | -5.320145 | 0.0001234 | 0.9999949 | 1.3988625 | |
| cg15856662 | chr22 | 38380501 |
|
chr22:38379093-38379964 | S_Shore | SOX10;SOX10;POLR2F;POLR2F | 5’UTR;1stExon;Body;Body | -0.0731822 | 2.4710208 | -5.213569 | 0.0001491 | 0.9999949 | 1.2194325 |
| cg00183468 | chr3 | 142926178 |
|
OpenSea | 0.0609202 | 1.9081959 | 5.151478 | 0.0001665 | 0.9999949 | 1.1141210 | |||
| cg11355135 | chr2 | 108602860 |
|
chr2:108602824-108603467 | Island | SLC5A7 | TSS200 | 0.1105003 | -3.8213288 | 5.120093 | 0.0001762 | 0.9999949 | 1.0606732 |
| cg15165643 | chr8 | 145518600 |
|
chr8:145514719-145515959 | S_Shelf | HSF1 | Body | -0.0813514 | 3.0891873 | -5.119631 | 0.0001763 | 0.9999949 | 1.0598861 |
| cg17144383 | chr9 | 96868748 |
|
OpenSea | PTPDC1;PTPDC1;PTPDC1;PTPDC1 | Body;Body;Body;Body | 0.0776330 | 1.6439371 | 5.102222 | 0.0001819 | 0.9999949 | 1.0301762 | |
| cg07375358 | chr11 | 114056431 |
|
OpenSea | ZBTB16;ZBTB16 | Body;Body | 0.0903823 | 0.3404948 | 5.073839 | 0.0001914 | 0.9999949 | 0.9816468 | |
| cg17852290 | chr16 | 83762355 |
|
OpenSea | CDH13;CDH13;CDH13;CDH13 | Body;Body;Body;Body | 0.0854262 | 1.3389770 | 5.058992 | 0.0001966 | 0.9999949 | 0.9562140 | |
| cg20740133 | chr5 | 106792666 |
|
OpenSea | EFNA5 | Body | 0.0651929 | 1.9314424 | 5.042657 | 0.0002025 | 0.9999949 | 0.9281973 |
Other models:
motor skills impairment
diagnosis
SRS
inverse IQ
inverse langguage
Regression analysis - within twin pair on motor skills impairment
#motor:820,1.053055e-12
design_tw <- model.matrix(~ Twin_Group + motor )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))## (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down 222760 3101 3113 1626 1500 2520
## NotSig 109515 784537 785845 787951 788096 786863
## Up 458383 3020 1700 1081 1062 1275
## Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13 motor
## Down 153258 2952 3106 1384 1889 444
## NotSig 607527 784694 785576 788169 787716 789820
## Up 29873 3012 1976 1105 1053 394top <- 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] 54093output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_guthrie_motor.csv",row.names=FALSE)
head(output,30) %>% kbl(caption="motor") %>% kable_paper("hover", full_width = F)| Name | chr | pos | strand | Islands_Name | Relation_to_Island | UCSC_RefGene_Name | UCSC_RefGene_Group | logFC | AveExpr | t | P.Value | adj.P.Val | B |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| cg22735331 | chr11 | 2369390 |
|
OpenSea | -4.484327 | 2.7410066 | -24.71697 | 0 | 0.0000010 | 11.327823 | |||
| cg08840551 | chr5 | 56718349 |
|
OpenSea | -3.005235 | -0.2449499 | -23.15409 | 0 | 0.0000012 | 11.117561 | |||
| cg17121412 | chr15 | 22833400 |
|
chr15:22833282-22833911 | Island | TUBGCP5;TUBGCP5;TUBGCP5;TUBGCP5 | 1stExon;5’UTR;5’UTR;1stExon | -2.164138 | -5.4037488 | -18.92657 | 0 | 0.0000112 | 10.333772 |
| cg09002832 | chr9 | 4300181 |
|
chr9:4297817-4300182 | Island | GLIS3 | TSS200 | 2.574426 | -4.7920560 | 17.58274 | 0 | 0.0000181 | 9.992929 |
| cg06741568 | chr12 | 129028641 |
|
OpenSea | TMEM132C | Body | 2.952675 | 4.8924071 | 17.52976 | 0 | 0.0000181 | 9.978313 | |
| cg05098590 | chr10 | 111767379 |
|
chr10:111767087-111768355 | Island | ADD3;ADD3;ADD3 | 5’UTR;TSS1500;TSS1500 | 2.376842 | -5.5089331 | 16.70940 | 0 | 0.0000280 | 9.739318 |
| cg25259284 | chr1 | 52521812 |
|
chr1:52521780-52522296 | Island | TXNDC12;TXNDC12;BTF3L4;BTF3L4;BTF3L4;TXNDC12;TXNDC12 | 1stExon;5’UTR;TSS200;TSS200;TSS200;Body;Body | 2.467156 | -4.1947802 | 16.15712 | 0 | 0.0000309 | 9.563949 |
| cg09804439 | chr17 | 40540457 |
|
chr17:40539837-40540775 | Island | STAT3;STAT3;STAT3;STAT3;STAT3 | 5’UTR;1stExon;1stExon;5’UTR;TSS200 | 2.102328 | -5.6125956 | 16.11285 | 0 | 0.0000309 | 9.549348 |
| cg23199907 | chr13 | 33305966 |
|
OpenSea | PDS5B | Body | -3.694601 | 2.7345231 | -16.06816 | 0 | 0.0000309 | 9.534528 | |
| cg27533482 | chr3 | 28390637 |
|
chr3:28389975-28390855 | Island | AZI2;AZI2;AZI2 | TSS200;TSS200;TSS200 | -1.831110 | -5.5148828 | -15.51524 | 0 | 0.0000408 | 9.343893 |
| cg04055615 | chr6 | 31697723 |
|
chr6:31695894-31698245 | Island | DDAH2 | 5’UTR | 2.216171 | -4.9480845 | 15.40152 | 0 | 0.0000408 | 9.302958 |
| cg07241733 | chr8 | 68255213 |
|
chr8:68255129-68256274 | Island | ARFGEF1 | Body | 2.255660 | -3.1419951 | 15.37574 | 0 | 0.0000408 | 9.293594 |
| cg01593834 | chr5 | 145718401 |
|
chr5:145718289-145720095 | Island | POU4F3 | TSS200 | -2.706083 | -1.8199821 | -15.25612 | 0 | 0.0000417 | 9.249729 |
| cg04802271 | chr22 | 39930989 |
|
chr22:39930803-39931115 | Island | 1.649464 | -6.3475686 | 14.78585 | 0 | 0.0000577 | 9.070416 | ||
| cg15592945 | chr1 | 226736711 |
|
chr1:226736354-226737412 | Island | C1orf95 | 1stExon | 1.737831 | -5.6952786 | 14.22197 | 0 | 0.0000883 | 8.840089 |
| cg07220356 | chr6 | 33263442 |
|
chr6:33266302-33267582 | N_Shelf | RGL2;RGL2 | Body;Body | -1.842579 | 3.5862188 | -13.97553 | 0 | 0.0000941 | 8.733810 |
| cg14883700 | chr2 | 173096328 |
|
chr2:173099482-173100098 | N_Shelf | -1.963418 | 3.6911719 | -13.95601 | 0 | 0.0000941 | 8.725239 | ||
| cg23422659 | chr17 | 44929215 |
|
chr17:44928287-44929690 | Island | WNT9B | Body | -1.862823 | -5.2103189 | -13.94869 | 0 | 0.0000941 | 8.722020 |
| cg09669099 | chr16 | 69419557 |
|
chr16:69419316-69420086 | Island | TERF2 | 1stExon | 1.609111 | -5.1850149 | 13.78537 | 0 | 0.0001015 | 8.649344 |
| cg04114460 | chr18 | 29239388 |
|
OpenSea | B4GALT6 | Body | -2.310542 | 2.8339004 | -13.75061 | 0 | 0.0001015 | 8.633664 | |
| cg01957330 | chr7 | 26904470 |
|
chr7:26903869-26904753 | Island | SKAP2 | TSS200 | 1.999595 | -4.4330191 | 13.55569 | 0 | 0.0001124 | 8.544358 |
| cg04426862 | chr14 | 37131008 |
|
chr14:37131181-37132785 | N_Shore | PAX9 | 5’UTR | 1.653746 | -4.4862989 | 13.53799 | 0 | 0.0001124 | 8.536129 |
| cg26516701 | chr12 | 53845683 |
|
chr12:53845668-53846606 | Island | PCBP2;PCBP2;PCBP2;PCBP2;PCBP2;PCBP2;PCBP2 | TSS1500;TSS1500;TSS1500;TSS1500;TSS1500;TSS1500;TSS1500 | 1.502632 | -4.5177605 | 13.32026 | 0 | 0.0001239 | 8.433260 |
| cg10195929 | chr12 | 58019223 |
|
chr12:58021294-58022037 | N_Shelf | SLC26A10 | Body | -1.985173 | 2.7720489 | -13.29862 | 0 | 0.0001239 | 8.422867 |
| cg26763618 | chr16 | 14726469 |
|
chr16:14726516-14727366 | N_Shore | BFAR | TSS200 | 1.907346 | -4.4003705 | 13.23363 | 0 | 0.0001239 | 8.391467 |
| cg13523875 | chr7 | 117823728 |
|
chr7:117823984-117824485 | N_Shore | LSM8 | TSS1500 | 1.508943 | -5.5300468 | 13.23203 | 0 | 0.0001239 | 8.390689 |
| cg21434530 | chr4 | 4870486 |
|
chr4:4868440-4869173 | S_Shore | 1.615972 | -4.2911102 | 13.21643 | 0 | 0.0001239 | 8.383109 | ||
| cg02005873 | chr6 | 10385839 |
|
chr6:10383525-10384114 | S_Shore | 1.833494 | -1.2496069 | 12.70945 | 0 | 0.0001952 | 8.127591 | ||
| cg04570283 | chr14 | 82000256 |
|
chr14:81999751-82000434 | Island | SEL1L | TSS200 | 1.710792 | -2.9036410 | 12.66672 | 0 | 0.0001966 | 8.105219 |
| cg07762054 | chr11 | 75273187 |
|
chr11:75272950-75273616 | Island | SERPINH1;SERPINH1 | 5’UTR;1stExon | 2.139041 | -5.6365570 | 12.35060 | 0 | 0.0002591 | 7.935489 |
saveRDS(design_tw, "gu_design_motor.rds")Regression analysis - within twin pair on diagnosis
#diagnosis:0,3.080605e-05
design_tw <- model.matrix(~ Twin_Group + diagnosis )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))## (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down 246826 2327 2873 2500 1993 4380
## NotSig 54341 786065 785865 786193 786773 783884
## Up 489491 2266 1920 1965 1892 2394
## Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13 diagnosis
## Down 44847 2593 2859 2007 2569 0
## NotSig 733601 785482 785920 786577 786393 790658
## Up 12210 2583 1879 2074 1696 0top <- 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] 17072output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_guthrie_diagnosis.csv",row.names=FALSE)
head(output,30) %>% kbl(caption="diagnosis") %>% kable_paper("hover", full_width = F)| Name | chr | pos | strand | Islands_Name | Relation_to_Island | UCSC_RefGene_Name | UCSC_RefGene_Group | logFC | AveExpr | t | P.Value | adj.P.Val | B |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| cg05281810 | chr18 | 77542644 |
|
chr18:77542920-77543323 | N_Shore | -0.3454334 | 3.0991938 | -6.251185 | 0.0000251 | 0.986885 | -1.507792 | ||
| cg14507658 | chr21 | 46408142 |
|
chr21:46407925-46408143 | Island | 0.5668544 | 3.3429264 | 6.114589 | 0.0000315 | 0.986885 | -1.558126 | ||
| cg11355135 | chr2 | 108602860 |
|
chr2:108602824-108603467 | Island | SLC5A7 | TSS200 | 0.6331678 | -3.8213288 | 5.890972 | 0.0000457 | 0.986885 | -1.644564 |
| cg08008144 | chr8 | 1983127 |
|
OpenSea | -0.5202509 | 2.5733842 | -5.741776 | 0.0000589 | 0.986885 | -1.705144 | |||
| cg02344606 | chr7 | 112062297 |
|
OpenSea | IFRD1 | TSS1500 | 0.4802282 | -2.0952319 | 5.709395 | 0.0000623 | 0.986885 | -1.718609 | |
| cg00633277 | chr12 | 7250385 |
|
OpenSea | C1RL;C1RL;C1RL | Body;Body;Body | -0.4824283 | 0.7766609 | -5.478184 | 0.0000929 | 0.986885 | -1.818150 | |
| cg15740041 | chr16 | 87781336 |
|
OpenSea | KLHDC4;KLHDC4;KLHDC4 | Body;Body;Body | 0.6792984 | 1.5969029 | 5.448471 | 0.0000979 | 0.986885 | -1.831383 | |
| cg25272125 | chr1 | 31538096 |
|
chr1:31538412-31538938 | N_Shore | PUM1;PUM1 | 5’UTR;5’UTR | 0.4307488 | -5.3095399 | 5.381241 | 0.0001101 | 0.986885 | -1.861703 |
| cg24267485 | chr18 | 74766743 |
|
chr18:74770351-74770590 | N_Shelf | MBP;MBP | Body;Body | 0.5554228 | -2.8111305 | 5.379639 | 0.0001104 | 0.986885 | -1.862433 |
| cg14529256 | chr1 | 245843607 |
|
OpenSea | KIF26B | Body | 0.5951160 | 2.2974986 | 5.350759 | 0.0001162 | 0.986885 | -1.875625 | |
| cg11712265 | chr19 | 7968106 |
|
chr19:7968325-7969421 | N_Shore | MAP2K7;MAP2K7;MAP2K7 | TSS1500;TSS1500;TSS1500 | 0.3440832 | -3.7469755 | 5.309199 | 0.0001251 | 0.986885 | -1.894784 |
| cg27149150 | chr1 | 7449851 |
|
OpenSea | CAMTA1 | Body | -0.4878662 | -0.3789890 | -5.251698 | 0.0001385 | 0.986885 | -1.921632 | |
| cg16239628 | chr13 | 46626330 |
|
chr13:46626454-46627019 | N_Shore | ZC3H13 | 5’UTR | -0.5714673 | -4.6134867 | -5.238885 | 0.0001417 | 0.986885 | -1.927668 |
| cg08277291 | chr6 | 32037426 |
|
OpenSea | TNXB | Body | -0.4220068 | 3.0312857 | -5.220157 | 0.0001465 | 0.986885 | -1.936527 | |
| cg04900486 | chr1 | 151320117 |
|
chr1:151319326-151319545 | S_Shore | RFX5;RFX5 | TSS1500;TSS1500 | 0.4805099 | -4.3378613 | 5.210133 | 0.0001491 | 0.986885 | -1.941286 |
| cg23404427 | chr19 | 1653402 |
|
OpenSea | TCF3 | TSS1500 | -0.4918227 | -5.0274142 | -5.198232 | 0.0001523 | 0.986885 | -1.946952 | |
| cg17904575 | chr14 | 102290370 |
|
OpenSea | PPP2R5C;PPP2R5C;PPP2R5C;PPP2R5C;PPP2R5C | Body;Body;Body;Body;Body | -0.3908957 | 1.0365829 | -5.193349 | 0.0001536 | 0.986885 | -1.949282 | |
| cg08471713 | chr17 | 41738893 |
|
OpenSea | MEOX1;MEOX1;MEOX1 | 1stExon;1stExon;5’UTR | -0.5405729 | -2.6815531 | -5.087137 | 0.0001858 | 0.986885 | -2.000680 | |
| cg05116656 | chr5 | 94417892 |
|
OpenSea | MCTP1;MCTP1;MCTP1 | TSS1500;TSS1500;Body | 0.4626864 | -3.3779970 | 5.077437 | 0.0001891 | 0.986885 | -2.005444 | |
| cg06595448 | chr7 | 50727315 |
|
OpenSea | GRB10;GRB10;GRB10;GRB10 | Body;Body;Body;Body | 0.3855961 | 2.2886258 | 5.015083 | 0.0002116 | 0.986885 | -2.036343 | |
| cg01497262 | chr2 | 43635730 |
|
OpenSea | THADA;THADA;THADA | Body;Body;Body | 0.5221124 | 0.3162170 | 4.988724 | 0.0002220 | 0.986885 | -2.049552 | |
| cg25481455 | chr6 | 22147551 |
|
OpenSea | NBAT1;CASC15 | TSS200;Body | 0.4057161 | -4.6752568 | 4.979172 | 0.0002259 | 0.986885 | -2.054360 | |
| cg06753153 | chr14 | 65186982 |
|
OpenSea | PLEKHG3 | 5’UTR | -0.3736818 | 3.1606283 | -4.972126 | 0.0002288 | 0.986885 | -2.057914 | |
| cg25022405 | chr22 | 26137668 |
|
OpenSea | MYO18B | TSS1500 | 0.3428876 | 3.6893419 | 4.969355 | 0.0002299 | 0.986885 | -2.059313 | |
| cg11250279 | chr13 | 24463401 |
|
chr13:24463119-24463677 | Island | C1QTNF9B-AS1;MIPEP;C1QTNF9B-AS1 | 5’UTR;1stExon;TSS200 | 0.3385965 | -4.5582319 | 4.965330 | 0.0002316 | 0.986885 | -2.061347 |
| cg26684767 | chr2 | 231000057 |
|
OpenSea | 0.4303030 | -1.9273125 | 4.957122 | 0.0002351 | 0.986885 | -2.065502 | |||
| cg27315633 | chr21 | 46056170 |
|
OpenSea | KRTAP10-10;C21orf29 | TSS1500;Body | 0.3085857 | 2.7537308 | 4.948176 | 0.0002389 | 0.986885 | -2.070041 | |
| cg09761351 | chr10 | 70092558 |
|
chr10:70091055-70092573 | Island | PBLD;PBLD;HNRNPH3;PBLD;PBLD;HNRNPH3 | 1stExon;5’UTR;5’UTR;1stExon;5’UTR;5’UTR | 0.3265333 | -3.8393948 | 4.944569 | 0.0002405 | 0.986885 | -2.071873 |
| cg00010078 | chr2 | 109967172 |
|
OpenSea | SH3RF3 | Body | 0.3986020 | -0.0707169 | 4.910082 | 0.0002561 | 0.986885 | -2.089480 | |
| cg21175642 | chr3 | 48701770 |
|
chr3:48698335-48701667 | S_Shore | CELSR3 | TSS1500 | 0.5206228 | -3.2141028 | 4.900505 | 0.0002606 | 0.986885 | -2.094395 |
saveRDS(design_tw, "gu_design_diagnosis.rds")Regression analysis - within twin pair on SRS
#SRS:0,1.895198e-07
design_tw <- model.matrix(~ Twin_Group + SRS )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))## (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down 216132 2056 3244 1608 1989 53519
## NotSig 131030 786807 785571 787932 786964 728796
## Up 443496 1795 1843 1118 1705 8343
## Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13 SRS
## Down 203966 2976 3114 1352 1056 0
## NotSig 552849 784673 785492 788198 789016 790658
## Up 33843 3009 2052 1108 586 0top <- topTable(fit2_tw,coef="SRS",num=Inf, sort.by = "P")
head(top)## logFC AveExpr t P.Value adj.P.Val B
## cg19720702 -0.03838474 3.067023 -9.607598 1.905809e-07 0.1506843 7.225640
## cg23884297 0.02678726 2.566698 7.562793 3.067352e-06 0.4590059 4.317512
## cg15865511 0.02908070 2.877240 7.425775 3.760461e-06 0.4590059 4.104410
## cg08281515 -0.04464912 -3.210691 -7.408433 3.859340e-06 0.4590059 4.077265
## cg13031847 -0.03686804 -4.802769 -7.128509 5.899114e-06 0.4590059 3.633617
## cg03767531 -0.03513403 -4.924137 -6.963021 7.617575e-06 0.4590059 3.366441nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)## [1] 85176output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_guthrie_SRS.csv",row.names=FALSE)
head(output,30) %>% kbl(caption="SRS score") %>% kable_paper("hover", full_width = F)| Name | chr | pos | strand | Islands_Name | Relation_to_Island | UCSC_RefGene_Name | UCSC_RefGene_Group | logFC | AveExpr | t | P.Value | adj.P.Val | B |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| cg19720702 | chr1 | 53856531 |
|
OpenSea | -0.0383847 | 3.0670226 | -9.607598 | 2.00e-07 | 0.1506843 | 7.225640 | |||
| cg23884297 | chr20 | 30124811 |
|
OpenSea | HM13;HM13;HM13;HM13 | Body;Body;Body;Body | 0.0267873 | 2.5666976 | 7.562793 | 3.10e-06 | 0.4590059 | 4.317512 | |
| cg15865511 | chr2 | 43252424 |
|
OpenSea | 0.0290807 | 2.8772398 | 7.425775 | 3.80e-06 | 0.4590059 | 4.104410 | |||
| cg08281515 | chr11 | 111956885 |
|
chr11:111957352-111957681 | N_Shore | TIMM8B;SDHD;TIMM8B | Body;TSS1500;Body | -0.0446491 | -3.2106907 | -7.408433 | 3.90e-06 | 0.4590059 | 4.077265 |
| cg13031847 | chr12 | 124873515 |
|
chr12:124873241-124874008 | Island | NCOR2;NCOR2;NCOR2 | Body;Body;Body | -0.0368680 | -4.8027693 | -7.128509 | 5.90e-06 | 0.4590059 | 3.633617 |
| cg03767531 | chr2 | 182545572 |
|
chr2:182547873-182549177 | N_Shelf | NEUROD1 | TSS200 | -0.0351340 | -4.9241373 | -6.963021 | 7.60e-06 | 0.4590059 | 3.366441 |
| cg13534438 | chr3 | 197229854 |
|
chr3:197225631-197226150 | S_Shelf | -0.0254258 | 2.6981194 | -6.922437 | 8.10e-06 | 0.4590059 | 3.300358 | ||
| cg08748969 | chr12 | 69327779 |
|
chr12:69327021-69327532 | S_Shore | CPM;CPM;CPM | TSS1500;TSS1500;5’UTR | 0.0253039 | -1.2281672 | 6.721795 | 1.11e-05 | 0.4590059 | 2.970388 |
| cg13206794 | chr2 | 242010703 |
|
OpenSea | SNED1 | Body | -0.0253873 | 2.0189033 | -6.690690 | 1.17e-05 | 0.4590059 | 2.918747 | |
| cg19868495 | chr13 | 61995249 |
|
OpenSea | 0.0325162 | 0.8837211 | 6.667391 | 1.21e-05 | 0.4590059 | 2.879978 | |||
| cg17916960 | chr15 | 79447300 |
|
OpenSea | 0.0207600 | -2.1672398 | 6.578462 | 1.40e-05 | 0.4590059 | 2.731325 | |||
| cg03013237 | chr16 | 1802436 |
|
OpenSea | MAPK8IP3;MAPK8IP3 | Body;Body | -0.0338807 | 2.3924811 | -6.536434 | 1.50e-05 | 0.4590059 | 2.660697 | |
| cg01201782 | chr16 | 1652449 |
|
OpenSea | IFT140 | Body | -0.0238570 | 3.6739107 | -6.536178 | 1.50e-05 | 0.4590059 | 2.660267 | |
| cg06874119 | chr5 | 1892431 |
|
OpenSea | CTD-2194D22.4 | Body | 0.0246989 | -1.0880647 | 6.532912 | 1.51e-05 | 0.4590059 | 2.654769 | |
| cg15616400 | chr1 | 16062894 |
|
OpenSea | SLC25A34;SLC25A34 | 5’UTR;1stExon | -0.0345763 | 2.2462683 | -6.505202 | 1.57e-05 | 0.4590059 | 2.608056 | |
| cg17497176 | chr7 | 150763829 |
|
OpenSea | SLC4A2 | Body | -0.0321668 | 2.6216054 | -6.384151 | 1.92e-05 | 0.4590059 | 2.402772 | |
| cg21401292 | chr6 | 134570669 |
|
OpenSea | SGK1 | Body | 0.0262834 | 1.8889371 | 6.314084 | 2.15e-05 | 0.4590059 | 2.283034 | |
| cg16134349 | chr10 | 61900552 |
|
OpenSea | ANK3;ANK3;ANK3;ANK3;ANK3 | 1stExon;5’UTR;Body;Body;Body | -0.0194964 | 2.8568108 | -6.274146 | 2.30e-05 | 0.4590059 | 2.214485 | |
| cg22593006 | chr6 | 170603682 |
|
chr6:170604411-170606479 | N_Shore | FAM120B;FAM120B | TSS1500;Body | 0.0301467 | -0.7098248 | 6.268670 | 2.32e-05 | 0.4590059 | 2.205070 |
| cg13385114 | chr1 | 9656493 |
|
OpenSea | TMEM201;TMEM201 | Body;Body | -0.0318449 | 2.3799036 | -6.233960 | 2.45e-05 | 0.4590059 | 2.145291 | |
| cg09487134 | chr10 | 35195634 |
|
OpenSea | -0.0193997 | 0.2799785 | -6.162045 | 2.77e-05 | 0.4590059 | 2.020915 | |||
| cg11920999 | chr22 | 37417225 |
|
chr22:37414319-37416111 | S_Shore | MPST;MPST;TST;MPST | 5’UTR;5’UTR;TSS1500;Body | 0.0254532 | 3.2047572 | 6.161423 | 2.77e-05 | 0.4590059 | 2.019836 |
| cg14266530 | chr13 | 52511468 |
|
OpenSea | ATP7B;ATP7B;ATP7B | Body;Body;Body | -0.0285632 | 3.6982464 | -6.155577 | 2.80e-05 | 0.4590059 | 2.009694 | |
| cg13709529 | chr6 | 54001895 |
|
OpenSea | C6orf142 | Body | 0.0350504 | 2.9909831 | 6.108476 | 3.02e-05 | 0.4590059 | 1.927812 | |
| cg01038640 | chr5 | 157097845 |
|
chr5:157098263-157099041 | N_Shore | C5orf52 | TSS1500 | -0.0344498 | -3.1668696 | -6.106252 | 3.04e-05 | 0.4590059 | 1.923938 |
| cg23657393 | chr20 | 42790697 |
|
chr20:42788216-42789168 | S_Shore | JPH2 | Body | -0.0279981 | 2.5606519 | -6.104255 | 3.05e-05 | 0.4590059 | 1.920458 |
| cg12446246 | chr1 | 208418552 |
|
chr1:208416489-208418016 | S_Shore | PLXNA2 | TSS1500 | 0.0164020 | -0.4125177 | 6.093169 | 3.10e-05 | 0.4590059 | 1.901135 |
| cg18764121 | chr15 | 25017456 |
|
chr15:25018174-25018533 | N_Shore | -0.0261391 | 0.1134955 | -6.075099 | 3.20e-05 | 0.4590059 | 1.869605 | ||
| cg00647512 | chr8 | 42669499 |
|
OpenSea | -0.0350431 | 2.2195348 | -6.069488 | 3.23e-05 | 0.4590059 | 1.859805 | |||
| cg17594278 | chr12 | 112192763 |
|
OpenSea | ACAD10;ACAD10 | Body;Body | -0.0280782 | 2.5525592 | -6.068771 | 3.23e-05 | 0.4590059 | 1.858553 |
saveRDS(design_tw, "gu_design_srs.rds")Regression analysis - within twin pair on inverse IQ
#iiq:0,1.326967e-05
design_tw <- model.matrix(~ Twin_Group + iiq )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))## (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down 247669 2506 3404 2830 1994 14683
## NotSig 47612 786047 785459 785784 786656 771933
## Up 495377 2105 1795 2044 2008 4042
## Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13 iiq
## Down 183563 2909 2798 2032 3048 0
## NotSig 574299 784699 786188 786643 785864 790658
## Up 32796 3050 1672 1983 1746 0top <- topTable(fit2_tw,coef="iiq",num=Inf, sort.by = "P")
head(top)## logFC AveExpr t P.Value adj.P.Val B
## cg13470919 0.4687461 -3.022728 6.633508 1.302655e-05 0.3987298 -0.01425388
## cg09714307 -0.6297087 3.213380 -6.540738 1.511421e-05 0.3987298 -0.06379907
## cg08210568 0.4942110 -4.367803 6.499328 1.615722e-05 0.3987298 -0.08626969
## cg15934248 0.4696305 -4.797185 6.428510 1.812027e-05 0.3987298 -0.12521213
## cg26869412 -0.5128171 5.367107 -6.250046 2.426541e-05 0.3987298 -0.22629090
## cg07241312 0.5235001 -1.255773 6.077450 3.231910e-05 0.3987298 -0.32817023nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)## [1] 99076output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_guthrie_iIQ.csv",row.names=FALSE)
head(output,30) %>% kbl(caption="inverse IQ") %>% kable_paper("hover", full_width = F)| Name | chr | pos | strand | Islands_Name | Relation_to_Island | UCSC_RefGene_Name | UCSC_RefGene_Group | logFC | AveExpr | t | P.Value | adj.P.Val | B |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| cg13470919 | chr12 | 108955320 |
|
chr12:108956188-108956759 | N_Shore | ISCU;ISCU;SART3 | TSS1500;TSS1500;TSS200 | 0.4687461 | -3.0227285 | 6.633508 | 1.30e-05 | 0.3987298 | -0.0142539 |
| cg09714307 | chr5 | 10483440 |
|
OpenSea | -0.6297087 | 3.2133805 | -6.540738 | 1.51e-05 | 0.3987298 | -0.0637991 | |||
| cg08210568 | chr3 | 15141216 |
|
chr3:15140035-15140890 | S_Shore | ZFYVE20 | TSS1500 | 0.4942110 | -4.3678027 | 6.499328 | 1.62e-05 | 0.3987298 | -0.0862697 |
| cg15934248 | chr2 | 86332749 |
|
chr2:86332659-86333650 | Island | POLR1A;PTCD3 | Body;TSS1500 | 0.4696305 | -4.7971845 | 6.428510 | 1.81e-05 | 0.3987298 | -0.1252121 |
| cg26869412 | chr10 | 3023801 |
|
OpenSea | -0.5128171 | 5.3671068 | -6.250045 | 2.43e-05 | 0.3987298 | -0.2262909 | |||
| cg07241312 | chr16 | 87051295 |
|
OpenSea | 0.5235001 | -1.2557727 | 6.077450 | 3.23e-05 | 0.3987298 | -0.3281702 | |||
| cg00758881 | chr16 | 58534681 |
|
chr16:58535040-58535596 | N_Shore | NDRG4;NDRG4;NDRG4 | Body;Body;Body | -0.3803378 | -4.2022002 | -6.028860 | 3.51e-05 | 0.3987298 | -0.3576032 |
| cg00511674 | chr16 | 78080068 |
|
chr16:78079752-78080166 | Island | 0.3844774 | -5.2143180 | 6.011863 | 3.61e-05 | 0.3987298 | -0.3679777 | ||
| cg05845592 | chr16 | 28634766 |
|
chr16:28634440-28635031 | Island | SULT1A1;SULT1A1 | 5’UTR;1stExon | 0.6165195 | -5.6347444 | 5.952190 | 3.99e-05 | 0.3987298 | -0.4047294 |
| cg15957908 | chr6 | 158923478 |
|
OpenSea | TULP4;TULP4 | Body;Body | -0.3860089 | 3.6925825 | -5.944325 | 4.04e-05 | 0.3987298 | -0.4096113 | |
| cg13709580 | chr5 | 10668678 |
|
OpenSea | -0.3633339 | 2.2322906 | -5.883964 | 4.48e-05 | 0.3987298 | -0.4473784 | |||
| cg16336665 | chr17 | 73028495 |
|
chr17:73030676-73031160 | N_Shelf | KCTD2 | TSS200 | 0.7375034 | -2.2349760 | 5.873733 | 4.56e-05 | 0.3987298 | -0.4538322 |
| cg07943548 | chr15 | 86304357 |
|
chr15:86301949-86302213 | S_Shelf | KLHL25 | 3’UTR | -0.4948093 | 4.2305474 | -5.835656 | 4.86e-05 | 0.3987298 | -0.4779868 |
| cg09375488 | chr7 | 95115026 |
|
OpenSea | ASB4;ASB4 | TSS1500;TSS1500 | -0.4622716 | 0.2950509 | -5.825875 | 4.94e-05 | 0.3987298 | -0.4842258 | |
| cg23905944 | chr20 | 62571879 |
|
chr20:62571738-62572556 | Island | UCKL1 | Body | -0.5873078 | 3.4203253 | -5.805162 | 5.12e-05 | 0.3987298 | -0.4974851 |
| cg21419752 | chr3 | 139514052 |
|
OpenSea | -0.4104916 | 1.0266672 | -5.800482 | 5.16e-05 | 0.3987298 | -0.5004897 | |||
| cg25849262 | chr3 | 19189070 |
|
chr3:19189688-19190100 | N_Shore | KCNH8 | TSS1500 | 0.4570049 | -4.5601102 | 5.773415 | 5.41e-05 | 0.3987298 | -0.5179308 |
| cg14570121 | chr19 | 45874153 |
|
chr19:45873339-45874033 | S_Shore | ERCC2;ERCC2 | TSS1500;TSS1500 | 0.3972660 | -4.7537124 | 5.761779 | 5.52e-05 | 0.3987298 | -0.5254621 |
| cg07690326 | chr8 | 145019032 |
|
chr8:145018815-145019214 | Island | PLEC1;PLEC1;PLEC1;PLEC1;PLEC1;PLEC1 | Body;Body;Body;TSS1500;Body;TSS200 | -0.5396969 | 5.3541477 | -5.761448 | 5.52e-05 | 0.3987298 | -0.5256764 |
| cg19686637 | chr15 | 85427743 |
|
OpenSea | SLC28A1;SLC28A1 | TSS200;TSS200 | -0.3461389 | 3.2554165 | -5.755249 | 5.58e-05 | 0.3987298 | -0.5296973 | |
| cg26362488 | chr12 | 14924051 |
|
chr12:14922877-14923974 | S_Shore | HIST4H4;HIST4H4 | 1stExon;5’UTR | 0.8217243 | -4.6311437 | 5.753558 | 5.60e-05 | 0.3987298 | -0.5307954 |
| cg03873049 | chr5 | 149625194 |
|
OpenSea | CAMK2A;CAMK2A | Body;Body | -0.3790098 | 3.3047912 | -5.733738 | 5.79e-05 | 0.3987298 | -0.5436944 | |
| cg23644589 | chr1 | 228404256 |
|
chr1:228404148-228404397 | Island | OBSCN;OBSCN | Body;Body | -0.6045446 | 3.8354376 | -5.715842 | 5.97e-05 | 0.3987298 | -0.5553921 |
| cg21645826 | chr19 | 19313533 |
|
chr19:19314137-19314416 | N_Shore | NR2C2AP | Body | -0.4030073 | -2.4552364 | -5.710733 | 6.02e-05 | 0.3987298 | -0.5587407 |
| cg19166616 | chr20 | 62259876 |
|
chr20:62257589-62259476 | S_Shore | GMEB2 | TSS1500 | -0.4375938 | -3.4141333 | -5.710294 | 6.03e-05 | 0.3987298 | -0.5590287 |
| cg07019009 | chr15 | 102265803 |
|
chr15:102264097-102264750 | S_Shore | TARSL2 | TSS1500 | -0.5622911 | 1.0209687 | -5.703379 | 6.10e-05 | 0.3987298 | -0.5635668 |
| cg07502168 | chr1 | 214725311 |
|
chr1:214725192-214725528 | Island | PTPN14 | TSS1500 | 0.5791304 | -4.8125449 | 5.683416 | 6.31e-05 | 0.3987298 | -0.5767098 |
| cg04465078 | chr16 | 15766942 |
|
OpenSea | NDE1;NDE1 | Body;Body | 0.9766506 | -2.3211134 | 5.672261 | 6.44e-05 | 0.3987298 | -0.5840807 | |
| cg21353154 | chr13 | 99381820 |
|
OpenSea | SLC15A1 | Body | -0.4790397 | -0.3292831 | -5.670635 | 6.46e-05 | 0.3987298 | -0.5851567 | |
| cg03734595 | chr18 | 9333974 |
|
chr18:9333896-9334244 | Island | TWSG1 | TSS1500 | 0.4850929 | -4.4676286 | 5.658297 | 6.60e-05 | 0.3987298 | -0.5933334 |
saveRDS(design_tw, "gu_design_iiq.rds")Regression analysis - within twin pair on inverse language skills
#ilanguage:0,9.149182e-08
design_tw <- model.matrix(~ Twin_Group + ilanguage )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))## (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down 247378 1824 2268 2296 2004 5873
## NotSig 52181 787219 786885 786623 786779 782005
## Up 491099 1615 1505 1739 1875 2780
## Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13 ilanguage
## Down 138798 2733 2464 1795 2794 0
## NotSig 624400 785112 786597 787118 786104 790658
## Up 27460 2813 1597 1745 1760 0top <- topTable(fit2_tw,coef="ilanguage",num=Inf, sort.by = "P")
head(top)## logFC AveExpr t P.Value adj.P.Val B
## cg22165507 -0.7962383 -4.455296 -10.326679 8.885287e-08 0.07025224 -0.6822315
## cg13435526 -0.5163639 3.222080 -7.544551 3.388991e-06 0.99011776 -1.1824312
## cg25679269 0.9220801 -4.314887 6.932603 8.514489e-06 0.99011776 -1.3486002
## cg22533406 -0.6174507 4.882477 -6.902615 8.918888e-06 0.99011776 -1.3574555
## cg19107578 0.5399238 2.688710 6.731461 1.164995e-05 0.99011776 -1.4093703
## cg18698884 0.7116987 -3.398486 6.560063 1.528340e-05 0.99011776 -1.4637874nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)## [1] 28188output <-merge(ann_sub,top,by.x="Name",by.y="row.names")
output <- output[order(output$P.Value),]
write.csv(output,file="limma_guthrie_ilanguage.csv",row.names=FALSE)
head(output,30) %>% kbl(caption="inverse language") %>% kable_paper("hover", full_width = F)| Name | chr | pos | strand | Islands_Name | Relation_to_Island | UCSC_RefGene_Name | UCSC_RefGene_Group | logFC | AveExpr | t | P.Value | adj.P.Val | B |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| cg22165507 | chr1 | 165601007 |
|
chr1:165599563-165600574 | S_Shore | MGST3 | 5’UTR | -0.7962383 | -4.455296 | -10.326679 | 1.00e-07 | 0.0702522 | -0.6822315 |
| cg13435526 | chr8 | 22472050 |
|
OpenSea | CCAR2;CCAR2 | Body;Body | -0.5163639 | 3.222080 | -7.544551 | 3.40e-06 | 0.9901178 | -1.1824312 | |
| cg25679269 | chr11 | 118747550 |
|
OpenSea | 0.9220801 | -4.314887 | 6.932603 | 8.50e-06 | 0.9901178 | -1.3486002 | |||
| cg22533406 | chr19 | 49106849 |
|
chr19:49106882-49107178 | N_Shore | FAM83E | Body | -0.6174507 | 4.882477 | -6.902615 | 8.90e-06 | 0.9901178 | -1.3574555 |
| cg19107578 | chr5 | 493262 |
|
chr5:494803-495306 | N_Shore | SLC9A3 | Body | 0.5399238 | 2.688710 | 6.731461 | 1.16e-05 | 0.9901178 | -1.4093703 |
| cg18698884 | chr11 | 16632724 |
|
chr11:16632508-16632725 | Island | 0.7116987 | -3.398486 | 6.560063 | 1.53e-05 | 0.9901178 | -1.4637874 | ||
| cg10021749 | chr7 | 2557358 |
|
chr7:2558371-2559967 | N_Shore | LFNG;LFNG | Body;TSS200 | -0.5743078 | 2.793587 | -6.491209 | 1.71e-05 | 0.9901178 | -1.4863586 |
| cg12843467 | chr15 | 99190301 |
|
chr15:99190446-99194559 | N_Shore | IRAIN;IGF1R;IGF1R | Body;TSS1500;TSS1500 | -0.4933420 | -4.969983 | -6.413177 | 1.93e-05 | 0.9901178 | -1.5124454 |
| cg13938969 | chr15 | 29409443 |
|
chr15:29407695-29408229 | S_Shore | APBA2;APBA2 | 3’UTR;3’UTR | -0.5885312 | 2.757445 | -6.379493 | 2.04e-05 | 0.9901178 | -1.5238751 |
| cg13724550 | chr5 | 142001033 |
|
OpenSea | FGF1;FGF1;FGF1;FGF1;FGF1;FGF1;FGF1 | 5’UTR;TSS200;5’UTR;5’UTR;5’UTR;Body;Body | -0.5420348 | 3.488093 | -6.340121 | 2.18e-05 | 0.9901178 | -1.5373664 | |
| cg23404435 | chr7 | 1588396 |
|
chr7:1590387-1591020 | N_Shore | TMEM184A | Body | -0.5821756 | 2.408359 | -6.326735 | 2.23e-05 | 0.9901178 | -1.5419853 |
| cg20781140 | chr2 | 175352216 |
|
chr2:175351364-175352154 | S_Shore | GPR155;GPR155 | TSS1500;TSS1500 | 0.5758823 | -4.024754 | 6.302035 | 2.32e-05 | 0.9901178 | -1.5505524 |
| cg01369908 | chr11 | 34375707 |
|
chr11:34378229-34379993 | N_Shelf | ABTB2 | Body | -0.4642938 | 3.183633 | -6.187706 | 2.79e-05 | 0.9901178 | -1.5909508 |
| cg03547220 | chr19 | 7069348 |
|
chr19:7069396-7069921 | N_Shore | ZNF557;ZNF557;ZNF557 | TSS200;TSS200;TSS200 | 0.5613098 | -2.014080 | 6.179544 | 2.83e-05 | 0.9901178 | -1.5938823 |
| cg01231620 | chr14 | 23298884 |
|
chr14:23298984-23299313 | N_Shore | MRPL52;MRPL52;MRPL52;MRPL52;MRPL52;MRPL52 | TSS1500;TSS1500;TSS1500;TSS1500;TSS1500;TSS1500 | 0.5292737 | -4.349346 | 6.101669 | 3.22e-05 | 0.9901178 | -1.6221757 |
| cg11913416 | chr1 | 13909012 |
|
chr1:13910137-13910868 | N_Shore | PDPN;PDPN | TSS1500;TSS1500 | 0.4991428 | -0.702180 | 6.055432 | 3.48e-05 | 0.9901178 | -1.6392545 |
| cg09219421 | chr19 | 13779879 |
|
OpenSea | -0.4056121 | 2.949086 | -6.007212 | 3.77e-05 | 0.9901178 | -1.6572913 | |||
| cg02882857 | chr5 | 149381215 |
|
chr5:149379964-149380802 | S_Shore | HMGXB3;TIGD6;TIGD6 | 5’UTR;TSS1500;TSS1500 | 0.6416605 | -4.191304 | 6.005829 | 3.78e-05 | 0.9901178 | -1.6578122 |
| cg27594834 | chr19 | 11708137 |
|
chr19:11708277-11708683 | N_Shore | ZNF627 | TSS200 | -0.5826615 | -4.511569 | -5.861055 | 4.82e-05 | 0.9901178 | -1.7133981 |
| cg09459188 | chr3 | 158520176 |
|
chr3:158519509-158520329 | Island | MFSD1;MFSD1 | Body;Body | 0.8660146 | -4.773250 | 5.815394 | 5.21e-05 | 0.9901178 | -1.7313778 |
| cg17506328 | chr11 | 46868210 |
|
chr11:46867194-46868077 | S_Shore | CKAP5;CKAP5 | TSS1500;TSS1500 | -0.3988351 | -3.441229 | -5.780347 | 5.53e-05 | 0.9901178 | -1.7453263 |
| cg27583671 | chr6 | 151670624 |
|
OpenSea | AKAP12;AKAP12 | Body;Body | -0.5433919 | 2.656961 | -5.767792 | 5.65e-05 | 0.9901178 | -1.7503550 | |
| cg27522733 | chr1 | 193448591 |
|
OpenSea | 0.5176401 | -2.434439 | 5.731087 | 6.01e-05 | 0.9901178 | -1.7651516 | |||
| cg19572637 | chr2 | 66660175 |
|
chr2:66660452-66660794 | N_Shore | 0.3264399 | -4.483067 | 5.718716 | 6.14e-05 | 0.9901178 | -1.7701716 | ||
| cg06848865 | chr14 | 82000050 |
|
chr14:81999751-82000434 | Island | SEL1L | 1stExon | -0.7715305 | -4.890653 | -5.696292 | 6.38e-05 | 0.9901178 | -1.7793122 |
| cg09733297 | chr16 | 54155568 |
|
chr16:54155953-54156180 | N_Shore | -0.5844263 | 3.384625 | -5.675573 | 6.61e-05 | 0.9901178 | -1.7878062 | ||
| cg17427349 | chr15 | 38988004 |
|
OpenSea | C15orf53 | TSS1500 | -0.4387749 | 3.003292 | -5.673513 | 6.64e-05 | 0.9901178 | -1.7886529 | |
| cg05792660 | chr14 | 94929320 |
|
OpenSea | SERPINA9;SERPINA9;SERPINA9;SERPINA9 | 3’UTR;3’UTR;3’UTR;3’UTR | -0.3724785 | 2.242181 | -5.649184 | 6.92e-05 | 0.9901178 | -1.7986917 | |
| cg01456296 | chr12 | 54971406 |
|
OpenSea | PDE1B;PDE1B;PDE1B;PDE1B | 3’UTR;3’UTR;3’UTR;3’UTR | -0.3758029 | 3.390950 | -5.638268 | 7.05e-05 | 0.9901178 | -1.8032171 | |
| cg10078645 | chr12 | 124903668 |
|
chr12:124901554-124901843 | S_Shore | NCOR2;NCOR2 | Body;Body | -0.4205931 | 3.786148 | -5.615006 | 7.34e-05 | 0.9901178 | -1.8129036 |
saveRDS(design_tw, "gu_design_ilanguage.rds")Converting to beta from M values
bDat = ilogit2(mvals)
bDat_new = getBeta(gmset_flt)
#View(bDat)c
write.csv(bDat,file="ASD_guthrie_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 ADOSmylm <- lm(gwam~ADOS)
summary(mylm)##
## Call:
## lm(formula = gwam ~ ADOS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.123842 -0.021378 -0.009516 0.044810 0.074402
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7478296 0.0266803 28.029 <2e-16 ***
## ADOS -0.0001873 0.0021538 -0.087 0.932
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05386 on 20 degrees of freedom
## Multiple R-squared: 0.0003778, Adjusted R-squared: -0.0496
## F-statistic: 0.00756 on 1 and 20 DF, p-value: 0.9316summary(lm(gwam ~ Twin_Group + ADOS ))##
## Call:
## lm(formula = gwam ~ Twin_Group + ADOS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.04877 -0.02110 0.00000 0.02110 0.04877
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7850380 0.0323414 24.273 3.21e-10 ***
## Twin_Group2 -0.0188145 0.0393075 -0.479 0.6425
## Twin_Group3 -0.0473057 0.0409223 -1.156 0.2746
## Twin_Group4 -0.0417615 0.0390814 -1.069 0.3104
## Twin_Group5 -0.0058695 0.0390814 -0.150 0.8836
## Twin_Group6 -0.0881363 0.0409223 -2.154 0.0567 .
## Twin_Group7 -0.1214234 0.0413367 -2.937 0.0149 *
## Twin_Group8 0.0290376 0.0469005 0.619 0.5497
## Twin_Group9 -0.0139663 0.0384659 -0.363 0.7241
## Twin_Group12 0.0322242 0.0417797 0.771 0.4584
## Twin_Group13 -0.0490928 0.0384659 -1.276 0.2307
## ADOS -0.0008717 0.0023356 -0.373 0.7168
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03845 on 10 degrees of freedom
## Multiple R-squared: 0.7453, Adjusted R-squared: 0.4651
## F-statistic: 2.66 on 11 and 10 DF, p-value: 0.06726message("GWAM assocation with SRS")## GWAM assocation with SRSmylm <- lm(gwam~SRS)
summary(mylm)##
## Call:
## lm(formula = gwam ~ SRS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.12381 -0.02376 -0.00655 0.04623 0.07058
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7563524 0.0304976 24.800 <2e-16 ***
## SRS -0.0001336 0.0003558 -0.376 0.711
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05368 on 20 degrees of freedom
## Multiple R-squared: 0.007003, Adjusted R-squared: -0.04265
## F-statistic: 0.141 on 1 and 20 DF, p-value: 0.7112plot(gwam~SRS)
abline(mylm)
summary(lm(gwam ~ Twin_Group + SRS ))##
## Call:
## lm(formula = gwam ~ Twin_Group + SRS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.03984 -0.01787 0.00000 0.01787 0.03984
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8858552 0.0649284 13.644 8.66e-08 ***
## Twin_Group2 -0.0664369 0.0420404 -1.580 0.14512
## Twin_Group3 -0.0793746 0.0369698 -2.147 0.05735 .
## Twin_Group4 -0.1004921 0.0482648 -2.082 0.06397 .
## Twin_Group5 -0.0324477 0.0363394 -0.893 0.39289
## Twin_Group6 -0.1058273 0.0344674 -3.070 0.01183 *
## Twin_Group7 -0.1491355 0.0359534 -4.148 0.00199 **
## Twin_Group8 0.0146998 0.0338343 0.434 0.67318
## Twin_Group9 -0.0120059 0.0337741 -0.355 0.72962
## Twin_Group12 -0.0347440 0.0480725 -0.723 0.48640
## Twin_Group13 -0.1377999 0.0604396 -2.280 0.04579 *
## SRS -0.0009585 0.0005391 -1.778 0.10579
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03375 on 10 degrees of freedom
## Multiple R-squared: 0.8037, Adjusted R-squared: 0.5879
## F-statistic: 3.723 on 11 and 10 DF, p-value: 0.02372message("GWAM assocation with motor impairment")## GWAM assocation with motor impairmentmylm <- lm(gwam~motor)
summary(mylm)##
## Call:
## lm(formula = gwam ~ motor)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.120972 -0.021477 -0.009074 0.046560 0.072041
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.747943 0.015886 47.080 <2e-16 ***
## motor -0.002429 0.012088 -0.201 0.843
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05381 on 20 degrees of freedom
## Multiple R-squared: 0.002014, Adjusted R-squared: -0.04789
## F-statistic: 0.04037 on 1 and 20 DF, p-value: 0.8428summary(lm(gwam ~ Twin_Group + motor ))##
## Call:
## lm(formula = gwam ~ Twin_Group + motor)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.05094 -0.01497 0.00000 0.01497 0.05094
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.84321 0.05678 14.851 3.85e-08 ***
## Twin_Group2 -0.05422 0.04398 -1.233 0.24584
## Twin_Group3 -0.08489 0.04398 -1.930 0.08241 .
## Twin_Group4 -0.10385 0.06220 -1.670 0.12592
## Twin_Group5 -0.07319 0.06220 -1.177 0.26653
## Twin_Group6 -0.15807 0.06220 -2.542 0.02929 *
## Twin_Group7 -0.12709 0.03591 -3.539 0.00536 **
## Twin_Group8 0.01901 0.03591 0.529 0.60803
## Twin_Group9 -0.01440 0.03591 -0.401 0.69681
## Twin_Group12 -0.03859 0.06220 -0.620 0.54888
## Twin_Group13 -0.11337 0.06220 -1.823 0.09835 .
## motor -0.03235 0.02539 -1.274 0.23141
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03591 on 10 degrees of freedom
## Multiple R-squared: 0.7778, Adjusted R-squared: 0.5334
## F-statistic: 3.182 on 11 and 10 DF, p-value: 0.03936message("GWAM assocation with diagnosis")## GWAM assocation with diagnosismylm <- lm(gwam~diagnosis)
summary(mylm)##
## Call:
## lm(formula = gwam ~ diagnosis)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.113165 -0.023960 -0.006969 0.051259 0.070328
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.76403 0.02087 36.606 <2e-16 ***
## diagnosis -0.01438 0.01384 -1.038 0.311
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05247 on 20 degrees of freedom
## Multiple R-squared: 0.05116, Adjusted R-squared: 0.003721
## F-statistic: 1.078 on 1 and 20 DF, p-value: 0.3114summary(lm(gwam ~ Twin_Group + diagnosis ))##
## Call:
## lm(formula = gwam ~ Twin_Group + diagnosis)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.04545 -0.02117 0.00000 0.02117 0.04545
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.789484 0.029148 27.085 1.09e-10 ***
## Twin_Group2 -0.016374 0.037801 -0.433 0.6741
## Twin_Group3 -0.047044 0.037801 -1.245 0.2417
## Twin_Group4 -0.044638 0.037801 -1.181 0.2650
## Twin_Group5 -0.008485 0.037286 -0.228 0.8246
## Twin_Group6 -0.087875 0.037801 -2.325 0.0424 *
## Twin_Group7 -0.116106 0.039303 -2.954 0.0144 *
## Twin_Group8 0.029997 0.039303 0.763 0.4630
## Twin_Group9 -0.019894 0.037801 -0.526 0.6102
## Twin_Group12 0.031614 0.037801 0.836 0.4225
## Twin_Group13 -0.048657 0.037286 -1.305 0.2211
## diagnosis -0.010983 0.012429 -0.884 0.3976
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03729 on 10 degrees of freedom
## Multiple R-squared: 0.7604, Adjusted R-squared: 0.4969
## F-statistic: 2.886 on 11 and 10 DF, p-value: 0.05303message("GWAM assocation with iIQ")## GWAM assocation with iIQmylm <- lm(gwam~iiq)
summary(mylm)##
## Call:
## lm(formula = gwam ~ iiq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.121210 -0.027717 0.001123 0.037264 0.071896
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.74574 0.01118 66.716 <2e-16 ***
## iiq -0.01207 0.01144 -1.055 0.304
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05243 on 20 degrees of freedom
## Multiple R-squared: 0.05268, Adjusted R-squared: 0.005316
## F-statistic: 1.112 on 1 and 20 DF, p-value: 0.3042plot(gwam~iiq)
abline(mylm)
summary(lm(gwam ~ Twin_Group + iiq ))##
## Call:
## lm(formula = gwam ~ Twin_Group + iiq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.03422 -0.01737 0.00000 0.01737 0.03422
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.788777 0.023785 33.163 1.47e-11 ***
## Twin_Group2 -0.048601 0.035499 -1.369 0.20093
## Twin_Group3 -0.081506 0.035944 -2.268 0.04676 *
## Twin_Group4 -0.051835 0.033445 -1.550 0.15222
## Twin_Group5 0.001069 0.033178 0.032 0.97493
## Twin_Group6 -0.094183 0.032827 -2.869 0.01669 *
## Twin_Group7 -0.128585 0.032833 -3.916 0.00288 **
## Twin_Group8 0.022520 0.032872 0.685 0.50886
## Twin_Group9 -0.041629 0.035594 -1.170 0.26930
## Twin_Group12 0.009342 0.033903 0.276 0.78850
## Twin_Group13 -0.060053 0.033326 -1.802 0.10172
## iiq -0.019321 0.009769 -1.978 0.07616 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03282 on 10 degrees of freedom
## Multiple R-squared: 0.8143, Adjusted R-squared: 0.6101
## F-statistic: 3.987 on 11 and 10 DF, p-value: 0.01883message("GWAM assocation with ilanguage")## GWAM assocation with ilanguagemylm <- lm(gwam~ilanguage)
summary(mylm)##
## Call:
## lm(formula = gwam ~ ilanguage)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.122403 -0.026473 -0.004192 0.038588 0.071380
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.745736 0.011352 65.693 <2e-16 ***
## ilanguage -0.007961 0.011619 -0.685 0.501
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05325 on 20 degrees of freedom
## Multiple R-squared: 0.02293, Adjusted R-squared: -0.02592
## F-statistic: 0.4694 on 1 and 20 DF, p-value: 0.5011plot(gwam~ilanguage)
abline(mylm)
summary(lm(gwam ~ Twin_Group + ilanguage ))##
## Call:
## lm(formula = gwam ~ Twin_Group + ilanguage)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.04184 -0.01944 0.00000 0.01944 0.04184
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.786167 0.027610 28.474 6.64e-11 ***
## Twin_Group2 -0.042759 0.043673 -0.979 0.35064
## Twin_Group3 -0.071806 0.042768 -1.679 0.12408
## Twin_Group4 -0.051054 0.039417 -1.295 0.22434
## Twin_Group5 -0.009145 0.037203 -0.246 0.81079
## Twin_Group6 -0.096835 0.037390 -2.590 0.02695 *
## Twin_Group7 -0.127220 0.037196 -3.420 0.00654 **
## Twin_Group8 0.024651 0.037705 0.654 0.52801
## Twin_Group9 -0.030154 0.041003 -0.735 0.47898
## Twin_Group12 0.011249 0.040608 0.277 0.78741
## Twin_Group13 -0.051676 0.037343 -1.384 0.19652
## ilanguage -0.011556 0.012659 -0.913 0.38277
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0372 on 10 degrees of freedom
## Multiple R-squared: 0.7616, Adjusted R-squared: 0.4993
## F-statistic: 2.904 on 11 and 10 DF, p-value: 0.05204probes_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("Shelf",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_guthrie_ADOS.csv")
TOPPROBES <- head(res$Name,9)
par(mfrow=c(3,3))
sapply(TOPPROBES,effect_plot)
## $cg00546248
## NULL
##
## $cg15740041
## NULL
##
## $cg05116656
## NULL
##
## $cg11632273
## NULL
##
## $cg01656620
## NULL
##
## $cg10881749
## NULL
##
## $cg13699963
## NULL
##
## $cg07162704
## NULL
##
## $cg13277464
## NULLTOPPROBES <- head(res$Name,4)
par(mfrow=c(2,2))
sapply(TOPPROBES,effect_plot)
## $cg00546248
## NULL
##
## $cg15740041
## NULL
##
## $cg05116656
## NULL
##
## $cg11632273
## NULLsig <- subset(res,P.Value<1e-4)
sig <- sig[order(-abs(sig$logFC)),]
TOPPROBES <- head(sig$Name,4)
sapply(TOPPROBES,effect_plot)
## $cg23157290
## NULL
##
## $cg01656620
## NULL
##
## $cg11632273
## NULL
##
## $cg15740041
## NULLpdf("effect_guthrie.pdf")
par(mfrow=c(2,2))
sapply(TOPPROBES,effect_plot)## $cg23157290
## NULL
##
## $cg01656620
## NULL
##
## $cg11632273
## NULL
##
## $cg15740041
## NULLdev.off()## png
## 2par(mfrow=c(1,1))Gene ontology of top DMPs LIMMA data
Gene Ontology analysis (gometh): top 1000 probes (limma data).
#res <- read.csv("ASD_guthrie_top_dmps_onADOS_withintw_Nov27_limma.csv", header = TRUE)
sigCpGs_1k = res$Name[1:1000]
sigCpGs_1k = as.character(sigCpGs_1k)
all = res$Name
length(all)## [1] 790658# kegg
gometh_kegg <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "KEGG", prior.prob=TRUE)## All input CpGs are used for testing.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)## Description N DE P.DE FDR
## hsa04974 Protein digestion and absorption 92 6 0.04927660 1
## hsa05217 Basal cell carcinoma 63 5 0.06384917 1
## hsa04350 TGF-beta signaling pathway 106 6 0.08337315 1
## hsa04930 Type II diabetes mellitus 45 4 0.08569469 1
## hsa05165 Human papillomavirus infection 314 14 0.11804564 1
## hsa04922 Glucagon signaling pathway 99 5 0.12379190 1
## hsa03420 Nucleotide excision repair 57 3 0.12440007 1
## hsa04658 Th1 and Th2 cell differentiation 88 5 0.12923937 1
## hsa04514 Cell adhesion molecules 141 7 0.13258477 1
## hsa00514 Other types of O-glycan biosynthesis 43 3 0.14727719 1
## hsa04211 Longevity regulating pathway 88 5 0.15136069 1
## hsa04927 Cortisol synthesis and secretion 63 4 0.15856029 1
## hsa01240 Biosynthesis of cofactors 147 5 0.17288744 1
## hsa04910 Insulin signaling pathway 131 6 0.18441491 1
## hsa04931 Insulin resistance 105 5 0.19259141 1
## hsa04934 Cushing syndrome 154 7 0.24014001 1
## hsa04120 Ubiquitin mediated proteolysis 134 5 0.25206213 1
## hsa04714 Thermogenesis 213 7 0.25536667 1
## hsa05146 Amoebiasis 98 4 0.26413455 1
## hsa04340 Hedgehog signaling pathway 55 3 0.29100761 1# GO terms
gometh_go <- gometh(sig.cpg = sigCpGs_1k, all.cpg = all, collection = "GO" , prior.prob=TRUE)## All input CpGs are used for testing.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)## ONTOLOGY TERM N
## GO:0006517 BP protein deglycosylation 15
## GO:0048625 BP myoblast fate commitment 6
## GO:0033089 BP positive regulation of T cell differentiation in thymus 10
## GO:0031526 CC brush border membrane 57
## GO:0033077 BP T cell differentiation in thymus 84
## GO:0090160 BP Golgi to lysosome transport 12
## GO:0044319 BP wound healing, spreading of cells 40
## GO:0090505 BP epiboly involved in wound healing 40
## GO:0001837 BP epithelial to mesenchymal transition 182
## GO:0005903 CC brush border 100
## GO:0090504 BP epiboly 41
## GO:0045061 BP thymic T cell selection 24
## GO:0035313 BP wound healing, spreading of epidermal cells 26
## GO:0045059 BP positive thymic T cell selection 15
## GO:0007566 BP embryo implantation 54
## GO:0072578 BP neurotransmitter-gated ion channel clustering 16
## GO:0019827 BP stem cell population maintenance 175
## GO:0045582 BP positive regulation of T cell differentiation 115
## GO:0098727 BP maintenance of cell number 179
## GO:0098644 CC complex of collagen trimers 19
## DE P.DE FDR
## GO:0006517 4 0.0006180811 1
## GO:0048625 3 0.0011160042 1
## GO:0033089 3 0.0028081873 1
## GO:0031526 6 0.0029422814 1
## GO:0033077 8 0.0032663249 1
## GO:0090160 3 0.0040043647 1
## GO:0044319 5 0.0069659613 1
## GO:0090505 5 0.0069659613 1
## GO:0001837 12 0.0072278397 1
## GO:0005903 8 0.0079006136 1
## GO:0090504 5 0.0082426279 1
## GO:0045061 4 0.0092813761 1
## GO:0035313 4 0.0094181597 1
## GO:0045059 3 0.0113805293 1
## GO:0007566 5 0.0123428913 1
## GO:0072578 3 0.0138685763 1
## GO:0019827 11 0.0139224099 1
## GO:0045582 8 0.0148954486 1
## GO:0098727 11 0.0161536542 1
## GO:0098644 3 0.0197022909 1Now 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 790658 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 cachehead(results.ranges, 20)## seqnames start end width strand no.cpgs min_smoothed_fdr
## 1 chr22 38243770 38244199 430 * 3 3.003591e-08
## 2 chr9 24545490 24546166 677 * 4 7.545068e-06
## 3 chr1 46955557 46955679 123 * 2 1.600425e-05
## 4 chr11 22687582 22687660 79 * 5 3.778833e-05
## 5 chr13 113824335 113824550 216 * 2 7.287503e-04
## Stouffer HMFDR Fisher maxdiff meandiff
## 1 0.0019261661 0.0001088322 0.0001388725 0.013216983 4.805928e-03
## 2 0.0041498334 0.0001088322 0.0002364104 0.005518678 -3.043793e-05
## 3 0.0002470517 0.0001911644 0.0001458473 0.012211160 3.655512e-03
## 4 0.0298022656 0.0001088322 0.0007223910 0.015894592 5.082461e-03
## 5 0.0001666470 0.0002415720 0.0001388725 0.013864206 8.027718e-03
## overlapping.genes
## 1 ANKRD54, MIR659
## 2 RP11-20A20.2, IZUMO3
## 3 <NA>
## 4 GAS2
## 5 PROZwrite.csv(results.ranges, file = "dmrcoutput_guthrie.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