Hello,

While there is a detailed, step by step procedure for using MEDIPS for the methylated DNA immunoprecipitation ("MEDIPS: genome-wide differential coverage analysis of sequencing data derived from DNA enrichment experiments.") I am so far unable to find a similar recommended and tested procedure for the ChIP-Seq data.

For example the recent review: Eder, T., Grebien, F. Comprehensive assessment of differential ChIP-seq tools guides optimal algorithm selection. Genome Biol 23, 119 (2022). https://doi.org/10.1186/s13059-022-02686-y does list MEDIPS as one of the best tools for ChIP-seq broad peak detection but AFAIK there is no code used to call the peaks.

Any hints as where such procedures are described will be greatly appreciated.

Many thanks for your help in advance

Darek Kedra

Update:

I was able to transform part of the bash shell script: https://github.com/Edert/DCS_predictions/blob/main/run_pi_medips.sh

to the R script below using some fastq files generated by Eder et al. Good news: it does run and produces the output. Not so good news: A lot of peaks/differentially covered regions (31k) for just one mouse chromosome. I guess I still miss some steps to get proper broad peaks.

#!/usr/bin/env Rscript

library(MEDIPS)
library(BSgenome.Mmusculus.UCSC.mm10)


get_medips_set <- function(file_name) {
    BSgenome='BSgenome.Mmusculus.UCSC.mm10'
    uniq <- 1
    extend <- 300
    shift <- 0
    ws <- 100
    chr.select <- 'chr19'

    result = MEDIPS.createSet(file = file_name, 
        BSgenome = BSgenome, 
        extend = extend, 
        shift = shift, 
        uniq = uniq, 
        window_size = ws, 
        chr.select = chr.select)

return(result)
}


fn_s1_inpt <- "mm10_chr19_broad_100u0d_sample1-INPUT.mm10.bwamem2.markdup.bam" 
fn_s1_rep1 <- "mm10_chr19_broad_100u0d_sample1-rep1.mm10.bwamem2.markdup.bam"
fn_s1_rep2 <- "mm10_chr19_broad_100u0d_sample1-rep2.mm10.bwamem2.markdup.bam"

fn_s2_inpt <- "mm10_chr19_broad_100u0d_sample2-INPUT.mm10.bwamem2.markdup.bam"
fn_s2_rep1 <- "mm10_chr19_broad_100u0d_sample2-rep1.mm10.bwamem2.markdup.bam"
fn_s2_rep2 <- "mm10_chr19_broad_100u0d_sample2-rep2.mm10.bwamem2.markdup.bam"

s1_inpt <- get_medips_set(fn_s1_inpt)
s1_rep1 <- get_medips_set(fn_s1_rep1)
s1_rep2 <- get_medips_set(fn_s1_rep2)

set_1 <- c(s1_rep1, s1_rep2)

s2_inpt <- get_medips_set(fn_s2_inpt)
s2_rep1 <- get_medips_set(fn_s2_rep1)
s2_rep2 <- get_medips_set(fn_s2_rep2)

set_2 <- c(s2_rep1, s2_rep2)


CS = MEDIPS.couplingVector(pattern = 'CG', refObj = set_1[[1]])

print("after CS")

resultTable = MEDIPS.meth(MSet1 = set_1, 
MSet2 = set_2, 
CSet = CS, 
ISet1 = s1_inpt, 
ISet2 = s2_inpt, 
chr = 'chr19', 
p.adj='BH', 
diff.method='edgeR', 
CNV=FALSE, 
MeDIP=FALSE, 
diffnorm='tmm')

print("after resultTable")

sig = MEDIPS.selectSig(results=resultTable, 
p.value=0.05, 
adj=TRUE, 
ratio=NULL, 
bg.counts=NULL, 
CNV=FALSE)

print("after selectSig")

options(scipen = 999)
write.table(file='results_medips_tmm_ws100.txt', 
sig, 
quote = FALSE, 
sep = "\t",
row.names = FALSE, 
col.names = TRUE)