Autism Spectrum Disorder Methylation analysis - Guthrie card analysis
Namitha Mohandras, Mark Ziemann
2023-11-10
Source: https://github.com/markziemann/asd_meth
Introduction
Here we are analysing Guthrie card methylation in 9 twin pairs.
baseDir=getwd()
dataDir=paste(baseDir,"/ASD_EPIC_DATA",sep="")
suppressPackageStartupMessages({
library("missMethyl")
library("limma")
library("minfi")
library("IlluminaHumanMethylation450kanno.ilmn12.hg19")
library("IlluminaHumanMethylationEPICanno.ilm10b2.hg19")
library("ruv")
library("RColorBrewer")
library("matrixStats")
library("gplots")
library("WGCNA")
library("FlowSorted.Blood.450k")
library("reshape2")
library("ggplot2")
library("missMethyl")
library("DMRcate")
library("FlowSorted.Blood.EPIC")
library("mitch")
library("kableExtra")
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 = "ASD_guthrie_summarysheet.csv")## [read.metharray.sheet] Found the following CSV files:## [1] "/home/mdz/projects/asd_meth/ASD_EPIC_DATA/ASD_guthrie_summarysheet.csv"#targets$ID = paste(targets$Sample_Group,targets_gen$Sample_Name,sep=".")
rgSet = read.metharray.exp(targets = targets_gen)
sampleNames(rgSet) = targets_gen$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## X11cairo
## 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 = 22, p-value = 0.491
## alternative hypothesis: true location shift is not equal to 0
##
##
## $Bcell
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 36, p-value = 0.3773
## alternative hypothesis: true location shift is not equal to 0
##
##
## $Mono
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 23, p-value = 0.6612
## alternative hypothesis: true location shift is not equal to 0
##
##
## $Gran
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 26, p-value = 0.913
## alternative hypothesis: true location shift is not equal to 0Try estimatecellcounts2.
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.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()## X11cairo
## 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 = 32, p-value = 0.6612
## alternative hypothesis: true location shift is not equal to 0
##
##
## $CD4T
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 36, p-value = 0.3773
## alternative hypothesis: true location shift is not equal to 0
##
##
## $NK
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 19, p-value = 0.3773
## alternative hypothesis: true location shift is not equal to 0
##
##
## $Bcell
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 39, p-value = 0.2212
## alternative hypothesis: true location shift is not equal to 0
##
##
## $Mono
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 20, p-value = 0.4409
## alternative hypothesis: true location shift is not equal to 0
##
##
## $Neu
##
## Wilcoxon rank sum exact test
##
## data: a[, i] and c[, i]
## W = 29, p-value = 0.913
## alternative hypothesis: true location shift is not equal to 0
Filtering
None of the p-values are less than 0.05.
Now filter for high detection p-value and overlap with SNP.
Need to get a copy of the Xreact probes from Namitha. In the mean time I have downloaded from github - link below.
detP <- detectionP(rgSet)
keep <- rowSums(detP < 0.01) == ncol(rgSet)
mset <- mset[keep,]
gmset <- mapToGenome(mset)
#remove SNPs
gmset_flt = dropLociWithSnps(gmset, snps = c("CpG", "SBE"))
#Removing cross-reactive probes
XURL="https://raw.githubusercontent.com/sirselim/illumina450k_filtering/master/EPIC/13059_2016_1066_MOESM1_ESM.csv"
Xreact <- read.csv(XURL)
#Xreact = read.csv(file="/group/canc2/puumba/Data/InfiniumData/NamithaData/Rprojects/Autism/Analysis_Sept11/EPIC_850k_crossreactiveProbes.csv", stringsAsFactors=FALSE)
#Xreact = read.csv(file="~/48639-non-specific-probes-Illumina450k.csv", stringsAsFactors=FALSE)
noXreact <- !(featureNames(gmset) %in% Xreact$X)
gmset <- gmset[noXreact,]
#Removing probes on X and Y chromosomes
autosomes <- !(featureNames(gmset) %in% ann$Name[ann$chr %in% c("chrX","chrY")])
gmset_flt <- gmset[autosomes,]
#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()## X11cairo
## 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 245993 2344 2566 2439 1909 3492
## NotSig 59114 786047 786358 786292 786913 785103
## Up 485551 2267 1734 1927 1836 2063
## Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13 ADOS
## Down 48678 2099 2939 1854 2557 0
## NotSig 728154 786567 785798 786943 786409 790658
## Up 13826 1992 1921 1861 1692 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] 16output <-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 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| cg12565057 | chr7 | 100482617 |
|
chr7:100484425-100484651 | N_Shore | SRRT;SRRT;SRRT;SRRT | Body;Body;Body;Body | 0.0826621 | 3.0460140 | 7.387285 | 0.0000040 | 0.9999951 | 4.581264 |
| cg00546248 | chr1 | 2412743 |
|
chr1:2411125-2411404 | S_Shore | PLCH2 | Body | 0.1034780 | 1.8014843 | 6.865500 | 0.0000089 | 0.9999951 | 3.850448 |
| cg15740041 | chr16 | 87781336 |
|
OpenSea | KLHDC4;KLHDC4;KLHDC4 | Body;Body;Body | 0.1317611 | 1.6008461 | 6.823044 | 0.0000095 | 0.9999951 | 3.789115 | |
| cg05116656 | chr5 | 94417892 |
|
OpenSea | MCTP1;MCTP1;MCTP1 | TSS1500;TSS1500;Body | 0.0915714 | -3.3750586 | 6.719336 | 0.0000112 | 0.9999951 | 3.638100 | |
| cg11632273 | chr22 | 20102580 |
|
chr22:20104374-20105729 | N_Shore | TRMT2A;TRMT2A;TRMT2A;RANBP1;MIR6816 | Body;Body;Body;TSS1500;TSS1500 | -0.1332566 | 2.2681145 | -6.338476 | 0.0000208 | 0.9999951 | 3.068838 |
| cg01656620 | chr12 | 124788698 |
|
chr12:124788559-124788884 | Island | FAM101A | 5’UTR | -0.1330939 | 3.0511986 | -6.203748 | 0.0000259 | 0.9999951 | 2.861912 |
| cg10881749 | chr1 | 234792390 |
|
OpenSea | 0.0804887 | 0.9589997 | 6.114550 | 0.0000301 | 0.9999951 | 2.723318 | |||
| cg13699963 | chr17 | 71548861 |
|
OpenSea | SDK2 | Body | -0.0663828 | 2.9801617 | -5.750106 | 0.0000558 | 0.9999951 | 2.143861 | |
| cg22656126 | chr17 | 1637206 |
|
chr17:1637053-1637492 | Island | WDR81;WDR81;WDR81;WDR81 | Body;Body;Body;Body | 0.1395930 | 5.7130309 | 5.698989 | 0.0000609 | 0.9999951 | 2.060902 |
| cg07162704 | chr6 | 91113687 |
|
OpenSea | 0.1009010 | 0.5704454 | 5.695702 | 0.0000613 | 0.9999951 | 2.055553 | |||
| cg13277464 | chr2 | 43393040 |
|
OpenSea | 0.1211560 | 2.3025832 | 5.606707 | 0.0000715 | 0.9999951 | 1.910100 | |||
| cg18731055 | chr13 | 49108131 |
|
chr13:49106411-49107463 | S_Shore | RCBTB2 | TSS1500 | 0.1016068 | 0.0882235 | 5.584670 | 0.0000743 | 0.9999951 | 1.873891 |
| cg03188919 | chr10 | 31372663 |
|
OpenSea | -0.1030369 | 2.0662356 | -5.526631 | 0.0000823 | 0.9999951 | 1.778167 | |||
| cg01128603 | chr11 | 65820391 |
|
chr11:65819367-65820224 | S_Shore | SF3B2 | Body | 0.0946490 | -3.2440188 | 5.486222 | 0.0000883 | 0.9999951 | 1.711213 |
| cg03177876 | chr4 | 71701931 |
|
chr4:71704664-71705736 | N_Shelf | GRSF1;GRSF1 | 5’UTR;Body | 0.0982965 | -1.2874739 | 5.481370 | 0.0000891 | 0.9999951 | 1.703158 |
| cg05791256 | chr9 | 96077420 |
|
chr9:96080113-96080330 | N_Shelf | WNK2;WNK2 | Body;Body | -0.0749658 | 3.6607111 | -5.477688 | 0.0000896 | 0.9999951 | 1.697042 |
| cg25723331 | chr11 | 44613980 |
|
OpenSea | CD82;CD82 | 5’UTR;5’UTR | 0.1123597 | 0.8048882 | 5.413638 | 0.0001004 | 0.9999951 | 1.590320 | |
| cg01016829 | chr16 | 46682692 |
|
OpenSea | 0.0860670 | 2.2697378 | 5.369826 | 0.0001084 | 0.9999951 | 1.516959 | |||
| cg23157290 | chr7 | 138602695 |
|
OpenSea | KIAA1549;KIAA1549 | Body;Body | 0.1583347 | 3.2345901 | 5.356147 | 0.0001111 | 0.9999951 | 1.493994 | |
| cg25652742 | chr2 | 240152334 |
|
chr2:240153476-240153792 | N_Shore | HDAC4 | Body | 0.1179861 | -1.9251649 | 5.354083 | 0.0001115 | 0.9999951 | 1.490528 |
| cg23061257 | chr1 | 54355272 |
|
chr1:54355354-54355656 | N_Shore | YIPF1 | 5’UTR | -0.0648082 | -3.3673834 | -5.351267 | 0.0001121 | 0.9999951 | 1.485795 |
| cg08008144 | chr8 | 1983127 |
|
OpenSea | -0.0927890 | 2.5786632 | -5.319054 | 0.0001187 | 0.9999951 | 1.431580 | |||
| cg23404427 | chr19 | 1653402 |
|
OpenSea | TCF3 | TSS1500 | -0.0903479 | -5.0294197 | -5.300923 | 0.0001226 | 0.9999951 | 1.400998 | |
| cg10619752 | chr2 | 233863457 |
|
OpenSea | NGEF | 5’UTR | -0.0820301 | 3.0526878 | -5.290482 | 0.0001249 | 0.9999951 | 1.383363 | |
| cg08225137 | chr3 | 54915807 |
|
OpenSea | CACNA2D3 | Body | -0.0561511 | 2.6468160 | -5.227687 | 0.0001397 | 0.9999951 | 1.276965 | |
| cg15856662 | chr22 | 38380501 |
|
chr22:38379093-38379964 | S_Shore | SOX10;SOX10;POLR2F;POLR2F | 5’UTR;1stExon;Body;Body | -0.0734815 | 2.4748310 | -5.180008 | 0.0001522 | 0.9999951 | 1.195788 |
| cg11355135 | chr2 | 108602860 |
|
chr2:108602824-108603467 | Island | SLC5A7 | TSS200 | 0.1104571 | -3.8226413 | 5.113981 | 0.0001714 | 0.9999951 | 1.082825 |
| cg07873041 | chr1 | 91301840 |
|
chr1:91300979-91301891 | Island | 0.0679299 | -5.4645582 | 5.090814 | 0.0001788 | 0.9999951 | 1.043041 | ||
| cg15165643 | chr8 | 145518600 |
|
chr8:145514719-145515959 | S_Shelf | HSF1 | Body | -0.0814095 | 3.0943937 | -5.080871 | 0.0001820 | 0.9999951 | 1.025942 |
| cg00183468 | chr3 | 142926178 |
|
OpenSea | 0.0609273 | 1.9139298 | 5.080236 | 0.0001822 | 0.9999951 | 1.024849 |
Other models:
motor skills impairment
diagnosis
SRS
inverse IQ
inverse langguage
Regression analysis - within twin pair on motor skills impairment
#motor:820,1.053055e-12
design_tw <- model.matrix(~ Twin_Group + motor )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))## (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down 222835 3100 3073 1578 1451 2457
## NotSig 109032 784572 785915 788025 788150 786947
## Up 458791 2986 1670 1055 1057 1254
## Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13 motor
## Down 152974 2956 3138 1373 1855 433
## NotSig 607875 784741 785529 788188 787768 789838
## Up 29809 2961 1991 1097 1035 387top <- 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] 53567output <-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.487554 | 2.7456613 | -24.58356 | 0 | 0.0000008 | 11.331647 | |||
| cg01593834 | chr5 | 145718401 |
|
chr5:145718289-145720095 | Island | POU4F3 | TSS200 | -3.047554 | -1.4836163 | -19.69127 | 0 | 0.0000080 | 10.529614 |
| cg08840551 | chr5 | 56718349 |
|
OpenSea | -2.630660 | -0.3473951 | -18.30316 | 0 | 0.0000140 | 10.208855 | |||
| cg09002832 | chr9 | 4300181 |
|
chr9:4297817-4300182 | Island | GLIS3 | TSS200 | 2.576547 | -4.7977912 | 17.43241 | 0 | 0.0000198 | 9.978451 |
| cg01059824 | chr15 | 69706539 |
|
chr15:69706516-69707352 | Island | KIF23;KIF23;KIF23 | TSS200;TSS200;TSS200 | 2.125246 | -4.7737263 | 16.53201 | 0 | 0.0000281 | 9.712490 |
| cg05098590 | chr10 | 111767379 |
|
chr10:111767087-111768355 | Island | ADD3;ADD3;ADD3 | 5’UTR;TSS1500;TSS1500 | 2.378836 | -5.5155179 | 16.45865 | 0 | 0.0000281 | 9.689455 |
| cg23199907 | chr13 | 33305966 |
|
OpenSea | PDS5B | Body | -3.696135 | 2.7389798 | -16.07868 | 0 | 0.0000312 | 9.566639 | |
| cg25259284 | chr1 | 52521812 |
|
chr1:52521780-52522296 | Island | TXNDC12;TXNDC12;BTF3L4;BTF3L4;BTF3L4;TXNDC12;TXNDC12 | 1stExon;5’UTR;TSS200;TSS200;TSS200;Body;Body | 2.467314 | -4.1955426 | 15.94548 | 0 | 0.0000312 | 9.522146 |
| cg09804439 | chr17 | 40540457 |
|
chr17:40539837-40540775 | Island | STAT3;STAT3;STAT3;STAT3;STAT3 | 5’UTR;1stExon;1stExon;5’UTR;TSS200 | 2.103803 | -5.6192796 | 15.82500 | 0 | 0.0000312 | 9.481235 |
| cg04055615 | chr6 | 31697723 |
|
chr6:31695894-31698245 | Island | DDAH2 | 5’UTR | 2.215811 | -4.9498235 | 15.15245 | 0 | 0.0000449 | 9.240656 |
| cg07241733 | chr8 | 68255213 |
|
chr8:68255129-68256274 | Island | ARFGEF1 | Body | 2.256471 | -3.1368420 | 15.14960 | 0 | 0.0000449 | 9.239591 |
| cg19571293 | chr19 | 1028488 |
|
chr19:1028319-1028641 | Island | CNN2;CNN2 | Body;Body | 2.420525 | -1.9882125 | 15.03487 | 0 | 0.0000455 | 9.196366 |
| cg15592945 | chr1 | 226736711 |
|
chr1:226736354-226737412 | Island | C1orf95 | 1stExon | 1.738758 | -5.7015620 | 13.92959 | 0 | 0.0000939 | 8.744156 |
| cg17121412 | chr15 | 22833400 |
|
chr15:22833282-22833911 | Island | TUBGCP5;TUBGCP5;TUBGCP5;TUBGCP5 | 1stExon;5’UTR;5’UTR;1stExon | -1.746412 | -4.7364030 | -13.85471 | 0 | 0.0000939 | 8.711002 |
| cg07762054 | chr11 | 75273187 |
|
chr11:75272950-75273616 | Island | SERPINH1;SERPINH1 | 5’UTR;1stExon | 1.789843 | -5.0227036 | 13.78260 | 0 | 0.0000939 | 8.678756 |
| cg14883700 | chr2 | 173096328 |
|
chr2:173099482-173100098 | N_Shelf | -1.964926 | 3.6945378 | -13.77481 | 0 | 0.0000939 | 8.675253 | ||
| cg07220356 | chr6 | 33263442 |
|
chr6:33266302-33267582 | N_Shelf | RGL2;RGL2 | Body;Body | -1.844915 | 3.5902570 | -13.75131 | 0 | 0.0000939 | 8.664668 |
| cg22157239 | chr1 | 18972096 |
|
chr1:18971729-18972097 | Island | PAX7;PAX7;PAX7 | Body;Body;Body | 2.008880 | -4.8748956 | 13.71440 | 0 | 0.0000939 | 8.647966 |
| cg23422659 | chr17 | 44929215 |
|
chr17:44928287-44929690 | Island | WNT9B | Body | -1.862788 | -5.2099662 | -13.70637 | 0 | 0.0000939 | 8.644320 |
| cg04114460 | chr18 | 29239388 |
|
OpenSea | B4GALT6 | Body | -2.312635 | 2.8380218 | -13.65741 | 0 | 0.0000939 | 8.622013 | |
| cg09669099 | chr16 | 69419557 |
|
chr16:69419316-69420086 | Island | TERF2 | 1stExon | 1.608697 | -5.1869754 | 13.44182 | 0 | 0.0001095 | 8.521973 |
| cg01957330 | chr7 | 26904470 |
|
chr7:26903869-26904753 | Island | SKAP2 | TSS200 | 1.999471 | -4.4333864 | 13.35145 | 0 | 0.0001140 | 8.479147 |
| cg04426862 | chr14 | 37131008 |
|
chr14:37131181-37132785 | N_Shore | PAX9 | 5’UTR | 1.655081 | -4.4928784 | 13.20356 | 0 | 0.0001256 | 8.407898 |
| cg26763618 | chr16 | 14726469 |
|
chr16:14726516-14727366 | N_Shore | BFAR | TSS200 | 1.907022 | -4.4007952 | 13.04177 | 0 | 0.0001408 | 8.328258 |
| cg26516701 | chr12 | 53845683 |
|
chr12:53845668-53846606 | Island | PCBP2;PCBP2;PCBP2;PCBP2;PCBP2;PCBP2;PCBP2 | TSS1500;TSS1500;TSS1500;TSS1500;TSS1500;TSS1500;TSS1500 | 1.502112 | -4.5229391 | 12.94985 | 0 | 0.0001479 | 8.282206 |
| cg13523875 | chr7 | 117823728 |
|
chr7:117823984-117824485 | N_Shore | LSM8 | TSS1500 | 1.508651 | -5.5317694 | 12.88476 | 0 | 0.0001516 | 8.249235 |
| cg21434530 | chr4 | 4870486 |
|
chr4:4868440-4869173 | S_Shore | 1.814807 | -4.8330389 | 12.78776 | 0 | 0.0001607 | 8.199539 | ||
| cg20221020 | chr11 | 67209257 |
|
OpenSea | CORO1B;CORO1B | Body;Body | -1.974616 | 4.9962484 | -12.63877 | 0 | 0.0001797 | 8.121881 | |
| cg02005873 | chr6 | 10385839 |
|
chr6:10383525-10384114 | S_Shore | 1.833420 | -1.2450526 | 12.55614 | 0 | 0.0001885 | 8.078111 | ||
| cg20729301 | chr2 | 233368115 |
|
chr2:233367780-233368577 | Island | 1.782130 | -6.2407714 | 12.46508 | 0 | 0.0001998 | 8.029276 |
saveRDS(design_tw, "gu_design_motor.rds")Regression analysis - within twin pair on diagnosis
#diagnosis:0,3.080605e-05
design_tw <- model.matrix(~ Twin_Group + diagnosis )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))## (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down 246805 2371 2834 2476 1954 4329
## NotSig 54218 786027 785932 786232 786824 783943
## Up 489635 2260 1892 1950 1880 2386
## Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13 diagnosis
## Down 45182 2607 2861 1993 2546 0
## NotSig 733235 785493 785921 786617 786426 790658
## Up 12241 2558 1876 2048 1686 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] 16978output <-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 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| cg14507658 | chr21 | 46408142 |
|
chr21:46407925-46408143 | Island | 0.5671208 | 3.3489041 | 6.094954 | 0.0000308 | 0.9905792 | -1.590186 | ||
| cg05281810 | chr18 | 77542644 |
|
chr18:77542920-77543323 | N_Shore | -0.3448119 | 3.1050953 | -6.089657 | 0.0000311 | 0.9905792 | -1.592175 | ||
| cg11355135 | chr2 | 108602860 |
|
chr2:108602824-108603467 | Island | SLC5A7 | TSS200 | 0.6328982 | -3.8226413 | 5.874010 | 0.0000448 | 0.9905792 | -1.675536 |
| cg08008144 | chr8 | 1983127 |
|
OpenSea | -0.5206780 | 2.5786632 | -5.716288 | 0.0000587 | 0.9905792 | -1.739575 | |||
| cg02344606 | chr7 | 112062297 |
|
OpenSea | IFRD1 | TSS1500 | 0.4806032 | -2.0898626 | 5.667407 | 0.0000639 | 0.9905792 | -1.759966 | |
| cg18731055 | chr13 | 49108131 |
|
chr13:49106411-49107463 | S_Shore | RCBTB2 | TSS1500 | 0.5571933 | 0.0882235 | 5.572878 | 0.0000753 | 0.9905792 | -1.800145 |
| cg00633277 | chr12 | 7250385 |
|
OpenSea | C1RL;C1RL;C1RL | Body;Body;Body | -0.4834382 | 0.7795169 | -5.473622 | 0.0000896 | 0.9905792 | -1.843412 | |
| cg15740041 | chr16 | 87781336 |
|
OpenSea | KLHDC4;KLHDC4;KLHDC4 | Body;Body;Body | 0.6813155 | 1.6008461 | 5.469788 | 0.0000902 | 0.9905792 | -1.845106 | |
| cg24267485 | chr18 | 74766743 |
|
chr18:74770351-74770590 | N_Shelf | MBP;MBP | Body;Body | 0.5571721 | -2.8069857 | 5.386618 | 0.0001045 | 0.9905792 | -1.882264 |
| cg14529256 | chr1 | 245843607 |
|
OpenSea | KIF26B | Body | 0.5967061 | 2.3027143 | 5.337372 | 0.0001141 | 0.9905792 | -1.904644 | |
| cg25272125 | chr1 | 31538096 |
|
chr1:31538412-31538938 | N_Shore | PUM1;PUM1 | 5’UTR;5’UTR | 0.4303036 | -5.3122575 | 5.324817 | 0.0001166 | 0.9905792 | -1.910395 |
| cg16239628 | chr13 | 46626330 |
|
chr13:46626454-46627019 | N_Shore | ZC3H13 | 5’UTR | -0.5717242 | -4.6148948 | -5.229213 | 0.0001384 | 0.9905792 | -1.954801 |
| cg11712265 | chr19 | 7968106 |
|
chr19:7968325-7969421 | N_Shore | MAP2K7;MAP2K7;MAP2K7 | TSS1500;TSS1500;TSS1500 | 0.3435839 | -3.7499297 | 5.220197 | 0.0001406 | 0.9905792 | -1.959045 |
| cg27149150 | chr1 | 7449851 |
|
OpenSea | CAMTA1 | Body | -0.4871695 | -0.3731033 | -5.217092 | 0.0001414 | 0.9905792 | -1.960509 | |
| cg23404427 | chr19 | 1653402 |
|
OpenSea | TCF3 | TSS1500 | -0.4921666 | -5.0294197 | -5.184932 | 0.0001499 | 0.9905792 | -1.975739 | |
| cg04900486 | chr1 | 151320117 |
|
chr1:151319326-151319545 | S_Shore | RFX5;RFX5 | TSS1500;TSS1500 | 0.4810737 | -4.3365457 | 5.163718 | 0.0001557 | 0.9905792 | -1.985854 |
| cg08277291 | chr6 | 32037426 |
|
OpenSea | TNXB | Body | -0.4216451 | 3.0373254 | -5.158534 | 0.0001572 | 0.9905792 | -1.988333 | |
| cg17904575 | chr14 | 102290370 |
|
OpenSea | PPP2R5C;PPP2R5C;PPP2R5C;PPP2R5C;PPP2R5C | Body;Body;Body;Body;Body | -0.3901626 | 1.0410297 | -5.134116 | 0.0001642 | 0.9905792 | -2.000058 | |
| cg05116656 | chr5 | 94417892 |
|
OpenSea | MCTP1;MCTP1;MCTP1 | TSS1500;TSS1500;Body | 0.4638695 | -3.3750586 | 5.081264 | 0.0001807 | 0.9905792 | -2.025685 | |
| cg08471713 | chr17 | 41738893 |
|
OpenSea | MEOX1;MEOX1;MEOX1 | 1stExon;1stExon;5’UTR | -0.5408241 | -2.6803562 | -5.078389 | 0.0001817 | 0.9905792 | -2.027089 | |
| cg01497262 | chr2 | 43635730 |
|
OpenSea | THADA;THADA;THADA | Body;Body;Body | 0.5232107 | 0.3214082 | 5.007192 | 0.0002068 | 0.9905792 | -2.062179 | |
| cg06595448 | chr7 | 50727315 |
|
OpenSea | GRB10;GRB10;GRB10;GRB10 | Body;Body;Body;Body | 0.3863605 | 2.2935124 | 4.959104 | 0.0002258 | 0.9905792 | -2.086235 | |
| cg26684767 | chr2 | 231000057 |
|
OpenSea | 0.4316435 | -1.9244246 | 4.932767 | 0.0002370 | 0.9905792 | -2.099533 | |||
| cg25481455 | chr6 | 22147551 |
|
OpenSea | NBAT1;CASC15 | TSS200;Body | 0.4057296 | -4.6757321 | 4.928824 | 0.0002387 | 0.9905792 | -2.101531 | |
| cg25022405 | chr22 | 26137668 |
|
OpenSea | MYO18B | TSS1500 | 0.3427957 | 3.6944563 | 4.919226 | 0.0002429 | 0.9905792 | -2.106403 | |
| cg06753153 | chr14 | 65186982 |
|
OpenSea | PLEKHG3 | 5’UTR | -0.3732891 | 3.1661977 | -4.914021 | 0.0002453 | 0.9905792 | -2.109050 | |
| cg11250279 | chr13 | 24463401 |
|
chr13:24463119-24463677 | Island | C1QTNF9B-AS1;MIPEP;C1QTNF9B-AS1 | 5’UTR;1stExon;TSS200 | 0.3386238 | -4.5577692 | 4.896740 | 0.0002532 | 0.9905792 | -2.117864 |
| cg21175642 | chr3 | 48701770 |
|
chr3:48698335-48701667 | S_Shore | CELSR3 | TSS1500 | 0.5201876 | -3.2165700 | 4.887731 | 0.0002574 | 0.9905792 | -2.122473 |
| cg20492597 | chr5 | 161274081 |
|
OpenSea | GABRA1;GABRA1;GABRA1;GABRA1 | TSS1500;TSS1500;TSS200;TSS1500 | 0.5847567 | -2.6606980 | 4.885151 | 0.0002586 | 0.9905792 | -2.123795 | |
| cg00010078 | chr2 | 109967172 |
|
OpenSea | SH3RF3 | Body | 0.4001143 | -0.0677274 | 4.882349 | 0.0002600 | 0.9905792 | -2.125231 |
saveRDS(design_tw, "gu_design_diagnosis.rds")Regression analysis - within twin pair on SRS
#SRS:0,1.895198e-07
design_tw <- model.matrix(~ Twin_Group + SRS )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))## (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down 216305 2079 3241 1609 1985 54923
## NotSig 130496 786758 785600 787930 786968 727301
## Up 443857 1821 1817 1119 1705 8434
## Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13 SRS
## Down 203038 2981 3165 1345 1041 0
## NotSig 553849 784700 785433 788210 789043 790658
## Up 33771 2977 2060 1103 574 0top <- topTable(fit2_tw,coef="SRS",num=Inf, sort.by = "P")
head(top)## logFC AveExpr t P.Value adj.P.Val B
## cg19720702 -0.03841315 3.071925 -9.525594 1.895198e-07 0.1498453 7.225728
## cg23884297 0.02680846 2.572387 7.478187 3.217331e-06 0.4586088 4.265063
## cg08281515 -0.04465879 -3.213167 -7.424466 3.488795e-06 0.4586088 4.180428
## cg15865511 0.02908008 2.876906 7.373374 3.769491e-06 0.4586088 4.099582
## cg03767531 -0.03515160 -4.926828 -6.952022 7.228624e-06 0.4586088 3.419673
## cg13534438 -0.02551836 2.702726 -6.877745 8.127541e-06 0.4586088 3.297352nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)## [1] 85217output <-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.0384131 | 3.0719249 | -9.525594 | 2.00e-07 | 0.1498453 | 7.225728 | |||
| cg23884297 | chr20 | 30124811 |
|
OpenSea | HM13;HM13;HM13;HM13 | Body;Body;Body;Body | 0.0268085 | 2.5723870 | 7.478187 | 3.20e-06 | 0.4586088 | 4.265063 | |
| cg08281515 | chr11 | 111956885 |
|
chr11:111957352-111957681 | N_Shore | TIMM8B;SDHD;TIMM8B | Body;TSS1500;Body | -0.0446588 | -3.2131669 | -7.424466 | 3.50e-06 | 0.4586088 | 4.180428 |
| cg15865511 | chr2 | 43252424 |
|
OpenSea | 0.0290801 | 2.8769063 | 7.373374 | 3.80e-06 | 0.4586088 | 4.099582 | |||
| cg03767531 | chr2 | 182545572 |
|
chr2:182547873-182549177 | N_Shelf | NEUROD1 | TSS200 | -0.0351516 | -4.9268280 | -6.952022 | 7.20e-06 | 0.4586088 | 3.419673 |
| cg13534438 | chr3 | 197229854 |
|
chr3:197225631-197226150 | S_Shelf | -0.0255184 | 2.7027261 | -6.877745 | 8.10e-06 | 0.4586088 | 3.297352 | ||
| cg08748969 | chr12 | 69327779 |
|
chr12:69327021-69327532 | S_Shore | CPM;CPM;CPM | TSS1500;TSS1500;5’UTR | 0.0253835 | -1.2268838 | 6.696612 | 1.09e-05 | 0.4586088 | 2.995932 |
| cg19868495 | chr13 | 61995249 |
|
OpenSea | 0.0325808 | 0.8882802 | 6.667498 | 1.14e-05 | 0.4586088 | 2.947068 | |||
| cg13206794 | chr2 | 242010703 |
|
OpenSea | SNED1 | Body | -0.0254861 | 2.0212841 | -6.654190 | 1.16e-05 | 0.4586088 | 2.924695 | |
| cg03013237 | chr16 | 1802436 |
|
OpenSea | MAPK8IP3;MAPK8IP3 | Body;Body | -0.0339334 | 2.3972206 | -6.523127 | 1.44e-05 | 0.4586088 | 2.703060 | |
| cg15616400 | chr1 | 16062894 |
|
OpenSea | SLC25A34;SLC25A34 | 5’UTR;1stExon | -0.0347249 | 2.2495978 | -6.514356 | 1.46e-05 | 0.4586088 | 2.688143 | |
| cg06874119 | chr5 | 1892431 |
|
OpenSea | CTD-2194D22.4 | Body | 0.0247034 | -1.0871226 | 6.487037 | 1.52e-05 | 0.4586088 | 2.641619 | |
| cg17916960 | chr15 | 79447300 |
|
OpenSea | 0.0207799 | -2.1616055 | 6.463602 | 1.58e-05 | 0.4586088 | 2.601628 | |||
| cg01201782 | chr16 | 1652449 |
|
OpenSea | IFT140 | Body | -0.0237944 | 3.6792856 | -6.460084 | 1.59e-05 | 0.4586088 | 2.595618 | |
| cg17497176 | chr7 | 150763829 |
|
OpenSea | SLC4A2 | Body | -0.0322948 | 2.6257672 | -6.381261 | 1.81e-05 | 0.4586088 | 2.460518 | |
| cg21401292 | chr6 | 134570669 |
|
OpenSea | SGK1 | Body | 0.0263213 | 1.8937236 | 6.278702 | 2.15e-05 | 0.4586088 | 2.283470 | |
| cg22593006 | chr6 | 170603682 |
|
chr6:170604411-170606479 | N_Shore | FAM120B;FAM120B | TSS1500;Body | 0.0302294 | -0.7067929 | 6.253230 | 2.24e-05 | 0.4586088 | 2.239277 |
| cg13385114 | chr1 | 9656493 |
|
OpenSea | TMEM201;TMEM201 | Body;Body | -0.0318501 | 2.3856163 | -6.245681 | 2.27e-05 | 0.4586088 | 2.226162 | |
| cg15697902 | chr19 | 14676213 |
|
chr19:14676201-14676542 | Island | TECR | Body | 0.0383106 | 3.4063584 | 6.206997 | 2.42e-05 | 0.4586088 | 2.158836 |
| cg16134349 | chr10 | 61900552 |
|
OpenSea | ANK3;ANK3;ANK3;ANK3;ANK3 | 1stExon;5’UTR;Body;Body;Body | -0.0194575 | 2.8615804 | -6.178204 | 2.54e-05 | 0.4586088 | 2.108592 | |
| cg09451572 | chr22 | 46760120 |
|
chr22:46759899-46760121 | Island | CELSR1 | Body | -0.0263675 | 2.3694795 | -6.129153 | 2.76e-05 | 0.4586088 | 2.022739 |
| cg14266530 | chr13 | 52511468 |
|
OpenSea | ATP7B;ATP7B;ATP7B | Body;Body;Body | -0.0285055 | 3.7038678 | -6.128550 | 2.76e-05 | 0.4586088 | 2.021682 | |
| cg11920999 | chr22 | 37417225 |
|
chr22:37414319-37416111 | S_Shore | MPST;MPST;TST;MPST | 5’UTR;5’UTR;TSS1500;Body | 0.0254667 | 3.2107679 | 6.122789 | 2.79e-05 | 0.4586088 | 2.011577 |
| cg13709529 | chr6 | 54001895 |
|
OpenSea | C6orf142 | Body | 0.0351276 | 2.9963215 | 6.106824 | 2.87e-05 | 0.4586088 | 1.983550 | |
| cg01038640 | chr5 | 157097845 |
|
chr5:157098263-157099041 | N_Shore | C5orf52 | TSS1500 | -0.0344627 | -3.1684180 | -6.100173 | 2.90e-05 | 0.4586088 | 1.971864 |
| cg23657393 | chr20 | 42790697 |
|
chr20:42788216-42789168 | S_Shore | JPH2 | Body | -0.0280709 | 2.5652689 | -6.097772 | 2.91e-05 | 0.4586088 | 1.967644 |
| cg09487134 | chr10 | 35195634 |
|
OpenSea | -0.0194752 | 0.2801082 | -6.079980 | 3.00e-05 | 0.4586088 | 1.936345 | |||
| cg00647512 | chr8 | 42669499 |
|
OpenSea | -0.0351548 | 2.2233507 | -6.071817 | 3.04e-05 | 0.4586088 | 1.921973 | |||
| cg18764121 | chr15 | 25017456 |
|
chr15:25018174-25018533 | N_Shore | -0.0261292 | 0.1184601 | -6.052347 | 3.14e-05 | 0.4586088 | 1.887653 | ||
| cg17594278 | chr12 | 112192763 |
|
OpenSea | ACAD10;ACAD10 | Body;Body | -0.0281927 | 2.5565854 | -6.043254 | 3.19e-05 | 0.4586088 | 1.871607 |
saveRDS(design_tw, "gu_design_srs.rds")Regression analysis - within twin pair on inverse IQ
#iiq:0,1.326967e-05
design_tw <- model.matrix(~ Twin_Group + iiq )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))## (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down 247642 2515 3381 2813 1956 14775
## NotSig 47504 786029 785507 785809 786712 771811
## Up 495512 2114 1770 2036 1990 4072
## Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13 iiq
## Down 182628 2926 2792 2019 3029 0
## NotSig 575220 784703 786190 786679 785893 790658
## Up 32810 3029 1676 1960 1736 0top <- topTable(fit2_tw,coef="iiq",num=Inf, sort.by = "P")
head(top)## logFC AveExpr t P.Value adj.P.Val B
## cg13470919 0.4700907 -3.017373 6.591035 1.326967e-05 0.3977294 -0.05845352
## cg09714307 -0.6303321 3.218485 -6.556696 1.402703e-05 0.3977294 -0.07677677
## cg08210568 0.4937255 -4.368887 6.459055 1.643964e-05 0.3977294 -0.12970089
## cg15934248 0.4688290 -4.799787 6.399469 1.812315e-05 0.3977294 -0.16260397
## cg03922146 0.4240616 2.534101 6.378132 1.876925e-05 0.3977294 -0.17449911
## cg26869412 -0.5128173 5.366690 -6.229826 2.398524e-05 0.3977294 -0.25885486nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)## [1] 99359output <-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.4700907 | -3.0173732 | 6.591035 | 1.33e-05 | 0.3977294 | -0.0584535 |
| cg09714307 | chr5 | 10483440 |
|
OpenSea | -0.6303321 | 3.2184848 | -6.556697 | 1.40e-05 | 0.3977294 | -0.0767768 | |||
| cg08210568 | chr3 | 15141216 |
|
chr3:15140035-15140890 | S_Shore | ZFYVE20 | TSS1500 | 0.4937255 | -4.3688871 | 6.459054 | 1.64e-05 | 0.3977294 | -0.1297009 |
| cg15934248 | chr2 | 86332749 |
|
chr2:86332659-86333650 | Island | POLR1A;PTCD3 | Body;TSS1500 | 0.4688290 | -4.7997874 | 6.399469 | 1.81e-05 | 0.3977294 | -0.1626040 |
| cg03922146 | chr6 | 30656499 |
|
chr6:30654392-30654934 | S_Shore | KIAA1949;NRM;KIAA1949 | TSS1500;Body;TSS1500 | 0.4240616 | 2.5341006 | 6.378132 | 1.88e-05 | 0.3977294 | -0.1744991 |
| cg26869412 | chr10 | 3023801 |
|
OpenSea | -0.5128173 | 5.3666898 | -6.229826 | 2.40e-05 | 0.3977294 | -0.2588549 | |||
| cg07241312 | chr16 | 87051295 |
|
OpenSea | 0.5246257 | -1.2528635 | 6.078011 | 3.09e-05 | 0.3977294 | -0.3483084 | |||
| cg05845592 | chr16 | 28634766 |
|
chr16:28634440-28635031 | Island | SULT1A1;SULT1A1 | 5’UTR;1stExon | 0.6166146 | -5.6410886 | 5.938857 | 3.92e-05 | 0.3977294 | -0.4331377 |
| cg15957908 | chr6 | 158923478 |
|
OpenSea | TULP4;TULP4 | Body;Body | -0.3848645 | 3.6984208 | -5.887559 | 4.27e-05 | 0.3977294 | -0.4651084 | |
| cg16336665 | chr17 | 73028495 |
|
chr17:73030676-73031160 | N_Shelf | KCTD2 | TSS200 | 0.7384996 | -2.2320450 | 5.869774 | 4.41e-05 | 0.3977294 | -0.4762820 |
| cg07943548 | chr15 | 86304357 |
|
chr15:86301949-86302213 | S_Shelf | KLHL25 | 3’UTR | -0.4935287 | 4.2351942 | -5.824888 | 4.76e-05 | 0.3977294 | -0.5046867 |
| cg13709580 | chr5 | 10668678 |
|
OpenSea | -0.3628510 | 2.2360219 | -5.820586 | 4.79e-05 | 0.3977294 | -0.5074243 | |||
| cg23905944 | chr20 | 62571879 |
|
chr20:62571738-62572556 | Island | UCKL1 | Body | -0.5878267 | 3.4256349 | -5.807980 | 4.90e-05 | 0.3977294 | -0.5154632 |
| cg09375488 | chr7 | 95115026 |
|
OpenSea | ASB4;ASB4 | TSS1500;TSS1500 | -0.4619154 | 0.2983563 | -5.806454 | 4.91e-05 | 0.3977294 | -0.5164376 | |
| cg08848154 | chr10 | 6278954 |
|
OpenSea | -0.5018118 | 5.9889411 | -5.805757 | 4.92e-05 | 0.3977294 | -0.5168831 | |||
| cg21419752 | chr3 | 139514052 |
|
OpenSea | -0.4099476 | 1.0313638 | -5.766283 | 5.26e-05 | 0.3977294 | -0.5422187 | |||
| cg26362488 | chr12 | 14924051 |
|
chr12:14922877-14923974 | S_Shore | HIST4H4;HIST4H4 | 1stExon;5’UTR | 0.8214279 | -4.6320444 | 5.765836 | 5.27e-05 | 0.3977294 | -0.5425070 |
| cg07690326 | chr8 | 145019032 |
|
chr8:145018815-145019214 | Island | PLEC1;PLEC1;PLEC1;PLEC1;PLEC1;PLEC1 | Body;Body;Body;TSS1500;Body;TSS200 | -0.5396723 | 5.3537324 | -5.753420 | 5.38e-05 | 0.3977294 | -0.5505243 |
| cg25849262 | chr3 | 19189070 |
|
chr3:19189688-19190100 | N_Shore | KCNH8 | TSS1500 | 0.4565090 | -4.5611537 | 5.750801 | 5.41e-05 | 0.3977294 | -0.5522186 |
| cg14570121 | chr19 | 45874153 |
|
chr19:45873339-45874033 | S_Shore | ERCC2;ERCC2 | TSS1500;TSS1500 | 0.3972347 | -4.7541665 | 5.725689 | 5.65e-05 | 0.3977294 | -0.5685141 |
| cg23644589 | chr1 | 228404256 |
|
chr1:228404148-228404397 | Island | OBSCN;OBSCN | Body;Body | -0.6047790 | 3.8404119 | -5.712356 | 5.78e-05 | 0.3977294 | -0.5772043 |
| cg27325462 | chr11 | 129903118 |
|
chr11:129902811-129903119 | Island | -0.6266766 | 4.3295504 | -5.698560 | 5.92e-05 | 0.3977294 | -0.5862246 | ||
| cg07019009 | chr15 | 102265803 |
|
chr15:102264097-102264750 | S_Shore | TARSL2 | TSS1500 | -0.5618290 | 1.0241536 | -5.697717 | 5.93e-05 | 0.3977294 | -0.5867766 |
| cg03873049 | chr5 | 149625194 |
|
OpenSea | CAMK2A;CAMK2A | Body;Body | -0.3784796 | 3.3104549 | -5.692352 | 5.98e-05 | 0.3977294 | -0.5902926 | |
| cg04465078 | chr16 | 15766942 |
|
OpenSea | NDE1;NDE1 | Body;Body | 0.9780949 | -2.3163770 | 5.691853 | 5.99e-05 | 0.3977294 | -0.5906199 | |
| cg02025435 | chr16 | 89986317 |
|
chr16:89985913-89986356 | Island | MC1R | 1stExon | -0.5241210 | 5.4439103 | -5.686550 | 6.04e-05 | 0.3977294 | -0.5941000 |
| cg19686637 | chr15 | 85427743 |
|
OpenSea | SLC28A1;SLC28A1 | TSS200;TSS200 | -0.3449407 | 3.2605496 | -5.685110 | 6.06e-05 | 0.3977294 | -0.5950455 | |
| cg19166616 | chr20 | 62259876 |
|
chr20:62257589-62259476 | S_Shore | GMEB2 | TSS1500 | -0.4381218 | -3.4166360 | -5.684392 | 6.07e-05 | 0.3977294 | -0.5955177 |
| cg21645826 | chr19 | 19313533 |
|
chr19:19314137-19314416 | N_Shore | NR2C2AP | Body | -0.4032165 | -2.4565009 | -5.673476 | 6.18e-05 | 0.3977294 | -0.6026977 |
| cg21353154 | chr13 | 99381820 |
|
OpenSea | SLC15A1 | Body | -0.4783244 | -0.3243186 | -5.654775 | 6.39e-05 | 0.3977294 | -0.6150408 |
saveRDS(design_tw, "gu_design_iiq.rds")Regression analysis - within twin pair on inverse language skills
#ilanguage:0,9.149182e-08
design_tw <- model.matrix(~ Twin_Group + ilanguage )
fit_tw <- lmFit(mvals, design_tw)
fit2_tw <- eBayes(fit_tw)
summary(decideTests(fit2_tw))## (Intercept) Twin_Group2 Twin_Group3 Twin_Group4 Twin_Group5 Twin_Group6
## Down 247342 1826 2269 2268 1972 5818
## NotSig 52060 787197 786886 786664 786825 782068
## Up 491256 1635 1503 1726 1861 2772
## Twin_Group7 Twin_Group8 Twin_Group9 Twin_Group12 Twin_Group13 ilanguage
## Down 138499 2740 2456 1792 2752 0
## NotSig 624750 785140 786604 787138 786170 790658
## Up 27409 2778 1598 1728 1736 0top <- topTable(fit2_tw,coef="ilanguage",num=Inf, sort.by = "P")
head(top)## logFC AveExpr t P.Value adj.P.Val B
## cg22165507 -0.7964738 -4.457186 -10.209898 9.149182e-08 0.07233874 -0.7254298
## cg13435526 -0.5169839 3.227245 -7.450780 3.624389e-06 0.99022496 -1.2323101
## cg25679269 0.9219619 -4.314010 6.921194 8.134347e-06 0.99022496 -1.3773828
## cg22533406 -0.6174444 4.882065 -6.864434 8.890001e-06 0.99022496 -1.3941709
## cg09476130 -0.4752016 -3.483859 -6.629915 1.289124e-05 0.99022496 -1.4662838
## cg18698884 0.7123973 -3.397628 6.556641 1.450054e-05 0.99022496 -1.4897535nsig <- sum(top$adj.P.Val < 0.05)
sum(top$P.Value< 0.05)## [1] 27881output <-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.7964738 | -4.457186 | -10.209898 | 1.00e-07 | 0.0723387 | -0.7254298 |
| cg13435526 | chr8 | 22472050 |
|
OpenSea | CCAR2;CCAR2 | Body;Body | -0.5169839 | 3.227245 | -7.450780 | 3.60e-06 | 0.9902250 | -1.2323101 | |
| cg25679269 | chr11 | 118747550 |
|
OpenSea | 0.9219619 | -4.314011 | 6.921194 | 8.10e-06 | 0.9902250 | -1.3773828 | |||
| cg22533406 | chr19 | 49106849 |
|
chr19:49106882-49107178 | N_Shore | FAM83E | Body | -0.6174444 | 4.882065 | -6.864434 | 8.90e-06 | 0.9902250 | -1.3941709 |
| cg09476130 | chr1 | 159870086 |
|
chr1:159869901-159870143 | Island | CCDC19 | TSS200 | -0.4752016 | -3.483859 | -6.629915 | 1.29e-05 | 0.9902250 | -1.4662838 |
| cg18698884 | chr11 | 16632724 |
|
chr11:16632508-16632725 | Island | 0.7123973 | -3.397628 | 6.556641 | 1.45e-05 | 0.9902250 | -1.4897535 | ||
| cg10021749 | chr7 | 2557358 |
|
chr7:2558371-2559967 | N_Shore | LFNG;LFNG | Body;TSS200 | -0.5742813 | 2.799210 | -6.470328 | 1.67e-05 | 0.9902250 | -1.5179938 |
| cg12843467 | chr15 | 99190301 |
|
chr15:99190446-99194559 | N_Shore | IRAIN;IGF1R;IGF1R | Body;TSS1500;TSS1500 | -0.4936995 | -4.971481 | -6.364014 | 1.98e-05 | 0.9902250 | -1.5536832 |
| cg13938969 | chr15 | 29409443 |
|
chr15:29407695-29408229 | S_Shore | APBA2;APBA2 | 3’UTR;3’UTR | -0.5900858 | 2.762155 | -6.346409 | 2.04e-05 | 0.9902250 | -1.5596911 |
| cg23404435 | chr7 | 1588396 |
|
chr7:1590387-1591020 | N_Shore | TMEM184A | Body | -0.5840304 | 2.412321 | -6.306655 | 2.18e-05 | 0.9902250 | -1.5733623 |
| cg13724550 | chr5 | 142001033 |
|
OpenSea | FGF1;FGF1;FGF1;FGF1;FGF1;FGF1;FGF1 | 5’UTR;TSS200;5’UTR;5’UTR;5’UTR;Body;Body | -0.5416578 | 3.493387 | -6.289624 | 2.24e-05 | 0.9902250 | -1.5792638 | |
| cg20781140 | chr2 | 175352216 |
|
chr2:175351364-175352154 | S_Shore | GPR155;GPR155 | TSS1500;TSS1500 | 0.5757550 | -4.026189 | 6.268960 | 2.32e-05 | 0.9902250 | -1.5864597 |
| cg19107578 | chr5 | 493262 |
|
chr5:494803-495306 | N_Shore | SLC9A3 | Body | 0.5244209 | 3.194001 | 6.174053 | 2.71e-05 | 0.9902250 | -1.6200247 |
| cg03547220 | chr19 | 7069348 |
|
chr19:7069396-7069921 | N_Shore | ZNF557;ZNF557;ZNF557 | TSS200;TSS200;TSS200 | 0.5614201 | -2.007877 | 6.150281 | 2.82e-05 | 0.9902250 | -1.6285656 |
| cg01369908 | chr11 | 34375707 |
|
chr11:34378229-34379993 | N_Shelf | ABTB2 | Body | -0.4638439 | 3.189187 | -6.143036 | 2.86e-05 | 0.9902250 | -1.6311793 |
| cg01231620 | chr14 | 23298884 |
|
chr14:23298984-23299313 | N_Shore | MRPL52;MRPL52;MRPL52;MRPL52;MRPL52;MRPL52 | TSS1500;TSS1500;TSS1500;TSS1500;TSS1500;TSS1500 | 0.5300644 | -4.348508 | 6.073912 | 3.21e-05 | 0.9902250 | -1.6563723 |
| cg11913416 | chr1 | 13909012 |
|
chr1:13910137-13910868 | N_Shore | PDPN;PDPN | TSS1500;TSS1500 | 0.4995019 | -0.696540 | 6.009515 | 3.57e-05 | 0.9902250 | -1.6802627 |
| cg02882857 | chr5 | 149381215 |
|
chr5:149379964-149380802 | S_Shore | HMGXB3;TIGD6;TIGD6 | 5’UTR;TSS1500;TSS1500 | 0.6426011 | -4.189717 | 5.997624 | 3.65e-05 | 0.9902250 | -1.6847187 |
| cg09219421 | chr19 | 13779879 |
|
OpenSea | -0.4059665 | 2.954756 | -5.922746 | 4.14e-05 | 0.9902250 | -1.7131055 | |||
| cg27594834 | chr19 | 11708137 |
|
chr19:11708277-11708683 | N_Shore | ZNF627 | TSS200 | -0.5832143 | -4.512216 | -5.851430 | 4.67e-05 | 0.9902250 | -1.7406710 |
| cg09459188 | chr3 | 158520176 |
|
chr3:158519509-158520329 | Island | MFSD1;MFSD1 | Body;Body | 0.8657676 | -4.775907 | 5.837325 | 4.79e-05 | 0.9902250 | -1.7461847 |
| cg03081615 | chr4 | 8007156 |
|
OpenSea | ABLIM2;ABLIM2;ABLIM2;ABLIM2;ABLIM2;ABLIM2 | Body;Body;Body;Body;Body;Body | -0.4584577 | 4.502980 | -5.804816 | 5.06e-05 | 0.9902250 | -1.7589717 | |
| cg27583671 | chr6 | 151670624 |
|
OpenSea | AKAP12;AKAP12 | Body;Body | -0.5434487 | 2.662483 | -5.714514 | 5.91e-05 | 0.9902250 | -1.7950719 | |
| cg17506328 | chr11 | 46868210 |
|
chr11:46867194-46868077 | S_Shore | CKAP5;CKAP5 | TSS1500;TSS1500 | -0.3991032 | -3.443528 | -5.703255 | 6.03e-05 | 0.9902250 | -1.7996333 |
| cg27522733 | chr1 | 193448591 |
|
OpenSea | 0.5184504 | -2.431826 | 5.680984 | 6.26e-05 | 0.9902250 | -1.8086967 | |||
| cg09733297 | chr16 | 54155568 |
|
chr16:54155953-54156180 | N_Shore | -0.5840849 | 3.388928 | -5.658058 | 6.52e-05 | 0.9902250 | -1.8180817 | ||
| cg17427349 | chr15 | 38988004 |
|
OpenSea | C15orf53 | TSS1500 | -0.4379188 | 3.008786 | -5.618339 | 6.98e-05 | 0.9902250 | -1.8344761 | |
| cg19572637 | chr2 | 66660175 |
|
chr2:66660452-66660794 | N_Shore | 0.3265461 | -4.483144 | 5.594407 | 7.28e-05 | 0.9902250 | -1.8444370 | ||
| cg00637015 | chr1 | 31468126 |
|
OpenSea | PUM1;PUM1 | Body;Body | 0.6517028 | 1.224266 | 5.593356 | 7.29e-05 | 0.9902250 | -1.8448759 | |
| cg04263436 | chr6 | 33245467 |
|
chr6:33244677-33245554 | Island | B3GALT4 | 1stExon | 0.4051744 | -5.021011 | 5.568245 | 7.62e-05 | 0.9902250 | -1.8553980 |
saveRDS(design_tw, "gu_design_ilanguage.rds")Converting to beta from M values
bDat = ilogit2(mvals)
bDat_new = getBeta(gmset_flt)
#View(bDat)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.122830 -0.021397 -0.009488 0.044541 0.074036
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.748293 0.026526 28.210 <2e-16 ***
## ADOS -0.000183 0.002141 -0.085 0.933
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05354 on 20 degrees of freedom
## Multiple R-squared: 0.0003651, Adjusted R-squared: -0.04962
## F-statistic: 0.007305 on 1 and 20 DF, p-value: 0.9327summary(lm(gwam ~ Twin_Group + ADOS ))##
## Call:
## lm(formula = gwam ~ Twin_Group + ADOS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.04849 -0.02106 0.00000 0.02106 0.04849
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7854076 0.0321787 24.408 3.04e-10 ***
## Twin_Group2 -0.0187093 0.0391098 -0.478 0.643
## Twin_Group3 -0.0471488 0.0407165 -1.158 0.274
## Twin_Group4 -0.0415169 0.0388848 -1.068 0.311
## Twin_Group5 -0.0058929 0.0388848 -0.152 0.883
## Twin_Group6 -0.0875574 0.0407165 -2.150 0.057 .
## Twin_Group7 -0.1205572 0.0411289 -2.931 0.015 *
## Twin_Group8 0.0289974 0.0466647 0.621 0.548
## Twin_Group9 -0.0140098 0.0382724 -0.366 0.722
## Twin_Group12 0.0321116 0.0415696 0.772 0.458
## Twin_Group13 -0.0490648 0.0382724 -1.282 0.229
## ADOS -0.0008733 0.0023238 -0.376 0.715
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03825 on 10 degrees of freedom
## Multiple R-squared: 0.7449, Adjusted R-squared: 0.4642
## F-statistic: 2.654 on 11 and 10 DF, p-value: 0.06766message("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.122806 -0.023720 -0.006553 0.045959 0.070253
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7567609 0.0303220 24.957 <2e-16 ***
## SRS -0.0001323 0.0003537 -0.374 0.712
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05337 on 20 degrees of freedom
## Multiple R-squared: 0.006948, Adjusted R-squared: -0.0427
## F-statistic: 0.1399 on 1 and 20 DF, p-value: 0.7123plot(gwam~SRS)
abline(mylm)
summary(lm(gwam ~ Twin_Group + SRS ))##
## Call:
## lm(formula = gwam ~ Twin_Group + SRS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.03950 -0.01778 0.00000 0.01778 0.03950
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8859084 0.0645649 13.721 8.21e-08 ***
## Twin_Group2 -0.0662112 0.0418050 -1.584 0.14432
## Twin_Group3 -0.0791515 0.0367628 -2.153 0.05677 .
## Twin_Group4 -0.1000688 0.0479946 -2.085 0.06366 .
## Twin_Group5 -0.0324082 0.0361359 -0.897 0.39088
## Twin_Group6 -0.1052230 0.0342744 -3.070 0.01184 *
## Twin_Group7 -0.1482176 0.0357521 -4.146 0.00199 **
## Twin_Group8 0.0146529 0.0336449 0.436 0.67244
## Twin_Group9 -0.0120570 0.0335851 -0.359 0.72706
## Twin_Group12 -0.0346958 0.0478033 -0.726 0.48460
## Twin_Group13 -0.1375186 0.0601013 -2.288 0.04516 *
## SRS -0.0009558 0.0005361 -1.783 0.10495
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03356 on 10 degrees of freedom
## Multiple R-squared: 0.8037, Adjusted R-squared: 0.5877
## F-statistic: 3.721 on 11 and 10 DF, p-value: 0.02376message("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.119997 -0.021356 -0.009054 0.046267 0.071705
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.748427 0.015795 47.38 <2e-16 ***
## motor -0.002399 0.012018 -0.20 0.844
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0535 on 20 degrees of freedom
## Multiple R-squared: 0.001988, Adjusted R-squared: -0.04791
## F-statistic: 0.03984 on 1 and 20 DF, p-value: 0.8438summary(lm(gwam ~ Twin_Group + motor ))##
## Call:
## lm(formula = gwam ~ Twin_Group + motor)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.05067 -0.01489 0.00000 0.01489 0.05067
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.84340 0.05648 14.934 3.65e-08 ***
## Twin_Group2 -0.05404 0.04375 -1.235 0.24499
## Twin_Group3 -0.08466 0.04375 -1.935 0.08173 .
## Twin_Group4 -0.10344 0.06187 -1.672 0.12549
## Twin_Group5 -0.07305 0.06187 -1.181 0.26501
## Twin_Group6 -0.15734 0.06187 -2.543 0.02921 *
## Twin_Group7 -0.12623 0.03572 -3.534 0.00541 **
## Twin_Group8 0.01895 0.03572 0.531 0.60725
## Twin_Group9 -0.01445 0.03572 -0.404 0.69440
## Twin_Group12 -0.03854 0.06187 -0.623 0.54725
## Twin_Group13 -0.11317 0.06187 -1.829 0.09730 .
## motor -0.03227 0.02526 -1.278 0.23024
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03572 on 10 degrees of freedom
## Multiple R-squared: 0.7776, Adjusted R-squared: 0.5329
## F-statistic: 3.178 on 11 and 10 DF, p-value: 0.03951message("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.112240 -0.023919 -0.007083 0.051046 0.069996
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.76440 0.02075 36.833 <2e-16 ***
## diagnosis -0.01426 0.01377 -1.036 0.312
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05217 on 20 degrees of freedom
## Multiple R-squared: 0.05095, Adjusted R-squared: 0.003499
## F-statistic: 1.074 on 1 and 20 DF, p-value: 0.3125summary(lm(gwam ~ Twin_Group + diagnosis ))##
## Call:
## lm(formula = gwam ~ Twin_Group + diagnosis)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.04520 -0.02109 0.00000 0.02109 0.04520
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.789789 0.029004 27.231 1.03e-10 ***
## Twin_Group2 -0.016300 0.037614 -0.433 0.6740
## Twin_Group3 -0.046923 0.037614 -1.248 0.2406
## Twin_Group4 -0.044363 0.037614 -1.179 0.2655
## Twin_Group5 -0.008513 0.037102 -0.229 0.8231
## Twin_Group6 -0.087332 0.037614 -2.322 0.0426 *
## Twin_Group7 -0.115302 0.039109 -2.948 0.0146 *
## Twin_Group8 0.029886 0.039109 0.764 0.4624
## Twin_Group9 -0.019912 0.037614 -0.529 0.6081
## Twin_Group12 0.031464 0.037614 0.837 0.4224
## Twin_Group13 -0.048628 0.037102 -1.311 0.2193
## diagnosis -0.010932 0.012367 -0.884 0.3975
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0371 on 10 degrees of freedom
## Multiple R-squared: 0.76, Adjusted R-squared: 0.4961
## F-statistic: 2.879 on 11 and 10 DF, p-value: 0.05338message("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.120219 -0.027573 0.001011 0.037086 0.071548
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.74625 0.01111 67.148 <2e-16 ***
## iiq -0.01199 0.01137 -1.054 0.304
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05213 on 20 degrees of freedom
## Multiple R-squared: 0.05263, Adjusted R-squared: 0.005262
## F-statistic: 1.111 on 1 and 20 DF, p-value: 0.3044plot(gwam~iiq)
abline(mylm)
summary(lm(gwam ~ Twin_Group + iiq ))##
## Call:
## lm(formula = gwam ~ Twin_Group + iiq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.03390 -0.01728 0.00000 0.01728 0.03390
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.789099 0.023654 33.361 1.38e-11 ***
## Twin_Group2 -0.048409 0.035303 -1.371 0.20028
## Twin_Group3 -0.081258 0.035746 -2.273 0.04632 *
## Twin_Group4 -0.051541 0.033261 -1.550 0.15228
## Twin_Group5 0.001008 0.032995 0.031 0.97624
## Twin_Group6 -0.093611 0.032646 -2.867 0.01674 *
## Twin_Group7 -0.127724 0.032652 -3.912 0.00291 **
## Twin_Group8 0.022448 0.032691 0.687 0.50789
## Twin_Group9 -0.041579 0.035398 -1.175 0.26736
## Twin_Group12 0.009276 0.033716 0.275 0.78883
## Twin_Group13 -0.059985 0.033142 -1.810 0.10042
## iiq -0.019254 0.009715 -1.982 0.07564 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03264 on 10 degrees of freedom
## Multiple R-squared: 0.8142, Adjusted R-squared: 0.6099
## F-statistic: 3.985 on 11 and 10 DF, p-value: 0.01887message("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.121403 -0.026341 -0.004373 0.038349 0.071034
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.746247 0.011286 66.120 <2e-16 ***
## ilanguage -0.007916 0.011552 -0.685 0.501
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05294 on 20 degrees of freedom
## Multiple R-squared: 0.02294, Adjusted R-squared: -0.02592
## F-statistic: 0.4695 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.04157 -0.01939 0.00000 0.01939 0.04157
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.786520 0.027464 28.639 6.27e-11 ***
## Twin_Group2 -0.042647 0.043442 -0.982 0.3494
## Twin_Group3 -0.071648 0.042542 -1.684 0.1231
## Twin_Group4 -0.050797 0.039208 -1.296 0.2242
## Twin_Group5 -0.009173 0.037006 -0.248 0.8092
## Twin_Group6 -0.096264 0.037192 -2.588 0.0270 *
## Twin_Group7 -0.126365 0.036999 -3.415 0.0066 **
## Twin_Group8 0.024588 0.037505 0.656 0.5269
## Twin_Group9 -0.030189 0.040786 -0.740 0.4762
## Twin_Group12 0.011134 0.040393 0.276 0.7884
## Twin_Group13 -0.051645 0.037145 -1.390 0.1946
## ilanguage -0.011549 0.012592 -0.917 0.3806
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.037 on 10 degrees of freedom
## Multiple R-squared: 0.7614, Adjusted R-squared: 0.4988
## F-statistic: 2.9 on 11 and 10 DF, p-value: 0.05224probes_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)
## $cg12565057
## NULL
##
## $cg00546248
## NULL
##
## $cg15740041
## NULL
##
## $cg05116656
## NULL
##
## $cg11632273
## NULL
##
## $cg01656620
## NULL
##
## $cg10881749
## NULL
##
## $cg13699963
## NULL
##
## $cg22656126
## NULLTOPPROBES <- head(res$Name,4)
par(mfrow=c(2,2))
sapply(TOPPROBES,effect_plot)
## $cg12565057
## NULL
##
## $cg00546248
## NULL
##
## $cg15740041
## NULL
##
## $cg05116656
## NULLsig <- subset(res,P.Value<1e-4)
sig <- sig[order(-abs(sig$logFC)),]
TOPPROBES <- head(sig$Name,4)
sapply(TOPPROBES,effect_plot)
## $cg22656126
## NULL
##
## $cg11632273
## NULL
##
## $cg01656620
## NULL
##
## $cg15740041
## NULLpdf("effect_guthrie.pdf")
par(mfrow=c(2,2))
sapply(TOPPROBES,effect_plot)## $cg22656126
## NULL
##
## $cg11632273
## NULL
##
## $cg01656620
## NULL
##
## $cg15740041
## NULLdev.off()## X11cairo
## 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
## hsa04974 Protein digestion and absorption 92 6 0.04575168
## hsa05222 Small cell lung cancer 88 6 0.05318918
## hsa05217 Basal cell carcinoma 63 5 0.06104022
## hsa05165 Human papillomavirus infection 314 15 0.06331859
## hsa04928 Parathyroid hormone synthesis, secretion and action 105 7 0.06457119
## hsa04350 TGF-beta signaling pathway 106 6 0.07901649
## hsa04930 Type II diabetes mellitus 45 4 0.08183071
## hsa05216 Thyroid cancer 37 3 0.11608492
## hsa03420 Nucleotide excision repair 57 3 0.11993290
## hsa04658 Th1 and Th2 cell differentiation 88 5 0.12312867
## hsa04976 Bile secretion 86 4 0.12597522
## hsa04512 ECM-receptor interaction 87 5 0.14250289
## hsa04927 Cortisol synthesis and secretion 63 4 0.15138068
## hsa01240 Biosynthesis of cofactors 147 5 0.16380357
## hsa04934 Cushing syndrome 154 7 0.22946999
## hsa04514 Cell adhesion molecules 142 6 0.24322446
## hsa04120 Ubiquitin mediated proteolysis 134 5 0.24361552
## hsa05200 Pathways in cancer 507 18 0.24738145
## hsa05146 Amoebiasis 98 4 0.25389948
## hsa04922 Glucagon signaling pathway 99 4 0.25989098
## FDR
## hsa04974 1
## hsa05222 1
## hsa05217 1
## hsa05165 1
## hsa04928 1
## hsa04350 1
## hsa04930 1
## hsa05216 1
## hsa03420 1
## hsa04658 1
## hsa04976 1
## hsa04512 1
## hsa04927 1
## hsa01240 1
## hsa04934 1
## hsa04514 1
## hsa04120 1
## hsa05200 1
## hsa05146 1
## hsa04922 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:0048625 BP myoblast fate commitment 6
## GO:0031526 CC brush border membrane 54
## GO:0033077 BP T cell differentiation in thymus 83
## GO:0033089 BP positive regulation of T cell differentiation in thymus 10
## GO:0033157 BP regulation of intracellular protein transport 214
## GO:0006517 BP protein deglycosylation 25
## GO:0090160 BP Golgi to lysosome transport 12
## GO:0002076 BP osteoblast development 17
## GO:0005903 CC brush border 98
## GO:0045061 BP thymic T cell selection 22
## GO:0045059 BP positive thymic T cell selection 14
## GO:0071218 BP cellular response to misfolded protein 21
## GO:0019827 BP stem cell population maintenance 170
## GO:0045582 BP positive regulation of T cell differentiation 112
## GO:0038065 BP collagen-activated signaling pathway 14
## GO:0098727 BP maintenance of cell number 174
## GO:0061101 BP neuroendocrine cell differentiation 16
## GO:0090316 BP positive regulation of intracellular protein transport 147
## GO:0005902 CC microvillus 92
## GO:0051788 BP response to misfolded protein 23
## DE P.DE FDR
## GO:0048625 3 0.001076981 1
## GO:0031526 6 0.001785434 1
## GO:0033077 8 0.002671078 1
## GO:0033089 3 0.002673410 1
## GO:0033157 13 0.003118817 1
## GO:0006517 4 0.003283324 1
## GO:0090160 3 0.003837176 1
## GO:0002076 4 0.004035289 1
## GO:0005903 8 0.006108216 1
## GO:0045061 4 0.006922502 1
## GO:0045059 3 0.009118887 1
## GO:0071218 3 0.010242147 1
## GO:0019827 11 0.010273633 1
## GO:0045582 8 0.010693065 1
## GO:0038065 3 0.010725563 1
## GO:0098727 11 0.012030153 1
## GO:0061101 3 0.013676681 1
## GO:0090316 9 0.014310983 1
## GO:0005902 7 0.014552391 1
## GO:0051788 3 0.015264279 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 = "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 cachehead(results.ranges, 20)## seqnames start end width strand no.cpgs min_smoothed_fdr Stouffer
## 1 chr1 234792390 234792752 363 * 2 8.755567e-07 1
## 2 chr2 46525649 46525678 30 * 2 9.543493e-04 1
## 3 chr6 91113506 91113687 182 * 2 8.564702e-05 1
## 4 chr9 96077410 96077420 11 * 2 9.668722e-04 1
## 5 chr11 44613980 44614255 276 * 2 4.084688e-05 1
## 6 chr22 20102216 20102580 365 * 3 1.418446e-06 1
## 7 chr12 75784541 75785295 755 * 11 8.243423e-09 1
## 8 chr2 33359198 33359688 491 * 10 9.097454e-04 1
## 9 chr22 31318147 31318444 298 * 7 9.447998e-04 1
## 10 chr7 100481990 100482960 971 * 5 3.223939e-08 1
## 11 chr16 87780848 87781501 654 * 5 5.201435e-07 1
## 12 chr11 65820391 65820516 126 * 4 4.325928e-04 1
## 13 chr12 124788627 124789036 410 * 4 2.317340e-05 1
## HMFDR Fisher maxdiff meandiff overlapping.genes
## 1 0.9999951 1 0.011843398 0.007349988 RP4-781K5.4
## 2 0.9999951 1 -0.003715122 -0.002901168 EPAS1
## 3 0.9999951 1 0.015641676 0.011832067 <NA>
## 4 0.9999951 1 -0.004274340 -0.001878717 WNK2
## 5 0.9999951 1 0.017079330 0.015007099 CD82
## 6 0.9999951 1 -0.012596545 -0.003140646 TRMT2A
## 7 0.9999951 1 0.006045794 0.001967226 GLIPR1L2, CAPS2
## 8 0.9999951 1 0.015386270 0.008549527 LTBP1
## 9 0.9999951 1 0.009559386 0.003625967 MORC2-AS1
## 10 0.9999951 1 0.005702107 0.002285903 SRRT
## 11 0.9999951 1 0.015904041 0.006726073 KLHDC4, RP11-278A23.2
## 12 0.9999951 1 0.006046969 -0.001162702 SF3B2
## 13 0.9999951 1 -0.009740326 -0.002426397 FAM101Awrite.csv(results.ranges, file = "dmrcoutput_guthrie.csv", row.names = TRUE)Session information
sessionInfo()## R version 4.3.2 (2023-10-31)
## 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] IlluminaHumanMethylation450kmanifest_0.4.0
## [3] IlluminaHumanMethylationEPICmanifest_0.3.0
## [4] vioplot_0.4.0
## [5] zoo_1.8-12
## [6] sm_2.2-5.7.1
## [7] kableExtra_1.3.4
## [8] mitch_1.12.0
## [9] FlowSorted.Blood.EPIC_2.4.2
## [10] ExperimentHub_2.8.1
## [11] AnnotationHub_3.8.0
## [12] BiocFileCache_2.8.0
## [13] dbplyr_2.3.3
## [14] DMRcate_2.14.1
## [15] ggplot2_3.4.3
## [16] reshape2_1.4.4
## [17] FlowSorted.Blood.450k_1.38.0
## [18] WGCNA_1.72-1
## [19] fastcluster_1.2.3
## [20] dynamicTreeCut_1.63-1
## [21] gplots_3.1.3
## [22] RColorBrewer_1.1-3
## [23] ruv_0.9.7.1
## [24] IlluminaHumanMethylationEPICanno.ilm10b2.hg19_0.6.0
## [25] limma_3.56.2
## [26] missMethyl_1.34.0
## [27] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [28] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1
## [29] minfi_1.46.0
## [30] bumphunter_1.42.0
## [31] locfit_1.5-9.8
## [32] iterators_1.0.14
## [33] foreach_1.5.2
## [34] Biostrings_2.68.1
## [35] XVector_0.40.0
## [36] SummarizedExperiment_1.30.2
## [37] Biobase_2.60.0
## [38] MatrixGenerics_1.12.3
## [39] matrixStats_1.0.0
## [40] GenomicRanges_1.52.0
## [41] GenomeInfoDb_1.36.3
## [42] IRanges_2.34.1
## [43] S4Vectors_0.38.1
## [44] 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.2 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 labeling_0.4.3
## [25] sass_0.4.7 readr_2.1.4
## [27] genefilter_1.82.1 askpass_1.2.0
## [29] systemfonts_1.0.4 Rsamtools_2.16.0
## [31] foreign_0.8-85 svglite_2.1.1
## [33] siggenes_1.74.0 illuminaio_0.42.0
## [35] R.utils_2.12.2 dichromat_2.0-0.1
## [37] scrime_1.3.5 BSgenome_1.68.0
## [39] readxl_1.4.3 rstudioapi_0.15.0
## [41] impute_1.74.1 RSQLite_2.3.1
## [43] generics_0.1.3 BiocIO_1.10.0
## [45] gtools_3.9.4 dplyr_1.1.3
## [47] GO.db_3.17.0 Matrix_1.6-1.1
## [49] interp_1.1-4 fansi_1.0.4
## [51] abind_1.4-5 R.methodsS3_1.8.2
## [53] lifecycle_1.0.3 edgeR_3.42.4
## [55] yaml_2.3.7 rhdf5_2.44.0
## [57] grid_4.3.2 blob_1.2.4
## [59] promises_1.2.1 crayon_1.5.2
## [61] lattice_0.22-5 echarts4r_0.4.5
## [63] GenomicFeatures_1.52.2 annotate_1.78.0
## [65] KEGGREST_1.40.0 pillar_1.9.0
## [67] knitr_1.44 beanplot_1.3.1
## [69] rjson_0.2.21 codetools_0.2-19
## [71] glue_1.6.2 data.table_1.14.8
## [73] vctrs_0.6.3 png_0.1-8
## [75] cellranger_1.1.0 gtable_0.3.4
## [77] cachem_1.0.8 xfun_0.40
## [79] mime_0.12 S4Arrays_1.0.6
## [81] survival_3.5-7 statmod_1.5.0
## [83] ellipsis_0.3.2 interactiveDisplayBase_1.38.0
## [85] nlme_3.1-163 bit64_4.0.5
## [87] bsseq_1.36.0 progress_1.2.2
## [89] filelock_1.0.2 bslib_0.5.1
## [91] nor1mix_1.3-0 KernSmooth_2.23-22
## [93] rpart_4.1.21 colorspace_2.1-0
## [95] DBI_1.1.3 Hmisc_5.1-1
## [97] nnet_7.3-19 tidyselect_1.2.0
## [99] bit_4.0.5 compiler_4.3.2
## [101] curl_5.0.2 rvest_1.0.3
## [103] htmlTable_2.4.1 BiasedUrn_2.0.11
## [105] xml2_1.3.5 DelayedArray_0.26.7
## [107] rtracklayer_1.60.1 checkmate_2.2.0
## [109] scales_1.2.1 caTools_1.18.2
## [111] quadprog_1.5-8 rappdirs_0.3.3
## [113] stringr_1.5.0 digest_0.6.33
## [115] rmarkdown_2.24 GEOquery_2.68.0
## [117] htmltools_0.5.6 pkgconfig_2.0.3
## [119] jpeg_0.1-10 base64enc_0.1-3
## [121] sparseMatrixStats_1.12.2 highr_0.10
## [123] fastmap_1.1.1 ensembldb_2.24.0
## [125] rlang_1.1.1 htmlwidgets_1.6.2
## [127] shiny_1.7.5 DelayedMatrixStats_1.22.6
## [129] farver_2.1.1 jquerylib_0.1.4
## [131] jsonlite_1.8.7 BiocParallel_1.34.2
## [133] mclust_6.0.0 R.oo_1.25.0
## [135] VariantAnnotation_1.46.0 RCurl_1.98-1.12
## [137] magrittr_2.0.3 Formula_1.2-5
## [139] GenomeInfoDbData_1.2.10 Rhdf5lib_1.22.1
## [141] munsell_0.5.0 Rcpp_1.0.11
## [143] stringi_1.7.12 zlibbioc_1.46.0
## [145] MASS_7.3-60 plyr_1.8.8
## [147] org.Hs.eg.db_3.17.0 deldir_1.0-9
## [149] splines_4.3.2 multtest_2.56.0
## [151] hms_1.1.3 rngtools_1.5.2
## [153] biomaRt_2.56.1 BiocVersion_3.17.1
## [155] XML_3.99-0.14 evaluate_0.21
## [157] latticeExtra_0.6-30 biovizBase_1.48.0
## [159] BiocManager_1.30.22 httpuv_1.6.11
## [161] tzdb_0.4.0 tidyr_1.3.0
## [163] openssl_2.1.0 purrr_1.0.2
## [165] reshape_0.8.9 xtable_1.8-4
## [167] restfulr_0.0.15 AnnotationFilter_1.24.0
## [169] later_1.3.1 viridisLite_0.4.2
## [171] tibble_3.2.1 beeswarm_0.4.0
## [173] memoise_2.0.1 AnnotationDbi_1.62.2
## [175] GenomicAlignments_1.36.0 cluster_2.1.4