<commit_before>#!/usr/bin/env python2

'''
This script takes two arguments:
- an input filename (the output of epacts with a single phenotype)
- an output filename (please end it with `.json`)

It creates a json file which can be used to render a QQ plot.
'''

from __future__ import print_function, division, absolute_import

import os.path
import sys
import gzip
import re
import json
import math
import collections
import scipy.stats
from itertools import islice
from extreme_collection import ExtremeCollection

NEGLOG10_PVAL_BIN_SIZE = 0.05 # Use 0.05, 0.1, 0.15, etc
NEGLOG10_PVAL_BIN_DIGITS = 2 # Then keep 2 digits after the decimal

def round_sig(x, digits):
    return 0 if x==0 else round(x, digits-1-int(math.floor(math.log10(abs(x)))))
assert round_sig(0.00123, 2) == 0.0012
assert round_sig(1.59e-10, 2) == 1.6e-10

def median(x):
    """returns median assuming values are sorted"""
    length = len(x)
    middle = (length - 1) // 2

    if (length %2):
        return x[middle]
    else:
        return (x[middle] + x[middle+1]) / 2

def approx_equal(a, b, tolerance=1e-4):
    return abs(a-b) <= max(abs(a), abs(b)) * tolerance
assert approx_equal(42, 42.0000001)
assert not approx_equal(42, 42.01)


def gc_value(median_pval):
    # This should be equivalent to this R: `qchisq(p, df=1, lower.tail=F) / qchisq(.5, df=1, lower.tail=F)`
    return scipy.stats.chi2.ppf(1 - median_pval, 1) / scipy.stats.chi2.ppf(0.5, 1)
assert approx_equal(gc_value(0.49), 1.047457) # I computed these using that R code.
assert approx_equal(gc_value(0.5), 1)
assert approx_equal(gc_value(0.50001), 0.9999533)
assert approx_equal(gc_value(0.6123), 0.5645607)

