<gh_stars>0

'''
calc_rpvg_gibbs_stats.py 
Calculates rpvg Gibbs samples statistics 
(mean, median and credibility intervals). 
'''

import sys
import os
import subprocess
import gzip
import numpy

from utils import *


printScriptHeader()

if len(sys.argv) != 3:

  print("Usage: python calc_rpvg_gibbs_stats.py <rpvg_gibbs_gz_name> <output_file_name>\n")
  sys.exit(1)


read_count_sums = []

gibbs_file = gzip.open(sys.argv[1], "rb")

for line in gibbs_file:

    line_split = line.decode().strip().split("\t")

    if line_split[0] == "Name":

        assert(len(line_split) > 2)
        assert(len(read_count_sums) == 0)      

        for i in range(0, len(line_split) - 2):

            read_count_sums.append(0)

        continue

    assert(len(line_split) == len(read_count_sums) + 2)

    for i in range(0, len(read_count_sums)):

        read_count_sums[i] = read_count_sums[i] + float(line_split[i + 2])

gibbs_file.close()

print(min(read_count_sums))
print(numpy.mean(read_count_sums))
print(numpy.median(read_count_sums))
print(max(read_count_sums))

out_file = open(sys.argv[2], "w")
out_file.write("Name\tClusterID\tMean\tMedian\tMin\tMax\tCILower95\tCIUpper95\tCILower90\tCIUpper90\n")

gibbs_file = gzip.open(sys.argv[1], "rb")

for line in gibbs_file:

    line_split = line.decode().strip().split("\t")

    if line_split[0] == "Name":

        continue

    assert(len(line_split) == len(read_count_sums) + 2)

    samples = []

    for i in range(0, len(read_count_sums)):

        samples.append(float(line_split[i + 2]) / read_count_sums[i] * 1000000)

    samples.sort()

    cil95 = samples[int(round(0.025  * len(samples)))]
    ciu95 = samples[int(round(0.975 * len(samples)))]
    cil90 = samples[int(round(0.05 * len(samples)))]
    ciu90 = samples[int(round(0.95 * len(samples)))]

    out_line = [line_split[0], line_split[1], numpy.mean(samples), numpy.median(samples), min(samples), max(samples), cil95, ciu95, cil90, ciu90]
    out_file.write("\t".join([str(x) for x in out_line]) + "\n")

gibbs_file.close()
out_file.close()

print("\nDone")
