Kathleen Durkin 2026-09-04
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).
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)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 -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 ...## 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")
}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 ,## [1] "(Intercept)" "treatmentExposed"
## [3] "sexMale" "treatmentExposed:sexMale"Fit an additive model (no interaction) to test for signal of parental treatment effect while controlling for 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: 1871A 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): 7217save_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")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.6## 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## Female: DMLs = 10584 DMRs = 8489## Larvae: DMLs = 67461 DMRs = 22773# 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")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”)