galaxy/tools/solid_tools/solid_qual_stats.py

#!/usr/bin/env python
# Guruprasad Ananda
from __future__ import print_function

import sys
import tempfile
import zipfile

QUAL_UPPER_BOUND = 41
QUAL_LOWER_BOUND = 1


def stop_err(msg):
    sys.stderr.write("%s\n" % msg)
    sys.exit()


def unzip(filename):
    zip_file = zipfile.ZipFile(filename, "r")
    tmpfilename = tempfile.NamedTemporaryFile().name
    for name in zip_file.namelist():
        open(tmpfilename, "a").write(zip_file.read(name))
    zip_file.close()
    return tmpfilename


def __main__():
    infile_score_name = sys.argv[1].strip()
    with open(sys.argv[2].strip(), "w") as fout:
        if zipfile.is_zipfile(infile_score_name):
            infile_name = unzip(infile_score_name)
        else:
            infile_name = infile_score_name

        readlen = None
        invalid_lines = 0
        j = 0
        with open(infile_name) as fin:
            for line in fin:
                line = line.strip()
                if not (line) or line.startswith("#") or line.startswith(">"):
                    continue
                elems = line.split()
                try:
                    for item in elems:
                        int(item)
                    if not readlen:
                        readlen = len(elems)
                    if len(elems) != readlen:
                        print("Note: Reads in the input dataset are of variable lengths.")
                    j += 1
                except ValueError:
                    invalid_lines += 1
                if j > 10:
                    break

            position_dict = {}
            print("column\tcount\tmin\tmax\tsum\tmean\tQ1\tmed\tQ3\tIQR\tlW\trW", file=fout)
            for line in open(infile_name):
                line = line.strip()
                if not (line) or line.startswith("#") or line.startswith(">"):
                    continue
                elems = line.split()
                if position_dict == {}:
                    for pos in range(readlen):
                        position_dict[pos] = [0] * QUAL_UPPER_BOUND
                if len(elems) != readlen:
                    invalid_lines += 1
                    continue
                for ind, item in enumerate(elems):
                    try:
                        item = int(item)
                        position_dict[ind][item] += 1
                    except Exception:
                        pass

            invalid_positions = 0
            for pos in position_dict:
                carr = position_dict[pos]  # count array for position pos
                total = sum(carr)  # number of bases found in this column.
                med_elem = int(round(total / 2.0))
                lowest = None  # Lowest quality score value found in this column.
                highest = None  # Highest quality score value found in this column.
                median = None  # Median quality score value found in this column.
                qsum = 0.0  # Sum of quality score values for this column.
                q1 = None  # 1st quartile quality score.
                q3 = None  # 3rd quartile quality score.
                q1_elem = int(round((total + 1) / 4.0))
                q3_elem = int(round((total + 1) * 3 / 4.0))

                try:
                    for ind, cnt in enumerate(carr):
                        qsum += ind * cnt

                        if cnt != 0:
                            highest = ind

                        if lowest is None and cnt != 0:  # first non-zero count
                            lowest = ind

                        if q1 is None:
                            if sum(carr[: ind + 1]) >= q1_elem:
                                q1 = ind

                        if median is None:
                            if sum(carr[: ind + 1]) < med_elem:
                                continue
                            median = ind
                            if total % 2 == 0:  # even number of elements
                                median2 = median
                                if sum(carr[: ind + 1]) < med_elem + 1:
                                    for ind2, elem in enumerate(carr[ind + 1 :]):
                                        if elem != 0:
                                            median2 = ind + ind2 + 1
                                            break
                                median = (median + median2) / 2.0

                        if q3 is None:
                            if sum(carr[: ind + 1]) >= q3_elem:
                                q3 = ind

                    mean = qsum / total  # Mean quality score value for this column.
                    iqr = q3 - q1
                    left_whisker = max(q1 - 1.5 * iqr, lowest)
                    right_whisker = min(q3 + 1.5 * iqr, highest)

                    print(
                        "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s"
                        % (
                            pos + 1,
                            total,
                            lowest,
                            highest,
                            qsum,
                            mean,
                            q1,
                            median,
                            q3,
                            iqr,
                            left_whisker,
                            right_whisker,
                        ),
                        file=fout,
                    )
                except Exception:
                    invalid_positions += 1
                    nullvals = ["NA"] * 11
                    print("%s\t%s" % (pos + 1, "\t".join(nullvals)), file=fout)

    if invalid_lines:
        print("Skipped %d reads as invalid." % invalid_lines)
    if invalid_positions:
        print("Skipped stats computation for %d read positions." % invalid_positions)


if __name__ == "__main__":
    __main__()