annotate get_contigs_from_re_archive.py @ 8:2557cad81607 draft

Uploaded
author petr-novak
date Wed, 03 May 2023 11:12:41 +0000
parents 814cba36e435
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
1 #!/usr/bin/env python
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
2 '''
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
3 parse .aln file - output from cap3 program. Output is fasta file and
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
4 profile file
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
5 '''
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
6 import argparse
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
7 import re
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
8 import zipfile
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
9 import tempfile
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
10 import textwrap
1
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
11 import os
0
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
12
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
13 def parse_args():
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
14 '''Argument parsin'''
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
15 description = """
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
16 parsing cap3 assembly aln output
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
17 """
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
18
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
19 parser = argparse.ArgumentParser(
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
20 description=description, formatter_class=argparse.RawTextHelpFormatter)
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
21 parser.add_argument('-re',
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
22 '--re_file',
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
23 default=None,
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
24 required=True,
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
25 help="RepeatExlorer archive or directory",
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
26 type=str,
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
27 action='store')
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
28 parser.add_argument('-f',
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
29 '--fasta',
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
30 default=None,
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
31 required=True,
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
32 help="fasta output file name",
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
33 type=str,
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
34 action='store')
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
35 parser.add_argument('-m',
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
36 '--min_coverage',
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
37 default=5,
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
38 required=False,
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
39 help="minimum contig coverage",
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
40 type=int,
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
41 action="store")
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
42 parser.add_argument('-L',
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
43 '--min_contig_length',
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
44 default=50,
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
45 required=False,
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
46 help="minimum contig length",
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
47 type=int,
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
48 action="store")
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
49 return parser.parse_args()
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
50
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
51
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
52 def get_header(f):
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
53 aln_header = " . : . : . : . : . : . :"
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
54 contig_lead = "******************"
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
55 aln_start = -1
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
56 while True:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
57 line = f.readline()
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
58 if not line:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
59 return None, None
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
60 if line[0:18] == contig_lead:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
61 line2 = f.readline()
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
62 else:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
63 continue
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
64 if aln_header in line2:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
65 aln_start = line2.index(aln_header)
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
66 break
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
67 contig_name = line.split()[1] + line.split()[2]
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
68 return contig_name, aln_start
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
69
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
70
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
71 def segment_start(f):
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
72 pos = f.tell()
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
73 line = f.readline()
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
74 # detect next contig or end of file
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
75 if "********" in line or line == "" or "Number of segment pairs = " in line:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
76 segment = False
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
77 else:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
78 segment = True
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
79 f.seek(pos)
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
80 return segment
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
81
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
82
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
83 def get_segment(f, seq_start):
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
84 if not segment_start(f):
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
85 return None, None
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
86 aln = []
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
87 while True:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
88 line = f.readline()
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
89 if ". : . :" in line:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
90 continue
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
91 if "__________" in line:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
92 consensus = f.readline().rstrip('\n')[seq_start:]
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
93 f.readline() # empty line
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
94 break
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
95 else:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
96 aln.append(line.rstrip('\n')[seq_start:])
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
97 return aln, consensus
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
98
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
99
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
100 def aln2coverage(aln):
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
101 coverage = [0] * len(aln[0])
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
102 for a in aln:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
103 for i, c in enumerate(a):
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
104 if c not in " -":
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
105 coverage[i] += 1
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
106 return coverage
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
107
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
108
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
109 def read_contig(f, seq_start):
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
110 contig = ""
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
111 coverage = []
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
112 while True:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
113 aln, consensus = get_segment(f, seq_start)
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
114 if aln:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
115 contig += consensus
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
116 coverage += aln2coverage(aln)
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
117 else:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
118 break
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
119 return contig, coverage
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
120
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
121
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
122 def remove_gaps(consensus, coverage):
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
123 if "-" not in consensus:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
124 return consensus, coverage
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
125 new_coverage = [
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
126 cov for cons, cov in zip(consensus, coverage) if cons != "-"
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
127 ]
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
128 new_consensus = consensus.replace("-", "")
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
129 return new_consensus, new_coverage
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
130
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
131
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
132 def extract_contigs_from_re_archive(archive, aln_output):
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
133 with zipfile.ZipFile(archive, 'r') as zip_object, open(aln_output,
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
134 'w') as fout:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
135 flist = zip_object.infolist()
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
136 for fn in flist:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
137 if re.match('seqclust.+[.]aln$', fn.filename):
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
138 with zip_object.open(fn.filename) as aln:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
139 for l in aln:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
140 fout.write(l.decode('utf-8'))
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
141 return aln_output
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
142
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
143 def read_tarean_fasta(fobj):
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
144 ids = []
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
145 s = []
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
146 for i in fobj:
1
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
147 if isinstance(i, str):
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
148 ii = i
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
149 else:
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
150 ii = i.decode('utf-8')
0
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
151 if ii[0] == ">":
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
152 ids.append(ii)
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
153 s.append("")
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
154 else:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
155 s[-1] = s[-1] + ii.strip()
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
156 return ids, s
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
157
1
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
158
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
159 def extract_contigs_from_re_directory(re_dir, aln_output):
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
160 with open(aln_output, 'w') as fout:
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
161 for subdir, dirs, files in os.walk(re_dir):
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
162 for fn in files:
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
163 fn_full = subdir + os.sep + fn
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
164 if re.match('^.+seqclust.+[.]aln$', fn_full):
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
165 print(fn_full)
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
166 with open(fn_full) as aln:
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
167 for l in aln:
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
168 fout.write(l)
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
169 return aln_output
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
170
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
171
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
172
0
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
173 def extract_tarean_contigs_from_re_archive(archive):
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
174 with zipfile.ZipFile(archive, 'r') as zip_object:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
175 flist = zip_object.infolist()
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
176 seqs_all = []
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
177 ids_all = []
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
178 for fn in flist:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
179 if re.match("seqclust.+dir_CL[0-9]+[/]tarean_contigs.fasta", fn.filename):
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
180 print(fn.filename)
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
181 with zip_object.open(fn.filename) as fobj:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
182 ids, seqs = read_tarean_fasta(fobj)
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
183 # wrap sequences
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
184 seqs = ["\n".join(textwrap.wrap(s + s, 80)) for s in seqs]
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
185 seqs_all += seqs
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
186 ids_all += ids
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
187 return ids_all, seqs_all
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
188
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
189
1
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
190 def extract_tarean_contigs_from_re_directory(re_dir):
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
191 seqs_all = []
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
192 ids_all = []
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
193 for subdir, dirs, files in os.walk(re_dir):
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
194 for fn in files:
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
195 fn_full = subdir + os.sep + fn
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
196 if re.match("^.+seqclust.+dir_CL[0-9]+[/]tarean_contigs.fasta", fn_full):
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
197 print (fn_full)
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
198 with open(fn_full) as fobj:
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
199 ids, seqs = read_tarean_fasta(fobj)
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
200 # wrap sequences
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
201 seqs = ["\n".join(textwrap.wrap(s + s, 80)) for s in seqs]
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
202 seqs_all += seqs
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
203 ids_all += ids
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
204 return ids_all, seqs_all
0
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
205
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
206 def filter_contigs(consensus, coverage, min_coverage=5):
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
207 x = "".join([
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
208 s if cov >= min_coverage else " "
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
209 for s, cov in zip(consensus, coverage)
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
210 ]).strip()
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
211 consensus_N = "\n".join(textwrap.wrap(x.replace(" ", "N"),80))
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
212 return consensus_N
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
213
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
214
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
215 def main():
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
216 args = parse_args()
1
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
217 if os.path.isdir(args.re_file):
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
218 # extract from directory
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
219 ids, seqs = extract_tarean_contigs_from_re_directory(args.re_file)
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
220 aln_file = extract_contigs_from_re_directory(
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
221 args.re_file,
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
222 tempfile.NamedTemporaryFile().name)
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
223 else:
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
224 # extract aln from archive
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
225 ids, seqs = extract_tarean_contigs_from_re_archive(args.re_file)
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
226 aln_file = extract_contigs_from_re_archive(
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
227 args.re_file,
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
228 tempfile.NamedTemporaryFile().name)
814cba36e435 Uploaded
mvdbeek
parents: 0
diff changeset
229
0
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
230 with open(aln_file, 'r') as f1, open(args.fasta, 'w') as ffasta:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
231 while True:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
232 contig_name, seq_start = get_header(f1)
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
233 if contig_name:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
234 consensus, coverage = remove_gaps(*read_contig(f1, seq_start))
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
235 clean_consensus = filter_contigs(consensus, coverage,
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
236 args.min_coverage)
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
237 if len(clean_consensus) >= args.min_contig_length:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
238 ffasta.write(">{}\n".format(contig_name))
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
239 ffasta.write("{}\n".format(clean_consensus))
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
240 else:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
241 break
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
242
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
243 # write tarean sequences:
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
244 for i, s in zip(ids, seqs):
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
245 ffasta.write(i)
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
246 ffasta.write(s + "\n")
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
247
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
248
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
249
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
250 if __name__ == "__main__":
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
251
ea6a3059a6af Uploaded
petr-novak
parents:
diff changeset
252 main()