WGBS_pipe_lite.py~

__version__="v1.2.1"


import os
import sys
import time
import argparse

#to-do: reduce the imports, take too long
import zipfile
import re
import subprocess
import pandas
import numpy
import commands
import math
import time
from collections import OrderedDict
import logging
from operator import is_not
from functools import partial
import statistics
import pysam
import fnmatch

sys.path.insert(0, "/data/manke/repository/scripts/DNA_methylation/ruffus")
from ruffus import *
from ruffus.proxy_logger import *
from ruffus.combinatorics import *
from ruffus.drmaa_wrapper import run_job, run_job_using_drmaa, error_drmaa_job

from PIL import Image
import string

parser = argparse.ArgumentParser(prog="methPipe", version="1.2.0", description="runs complete CpG methylation pipeline", formatter_class=argparse.RawTextHelpFormatter)
parser.add_argument("-ri", "--readIn", dest="readdir", action="store", default=False, help="input read folder")
parser.add_argument("-qi", "--fqcIn", dest="fqcdir", action="store", default=False, help="folder with fastqc.zip results for raw reads")
parser.add_argument("-g", "--ref", dest="refpath", action="store", default=False, help="path to indexed reference genome")
parser.add_argument("-cg", "--cref", dest="convrefpath", action="store", default=False, help="path to converted and indexed reference genome")
parser.add_argument("-il", "--intList", dest="intList", action="append", default=None, help="target interval file(s)")
parser.add_argument("-bl", "--blackList", dest="blackList", action="store", default=None, help="SNP black list")
parser.add_argument("-si", "--sampleInfo", dest="sampleInfo", action="store", default=None, help="sample sheet")
parser.add_argument("-wd", "--wdir", dest="wdir", action="store", default=False, help="output folder")
parser.add_argument("-vb", "--verbose", dest="verbose", action="store_true", help="more detailed log")
parser.add_argument("-bs", "--batchSize", dest="bsize", action="store",metavar="INT",type=int,default=10, help="number of samples to process in parallel")
parser.add_argument("-nt", "--numThr", dest="nthreads", action="store",metavar="INT",type=int,default=8, help="number of threads to use per sample")
parser.add_argument("-tR", "--trimReads", dest="trimReads", action="store",default=None,choices=['auto','user',None], help="adapter-trim and hardclip the reads")
parser.add_argument("--nextera", action="store_true", help="Trim Nextera adapters, rather than the default TruSeq adapters. Requires --trimReads")
parser.add_argument("--trimThreshold", default=10, type=int, help="If --trimReads is specified, this sets the phred score threshold.")
parser.add_argument("--trimOtherArgs", default="", help="Other arguments you would like passed to cutadapt.")
parser.add_argument("-cr", "--convRef", dest="convRef", action="store_true", help="BS-convert reference genome")
parser.add_argument('--mbias', dest="mbias_ignore",action="store",default="auto",help="number of nucleotides with mbias to ignore during methylation extraction")
parser.add_argument("--touchOnly", dest="touchOnly", action="store_true", help="only touch files")
parser.add_argument("--target_tasks", dest="target_tasks", action="append",default=[], help="target tasks")
parser.add_argument("--forcedtorun_tasks", dest="forcedtorun_tasks", action="append",default=[], help="forced to run tasks")

args = parser.parse_args()
pipev=__version__

#setup central working directory
wdir=args.wdir
if not os.path.exists(wdir):
    os.makedirs(wdir)
    
#setup logging
logger = logging.getLogger(__name__)
fhandler = logging.FileHandler(filename=os.path.join(wdir,'pipeline.log'), mode='a')
formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')
fhandler.setFormatter(formatter)
logger.addHandler(fhandler)
logger.setLevel(logging.DEBUG)

logger.debug(subprocess.check_output('echo $DRMAA_LIBRARY_PATH',shell=True))

import drmaa

#initiate 1 drmaa session for the whole pipeline
mySession=drmaa.Session()
mySession.initialize()

