Testing base methylKit code

Alright, it’s time to run methylKit on mox! Steven and Sam suggested I run my methylKit code in this discussion, so I’ll finish testing my code, make some modifications, and create a script to run on mox.

Percent difference for DML

The first thing I needed to make a decision about was the percent difference I wanted to use as a threshold for DML identification. In this discussion, I pointed out that some papers use a 25% difference, while others (including my C. virginica paper) use a 50% difference. Initially, I set up my code to include 25%, 50%, 75%, and 99% differences in methylation between treatments. After loading in my saved .RData, I ran getMethylDiff using objects generated from calculateDiffMeth with no additional modifications and just the overdispersion correction.

Using a 25% or 50% methylation difference, I got > 1000 DML when not applying any covariate matrix or overdispersion correction! The number of DML identified were in the 100s when using a 75% or 99% difference. When I used the output from calculateDiffMeth with the overdispersion correcrtion, I identified 2119 DML that were 25% different, 564 that were 50% different, and 80 that were 75% different. Clearly making the model more complicated reduces the amount of DML identified. I think if I used output with a covariate matrix and overdispersion correction, a 25% or 50% methylation difference wouldn’t get me more than ~500 DML. I set up the DML identification code to use both a 25% and 50% methylation difference so I can compare the output from mox.

Testing randomization code

The main parts of my R Markdown script I needed to still test were the randomization trial and subsequent statistical testing. I started by running a modified version of the randomization trial that did not include any overdispersion or covariate information for calculateDiffMeth:

for (i in 1:1000) { #For each iteration
  pHRandTreat <- sample(sampleMetadata$pHTreatment,
                        size = length(sampleMetadata$pHTreatment),
                        replace = FALSE) #Randomize treatment information
  pHRandDML <- methylKit::reorganize(processedFiles,
                                     sample.ids = sampleMetadata$sampleID,
                                     treatment = pHRandTreat) %>%
    methylKit::filterByCoverage(., lo.count = 5, lo.perc = NULL,
                                hi.count = NULL, hi.perc = 99.9) %>%
    methylKit::normalizeCoverage(.) %>%
    methylKit::unite(., destrand = FALSE) %>%
    methylKit::calculateDiffMeth(.) %>%
    methylKit::getMethylDiff(., difference = 50, qvalue = 0.01) %>%
    nrow() -> pHRandTest50[i] ##Reorganize, filter, normalize, and unite samples. Calculate diffMeth and obtain DML that are 50% and have a q-value of 0.01. Save the number of DML identified
}

I thought the code would finish running in a few hours, but since it was taking a while and I didn’t get any error messages, I interrupted the code chunk and commented it out. I then set up two randomization tests: one for a 25% methylation difference, and another for a 50% methylation difference:

pHRandTest25 <- NULL #Create an empty matrix to store randomization trial results for 25% difference

for (i in 1:1000) { #For each iteration
  pHRandTreat <- sample(sampleMetadata$pHTreatment,
                        size = length(sampleMetadata$pHTreatment),
                        replace = FALSE) #Randomize treatment information
  pHRandDML <- methylKit::reorganize(processedFiles,
                                     sample.ids = sampleMetadata$sampleID,
                                     treatment = pHRandTreat) %>%
    methylKit::filterByCoverage(., lo.count = 5, lo.perc = NULL,
                                hi.count = NULL, hi.perc = 99.9) %>%
    methylKit::normalizeCoverage(.) %>%
    methylKit::unite(., destrand = FALSE) %>%
    methylKit::calculateDiffMeth(., covariates = covariateMetadata,
                                 overdispersion = "MN",
                                 test = "Chisq") %>%
    methylKit::getMethylDiff(., difference = 25, qvalue = 0.01) %>%
    nrow() -> pHRandTest25[i] ##Reorganize, filter, normalize, and unite samples. Calculate diffMeth and obtain DML that are 50% and have a q-value of 0.01. Save the number of DML identified
}

pHRandTest50 <- NULL #Create an empty matrix to store randomization trial results for 50% difference

for (i in 1:1000) { #For each iteration
  pHRandTreat <- sample(sampleMetadata$pHTreatment,
                        size = length(sampleMetadata$pHTreatment),
                        replace = FALSE) #Randomize treatment information
  pHRandDML <- methylKit::reorganize(processedFiles,
                                     sample.ids = sampleMetadata$sampleID,
                                     treatment = pHRandTreat) %>%
    methylKit::filterByCoverage(., lo.count = 5, lo.perc = NULL,
                                hi.count = NULL, hi.perc = 99.9) %>%
    methylKit::normalizeCoverage(.) %>%
    methylKit::unite(., destrand = FALSE) %>%
    methylKit::calculateDiffMeth(., covariates = covariateMetadata,
                                 overdispersion = "MN",
                                 test = "Chisq") %>%
    methylKit::getMethylDiff(., difference = 50, qvalue = 0.01) %>%
    nrow() -> pHRandTest50[i] ##Reorganize, filter, normalize, and unite samples. Calculate diffMeth and obtain DML that are 50% and have a q-value of 0.01. Save the number of DML identified
}

The output of the randomization test is a vector with the number of DML identified in each trial. I needed a way to determine if the number of DML I identified with getMethylDiff was significantly greater than what may be identified by chance. I dug deep into my stats brain (you know, that part of my brain that barely turns on) and decided to run a one-tailed t-test:

pHRandTestResults50 <- t.test(pHRandTest50, alternative = "less", mu = nrow(diffMethStatsTreatment50)) #Conduct a one-sample t-test to determine the probability of identifying more DML by chance. mu is the number of DML identified with a 50% methylation difference
pHRandTestResults50$statistic #t
pHRandTestResults50$parameter #df
pHRandTestResults50$p.value #P-value

I also wrote code to create a histogram with the random distribution of DML identified, and include a vertical line for the number of DML obtained with getMethylDiff:

jpeg(filename = "../output/06-methylKit/rand-test/Random-Distribution-diff50.jpeg", height = 1000, width = 1000) #Save file with designated name
hist(pHRandTest50, plot = TRUE, main = "", xlab = "Number of DML (50% difference)") #Plot a histogram with randomization trial results
dev.off()

Now that I’ve tested my base code, I will work on separating female and indeterminate samples and creating the actual mox script!

Going forward

1. Separate female and indeterminate samples
2. Write a mox script
3. Run R script on mox
4. Update methods
5. Obtain relatedness matrix and SNPs with EpiDiverse/snp
6. Revise methylKit study design: separate tests for indeterminate and female oysters
7. Write results
8. Identify genomic location of DML
9. Determine if RNA should be extracted
10. Determine if larval DNA/RNA should be extracted