Mercurial > repos > iuc > extract_genomic_dna

comparison extract_genomic_dna.py @ 0:8dd8e89c0603 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/extract_genomic_dna commit b'67cff25a50ba173b0468819204d0999496f68ea9'
author iuc
date Tue, 19 Jan 2016 09:34:23 -0500
parents
children 702970e4a134
comparison
equal deleted inserted replaced
-1:000000000000 0:8dd8e89c0603
1 #!/usr/bin/env python
2 import argparse
3 import os
4
5 import extract_genomic_dna_utils as egdu
6 import bx.seq.nib
7 import bx.seq.twobit
8 from bx.intervals.io import Header, Comment
9
10
11 parser = argparse.ArgumentParser()
12 parser.add_argument('--input_format', dest='input_format', help="Input dataset format")
13 parser.add_argument('--input', dest='input', help="Input dataset")
14 parser.add_argument('--genome', dest='genome', help="Input dataset genome build")
15 parser.add_argument('--interpret_features', dest='interpret_features', default=None, help="Interpret features if input format is gff")
16 parser.add_argument('--columns', dest='columns', help="Columns to use in input file")
17 parser.add_argument('--reference_genome_source', dest='reference_genome_source', help="Source of reference genome file")
18 parser.add_argument('--reference_genome', dest='reference_genome', help="Reference genome file")
19 parser.add_argument('--output_format', dest='output_format', help="Output format")
20 parser.add_argument('--output', dest='output', help="Output dataset")
21 args = parser.parse_args()
22
23 input_is_gff = args.input_format == 'gff'
24 interpret_features = input_is_gff and args.interpret_features == "yes"
25 if len(args.columns.split(',')) == 5:
26 # Bed file.
27 chrom_col, start_col, end_col, strand_col, name_col = egdu.parse_cols_arg(args.columns)
28 else:
29 # Gff file.
30 chrom_col, start_col, end_col, strand_col = egdu.parse_cols_arg(args.columns)
31 name_col = False
32
33 if args.reference_genome_source == "history":
34 seq_path = egdu.convert_to_twobit(args.reference_genome)
35 else:
36 seq_path = args.reference_genome
37 seq_dir = os.path.split(seq_path)[0]
38
39 includes_strand_col = strand_col >= 0
40 strand = None
41 nibs = {}
42 skipped_lines = 0
43 first_invalid_line = 0
44 invalid_lines = []
45 warnings = []
46 warning = ''
47 twobitfile = None
48 line_count = 1
49 file_iterator = open(args.input)
50 if interpret_features:
51 file_iterator = egdu.GFFReaderWrapper(file_iterator, fix_strand=False)
52 out = open(args.output, 'wt')
53
54 for feature in file_iterator:
55 # Ignore comments, headers.
56 if isinstance(feature, (Header, Comment)):
57 line_count += 1
58 continue
59 name = ""
60 if interpret_features:
61 # Processing features.
62 egdu.convert_gff_coords_to_bed(feature)
63 chrom = feature.chrom
64 start = feature.start
65 end = feature.end
66 strand = feature.strand
67 else:
68 # Processing lines, either interval or GFF format.
69 line = feature.rstrip('\r\n')
70 if line and not line.startswith("#"):
71 fields = line.split('\t')
72 try:
73 chrom = fields[chrom_col]
74 start = int(fields[start_col])
75 end = int(fields[end_col])
76 if name_col:
77 name = fields[name_col]
78 if input_is_gff:
79 start, end = egdu.convert_gff_coords_to_bed([start, end])
80 if includes_strand_col:
81 strand = fields[strand_col]
82 except:
83 warning = "Invalid chrom, start or end column values. "
84 warnings.append(warning)
85 if not invalid_lines:
86 invalid_lines = egdu.get_lines(feature)
87 first_invalid_line = line_count
88 skipped_lines += len(invalid_lines)
89 continue
90 if start > end:
91 warning = "Invalid interval, start '%d' > end '%d'. " % (start, end)
92 warnings.append(warning)
93 if not invalid_lines:
94 invalid_lines = egdu.get_lines(feature)
95 first_invalid_line = line_count
96 skipped_lines += len(invalid_lines)
97 continue
98 if strand not in ['+', '-']:
99 strand = '+'
100 sequence = ''
101 else:
102 continue
103 # Open sequence file and get sequence for feature/interval.
104 if os.path.exists("%s/%s.nib" % (seq_dir, chrom)):
105 if chrom in nibs:
106 nib = nibs[chrom]
107 else:
108 nibs[chrom] = nib = bx.seq.nib.NibFile(open("%s/%s.nib" % (seq_path, chrom)))
109 try:
110 sequence = nib.get(start, end - start)
111 except Exception, e:
112 warning = "Unable to fetch the sequence from '%d' to '%d' for build '%s'. " % (start, end - start, args.genome)
113 warnings.append(warning)
114 if not invalid_lines:
115 invalid_lines = egdu.get_lines(feature)
116 first_invalid_line = line_count
117 skipped_lines += len(invalid_lines)
118 continue
119 elif os.path.isfile(seq_path):
120 if not(twobitfile):
121 twobitfile = bx.seq.twobit.TwoBitFile(open(seq_path))
122 try:
123 if interpret_features:
124 # Create sequence from intervals within a feature.
125 sequence = ''
126 for interval in feature.intervals:
127 sequence += twobitfile[interval.chrom][interval.start:interval.end]
128 else:
129 sequence = twobitfile[chrom][start:end]
130 except:
131 warning = "Unable to fetch the sequence from '%d' to '%d' for chrom '%s'. " % (start, end - start, chrom)
132 warnings.append(warning)
133 if not invalid_lines:
134 invalid_lines = egdu.get_lines(feature)
135 first_invalid_line = line_count
136 skipped_lines += len(invalid_lines)
137 continue
138 else:
139 warning = "Chromosome by name '%s' was not found for build '%s'. " % (chrom, args.genome)
140 warnings.append(warning)
141 if not invalid_lines:
142 invalid_lines = egdu.get_lines(feature)
143 first_invalid_line = line_count
144 skipped_lines += len(invalid_lines)
145 continue
146 if sequence == '':
147 warning = "Chrom: '%s', start: '%d', end: '%d' is either invalid or not present in build '%s'. " % (chrom, start, end, args.genome)
148 warnings.append(warning)
149 if not invalid_lines:
150 invalid_lines = egdu.get_lines(feature)
151 first_invalid_line = line_count
152 skipped_lines += len(invalid_lines)
153 continue
154 if includes_strand_col and strand == "-":
155 sequence = egdu.reverse_complement(sequence)
156 if args.output_format == "fasta":
157 l = len(sequence)
158 c = 0
159 if input_is_gff:
160 start, end = egdu.convert_bed_coords_to_gff([start, end])
161 fields = [args.genome, str(chrom), str(start), str(end), strand]
162 meta_data = "_".join(fields)
163 if name.strip():
164 out.write(">%s %s\n" % (meta_data, name))
165 else:
166 out.write(">%s\n" % meta_data)
167 while c < l:
168 b = min(c + 50, l)
169 out.write("%s\n" % str(sequence[c:b]))
170 c = b
171 else:
172 # output_format == "interval".
173 if interpret_features:
174 meta_data = "\t".join([feature.chrom,
175 "galaxy_extract_genomic_dna",
176 "interval",
177 str(feature.start),
178 str(feature.end),
179 feature.score,
180 feature.strand,
181 ".",
182 egdu.gff_attributes_to_str(feature.attributes, "GTF")])
183 else:
184 # Where is fields being set here?
185 meta_data = "\t".join(fields)
186 if input_is_gff:
187 format_str = "%s seq \"%s\";\n"
188 else:
189 format_str = "%s\t%s\n"
190 out.write(format_str % (meta_data, str(sequence)))
191 # Update line count.
192 if isinstance(feature, egdu.GFFFeature):
193 line_count += len(feature.intervals)
194 else:
195 line_count += 1
196 out.close()
197
198 if warnings:
199 warn_msg = "%d warnings, 1st is: " % len(warnings)
200 warn_msg += warnings[0]
201 print warn_msg
202 if skipped_lines:
203 # Error message includes up to the first 10 skipped lines.
204 print 'Skipped %d invalid lines, 1st is #%d, "%s"' % (skipped_lines, first_invalid_line, '\n'.join(invalid_lines[:10]))
205
206 if args.reference_genome_source == "history":
207 os.remove(seq_path)