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