view tools/rgenetics/rgQC.py @ 1:cdcb0ce84a1b

Uploaded

# oct 15 rpy replaced - temp fix until we get gnuplot working 
# rpy deprecated - replace with RRun
# fixes to run functional test! oct1 2009
# needed to expand our path with os.path.realpath to get newpath working with
# all the fancy pdfnup stuff
# and a fix to pruneld to write output to where it should be
# smallish data in test-data/smallwga in various forms
# python ../tools/rgenetics/rgQC.py -i smallwga -o smallwga -s smallwga/smallwga.html -p smallwga
# child files are deprecated and broken as at july 19 2009
# need to move them to the html file extrafiles path
# found lots of corner cases with some illumina data where cnv markers were
# included
# and where affection status was all missing !
# added links to tab files showing worst 1/keepfrac markers and subjects
# ross lazarus january 2008
#
# added named parameters
# to ensure no silly slippages if non required parameter in the most general case
# some potentially useful things here reusable in complex scripts
# with lots'o'html (TM)
# aug 17 2007 rml
#
# added marker and subject and parenting april 14 rml
# took a while to get the absolute paths right for all the file munging
# as of april 16 seems to work..
# getting galaxy to serve images in html reports is a little tricky
# we don't want QC reports to be dozens of individual files, so need
# to use the url /static/rg/... since galaxy's web server will happily serve images
# from there
# galaxy passes output files as relative paths
# these have to be munged by rgQC.py before calling this
# galaxy will pass in 2 file names - one for the log
# and one for the final html report
# of the form './database/files/dataset_66.dat'
# we need to be working in that directory so our plink output files are there
# so these have to be munged by rgQC.py before calling this
# note no ped file passed so had to remove the -l option
# for plinkParse.py that makes a heterozygosity report from the ped
# file - needs fixing...
# new: importing manhattan/qqplot plotter
# def doManQQ(input_fname,chrom_col,offset_col,pval_cols,title,grey,ctitle,outdir):
#    """ draw a qq for pvals and a manhattan plot if chrom/offset <> 0
#    contains some R scripts as text strings - we substitute defaults into the calls
#    to make them do our bidding - and save the resulting code for posterity
#    this can be called externally, I guess...for QC eg?
#    """
#
#    rcmd = '%s%s' % (rcode,rcode2 % (input_fname,chrom_col,offset_col,pval_cols,title,grey))
#    rlog,flist = RRun(rcmd=rcmd,title=ctitle,outdir=outdir)
#    return rlog,flist
  

from optparse import OptionParser

import sys,os,shutil, glob, math, subprocess, time, operator, random, tempfile, copy, string
from os.path import abspath
from rgutils import galhtmlprefix, galhtmlpostfix, RRun, timenow, plinke, rexe, runPlink, pruneLD
import rgManQQ

prog = os.path.split(sys.argv[0])[1]
vers = '0.4 april 2009 rml'
idjoiner = '_~_~_' # need something improbable..
# many of these may need fixing for a new install

myversion = vers
keepfrac = 20 # fraction to keep after sorting by each interesting value

missvals = {'0':'0','N':'N','-9':'-9','-':'-'} # fix me if these change!

