Hello Fellas we are running a set of MeDip-seq data using the MeQa pipeline. the pipeline breaks when it tries to read the BAM file for MEDIPS analysis and give this error
Error: could not find function "MEDIPS.readAlignedSequences" Execution halted
i did some reading material and it seems that the MEDIPS.readAlignedSequences function in the MEDIPS package has been deprecated. and we should use the following function instead
MEDIPS.createSet(file=NULL, extend=0, shift=0, window_size=300, BSgenome=NULL, uniq=1e-3, chr.select=NULL, paired = F, sample_name=NULL, bwa=FALSE.
my questions is if we change the function and set the parameters for it in the config file that runs the pipeline it is gonna be a problem for the other steps of the pipeline
here is the part of the python script we use and the part of the config file containing the parameters for the function
###For treat bed file
fw.write('treatBed="%s"\n' % (treatBed))
**fw.write('TREAT.SET=MEDIPS.readAlignedSequences(BSgenome = \"%s\",file = treatBed, numrows = %s)\n' % (BS,mapcount(treatBed)))**
fw.write('TREAT.SET = MEDIPS.genomeVector(data = TREAT.SET, bin_size = 50,extend = 400)\n')
fw.write('TREAT.SET = MEDIPS.getPositions(data = TREAT.SET, pattern = \"CG\")\n')
fw.write('TREAT.SET = MEDIPS.couplingVector(data = TREAT.SET, fragmentLength = 700,func = \"count\")\n')
fw.write('TREAT.SET = MEDIPS.calibrationCurve(data = TREAT.SET)\n')
fw.write('TREAT.SET = MEDIPS.normalize(data = TREAT.SET)\n')
###For control bed file
fw.write('controlBed="%s"\n' % (controlBed))
**fw.write('CONTROL.SET=MEDIPS.readAlignedSequences(BSgenome = \"%s\",file = controlBed, numrows = %s)\n' % (BS,mapcount(controlBed))**)
fw.write('CONTROL.SET = MEDIPS.genomeVector(data = CONTROL.SET, bin_size = 50,extend = 400)\n')
fw.write('CONTROL.SET = MEDIPS.getPositions(data = CONTROL.SET, pattern = \"CG\")\n')
fw.write('CONTROL.SET = MEDIPS.couplingVector(data = CONTROL.SET, fragmentLength = 700,func = \"count\")\n')
fw.write('CONTROL.SET = MEDIPS.calibrationCurve(data = CONTROL.SET)\n')
fw.write('CONTROL.SET = MEDIPS.normalize(data = CONTROL.SET)\n')
###For control bed file
if inputBed and len(inputBed)>1:
fw.write('inputBed="%s"\n' % (inputBed))
**fw.write('INPUT.SET=MEDIPS.readAlignedSequences(BSgenome = \"%s\",file = inputBed, numrows = %s)\n' % (BS,mapcount(inputBed)))**
fw.write('INPUT.SET = MEDIPS.genomeVector(data = INPUT.SET, bin_size = 50,extend = 400)\n')
fw.write('INPUT.SET = MEDIPS.getPositions(data = INPUT.SET, pattern = \"CG\")\n')
fw.write('INPUT.SET = MEDIPS.couplingVector(data = INPUT.SET, fragmentLength = 700,func = \"count\")\n')
fw.write('INPUT.SET = MEDIPS.calibrationCurve(data = INPUT.SET)\n')
fw.write('INPUT.SET = MEDIPS.normalize(data = INPUT.SET)\n')
if inputBed and len(inputBed)>1:
fw.write('diff.meth = MEDIPS.methylProfiling(data1 = CONTROL.SET, data2 = TREAT.SET, input = INPUT.SET, frame_size = 500, select = 2)\n')
else:
fw.write('diff.meth = MEDIPS.methylProfiling(data1 = CONTROL.SET, data2 = TREAT.SET, frame_size = 500, select = 2)\n')
fw.write('write.table(diff.meth, file = "%s%s_diff.meth.txt", sep =\"\t\", quote = F, col.names = T, row.names = F)\n' % (dout,name))
fw.close()
runDMR(fo,dout+name)
def runDMR(filename,name):
cmd="R CMD BATCH " + filename +" "+name+"-DMR.log"
info(cmd)
info("This will take takes a comparable long processing time.")
try:
p = subprocess.Popen(['/bin/bash', '-c', cmd])
sts = os.waitpid(p.pid, 0)
info('#Congratulations! Dierentially methylated regions have been output!')
except:
info('#Congratulations! Run %s using R for the dierentially methylated regions! DMR could not run R directly.' % (filename))
info('Please run %s in your command line.', cmd)
here is the config file
#Here set the parameters for MEDIPS
#Full parameters are described in: Computational analysis of genome-wide DNA methylation during the differentiation of human embryonic stem cells along the endodermal lineage, Genome Res. 2010 Oct;20(10):1441-50. Epub 2010 Aug 27 (PMID: 20802089).
#http://bioconductor.org/packages/2.9/bioc/vignettes/MEDIPS/inst/doc/MEDIPS.pdf
[Calibration]
#[Calibration]:Calibration Curve and Linear Regression
#chr should set as chr=
chr=
xrange=50
[methylProfiling]
#chr should set as chr=
chr=
select=2
transf=
math=mean
frame_size=500
step=250
ROI_file=
[genomeVector]
bin_size=50
extend=300
[couplingVector]
fragmentLength=700
func=count
[Output]
#The result output directory
OutDirectory=.
thank you , any help is valuable