annotate fasta_interlacer.py @ 13:f12eb3896842 draft

Uploaded
author petr-novak
date Fri, 07 Feb 2020 02:25:47 -0500
parents d397f5a85464
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
1 #!/usr/bin/env python
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
2 ''' interlacing two fastq sequences'''
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
3 import sys
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
4
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
5 def readSingleSeq(file):
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
6 ''' read single seq from fasta file'''
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
7 line = file.readline()
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
8 if not line:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
9 return False # end of file
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
10 if line[0] != ">":
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
11 raise Exception("no header on the first line")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
12 seqname = line[1:].strip()
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
13 seq = ""
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
14 # read sequences
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
15 while True:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
16 last_pos = file.tell()
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
17 line = file.readline()
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
18 if not line:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
19 break
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
20 if line[0] == ">":
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
21 file.seek(last_pos)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
22 break
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
23 seq = seq + line.strip()
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
24 return {'name': seqname, 'sequence': seq}
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
25
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
26
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
27 def writeSingleSeq(fileobject, seq):
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
28 ''' write single sequence to fasta file'''
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
29 fileobject.write(">")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
30 fileobject.write(seq['name'] + "\n")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
31 fileobject.write(seq['sequence'] + "\n")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
32
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
33
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
34 def main():
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
35 from optparse import OptionParser
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
36 parser = OptionParser()
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
37 parser.add_option("-a",
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
38 "--fasta_file_A",
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
39 dest="seqfileA",
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
40 help="input sequences in fasta format")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
41 parser.add_option("-b",
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
42 "--fasta_file_B",
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
43 dest="seqfileB",
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
44 help="input sequences in fasta format")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
45 parser.add_option("-p",
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
46 "--fasta_file_pairs",
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
47 dest="seqfile_pairs",
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
48 help="output file with paired sequences")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
49 parser.add_option("-x",
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
50 "--fasta_file_singles",
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
51 dest="seqfile_singles",
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
52 help="output file with single sequences")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
53 options, args = parser.parse_args()
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
54
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
55 # Input files
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
56 fA = open(options.seqfileA, 'r')
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
57 fB = open(options.seqfileB, 'r')
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
58 # Output files
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
59 if options.seqfile_pairs:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
60 fPairs = open(options.seqfile_pairs, 'w')
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
61 else:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
62 fPairs = open(options.seqfileA + ".pairs", 'w')
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
63
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
64 if options.seqfile_singles:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
65 single = open(options.seqfile_singles, "w")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
66 else:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
67 single = open(options.seqfileA + ".single", "w")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
68
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
69
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
70 sA1 = readSingleSeq(fA)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
71 sB1 = readSingleSeq(fB)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
72 if not sA1 or not sB1:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
73 raise Exception("\nEmpty sequence on input, nothing to interlace!\n")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
74 charA = sA1['name'][-1]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
75 charB = sB1['name'][-1]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
76 # validate sequence names
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
77 if charA == charB:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
78 sys.stderr.write(
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
79 "last character of sequence id must be used for distinguishing pairs!")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
80 exit(1)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
81 # check first thousand!
4
d397f5a85464 Uploaded
petr-novak
parents: 0
diff changeset
82 for i in range(3):
0
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
83 seqA = readSingleSeq(fA)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
84 seqB = readSingleSeq(fB)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
85 if (not seqA) or (not seqB):
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
86 # end of file:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
87 if i == 0:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
88 sys.stderr.write("input file is empty")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
89 exit(1)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
90 else:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
91 break
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
92 if seqA['name'][-1] == charA and seqB['name'][-1] == charB:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
93 continue
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
94 else:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
95 sys.stderr.write(
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
96 "last character of sequence id must be used for distinguishing pairs!")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
97 exit(1)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
98
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
99 fA.seek(0)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
100 fB.seek(0)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
101
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
102 buffA = {}
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
103 buffB = {}
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
104 buffA_names = []
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
105 buffB_names = []
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
106
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
107 while True:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
108 seqA = readSingleSeq(fA)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
109 seqB = readSingleSeq(fB)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
110 if not seqA and not seqB:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
111 break # end of file
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
112
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
113 ## validation and direct checking only if not end of files
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
114 if seqA and seqB:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
115 #validate:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
116 if not (seqA['name'][-1] == charA and seqB['name'][-1] == charB):
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
117 sys.stderr.write(
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
118 "last character of sequence id must be used for distinguishing pairs!")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
119 exit(1)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
120
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
121 # check if current seqs are pairs
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
122 if seqA['name'][:-1] == seqB['name'][:-1]:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
123 writeSingleSeq(fPairs, seqA)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
124 writeSingleSeq(fPairs, seqB)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
125 continue
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
126
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
127 ### compare whith buffers
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
128 ### seqA vs buffB
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
129 if seqA:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
130 if seqA["name"][:-1] in buffB:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
131 writeSingleSeq(fPairs, seqA)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
132 seqtmp = {"name": seqA["name"][:-1] + charB,
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
133 "sequence": buffB[seqA["name"][:-1]]}
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
134 writeSingleSeq(fPairs, seqtmp)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
135 # can I empty buffA ???
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
136 for i in buffA_names:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
137 seqtmp = {"name": i + charA, "sequence": buffA[i]}
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
138 writeSingleSeq(single, seqtmp)
4
d397f5a85464 Uploaded
petr-novak
parents: 0
diff changeset
139 buffA = {}
d397f5a85464 Uploaded
petr-novak
parents: 0
diff changeset
140 buffA_names = []
0
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
141
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
142 j = 0
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
143 for i in buffB_names:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
144 seqtmp = {"name": i + charB, "sequence": buffB[i]}
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
145 del buffB[i]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
146 j += 1
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
147 if i == seqA["name"][:-1]:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
148 del buffB_names[0:j]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
149 break
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
150 else:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
151 writeSingleSeq(single, seqtmp)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
152 else:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
153 buffA[seqA["name"][:-1]] = seqA['sequence']
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
154 buffA_names.append(seqA["name"][:-1])
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
155
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
156 ### seqA vs buffB
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
157 if seqB:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
158 if seqB["name"][:-1] in buffA:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
159 seqtmp = {"name": seqB["name"][:-1] + charA,
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
160 "sequence": buffA[seqB["name"][:-1]]}
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
161 writeSingleSeq(fPairs, seqtmp)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
162 writeSingleSeq(fPairs, seqB)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
163 # can I empty buffB ???
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
164 for i in buffB_names:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
165 seqtmp = {"name": i + charB, "sequence": buffB[i]}
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
166 writeSingleSeq(single, seqtmp)
4
d397f5a85464 Uploaded
petr-novak
parents: 0
diff changeset
167 buffB = {}
d397f5a85464 Uploaded
petr-novak
parents: 0
diff changeset
168 buffB_names = []
0
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
169
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
170 j = 0
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
171 for i in buffA_names:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
172 seqtmp = {"name": i + charA, "sequence": buffA[i]}
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
173 del buffA[i]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
174 j += 1
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
175 if i == seqB["name"][:-1]:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
176 del buffA_names[0:j]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
177 break
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
178 else:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
179 writeSingleSeq(single, seqtmp)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
180
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
181 else:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
182 buffB[seqB["name"][:-1]] = seqB['sequence']
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
183 buffB_names.append(seqB["name"][:-1])
4
d397f5a85464 Uploaded
petr-novak
parents: 0
diff changeset
184
d397f5a85464 Uploaded
petr-novak
parents: 0
diff changeset
185 fA.close()
d397f5a85464 Uploaded
petr-novak
parents: 0
diff changeset
186 fB.close()
d397f5a85464 Uploaded
petr-novak
parents: 0
diff changeset
187 fPairs.close()
d397f5a85464 Uploaded
petr-novak
parents: 0
diff changeset
188
d397f5a85464 Uploaded
petr-novak
parents: 0
diff changeset
189 # write rest of singles:
0
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
190 for i in buffA:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
191 seqtmp = {"name": i + charA, "sequence": buffA[i]}
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
192 writeSingleSeq(single, seqtmp)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
193 for i in buffB:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
194 seqtmp = {"name": i + charB, "sequence": buffB[i]}
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
195 writeSingleSeq(single, seqtmp)
4
d397f5a85464 Uploaded
petr-novak
parents: 0
diff changeset
196 single.close()
0
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
197
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
198
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
199 if __name__ == "__main__":
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
200 main()