mogresize = "x300" # this controls the width for jpeg thumbnails



            
def makePlots(markers=[],subjects=[],newfpath='.',basename='test',nbreaks='20',nup=3,height=10,width=8,rgbin=''):
    """
    marker rhead = ['snp','chrom','maf','a1','a2','missfrac',
    'p_hwe_all','logp_hwe_all','p_hwe_unaff','logp_hwe_unaff','N_Mendel']
    subject rhead = ['famId','iId','FracMiss','Mendel_errors','Ped_sex','SNP_sex','Status','Fest']
    """

        
    def rHist(plotme=[],outfname='',xlabname='',title='',basename='',nbreaks=50):
        """   rHist <- function(plotme,froot,plotname,title,mfname,nbreaks=50)
        # generic histogram and vertical boxplot in a 3:1 layout
        # returns the graphic file name for inclusion in the web page
        # broken out here for reuse
        # ross lazarus march 2007
        """
        screenmat = (1,2,1,2) # create a 2x2 cabvas
        widthlist = (80,20) # change to 4:1 ratio for histo and boxplot
        rpy.r.pdf( outfname, height , width  )
        #rpy.r.layout(rpy.r.matrix(rpy.r.c(1,1,1,2), 1, 4, byrow = True)) # 3 to 1 vertical plot
        m = rpy.r.matrix((1,1,1,2),nrow=1,ncol=4,byrow=True)
        # in R, m = matrix(c(1,2),nrow=1,ncol=2,byrow=T)
        rpy.r("layout(matrix(c(1,1,1,2),nrow=1,ncol=4,byrow=T))") # 4 to 1 vertical plot
        maint = 'QC for %s - %s' % (basename,title)
        rpy.r.hist(plotme,main=maint, xlab=xlabname,breaks=nbreaks,col="maroon",cex=0.8)
        rpy.r.boxplot(plotme,main='',col="maroon",outline=False)
        rpy.r.dev_off()

    def rCum(plotme=[],outfname='',xlabname='',title='',basename='',nbreaks=100):
        """
        Useful to see what various cutoffs yield - plot percentiles
        """
        n = len(plotme)
        maxveclen = 1000.0 # for reasonable pdf sizes!
        yvec = copy.copy(plotme)
        # arrives already in decending order of importance missingness or mendel count by subj or marker
        xvec = range(n)
        xvec = [100.0*(n-x)/n for x in xvec] # convert to centile
        # now have half a million markers eg - too many to plot all for a pdf - sample to get 10k or so points
        if n > maxveclen: # oversample part of the distribution
            always = min(1000,n/20) # oversample smaller of lowest few hundred items or 5%
            skip = int(n/maxveclen) # take 1 in skip to get about maxveclen points
            samplei = [i for i in range(n) if (i % skip == 0) or (i < always)] # always oversample first sorted values
            yvec = [yvec[i] for i in samplei] # always get first and last
            xvec = [xvec[i] for i in samplei] # always get first and last
        # need to supply the x axis or else rpy prints the freaking vector on the pdf - go figure
        rpy.r.pdf( outfname, height , width  )
        maint = 'QC for %s - %s' % (basename,title)
        rpy.r("par(lab=c(10,10,10))") # so our grid is denser than the default 5
        rpy.r.plot(xvec,yvec,type='p',main=maint, ylab=xlabname, xlab='Sample Percentile',col="maroon",cex=0.8)
        rpy.r.grid(nx = None, ny = None, col = "lightgray", lty = "dotted")
        rpy.r.dev_off()

    def rQQ(plotme=[], outfname='fname',title='title',xlabname='Sample',basename=''):
        """
        y is data for a qq plot and ends up on the x axis go figure
        if sampling, oversample low values - all the top 1% ?
        this version called with -log10 transformed hwe
        """
        nrows = len(plotme)
        fn = float(nrows)
        xvec = [-math.log10(x/fn) for x in range(1,(nrows+1))]
        mx = [0,math.log10(fn)] # if 1000, becomes 3 for the null line
        maxveclen = 3000
        yvec = copy.copy(plotme)
        if nrows > maxveclen:
            # now have half a million markers eg - too many to plot all for a pdf - sample to get 10k or so points
            # oversample part of the distribution
            always = min(1000,nrows/20) # oversample smaller of lowest few hundred items or 5%
            skip = int(nrows/float(maxveclen)) # take 1 in skip to get about maxveclen points
            samplei = [i for i in range(nrows) if (i < always) or (i % skip == 0)]
            # always oversample first sorted (here lowest) values
            yvec = [yvec[i] for i in samplei] # always get first and last
            xvec = [xvec[i] for i in samplei] # and sample xvec same way
            maint='Log QQ Plot (random %d of %d)' % (len(yvec),nrows)
        else:
            maint='Log QQ Plot(n=%d)' % (nrows)
        mx = [0,math.log10(fn)] # if 1000, becomes 3 for the null line
        ylab = '%s' % xlabname
        xlab = '-log10(Uniform 0-1)'
        # need to supply the x axis or else rpy prints the freaking vector on the pdf - go figure
        rpy.r.pdf( outfname, height , width  )
        rpy.r("par(lab=c(10,10,10))") # so our grid is denser than the default 5
        rpy.r.qqplot(xvec,yvec,xlab=xlab,ylab=ylab,main=maint,sub=title,pch=19,col="maroon",cex=0.8)
        rpy.r.points(mx,mx,type='l')
        rpy.r.grid(nx = None, ny = None, col = "lightgray", lty = "dotted")
        rpy.r.dev_off()

    def rMultiQQ(plotme = [],nsplits=5, outfname='fname',title='title',xlabname='Sample',basename=''):
        """
        data must contain p,x,y as data for a qq plot, quantiles of x and y axis used to create a
        grid of qq plots to show departure from null at extremes of data quality
        Need to plot qqplot(p,unif) where the p's come from one x and one y quantile
        and ends up on the x axis go figure
        if sampling, oversample low values - all the top 1% ?
        """
        data = copy.copy(plotme)
        nvals = len(data)
        stepsize = nvals/nsplits
        logstep = math.log10(stepsize) # so is 3 for steps of 1000
        quints = range(0,nvals,stepsize) # quintile cutpoints for each layer
        data.sort(key=itertools.itemgetter(1)) # into x order
        rpy.r.pdf( outfname, height , width  )
        rpy.r("par(mfrow = c(%d,%d))" % (nsplits,nsplits))
        yvec = [-math.log10(random.random()) for x in range(stepsize)]
        yvec.sort() # size of each step is expected range for xvec under null?!
        for rowstart in quints:
            rowend = rowstart + stepsize
            if nvals - rowend < stepsize: # finish last split
                rowend = nvals
            row = data[rowstart:rowend]
            row.sort(key=itertools.itemgetter(2)) # into y order
            for colstart in quints:
                colend = colstart + stepsize
                if nvals - colend < stepsize: # finish last split
                    colend = nvals
                cell = row[colstart:colend]
                xvec = [-math.log10(x[0]) for x in cell] # all the pvalues for this cell
                rpy.r.qqplot(xvec,yvec,xlab=xlab,ylab=ylab,pch=19,col="maroon",cex=0.8)
                rpy.r.points(c(0,logstep),c(0,logstep),type='l')
        rpy.r.dev_off()


    def rQQNorm(plotme=[], outfname='fname',title='title',xlabname='Sample',basename=''):
        """
        y is data for a qqnorm plot
        if sampling, oversample low values - all the top 1% ?
        """
        rangeunif = len(plotme)
        nunif = 1000
        maxveclen = 3000
        nrows = len(plotme)
        data = copy.copy(plotme)
        if nrows > maxveclen:
            # now have half a million markers eg - too many to plot all for a pdf - sample to get 10k or so points
            # oversample part of the distribution
            always = min(1000,nrows/20) # oversample smaller of lowest few hundred items or 5%
            skip = int((nrows-always)/float(maxveclen)) # take 1 in skip to get about maxveclen points
            samplei = [i for i in range(nrows) if (i % skip == 0) or (i < always)]
            # always oversample first sorted (here lowest) values
            yvec = [data[i] for i in samplei] # always get first and last
            maint='Log QQ Plot (random %d of %d)' % (len(yvec),nrows)
        else:
            yvec = data
            maint='Log QQ Plot(n=%d)' % (nrows)
        n = 1000
        ylab = '%s' % xlabname
        xlab = 'Normal'
        # need to supply the x axis or else rpy prints the freaking vector on the pdf - go figure
        rpy.r.pdf( outfname, height , width  )
        rpy.r("par(lab=c(10,10,10))") # so our grid is denser than the default 5
        rpy.r.qqnorm(yvec,xlab=xlab,ylab=ylab,main=maint,sub=title,pch=19,col="maroon",cex=0.8)
        rpy.r.grid(nx = None, ny = None, col = "lightgray", lty = "dotted")
        rpy.r.dev_off()

    def rMAFMissqq(plotme=[], outfname='fname',title='title',xlabname='Sample',basename=''):
        """
        layout qq plots for pvalues within rows of increasing MAF and columns of increasing missingness
        like the GAIN qc tools
        y is data for a qq plot and ends up on the x axis go figure
        if sampling, oversample low values - all the top 1% ?
        """
        rangeunif = len(plotme)
        nunif = 1000
        fn = float(rangeunif)
        xvec = [-math.log10(x/fn) for x in range(1,(rangeunif+1))]
        skip = max(int(rangeunif/fn),1)
        # force include last points
        mx = [0,math.log10(fn)] # if 1000, becomes 3 for the null line
        maxveclen = 2000
        nrows = len(plotme)
        data = copy.copy(plotme)
        data.sort() # low to high - oversample low values
        if nrows > maxveclen:
            # now have half a million markers eg - too many to plot all for a pdf - sample to get 10k or so points
            # oversample part of the distribution
            always = min(1000,nrows/20) # oversample smaller of lowest few hundred items or 5%
            skip = int(nrows/float(maxveclen)) # take 1 in skip to get about maxveclen points
            samplei = [i for i in range(nrows) if (i % skip == 0) or (i < always)]
            # always oversample first sorted (here lowest) values
            yvec = [data[i] for i in samplei] # always get first and last
            xvec = [xvec[i] for i in samplei] # and sample xvec same way
            maint='Log QQ Plot (random %d of %d)' % (len(yvec),nrows)
        else:
            yvec = data
            maint='Log QQ Plot(n=%d)' % (nrows)
        n = 1000
        mx = [0,log10(fn)] # if 1000, becomes 3 for the null line
        ylab = '%s' % xlabname
        xlab = '-log10(Uniform 0-1)'
        # need to supply the x axis or else rpy prints the freaking vector on the pdf - go figure
        rpy.r.pdf( outfname, height , width  )
        rpy.r("par(lab=c(10,10,10))") # so our grid is denser than the default 5
        rpy.r.qqplot(xvec,yvec,xlab=xlab,ylab=ylab,main=maint,sub=title,pch=19,col="maroon",cex=0.8)
        rpy.r.points(mx,mx,type='l')
        rpy.r.grid(nx = None, ny = None, col = "lightgray", lty = "dotted")
        rpy.r.dev_off()


    fdsto,stofile = tempfile.mkstemp()
    sto = open(stofile,'w')
    import rpy # delay to avoid rpy stdout chatter replacing galaxy file blurb
    mog = 'mogrify'
    pdfnup = 'pdfnup'
    pdfjoin = 'pdfjoin'
    shead = subjects.pop(0) # get rid of head
    mhead = markers.pop(0)
    maf = mhead.index('maf')
    missfrac = mhead.index('missfrac')
    logphweall = mhead.index('logp_hwe_all')
    logphweunaff = mhead.index('logp_hwe_unaff')
    # check for at least some unaffected rml june 2009
    m_mendel = mhead.index('N_Mendel')
    fracmiss = shead.index('FracMiss')
    s_mendel = shead.index('Mendel_errors')
    s_het = shead.index('F_Stat')
    params = {}
    hweres = [float(x[logphweunaff]) for x in markers if len(x[logphweunaff]) >= logphweunaff
         and x[logphweunaff].upper() <> 'NA']
    if len(hweres) <> 0:
        xs = [logphweunaff, missfrac, maf, m_mendel, fracmiss, s_mendel, s_het]
        # plot for each of these cols
    else: # try hwe all instead - maybe no affection status available
        xs = [logphweall, missfrac, maf, m_mendel, fracmiss, s_mendel, s_het]
    ordplotme = [1,1,1,1,1,1,1] # ordered plots for everything!
    oreverseme = [1,1,0,1,1,1,0] # so larger values are oversampled
    qqplotme = [1,0,0,0,0,0,0] #
    qnplotme = [0,0,0,0,0,0,1] #
    nplots = len(xs)
    xlabnames = ['log(p) HWE (unaff)', 'Missing Rate: Markers', 'Minor Allele Frequency',
                 'Marker Mendel Error Count', 'Missing Rate: Subjects',
                 'Subject Mendel Error Count','Subject Inbreeding (het) F Statistic']
    plotnames = ['logphweunaff', 'missfrac', 'maf', 'm_mendel', 'fracmiss', 's_mendel','s_het']
    ploturls = ['%s_%s.pdf' % (basename,x) for x in plotnames] # real plotnames
    ordplotnames = ['%s_cum' % x for x in plotnames]
    ordploturls = ['%s_%s.pdf' % (basename,x) for x in ordplotnames] # real plotnames
    outfnames = [os.path.join(newfpath,ploturls[x]) for x in range(nplots)]
    ordoutfnames = [os.path.join(newfpath,ordploturls[x]) for x in range(nplots)]
    datasources = [markers,markers,markers,markers,subjects,subjects,subjects] # use this table
    titles = ["Marker HWE","Marker Missing Genotype", "Marker MAF","Marker Mendel",
        "Subject Missing Genotype","Subject Mendel",'Subject F Statistic']
    html = []
    pdflist = []
    for n,column in enumerate(xs):
        dat = [float(x[column]) for x in datasources[n] if len(x) >= column
               and x[column][:2].upper() <> 'NA'] # plink gives both!
        if sum(dat) <> 0: # eg nada for mendel if case control?
            rHist(plotme=dat,outfname=outfnames[n],xlabname=xlabnames[n],
              title=titles[n],basename=basename,nbreaks=nbreaks)
            row = [titles[n],ploturls[n],outfnames[n] ]
            html.append(row)
            pdflist.append(outfnames[n])
            if ordplotme[n]: # for missingness, hwe - plots to see where cutoffs will end up
                otitle = 'Ranked %s' % titles[n]
                dat.sort()
                if oreverseme[n]:
                    dat.reverse()
                rCum(plotme=dat,outfname=ordoutfnames[n],xlabname='Ordered %s' % xlabnames[n],
                  title=otitle,basename=basename,nbreaks=1000)
                row = [otitle,ordploturls[n],ordoutfnames[n]]
                html.append(row)
                pdflist.append(ordoutfnames[n])
            if qqplotme[n]: #
                otitle = 'LogQQ plot %s' % titles[n]
                dat.sort()
                dat.reverse()
                ofn = os.path.split(ordoutfnames[n])
                ofn = os.path.join(ofn[0],'QQ%s' % ofn[1])
                ofu = os.path.split(ordploturls[n])
                ofu = os.path.join(ofu[0],'QQ%s' % ofu[1])
                rQQ(plotme=dat,outfname=ofn,xlabname='QQ %s' % xlabnames[n],
                  title=otitle,basename=basename)
                row = [otitle,ofu,ofn]
                html.append(row)
                pdflist.append(ofn)
            elif qnplotme[n]:
                otitle = 'F Statistic %s' % titles[n]
                dat.sort()
                dat.reverse()
                ofn = os.path.split(ordoutfnames[n])
                ofn = os.path.join(ofn[0],'FQNorm%s' % ofn[1])
                ofu = os.path.split(ordploturls[n])
                ofu = os.path.join(ofu[0],'FQNorm%s' % ofu[1])
                rQQNorm(plotme=dat,outfname=ofn,xlabname='F QNorm %s' % xlabnames[n],
                  title=otitle,basename=basename)
                row = [otitle,ofu,ofn]
                html.append(row)
                pdflist.append(ofn)
        else:
            print '#$# no data for # %d - %s, data[:10]=%s' % (n,titles[n],dat[:10])
    if nup>0:
        # pdfjoin --outfile chr1test.pdf `ls database/files/dataset_396_files/*.pdf`
        # pdfnup chr1test.pdf --nup 3x3 --frame true --outfile chr1test3.pdf
        filestojoin = ' '.join(pdflist) # all the file names so far
        afname = '%s_All_Paged.pdf' % (basename)
        fullafname = os.path.join(newfpath,afname)
        expl = 'All %s QC Plots joined into a single pdf' % basename
        vcl = '%s %s --outfile %s ' % (pdfjoin,filestojoin, fullafname)
        # make single page pdf
        x=subprocess.Popen(vcl,shell=True,cwd=newfpath,stderr=sto,stdout=sto)
        retval = x.wait()
        row = [expl,afname,fullafname]
        html.insert(0,row) # last rather than second
        nfname = '%s_All_%dx%d.pdf' % (basename,nup,nup)
        fullnfname = os.path.join(newfpath,nfname)
        expl = 'All %s QC Plots %d by %d to a page' % (basename,nup,nup)
        vcl = '%s %s --nup %dx%d --frame true --outfile %s' % (pdfnup,afname,nup,nup,fullnfname)
        # make thumbnail images
        x=subprocess.Popen(vcl,shell=True,cwd=newfpath,stderr=sto,stdout=sto)
        retval = x.wait()
        row = [expl,nfname,fullnfname]
        html.insert(1,row) # this goes second
    vcl = '%s -format jpg -resize %s %s' % (mog, mogresize, os.path.join(newfpath,'*.pdf'))
    # make thumbnail images
    x=subprocess.Popen(vcl,shell=True,cwd=newfpath,stderr=sto,stdout=sto)
    retval = x.wait()
    sto.close()
    cruft = open(stofile,'r').readlines()
    return html,cruft # elements for an ordered list of urls or whatever..  


