##################################
#                                #
# Last modified 03/11/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>'
        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]
        outfile.write('library(spp);'+'\n')
        if doBAM:
            outfile.write('chip.data <- read.bam.tags("' + bowtie + '");'+'\n')
        else:
            outfile.write('chip.data <- read.bowtie.tags("' + bowtie + '");'+'\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=")
        outfile.write('"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('\n')

    outfile.close()
   
run()