|
0
|
1 import argparse
|
|
|
2 import csv
|
|
|
3 import os
|
|
|
4 import sys
|
|
|
5
|
|
|
6 from fasta_extract_utils import Fasta
|
|
|
7
|
|
|
8
|
|
|
9 def reverse_complement(bases):
|
|
|
10 complements = {'A': 'T', 'T': 'A', 'C': 'G', 'G': 'C'}
|
|
|
11 return ''.join(complements[b.upper()] for b in reversed(bases))
|
|
|
12
|
|
|
13
|
|
|
14 def get_output_path(hid, subtract_from_start, add_to_end, extend_existing, consider_strand, orphan=False):
|
|
|
15 attrs = 'u%dd%d' % (subtract_from_start, add_to_end)
|
|
|
16 if extend_existing:
|
|
|
17 attrs += 'x'
|
|
|
18 if consider_strand:
|
|
|
19 attrs += '_s'
|
|
|
20 if orphan:
|
|
|
21 attrs += '_orphan'
|
|
|
22 format = 'gff'
|
|
7
|
23 output_dir = 'output_orphan_dir'
|
|
0
|
24 else:
|
|
|
25 format = 'fasta'
|
|
7
|
26 output_dir = 'output_dir'
|
|
|
27 return os.path.join(output_dir, '%s_on_data_%d.%s' % (attrs, hid, format))
|
|
0
|
28
|
|
|
29
|
|
|
30 def stop_err(msg):
|
|
|
31 sys.stderr.write(msg)
|
|
|
32 sys.exit(1)
|
|
|
33
|
|
|
34
|
|
|
35 parser = argparse.ArgumentParser()
|
|
1
|
36 parser.add_argument('--input', dest='inputs', action='append', nargs=2, help="Input datasets")
|
|
0
|
37 parser.add_argument('--genome_file', dest='genome_file', help='Reference genome fasta index file.')
|
|
|
38 parser.add_argument('--subtract_from_start', dest='subtract_from_start', type=int, help='Distance to subtract from start.')
|
|
|
39 parser.add_argument('--add_to_end', dest='add_to_end', type=int, help='Distance to add to end.')
|
|
|
40 parser.add_argument('--extend_existing', dest='extend_existing', help='Extend existing start/end rather or from computed midpoint.')
|
|
|
41 parser.add_argument('--strand', dest='strand', help='Consider strandedness: reverse complement extracted sequence on reverse strand.')
|
|
|
42 args = parser.parse_args()
|
|
|
43
|
|
|
44 fasta = Fasta(args.genome_file)
|
|
|
45
|
|
|
46 for (input_filename, hid) in args.inputs:
|
|
|
47 extend_existing = args.extend_existing == 'existing'
|
|
|
48 consider_strand = args.strand == 'yes'
|
|
|
49 reader = csv.reader(open(input_filename, 'rU'), delimiter='\t')
|
|
3
|
50 fasta_output_path = get_output_path(hid,
|
|
0
|
51 args.subtract_from_start,
|
|
|
52 args.add_to_end,
|
|
|
53 extend_existing,
|
|
|
54 consider_strand)
|
|
|
55 output = open(fasta_output_path, 'wb')
|
|
3
|
56 gff_output_path = get_output_path(hid,
|
|
0
|
57 args.subtract_from_start,
|
|
|
58 args.add_to_end,
|
|
|
59 extend_existing,
|
|
|
60 consider_strand,
|
|
|
61 orphan=True)
|
|
|
62 orphan_writer = csv.writer(open(gff_output_path, 'wb'), delimiter='\t')
|
|
|
63 for row in reader:
|
|
|
64 if len(row) != 9 or row[0].startswith('#'):
|
|
|
65 continue
|
|
|
66 try:
|
|
|
67 cname = row[0]
|
|
|
68 start = int(row[3])
|
|
|
69 end = int(row[4])
|
|
|
70 strand = row[6]
|
|
|
71 if extend_existing:
|
|
|
72 start -= args.subtract_from_start
|
|
|
73 end += args.add_to_end
|
|
|
74 else:
|
|
|
75 midpoint = (start + end) // 2
|
|
|
76 start = midpoint - args.subtract_from_start
|
|
|
77 end = midpoint + args.add_to_end
|
|
|
78 if 1 <= start and end <= len(fasta[cname]):
|
|
|
79 output.write('>%s:%s-%s_%s\n' % (cname, start, end, strand))
|
|
|
80 bases = fasta[cname][start-1:end]
|
|
|
81 if consider_strand and strand == '-':
|
|
|
82 bases = reverse_complement(bases)
|
|
|
83 output.write('%s\n' % bases)
|
|
|
84 else:
|
|
|
85 orphan_writer.writerow(row)
|
|
|
86 except Exception, e:
|
|
|
87 stop_err(str(e))
|
|
7
|
88 finally:
|
|
|
89 output.clode()
|