logger.info('Pipeline version is '+pipev)
logger.info('Working directory is : ' + wdir)
logger.info(args)

readdir=args.readdir
os.chdir(readdir)

libset2 = []
# Walk through directory
for dName, sdName, fList in os.walk(readdir):
    for fileName in fList:
        if fnmatch.fnmatch(fileName, "*fastq.gz"): # Match search string
            libset2.append(os.path.join(dName, fileName))

libset2_R1=filter(lambda x:'_R1.fastq.gz' in x, libset2)
libset2_R1.sort()
libset2_R2=filter(lambda x:'_R2.fastq.gz' in x, libset2)
libset2_R2.sort()
read_root=[ re.sub('_R1.fastq.gz','',R1f) for R1f in libset2_R1 ]
INfiles=list(zip(libset2_R1,libset2_R2))
    
logger.debug(INfiles)    

######################## paths to unconverted reference genome #######################################
if os.path.isfile(args.refpath):
    refG=args.refpath
    refGpath=os.path.dirname(refG)
    
elif os.path.isfile(os.path.join('/data/repository/organisms',(args.refpath+'_ensembl'),'genome_fasta/genome.fa')):
    refG=os.path.join('/data/repository/organisms',(args.refpath+'_ensembl'),'genome_fasta/genome.fa')
    refGpath=os.path.dirname(refG)
    
else:
    logger.info('Reference genome not recognized. Please check spelling and/or path.')
        
######################## paths to converted reference genome #######################################

if os.path.isfile(str(args.convrefpath)):
    crefG=args.convrefpath
    crefGpath=os.path.dirname(crefG)
elif os.path.isfile(os.path.join('/data/repository/organisms',(str(args.refpath)+'_ensembl'),'BWAmethIndex','genome.fa')):
    
    crefG=os.path.join('/data/repository/organisms',(args.refpath+'_ensembl'),'BWAmethIndex','genome.fa')
    crefGpath=os.path.dirname(crefG)
elif not args.convrefpath:
    print('Converted reference genome not specified, will be generated by the pipeline.')
    args.convRef=True
    crefG='not available.'
else:
    logger.error('Converted reference genome not recognized. Please check spelling and/or path.')    
    
logger.debug('Reference genome is '+refG)
logger.debug('Converted reference genome is '+crefG)
logger.debug('Convert genome? ' + str(args.convRef))

##################PATHS TO EXECUTABLES###############################################################
FQCpath='/package/FastQC-0.11.3'
#cutpath='/package/cutadapt-1.9.1/bin'
cutpath='module load cutadapt/1.16 ;'
mCTpath='/data/manke/repository/scripts/DNA_methylation/methylCtools'
tabpath='/package/tabix-1.2.1'
bwapath='/package/bwa-0.7.4/bin'
#bmethpath='/data/manke/repository/scripts/DNA_methylation/bwa-meth-master_2016'
bmethpath='module load bwameth; '
sampath='/package/samtools-1.3/bin'
#Picpath='/package/picard-tools-1.136'
sambambapath='module load sambamba/0.6.6;'
GATKpath='/package/GenomeAnalysisTK-3.5'
POMpath='/package/MethylDackel-0.3.0/bin'
bedpath='module load bedtools2/2.27.0;'
metipath='/data/manke/repository/scripts/DNA_methylation/metilene_v0.2-6'
Rpath='/package/R-3.3.1/bin'
Rlib='/data/manke/repository/scripts/DNA_methylation/Rlibs.3.3.1'

######################################################################################################

#add folder with custom modules to path            
#sys.path.append('/data/manke/repository/scripts/DNA_methylation/WGBS_pipe/v0.02')
sys.path.append(os.path.dirname(os.path.realpath(__file__)))


