|
3
|
1 #!/usr/bin/python -tt
|
|
4
|
2 """
|
|
|
3 pal_filter
|
|
|
4 https://github.com/graemefox/pal_filter
|
|
|
5
|
|
|
6 Graeme Fox - 03/03/2016 - graeme.fox@manchester.ac.uk
|
|
|
7 Tested on 64-bit Ubuntu, with Python 2.7
|
|
|
8
|
|
|
9 ~~~~~~~~~~~~~~~~~~~
|
|
|
10 PROGRAM DESCRIPTION
|
|
|
11
|
|
|
12 Program to pick optimum loci from the output of pal_finder_v0.02.04
|
|
|
13
|
|
|
14 This program can be used to filter output from pal_finder and choose the
|
|
|
15 'optimum' loci.
|
|
|
16
|
|
|
17 For the paper referncing this workflow, see Griffiths et al.
|
|
|
18 (unpublished as of 15/02/2016) (sarah.griffiths-5@postgrad.manchester.ac.uk)
|
|
|
19
|
|
|
20 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
21 This program also contains a quality-check method to improve the rate of PCR
|
|
|
22 success. For this QC method, paired end reads are assembled using
|
|
|
23 PANDAseq so you must have PANDAseq installed.
|
|
|
24
|
|
|
25 For the paper referencing this assembly-QC method see Fox et al.
|
|
|
26 (unpublished as of 15/02/2016) (graeme.fox@manchester.ac.uk)
|
|
|
27
|
|
|
28 For best results in PCR for marker development, I suggest enabling all the
|
|
|
29 filter options AND the assembly based QC
|
|
|
30
|
|
|
31 ~~~~~~~~~~~~
|
|
|
32 REQUIREMENTS
|
|
|
33
|
|
|
34 Must have Biopython installed (www.biopython.org).
|
|
|
35
|
|
|
36 If you with to perform the assembly QC step, you must have:
|
|
|
37 PandaSeq (https://github.com/neufeld/pandaseq)
|
|
|
38 PandaSeq must be in your $PATH / able to run from anywhere
|
|
|
39
|
|
|
40 ~~~~~~~~~~~~~~~~
|
|
|
41 REQUIRED OPTIONS
|
|
|
42
|
|
|
43 -i forward_paired_ends.fastQ
|
|
|
44 -j reverse_paired_ends.fastQ
|
|
|
45 -p pal_finder output - by default pal_finder names this "x_PAL_summary.txt"
|
|
|
46
|
|
|
47 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
48 BY DEFAULT THIS PROGRAM DOES NOTHING. ENABLE SOME OF THE OPTIONS BELOW.
|
|
|
49 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
50 NON-REQUIRED OPTIONS
|
|
|
51
|
|
|
52 -assembly: turn on the pandaseq assembly QC step
|
|
|
53
|
|
|
54 -primers: filter microsatellite loci to just those which have primers designed
|
|
|
55
|
|
|
56 -occurrences: filter microsatellite loci to those with primers
|
|
|
57 which appear only once in the dataset
|
|
|
58
|
|
|
59 -rankmotifs: filter microsatellite loci to just those with perfect motifs.
|
|
|
60 Rank the output by size of motif (largest first)
|
|
|
61
|
|
|
62 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
63 For repeat analysis, the following extra non-required options may be useful:
|
|
|
64
|
|
|
65 Since PandaSeq Assembly, and fastq -> fasta conversion are slow, do them the
|
|
|
66 first time, generate the files and then skip either, or both steps with
|
|
|
67 the following:
|
|
|
68
|
|
|
69 -a: skip assembly step
|
|
|
70 -c: skip fastq -> fasta conversion step
|
|
|
71
|
|
|
72 Just make sure to keep the assembled/converted files in the correct directory
|
|
|
73 with the correct filename(s)
|
|
|
74
|
|
|
75 ~~~~~~~~~~~~~~~
|
|
|
76 EXAMPLE USAGE:
|
|
|
77
|
|
|
78 pal_filtery.py -i R1.fastq -j R2.fastq
|
|
|
79 -p pal_finder_output.tabular -primers -occurrences -rankmotifs -assembly
|
|
|
80
|
|
|
81 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
82 """
|
|
3
|
83 import Bio, subprocess, argparse, csv, os, re, time
|
|
|
84 from Bio import SeqIO
|
|
4
|
85 __version__ = "1.0.0"
|
|
|
86 ############################################################
|
|
|
87 # Function List #
|
|
|
88 ############################################################
|
|
|
89 def ReverseComplement1(seq):
|
|
|
90 """
|
|
|
91 take a nucleotide sequence and reverse-complement it
|
|
|
92 """
|
|
|
93 seq_dict = {'A':'T','T':'A','G':'C','C':'G'}
|
|
|
94 return "".join([seq_dict[base] for base in reversed(seq)])
|
|
|
95
|
|
|
96 def fastq_to_fasta(input_file, wanted_set):
|
|
|
97 """
|
|
|
98 take a file in fastq format, convert to fasta format and filter on
|
|
|
99 the set of sequences that we want to keep
|
|
|
100 """
|
|
|
101 file_name = os.path.splitext(os.path.basename(input_file))[0]
|
|
|
102 with open(file_name + "_filtered.fasta", "w") as out:
|
|
|
103 for record in SeqIO.parse(input_file, "fastq"):
|
|
|
104 ID = str(record.id)
|
|
|
105 SEQ = str(record.seq)
|
|
|
106 if ID in wanted_set:
|
|
|
107 out.write(">" + ID + "\n" + SEQ + "\n")
|
|
|
108
|
|
|
109 def strip_barcodes(input_file, wanted_set):
|
|
|
110 """
|
|
|
111 take fastq data containing sequencing barcodes and strip the barcode
|
|
|
112 from each sequence. Filter on the set of sequences that we want to keep
|
|
|
113 """
|
|
|
114 file_name = os.path.splitext(os.path.basename(input_file))[0]
|
|
|
115 with open(file_name + "_adapters_removed.fasta", "w") as out:
|
|
|
116 for record in SeqIO.parse(input_file, "fasta"):
|
|
|
117 match = re.search(r'\S*:', record.id)
|
|
|
118 if match:
|
|
|
119 correct = match.group().rstrip(":")
|
|
|
120 else:
|
|
|
121 correct = str(record.id)
|
|
|
122 SEQ = str(record.seq)
|
|
|
123 if correct in wanted_set:
|
|
|
124 out.write(">" + correct + "\n" + SEQ + "\n")
|
|
|
125
|
|
|
126 ############################################################
|
|
|
127 # MAIN PROGRAM #
|
|
|
128 ############################################################
|
|
|
129 print "\n~~~~~~~~~~"
|
|
|
130 print "pal_filter"
|
|
|
131 print "~~~~~~~~~~"
|
|
|
132 print "Version: " + __version__
|
|
|
133 time.sleep(1)
|
|
|
134 print "\nFind the optimum loci in your pal_finder output and increase "\
|
|
|
135 "the rate of successful microsatellite marker development"
|
|
|
136 print "\nSee Griffiths et al. (currently unpublished) for more details......"
|
|
|
137 time.sleep(2)
|
|
3
|
138
|
|
|
139 # Get values for all the required and optional arguments
|
|
|
140
|
|
|
141 parser = argparse.ArgumentParser(description='pal_filter')
|
|
|
142 parser.add_argument('-i','--input1', help='Forward paired-end fastq file', \
|
|
|
143 required=True)
|
|
|
144
|
|
|
145 parser.add_argument('-j','--input2', help='Reverse paired-end fastq file', \
|
|
|
146 required=True)
|
|
|
147
|
|
|
148 parser.add_argument('-p','--pal_finder', help='Output from pal_finder ', \
|
|
|
149 required=True)
|
|
|
150
|
|
4
|
151 parser.add_argument('-assembly', help='Perform the PandaSeq based QC', \
|
|
|
152 action='store_true')
|
|
3
|
153
|
|
|
154 parser.add_argument('-a','--skip_assembly', help='If the assembly has already \
|
|
4
|
155 been run, skip it with -a', action='store_true')
|
|
3
|
156
|
|
|
157 parser.add_argument('-c','--skip_conversion', help='If the fastq to fasta \
|
|
|
158 conversion has already been run, skip it with -c', \
|
|
4
|
159 action='store_true')
|
|
3
|
160
|
|
4
|
161 parser.add_argument('-primers', help='Filter \
|
|
3
|
162 pal_finder output to just those loci which have primers \
|
|
4
|
163 designed', action='store_true')
|
|
3
|
164
|
|
4
|
165 parser.add_argument('-occurrences', \
|
|
3
|
166 help='Filter pal_finder output to just loci with primers \
|
|
4
|
167 which only occur once in the dataset', action='store_true')
|
|
3
|
168
|
|
4
|
169 parser.add_argument('-rankmotifs', \
|
|
3
|
170 help='Filter pal_finder output to just loci which are a \
|
|
|
171 perfect repeat unit. Also, rank the loci by motif size \
|
|
4
|
172 (largest first)', action='store_true')
|
|
3
|
173
|
|
4
|
174 parser.add_argument('-v', '--get_version', help='Print the version number of \
|
|
|
175 this pal_filter script', action='store_true')
|
|
3
|
176
|
|
4
|
177 args = parser.parse_args()
|
|
3
|
178
|
|
4
|
179 if not args.assembly and not args.primers and not args.occurrences \
|
|
|
180 and not args.rankmotifs:
|
|
3
|
181 print "\nNo optional arguments supplied."
|
|
|
182 print "\nBy default this program does nothing."
|
|
|
183 print "\nNo files produced and no modifications made."
|
|
|
184 print "\nFinished.\n"
|
|
|
185 exit()
|
|
|
186 else:
|
|
|
187 print "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
|
|
|
188 print "Checking supplied filtering parameters:"
|
|
|
189 print "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
|
|
|
190 time.sleep(2)
|
|
4
|
191 if args.get_version:
|
|
|
192 print "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
|
|
|
193 print "pal_filter version is " + __version__ + " (03/03/2016)"
|
|
|
194 print "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"
|
|
|
195 if args.primers:
|
|
3
|
196 print "-primers flag supplied."
|
|
|
197 print "Filtering pal_finder output on the \"Primers found (1=y,0=n)\"" \
|
|
|
198 " column."
|
|
|
199 print "Only rows where primers have successfully been designed will"\
|
|
|
200 " pass the filter.\n"
|
|
|
201 time.sleep(2)
|
|
4
|
202 if args.occurrences:
|
|
3
|
203 print "-occurrences flag supplied."
|
|
4
|
204 print "Filtering pal_finder output on the \"Occurrences of Forward" \
|
|
|
205 " Primer in Reads\" and \"Occurrences of Reverse Primer" \
|
|
3
|
206 " in Reads\" columns."
|
|
|
207 print "Only rows where both primers occur only a single time in the"\
|
|
|
208 " reads pass the filter.\n"
|
|
|
209 time.sleep(2)
|
|
4
|
210 if args.rankmotifs:
|
|
3
|
211 print "-rankmotifs flag supplied."
|
|
|
212 print "Filtering pal_finder output on the \"Motifs(bases)\" column to" \
|
|
|
213 " just those with perfect repeats."
|
|
|
214 print "Only rows containing 'perfect' repeats will pass the filter."
|
|
|
215 print "Also, ranking output by size of motif (largest first).\n"
|
|
|
216 time.sleep(2)
|
|
4
|
217
|
|
3
|
218 # index the raw fastq files so that the sequences can be pulled out and
|
|
|
219 # added to the filtered output file
|
|
|
220 print "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
|
|
|
221 print "Indexing FastQ files....."
|
|
|
222 print "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
|
|
4
|
223 R1fastq_sequences_index = SeqIO.index(args.input1,'fastq')
|
|
|
224 R2fastq_sequences_index = SeqIO.index(args.input2,'fastq')
|
|
3
|
225 print "Indexing complete."
|
|
|
226
|
|
|
227 # create a set to hold the filtered output
|
|
|
228 wanted_lines = set()
|
|
|
229
|
|
|
230 # get lines from the pal_finder output which meet filter settings
|
|
|
231 # read the pal_finder output file into a csv reader
|
|
4
|
232 with open (args.pal_finder) as csvfile_infile:
|
|
3
|
233 csv_f = csv.reader(csvfile_infile, delimiter='\t')
|
|
4
|
234 header = csv_f.next()
|
|
|
235 header.extend(("R1_Sequence_ID", \
|
|
|
236 "R1_Sequence", \
|
|
|
237 "R2_Sequence_ID", \
|
|
|
238 "R2_Sequence" + "\n"))
|
|
|
239 with open( \
|
|
|
240 os.path.splitext(os.path.basename(args.pal_finder))[0] + \
|
|
|
241 ".filtered", 'w') as csvfile_outfile:
|
|
|
242 # write the header line for the output file
|
|
|
243 csvfile_outfile.write('\t'.join(header))
|
|
3
|
244 for row in csv_f:
|
|
4
|
245 # get the sequence ID
|
|
3
|
246 seq_ID = row[0]
|
|
4
|
247 # get the raw sequence reads and convert to a format that can
|
|
|
248 # go into a tsv file
|
|
3
|
249 R1_sequence = R1fastq_sequences_index[seq_ID].format("fasta").\
|
|
|
250 replace("\n","\t",1).replace("\n","")
|
|
|
251 R2_sequence = R2fastq_sequences_index[seq_ID].format("fasta").\
|
|
|
252 replace("\n","\t",1).replace("\n","")
|
|
|
253 seq_info = "\t" + R1_sequence + "\t" + R2_sequence + "\n"
|
|
4
|
254 # navigate through all different combinations of filter options
|
|
|
255 # if the primer filter is switched on
|
|
|
256 if args.primers:
|
|
|
257 # check the occurrences of primers field
|
|
3
|
258 if row[5] == "1":
|
|
4
|
259 # if filter occurrences of primers is switched on
|
|
|
260 if args.occurrences:
|
|
|
261 # check the occurrences of primers field
|
|
3
|
262 if (row[15] == "1" and row[16] == "1"):
|
|
4
|
263 # if rank by motif is switched on
|
|
|
264 if args.rankmotifs:
|
|
|
265 # check for perfect motifs
|
|
3
|
266 if row[1].count('(') == 1:
|
|
4
|
267 # all 3 filter switched on
|
|
|
268 # write line out to output
|
|
3
|
269 csvfile_outfile.write('\t'.join(row) + \
|
|
|
270 seq_info)
|
|
|
271 else:
|
|
|
272 csvfile_outfile.write('\t'.join(row) + seq_info)
|
|
4
|
273 elif args.rankmotifs:
|
|
3
|
274 if row[1].count('(') == 1:
|
|
|
275 csvfile_outfile.write('\t'.join(row) + seq_info)
|
|
|
276 else:
|
|
|
277 csvfile_outfile.write('\t'.join(row) + seq_info)
|
|
4
|
278 elif args.occurrences:
|
|
3
|
279 if (row[15] == "1" and row[16] == "1"):
|
|
4
|
280 if args.rankmotifs:
|
|
3
|
281 if row[1].count('(') == 1:
|
|
|
282 csvfile_outfile.write('\t'.join(row) + seq_info)
|
|
|
283 else:
|
|
|
284 csvfile_outfile.write('\t'.join(row) + seq_info)
|
|
4
|
285 elif args.rankmotifs:
|
|
3
|
286 if row[1].count('(') == 1:
|
|
|
287 csvfile_outfile.write('\t'.join(row) + seq_info)
|
|
|
288 else:
|
|
|
289 csvfile_outfile.write('\t'.join(row) + seq_info)
|
|
|
290
|
|
|
291 # if filter_rank_motifs is active, order the file by the size of the motif
|
|
4
|
292 if args.rankmotifs:
|
|
3
|
293 rank_motif = []
|
|
|
294 ranked_list = []
|
|
|
295 # read in the non-ordered file and add every entry to rank_motif list
|
|
4
|
296 with open( \
|
|
|
297 os.path.splitext(os.path.basename(args.pal_finder))[0] + \
|
|
|
298 ".filtered") as csvfile_ranksize:
|
|
3
|
299 csv_rank = csv.reader(csvfile_ranksize, delimiter='\t')
|
|
|
300 header = csv_rank.next()
|
|
|
301 for line in csv_rank:
|
|
|
302 rank_motif.append(line)
|
|
|
303
|
|
4
|
304 # open the file ready to write the ordered list
|
|
|
305 with open( \
|
|
|
306 os.path.splitext(os.path.basename(args.pal_finder))[0] + \
|
|
|
307 ".filtered", 'w') as rank_outfile:
|
|
3
|
308 rankwriter = csv.writer(rank_outfile, delimiter='\t', \
|
|
|
309 lineterminator='\n')
|
|
|
310 rankwriter.writerow(header)
|
|
|
311 count = 2
|
|
|
312 while count < 10:
|
|
|
313 for row in rank_motif:
|
|
4
|
314 # count size of motif
|
|
3
|
315 motif = re.search(r'[ATCG]*',row[1])
|
|
|
316 if motif:
|
|
|
317 the_motif = motif.group()
|
|
4
|
318 # rank it and write into ranked_list
|
|
3
|
319 if len(the_motif) == count:
|
|
|
320 ranked_list.insert(0, row)
|
|
|
321 count = count + 1
|
|
4
|
322 # write out the ordered list, into the .filtered file
|
|
3
|
323 for row in ranked_list:
|
|
|
324 rankwriter.writerow(row)
|
|
|
325
|
|
|
326 print "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
|
|
|
327 print "Checking assembly flags supplied:"
|
|
|
328 print "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
|
|
4
|
329 if not args.assembly:
|
|
3
|
330 print "Assembly flag not supplied. Not performing assembly QC.\n"
|
|
4
|
331 if args.assembly:
|
|
|
332 print "-assembly flag supplied: Perform PandaSeq assembly quality checks."
|
|
3
|
333 print "See Fox et al. (currently unpublished) for full details on the"\
|
|
|
334 " quality-check process.\n"
|
|
|
335 time.sleep(5)
|
|
|
336
|
|
|
337 # Get readID, F primers, R primers and motifs from filtered pal_finder output
|
|
|
338 seqIDs = []
|
|
|
339 motif = []
|
|
|
340 F_primers = []
|
|
|
341 R_primers = []
|
|
4
|
342 with open( \
|
|
|
343 os.path.splitext(os.path.basename(args.pal_finder))[0] + \
|
|
|
344 ".filtered") as input_csv:
|
|
3
|
345 pal_finder_csv = csv.reader(input_csv, delimiter='\t')
|
|
|
346 header = pal_finder_csv.next()
|
|
|
347 for row in pal_finder_csv:
|
|
|
348 seqIDs.append(row[0])
|
|
|
349 motif.append(row[1])
|
|
|
350 F_primers.append(row[7])
|
|
|
351 R_primers.append(row[9])
|
|
|
352
|
|
|
353 # Get a list of just the unique IDs we want
|
|
|
354 wanted = set()
|
|
|
355 for line in seqIDs:
|
|
|
356 wanted.add(line)
|
|
4
|
357 """
|
|
|
358 Assemble the paired end reads into overlapping contigs using PandaSeq
|
|
|
359 (can be skipped with the -a flag if assembly has already been run
|
|
|
360 and the appropriate files are in the same directory as the script,
|
|
|
361 and named "Assembly.fasta" and "Assembly_adapters_removed.fasta")
|
|
3
|
362
|
|
4
|
363 The first time you riun the script you MUST not enable the -a flag.t
|
|
|
364 but you are able to skip the assembly in subsequent analysis using the
|
|
|
365 -a flag.
|
|
|
366 """
|
|
|
367 if not args.skip_assembly:
|
|
|
368 pandaseq_command = 'pandaseq -A pear -f ' + args.input1 + ' -r ' + \
|
|
|
369 args.input2 + ' -o 25 -t 0.95 -w Assembly.fasta'
|
|
3
|
370 subprocess.call(pandaseq_command, shell=True)
|
|
4
|
371 strip_barcodes("Assembly.fasta", wanted)
|
|
3
|
372 print "\nPaired end reads been assembled into overlapping reads."
|
|
4
|
373 print "\nFor future analysis, you can skip this assembly step using" \
|
|
|
374 " the -a flag, provided that the assembly.fasta file" \
|
|
|
375 " is intact and in the same location."
|
|
3
|
376 else:
|
|
|
377 print "\n(Skipping the assembly step as you provided the -a flag)"
|
|
4
|
378 """
|
|
|
379 Fastq files need to be converted to fasta. The first time you run the script
|
|
|
380 you MUST not enable the -c flag, but you are able to skip the conversion
|
|
|
381 later using the -c flag. Make sure the fasta files are in the same location
|
|
|
382 and do not change the filenames
|
|
|
383 """
|
|
|
384 if not args.skip_conversion:
|
|
|
385 fastq_to_fasta(args.input1, wanted)
|
|
|
386 fastq_to_fasta(args.input2, wanted)
|
|
3
|
387 print "\nThe input fastq files have been converted to the fasta format."
|
|
4
|
388 print "\nFor any future analysis, you can skip this conversion step" \
|
|
|
389 " using the -c flag, provided that the fasta files" \
|
|
|
390 " are intact and in the same location."
|
|
3
|
391 else:
|
|
|
392 print "\n(Skipping the fastq -> fasta conversion as you provided the" \
|
|
|
393 " -c flag).\n"
|
|
|
394
|
|
4
|
395 # get the files and everything else needed
|
|
|
396 # Assembled fasta file
|
|
3
|
397 assembly_file = "Assembly_adapters_removed.fasta"
|
|
4
|
398 # filtered R1 reads
|
|
|
399 R1_fasta = os.path.splitext( \
|
|
|
400 os.path.basename(args.input1))[0] + "_filtered.fasta"
|
|
|
401 # filtered R2 reads
|
|
|
402 R2_fasta = os.path.splitext( \
|
|
|
403 os.path.basename(args.input2))[0] + "_filtered.fasta"
|
|
|
404 outputfilename = os.path.splitext(os.path.basename(args.input1))[0]
|
|
|
405 # parse the files with SeqIO
|
|
3
|
406 assembly_sequences = SeqIO.parse(assembly_file,'fasta')
|
|
|
407 R1fasta_sequences = SeqIO.parse(R1_fasta,'fasta')
|
|
|
408
|
|
4
|
409 # create some empty lists to hold the ID tags we are interested in
|
|
3
|
410 assembly_IDs = []
|
|
|
411 fasta_IDs = []
|
|
|
412
|
|
4
|
413 # populate the above lists with sequence IDs
|
|
3
|
414 for sequence in assembly_sequences:
|
|
|
415 assembly_IDs.append(sequence.id)
|
|
|
416 for sequence in R1fasta_sequences:
|
|
|
417 fasta_IDs.append(sequence.id)
|
|
|
418
|
|
4
|
419 # Index the assembly fasta file
|
|
3
|
420 assembly_sequences_index = SeqIO.index(assembly_file,'fasta')
|
|
|
421 R1fasta_sequences_index = SeqIO.index(R1_fasta,'fasta')
|
|
|
422 R2fasta_sequences_index = SeqIO.index(R2_fasta,'fasta')
|
|
|
423
|
|
4
|
424 # prepare the output file
|
|
|
425 with open ( \
|
|
|
426 outputfilename + "_pal_filter_assembly_output.txt", 'w') \
|
|
|
427 as outputfile:
|
|
|
428 # write the headers for the output file
|
|
|
429 output_header = ("readPairID", \
|
|
|
430 "Forward Primer",\
|
|
|
431 "F Primer Position in Assembled Read", \
|
|
|
432 "Reverse Primer", \
|
|
|
433 "R Primer Position in Assembled Read", \
|
|
|
434 "Motifs(bases)", \
|
|
|
435 "Assembled Read ID", \
|
|
|
436 "Assembled Read Sequence", \
|
|
|
437 "Raw Forward Read ID", \
|
|
|
438 "Raw Forward Read Sequence", \
|
|
|
439 "Raw Reverse Read ID", \
|
|
|
440 "Raw Reverse Read Sequence\n")
|
|
|
441 outputfile.write("\t".join(output_header))
|
|
3
|
442
|
|
4
|
443 # cycle through parameters from the pal_finder output
|
|
3
|
444 for x, y, z, a in zip(seqIDs, F_primers, R_primers, motif):
|
|
|
445 if str(x) in assembly_IDs:
|
|
|
446 # get the raw sequences ready to go into the output file
|
|
|
447 assembly_seq = (assembly_sequences_index.get_raw(x).decode())
|
|
4
|
448 # fasta entries need to be converted to single line so sit
|
|
|
449 # nicely in the output
|
|
|
450 assembly_output = assembly_seq.replace("\n","\t").strip('\t')
|
|
3
|
451 R1_fasta_seq = (R1fasta_sequences_index.get_raw(x).decode())
|
|
|
452 R1_output = R1_fasta_seq.replace("\n","\t",1).replace("\n","")
|
|
|
453 R2_fasta_seq = (R2fasta_sequences_index.get_raw(x).decode())
|
|
|
454 R2_output = R2_fasta_seq.replace("\n","\t",1).replace("\n","")
|
|
|
455 assembly_no_id = '\n'.join(assembly_seq.split('\n')[1:])
|
|
|
456
|
|
4
|
457 # check that both primer sequences can be seen in the
|
|
|
458 # assembled contig
|
|
3
|
459 if y or ReverseComplement1(y) in assembly_no_id and z or \
|
|
|
460 ReverseComplement1(z) in assembly_no_id:
|
|
|
461 if y in assembly_no_id:
|
|
4
|
462 # get the positions of the primers in the assembly
|
|
|
463 # (can be used to predict fragment length)
|
|
3
|
464 F_position = assembly_no_id.index(y)+len(y)+1
|
|
|
465 if ReverseComplement1(y) in assembly_no_id:
|
|
4
|
466 F_position = assembly_no_id.index( \
|
|
|
467 ReverseComplement1(y))+len(ReverseComplement1(y))+1
|
|
3
|
468 if z in assembly_no_id:
|
|
|
469 R_position = assembly_no_id.index(z)+1
|
|
|
470 if ReverseComplement1(z) in assembly_no_id:
|
|
4
|
471 R_position = assembly_no_id.index( \
|
|
|
472 ReverseComplement1(z))+1
|
|
|
473 output = (str(x),
|
|
|
474 str(y),
|
|
|
475 str(F_position),
|
|
|
476 str(z),
|
|
|
477 str(R_position),
|
|
|
478 str(a),
|
|
|
479 str(assembly_output),
|
|
|
480 str(R1_output),
|
|
|
481 str(R2_output + "\n"))
|
|
|
482 outputfile.write("\t".join(output))
|
|
3
|
483 print "\nPANDAseq quality check complete."
|
|
|
484 print "Results from PANDAseq quality check (and filtering, if any" \
|
|
|
485 " any filters enabled) written to output file" \
|
|
4
|
486 " ending \"_pal_filter_assembly_output.txt\".\n"
|
|
3
|
487 print "Filtering of pal_finder results complete."
|
|
|
488 print "Filtered results written to output file ending \".filtered\"."
|
|
|
489 print "\nFinished\n"
|
|
|
490 else:
|
|
4
|
491 if args.skip_assembly or args.skip_conversion:
|
|
3
|
492 print "\nERROR: You cannot supply the -a flag or the -c flag without \
|
|
|
493 also supplying the -assembly flag.\n"
|
|
|
494 print "\nProgram Finished\n"
|
