def get_calpha(args):
    db = args.db
    if not (args.controls and args.cases):
        case, control = _get_case_and_control_samples(args)
    else:
        case = args.cases
        control = args.controls
    assert (case and control), ("Phenotypes not found in the database and "
                                "--cases and --controls are not set.")

    samples = control + case
    # p_0 = the fraction of samples that are cases (used for weighting)
    p_0 = float(len(case)) / float(len(samples))

    if args.nonsynonymous:
        ns = _nonsynonymous_variants(args)
    else:
        ns = _medium_or_high_impact_variants(args)

    variants_in_gene, variants = _calculate_counts(ns, samples)
    header = ["gene", "T", "c", "Z", "p_value"]
    print "\t".join(header)

    if args.permutations > 0:
        perms = permute_cases(samples, args.permutations, case)

    for gene in variants_in_gene:
        vig = variants_in_gene[gene]

        # m = the number of variants observed for this gene
        m = len(vig.keys())

        # m_n is the number of variants with n copies (i.e., samples with the variant)
        #m_n = Counter([len(x) for x in vig.values()])

        # n_i is a list reflecting the total number of samples
        # having each variant
        n_i = [len(x) for x in vig.values()]

        # y_i is a list reflecting the total number of __cases__
        # having each variant
        y_i = [len(filter(lambda a: a in case, x)) for x in vig.values()]

        # "The C-alpha test statistic T contrasts the variance of each observed
        # count with the expected variance, assuming the binomial distribution."
        # In other words, given that we have n total samples and p_0 * n of them
        # are cases, we _expect_ the variant copies to be distributed among the
        # samples following a binomal distribution.  The T statistic contrasts
        # the observed count distributions with the expected:
        #
        # T = SUM{i=(1,m)} [(y_i - n_i*p_0)^2 - n_i*p_0(1 - p_0)]
        #
        T = _calculate_T(m, p_0, n_i, y_i)

        # Calculate the variance of T in order to normalize it
        c = _calculate_c(n_i, p_0)

        # The final test statistic, Z, id just the original test statistic divided
        # by its standard deviation. "We reject the null when Z is larger than expected
        # using a one-tailed standard normal distribution for reference.
        if c == 0:
            Z = np.NaN
            p_value = np.NaN
            print "\t".join([gene, str(T), str(c), str(Z), str(p_value)])
            continue
        else:
            Z = T / math.sqrt(c)

        if args.permutations == 0:
            # sf is the survival function ... same as 1 - CDF.
            p_value = norm.sf(Z)
        else:
            # this permutes the cases without replacement, important for
            # calculating an exact p-value
            T_scores = []
            for perm_case in perms:
                y_i = [len(filter(lambda a: a in perm_case, x)) for x in vig.values()]
                T_permuted = _calculate_T(m, p_0, n_i, y_i)
                T_scores.append(T_permuted)
            if args.save_tscores:
                with open("permutated_t_scores.txt", "a") as out_handle:
                    out_handle.write("\t".join([gene] + map(str, T_scores)) + "\n")
            false_hits = sum([x >= T for x in T_scores])
            # the + 1 to make it an unbiased estimator
            # Permutation P-values Should Never Be Zero: Calculating Exact
            # P-values When Permutations Are Randomly Drawn
            # http://www.degruyter.com/view/j/sagmb.2010.9.1/sagmb.2010.9.1.1585/sagmb.2010.9.1.1585.xml
            p_value = (float(false_hits) + 1) / (float(args.permutations + 1))

        print "\t".join([gene, str(T), str(c), str(Z), str(p_value)])
Calculate the C-alpha statistic for each gene based on the observed
dataset/ETHPy150Open arq5x/gemini/gemini/tool_burden_tests.py/get_calpha