|
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'
|
|
|
23 else:
|
|
|
24 format = 'fasta'
|
|
5
|
25 if orphan:
|
|
|
26 return os.path.join('output_orphan_dir', '%s_on_data_%d.%s' % (attrs, hid, format))
|
|
|
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:
|
|
4
|
47 hid = int(hid)
|
|
0
|
48 extend_existing = args.extend_existing == 'existing'
|
|
|
49 consider_strand = args.strand == 'yes'
|
|
|
50 reader = csv.reader(open(input_filename, 'rU'), delimiter='\t')
|
|
3
|
51 fasta_output_path = get_output_path(hid,
|
|
0
|
52 args.subtract_from_start,
|
|
|
53 args.add_to_end,
|
|
|
54 extend_existing,
|
|
|
55 consider_strand)
|
|
|
56 output = open(fasta_output_path, 'wb')
|
|
3
|
57 gff_output_path = get_output_path(hid,
|
|
0
|
58 args.subtract_from_start,
|
|
|
59 args.add_to_end,
|
|
|
60 extend_existing,
|
|
|
61 consider_strand,
|
|
|
62 orphan=True)
|
|
|
63 orphan_writer = csv.writer(open(gff_output_path, 'wb'), delimiter='\t')
|
|
|
64 for row in reader:
|
|
|
65 if len(row) != 9 or row[0].startswith('#'):
|
|
|
66 continue
|
|
|
67 try:
|
|
|
68 cname = row[0]
|
|
|
69 start = int(row[3])
|
|
|
70 end = int(row[4])
|
|
|
71 strand = row[6]
|
|
|
72 if extend_existing:
|
|
|
73 start -= args.subtract_from_start
|
|
|
74 end += args.add_to_end
|
|
|
75 else:
|
|
|
76 midpoint = (start + end) // 2
|
|
|
77 start = midpoint - args.subtract_from_start
|
|
|
78 end = midpoint + args.add_to_end
|
|
|
79 if 1 <= start and end <= len(fasta[cname]):
|
|
|
80 output.write('>%s:%s-%s_%s\n' % (cname, start, end, strand))
|
|
|
81 bases = fasta[cname][start-1:end]
|
|
|
82 if consider_strand and strand == '-':
|
|
|
83 bases = reverse_complement(bases)
|
|
|
84 output.write('%s\n' % bases)
|
|
|
85 else:
|
|
|
86 orphan_writer.writerow(row)
|
|
|
87 except Exception, e:
|
|
|
88 stop_err(str(e))
|
