##################################
#                                #
# Last modified 04/20/2011       # 
#                                #
# Georgi Marinov                 #
#                                # 
##################################

import sys
import string

def run():

    if len(sys.argv) < 2:
        print 'usage: python %s inputfile outputfilename [-BAM]' % sys.argv[0]
        print 'inputfile format: <output prefix> <tab> <bowtie file name> <tab> <control bowtie file name>'
        sys.exit(1)
    
    input = sys.argv[1]
    outfilename = sys.argv[2]

    outfile = open(outfilename, 'w')

    doBAM=False
    if '-BAM' in sys.argv:
        doBAM=True

    inputdatafile = open(input)
    for line in inputdatafile:
        fields=line.strip().split('\t')
        prefix=fields[0]
        bowtie=fields[1]
        control=fields[2]
        outfile.write('library(spp);'+'\n')
        if doBAM:
            outfile.write('chip.data <- read.bam.tags("' + bowtie + '");'+'\n')
            outfile.write('input.data <- read.bam.tags("' + control + '");'+'\n')
        else:
            outfile.write('chip.data <- read.bowtie.tags("' + bowtie + '");'+'\n')
            outfile.write('input.data <- read.bowtie.tags("' + control + '");'+'\n')
        outfile.write('binding.characteristics <- get.binding.characteristics(chip.data,srange=c(0,400),bin=2);'+'\n')
        outfile.write('print(paste("binding peak separation distance =",binding.characteristics$peak$x));'+'\n')
        outfile.write('pdf(file="' + prefix + '.crosscorrelation.pdf",width=5,height=5);'+'\n')
        outfile.write('par(mar = c(3.5,3.5,1.0,0.5), mgp = c(2,0.65,0), cex = 0.8);'+'\n')
        outfile.write("plot(binding.characteristics$cross.correlation,type='l',xlab=" + '"strand shift"' + ",ylab=" + '"cross-correlation");\n')
        outfile.write('abline(v=binding.characteristics$peak$x,lty=2,col=2);'+'\n')
        outfile.write('dev.off();'+'\n')
        outfile.write('chip.data <- select.informative.tags(chip.data,binding.characteristics);'+'\n')
        outfile.write('input.data <- select.informative.tags(input.data,binding.characteristics);'+'\n')
        outfile.write('chip.data <- remove.local.tag.anomalies(chip.data);'+'\n')
        outfile.write('input.data <- remove.local.tag.anomalies(input.data);'+'\n')
        outfile.write('tag.shift <- round(binding.characteristics$peak$x/2)'+'\n')
        outfile.write('smoothed.density <- get.smoothed.tag.density(chip.data,control.tags=input.data,bandwidth=200,step=100,tag.shift=tag.shift);'+'\n')
        outfile.write('writewig(smoothed.density,"' + prefix + '.density.wig","Example smoothed, background-subtracted tag density");'+'\n')
        outfile.write('rm(smoothed.density);'+'\n')
        outfile.write('enrichment.estimates <- get.conservative.fold.enrichment.profile(chip.data,input.data,fws=500,step=100,alpha=0.01);'+'\n')
        outfile.write('writewig(enrichment.estimates,"' + prefix + '.enrichment.estimates.wig","Example conservative fold-enrichment/depletion estimates shown on log2 scale");'+'\n')
        outfile.write('rm(enrichment.estimates);'+'\n')
        outfile.write('broad.clusters <- get.broad.enrichment.clusters(chip.data,input.data,window.size=1e3,z.thr=3,tag.shift=round(binding.characteristics$peak$x/2))'+'\n')
        outfile.write('write.broadpeak.info(broad.clusters,"' + prefix + '.broadPeak")'+'\n')
        outfile.write('fdr <- 1e-2; '+'\n')
        outfile.write('detection.window.halfsize <- binding.characteristics$whs;'+'\n')
        outfile.write('bp <- find.binding.positions(signal.data=chip.data,control.data=input.data,fdr=fdr,whs=detection.window.halfsize)'+'\n')
        outfile.write('print(paste("detected",sum(unlist(lapply(bp$npl,function(d) length(d$x)))),"peaks"));'+'\n')
        outfile.write('output.binding.results(bp,"' + prefix + '.FDR-1e-2.binding.positions.txt");'+'\n')
        outfile.write('bp <- add.broad.peak.regions(chip.data,input.data,bp,window.size=100,z.thr=3)'+'\n')
        outfile.write('write.narrowpeak.binding(bp,"' + prefix + '-FRD1e-2.narrowPeak")'+'\n')
        outfile.write('fdr <- 1e-3; '+'\n')
        outfile.write('detection.window.halfsize <- binding.characteristics$whs;'+'\n')
        outfile.write('bp <- find.binding.positions(signal.data=chip.data,control.data=input.data,fdr=fdr,whs=detection.window.halfsize)'+'\n')
        outfile.write('print(paste("detected",sum(unlist(lapply(bp$npl,function(d) length(d$x)))),"peaks"));'+'\n')
        outfile.write('output.binding.results(bp,"' + prefix + '.FDR-1e-3.binding.positions.txt");'+'\n')
        outfile.write('bp <- add.broad.peak.regions(chip.data,input.data,bp,window.size=100,z.thr=3)'+'\n')
        outfile.write('write.narrowpeak.binding(bp,"' + prefix + '-FRD1e-3.narrowPeak")'+'\n')
        outfile.write('\n')

    outfile.close()
   
run()