def RmakePlots(markers=[],subjects=[],newfpath='.',basename='test',nbreaks='100',nup=3,height=8,width=10,rexe=''):
    """
    nice try but the R scripts are huge and take forever to run if there's a lot of data
    marker rhead = ['snp','chrom','maf','a1','a2','missfrac',
    'p_hwe_all','logp_hwe_all','p_hwe_unaff','logp_hwe_unaff','N_Mendel']
    subject rhead = ['famId','iId','FracMiss','Mendel_errors','Ped_sex','SNP_sex','Status','Fest']
    """
    colour = "maroon"
        
    def rHist(plotme='',outfname='',xlabname='',title='',basename='',nbreaks=nbreaks):
        """   rHist <- function(plotme,froot,plotname,title,mfname,nbreaks=50)
        # generic histogram and vertical boxplot in a 3:1 layout
        # returns the graphic file name for inclusion in the web page
        # broken out here for reuse
        # ross lazarus march 2007
        """
        R = []
        maint = 'QC for %s - %s' % (basename,title)
        screenmat = (1,2,1,2) # create a 2x2 canvas
        widthlist = (80,20) # change to 4:1 ratio for histo and boxplot
        R.append('pdf("%s",h=%d,w=%d)' % (outfname,height,width))
        R.append("layout(matrix(c(1,1,1,2),nrow=1,ncol=4,byrow=T))")
        R.append("plotme = read.table(file='%s',head=F,sep='\t')" % plotme)
        R.append('hist(plotme, main="%s",xlab="%s",breaks=%d,col="%s")' % (maint,xlabname,nbreaks,colour))
        R.append('boxplot(plotme,main="",col="%s",outline=F)' % (colour) )
        R.append('dev.off()')
        return R
        
    def rCum(plotme='',outfname='',xlabname='',title='',basename='',nbreaks=100):
        """
        Useful to see what various cutoffs yield - plot percentiles
        """
        R = []
        n = len(plotme)
        R.append("plotme = read.table(file='%s',head=T,sep='\t')" % plotme)
        # arrives already in decending order of importance missingness or mendel count by subj or marker
        maint = 'QC for %s - %s' % (basename,title)
        R.append('pdf("%s",h=%d,w=%d)' % (outfname,height,width))
        R.append("par(lab=c(10,10,10))")
        R.append('plot(plotme$xvec,plotme$yvec,type="p",main="%s",ylab="%s",xlab="Sample Percentile",col="%s")' % (maint,xlabname,colour))
        R.append('dev.off()')
        return R

    def rQQ(plotme='', outfname='fname',title='title',xlabname='Sample',basename=''):
        """
        y is data for a qq plot and ends up on the x axis go figure
        if sampling, oversample low values - all the top 1% ?
        this version called with -log10 transformed hwe
        """
        R = []
        nrows = len(plotme)
        fn = float(nrows)
        xvec = [-math.log10(x/fn) for x in range(1,(nrows+1))]
        #mx = [0,math.log10(fn)] # if 1000, becomes 3 for the null line
        maxveclen = 3000
        yvec = copy.copy(plotme)
        if nrows > maxveclen:
            # now have half a million markers eg - too many to plot all for a pdf - sample to get 10k or so points
            # oversample part of the distribution
            always = min(1000,nrows/20) # oversample smaller of lowest few hundred items or 5%
            skip = int(nrows/float(maxveclen)) # take 1 in skip to get about maxveclen points
            if skip < 2:
                skip = 2
            samplei = [i for i in range(nrows) if (i < always) or (i % skip == 0)]
            # always oversample first sorted (here lowest) values
            yvec = [yvec[i] for i in samplei] # always get first and last
            xvec = [xvec[i] for i in samplei] # and sample xvec same way
            maint='Log QQ Plot (random %d of %d)' % (len(yvec),nrows)
        else:
            maint='Log QQ Plot(n=%d)' % (nrows)
        mx = [0,math.log10(fn)] # if 1000, becomes 3 for the null line
        ylab = '%s' % xlabname
        xlab = '-log10(Uniform 0-1)'
        # need to supply the x axis or else rpy prints the freaking vector on the pdf - go figure
        x = ['%f' % x for x in xvec]
        R.append('xvec = c(%s)' % ','.join(x))
        y = ['%f' % x for x in yvec]
        R.append('yvec = c(%s)' % ','.join(y))
        R.append('mx = c(0,%f)' % (math.log10(fn)))
        R.append('pdf("%s",h=%d,w=%d)' % (outfname,height,width))
        R.append("par(lab=c(10,10,10))")
        R.append('qqplot(xvec,yvec,xlab="%s",ylab="%s",main="%s",sub="%s",pch=19,col="%s",cex=0.8)' % (xlab,ylab,maint,title,colour))
        R.append('points(mx,mx,type="l")')
        R.append('grid(col="lightgray",lty="dotted")')
        R.append('dev.off()')
        return R

    def rMultiQQ(plotme = '',nsplits=5, outfname='fname',title='title',xlabname='Sample',basename=''):
        """
        data must contain p,x,y as data for a qq plot, quantiles of x and y axis used to create a
        grid of qq plots to show departure from null at extremes of data quality
        Need to plot qqplot(p,unif) where the p's come from one x and one y quantile
        and ends up on the x axis go figure
        if sampling, oversample low values - all the top 1% ?
        """
        data = copy.copy(plotme)
        nvals = len(data)
        stepsize = nvals/nsplits
        logstep = math.log10(stepsize) # so is 3 for steps of 1000
        R.append('mx = c(0,%f)' % logstep)
        quints = range(0,nvals,stepsize) # quintile cutpoints for each layer
        data.sort(key=itertools.itemgetter(1)) # into x order
        #rpy.r.pdf( outfname, h , w  )
        #rpy.r("par(mfrow = c(%d,%d))" % (nsplits,nsplits))
        R.append('pdf("%s",h=%d,w=%d)' % (outfname,height,width))
        yvec = [-math.log10(random.random()) for x in range(stepsize)]
        yvec.sort() # size of each step is expected range for xvec under null?!
        y = ['%f' % x for x in yvec]
        R.append('yvec = c(%s)' % ','.join(y))
        for rowstart in quints:
            rowend = rowstart + stepsize
            if nvals - rowend < stepsize: # finish last split
                rowend = nvals
            row = data[rowstart:rowend]
            row.sort(key=itertools.itemgetter(2)) # into y order
            for colstart in quints:
                colend = colstart + stepsize
                if nvals - colend < stepsize: # finish last split
                    colend = nvals
                cell = row[colstart:colend]
                xvec = [-math.log10(x[0]) for x in cell] # all the pvalues for this cell
                x = ['%f' % x for x in xvec]
                R.append('xvec = c(%s)' % ','.join(x))
                R.append('qqplot(xvec,yvec,xlab="%s",ylab="%s",main="%s",sub="%s",pch=19,col="%s",cex=0.8)' % (xlab,ylab,maint,title,colour))
                R.append('points(mx,mx,type="l")')
                R.append('grid(col="lightgray",lty="dotted")')
                #rpy.r.qqplot(xvec,yvec,xlab=xlab,ylab=ylab,pch=19,col="maroon",cex=0.8)
                #rpy.r.points(c(0,logstep),c(0,logstep),type='l')
        R.append('dev.off()')
        #rpy.r.dev_off()
        return R


    def rQQNorm(plotme=[], outfname='fname',title='title',xlabname='Sample',basename=''):
        """
        y is data for a qqnorm plot
        if sampling, oversample low values - all the top 1% ?
        """
        rangeunif = len(plotme)
        nunif = 1000
        maxveclen = 3000
        nrows = len(plotme)
        data = copy.copy(plotme)
        if nrows > maxveclen:
            # now have half a million markers eg - too many to plot all for a pdf - sample to get 10k or so points
            # oversample part of the distribution
            always = min(1000,nrows/20) # oversample smaller of lowest few hundred items or 5%
            skip = int((nrows-always)/float(maxveclen)) # take 1 in skip to get about maxveclen points
            samplei = [i for i in range(nrows) if (i % skip == 0) or (i < always)]
            # always oversample first sorted (here lowest) values
            yvec = [data[i] for i in samplei] # always get first and last
            maint='Log QQ Plot (random %d of %d)' % (len(yvec),nrows)
        else:
            yvec = data
            maint='Log QQ Plot(n=%d)' % (nrows)
        n = 1000
        ylab = '%s' % xlabname
        xlab = 'Normal'
        # need to supply the x axis or else rpy prints the freaking vector on the pdf - go figure
        #rpy.r.pdf( outfname, h , w  )
        #rpy.r("par(lab=c(10,10,10))") # so our grid is denser than the default 5
        #rpy.r.qqnorm(yvec,xlab=xlab,ylab=ylab,main=maint,sub=title,pch=19,col="maroon",cex=0.8)
        #rpy.r.grid(nx = None, ny = None, col = "lightgray", lty = "dotted")
        #rpy.r.dev_off()
        y = ['%f' % x for x in yvec]
        R.append('yvec = c(%s)' % ','.join(y))
        R.append('pdf("%s",h=%d,w=%d)' % (outfname,height,width))
        R.append("par(lab=c(10,10,10))")
        R.append('qqnorm(yvec,xlab="%s",ylab="%s",main="%s",sub="%s",pch=19,col="%s",cex=0.8)' % (xlab,ylab,maint,title,colour))
        R.append('grid(col="lightgray",lty="dotted")')
        R.append('dev.off()')
        return R

    def rMAFMissqq(plotme=[], outfname='fname',title='title',xlabname='Sample',basename=''):
        """
        layout qq plots for pvalues within rows of increasing MAF and columns of increasing missingness
        like the GAIN qc tools
        y is data for a qq plot and ends up on the x axis go figure
        if sampling, oversample low values - all the top 1% ?
        """
        rangeunif = len(plotme)
        nunif = 1000
        fn = float(rangeunif)
        xvec = [-math.log10(x/fn) for x in range(1,(rangeunif+1))]
        skip = max(int(rangeunif/fn),1)
        # force include last points
        mx = [0,math.log10(fn)] # if 1000, becomes 3 for the null line
        maxveclen = 2000
        nrows = len(plotme)
        data = copy.copy(plotme)
        data.sort() # low to high - oversample low values
        if nrows > maxveclen:
            # now have half a million markers eg - too many to plot all for a pdf - sample to get 10k or so points
            # oversample part of the distribution
            always = min(1000,nrows/20) # oversample smaller of lowest few hundred items or 5%
            skip = int(nrows/float(maxveclen)) # take 1 in skip to get about maxveclen points
            samplei = [i for i in range(nrows) if (i % skip == 0) or (i < always)]
            # always oversample first sorted (here lowest) values
            yvec = [data[i] for i in samplei] # always get first and last
            xvec = [xvec[i] for i in samplei] # and sample xvec same way
            maint='Log QQ Plot (random %d of %d)' % (len(yvec),nrows)
        else:
            yvec = data
            maint='Log QQ Plot(n=%d)' % (nrows)
        n = 1000
        mx = [0,log10(fn)] # if 1000, becomes 3 for the null line
        ylab = '%s' % xlabname
        xlab = '-log10(Uniform 0-1)'
        R.append('mx = c(0,%f)' % (math.log10(fn)))
        x = ['%f' % x for x in xvec]
        R.append('xvec = c(%s)' % ','.join(x))
        y = ['%f' % x for x in yvec]
        R.append('yvec = c(%s)' % ','.join(y))
        R.append('pdf("%s",h=%d,w=%d)' % (outfname,height,width))
        R.append("par(lab=c(10,10,10))")
        R.append('qqplot(xvec,yvec,xlab="%s",ylab="%s",main="%s",sub="%s",pch=19,col="%s",cex=0.8)' % (xlab,ylab,maint,title,colour))
        R.append('points(mx,mx,type="l")')
        R.append('grid(col="lightgray",lty="dotted")')
        R.append('dev.off()')


    shead = subjects.pop(0) # get rid of head
    mhead = markers.pop(0)
    maf = mhead.index('maf')
    missfrac = mhead.index('missfrac')
    logphweall = mhead.index('logp_hwe_all')
    logphweunaff = mhead.index('logp_hwe_unaff')
    # check for at least some unaffected rml june 2009
    m_mendel = mhead.index('N_Mendel')
    fracmiss = shead.index('FracMiss')
    s_mendel = shead.index('Mendel_errors')
    s_het = shead.index('F_Stat')
    params = {}
    h = [float(x[logphweunaff]) for x in markers if len(x[logphweunaff]) >= logphweunaff
         and x[logphweunaff].upper() <> 'NA']
    if len(h) <> 0:
        xs = [logphweunaff, missfrac, maf, m_mendel, fracmiss, s_mendel, s_het]
        # plot for each of these cols
    else: # try hwe all instead - maybe no affection status available
        xs = [logphweall, missfrac, maf, m_mendel, fracmiss, s_mendel, s_het]
    ordplotme = [1,1,1,1,1,1,1] # ordered plots for everything!
    oreverseme = [1,1,0,1,1,1,0] # so larger values are oversampled
    qqplotme = [1,0,0,0,0,0,0] #
    qnplotme = [0,0,0,0,0,0,1] #
    nplots = len(xs)
    xlabnames = ['log(p) HWE (unaff)', 'Missing Rate: Markers', 'Minor Allele Frequency',
                 'Marker Mendel Error Count', 'Missing Rate: Subjects',
                 'Subject Mendel Error Count','Subject Inbreeding (het) F Statistic']
    plotnames = ['logphweunaff', 'missfrac', 'maf', 'm_mendel', 'fracmiss', 's_mendel','s_het']
    ploturls = ['%s_%s.pdf' % (basename,x) for x in plotnames] # real plotnames
    ordplotnames = ['%s_cum' % x for x in plotnames]
    ordploturls = ['%s_%s.pdf' % (basename,x) for x in ordplotnames] # real plotnames
    outfnames = [os.path.join(newfpath,ploturls[x]) for x in range(nplots)]
    ordoutfnames = [os.path.join(newfpath,ordploturls[x]) for x in range(nplots)]
    datasources = [markers,markers,markers,markers,subjects,subjects,subjects] # use this table
    titles = ["Marker HWE","Marker Missing Genotype", "Marker MAF","Marker Mendel",
        "Subject Missing Genotype","Subject Mendel",'Subject F Statistic']
    html = []
    pdflist = []
    R = []
    for n,column in enumerate(xs):
        dfn = '%d_%s.txt' % (n,titles[n])
        dfilepath = os.path.join(newfpath,dfn)
        dat = [float(x[column]) for x in datasources[n] if len(x) >= column
               and x[column][:2].upper() <> 'NA'] # plink gives both!
        if sum(dat) <> 0: # eg nada for mendel if case control?
            plotme = file(dfilepath,'w')
            plotme.write('\n'.join(['%f' % x for x in dat])) # pass as a file - copout!
            tR = rHist(plotme=dfilepath,outfname=outfnames[n],xlabname=xlabnames[n],
              title=titles[n],basename=basename,nbreaks=nbreaks)
            R += tR
            row = [titles[n],ploturls[n],outfnames[n] ]
            html.append(row)
            pdflist.append(outfnames[n])
            if ordplotme[n]: # for missingness, hwe - plots to see where cutoffs will end up
                otitle = 'Ranked %s' % titles[n]
                dat.sort()
                if oreverseme[n]:
                    dat.reverse()
                    ndat = len(dat)
                    xvec = range(ndat)
                    xvec = [100.0*(n-x)/n for x in xvec] # convert to centile
                    # now have half a million markers eg - too many to plot all for a pdf - sample to get 10k or so points
                    maxveclen = 1000.0 # for reasonable pdf sizes!
                    if ndat > maxveclen: # oversample part of the distribution
                        always = min(1000,ndat/20) # oversample smaller of lowest few hundred items or 5%
                        skip = int(ndat/maxveclen) # take 1 in skip to get about maxveclen points
                        samplei = [i for i in range(ndat) if (i % skip == 0) or (i < always)] # always oversample first sorted values
                        yvec = [yvec[i] for i in samplei] # always get first and last
                        xvec = [xvec[i] for i in samplei] # always get first and last
                        plotme = file(dfilepath,'w')
                        plotme.write('xvec\tyvec\n')
                        plotme.write('\n'.join(['%f\t%f' % (xvec[i],y) for y in yvec])) # pass as a file - copout!
                tR = rCum(plotme=dat,outfname=ordoutfnames[n],xlabname='Ordered %s' % xlabnames[n],
                  title=otitle,basename=basename,nbreaks=nbreaks)
                R += tR
                row = [otitle,ordploturls[n],ordoutfnames[n]]
                html.append(row)
                pdflist.append(ordoutfnames[n])
            if qqplotme[n]: #
                otitle = 'LogQQ plot %s' % titles[n]
                dat.sort()
                dat.reverse()
                ofn = os.path.split(ordoutfnames[n])
                ofn = os.path.join(ofn[0],'QQ%s' % ofn[1])
                ofu = os.path.split(ordploturls[n])
                ofu = os.path.join(ofu[0],'QQ%s' % ofu[1])
                tR = rQQ(plotme=dat,outfname=ofn,xlabname='QQ %s' % xlabnames[n],
                  title=otitle,basename=basename)
                R += tR
                row = [otitle,ofu,ofn]
                html.append(row)
                pdflist.append(ofn)
            elif qnplotme[n]:
                otitle = 'F Statistic %s' % titles[n]
                dat.sort()
                dat.reverse()
                ofn = os.path.split(ordoutfnames[n])
                ofn = os.path.join(ofn[0],'FQNorm%s' % ofn[1])
                ofu = os.path.split(ordploturls[n])
                ofu = os.path.join(ofu[0],'FQNorm%s' % ofu[1])
                tR = rQQNorm(plotme=dat,outfname=ofn,xlabname='F QNorm %s' % xlabnames[n],
                  title=otitle,basename=basename)
                R += tR
                row = [otitle,ofu,ofn]
                html.append(row)
                pdflist.append(ofn)
        else:
            print '#$# no data for # %d - %s, data[:10]=%s' % (n,titles[n],dat[:10])
    rlog,flist = RRun(rcmd=R,title='makeQCplots',outdir=newfpath)
    if nup>0:
        # pdfjoin --outfile chr1test.pdf `ls database/files/dataset_396_files/*.pdf`
        # pdfnup chr1test.pdf --nup 3x3 --frame true --outfile chr1test3.pdf
        filestojoin = ' '.join(pdflist) # all the file names so far
        afname = '%s_All_Paged.pdf' % (basename)
        fullafname = os.path.join(newfpath,afname)
        expl = 'All %s QC Plots joined into a single pdf' % basename
        vcl = 'pdfjoin %s --outfile %s ' % (filestojoin, fullafname)
        # make single page pdf
        x=subprocess.Popen(vcl,shell=True,cwd=newfpath)
        retval = x.wait()
        row = [expl,afname,fullafname]
        html.insert(0,row) # last rather than second
        nfname = '%s_All_%dx%d.pdf' % (basename,nup,nup)
        fullnfname = os.path.join(newfpath,nfname)
        expl = 'All %s QC Plots %d by %d to a page' % (basename,nup,nup)
        vcl = 'pdfnup %s --nup %dx%d --frame true --outfile %s' % (afname,nup,nup,fullnfname)
        # make thumbnail images
        x=subprocess.Popen(vcl,shell=True,cwd=newfpath)
        retval = x.wait()
        row = [expl,nfname,fullnfname]
        html.insert(1,row) # this goes second
    vcl = 'mogrify -format jpg -resize %s %s' % (mogresize, os.path.join(newfpath,'*.pdf'))
    # make thumbnail images
    x=subprocess.Popen(vcl,shell=True,cwd=newfpath)
    retval = x.wait()
    return html # elements for an ordered list of urls or whatever..

