mirror of
https://github.com/galaxyproject/galaxy.git
synced 2026-09-24 16:30:27 +08:00
Fix solid_qual_stats.py
This commit is contained in:
mvdbeek
@@ -26,109 +26,110 @@ def unzip(filename):
|
||||
|
||||
def __main__():
|
||||
infile_score_name = sys.argv[1].strip()
|
||||
fout = open(sys.argv[2].strip(), 'r+w')
|
||||
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
|
||||
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
|
||||
for line in open(infile_name):
|
||||
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
|
||||
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 k, line in enumerate(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
|
||||
position_dict = {}
|
||||
print("column\tcount\tmin\tmax\tsum\tmean\tQ1\tmed\tQ3\tIQR\tlW\trW", file=fout)
|
||||
for k, line in enumerate(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))
|
||||
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
|
||||
try:
|
||||
for ind, cnt in enumerate(carr):
|
||||
qsum += ind * cnt
|
||||
|
||||
if cnt != 0:
|
||||
highest = ind
|
||||
if cnt != 0:
|
||||
highest = ind
|
||||
|
||||
if lowest is None and cnt != 0: # first non-zero count
|
||||
lowest = 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 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 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
|
||||
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)
|
||||
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)
|
||||
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)
|
||||
|
||||
Reference in New Issue
Block a user