"""
@summary: Get lists of eQTLs liked to significant hotspots

@author: nanette.coetzer@gmail.com
@version 0.1

"""
import optparse, sys
import subprocess
import tempfile
import os, re

def stop_err( msg ):
    sys.stderr.write( "%s\n" % msg )
    sys.exit()
 
def __main__():
    #Parse Command Line
    parser = optparse.OptionParser()
    parser.add_option("-a", "--rscript1", default=None, dest="rscript1", 
                      help="R script 1")
    parser.add_option("-b", "--rscript2", default=None, dest="rscript2", 
                      help="R script 2")
    parser.add_option("-i", "--input1", default=None, dest="input1", 
                      help="Gene and eQTL frequency per bin table")
    parser.add_option("-j", "--input2", default=None, dest="input2", 
                      help="Gene and eQTL sliding frequency")
    parser.add_option("-k", "--input3", default=None, dest="input3", 
                      help="Frequency summary file")
    parser.add_option("-l", "--input4", default=None, dest="input4", 
                      help="all eQTL classification file")
    parser.add_option("-m", "--input5", default=None, dest="input5", 
                      help="cis-eQTL classification file")
    parser.add_option("-n", "--input6", default=None, dest="input6", 
                      help="trans-eQTL classification file")
    parser.add_option("-g", "--input7", default=None, dest="input7", 
                      help="Sliding map file")
    
    parser.add_option("-o", "--output1", default=None, dest="output1", 
                      help="Unbiased all eQTL hotspots result file")
    parser.add_option("-p", "--output2", default=None, dest="output2", 
                      help="List of all eQTLs per hotspot")
    parser.add_option("-q", "--output3", default=None, dest="output3", 
                      help="Unbiased cis eQTL hotspots result file")
    parser.add_option("-r", "--output4", default=None, dest="output4", 
                      help="List of cis eQTLs per hotspot")
    parser.add_option("-s", "--output5", default=None, dest="output5", 
                      help="Unbiased trans eQTL hotspots result file")
    parser.add_option("-t", "--output6", default=None, dest="output6", 
                      help="List of trans eQTLs per hotspot")
    parser.add_option("-u", "--output7", default=None, dest="output7", 
                      help="Full summary file")
    parser.add_option("-v", "--output8", default=None, dest="output8", 
                      help="zip folder with R plots")
    (options, args) = parser.parse_args()

    try:
        open(options.input1, "r").close()
    except TypeError, e:
        stop_err("You need to supply the Gene and eQTL frequency file:\n" + str(e))
    except IOError, e:
        stop_err("Can not open the Gene and eQTL frequency file:\n" + str(e))
        
    try:
        open(options.input2, "r").close()
    except TypeError, e:
        stop_err("You need to supply the Gene and eQTL sliding window frequency file:\n" + str(e))
    except IOError, e:
        stop_err("Can not open the Gene and eQTL sliding window frequency file:\n" + str(e))

    try:
        open(options.input3, "r").close()
    except TypeError, e:
        stop_err("You need to supply the Frequency summary file:\n" + str(e))
    except IOError, e:
        stop_err("Can not open the Frequency summary file:\n" + str(e))

    try:
        open(options.input4, "r").close()
    except TypeError, e:
        stop_err("You need to supply the all eQTL classification file:\n" + str(e))
    except IOError, e:
        stop_err("Can not open the all eQTL classification file:\n" + str(e))

    try:
        open(options.input5, "r").close()
    except TypeError, e:
        stop_err("You need to supply the cis-eQTL classification file:\n" + str(e))
    except IOError, e:
        stop_err("Can not open the cis-eQTL classification file:\n" + str(e))

    try:
        open(options.input6, "r").close()
    except TypeError, e:
        stop_err("You need to supply the trans-eQTL classification file:\n" + str(e))
    except IOError, e:
        stop_err("Can not open the trans-eQTL classification file:\n" + str(e))
  
    try:
        open(options.input7, "r").close()
    except TypeError, e:
        stop_err("You need to supply the Sliding-Lookup file:\n" + str(e))
    except IOError, e:
        stop_err("Can not open the Sliding-Lookup file:\n" + str(e))

    

    ##########################################################
    
    # Create temp direcotry
    tempdir = tempfile.mkdtemp()

    # copy INPUT file to the temp directory
    s = "cp %s %s/frequency_summary.txt" %(options.input3, tempdir)
    subprocess.call(s, shell=True)
    s = "cp %s %s/sliding_frequency.txt" %(options.input2, tempdir)
    subprocess.call(s, shell=True)

    # create R script => save in temp directory
    # generate new header 
    new_script = open(tempdir+"/new_script.r","w")
    header = "setwd(\"%s\")" %tempdir
    new_script.write(header+"\n")
    # add script body
    script = open(options.rscript1,"r")
    for line in script:
        new_script.write(line.strip()+"\n")
    new_script.close()
    
    # run R script from temp directory
    s = "R CMD BATCH %s/new_script.r out.txt" %tempdir
    subprocess.call(s, shell=True)
    
    ###
    
    # get R output file and write to options.output --> all
    outfile = open(tempdir+"/chisq_out_all2.txt",'w')     # chisq output file
    r_outfile = open(tempdir+"/chisq_out_all.txt","r")
    i = 0
    for r_line in r_outfile:
        i += 1
        if i == 1:
            outfile.write("bin\tsliding.eQTL\tsliding.genes\texp.value.ok\tsliding.cM\tchr\tchisq.result\n")
        else:
            rl = r_line.strip().split(" ")
            new_rl = "\t".join(rl[1:])
            outfile.write(new_rl+"\n")
    outfile.close()
    r_outfile.close()

    # get R output file and write to options.output --> cis
    outfile = open(tempdir+"/chisq_out_cis2.txt",'w')     # chisq output file
    r_outfile = open(tempdir+"/chisq_out_cis.txt","r")
    i = 0
    for r_line in r_outfile:
        i += 1
        if i == 1:
            outfile.write("bin\tsliding.eQTL\tsliding.genes\texp.value.ok\tsliding.cM\tchr\tchisq.result\n")
        else:
            rl = r_line.strip().split(" ")
            new_rl = "\t".join(rl[1:])
            outfile.write(new_rl+"\n")
    outfile.close()
    r_outfile.close()
    
    # get R output file and write to options.output --> trans
    outfile = open(tempdir+"/chisq_out_trans2.txt",'w')     # chisq output file
    r_outfile = open(tempdir+"/chisq_out_trans.txt","r")
    i = 0
    for r_line in r_outfile:
        i += 1
        if i == 1:
            outfile.write("bin\tsliding.eQTL\tsliding.genes\texp.value.ok\tsliding.cM\tchr\tchisq.result\n")
        else:
            rl = r_line.strip().split(" ")
            new_rl = "\t".join(rl[1:])
            outfile.write(new_rl+"\n")
    outfile.close()
    r_outfile.close()

    #####
    
    # get R output file and write to options.output
    r_outfile2 = open(tempdir+"/threshold_chisq_summary.txt","r")   
    for r_line2 in r_outfile2:
        if r_line2.startswith("\"1\""):
            threshold_all = r_line2.strip().split(" ")[1]
        if r_line2.startswith("\"2\""):
            t_all = r_line2.strip().split(" ")[1]
        if r_line2.startswith("\"3\""):
            g_left_all = r_line2.strip().split(" ")[1]           
        if r_line2.startswith("\"4\""):
            e_right_all = r_line2.strip().split(" ")[1]
	
	if r_line2.startswith("\"5\""):
            threshold_cis = r_line2.strip().split(" ")[1]
        if r_line2.startswith("\"6\""):
            t_cis = r_line2.strip().split(" ")[1]
        if r_line2.startswith("\"7\""):
            g_left_cis = r_line2.strip().split(" ")[1]           
        if r_line2.startswith("\"8\""):
            e_right_cis = r_line2.strip().split(" ")[1]
	
	if r_line2.startswith("\"9\""):
            threshold_trans = r_line2.strip().split(" ")[1]
        if r_line2.startswith("\"10\""):
            t_trans = r_line2.strip().split(" ")[1]
        if r_line2.startswith("\"11\""):
            g_left_trans = r_line2.strip().split(" ")[1]           
        if r_line2.startswith("\"12\""):
            e_right_trans = r_line2.strip().split(" ")[1]

    outfile7 = open(options.output7,'w')
    infile3 = open(options.input3,'r')
    for line in infile3:
        l = line.strip()
        outfile7.write(l+"\n")
	#if line.startswith("eQTL"):
	    #eqtl_option = line.strip().split("\t")[1]
    outfile7.write("Calculated permutation threshold for all eQTL (eQTL/cM)\t"+str(round(float(threshold_all),2))+"\n")
    perm_threshold_all = round(float(threshold_all),2)
    outfile7.write("Chi-squared test population estimate for all eQTL (genes:eQTL)\t"+str(round(float(g_left_all),3))+":"+str(round(float(e_right_all),3))+"\n")
    outfile7.write("Chi-squared test assumption for all eQTL: min number of [eQTL + genes] in bin\t"+str(round(float(t_all),2))+"\n")
    
    outfile7.write("Calculated permutation threshold for all eQTL (eQTL/cM)\t"+str(round(float(threshold_cis),2))+"\n")
    perm_threshold_cis = round(float(threshold_cis),2)
    outfile7.write("Chi-squared test population estimate for cis-eQTL (genes:eQTL)\t"+str(round(float(g_left_cis),3))+":"+str(round(float(e_right_cis),3))+"\n")
    outfile7.write("Chi-squared test assumption for cis-eQTL: min number of [eQTL + genes] in bin\t"+str(round(float(t_cis),2))+"\n")
    
    outfile7.write("Calculated permutation threshold for all eQTL (eQTL/cM)\t"+str(round(float(threshold_trans),2))+"\n")
    perm_threshold_trans = round(float(threshold_trans),2)
    outfile7.write("Chi-squared test population estimate for trans-eQTL (genes:eQTL)\t"+str(round(float(g_left_trans),3))+":"+str(round(float(e_right_trans),3))+"\n")
    outfile7.write("Chi-squared test assumption for trans-eQTL: min number of [eQTL + genes] in bin\t"+str(round(float(t_trans),2))+"\n")
    
    outfile7.close()
    infile3.close()
    
    ##########################################################

    ############################################
    # All eQTL hotspots and lists
    hotspot_result_all = open(options.output1,'w')    # unbiased_hotspots_result.txt
    gene_lists_results_all = open(options.output2,'w')    # eQTL_list_per_hotspot.txt
    hotspot_result_cis = open(options.output3,'w')    # unbiased_hotspots_result.txt
    gene_lists_results_cis = open(options.output4,'w')    # eQTL_list_per_hotspot.txt
    hotspot_result_trans = open(options.output5,'w')    # unbiased_hotspots_result.txt
    gene_lists_results_trans = open(options.output6,'w')    # eQTL_list_per_hotspot.txt
    
    ############################################
    # freq_dict
    freq = open(options.input1,'r')  # lookup_freq.txt
    freq_dict_all = {}
    freq_dict_cis = {}
    freq_dict_trans = {}
    lookup_tmp_dict = {}
    header = 0
    for line in freq:
	l = line.strip().split("\t")
	if l[1].startswith("c"):
	    header = 1
	else:
	    header = 0
	if header == 0:
		freq_dict_all[l[0]] = int(float(l[7]))
		freq_dict_cis[l[0]] = int(float(l[8]))
		freq_dict_trans[l[0]] = int(float(l[9]))
		lookup_tmp_dict[l[0]] = [l[1],str(round(float(l[3]),4))]
    freq.close()
    
    ############################################
    # s_freq_dict
    # sliding.id: [chr, sliding.cM, sliding.eQTL, sliding.genes, eQTL/cM, sliding.eQTL/cM/10genes, threshold]
    s_freq = open(options.input2,'r')  
    s_freq_dict_all = {}
    s_freq_dict_cis = {}
    s_freq_dict_trans = {}
    j = 0
    for line in s_freq:
	j += 1
	if j > 1:
	    l = line.strip().split("\t")
	    s_freq_dict_all[l[0]] = [l[0],l[1],l[2],l[3],l[6],l[7]]
	    s_freq_dict_cis[l[0]] = [l[0],l[1],l[2],l[4],l[6],l[8]]
	    s_freq_dict_trans[l[0]] = [l[0],l[1],l[2],l[5],l[6],l[9]]
	    # add threshold columns
	    if float(l[7]) > float(perm_threshold_all):
		s_freq_dict_all[l[0]].append('1')
	    else:
		s_freq_dict_all[l[0]].append('0')
	    if float(l[8]) > float(perm_threshold_cis):
		s_freq_dict_cis[l[0]].append('1')
	    else:
		s_freq_dict_cis[l[0]].append('0')
	    if float(l[9]) > float(perm_threshold_trans):
		s_freq_dict_trans[l[0]].append('1')
	    else:
		s_freq_dict_trans[l[0]].append('0')
    s_freq.close()
    
    # sliding.id: [chr, sliding.cM, sliding.eQTL, sliding.genes, sliding.eQTL/cM/10genes, eQTL/cM, threshold, chisq.result, hotspot.nr]
    # WRITE hotspot_result --> all eQTLs
    chisq_results = open(tempdir+"/chisq_out_all2.txt",'r')
    i = 0
    hotspot_nr = 0
    prev_id = 0
    prev_chr = 0
    hotspot_slID_dict_all = {}
    #hotspot_result.write("sliding.id\tchr\tsliding.cM\tsliding.eQTL\tsliding.genes\tsliding.eQTL/cM\tsliding.eQTL/cM/10genes\tthreshold.result\tchisq.result\tsignificant\thotspot.number\n")
    hotspot_result_all.write("sliding.id\tchr\tsliding.cM\tsliding.all.eQTL\tsliding.genes\tsliding.all.eQTL/cM\tthreshold.result\tchisq.result\tsignificant\thotspot.number\n")
    for line in chisq_results:
        i += 1
        if i > 1:
            l = line.strip().split("\t")
            s_freq_dict_all[l[0]].append(l[6])
	    if (s_freq_dict_all[l[0]][6] == '1') and (s_freq_dict_all[l[0]][7] == '1'):
		id = s_freq_dict_all[l[0]][0]
		s_freq_dict_all[l[0]].append('*')
		if (prev_id == int(id)-1) and (s_freq_dict_all[l[0]][1] == prev_chr):
		    s_freq_dict_all[l[0]].append(str(hotspot_nr))
		    hotspot_slID_dict_all[hotspot_nr].append(int(id))
		else:
		    hotspot_nr += 1
		    s_freq_dict_all[l[0]].append(str(hotspot_nr))
		    hotspot_slID_dict_all[hotspot_nr] = [int(id)]
		prev_id = int(id)
		prev_chr = s_freq_dict_all[l[0]][1]
	    else:
		s_freq_dict_all[l[0]].append("\t")
	    joinline = "\t".join(s_freq_dict_all[l[0]])
	    hotspot_result_all.write(joinline+"\n")
    chisq_results.close()
    hotspot_result_all.close()
    
    # WRITE hotspot_result --> cis eQTLs
    chisq_results = open(tempdir+"/chisq_out_cis2.txt",'r')
    i = 0
    hotspot_nr = 0
    prev_id = 0
    prev_chr = 0
    hotspot_slID_dict_cis = {}
    #hotspot_result.write("sliding.id\tchr\tsliding.cM\tsliding.eQTL\tsliding.genes\tsliding.eQTL/cM\tsliding.eQTL/cM/10genes\tthreshold.result\tchisq.result\tsignificant\thotspot.number\n")
    hotspot_result_cis.write("sliding.id\tchr\tsliding.cM\tsliding.all.eQTL\tsliding.genes\tsliding.all.eQTL/cM\tthreshold.result\tchisq.result\tsignificant\thotspot.number\n")
    for line in chisq_results:
        i += 1
        if i > 1:
            l = line.strip().split("\t")
            s_freq_dict_cis[l[0]].append(l[6])
	    if (s_freq_dict_cis[l[0]][6] == '1') and (s_freq_dict_cis[l[0]][7] == '1'):
		id = s_freq_dict_cis[l[0]][0]
		s_freq_dict_cis[l[0]].append('*')
		if (prev_id == int(id)-1) and (s_freq_dict_cis[l[0]][1] == prev_chr):
		    s_freq_dict_cis[l[0]].append(str(hotspot_nr))
		    hotspot_slID_dict_cis[hotspot_nr].append(int(id))
		else:
		    hotspot_nr += 1
		    s_freq_dict_cis[l[0]].append(str(hotspot_nr))
		    hotspot_slID_dict_cis[hotspot_nr] = [int(id)]
		prev_id = int(id)
		prev_chr = s_freq_dict_cis[l[0]][1]
	    else:
		s_freq_dict_cis[l[0]].append("\t")
	    joinline = "\t".join(s_freq_dict_cis[l[0]])
	    hotspot_result_cis.write(joinline+"\n")
    chisq_results.close()
    hotspot_result_cis.close()
    
    # WRITE hotspot_result --> trans eQTLs
    chisq_results = open(tempdir+"/chisq_out_trans2.txt",'r')
    i = 0
    hotspot_nr = 0
    prev_id = 0
    prev_chr = 0
    hotspot_slID_dict_trans = {}
    #hotspot_result.write("sliding.id\tchr\tsliding.cM\tsliding.eQTL\tsliding.genes\tsliding.eQTL/cM\tsliding.eQTL/cM/10genes\tthreshold.result\tchisq.result\tsignificant\thotspot.number\n")
    hotspot_result_trans.write("sliding.id\tchr\tsliding.cM\tsliding.all.eQTL\tsliding.genes\tsliding.all.eQTL/cM\tthreshold.result\tchisq.result\tsignificant\thotspot.number\n")
    for line in chisq_results:
        i += 1
        if i > 1:
            l = line.strip().split("\t")
            s_freq_dict_trans[l[0]].append(l[6])
	    if (s_freq_dict_trans[l[0]][6] == '1') and (s_freq_dict_trans[l[0]][7] == '1'):
		id = s_freq_dict_trans[l[0]][0]
		s_freq_dict_trans[l[0]].append('*')
		if (prev_id == int(id)-1) and (s_freq_dict_trans[l[0]][1] == prev_chr):
		    s_freq_dict_trans[l[0]].append(str(hotspot_nr))
		    hotspot_slID_dict_trans[hotspot_nr].append(int(id))
		else:
		    hotspot_nr += 1
		    s_freq_dict_trans[l[0]].append(str(hotspot_nr))
		    hotspot_slID_dict_trans[hotspot_nr] = [int(id)]
		prev_id = int(id)
		prev_chr = s_freq_dict_trans[l[0]][1]
	    else:
		s_freq_dict_trans[l[0]].append("\t")
	    joinline = "\t".join(s_freq_dict_trans[l[0]])
	    hotspot_result_trans.write(joinline+"\n")
    chisq_results.close()
    hotspot_result_trans.close()

    ############################################
    #tot_in_hotspots_all = 0
    #tot_in_hotspots_cis = 0
    #tot_in_hotspots_trans = 0
    sl_lookup_dict = {}
    slidinglookup = open(options.input7,'r')  # threshold_output.txt
    for sline in slidinglookup:
    	sl = sline.strip().split("\t")
	sl_lookup_dict[sl[0]] = sl[1][1:-1].split(",")
    slidinglookup.close()
    
    # for each hotspot, for each sliding.id, get the int.id
    # WRITE gene_lists_results --> all
    for i in range(0,len(hotspot_slID_dict_all.keys())):
    	hid = hotspot_slID_dict_all.keys()[i]
	sl_ids = hotspot_slID_dict_all[int(i)+1]
	lu_ids = []
	for j in sl_ids:
	    for k in sl_lookup_dict[str(j)]:
	        if int(k.strip()) not in lu_ids:
		    lu_ids.append(int(k.strip()))
	chr = lookup_tmp_dict[str(lu_ids[0])][0]

	freq_count = 0
	for luid in lu_ids:
	    freq_count = freq_count + freq_dict_all[str(luid)]
	#tot_in_hotspots_all = tot_in_hotspots_all + freq_count
	gene_lists_results_all.write("\n\n= = = = =   Hotspot "+str(hid)+"  chr "+str(chr)+"	sliding.ids: "+str(sl_ids)+"	int.ids: "+str(lu_ids)+"	nr.eQTL: "+str(freq_count)+"   = = = = =\n")
	for l in lu_ids:
	    intpos = lookup_tmp_dict[str(l)]
	    eqtl = open(options.input4,'r')  # eqtls.txt 
	    for eqtlline in eqtl:
	    	el = eqtlline.strip().split("\t")
	    	if not el[2].startswith("c"):
	    	    pos = [el[2],str(round(float(el[8]),4))]
	    	    if intpos == pos:
			gene_lists_results_all.write("\t".join(el)+"\n")
	    eqtl.close()
    gene_lists_results_all.close()
    
    # WRITE gene_lists_results --> cis
    for i in range(0,len(hotspot_slID_dict_cis.keys())):
    	hid = hotspot_slID_dict_cis.keys()[i]
	sl_ids = hotspot_slID_dict_cis[int(i)+1]
	lu_ids = []
	for j in sl_ids:
	    for k in sl_lookup_dict[str(j)]:
	        if int(k.strip()) not in lu_ids:
		    lu_ids.append(int(k.strip()))
	chr = lookup_tmp_dict[str(lu_ids[0])][0]
	
	freq_count = 0
	for luid in lu_ids:
	    freq_count = freq_count + freq_dict_cis[str(luid)]
	#tot_in_hotspots = tot_in_hotspots + freq_count
	gene_lists_results_cis.write("\n\n= = = = =   Hotspot "+str(hid)+"  chr "+str(chr)+"	sliding.ids: "+str(sl_ids)+"	int.ids: "+str(lu_ids)+"	nr.eQTL: "+str(freq_count)+"   = = = = =\n")
	for l in lu_ids:
	    intpos = lookup_tmp_dict[str(l)]
	    eqtl = open(options.input5,'r')  # eqtls.txt 
	    for eqtlline in eqtl:
	    	el = eqtlline.strip().split("\t")
	    	if not el[2].startswith("c"):
	    	    pos = [el[2],str(round(float(el[8]),4))]
	    	    if intpos == pos:
			gene_lists_results_cis.write("\t".join(el)+"\n")
	    eqtl.close()
    gene_lists_results_cis.close()
    
    # WRITE gene_lists_results --> trans
    for i in range(0,len(hotspot_slID_dict_trans.keys())):
    	hid = hotspot_slID_dict_trans.keys()[i]
	sl_ids = hotspot_slID_dict_trans[int(i)+1]
	lu_ids = []
	for j in sl_ids:
	    for k in sl_lookup_dict[str(j)]:
	        if int(k.strip()) not in lu_ids:
		    lu_ids.append(int(k.strip()))
	chr = lookup_tmp_dict[str(lu_ids[0])][0]
	freq_count = 0
	for luid in lu_ids:
	    freq_count = freq_count + freq_dict_trans[str(luid)]
	#tot_in_hotspots = tot_in_hotspots + freq_count
	gene_lists_results_trans.write("\n\n= = = = =   Hotspot "+str(hid)+"  chr "+str(chr)+"	sliding.ids: "+str(sl_ids)+"	int.ids: "+str(lu_ids)+"	nr.eQTL: "+str(freq_count)+"   = = = = =\n")
	for l in lu_ids:
	    intpos = lookup_tmp_dict[str(l)]
	    eqtl = open(options.input6,'r')  # eqtls.txt 
	    for eqtlline in eqtl:
	    	el = eqtlline.strip().split("\t")
	    	if not el[2].startswith("c"):
	    	    pos = [el[2],str(round(float(el[8]),4))]
	    	    if intpos == pos:
			gene_lists_results_trans.write("\t".join(el)+"\n")
	    eqtl.close()
    gene_lists_results_trans.close()
    
    
    ########################################## R plots ##########################################
    
    # Create temp direcotry
    tempdir = tempfile.mkdtemp()

    # copy INPUT file to the temp directory
    s = "cp %s %s/all_hotspots_out.txt" %(options.output1, tempdir)
    subprocess.call(s, shell=True)
    s = "cp %s %s/cis_hotspots_out.txt" %(options.output3, tempdir)
    subprocess.call(s, shell=True)
    s = "cp %s %s/trans_hotspots_out.txt" %(options.output5, tempdir)
    subprocess.call(s, shell=True)
    s = "cp %s %s/full_summary.txt" %(options.output7, tempdir)
    subprocess.call(s, shell=True)

    # create R script => save in temp directory
    # generate new header
    new_script = open(tempdir+"/new_script.r","w")
    header = "setwd(\"%s\")" %tempdir
    new_script.write(header+"\n")
    # add script body
    
    script = open(options.rscript2,"r")
    for line in script:
        new_script.write(line.strip()+"\n")
    new_script.close()
    
    # run R script from temp directory
    s = "R CMD BATCH %s/new_script.r out.txt" %tempdir
    subprocess.call(s, shell=True)

    tempdir2 = tempfile.mkdtemp()
    #os.system("mv %s/genome_wide_EcM10G.pdf %s/genome_wide_EcM10G.pdf" %(tempdir,tempdir2))
    #os.system("mv %s/chr_plots_EcM10G.pdf %s/chr_plots_EcM10G.pdf" %(tempdir,tempdir2))
    os.system("mv %s/genome_wide_EcM_all.pdf %s/genome_wide_EcM_all.pdf" %(tempdir,tempdir2))
    os.system("mv %s/chr_plots_EcM_all.pdf %s/chr_plots_EcM_all.pdf" %(tempdir,tempdir2))
    os.system("mv %s/genome_wide_EcM_cis.pdf %s/genome_wide_EcM_cis.pdf" %(tempdir,tempdir2))
    os.system("mv %s/chr_plots_EcM_cis.pdf %s/chr_plots_EcM_cis.pdf" %(tempdir,tempdir2))
    os.system("mv %s/genome_wide_EcM_trans.pdf %s/genome_wide_EcM_trans.pdf" %(tempdir,tempdir2))
    os.system("mv %s/chr_plots_EcM_trans.pdf %s/chr_plots_EcM_trans.pdf" %(tempdir,tempdir2))
    s = "zip -j %s/output.zip %s/*" %(tempdir,tempdir2)
    subprocess.call(s, shell=True)
    
    # move the zipped temp directory to options.output
    os.system("mv %s/output.zip %s" %(tempdir,options.output8))
    
    

    ##############################################
    
if __name__=="__main__": 
    __main__()