#############################PIPELINE#################################################################################
#TRIM READS############################################################################################################
if not args.trimReads is None :
    fqcout=os.path.join(wdir,'fastqc_cut')
    cutout=os.path.join(wdir,'reads_cut')
    #setup shared logging 
    log_args={}
    log_args["file_name"] = os.path.join(fqcout,"LOG")
    log_args["formatter"] = "%(asctime)s - %(name)s - %(levelname)6s - %(message)s"
    log_args["delay"]     = True
    log_args["level"]     = logging.DEBUG
    (shared_logger,logging_mutex) = make_shared_logger_and_proxy (setup_std_shared_logger,
                                               "shared_logger",
                                               log_args)  
    if args.fqcdir and args.trimReads == "auto":    
        fqcdir=args.fqcdir
        fqcL=[ os.path.join(fqcdir,z ) for z in os.listdir(fqcdir) ]
        zipR1=filter(lambda x:'_R1_fastqc.zip' in x, fqcL)
        zipR1.sort()
        zipR2=filter(lambda x:'_R2_fastqc.zip' in x, fqcL)
        zipR2.sort()
        zipL=list(zip(zipR1,zipR2))
        zipL2=[ list(z) for z in zipL ]
        os.chdir(fqcdir)
        #from BSprecut import calc_cutThd
        @mkdir(fqcout)
        @transform(zipL2,suffix('_R1_fastqc.zip'),'.R12.ct.txt',output_dir=fqcout)
        def get_cutThd(input_files,output_file):
            ii1 = input_files[0]
            ii2 = input_files[1]
            from BSprecut import calc_cutThd
            with open(output_file, "w") as oo:
                cutThdRes=calc_cutThd([ii1,ii2],args.fqcdir,fqcout,shared_logger,logging_mutex)
                oo.write('\n'.join('%s\t%s\n' % x for x in cutThdRes))
                    
        
    os.chdir(readdir)
    
if args.trimReads == 'auto':
    @follows(get_cutThd)
    @mkdir(cutout,os.path.join(cutout,'logs'))
    @transform(input=INfiles,filter=suffix('_R1.fastq.gz'),output=['_R1.fastq.gz','_R2.fastq.gz'],output_dir=cutout) 
    def trim_reads(input_files,output_files):
        ii1 = input_files[0]
        ii2 = input_files[1]
        ii3 = os.path.join(fqcout,re.sub('_R1.fastq.gz','.R12.ct.txt',os.path.basename(ii1)))
        print(ii3)
        oo1 = output_files[0]
        oo2 = output_files[1]
        #prepare threshold values
        ctfile = open(ii3)
        for line in ctfile:
            fields = line.strip().split()
            ct1=fields[0]
            ct2=fields[1]
    
        from BSprecut import cut_reads_auto
        cut_reads_auto(ii1,ii2,oo1,oo2,ct1,ct2,cutpath,mySession,cutout,logger,args)

elif args.trimReads=='user':
    @mkdir(cutout,os.path.join(cutout,'logs'))
    @transform(input=INfiles,filter=suffix('_R1.fastq.gz'),output=['_R1.fastq.gz','_R2.fastq.gz'],output_dir=cutout) 
    def trim_reads(input_files,output_files):
        ii1 = input_files[0]
        ii2 = input_files[1]
        oo1 = output_files[0]
        oo2 = output_files[1]
        from BSprecut import cut_reads_user
        cut_reads_user(ii1,ii2,oo1,oo2,cutpath,mySession,cutout,logger,args)

    @mkdir(fqcout,os.path.join(fqcout,'logs'))
    @transform(input=trim_reads,filter=suffix('_R1.fastq.gz'),output=['_R1.zip','_R1.html','_R2.zip','_R2.html'],output_dir=fqcout)
    def postTrim_fastqc(input_files,output_files):
        ii1 = input_files[0]
        ii2 = input_files[1]
        from BSprecut import post_trim_fqc
        post_trim_fqc(ii1,ii2,fqcout,FQCpath,mySession,logger)   
        
#CONVERT REFERENCE GENOME##########################################################################################
if args.convRef:
    crefGpath=os.path.join(wdir,'conv_ref')
    if not os.path.exists(crefGpath):
        os.makedirs(crefGpath)
    
