10.1-diff-methyl-DSS-parents ================ Kathleen Durkin 2026-09-04 - [1 Setup](#1-setup) - [2 Download methylation calls](#2-download-methylation-calls) - [3 Read Bismark coverage into a BSseq object](#3-read-bismark-coverage-into-a-bsseq-object) - [4 Multi-factor model (design-level tests)](#4-multi-factor-model-design-level-tests) - [4.1 Parental treatment effect](#41-parental-treatment-effect) - [4.2 sex-dependence of the treatment effect](#42-sex-dependence-of-the-treatment-effect) - [5 Save outputs](#5-save-outputs) - [6 Two-group smoothed test within each sex](#6-two-group-smoothed-test-within-each-sex) - [7 Save outputs](#7-save-outputs) - [7.1 Directional concordance with parental gamete DMLs – UPDATED parental DMLs](#71-directional-concordance-with-parental-gamete-dmls--updated-parental-dmls) Differential methylation analysis using the `DSS` package (Dispersion Shrinkage for Sequencing data). This is essentially an identical copy of `10-diff-methyl-DSS.Rmd`, BUT it will be applied to the data from parental gametes (analyzed in more depth in Venkataraman et al. 2024). # 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. ``` bash #!/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\")"' ``` ``` r # 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: 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. ``` bash # 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\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. ``` bash 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. ``` r 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 | ``` r 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 ... ``` r 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 ``` r # 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: 19684036 EDIT: 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: ``` r # 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) ``` r 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 , ``` r # 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*. ``` r 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 , ``` r 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 ``` r # 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 ``` r # 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: 1871 ## 4.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 ``` r 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-17 ``` r dml_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): 7217 # 5 Save outputs ``` r 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. ``` r 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.6 ``` r res_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.6 ``` r cat("Female: DMLs =", nrow(res_fem$dml), " DMRs =", nrow(res_fem$dmr), "\n") ``` ## Female: DMLs = 10584 DMRs = 8489 ``` r cat("Larvae: DMLs =", nrow(res_male$dml), " DMRs =", nrow(res_male$dmr), "\n") ``` ## Larvae: DMLs = 67461 DMRs = 22773 # 7 Save outputs ``` r # 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 ``` r 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. ``` r 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. ``` r 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. ``` r 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. ``` r 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 ## ## 1 100 94 94 ``` r summarise(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 ## ## 1 262 254 96.9 ``` r shared_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 ## ## 1 1120 1113 99.4 ``` r summarise(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 ## ## 1 1022 1016 99.4 ``` r write_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”)