def countHet(bedf='fakeped_500000',linkageped=True,froot='fake500k',outfname="ahetf",logf='rgQC.log'):
    """
    NO LONGER USED - historical interest
    count het loci for each subject to look for outliers = ? contamination
    assume ped file is linkage format
    need to make a ped file from the bed file we were passed
    """
    vcl = [plinke,'--bfile',bedf,'--recode','--out','%s_recode' % froot] # write a recoded ped file from the real .bed file
    p=subprocess.Popen(' '.join(vcl),shell=True)
    retval = p.wait()
    f = open('%s_recode.recode.ped' % froot,'r')
    if not linkageped:
        head = f.next() # throw away header
    hets = [] # simple count of het loci per subject. Expect poisson?
    n = 1
    for l in f:
        n += 1
        ll = l.strip().split()
        if len(ll) > 6:
            iid = idjoiner.join(ll[:2]) # fam_iid
            gender = ll[4]
            alleles = ll[6:]
            nallele = len(alleles)
            nhet = 0
            for i in range(nallele/2):
                a1=alleles[2*i]
                a2=alleles[2*i+1]
                if alleles[2*i] <> alleles[2*i+1]: # must be het
                    if not missvals.get(a1,None) and not missvals.get(a2,None):
                        nhet += 1
            hets.append((nhet,iid,gender)) # for sorting later
    f.close()
    hets.sort()
    hets.reverse() # biggest nhet now on top
    f = open(outfname ,'w')
    res = ['%d\t%s\t%s' % (x) for x in hets] # I love list comprehensions
    res.insert(0,'nhetloci\tfamid_iid\tgender')
    res.append('')
    f.write('\n'.join(res))
    f.close()