@active_if(args.convRef)
@transform(refG,suffix('.fa'),'.fa.bwameth.c2t',output_dir=crefGpath)
def conv_ref(input_file,output_files):
    os.chdir(refGpath)
    ii = input_file
    from BSrefprep import bmeth_refprep
    bmeth_refprep(bmethpath,ii,crefGpath,logger)
    os.chdir(wdir)
    
#CONVERT AND MAP READS############################################################################################
bamoutO=os.path.join(wdir,'bams')
###################################################################################################################  
if args.convRef:
    if args.trimReads:
        @follows("conv_ref")
        @mkdir(bamoutO,os.path.join(bamoutO,'logs'))
        @transform(trim_reads,suffix('_R1.fastq.gz'),'.sorted.bam',output_dir=bamoutO)
        def map_reads(input_files,output_file):
            ii1 = input_files[0]
            ii2 = input_files[1]
            oo = output_file
            crefG=os.path.join(crefGpath,os.path.basename(refG))               
            from BSmapWGBS import bMeth_map_reads
            bMeth_map_reads(ii1,ii2,oo,bmethpath,sampath,bamoutO,crefG,args.nthreads,mySession,logger)
    else:
        @follows("conv_ref")
        @mkdir(bamoutO,os.path.join(bamoutO,'logs'))
        @transform(INfiles,suffix('_R1.fastq.gz'),'.sorted.bam',output_dir=bamoutO)
        def map_reads(input_files,output_file):
            ii1 = input_files[0]
            ii2 = input_files[1]
            oo = output_file
            crefG=os.path.join(crefGpath,os.path.basename(refG))
            from BSmapWGBS import bMeth_map_reads
            bMeth_map_reads(ii1,ii2,oo,bmethpath,sampath,bamoutO,crefG,args.nthreads,mySession,logger)
else:
    if args.trimReads:
        @mkdir(bamoutO,os.path.join(bamoutO,'logs'))
        @transform(trim_reads,suffix('_R1.fastq.gz'),'.sorted.bam',output_dir=bamoutO)
        def map_reads(input_files,output_file):
            ii1 = input_files[0]
            ii2 = input_files[1]
            oo = output_file
            from BSmapWGBS import bMeth_map_reads
            bMeth_map_reads(ii1,ii2,oo,bmethpath,sampath,bamoutO,crefG,args.nthreads,mySession,logger)
    else:
        @mkdir(bamoutO,os.path.join(bamoutO,'logs'))
        @transform(INfiles,suffix('_R1.fastq.gz'),'.sorted.bam',output_dir=bamoutO)
        def map_reads(input_files,output_file):
            ii1 = input_files[0]
            ii2 = input_files[1]
            oo = output_file
            from BSmapWGBS import bMeth_map_reads
            bMeth_map_reads(ii1,ii2,oo,bmethpath,sampath,bamoutO,crefG,args.nthreads,mySession,logger) 

################## BAM POSTPROCESSING #########################################################
#bam postprocessing; bwa-meth is sorted and indexed, Bismark and methylCtools require these steps

@transform(map_reads,suffix('.sorted.bam'),'.sorted.bam.bai') 
def index_bam(input_file,output_file):
    ii = input_file
    oo = output_file  
    from BSmapWGBS import BS_index_bam
    BS_index_bam(ii,sampath,bamoutO,mySession,logger)
        
@transform(index_bam,suffix('.sorted.bam.bai'),'.PCRrm.bam')
def PCRdup_rm(input_file,output_file):
    ii=input_file
    oo=output_file
    from BSmapWGBS import BS_rm_dupes
    BS_rm_dupes(ii,oo,sambambapath,bamoutO,args.nthreads,mySession,logger)

######################################################################################################################  
#RUN VARIOUS COVERAGE AND QUALITY METRICS#######################################################################
metout=os.path.join(wdir,'QC_metrics')
auxdir=os.path.join(wdir,'aux_files')
reads_downsampled=os.path.join(wdir,'downsampled_reads')

