1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
#!/usr/bin/env python3
#
# Copyright (C) 2010 Institute for Computational Biomedicine,
# Weill Medical College of Cornell University
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
import getopt
import os
import stat
import sys
from goby.Reads import CompactReads
from goby.utils import commify
def usage():
print("usage:", sys.argv[0], "[-h|--help] [-v|--verbose] <filename>")
def main():
verbose = False
try:
opts, args = getopt.getopt(sys.argv[1:], "hv", ["help", "verbose"])
except getopt.GetoptError as err:
print(str(err),file=sys.stderr)
usage()
sys.exit(1)
# Collect options
for opt, arg in opts:
if opt in ("-h", "--help"):
usage()
sys.exit()
elif opt in ("-v", "--verbose"):
verbose = True
if len(args) != 1:
usage()
sys.exit(2)
filename = args[0]
print("Compact reads filename = %s" % filename)
print()
filesize = os.stat(filename)[stat.ST_SIZE]
number_of_entries = 0
number_of_identifiers = 0
number_of_descriptions = 0
number_of_sequences = 0
number_of_sequence_pairs = 0
number_of_quality_scores = 0
number_of_quality_score_pairs = 0
min_read_length = sys.maxsize
max_read_length = -sys.maxsize - 1
total_read_length = 0
total_read_length_pair = 0
reads_reader = CompactReads(filename)
for entry in reads_reader:
number_of_entries += 1
read_length = entry.read_length
total_read_length += read_length
total_read_length_pair += entry.read_length_pair
if entry.HasField("read_identifier"):
number_of_identifiers += 1
if entry.HasField("description"):
number_of_descriptions += 1
if entry.HasField("sequence"):
number_of_sequences += 1
if entry.HasField("sequence_pair"):
number_of_sequence_pairs += 1
if entry.HasField("quality_scores"):
number_of_quality_scores += 1
if entry.HasField("quality_scores_pair"):
number_of_quality_score_pairs += 1
min_read_length = min(min_read_length, read_length)
max_read_length = max(max_read_length, read_length)
print("Average bytes per entry: %s" % commify(filesize / float(number_of_entries)))
print("Average bytes per base: %s" % commify(filesize / float(total_read_length)))
print("Has identifiers = %s (%s)" % (number_of_identifiers > 0, commify(number_of_identifiers)))
print("Has descriptions = %s (%s)" % (number_of_descriptions > 0, commify(number_of_descriptions)))
print("Has sequences = %s (%s)" % (number_of_sequences > 0, commify(number_of_sequences)))
print("Has sequence pairs = %s (%s)" % (number_of_sequence_pairs > 0, commify(number_of_sequence_pairs)))
print("Has quality scores = %s (%s)" % (number_of_quality_scores > 0, commify(number_of_quality_scores)))
print("Has quality score pairs = %s (%s)" % (
number_of_quality_score_pairs > 0, commify(number_of_quality_score_pairs)))
print("Number of entries = %s" % commify(number_of_entries))
print("Min read length = %s" % commify(min_read_length))
print("Max read length = %s" % commify(max_read_length))
print("Avg read length = %s" % commify(total_read_length / float(number_of_entries)))
print("Avg read pair length = %s" % commify(total_read_length_pair / float(number_of_entries)))
if __name__ == "__main__":
main()
|
