Where’d my DMR’s go?

So Steven noticed an interesting thing, DMRs that are significant in the 10x cutoff aren’t present in the 5x coverage. I’m thinking this is a modeling effect, but I wanted to check to make sure my data subsetting was actually correct first.

I read in the methylation count data I use in Macau for each of the cutoffs.

library(readr)
macau_meth10x <- read_delim("~/Documents/Genewiz_hdd/Day10/macau/macau.meth10x.txt", " ", escape_double = FALSE, trim_ws = TRUE)

macau_meth7x <- read_csv("~/Documents/Genewiz_hdd/Day10/macauv2/macau.meth7x.txt")

macau_meth5x <- read_csv("~/Documents/Genewiz_hdd/Day10/macauv2/macau.meth5x.txt")

macau_meth3x <- read_csv("~/Documents/Genewiz_hdd/Day10/macauv2/macau.meth3x.txt")

Then I check to see how many rows I’m dealing with, this is equivilant to the number of potential loci that meet our data cutoffs.

nrow(macau_meth10x)

[1] 1977

nrow(macau_meth7x)

[1] 4250

nrow(macau_meth5x)

[1] 8972

nrow(macau_meth3x)

[1] 23284

Then I check to see that all of my 10x coverage sites are in the 7x, 5x, and 3x sites.

The code asks what is the length of the vector returned by looking at which sites co-occur in both the 7x and 10x data sets.

This should equal the length of the smaller dataset, if they are all present in the larger data set, and it does.

nrow(macau_meth10x)

[1] 1977

length(which(macau_meth7x$Site %in% macau_meth10x$Site == TRUE))

[1] 1977

length(which(macau_meth5x$Site %in% macau_meth10x$Site == TRUE))

[1] 1977

length(which(macau_meth3x$Site %in% macau_meth10x$Site == TRUE))

[1] 1977

Now I make sure that all of the 7x are in the 5x and 3x sets, and they are.

nrow(macau_meth7x)

[1] 4250

length(which(macau_meth5x$Site %in% macau_meth7x$Site == TRUE))

[1] 4250

length(which(macau_meth3x$Site %in% macau_meth7x$Site == TRUE))

[1] 4250

And that the 5x is in the 3x set. And it is.

nrow(macau_meth5x)

[1] 8972

length(which(macau_meth3x$Site %in% macau_meth5x$Site == TRUE))

[1] 8972

Next. Lets make sure that something didn’t get wonky with the count values themselves.

test <- as.data.frame(macau_meth10x$Site)
test$EPI103_10x <- macau_meth10x$`EPI-103`
test$EPI103_5x <- macau_meth5x$EPI103Meth[which(macau_meth5x$Site %in% macau_meth10x$Site == TRUE)]
test$EPI104_10x <- macau_meth10x$`EPI-104`
test$EPI104_5x <- macau_meth5x$`EPI-104Meth`[which(macau_meth5x$Site %in% macau_meth10x$Site == TRUE)]
test$EPI111_10x <- macau_meth10x$`EPI-111`
test$EPI111_5x <- macau_meth5x$`EPI-111Meth`[which(macau_meth5x$Site %in% macau_meth10x$Site == TRUE)]
head(test)

It looks like count values also match between a subset of samples at the 10x and 5x range

What does this tell us? That it’s likely a modeling effect. Next I’ll look at the outputs from Macau for each coverage level, and see what the P-values for individual sites look like, maybe it’s just an effect of borderline close p-values getting cut off in different runs.

TenX_results <- read_delim("~/Documents/Genewiz_hdd/Day10/macau/output/result.assoc.txt", 
    "\t", escape_double = FALSE, trim_ws = TRUE)

SevenX_assoc <- read_delim("~/Documents/Genewiz_hdd/Day10/macauv2/output/7x.assoc.txt", 
    "\t", escape_double = FALSE, trim_ws = TRUE)

FiveX_assoc <- read_delim("~/Documents/Genewiz_hdd/Day10/macauv2/output/5x.assoc.txt", 
    "\t", escape_double = FALSE, trim_ws = TRUE)

ThreeX_assoc <- read_delim("~/Documents/Genewiz_hdd/Day10/macauv2/output/3x.assoc.txt", 
    "\t", escape_double = FALSE, trim_ws = TRUE)

Lets check to make sure that all of the loci deemed significant by the 10X test are still located in the 7x, 5x, and 3x files.

TenX.sig.results <- TenX_results[TenX_results$pvalue <= 0.05,]
length(which(TenX.sig.results$id %in% SevenX_assoc$id == TRUE))

[1] 41

length(which(TenX.sig.results$id %in% FiveX_assoc$id == TRUE))

[1] 41

length(which(TenX.sig.results$id %in% ThreeX_assoc$id == TRUE))

[1] 41

They are. Now lets look at p-values.

TenvFive.p <- as.data.frame(TenX.sig.results$id)
TenvFive.p$TenXp <- TenX.sig.results$pvalue
TenvFive.p$FiveXp <- FiveX_assoc$pvalue[which(TenX.sig.results$id %in% FiveX_assoc$id == TRUE)]
print(TenvFive.p)

TenvThree.p <- as.data.frame(TenX.sig.results$id)
TenvThree.p$TenXp <- TenX.sig.results$pvalue
TenvThree.p$ThreeXp <- ThreeX_assoc$pvalue[which(TenX.sig.results$id %in% ThreeX_assoc$id == TRUE)]
print(TenvThree.p)

So it looks like all of the DMR sites identified in the 10x coverage run still exist in the other runs, it’s just an effect of the model being run that removes them from the DMR list.