@transform( loadPolyphen, suffix(".load"), ".genestats")
def analysePolyphen( infile, outfile ):
    
    table = P.toTable( infile )
    tablename_map = "polyphen_map"

    dbhandle = sqlite3.connect( PARAMS["database"] )
    cc = dbhandle.cursor()

    statement = '''
        SELECT i.gene_id,
               COUNT(DISTINCT map.locus_id) as nsnps, 
               COUNT(DISTINCT case t.prediction when 'possiblydamaging' then map.locus_id when 'probablydamaging' then map.locus_id else NULL end) AS ndeleterious,
               MAX(s.length)
               FROM %(table)s as t, 
                    %(tablename_map)s as map, 
                    protein_stats as s,
                    transcript_info as i 
        WHERE map.snp_id = t.snp_id AND 
              i.transcript_id = map.transcript_id AND
              s.protein_id = map.protein_id
        GROUP BY i.gene_id
     ''' % locals()

    data = cc.execute(statement).fetchall()

    statement = '''SELECT DISTINCT i.gene_id, MAX(s.length) 
                   FROM transcript_info AS i, protein_stats AS s 
                   WHERE s.protein_id = i.protein_id 
                   GROUP BY i.gene_id'''
    gene_ids = cc.execute(statement).fetchall()

    total_nsnps = sum( [ x[1] for x in data ] )
    total_ndel = sum( [ x[2] for x in data ] )
    total_length = sum( [ x[1] for x in gene_ids ] )
    del_p = float(total_ndel) / total_nsnps
    len_p = float(total_nsnps) / total_length
    com_p = float(total_ndel) / total_length

    E.info( "del: background probability: %i/%i = %f" % (total_ndel, total_nsnps, del_p ) )
    E.info( "len: background probability: %i/%i = %f" % (total_nsnps, total_length, len_p ) )
    E.info( "com: background probability: %i/%i = %f" % (total_ndel, total_length, com_p ) )

    outf = open( outfile, "w" )
    outf.write( "\t".join( ("gene_id", "code", 
                            "length", "nsnps", "ndel", 
                            "del_p", "del_pvalue", "del_qvalue", 
                            "len_p", "len_pvalue", "len_qvalue",
                            "com_p", "com_pvalue", "com_qvalue", ) ) + "\n" )

    del_pvalues, len_pvalues, com_pvalues = [], [], []
    for gene_id, nsnps, ndel, length in data:

        # use -1, because I need P( x >= X)
        # sf = 1 - cdf and cdf = P( x <= X ), thus sf = 1 - P( x <= X ) = P (x > X ).
        del_pvalues.append( scipy.stats.binom.sf( ndel - 1, nsnps, del_p ) ) 
        len_pvalues.append( scipy.stats.binom.sf( nsnps - 1, int(round(length)), len_p ) )
        com_pvalues.append( scipy.stats.binom.sf( ndel - 1, int(round(length)), com_p ) )

    del_q = Stats.doFDR( del_pvalues )
    len_q = Stats.doFDR( len_pvalues )
    com_q = Stats.doFDR( com_pvalues )

    fdr = PARAMS["polyphen_fdr"]

    found = set()

    for a, del_pvalue, del_qvalue, len_pvalue, len_qvalue, com_pvalue, com_qvalue in \
            zip(data, 
                del_pvalues, del_q.mQValues, 
                len_pvalues, len_q.mQValues,
                com_pvalues, com_q.mQValues,
                ):
        gene_id, nsnps, ndel, length = a
        found.add(gene_id)

        del_p = float(ndel) / nsnps
        len_p = float(nsnps) / length

        code = "".join( [ str(int(x < fdr)) for x in (del_qvalue, len_qvalue, com_qvalue) ] )

        outf.write( "\t".join( (gene_id,
                                code,
                                "%i" % int(round(length)),
                                "%i" % int(nsnps),
                                "%i" % int(ndel),
                                "%6.4f" % del_p,
                                "%6.4g" % del_pvalue,
                                "%6.4g" % del_qvalue,
                                "%6.4f" % len_p,
                                "%6.4g" % len_pvalue,
                                "%6.4g" % len_qvalue,
                                "%6.4f" % com_p,
                                "%6.4g" % com_pvalue,
                                "%6.4g" % com_qvalue,
                                ) ) + "\n" )

    # add missing genes:
    code = "---"
    for gene_id, length in gene_ids:
        if gene_id in found: continue
        outf.write( "\t".join( (gene_id,
                                code,
                                "%i" % int(round(length)),
                                "%i" % 0,
                                "%i" % 0,
                                "%6.4f" % 0,
                                "%6.4g" % 1,
                                "%6.4g" % 1,
                                "%6.4f" % 0,
                                "%6.4g" % 1,
                                "%6.4g" % 1,
                                "%6.4f" % 0,
                                "%6.4g" % 1,
                                "%6.4g" % 1,
                                ) ) + "\n" )

    outf.close()
    
###################################################################
###################################################################
###################################################################
Returns the user ID
dataset/ETHPy150Open CGATOxford/cgat/obsolete/pipeline_snps.py/analysePolyphen--crossover with--dataset/ETHPy150Open jgilhutton/pyxiewps_WPShack-Python/Version-1.0/pyxiewps-EN.py/Config.get_uid