#! /opt/exp_soft/python27/bin/python import argparse import os import re import sys import string import math gp_term= {"pdf" : "set terminal pdfcairo colour font GPFONT_SIZE",\ "postscript":"set terminal postscript eps colour font LEGEND",\ "jpeg" : "set terminal jpeg medium",\ "png" : "set terminal png",\ "aqua" : "set terminal aqua",\ "x11" : "set terminal x11 lw 1"} PS=0.75 gp_eps_styles = [(0,9,PS*1.5),(3,25,PS),(4,7,PS),(2,28,PS),(8,11,PS),(5,5,PS),(1,19,PS),(2,4,PS),(3,3,PS),(8,18,PS),(5,2,PS),(9,20,PS),(10,19,PS),(11,18,PS),(12,17,PS),(13,16,PS),(14,15,PS),(13,13,PS),(0,1,PS),(1,2,PS),(2,3,PS),(3,4,PS)] gp_x11_styles = [(0,9,PS),(3,25,PS),(4,7,PS),(2,28,PS),(8,11,PS),(5,5,PS),(1,19,PS),(2,4,PS),(3,3,PS),(8,18,PS),(5,2,PS),(9,20,PS),(10,19,PS),(11,18,PS),(12,17,PS),(13,16,PS),(14,15,PS),(13,13,PS),(0,1,PS),(1,2,PS),(2,3,PS),(3,4,PS)] gp_aqua_styles = [(1,9,PS),(3,25,PS),(4,7,PS),(2,28,PS),(8,11,PS),(5,5,PS),(1,19,PS),(2,4,PS),(3,3,PS),(8,18,PS),(5,2,PS),(9,20,PS),(10,19,PS),(11,18,PS),(12,17,PS),(13,16,PS),(14,15,PS),(13,13,PS),(0,1,PS),(1,2,PS),(2,3,PS),(3,4,PS)] gp_pdf_styles = [(0,4,PS),(1,7,PS),(2,5,PS),(3,9,PS),(4,4,PS),(5,7,PS),(7,12,PS),(8,8,PS),(1,2,PS),(2,3,PS),(4,5,PS),(5,10,PS),(6,11,PS),(7,12,PS),(8,13,PS),(9,14,PS),(10,15,PS)] gp_x11_styles = [(0,9,PS*2),(3,25,PS*2),(4,7,PS*2),(2,28,PS*2),(8,11,PS*2),(5,5,PS*2),(1,19,PS*2),(2,4,PS*2),(3,3,PS*2),(8,18,PS*2),(5,2,PS*2),(9,20,PS*2),(10,19,PS*2),(11,18,PS*2),(12,17,PS*2),(13,16,PS*2),(14,15,PS*2),(13,13,PS*2),(0,1,PS*2),(1,2,PS*2),(2,3,PS*2),(3,4,PS*2)] gp_styles= {"pdf" : gp_pdf_styles,\ "postscript": gp_eps_styles, "x11": gp_x11_styles, "aqua": gp_aqua_styles} gp_prefs = """ GPFONT="%(FONT)s" GPSIZE="%(FONTSIZE)s" LEGSIZE="%(LEGSIZE)s" GPFONT_SIZE=GPFONT.",".GPSIZE LEGEND=GPFONT.",".LEGSIZE set terminal push set key bmargin horizontal %(TERMTYPE)s """ perspective = """ set terminal push %s set view %d,%d """ def extend_style(style): extension = [] u=1 for (i,j,k) in style: extension.append((i+1,u,k*1.5)) u=u+1 for e in extension: style.append(e) for k in gp_styles.keys(): extend_style(gp_styles[k]) def do_perspectives(gout,args): # when 3D graph produced can provide different 2D perspectives persps = [(60,30),(120,120),(150,60),(330,300)] for (x,y) in persps: gout.write(perspective%(gp_term[args.gp_term],x,y)) if args.gp_output: fname = args.gp_output.replace(".","-%d-%d."%(x,y)) gout.write("set output \"%s\"\n"%fname) gout.write("\nreplot\nset terminal pop\n\n") if args.gp_output and args.epstopdf: gout.write('\nsystem "epstopdf %s";\n'%fname) if not args.gp_output: output = "pause mouse" def get_default_term(): if os.environ["OSTYPE"]=="darwin": return "aqua" else: return "x11" def setup(): parser = argparse.ArgumentParser(description='Generate Gnuplot code from evec and related files') parser.add_argument('evec', metavar='N', type=str, help='evec file') parser.add_argument('--data-has-phe', dest='dhp', action='store_true',\ default=False,\ help='evec file has phenotype (default empty)') parser.add_argument('--phe', dest='phe', action='store', default="",\ help='phenotype file (default empty)') parser.add_argument('--phe-col', dest='phe_col', action='store',\ default=2,\ help='column in phenotype file (default 2)') parser.add_argument('--outbase',dest='outbase',action='store',default="",\ help='output file name base (default evec') parser.add_argument('--omit-groups',dest='omitg',action='store',default="",\ help='groups to omit') parser.add_argument('--gp-output', dest='gp_output', action='store',\ default="", \ help="where gnuplot's output should go (default none)") parser.add_argument('--gp-term', dest='gp_term', action='store',\ default=get_default_term(), \ help='gnuplot terminal (default x11)') parser.add_argument('--epstopdf', dest='epstopdf', action='store_true',\ default=False, \ help='gp program runs epstopdf (default False)') parser.add_argument('--gp-font', dest='gp_font', action='store',\ default="Helvetica", \ help='gnuplot font (default Helvetica)') parser.add_argument('--gp-size', dest='gp_size', action='store',\ default="14", \ help='gnuplot font size (default 16)') parser.add_argument('--add-all-labels', dest='add_all_labels', \ action='store_true',default=False,\ help='add text labels (default False)') parser.add_argument('--group-label', dest='group_label', \ action='store',default="",\ help='which groups to label (default none)') parser.add_argument('--indiv-label', dest='indiv_label', \ action='store',default="",\ help='which individuals to label (default none)') parser.add_argument('--labelformat', dest='labelformat', \ action='store',default="",\ help='which part of label to use') parser.add_argument('--with-dots', dest='with_dots', \ action='store_true',default=False,\ help='draw dots (default True)') parser.add_argument('--pc', dest='pcs', action='store',\ default="1:2", \ help='principal components (default 1:2 or 1:2:3)') parser.add_argument('-3', dest='three', action='store_true', default=False, help='plot in 3d (default False)') parser.add_argument('--perspective', dest='perspective', \ action='store_true', default=False, \ help='produce perspectives for 3D images (default No)') parser.add_argument('--eigen-ratio', dest='eigen_ratio', \ action='store_true', default=False, \ help='use eigenvalue for aspect ration for 2D (default No)') args = parser.parse_args() if args.group_label: args.group_label = args.group_label.split(",") if args.dhp and args.phe: sys.exit("Mismatch can have --phe and --data-has-phe set") args.phe_col = int(args.phe_col) if not args.outbase: mm = re.search("(.*)\.pca.evec",args.evec) if not mm: sys.exit("Valid output name can't be made from <%s> -- " "expeced a file name with pca.evec suffix"%args.evec) args.outbase = "%s"%mm.group(1) args.omitg = args.omitg.split(",") if args.omitg else [] return args def get_group(args): group = {} if args.dhp: def gkey(d):return ("%s"%d[0]) args.phe_col = -1 args.phe = args.evec regexp = "([-_.\w:]+)" else: def gkey(d):return "%s:%s"%(d[0],d[1]) regexp = "([-_.\w]+)" if not args.phe: return group phef = open(args.phe) for pline in phef: data = re.findall(regexp,pline) if data[0][0] == "#": continue group[gkey(data)]= data[args.phe_col] phef.close() return group def getPCs(args): # work out the PCs to be used pcs = map(int,args.pcs.split(":")) if args.three and len(pcs)==2: pcs.append(pcs[-1]+1) the_pcs = map(lambda x: str(x+1), pcs) using_label = string.join(the_pcs,":") return (using_label,pcs) def plotIndivids(args): if not args.indiv_label: return [] if args.indiv_label[0:5] == "asis:": names = args.indiv_label[5:].split(",") else: print args.indiv_label f = open(args.indiv_label) names = map(lambda x:x.rstrip("\n"),f.readlines()) names= map(lambda x:x.replace(" ",":"),names) f.close() return names def get_eigenvectors(group,using_label,evecname): def title(x): if set(group.values()) == set(["ALL"]): return " notitle " else: print x return " title \"%s\" "%x # read in the raw data global maxes, mines plotcmd = [] if args.labelformat=="1": colfmt = ":(substr((stringcolumn(%d)),1,strstrt((stringcolumn(%d)),\":\")))" label_label = using_label+colfmt%(label_column,label_column) elif args.labelformat=="2": colfmt = ":(substr((stringcolumn(%d)),strstrt((stringcolumn(%d))+1,\":\"),strlen((stringcolumn(%d)))))" label_label = using_label+colfmt%(label_column,label_column,label_column) else: label_label = using_label+":(stringcolumn(%d))"%label_column ef = open(evecname) evecdata = {} for ev in ef: if "#" in ev: continue data = re.findall("([-_A-Za-z0-9:\.]+)",ev) if not group.has_key(data[0]): print "%s not in data\n"%(data[0]) grp = group.setdefault(data[0],"ALL") if grp in args.omitg: continue zeroes = reduce (lambda x,y: (x and float(y)==0), data[1:-2]) if zeroes: continue for (i,d) in list(enumerate(data))[1:-2]: if (float(d)>maxes[i]): maxes[i]=float(d) if (float(d) delta: return b def gnuplot_put(args,pc_indices,cmd): (plotcmd,labelcmd) = cmd # open gnuplot file and put prefix gout=open(args.outbase+".gp","w") the_term = gp_term.get(args.gp_term,"set term %s\n"%args.gp_term) GPARGS = {"FONT":args.gp_font, "FONTSIZE":args.gp_size,\ "TERMTYPE" : the_term, "LEGSIZE":str(int(args.gp_size)-2) } gout.write(gp_prefs%GPARGS) labels = "xyz" all_mins = [] all_max = [] for i in range(len(pc_indices)): gout.write('\nset %slabel "PC%s" font GPFONT_SIZE;\n'%(labels[i],pc_indices[i])) gout.write('\nset %stics %f font GPFONT_SIZE\n'%(labels[i],incr(pc_indices[i]))) the_min = 0.9*mines[pc_indices[i]] the_max = 1.1*maxes[pc_indices[i]] all_mins.append(the_min) all_max.append(the_max) if args.eigen_ratio: ef=open(evecname) eigens = ef.readline() if not re.search("^\s*#",eigens): sys.exit("Can't find eigen value in line <%s>"%eigens) eigs = map(float,re.findall("([-\d.]+)",eigens)) ratio = eigs[pc_indices[1]]/eigs[pc_indices[0]] ratiocmd ="\n\nset size ratio %f\n\n"%ratio gout.write(ratiocmd) ef.close() gout.write('\nset %srange [%f:%f]'%(labels[i],min(all_mins),max(all_max))) # Now do the plotting if args.gp_output: gout.write("\n\n\nset output \"%s\"\n\n"%args.gp_output) gout.write("\n%s\\\n %s"%(plot,plotcmd[0])) for p in plotcmd[1:]: gout.write(",\\\n %s"%p) for (ind,i,label) in labelcmd: gout.write(",\\\n \"\" index %d every 1::%d::%d using %s with labels notitle"%\ (ind,i,i,label)) if not args.gp_output: gout.write(";\n\n pause mouse\n\n") if args.epstopdf: gout.write('\nsystem "epstopdf %s";\n'%args.gp_output) if args.perspective: do_perspectives(gout,args) gout.close() args = setup() (evecname,phename,threed,gnu_term) = \ (args.evec,args.phe,args.three,args.gp_term) label_column = 1 plot = "splot" if threed else "plot" group = get_group(args) (gp_using_label, pc_indices) = getPCs(args) maxes = [-100]*100 mines = [100]*100 cmd = get_eigenvectors(group,gp_using_label,evecname) gnuplot_put(args,pc_indices,cmd)