Ok, time to re-analyze the parental data using DSS. I grabbed the bismark .cov.gz files for the parents of our offspring samples, then replicated the 10-diff-methyl-DSS pipeline to (a) identify treatment-associated DMLs and DMRs while controlling for sex, (b) identify DMLs with sex-dependent signal (effect differs between M and F), (c) identify within-sex DMLs and DMRs (using smoothed test), and (d) compare these updated parental DMLs to the offspring DMLs identified in 10-diff-methyl-DSS to look for overlap (evidence of inheritance)
Results
n=26 parent samples:
| treatment | sex | n |
|---|---|---|
| Control | Female | 8 |
| Control | Male | 4 |
| Exposed | Female | 8 |
| Exposed | Male | 6 |
In intial setup, 34,343,214 methylation loci with 19,684,036 loci after filtering (this is compared to 31,450,950 loci filtered to 7,594,836 in the offspring data - difference is probably related to low-coverage issues in the offspring data).
MULTI-FACTOR:
Parental treatment effect (sex-adjusted): 2,237 DMLs and 1,871 DMRs
Sex-dependence of treatment effect: 7,217 loci show a significant sex-dependent treatment effect
SEX-SPECIFIC
Female: DMLs = 10,584 DMRs = 8,489
Male: DMLs = 67,461 DMRs = 22,773
PARENT-OFFSPRING OVERLAP
- Egg/Zygote overlap: 100 shared DMLs, 94% directionally concordant
- Egg/Larvae overlap: 262 shared DMLs, 96.9% directionally concordant
- Sperm/Zygote overlap: 1,120 shared DMLs, 99.4% directionally concordant
- Sperm/Larvae overlap: 1,022 shared DMLs, 99.4 directionally concordant
Super cool! Re-calling the parental data with DSS (and using the lowered methylation difference of 25%, as with the offspring DMLs) led to more than 10X the number of overlapping DMLs! Again, the overlap is largely with paternal contribution (sperm), though there are also more overlapping DMLs with the egg data, and overlap across all 4 comparisons is still almost exclusively directionally-concordant – great support for inheritance!
Also of interest:
Of the 7,287 and 5,330 DMLs identified in zygotes and larvae, respectively, only 335 are present in both lifestages. However, of those stage-shared DMLS, ~97% are dirtectionally concordant! This is really interesting potential evidence of partial epigenetic reprogramming!
Additionally, of the 10,584 and 67,461 DMLs identified in parent eggs and sperm, respectively, only 1,319 are present in both sexes. Furthermore, within these 1,319 sex-shared DMLs, only 36% are directionally concordant. In other words, The majority of parental DMLs are sex-specific, and of those present in both sexes, the majority of treatment responses point in opposite directions! While both male and female oysters underwent methylomic shifts, the location and direction of those shifts differed dramatically between the two sexes.
Next steps will be to perform an enrichment analysis, to check whether the observed overlap in parent-offspring differential methylation exceeds overlap that we would expect by chance.
Code
See full code at ceasmallr/code/10.1-diff-methyl-DSS-parents
1 Setup
As in 10-diff-methyl-DSS.Rmd, will need to run this .Rmd file from the identically-named .sh script, due to permissions and memory issues.
#!/bin/bash
#SBATCH --account=srlab
#SBATCH --partition=cpu-g2-mem2x
#SBATCH --cpus-per-task=8
#SBATCH --mem=350G
#SBATCH --time=1-00:00:00
#SBATCH --output=dss_parent_%j.out
#SBATCH --error=dss_parent_%j.err
#SBATCH --chdir=/gscratch/srlab/kdurkin1/ceasmallr/code
apptainer exec --bind /gscratch:/gscratch \
/gscratch/srlab/kdurkin1/srlab-R4.4-bioinformatics-container-703094b.sif \
bash -c 'source /srlab/programs/miniforge3-24.7.1-0/etc/profile.d/conda.sh && conda activate /gscratch/srlab/kdurkin1/.conda/envs/dss && Rscript -e "rmarkdown::render(\"10.1-diff-methyl-DSS-parents.Rmd\")"'# Install if needed:
# if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
# BiocManager::install(c("DSS", "bsseq"))
library(DSS)
library(bsseq)
library(readr)
library(dplyr)
library(stringr)
library(tibble)
library(ggplot2)2 Download methylation calls
I believe I’ve found the bismark .cov.gz coverage files for the ceabigr parent gametes stored at: https://gannet.fish.washington.edu/seashell/bu-mox/scrubbed/120321-cvBS/
Importantly, not all of these individuals were included in the crosses that generated our study’s offspring individuals. I need to ID the parents of all of our individuals, get their original sample names, and then download their respective files.
# For all offspring, grab their two parent IDs from filenames
parents=$(ls ../data/bismark-methyl-extraction/*.cov.gz \
| xargs -n1 basename \
| sed -E 's/^([A-Za-z0-9]+)-([A-Za-z0-9]+)-.*/\1\n\2/' \
| sort -u)
# from parent IDs, identify sex (M/F) and treatment (C/E)
# also, from metadata file, use Sample ID to get original Parent ID (used in parent data file names)
: > "../data/bismark-methyl-extraction-parents/parent_lookup.tsv"
printf 'Parent.ID\ttreatment\tsex\tSample.ID\n' >> "../data/bismark-methyl-extraction-parents/parent_lookup.tsv"
sample_ids=()
while read -r p; do
case ${p:0:1} in C) trt=Control;; E) trt=Exposed;; *) trt=NA;; esac
case ${p:1:1} in F) sex=Female;; M) sex=Male;; *) sex=NA;; esac
sid=$(awk -F, -v pid="$p" 'NR>1 && $6==pid {print $1; exit}' "../data/adult-meta.csv")
if [ -z "$sid" ]; then
printf '%s\t%s\t%s\t<NO MATCH>\n' "$p" "$trt" "$sex" >> "../data/bismark-methyl-extraction-parents/parent_lookup.tsv"
echo "WARNING: no metadata match for parent $p , skipping" >&2
else
printf '%s\t%s\t%s\t%s\n' "$p" "$trt" "$sex" "$sid" >> "../data/bismark-methyl-extraction-parents/parent_lookup.tsv"
sample_ids+=("$sid")
fi
done <<< "$parents"
# now that we have list of Parent IDs, download corresponding .cov.gz files from storage
alt=$(IFS='|'; echo "${sample_ids[*]}")
wget -r -np -nd -nc \
--accept-regex "/(${alt})[^/]*\.cov\.gz$" \
-P "../data/bismark-methyl-extraction-parents" "https://gannet.fish.washington.edu/seashell/bu-mox/scrubbed/120321-cvBS/"IMPORTANT: There doesn’t seem to be a record for the parent CM03? While I have two offspring labelled with that parent, CM03 isn’t included in the adult-meta.csv file…
Will need to go back to original repo to try to figure out what’s going on.
cd ../data/bismark-methyl-extraction-parents/
echo "How many methylation calling files were downloaded?"
ls *.deduplicated.bismark.cov.gz | wc -l
echo ""
echo "Check file checksums:"
grep '.deduplicated.bismark.cov.gz' checksums.md5 | md5sum -c -3 Read Bismark coverage into a BSseq object
NOTE: When working with the offspring data, I ended up dropping several low-yeild libraries. While I should still actually check the multiqc report for parental data, all of the bismark .cov.gz files are large enough (>100MB) that I’m not overly worried about low-yield issues. For now I’ll keep all of them in.
file.list <- list.files("../data/bismark-methyl-extraction-parents", pattern = "\\.bismark\\.cov\\.gz$", full.names = TRUE)
# As in 06.2, drop the low-yield libraries (<100M Cs)
# BUT, could keep all samples and let DSS handle the coverage differences -- come back to later.
# drop <- c("CF01-CM01-Zygote", "CF08-CM04-Larvae", "CF08-CM05-Larvae", "EF04-EM04-Zygote", "EF05-EM05-Zygote")
# file.list <- file.list[!str_detect(basename(file.list), str_c(drop, collapse = "|"))]
# format sample IDs
sample.ids <- basename(file.list) %>%
str_replace("_R1_val_1_bismark_bt2_pe\\.deduplicated\\.bismark\\.cov\\.gz$", "")
# Grab treatment and sex info
# Sex info is contained within Sample IDs (M/F)
# Treatment info will require cross-reference with metadata file
adult_meta <- read_csv("../data/adult-meta.csv", show_col_types = FALSE)
meta <- tibble(sample = sample.ids, file = file.list) %>%
# sex straight from the filename ID (trailing F / M)
mutate(sex = case_when(str_ends(sample, "F") ~ "Female",
str_ends(sample, "M") ~ "Male",
TRUE ~ NA_character_)) %>%
# treatment (and metadata Sex, for the cross-check) by Sample.ID
left_join(adult_meta %>% select(Sample.ID, Treatment, Sex),
by = c("sample" = "Sample.ID")) %>%
mutate(
# reference levels: Control = effect of exposure; Female = reference sex
treatment = factor(Treatment, levels = c("Control", "Exposed")),
sex = factor(sex, levels = c("Female", "Male"))
)
knitr::kable(count(meta, treatment, sex))| treatment | sex | n |
|---|---|---|
| Control | Female | 8 |
| Control | Male | 4 |
| Exposed | Female | 8 |
| Exposed | Male | 6 |
colData <- DataFrame(treatment = meta$treatment,
sex = meta$sex,
row.names = meta$sample)
BSobj <- read.bismark(files = meta$file,
colData = colData,
rmZeroCov = TRUE, # drops positions covered in 0 samples
strandCollapse = FALSE, # .cov isn't a cytosine report; leave off
BPPARAM = MulticoreParam(workers = 12, progressbar = TRUE),
verbose = TRUE)## [read.bismark] Parsing files and constructing valid loci ...
## | | | 0% | |=== | 4% | |====== | 8% | |======== | 12% | |=========== | 16% | |============== | 20% | |================= | 24% | |==================== | 28% | |====================== | 32% | |========================= | 36% | |============================ | 40% | |=============================== | 44% | |================================== | 48% | |==================================== | 52% | |======================================= | 56% | |========================================== | 60% | |============================================= | 64% | |================================================ | 68% | |================================================== | 72% | |===================================================== | 76% | |======================================================== | 80% | |=========================================================== | 84% | |============================================================== | 88% | |================================================================ | 92% | |=================================================================== | 96% | |======================================================================| 100%
## Done in 65.6 secs
## [read.bismark] Parsing files and constructing 'M' and 'Cov' matrices ...
## | | | 0% | |=== | 4% | |===== | 8% | |======== | 12% | |=========== | 15% | |============= | 19% | |================ | 23% | |=================== | 27% | |====================== | 31% | |======================== | 35% | |=========================== | 38% | |============================== | 42% | |================================ | 46% | |=================================== | 50% | |====================================== | 54% | |======================================== | 58% | |=========================================== | 62% | |============================================== | 65% | |================================================ | 69% | |=================================================== | 73% | |====================================================== | 77% | |========================================================= | 81% | |=========================================================== | 85% | |============================================================== | 88% | |================================================================= | 92% | |=================================================================== | 96% | |======================================================================| 100%
## Done in 38.8 secs
## [read.bismark] Constructing BSseq object ...BSobj # positions are the UNION across samples; missing positions get N = 0## An object of type 'BSseq' with
## 34343214 methylation loci
## 26 samples
## has not been smoothed
## All assays are in-memory# DSS tolerates missingness
# For now, require non-zero coverage in at least half the samples of each treatment x sex cell to reduce noise
# (MAY WANT TO ADJUST LATER)
cov_mat <- getCoverage(BSobj, type = "Cov") # loci x samples
cells <- interaction(meta$treatment, meta$sex, drop = TRUE)
keep <- Reduce(`&`, lapply(levels(cells), function(cl) {
idx <- which(cells == cl)
rowSums(cov_mat[, idx, drop = FALSE] > 0) >= ceiling(length(idx) / 2)
}))
BSobj <- BSobj[keep, ]
cat("Loci retained after coverage filter:", nrow(BSobj), "\n")## Loci retained after coverage filter: 19684036EDIT: All of the downstream model fitting calls are super intensive and take forever to run. Since I have to run this from Rscript, and thus cannot easily retain R objects between runs, this is making script troubleshooting super time-intensive and annoying. To facilitate troubleshooting, I’m going to add a test run option that will randomly subsample a tiny fraction of the data to proceed with:
# Keeps all samples but randomly thins loci genome-wide
# set test_run to FALSE for full-data runs
test_run <- FALSE
test_frac <- 0.02 # keep 2% of loci
if (test_run) {
set.seed(1)
idx <- sort(sample(nrow(BSobj), size = floor(nrow(BSobj) * test_frac)))
BSobj <- BSobj[idx, ]
cat("TEST RUN — loci subsampled to:", nrow(BSobj), "\n")
}4 Multi-factor model (design-level tests)
Fit an interaction model for tests of sex-dependence (does effect of parental exposure differ between the sexes)
design <- data.frame(treatment = meta$treatment, sex = meta$sex)
# Interaction model: does the parental-treatment effect differ by sex?
fit_int <- DMLfit.multiFactor(BSobj, design = design,
formula = ~ treatment + sex + treatment:sex)## Fitting DML model for CpG site: 100000 , 200000 , 300000 , 400000 , 500000 , 600000 , 700000 , 800000 , 900000 , 1000000 , 1100000 , 1200000 , 1300000 , 1400000 , 1500000 , 1600000 , 1700000 , 1800000 , 1900000 , 2000000 , 2100000 , 2200000 , 2300000 , 2400000 , 2500000 , 2600000 , 2700000 , 2800000 , 2900000 , 3000000 , 3100000 , 3200000 , 3300000 , 3400000 , 3500000 , 3600000 , 3700000 , 3800000 , 3900000 , 4000000 , 4100000 , 4200000 , 4300000 , 4400000 , 4500000 , 4600000 , 4700000 , 4800000 , 4900000 , 5000000 , 5100000 , 5200000 , 5300000 , 5400000 , 5500000 , 5600000 , 5700000 , 5800000 , 5900000 , 6000000 , 6100000 , 6200000 , 6300000 , 6400000 , 6500000 , 6600000 , 6700000 , 6800000 , 6900000 , 7000000 , 7100000 , 7200000 , 7300000 , 7400000 , 7500000 , 7600000 , 7700000 , 7800000 , 7900000 , 8000000 , 8100000 , 8200000 , 8300000 , 8400000 , 8500000 , 8600000 , 8700000 , 8800000 , 8900000 , 9000000 , 9100000 , 9200000 , 9300000 , 9400000 , 9500000 , 9600000 , 9700000 , 9800000 , 9900000 , 10000000 , 10100000 , 10200000 , 10300000 , 10400000 , 10500000 , 10600000 , 10700000 , 10800000 , 10900000 , 11000000 , 11100000 , 11200000 , 11300000 , 11400000 , 11500000 , 11600000 , 11700000 , 11800000 , 11900000 , 12000000 , 12100000 , 12200000 , 12300000 , 12400000 , 12500000 , 12600000 , 12700000 , 12800000 , 12900000 , 13000000 , 13100000 , 13200000 , 13300000 , 13400000 , 13500000 , 13600000 , 13700000 , 13800000 , 13900000 , 14000000 , 14100000 , 14200000 , 14300000 , 14400000 , 14500000 , 14600000 , 14700000 , 14800000 , 14900000 , 15000000 , 15100000 , 15200000 , 15300000 , 15400000 , 15500000 , 15600000 , 15700000 , 15800000 , 15900000 , 16000000 , 16100000 , 16200000 , 16300000 , 16400000 , 16500000 , 16600000 , 16700000 , 16800000 , 16900000 , 17000000 , 17100000 , 17200000 , 17300000 , 17400000 , 17500000 , 17600000 , 17700000 , 17800000 , 17900000 , 18000000 , 18100000 , 18200000 , 18300000 , 18400000 , 18500000 , 18600000 , 18700000 , 18800000 , 18900000 , 19000000 , 19100000 , 19200000 , 19300000 , 19400000 , 19500000 , 19600000 ,# Inspect coefficient names so we test the right columns
colnames(fit_int$X)## [1] "(Intercept)" "treatmentExposed"
## [3] "sexMale" "treatmentExposed:sexMale"4.1 Parental treatment effect
Fit an additive model (no interaction) to test for signal of parental treatment effect while controlling for sex.
fit_add <- DMLfit.multiFactor(BSobj, design = design, formula = ~ treatment + sex)## Fitting DML model for CpG site: 100000 , 200000 , 300000 , 400000 , 500000 , 600000 , 700000 , 800000 , 900000 , 1000000 , 1100000 , 1200000 , 1300000 , 1400000 , 1500000 , 1600000 , 1700000 , 1800000 , 1900000 , 2000000 , 2100000 , 2200000 , 2300000 , 2400000 , 2500000 , 2600000 , 2700000 , 2800000 , 2900000 , 3000000 , 3100000 , 3200000 , 3300000 , 3400000 , 3500000 , 3600000 , 3700000 , 3800000 , 3900000 , 4000000 , 4100000 , 4200000 , 4300000 , 4400000 , 4500000 , 4600000 , 4700000 , 4800000 , 4900000 , 5000000 , 5100000 , 5200000 , 5300000 , 5400000 , 5500000 , 5600000 , 5700000 , 5800000 , 5900000 , 6000000 , 6100000 , 6200000 , 6300000 , 6400000 , 6500000 , 6600000 , 6700000 , 6800000 , 6900000 , 7000000 , 7100000 , 7200000 , 7300000 , 7400000 , 7500000 , 7600000 , 7700000 , 7800000 , 7900000 , 8000000 , 8100000 , 8200000 , 8300000 , 8400000 , 8500000 , 8600000 , 8700000 , 8800000 , 8900000 , 9000000 , 9100000 , 9200000 , 9300000 , 9400000 , 9500000 , 9600000 , 9700000 , 9800000 , 9900000 , 10000000 , 10100000 , 10200000 , 10300000 , 10400000 , 10500000 , 10600000 , 10700000 , 10800000 , 10900000 , 11000000 , 11100000 , 11200000 , 11300000 , 11400000 , 11500000 , 11600000 , 11700000 , 11800000 , 11900000 , 12000000 , 12100000 , 12200000 , 12300000 , 12400000 , 12500000 , 12600000 , 12700000 , 12800000 , 12900000 , 13000000 , 13100000 , 13200000 , 13300000 , 13400000 , 13500000 , 13600000 , 13700000 , 13800000 , 13900000 , 14000000 , 14100000 , 14200000 , 14300000 , 14400000 , 14500000 , 14600000 , 14700000 , 14800000 , 14900000 , 15000000 , 15100000 , 15200000 , 15300000 , 15400000 , 15500000 , 15600000 , 15700000 , 15800000 , 15900000 , 16000000 , 16100000 , 16200000 , 16300000 , 16400000 , 16500000 , 16600000 , 16700000 , 16800000 , 16900000 , 17000000 , 17100000 , 17200000 , 17300000 , 17400000 , 17500000 , 17600000 , 17700000 , 17800000 , 17900000 , 18000000 , 18100000 , 18200000 , 18300000 , 18400000 , 18500000 , 18600000 , 18700000 , 18800000 , 18900000 , 19000000 , 19100000 , 19200000 , 19300000 , 19400000 , 19500000 , 19600000 ,test_trt <- DMLtest.multiFactor(fit_add, coef = "treatmentExposed")
head(test_trt[order(test_trt$pvals), ])## chr pos stat pvals fdrs
## 3685034 NC_035781.1 51608334 -14.83936 8.152725e-50 1.604785e-42
## 9719204 NC_035784.1 48597326 14.48492 1.509023e-47 1.485183e-40
## 18755155 NC_035788.1 102536511 -13.03014 8.245804e-39 5.410357e-32
## 18755156 NC_035788.1 102536517 -12.24875 1.706032e-34 8.395400e-28
## 18755152 NC_035788.1 102536498 -12.21941 2.448535e-34 9.639410e-28
## 18755151 NC_035788.1 102536494 -12.18658 3.664794e-34 1.202299e-27# DMLs = FDR threshold (multiFactor gives per-site fdrs, not a methylKit-style % difference)
dml_trt <- test_trt[which(test_trt$fdrs < 0.05), ]
cat("Treatment DMLs (FDR < 0.05):", nrow(dml_trt), "\n")## Treatment DMLs (FDR < 0.05): 2237# Regions
dmr_trt <- callDMR(test_trt, p.threshold = 0.01, minlen = 50, minCG = 3, dis.merge = 100)
cat("Treatment DMRs:", nrow(dmr_trt), "\n")## Treatment DMRs: 18714.2 sex-dependence of the treatment effect
A null interaction at treatment-associated loci = signal present in both sexes. A significant interaction = the effect differs between male and female
test_intx <- DMLtest.multiFactor(fit_int, coef = "treatmentExposed:sexMale")
head(test_intx[order(test_intx$pvals), ])## chr pos stat pvals fdrs
## 2730352 NC_035781.1 19515763 13.61346 3.331080e-42 6.556911e-35
## 11456770 NC_035785.1 4756863 12.22096 2.402437e-34 2.364482e-27
## 15788416 NC_035787.1 63244810 11.70675 1.177044e-31 7.722994e-25
## 2579761 NC_035781.1 14358343 11.00214 3.731753e-28 1.836399e-21
## 18932450 NC_035789.1 4530142 10.22312 1.562321e-24 6.150558e-18
## 7306374 NC_035783.1 32042835 10.12476 4.290175e-24 1.407466e-17dml_intx <- test_intx[which(test_intx$fdrs < 0.05), ]
cat("Loci with sex-dependent treatment effect (FDR < 0.05):", nrow(dml_intx), "\n")## Loci with sex-dependent treatment effect (FDR < 0.05): 72175 Save outputs
save_bed <- function(df, path, score_col) {
# callDMR() returns NULL when no regions pass; callDML() can return 0 rows.
# Write an empty .bed in that case rather than erroring, so a knit completes.
if (is.null(df) || nrow(df) == 0) {
file.create(path)
message("save_bed: no rows for ", basename(path), " — wrote empty file.")
return(invisible(NULL))
}
# Region tables (callDMR) carry chr/start/end; single-base DML tables
# (DMLtest/callDML) carry chr/pos — derive a 0-based, half-open interval.
if (all(c("chr","start","end") %in% names(df))) {
start <- df$start; end <- df$end
} else if (all(c("chr","pos") %in% names(df))) {
start <- df$pos - 1L; end <- df$pos
} else {
stop("save_bed: df lacks chr/start/end or chr/pos for ", basename(path),
" - columns are: ", paste(names(df), collapse=", "))
}
score <- if (score_col %in% names(df)) df[[score_col]] else NA_real_
bed <- data.frame(chr = df$chr, start = start, end = end, score = score)
write.table(bed, path, quote = FALSE, sep = "\t", row.names = FALSE, col.names = FALSE)
}
# Site-level tables
write_tsv(test_trt, "../output/10.1-diff-methyl-DSS-parents/multifactor_treatment_test.tsv")
write_tsv(test_intx, "../output/10.1-diff-methyl-DSS-parents/multifactor_interaction_test.tsv")
# Region-level BEDs (multifactor DMR uses `areaStat` as region-level test statistic)
save_bed(dmr_trt, "../output/10.1-diff-methyl-DSS-parents/treatment_DMR.bed", "areaStat")6 Two-group smoothed test within each sex
Also want to try within-sex tests with DSS. For non-multifactor tets, DSS makes use of smoothing and dispersion shrinkage to essentially “borrow” info across neighboring CpG sites. This can be useful for handling low-coverage libraries.
It also produces effect estimates of the same type as methylKit (% methylation difference), which could be useful for more direct comparisons to the conventions used in Rondon et al. 2017 and Venkataraman et al. 2024. It would also be useful for questions involving directional concordance (hyper- v hypo-methylation)
WARNING: Both the sex-specific DMLtest() runs are very memory-intensive. I’ve needed to up the node request to ~300G just to successfully complete them, and it takes a while to finish.
run_sex <- function(sex_label) {
s_idx <- which(meta$sex == sex_label)
BS_s <- BSobj[, s_idx]
g_exp <- meta$sample[s_idx][meta$treatment[s_idx] == "Exposed"]
g_ctrl <- meta$sample[s_idx][meta$treatment[s_idx] == "Control"]
dml <- DMLtest(BS_s, group1 = g_exp, group2 = g_ctrl, smoothing = TRUE)
out <- list(
dml = callDML(dml, delta = 0.25, p.threshold = 0.01),
dmr = callDMR(dml, delta = 0.10, p.threshold = 0.01,
minlen = 50, minCG = 3, dis.merge = 100, pct.sig = 0.5)
)
# write per-sex results to disk immediately, so a later crash can't lose them
saveRDS(dml, file.path("../output/10.1-diff-methyl-DSS-parents", paste0(sex_label, "_DMLtest.rds")))
rm(dml, BS_s); gc()
out
}
res_fem <- run_sex("Female"); gc()## Smoothing ...
## Estimating dispersion for each CpG site, this will take a while ...
## Computing test statistics ...
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 10600975 566.2 30246922 1615.4 46708381 2494.5
## Vcells 1837819927 14021.5 5040215958 38453.8 7871336716 60053.6res_male <- run_sex("Male"); gc()## Smoothing ...
## Estimating dispersion for each CpG site, this will take a while ...
## Warning in mclapply(1:nrow(X2), foo, mc.cores = ncores): scheduled cores 158,
## 167, 169, 175, 183, 186, 187, 188, 189 did not deliver results, all values of
## the jobs will be affected
## Warning in shrk.phi[ix] <- shrk.phi2: number of items to replace is not a
## multiple of replacement length
## Computing test statistics ...
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 10601120 566.2 30246922 1615.4 46708381 2494.5
## Vcells 1838846063 14029.3 6058701112 46224.3 7871336716 60053.6cat("Female: DMLs =", nrow(res_fem$dml), " DMRs =", nrow(res_fem$dmr), "\n")## Female: DMLs = 10584 DMRs = 8489cat("Larvae: DMLs =", nrow(res_male$dml), " DMRs =", nrow(res_male$dmr), "\n")## Larvae: DMLs = 67461 DMRs = 227737 Save outputs
# Per-sex DML BEDs (diff = mu1 - mu2 = exposed - control)
save_bed(res_fem$dml, "../output/10.1-diff-methyl-DSS-parents/female_DML.bed", "diff")
save_bed(res_male$dml, "../output/10.1-diff-methyl-DSS-parents/male_DML.bed", "diff")
# Region-level BEDs (diff.Methy is exposed - control at the region level)
save_bed(res_fem$dmr, "../output/10.1-diff-methyl-DSS-parents/female_DMR.bed", "diff.Methy")
save_bed(res_male$dmr, "../output/10.1-diff-methyl-DSS-parents/male_DMR.bed", "diff.Methy")7.1 Directional concordance with parental gamete DMLs – UPDATED parental DMLs
Same check as in 10-diff-methyl-DSS, but this time using parental DMLs derived using DSS
zyg_dmls <- read_tsv("../output/10-diff-methyl-DSS/zygote_DML.bed",
col_names = c("chr", "start", "end", "zyg_diff"))## Rows: 7287 Columns: 4
## -- Column specification --------------------------------------------------------
## Delimiter: "\t"
## chr (1): chr
## dbl (3): start, end, zyg_diff
##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.lar_dmls <- read_tsv("../output/10-diff-methyl-DSS/larvae_DML.bed",
col_names = c("chr", "start", "end", "lar_diff"))## Rows: 5330 Columns: 4
## -- Column specification --------------------------------------------------------
## Delimiter: "\t"
## chr (1): chr
## dbl (3): start, end, lar_diff
##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.sperm <- read_tsv("../output/10.1-diff-methyl-DSS-parents/male_DML.bed",
col_names = c("chr", "start", "end", "par_diff"))## Rows: 67461 Columns: 4
## -- Column specification --------------------------------------------------------
## Delimiter: "\t"
## chr (1): chr
## dbl (3): start, end, par_diff
##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.egg <- read_tsv("../output/10.1-diff-methyl-DSS-parents/female_DML.bed",
col_names = c("chr", "start", "end", "par_diff"))## Rows: 10584 Columns: 4
## -- Column specification --------------------------------------------------------
## Delimiter: "\t"
## chr (1): chr
## dbl (3): start, end, par_diff
##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.shared_zyg_egg <- inner_join(zyg_dmls, egg, by = c("chr", "start", "end")) %>%
mutate(concordant = sign(zyg_diff) == sign(par_diff))
shared_lar_egg<- inner_join(lar_dmls, egg, by = c("chr", "start", "end")) %>%
mutate(concordant = sign(lar_diff) == sign(par_diff))
summarise(shared_zyg_egg, n_shared = n(), n_concordant = sum(concordant),
pct_concordant = mean(concordant) * 100)## # A tibble: 1 x 3
## n_shared n_concordant pct_concordant
## <int> <int> <dbl>
## 1 100 94 94summarise(shared_lar_egg, n_shared = n(), n_concordant = sum(concordant),
pct_concordant = mean(concordant) * 100)## # A tibble: 1 x 3
## n_shared n_concordant pct_concordant
## <int> <int> <dbl>
## 1 262 254 96.9shared_zyg_sperm <- inner_join(zyg_dmls, sperm, by = c("chr", "start", "end")) %>%
mutate(concordant = sign(zyg_diff) == sign(par_diff))
shared_lar_sperm <- inner_join(lar_dmls, sperm, by = c("chr", "start", "end")) %>%
mutate(concordant = sign(lar_diff) == sign(par_diff))
summarise(shared_zyg_sperm, n_shared = n(), n_concordant = sum(concordant),
pct_concordant = mean(concordant) * 100)## # A tibble: 1 x 3
## n_shared n_concordant pct_concordant
## <int> <int> <dbl>
## 1 1120 1113 99.4summarise(shared_lar_sperm, n_shared = n(), n_concordant = sum(concordant),
pct_concordant = mean(concordant) * 100)## # A tibble: 1 x 3
## n_shared n_concordant pct_concordant
## <int> <int> <dbl>
## 1 1022 1016 99.4write_tsv(dplyr::select(shared_zyg_egg, -concordant), "../output/10.1-diff-methyl-DSS-parents/shared_zygote_egg_DML.bed", col_names=TRUE)
write_tsv(dplyr::select(shared_lar_egg, -concordant), "../output/10.1-diff-methyl-DSS-parents/shared_larvae_egg_DML.bed", col_names=TRUE)
write_tsv(dplyr::select(shared_zyg_sperm, -concordant), "../output/10.1-diff-methyl-DSS-parents/shared_zygote_sperm_DML.bed", col_names=TRUE)
write_tsv(dplyr::select(shared_lar_sperm, -concordant), "../output/10.1-diff-methyl-DSS-parents/shared_larvae_sperm_DML.bed", col_names=TRUE)NExt I want to test whether parent-offspring DML overlap exceeds chance (“enriched”)