def importMACS( infile, outfile, suffix = ".norm.bam" ):
    '''import MACS results.

    Imports only positive peaks. It filters peaks by p-value,
    q-value and fold change and imports the diagnostic data.

    Does re-counting of peakcenter, peakval, ... using
    the normalized tag counts.
    '''

    track = infile[:-len(".macs")]    
    infilename = infile + "_peaks.xls"
    filename_diag = infile + "_diag.xls"
    filename_r = infile + "_model.r"
    
    if not os.path.exists(infilename):
        E.warn("could not find %s" % infilename )
        outs = open(outfile,"w")
        outs.close()
        return

    shift = getPeakShift( infile )
    assert shift != None, "could not determine peak shift from MACS file %s" % infile

    samfiles = [ pysam.Samfile( track + suffix, "rb" ) ]
    offsets = [ shift / 2 ]
    outtemp = P.getTempFile()

    outtemp.write( "\t".join( ( \
                "interval_id", 
                "contig", "start", "end",
                "npeaks", "peakcenter", 
                "length", 
                "avgval", "peakval",
                "nprobes",
                "pvalue", "fold", "qvalue",
                "macs_summit", "macs_nprobes",
                )) + "\n" )
    id = 0

    ## get thresholds
    max_qvalue = float(PARAMS["macs_max_qvalue"])
    # min, as it is -10log10
    min_pvalue = float(PARAMS["macs_min_pvalue"])
    min_fold = float(PARAMS["macs_min_fold"])
    
    counter = E.Counter()
    with open( infilename, "r" ) as ins:
        for line in ins:
            if line.startswith("#"): continue
            if line.startswith( "chr\tstart"): continue
            counter.input += 1
            data = line[:-1].split("\t")
            if len(data) == 9:
                contig,start,end,length,summit,ntags,pvalue,fold,qvalue = data
            elif len(data) == 8:
                contig,start,end,length,summit,ntags,pvalue,fold = data
                qvalue = 1.0
            else:
                raise ValueError( "could not parse line %s" % line )
            
            pvalue, qvalue, summit, fold = float(pvalue), float(qvalue), int(summit), float(fold)

            if qvalue > max_qvalue or pvalue < min_pvalue or fold < min_fold: 
                counter.skipped += 1
                continue

            # these are 1-based coordinates
            start, end = int(start)-1, int(end)
            assert start < end
            # macs can have negative start coordinates
            start = max(start, 0)
            npeaks, peakcenter, length, avgval, peakval, nreads = countPeaks( contig, start, end, samfiles, offsets )

            outtemp.write ( "\t".join( map(str, ( \
                            id, contig, start, end, npeaks, peakcenter, length, avgval, peakval, nreads,
                            pvalue, fold, qvalue,
                            start + summit - 1, 
                            ntags) ) ) + "\n" )
            id += 1                        
            counter.output += 1
    outtemp.close()

    tablename = "%s_intervals" % track
    tmpfilename = outtemp.name

    statement = '''
   python %(scriptsdir)s/csv2db.py %(csv2db_options)s \
              --index=interval_id \
              --table=%(tablename)s \
    < %(tmpfilename)s > %(outfile)s
    '''

    P.run( **dict( locals().items() + PARAMS.items() ) )

    # import diagnotic data
    if os.path.exists( filename_diag ):

        tablename = "%s_macsdiag" % track

        statement = '''
        sed "s/FC range.*/fc\\tnpeaks\\tp90\\tp80\\tp70\\tp60\\tp50\\tp40\\tp30\\tp20/" < %(filename_diag)s |\
       python %(scriptsdir)s/csv2db.py %(csv2db_options)s \
                  --map=fc:str \
                  --table=%(tablename)s \
        > %(outfile)s
        '''

        P.run( **dict( locals().items() + PARAMS.items() ) )


    # create plot
    if os.path.exists( filename_r ):

        target_path = os.path.join( os.getcwd(), "doc", "_static", "MACS" )
        try:
            os.makedirs( target_path )
        except __HOLE__: 
            # ignore "file exists" exception
            pass

        statement = '''
        R --vanilla < %(track)s.macs_model.r > %(outfile)s
        '''
        
        P.run( **dict( locals().items() + PARAMS.items() ) )

        shutil.copyfile(
            "%s.macs_model.pdf" % track,
            os.path.join( target_path, "%s_model.pdf" % track) )
        
    os.unlink( tmpfilename )

    E.info("%s: %s" % (track, str(counter))) 

############################################################
############################################################
############################################################
dataset/ETHPy150Open CGATOxford/cgat/obsolete/pipeline_vitaminD_intervals.py/importMACS