def subjectRep(froot='cleantest',outfname="srep",newfpath='.',logf = None):
    """by subject (missingness = .imiss, mendel = .imendel)
    assume replicates have an underscore in family id for
    hapmap testing
    For sorting, we need floats and integers
    """
    isexfile = '%s.sexcheck' % froot
    imissfile = '%s.imiss' % froot
    imendfile = '%s.imendel' % froot
    ihetfile = '%s.het' % froot
    logf.write('## subject reports starting at %s\n' % timenow())
    outfile = os.path.join(newfpath,outfname)
    idlist = []
    imissdict = {}
    imenddict = {}
    isexdict = {}
    ihetdict = {}
    Tops = {}
    Tnames = ['Ranked Subject Missing Genotype', 'Ranked Subject Mendel',
              'Ranked Sex check', 'Ranked Inbreeding (het) F statistic']
    Tsorts = [2,3,6,8]
    Treverse = [True,True,True,False] # so first values are worser
    #rhead = ['famId','iId','FracMiss','Mendel_errors','Ped_sex','SNP_sex','Status','Fest']
    ##              FID            IID MISS_PHENO   N_MISS   N_GENO   F_MISS
    ##  1552042370_A   1552042370_A          N     5480   549883 0.009966
    ##  1552042410_A   1552042410_A          N     1638   549883 0.002979
 
    # ------------------missing--------------------
    # imiss has FID  IID MISS_PHENO N_MISS  F_MISS
    # we want F_MISS
    try:
        f = open(imissfile,'r')
    except:
        logf.write('# file %s is missing - talk about irony\n' % imissfile)
        f = None
    if f:
        for n,line in enumerate(f):
            ll = line.strip().split()
            if n == 0:
                head = [x.upper() for x in ll] # expect above                
                fidpos = head.index('FID')
                iidpos = head.index('IID')
                fpos = head.index('F_MISS')
            elif len(ll) >= fpos: # full line
                fid = ll[fidpos]
                #if fid.find('_') == -1: # not replicate! - icondb ids have these...
                iid = ll[iidpos]
                fmiss = ll[fpos]
                id = '%s%s%s' % (fid,idjoiner,iid)
                imissdict[id] = fmiss
                idlist.append(id)
        f.close()
    logf.write('## imissfile %s contained %d ids\n' % (imissfile,len(idlist)))
    # ------------------mend-------------------
    # *.imendel has FID  IID   N
    # we want N
    gotmend = True
    try:
        f = open(imendfile,'r')
    except:
        gotmend = False
        for id in idlist:
            imenddict[id] = '0'
    if gotmend:
        for n,line in enumerate(f):
            ll = line.strip().split()
            if n == 0:
                head = [x.upper() for x in ll] # expect above                
                npos = head.index('N')
                fidpos = head.index('FID')
                iidpos = head.index('IID')
            elif len(ll) >= npos: # full line
                fid = ll[fidpos]
                iid = ll[iidpos]
                id = '%s%s%s' % (fid,idjoiner,iid)
                nmend = ll[npos]
                imenddict[id] = nmend
        f.close()
    else:
        logf.write('## error No %s file - assuming not family data\n' % imendfile)
    # ------------------sex check------------------
    #[rerla@hg fresh]$ head /home/rerla/fresh/database/files/dataset_978_files/CAMP2007Dirty.sexcheck
    # sexcheck has FID IID PEDSEX SNPSEX STATUS F
    ##
    ##     FID     Family ID
    ##     IID     Individual ID
    ##     PEDSEX  Sex as determined in pedigree file (1=male, 2=female)
    ##     SNPSEX  Sex as determined by X chromosome
    ##     STATUS  Displays "PROBLEM" or "OK" for each individual
    ##     F       The actual X chromosome inbreeding (homozygosity) estimate
    ##
    ##    A PROBLEM arises if the two sexes do not match, or if the SNP data or pedigree data are
    ##    ambiguous with regard to sex.
    ##    A male call is made if F is more than 0.8; a femle call is made if F is less than 0.2.
    try:
        f = open(isexfile,'r')
        got_sexcheck = 1
    except:
        got_sexcheck = 0
    if got_sexcheck:
        for n,line in enumerate(f):
            ll = line.strip().split()
            if n == 0:
                head = [x.upper() for x in ll] # expect above                
                fidpos = head.index('FID')
                iidpos = head.index('IID')
                pedsexpos = head.index('PEDSEX')
                snpsexpos = head.index('SNPSEX')
                statuspos = head.index('STATUS')
                fpos = head.index('F')
            elif len(ll) >= fpos: # full line
                fid = ll[fidpos]
                iid = ll[iidpos]
                fest = ll[fpos]
                pedsex = ll[pedsexpos]
                snpsex = ll[snpsexpos]
                stat = ll[statuspos]
                id = '%s%s%s' % (fid,idjoiner,iid)
                isexdict[id] = (pedsex,snpsex,stat,fest)
        f.close()
    else:
        # this only happens if there are no subjects!
        logf.write('No %s file - assuming no sex errors' % isexfile)
    ##
    ##    FID  IID       O(HOM)       E(HOM)        N(NM)            F
    ##    457    2       490665    4.928e+05       722154    -0.009096
    ##    457    3       464519     4.85e+05       710986      -0.0908
    ##   1037    2       461632    4.856e+05       712025       -0.106
    ##   1037    1       491845    4.906e+05       719353     0.005577
    try:
        f = open(ihetfile,'r')
    except:
        f = None
        logf.write('## No %s file - did we run --het in plink?\n' % ihetfile)
    if f:
        for i,line in enumerate(f):
            ll = line.strip().split()
            if i == 0:
                head = [x.upper() for x in ll] # expect above                
                fidpos = head.index('FID')
                iidpos = head.index('IID')
                fpos = head.index('F')
                n = 0
            elif len(ll) >= fpos: # full line
                fid = ll[fidpos]            
                iid = ll[iidpos]
                fhet = ll[fpos]
                id = '%s%s%s' % (fid,idjoiner,iid)
                ihetdict[id] = fhet
        f.close()      # now assemble and output result list
    rhead = ['famId','iId','FracMiss','Mendel_errors','Ped_sex','SNP_sex','Status','XHomEst','F_Stat']
    res = []
    fres = [] # floats for sorting
    for id in idlist: # for each snp in found order
        fid,iid = id.split(idjoiner) # recover keys
        f_missing = imissdict.get(id,'0.0')
        nmend = imenddict.get(id,'0')
        (pedsex,snpsex,status,fest) = isexdict.get(id,('0','0','0','0.0'))
        fhet = ihetdict.get(id,'0.0')
        res.append([fid,iid,f_missing,nmend,pedsex,snpsex,status,fest,fhet])
        try:
            ff_missing = float(f_missing)
        except:
            ff_missing = 0.0
        try:
            inmend = int(nmend)
        except:
            inmend = 0
        try:
            ffest = float(fest)
        except:
            fest = 0.0
        try:
            ffhet = float(fhet)
        except:
            ffhet = 0.0
        fres.append([fid,iid,ff_missing,inmend,pedsex,snpsex,status,ffest,ffhet])
    ntokeep = max(20,len(res)/keepfrac)
    for i,col in enumerate(Tsorts):
        fres.sort(key=operator.itemgetter(col))
        if Treverse[i]:
            fres.reverse()
        repname = Tnames[i]
        Tops[repname] = fres[0:ntokeep]
        Tops[repname] = [map(str,x) for x in Tops[repname]]
        Tops[repname].insert(0,rhead)
    res.sort()
    res.insert(0,rhead)
    logf.write('### writing %s report with %s' % ( outfile,res[0]))  
    f = open(outfile,'w')
    f.write('\n'.join(['\t'.join(x) for x in res]))
    f.write('\n')
    f.close()
    return res,Tops

