Mercurial > repos > cstrittmatter > test_eurl_vtec_wgs_pt

comparison scripts/fastq_positional_quality_trimming.py @ 0:965517909457 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 15239f1674081ab51ab8dd75a9a40cf1bfaa93e8
author cstrittmatter
date Wed, 22 Jan 2020 08:41:44 -0500
parents
children 0cbed1c0a762
comparison
equal deleted inserted replaced
-1:000000000000 0:965517909457
1 # -*- coding: utf-8 -*-
2 """
3 This tool trims paired-end FASTQ files on the basis of quality score or left/right position, retaining mate integrity.
4 Reads without mate after filtering are saved in a separate output file.
5 """
6
7 import math
8 import optparse
9 import sys
10
11 def average(values):
12 """ Arithmetic mean of a list of values """
13 return math.fsum(values) / len(values) if len(values) else float('nan')
14
15
16 def phred2sanger(phred_scores):
17 """ Convert an array of Phred quality scores (integers in [0, 93]) to a Sanger-encoded quality string"""
18 return ''.join([chr(score + 33) for score in phred_scores])
19
20
21 def sanger2phred(sanger_string):
22 """ Convert a Sanger-encoded quality string to an array of Phred quality scores"""
23 return [ord(ch) - 33 for ch in sanger_string]
24
25
26 def trimming(sequ, qual, maxlengthtrim, lefttrim, righttrim, minqualtrim, avgqualtrim):
27 """ Trimming of sequence and quality of a read"""
28 # Maximum length trimming
29 if maxlengthtrim != -1:
30 if len(sequ) > maxlengthtrim:
31 sequ = sequ[: maxlengthtrim]
32 qual = qual[: maxlengthtrim]
33 # Left- and right-side trimming
34 if righttrim == 0:
35 sequ = sequ[lefttrim :]
36 qual = qual[lefttrim :]
37 else:
38 sequ = sequ[lefttrim : -righttrim]
39 qual = qual[lefttrim : -righttrim]
40 # Minimum quality right-side trimming
41 while len(sequ) and qual[-1] < minqualtrim:
42 qual = qual[:-1]
43 sequ = sequ[:-1]
44 # Average quality right-side trimming
45 while len(sequ) and average(qual) < avgqualtrim:
46 qual = qual[:-1]
47 sequ = sequ[:-1]
48 return sequ, qual
49
50
51 def __main__():
52 parser = optparse.OptionParser()
53 parser.add_option('-1', '--input1', dest='input1', help='forward or single-end reads file in Sanger FASTQ format')
54 parser.add_option('-2', '--input2', dest='input2', help='reverse reads file in Sanger FASTQ format')
55 parser.add_option('--maxlt', dest='maxlengthtrim', type='int', default=400, help='maximum length trimming')
56 parser.add_option('--lt', dest='lefttrim', type='int', default=0, help='left-side trimming')
57 parser.add_option('--rt', dest='righttrim', type='int', default=0, help='right-side trimming')
58 parser.add_option('--minqt', dest='minqualtrim', type='int', default=15, help='minimum quality right-side trimming')
59 parser.add_option('--avgqt', dest='avgqualtrim', type='float', default=20, help='average quality right-side trimming')
60 parser.add_option('--minlf', dest='minlen', type='int', default=25, help='minimum length filtering')
61 parser.add_option('--trimmed1', dest='trimmed1', help='trimmed forward FASTQ file')
62 parser.add_option('--trimmed2', dest='trimmed2', help='trimmed reverse FASTQ file')
63 parser.add_option('--trimmedunpaired', dest='trimmedunpaired', help='trimmed unpaired FASTQ file')
64 parser.add_option('--log', dest='logfile', help='log file')
65 (options, args) = parser.parse_args()
66 if len(args) > 0:
67 parser.error('Wrong number of arguments')
68
69 maxlengthtrim = options.maxlengthtrim
70 lefttrim = options.lefttrim
71 righttrim = options.righttrim
72 minqualtrim = options.minqualtrim
73 avgqualtrim = options.avgqualtrim
74 minlen = options.minlen
75
76 total_reads = 0
77 discarded_reads = 0
78 forward = open(options.input1)
79 if options.input2:
80 paired = True
81 passing_paired_reads = 0
82 passing_unpaired_reads = 0
83 reverse = open(options.input2)
84 trimmed_reverse = open(options.trimmed2, 'w')
85 trimmed_unpaired = open(options.trimmedunpaired, 'w')
86 else:
87 paired = False
88 passing_reads = 0
89 trimmed_forward = open(options.trimmed1, 'w')
90 log = open(options.logfile, 'w')
91 try:
92 while True:
93 headL = forward.next().rstrip()
94 sequL = forward.next().rstrip()
95 commL = forward.next().rstrip()
96 sangL = forward.next().rstrip()
97 qualL = sanger2phred(sangL)
98 trimmed_sequL, trimmed_qualL = trimming(sequL, qualL, maxlengthtrim, lefttrim, righttrim, minqualtrim, avgqualtrim)
99 if paired:
100 try:
101 headR = reverse.next().rstrip()
102 sequR = reverse.next().rstrip()
103 commR = reverse.next().rstrip()
104 sangR = reverse.next().rstrip()
105 except StopIteration:
106 sys.exit('Reverse FASTQ file contain less reads than forward FASTQ file.')
107 qualR = sanger2phred(sangR)
108 trimmed_sequR, trimmed_qualR = trimming(sequR, qualR, maxlengthtrim, lefttrim, righttrim, minqualtrim, avgqualtrim)
109 # Filter by residual length
110 if paired:
111 if len(trimmed_sequL) >= minlen and len(trimmed_sequR) >= minlen:
112 trimmed_forward.write(headL + '\n' + trimmed_sequL + '\n' + commL + '\n' + phred2sanger(trimmed_qualL) + '\n')
113 trimmed_reverse.write(headR + '\n' + trimmed_sequR + '\n' + commR + '\n' + phred2sanger(trimmed_qualR) + '\n')
114 passing_paired_reads += 1
115 elif len(trimmed_sequL) >= minlen and len(trimmed_sequR) < minlen:
116 trimmed_unpaired.write(headL + '\n' + trimmed_sequL + '\n' + commL + '\n' + phred2sanger(trimmed_qualL) + '\n')
117 passing_unpaired_reads += 1
118 elif len(trimmed_sequL) < minlen and len(trimmed_sequR) >= minlen:
119 trimmed_unpaired.write(headR + '\n' + trimmed_sequR + '\n' + commR + '\n' + phred2sanger(trimmed_qualR) + '\n')
120 passing_unpaired_reads += 1
121 else:
122 discarded_reads += 1
123 else:
124 if len(trimmed_sequL) >= minlen:
125 trimmed_forward.write(headL + '\n' + trimmed_sequL + '\n' + commL + '\n' + phred2sanger(trimmed_qualL) + '\n')
126 passing_reads += 1
127 else:
128 discarded_reads += 1
129 total_reads += 1
130 except StopIteration:
131 if paired:
132 try:
133 reverse.next()
134 except StopIteration:
135 log.write("Total paired reads : %d\n" % total_reads)
136 log.write("Passing paired reads : %d\n" % passing_paired_reads)
137 log.write("Passing unpaired reads : %d\n" % passing_unpaired_reads)
138 log.write("Discarded paired reads : %d\n" % discarded_reads)
139 else:
140 sys.exit('Forward FASTQ file contain less reads than reverse FASTQ file.')
141 else:
142 log.write("Total reads : %d\n" % total_reads)
143 log.write("Passing reads : %d\n" % passing_reads)
144 log.write("Discarded reads : %d\n" % discarded_reads)
145 finally:
146 forward.close()
147 trimmed_forward.close()
148 log.close()
149 if paired:
150 reverse.close()
151 trimmed_reverse.close()
152 trimmed_unpaired.close()
153
154
155 if __name__ == "__main__":
156 __main__()