MEDIPS for chIP-Seq
0
0
Entering edit mode
Darek ▴ 10
@c356bc88
Last seen 3 days ago
Spain

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)
ChIPSeq • 133 views
ADD COMMENT
2
Entering edit mode

The paper links to their GitHub with all code used: https://github.com/Edert/DCS_predictions

ADD REPLY
1
Entering edit mode

Thank you very much! I must suffered from a selective blindness: it was in the article I thought I scanned thoroughly...

ADD REPLY

Login before adding your answer.

Traffic: 231 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6