def markerRep(froot='cleantest',outfname="mrep",newfpath='.',logf=None,maplist=None ):
    """by marker (hwe = .hwe, missingness=.lmiss, freq = .frq)
    keep a list of marker order but keep all stats in dicts
    write out a fake xls file for R or SAS etc
    kinda clunky, but..
    TODO: ensure stable if any file not found?
    """
    mapdict = {}
    if maplist <> None:
       rslist = [x[1] for x in maplist]
       offset = [(x[0],x[3]) for x in maplist]
       mapdict = dict(zip(rslist,offset))
    hwefile = '%s.hwe' % froot
    lmissfile = '%s.lmiss' % froot
    freqfile = '%s.frq' % froot
    lmendfile = '%s.lmendel' % froot
    outfile = os.path.join(newfpath,outfname)
    markerlist = []
    chromlist = []
    hwedict = {}
    lmissdict = {}
    freqdict = {}
    lmenddict = {}
    Tops = {}
    Tnames = ['Ranked Marker MAF', 'Ranked Marker Missing Genotype', 'Ranked Marker HWE', 'Ranked Marker Mendel']
    Tsorts = [3,6,10,11]
    Treverse = [False,True,True,True] # so first values are worse(r)
    #res.append([rs,chrom,offset,maf,a1,a2,f_missing,hwe_all[0],hwe_all[1],hwe_unaff[0],hwe_unaff[1],nmend])
    #rhead = ['snp','chrom','maf','a1','a2','missfrac','p_hwe_all','logp_hwe_all','p_hwe_unaff','logp_hwe_unaff','N_Mendel']
    # -------------------hwe--------------------------
    #    hwe has SNP TEST  GENO   O(HET)   E(HET) P_HWD
    # we want all hwe where P_HWD <> NA
    # ah changed in 1.04 to
    ##  CHR         SNP     TEST   A1   A2                 GENO   O(HET)   E(HET)            P 
    ##   1   rs6671164      ALL    2    3           34/276/613    0.299   0.3032       0.6644
    ##   1   rs6671164      AFF    2    3                0/0/0      nan      nan           NA
    ##   1   rs6671164    UNAFF    2    3           34/276/613    0.299   0.3032       0.6644
    ##   1   rs4448553      ALL    2    3            8/176/748   0.1888   0.1848       0.5975
    ##   1   rs4448553      AFF    2    3                0/0/0      nan      nan           NA
    ##   1   rs4448553    UNAFF    2    3            8/176/748   0.1888   0.1848       0.5975
    ##   1   rs1990150      ALL    1    3           54/303/569   0.3272   0.3453       0.1067
    ##   1   rs1990150      AFF    1    3                0/0/0      nan      nan           NA
    ##   1   rs1990150    UNAFF    1    3           54/303/569   0.3272   0.3453       0.1067
    logf.write('## marker reports starting at %s\n' % timenow())
    try:
        f = open(hwefile,'r')
    except:
        f = None
        logf.write('## error - no hwefile %s found\n' % hwefile)
    if f:
        for i,line in enumerate(f):
            ll = line.strip().split()
            if i == 0: # head
                head = [x.upper() for x in ll] # expect above                
                try:
                    testpos = head.index('TEST')
                except:
                    testpos = 2 # patch for 1.04 which has b0rken headers - otherwise use head.index('TEST')
                try:
                    ppos = head.index('P')
                except:
                    ppos = 8 # patch - for head.index('P')
                snppos = head.index('SNP')
                chrpos = head.index('CHR')
                logf.write('hwe header testpos=%d,ppos=%d,snppos=%d\n' % (testpos,ppos,snppos))
            elif len(ll) >= ppos: # full line
                ps = ll[ppos].upper()
                rs = ll[snppos]
                chrom = ll[chrpos]
                test = ll[testpos]
                if not hwedict.get(rs,None):
                    hwedict[rs] = {}
                    markerlist.append(rs)
                chromlist.append(chrom) # one place to find it?
                lpvals = 0
                if ps.upper() <> 'NA' and ps.upper() <> 'NAN': # worth keeping
                    lpvals = '0'
                    if ps <> '1':
                        try:
                            pval = float(ps)
                            lpvals = '%f' % -math.log10(pval)
                        except:
                            pass
                    hwedict[rs][test] = (ps,lpvals)
            else:
                logf.write('short line #%d = %s\n' % (i,ll))
        f.close()
    # ------------------missing--------------------
    """lmiss has  
    CHR          SNP   N_MISS   N_GENO   F_MISS
   1   rs12354060        0       73        0
   1    rs4345758        1       73   0.0137
   1    rs2691310       73       73        1
   1    rs2531266       73       73        1
    # we want F_MISS"""
    try:
        f = open(lmissfile,'r')
    except:
        f = None
    if f:
        for i,line in enumerate(f):
            ll = line.strip().split()
            if i == 0:
                head = [x.upper() for x in ll] # expect above                
                fracpos = head.index('F_MISS')
                npos = head.index('N_MISS')
                snppos = head.index('SNP')
            elif len(ll) >= fracpos: # full line
                rs = ll[snppos]
                fracval = ll[fracpos]                
                lmissdict[rs] = fracval # for now, just that?
            else:
                lmissdict[rs] = 'NA'
        f.close()
    # ------------------freq-------------------
    # frq has CHR          SNP   A1   A2          MAF       NM
    # we want maf
    try:
        f = open(freqfile,'r')
    except:
        f = None
    if f:
        for i,line in enumerate(f):
            ll = line.strip().split()
            if i == 0:
                head = [x.upper() for x in ll] # expect above                
                mafpos = head.index('MAF')
                a1pos = head.index('A1')
                a2pos = head.index('A2')
                snppos = head.index('SNP')
            elif len(ll) >= mafpos: # full line
                rs = ll[snppos]
                a1 = ll[a1pos]
                a2 = ll[a2pos]
                maf = ll[mafpos]
                freqdict[rs] = (maf,a1,a2)
        f.close()
    # ------------------mend-------------------
    # lmend has CHR SNP   N
    # we want N
    gotmend = True
    try:
        f = open(lmendfile,'r')
    except:
        gotmend = False
        for rs in markerlist:
            lmenddict[rs] = '0'
    if gotmend:
        for i,line in enumerate(f):
            ll = line.strip().split()
            if i == 0:
                head = [x.upper() for x in ll] # expect above                
                npos = head.index('N')
                snppos = head.index('SNP')
            elif len(ll) >= npos: # full line
                rs = ll[snppos]
                nmend = ll[npos]
                lmenddict[rs] = nmend
        f.close()
    else:
        logf.write('No %s file - assuming not family data\n' % lmendfile)
    # now assemble result list
    rhead = ['snp','chromosome','offset','maf','a1','a2','missfrac','p_hwe_all','logp_hwe_all','p_hwe_unaff','logp_hwe_unaff','N_Mendel']
    res = []
    fres = []
    for rs in markerlist: # for each snp in found order
        f_missing = lmissdict.get(rs,'NA')
        maf,a1,a2 = freqdict.get(rs,('NA','NA','NA'))
        hwe_all = hwedict[rs].get('ALL',('NA','NA')) # hope this doesn't change...
        hwe_unaff = hwedict[rs].get('UNAFF',('NA','NA'))
        nmend = lmenddict.get(rs,'NA')
        (chrom,offset)=mapdict.get(rs,('?','0'))
        res.append([rs,chrom,offset,maf,a1,a2,f_missing,hwe_all[0],hwe_all[1],hwe_unaff[0],hwe_unaff[1],nmend])
    ntokeep = max(10,len(res)/keepfrac)

    def msortk(item=None):
        """
        deal with non numeric sorting
        """
        try:
           return float(item)
        except:
           return item

    for i,col in enumerate(Tsorts):
        res.sort(key=msortk(lambda x:x[col]))
        if Treverse[i]:
            res.reverse()
        repname = Tnames[i]
        Tops[repname] = res[0:ntokeep]
        Tops[repname].insert(0,rhead)
    res.sort(key=lambda x: '%s_%10d' % (x[1].ljust(4,'0'),int(x[2]))) # in chrom offset order
    res.insert(0,rhead)
    f = open(outfile,'w')
    f.write('\n'.join(['\t'.join(x) for x in res]))
    f.close()
    return res,Tops



  