def rounded(x):
    return round(x // NEGLOG10_PVAL_BIN_SIZE * NEGLOG10_PVAL_BIN_SIZE, NEGLOG10_PVAL_BIN_DIGITS)

def make_qq(variants, max_unbinned, num_bins):
    # smallest first
    sorted_variants = sorted((x for x in variants if x is not None), key=lambda x: x.pval)
    count = len(sorted_variants)

    max_exp_neglog10_pval = -math.log10(0.5 / count) #expected
    max_obs_neglog10_pval = -math.log10(sorted_variants[0].pval )

    current_bin = -1
    variants_in_bin = []
    binned_variants = []
    for i in xrange(max_unbinned+1, count):
        exp_neglog10_pval = -math.log10( (i+0.5) / count ) 
        exp_bin = int(exp_neglog10_pval / max_exp_neglog10_pval * num_bins)
        obs_neglog10_pval = -math.log10( sorted_variants[i].pval )
        if exp_bin != current_bin:
            if len(variants_in_bin)>0:
                binned_variants.append((
                    round_sig(current_bin / num_bins * max_exp_neglog10_pval, 3),
                    round_sig(sum(variants_in_bin) / float(len(variants_in_bin)), 3) 
                ))
                variants_in_bin = []
            current_bin = exp_bin
        variants_in_bin.append(obs_neglog10_pval)
    if len(variants_in_bin)>0:
        binned_variants.append((
            round_sig(current_bin / num_bins * max_exp_neglog10_pval, 3),
            round_sig(sum(variants_in_bin) / float(len(variants_in_bin)), 3) 
        ))

    unbinned_variants = []
    for i in xrange(0, min(max_unbinned, count)):
        variant = sorted_variants[i]
        obs_neglog10_pval = -math.log10( variant.pval )
        exp_neglog10_pval = -math.log10( (i+0.5) / count )
        rec = {
            'chrom': variant.chrom,
            'pos': variant.pos,
            'pval': variant.pval
        }
        rec.update(variant.other)
        unbinned_variants.append( (round_sig(exp_neglog10_pval,3), round_sig(obs_neglog10_pval,3), rec) )

    median_pval = sorted_variants[count//2].pval
    gc = {
        "50": round_sig(gc_value(median_pval),5) 
    }
        
    return {"variant_bins": binned_variants, "unbinned_variants": unbinned_variants, "gc": gc, "count": count} 


def process_file_unbinned(results, max_unbinned, num_bins):
    plot = {"description": "Overall Results", "layers": [make_qq(results, max_unbinned, num_bins)]}<commit_msg>Add conf int to QQ plots<commit_after>#!/usr/bin/env python2

'''
This script takes two arguments:
- an input filename (the output of epacts with a single phenotype)
- an output filename (please end it with `.json`)

It creates a json file which can be used to render a QQ plot.
'''

from __future__ import print_function, division, absolute_import

import os.path
import sys
import gzip
import re
import json
import math
import collections
import scipy.stats
from itertools import islice
from extreme_collection import ExtremeCollection

NEGLOG10_PVAL_BIN_SIZE = 0.05 # Use 0.05, 0.1, 0.15, etc
NEGLOG10_PVAL_BIN_DIGITS = 2 # Then keep 2 digits after the decimal

def round_sig(x, digits):
    return 0 if x==0 else round(x, digits-1-int(math.floor(math.log10(abs(x)))))
assert round_sig(0.00123, 2) == 0.0012
assert round_sig(1.59e-10, 2) == 1.6e-10

def median(x):
    """returns median assuming values are sorted"""
    length = len(x)
    middle = (length - 1) // 2

    if (length %2):
        return x[middle]
    else:
        return (x[middle] + x[middle+1]) / 2

def approx_equal(a, b, tolerance=1e-4):
    return abs(a-b) <= max(abs(a), abs(b)) * tolerance
assert approx_equal(42, 42.0000001)
assert not approx_equal(42, 42.01)


def gc_value(median_pval):
    # This should be equivalent to this R: `qchisq(p, df=1, lower.tail=F) / qchisq(.5, df=1, lower.tail=F)`
    return scipy.stats.chi2.ppf(1 - median_pval, 1) / scipy.stats.chi2.ppf(0.5, 1)
assert approx_equal(gc_value(0.49), 1.047457) # I computed these using that R code.
assert approx_equal(gc_value(0.5), 1)
assert approx_equal(gc_value(0.50001), 0.9999533)
assert approx_equal(gc_value(0.6123), 0.5645607)

def rounded(x):
    return round(x // NEGLOG10_PVAL_BIN_SIZE * NEGLOG10_PVAL_BIN_SIZE, NEGLOG10_PVAL_BIN_DIGITS)

def get_conf_int(nvar):
    slices = []
    for x in range(0, int(math.ceil(math.log(nvar,2)))):
        slices.append(2**x)
    slices.append(nvar-1);
    slices.reverse()

    points = []
    for slice in slices:
        rv = scipy.stats.beta(slice, nvar-slice)
        points.append((
            round(-math.log10((slice-0.5)/nvar),2),
            round(-math.log10(rv.ppf(0.05/2)),2), 
            round(-math.log10(rv.ppf(1-(0.05/2))),2)
        ))
    return points

def make_qq(variants, max_unbinned, num_bins):
    # smallest first
    sorted_variants = sorted((x for x in variants if x is not None), key=lambda x: x.pval)
    count = len(sorted_variants)

    max_exp_neglog10_pval = -math.log10(0.5 / count) #expected
    max_obs_neglog10_pval = -math.log10(sorted_variants[0].pval )

    current_bin = -1
    variants_in_bin = []
    binned_variants = []
    for i in xrange(max_unbinned+1, count):
        exp_neglog10_pval = -math.log10( (i+0.5) / count ) 
        exp_bin = int(exp_neglog10_pval / max_exp_neglog10_pval * num_bins)
        obs_neglog10_pval = -math.log10( sorted_variants[i].pval )
        if exp_bin != current_bin:
            if len(variants_in_bin)>0:
                binned_variants.append((
                    round_sig(current_bin / num_bins * max_exp_neglog10_pval, 3),
                    round_sig(sum(variants_in_bin) / float(len(variants_in_bin)), 3) 
                ))
                variants_in_bin = []
            current_bin = exp_bin
        variants_in_bin.append(obs_neglog10_pval)
    if len(variants_in_bin)>0:
        binned_variants.append((
            round_sig(current_bin / num_bins * max_exp_neglog10_pval, 3),
            round_sig(sum(variants_in_bin) / float(len(variants_in_bin)), 3) 
        ))

    unbinned_variants = []
    for i in xrange(0, min(max_unbinned, count)):
        variant = sorted_variants[i]
        obs_neglog10_pval = -math.log10( variant.pval )
        exp_neglog10_pval = -math.log10( (i+0.5) / count )
        rec = {
            'chrom': variant.chrom,
            'pos': variant.pos,
            'pval': variant.pval
        }
        rec.update(variant.other)
        unbinned_variants.append( (round_sig(exp_neglog10_pval,3), round_sig(obs_neglog10_pval,3), rec) )

    median_pval = sorted_variants[count//2].pval
    gc = {
        "50": round_sig(gc_value(median_pval),5) 
    }
    conf_int = get_conf_int(count)
        
    return {"variant_bins": binned_variants, "unbinned_variants": unbinned_variants, 
        "gc": gc, "count": count, "conf_int": conf_int} 


def process_file_unbinned(results, max_unbinned, num_bins):
    plot = {"description": "Overall Results", "layers": [make_qq(results, max_unbinned, num_bins)]}
{'fraction_non_alphanumeric': 0.12080536912751678, 'fraction_numerical': 0.0463911456493583, 'mean_word_length': 3.1393762183235867, 'pattern_counts': {'":': 15, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5917755', 'n_tokens_mistral': 3398, 'n_tokens_neox': 3140, 'n_words': 673}