#prepare bed file with random CpGs
@mkdir(auxdir,os.path.join(auxdir,'logs'))
@transform(refG,suffix('.fa'),'.poz.ran1M.sorted.bed',output_dir=auxdir)
def get_ran_CG(input_file,output_file):
    ii=input_file
    oo=output_file
    from BSmetricsWGBS import mCT_get_ranCG
    mCT_get_ranCG(ii,oo,refG,auxdir,mCTpath,logger)
    

#methylation bias
@mkdir(metout,os.path.join(metout,'logs'))
@transform(PCRdup_rm,suffix('.PCRrm.bam'),'.Mbias.txt',output_dir=metout)
def calc_Mbias(input_file,output_file):
    ii=input_file
    oo=output_file
    oos=re.sub('.txt','',oo)
    from BSmetricsWGBS import BS_Mbias
    BS_Mbias(ii,oos,POMpath,refG,metout,args.nthreads,mySession,logger)    


#GATK depth of coverage
if ( args.intList ): 
    int_dest=[re.sub('.bed','.mean.doc.sample_summary',os.path.basename(x)) for x in args.intList]
    int_dest[:0]=['mean.CG.doc.sample_summary']
    int_dest[:0]=['mean.genome.doc.sample_summary']
    @mkdir(metout,os.path.join(metout,'logs'))
    @follows(get_ran_CG)
    @transform(PCRdup_rm,suffix('PCRrm.bam'),int_dest,output_dir=metout)#
    def depth_of_cov(input_file,output_files):
        ii=input_file
        oos=output_files 
        oos2=[w.replace('.sample_summary', '') for w in oos]
        from BSmetricsWGBS import BS_doc_XT
        BS_doc_XT(ii,oos2,args.intList,auxdir,refG,GATKpath,metout,mySession,logger)
else:
    @mkdir(metout,os.path.join(metout,'logs'))
    @follows(get_ran_CG)
    @transform(PCRdup_rm,suffix('.PCRrm.bam'),['.mean.genome.doc.sample_summary','.mean.CG.doc.sample_summary'],output_dir=metout)#,
    def depth_of_cov(input_file,output_files):
        ii=input_file
        oos=output_files
        oos2=[w.replace('.sample_summary', '') for w in oos]
        from BSmetricsWGBS import BS_doc
        BS_doc(ii,oos2,refG,auxdir,GATKpath,metout,mySession,logger)

        
#conversion rate: currently from fastq, implement phiX control!
if  args.trimReads :
    @mkdir(reads_downsampled,os.path.join(reads_downsampled,'logs'))
    @transform(trim_reads,suffix('_R1.fastq.gz'),['_R1.fastq.gz','_R2.fastq.gz'],output_dir=reads_downsampled) #down to 5mln reads
    def dsmpl_reads(input_files,output_files):
        ii1=input_files[0]
        ii2=input_files[1]
        oo1=output_files[0]
        oo2=output_files[1]
        from BSmetricsWGBS import BS_downsample_reads
        BS_downsample_reads('5000000',args.nthreads,ii1,oo1,ii2,oo2,pipev,reads_downsampled,mySession,logger)

    @mkdir(metout,os.path.join(metout,'logs'))
    @transform(dsmpl_reads,suffix('_R1.fastq.gz'),'.conv.rate.txt',output_dir=metout) #it's ok to have both reads summarized in 1 file
    def conv_rate(input_files,output_file):
        ii1=input_files[0]
        ii1sub=re.sub('_R1.fastq.gz','',ii1)
        oo=output_file
        from BSmetricsWGBS import BS_conv_rate
        BS_conv_rate(ii1sub,oo,metout,mySession,logger)