def getfSize(fpath,outpath):
    """
    format a nice file size string
    """
    size = ''
    fp = os.path.join(outpath,fpath)
    if os.path.isfile(fp):
        n = float(os.path.getsize(fp))
        if n > 2**20:
            size = ' (%1.1f MB)' % (n/2**20)
        elif n > 2**10:
            size = ' (%1.1f KB)' % (n/2**10)
        elif n > 0:
            size = ' (%d B)' % (int(n))
    return size


if __name__ == "__main__":
    u = """ called in xml as
     <command interpreter="python">
        rgQC.py -i '$input_file.extra_files_path/$input_file.metadata.base_name' -o "$out_prefix" 
        -s '$html_file' -p '$html_file.files_path' -l '${GALAXY_DATA_INDEX_DIR}/rg/bin/plink' 
        -r '${GALAXY_DATA_INDEX_DIR}/rg/bin/R' 
    </command>

        Prepare a qc report - eg:
    print "%s %s -i birdlped -o birdlped -l qc.log -s htmlf -m marker.xls -s sub.xls -p ./" % (sys.executable,prog)

    """
    progname = os.path.basename(sys.argv[0])
    if len(sys.argv) < 9:
        print '%s requires 6 parameters - got %d = %s' % (progname,len(sys.argv),sys.argv)
        sys.exit(1)
    parser = OptionParser(usage=u, version="%prog 0.01")
    a = parser.add_option
    a("-i","--infile",dest="infile")
    a("-o","--oprefix",dest="opref")
    a("-l","--plinkexe",dest="plinke", default=plinke)
    a("-r","--rexe",dest="rexe", default=rexe)
    a("-s","--snps",dest="htmlf")
    #a("-m","--markerRaw",dest="markf")
    #a("-r","--rawsubject",dest="subjf")
    a("-p","--patho",dest="newfpath")
    (options,args) = parser.parse_args()
    basename = os.path.split(options.infile)[-1] # just want the file prefix to find the .xls files below
    killme = string.punctuation + string.whitespace
    trantab = string.maketrans(killme,'_'*len(killme))
    opref = options.opref.translate(trantab)
    alogh,alog = tempfile.mkstemp(suffix='.txt')
    plogh,plog = tempfile.mkstemp(suffix='.txt')
    alogf = open(alog,'w')
    plogf = open(plog,'w')
    ahtmlf = options.htmlf
    amarkf = 'MarkerDetails_%s.xls' % opref
    asubjf = 'SubjectDetails_%s.xls' % opref
    newfpath = options.newfpath
    newfpath = os.path.realpath(newfpath)
    try:
       os.makedirs(newfpath)
    except:
       pass
    ofn = basename
    bfn = options.infile
    try:
       mapf = '%s.bim' % bfn
       maplist = file(mapf,'r').readlines()
       maplist = [x.split() for x in maplist]
    except:
       maplist = None
       alogf.write('## error - cannot open %s to read map - no offsets will be available for output files')
    #rerla@beast galaxy]$ head test-data/tinywga.bim
    #22      rs2283802       0       21784722        4       2
    #22      rs2267000       0       21785366        4       2
    rgbin = os.path.split(rexe)[0] # get our rg bin path
    #plinktasks = [' --freq',' --missing',' --mendel',' --hardy',' --check-sex'] # plink v1 fixes that bug!
    # if we could, do all at once? Nope. Probably never.
    plinktasks = [['--freq',],['--hwe 0.0', '--missing','--hardy'],
    ['--mendel',],['--check-sex',]]
    vclbase = [options.plinke,'--noweb','--out',basename,'--bfile',bfn,'--mind','1.0','--geno','1.0','--maf','0.0']
    runPlink(logf=plogf,plinktasks=plinktasks,cd=newfpath, vclbase=vclbase)
    plinktasks = [['--bfile',bfn,'--indep-pairwise 40 20 0.5','--out %s' % basename],
    ['--bfile',bfn,'--extract %s.prune.in --make-bed --out ldp_%s' % (basename, basename)],
                  ['--bfile ldp_%s --het --out %s' % (basename,basename)]]
    # subset of ld independent markers for eigenstrat and other requirements
    vclbase = [options.plinke,'--noweb']
    plogout = pruneLD(plinktasks=plinktasks,cd=newfpath,vclbase = vclbase)
    plogf.write('\n'.join(plogout))
    plogf.write('\n')
    repout = os.path.join(newfpath,basename)
    subjects,subjectTops = subjectRep(froot=repout,outfname=asubjf,newfpath=newfpath,
                logf=alogf) # writes the subject_froot.xls file
    markers,markerTops = markerRep(froot=repout,outfname=amarkf,newfpath=newfpath,
                logf=alogf,maplist=maplist) # marker_froot.xls
    nbreaks = 100
    s = '## starting plotpage, newfpath=%s,m=%s,s=%s/n' % (newfpath,markers[:2],subjects[:2])
    alogf.write(s)
    print s
    plotpage,cruft = makePlots(markers=markers,subjects=subjects,newfpath=newfpath,
                         basename=basename,nbreaks=nbreaks,height=10,width=8,rgbin=rgbin)
    #plotpage = RmakePlots(markers=markers,subjects=subjects,newfpath=newfpath,basename=basename,nbreaks=nbreaks,rexe=rexe)

    # [titles[n],plotnames[n],outfnames[n] ]
    html = []
    html.append('<table cellpadding="5" border="0">')
    size = getfSize(amarkf,newfpath)
    html.append('<tr><td colspan="3"><a href="%s" type="application/vnd.ms-excel">%s</a>%s tab delimited</td></tr>' % \
                (amarkf,'Click here to download the Marker QC Detail report file',size))
    size = getfSize(asubjf,newfpath)
    html.append('<tr><td colspan="3"><a href="%s" type="application/vnd.ms-excel">%s</a>%s tab delimited</td></tr>' % \
                (asubjf,'Click here to download the Subject QC Detail report file',size))
    for (title,url,ofname) in plotpage:
        ttitle = 'Ranked %s' % title
        dat = subjectTops.get(ttitle,None)
        if not dat:
            dat = markerTops.get(ttitle,None)
        imghref = '%s.jpg' % os.path.splitext(url)[0] # removes .pdf
        thumbnail = os.path.join(newfpath,imghref)
        if not os.path.exists(thumbnail): # for multipage pdfs, mogrify makes multiple jpgs - fugly hack
            imghref = '%s-0.jpg' % os.path.splitext(url)[0] # try the first jpg
            thumbnail = os.path.join(newfpath,imghref)
        if not os.path.exists(thumbnail):
            html.append('<tr><td colspan="3"><a href="%s">%s</a></td></tr>' % (url,title))
        else:
            html.append('<tr><td><a href="%s"><img src="%s" alt="%s" hspace="10" align="middle">' \
                    % (url,imghref,title))
            if dat: # one or the other - write as an extra file and make a link here
                t = '%s.xls' % (ttitle.replace(' ','_'))
                fname = os.path.join(newfpath,t)
                f = file(fname,'w')
                f.write('\n'.join(['\t'.join(x) for x in dat])) # the report
                size = getfSize(t,newfpath)
                html.append('</a></td><td>%s</td><td><a href="%s">Worst data</a>%s</td></tr>' % (title,t,size))
            else:
                html.append('</a></td><td>%s</td><td>&nbsp;</td></tr>' % (title))
    html.append('</table><hr><h3>All output files from the QC run are available below</h3>')
    html.append('<table cellpadding="5" border="0">\n')
    flist = os.listdir(newfpath) # we want to catch 'em all
    flist.sort()
    for f in flist:
        fname = os.path.split(f)[-1]
        size = getfSize(fname,newfpath)
        html.append('<tr><td><a href="%s">%s</a>%s</td></tr>' % (fname,fname,size))
    html.append('</table>')
    alogf.close()
    plogf.close()
    llog = open(alog,'r').readlines()
    plogfile = open(plog,'r').readlines()
    os.unlink(alog)
    os.unlink(plog)
    llog += plogfile # add lines from pruneld log
    lf = file(ahtmlf,'w') # galaxy will show this as the default view
    lf.write(galhtmlprefix % progname)
    s = '\n<div>Output from Rgenetics QC report tool run at %s<br>\n' % (timenow())
    lf.write('<h4>%s</h4>\n' % s)
    lf.write('</div><div><h4>(Click any preview image to download a full sized PDF version)</h4><br><ol>\n')
    lf.write('\n'.join(html))
    lf.write('<h4>QC run log contents</h4>')
    lf.write('<pre>%s</pre>' % (''.join(llog))) # plink logs
    if len(cruft) > 0:
        lf.write('<h2>Blather from pdfnup follows:</h2><pre>%s</pre>' % (''.join(cruft))) # pdfnup
    lf.write('%s\n<hr>\n' % galhtmlpostfix)
    lf.close()