else: 
    @mkdir(reads_downsampled,os.path.join(reads_downsampled,'logs'))
    @transform(IN_files,suffix('_R1.fastq.gz'),['_R1.fastq.gz','_R2.fastq.gz'],output_dir=reads_downsampled) #down to 5mln reads
    def dsmpl_reads(input_files,output_files):
        ii1=input_files[0]
        ii2=input_files[1]
        oo1=output_files[0]
        oo2=output_files[1]
        from BSmetricsWGBS import BS_downsample_reads
        BS_downsample_reads('5000000',args.nthreads,ii1,oo1,ii2,oo2,pipev,reads_downsampled,mySession,logger)
    @mkdir(metout,os.path.join(metout,'logs'))
    @transform(dsmpl_reads,suffix('_R1.fastq.gz'),'.conv.rate.txt',output_dir=metout) #it's ok to have both reads summarized in 1 file
    def conv_rate(input_files,output_file):
        ii1=input_files[0][0]
        ii1sub=re.sub('_R1.fastq.gz','',ii1)
        oo=output_file
        from BSmetricsWGBS import BS_conv_rate
        BS_conv_rate(ii1sub,oo,metout,mySession,logger)
        
#flagstat -> mapping rate

@mkdir(metout,os.path.join(metout,'logs'))
@transform(PCRdup_rm,suffix('.PCRrm.bam'),'.flagstat',output_dir=metout)
def get_flagstat(input_file,output_file):
    ii=input_file
    oo=output_file
    from BSmetricsWGBS import BS_flagstat
    BS_flagstat(ii,oo,sampath,metout,mySession,logger)   
    
    
################################################################################################################## 
###final QC report: only if started from reads or bam files
if not os.path.exists(metout):
        os.makedirs(metout)
os.chdir(metout)
@merge(input=[depth_of_cov,conv_rate,calc_Mbias,get_flagstat],output='QC_report.pdf')
def produce_report(input_file,output_file):
    ii=input_file
    oo=output_file
    from BSmetricsWGBS import BS_QC_rep
    BS_QC_rep(Rpath,metout,logger)
    
#####EXTRACT METHYLATION COUNTS  ########################################################################
mextout=os.path.join(wdir,'methXT')
#auxdir=os.path.join(wdir,'aux_files') #defined previously
@follows('calc_Mbias',mkdir(mextout),mkdir(os.path.join(mextout,'logs')))
@transform(PCRdup_rm,suffix('.PCRrm.bam'),'_CpG.bedGraph',output_dir=mextout)
def methyl_extract(input_file,output_file):
    ii = input_file
    oo = output_file
    oos=re.sub('_CpG.bedGraph','',oo)
    from BSmethXT_WGBS import methXT_POM
    methXT_POM(ii,metout,oos,refG,POMpath,mextout,args.mbias_ignore,args.nthreads,mySession,logger)
@transform(methyl_extract,suffix('_CpG.bedGraph'),'.CpG.filt2.bed',output_dir=mextout) 
def CpG_filt(input_file,output_file):
    ii = input_file
    oo = output_file
    from BSmethXT_WGBS import filt_POM
    filt_POM(ii,bedpath,mextout,Rpath,Rlib,pipev,mySession,logger,args.blackList)    
#########################################################################################################    

     
####RUN SINGLE CYTOSINE STATISTICS IN R VIA LIMMA#########################################################
if args.sampleInfo :
    CpGstat_out=os.path.join(wdir,'singleCpG_stats_limma')
    if not os.path.exists(CpGstat_out) or not os.path.exists(os.path.join(CpGstat_out,'logs')):
        os.makedirs(CpGstat_out)
        os.makedirs(os.path.join(CpGstat_out,'logs'))
    os.chdir(CpGstat_out)
    @merge(input=CpG_filt,output=[os.path.join(CpGstat_out,'singleCpG.RData'),os.path.join(CpGstat_out,'limdat.LG.RData'),os.path.join(CpGstat_out,'metilene.IN.txt')])
    def CpG_stats(input_files,output_file):
        ii = input_files[0]
        oo = output_file
        from BSstats_WGBS import single_CpG_limma
        single_CpG_limma(os.path.dirname(ii),args.sampleInfo,CpGstat_out,Rpath,Rlib,pipev,mySession,logger)
########################################################################################################
        
    
####RUN INTERVAL AGGREGATE STATISTICS IN R ####################################################
###prepare bed files with all reference CpG positions within the genomic intervals
if (args.intList):
    int_dest=[re.sub('.bed','.CpGlist.bed',os.path.basename(x)) for x in args.intList]
    @follows(get_ran_CG)
    @transform(refG,suffix('.fa'),int_dest,output_dir=auxdir)
    def get_CG_per_int(input_file,output_files):
        ii=input_file
        oo=output_files
        from BSstats_WGBS import mCT_get_CpGxInt
        mCT_get_CpGxInt(ii,oo,refG,args.intList,auxdir,bedpath,mySession,logger)


if ( args.intList and args.sampleInfo ):
    int_dest=[re.sub('.bed','.aggCpG.RData',os.path.basename(x)) for x in args.intList]
    intStat_out=os.path.join(wdir,'aggregate_stats_limma')
    @mkdir(intStat_out,os.path.join(intStat_out,'logs'))
    @follows(get_CG_per_int)
    @transform(CpG_stats,suffix('.RData'),int_dest,output_dir=intStat_out)
    def intAgg_stats(input_files,output_files):
        ii = os.path.join(CpGstat_out,input_files[1])
        oo = output_files
        auxList=[os.path.join(auxdir,re.sub('.fa',re.sub('.bed','.CpGlist.bed',os.path.basename(x)),os.path.basename(refG))) for x in args.intList]
        from BSstats_WGBS import int_stats_limma
        int_stats_limma(ii,args.intList,auxList,args.sampleInfo,intStat_out,Rpath,Rlib,pipev,mySession,logger)
#####################################################################################################
        
####RUN DMR calling ################################################################################
if args.sampleInfo :
    DMRout=os.path.join(wdir,'metilene_out')
    @mkdir(DMRout,os.path.join(DMRout,'logs'))
    @transform(CpG_stats,suffix('.RData'),'.metilene.bed',output_dir=DMRout)
    def run_metilene (input_files,output_file):
        ii = input_files[2]
        oo = output_file
        from BS_DMR_WGBS import DMR_metilene
        DMR_metilene(ii,args.sampleInfo,oo,args.nthreads,metipath,mySession,logger)
    @follows(run_metilene)
    @transform(refG,suffix('.fa'),'.metilene.CpGlist.bed',output_dir=auxdir)
    def get_CG_metilene(input_file,output_file):
        ii=input_file
        oo=output_file
        from BSstats_WGBS import mCT_get_CpGxInt
        mCT_get_CpGxInt(ii,[oo],refG,[os.path.join(DMRout,'singleCpG.metilene.bed')],auxdir,bedpath,mySession,logger)

    os.chdir(CpGstat_out)
    @follows(get_CG_metilene)
    @merge(input=[run_metilene,CpG_stats],output=[os.path.join(DMRout,"singleCpG.metilene.limma.bed"),os.path.join(DMRout,"metilene.limma.annotated.txt")])
    def cleanup_metilene (input_files,output_files):
        ii1 = input_files[0]
        ii2 = input_files[1][1]
        oo = output_files
        auxList=os.path.join(auxdir,re.sub('.fa','.metilene.CpGlist.bed',os.path.basename(refG)))
        from BS_DMR_WGBS import clean_up_metilene
        clean_up_metilene(ii1,ii2,auxList,refG,args.sampleInfo,DMRout,Rpath,Rlib,bedpath,pipev,mySession,logger)
###################################################################################################

#####main

if __name__ == '__main__':
    with open(os.path.join(wdir,"pipelineGraph.png"),'w') as pipeGraph:
        pipeline_printout_graph(stream=pipeGraph,output_format='png',pipeline_name='WGBS',target_tasks=args.target_tasks)
    with open (os.path.join(wdir,"pipelinePrint.txt"),'w') as pipePrint:
        pipeline_printout(verbose_abbreviated_path=0,output_stream=pipePrint,target_tasks=args.target_tasks)    

    pipeline_run(touch_files_only=args.touchOnly,multiprocess=args.bsize,target_tasks=args.target_tasks,forcedtorun_tasks=args.forcedtorun_tasks,logger=logger)
    mySession.exit()