comparison filter.py @ 0:da1a6f33b504 draft default tip

Imported from capsule None
author devteam
date Mon, 27 Jan 2014 09:29:09 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:da1a6f33b504
1 #!/usr/bin/python
2
3 import os.path
4 import sys
5 import optparse
6
7 import vcfClass
8 from vcfClass import *
9
10 import tools
11 from tools import *
12
13 if __name__ == "__main__":
14 main()
15
16 def filterFail(text, file):
17 print >> sys.stderr, text
18 if file != None: os.remove(file)
19 exit(1)
20
21 def main():
22
23 # Parse the command line options
24 usage = "Usage: vcfPytools.py filter [options]"
25 parser = optparse.OptionParser(usage = usage)
26 parser.add_option("-i", "--in",
27 action="store", type="string",
28 dest="vcfFile", help="input vcf file")
29 parser.add_option("-o", "--out",
30 action="store", type="string",
31 dest="output", help="output vcf file")
32 parser.add_option("-q", "--quality",
33 action="store", type="int",
34 dest="quality", help="filter out SNPs with qualities lower than selected value")
35 parser.add_option("-n", "--info",
36 action="append", type="string", nargs=3,
37 dest="infoFilters", help="filter based on entries in the info string")
38 parser.add_option("-r", "--remove-genotypes",
39 action="store_true", default=False,
40 dest="removeGeno", help="remove the genotype strings from the vcf file")
41 parser.add_option("-m", "--mark-as-pass",
42 action="store_true", default=False,
43 dest="markPass", help="Mark all records as having passed filters")
44
45 (options, args) = parser.parse_args()
46
47 # Check that a single vcf file is given.
48 if options.vcfFile == None:
49 parser.print_help()
50 print >> sys.stderr, "\nInput vcf file (-i, --input) is required for vcf filtering."
51 exit(1)
52
53 # The --mark-as-pass option can only be used if no actual filters
54 # have been specified.
55 if options.markPass and options.infoFilters:
56 print >> sys.stderr, "--mark-as-pass cannot be used in conjunction with filters."
57 exit(1)
58
59 # Set the output file to stdout if no output file was specified.
60 outputFile, writeOut = setOutput(options.output) # tools.py
61
62 v = vcf() # Define vcf object.
63
64 # Open the vcf file.
65 v.openVcf(options.vcfFile)
66
67 # Read in the header information.
68 v.parseHeader(options.vcfFile, writeOut)
69 taskDescriptor = "##vcfPytools="
70 if options.infoFilters:
71 taskDescriptor += "filtered using the following filters: "
72 for filter, value, logic in options.infoFilters: taskDescriptor += str(filter) + str(value) + ","
73 taskDescriptor = taskDescriptor.rstrip(",")
74 if options.markPass: taskDescriptor += "marked all records as PASS"
75
76 writeHeader(outputFile, v, options.removeGeno, taskDescriptor)
77
78 # Check that specified filters from the info field are either integers or floats.
79 if options.infoFilters:
80 v.processInfo = True # Process the info string
81 filters = {}
82 filterValues = {}
83 filterLogic = {}
84 for filter, value, logic in options.infoFilters:
85 filterName = str(filter) + str(value)
86 if "-" in filter or "-" in value or "-" in logic:
87 print >> sys.stderr, "\n--info (-n) requires three arguments, for example:"
88 print >> sys.stderr, "\t--info DP 5 lt: filter records with DP less than (lt) 5.\n"
89 print >> sys.stderr, "allowed logic arguments:\n\tgt: greater than\n\tlt: less than."
90 print >> sys.stderr, "\nError in:", filter
91 exit(1)
92 if logic != "gt" and logic != "lt":
93 print >> sys.stderr, "\nfilter logic not recognised."
94 print >> sys.stderr, "allowed logic arguments:\n\tgt: greater than\n\tlt: less than."
95 print >> sys.stderr, "\nError in:", filter
96 exit(1)
97 if v.infoHeaderTags.has_key(filter):
98 if v.infoHeaderTags[filter][1].lower() == "integer":
99 try:
100 filters[filterName] = filter
101 filterValues[filterName] = int(value)
102 filterLogic[filterName] = logic
103 #filterLogic[filterName] = logic
104 except ValueError:
105 text = "Filter " + filter + " requires an integer entry, not " + str(type(value))
106 filterFail(text, options.output)
107
108 if v.infoHeaderTags[filter][1].lower() == "float":
109 try:
110 filters[filterName] = filter
111 filterValues[filterName] = float(value)
112 filterLogic[filterName] = logic
113 #filters[filterName] = float(value)
114 #filterLogic[filterName] = logic
115 except ValueError:
116 text = "Filter " + filter + " requires an float entry, not " + str(type(value))
117 filterFail(text, options.output)
118
119 else:
120 text = "Filter " + filter + " has no explanation in the header. Unknown type for the entry."
121 filterFail(text, options.output)
122
123 # Parse the vcf file and check if any of the filters are failed. If
124 # so, build up a string of failed filters.
125 while v.getRecord():
126 filterString = ""
127
128 # Mark the record as "PASS" if --mark-as-pass was applied.
129 if options.markPass: v.filters = "PASS"
130
131 # Check for quality filtering.
132 if options.quality != None:
133 if v.quality < options.quality:
134 filterString = filterString + ";" + "Q" + str(options.quality) if filterString != "" else "Q" + str(options.quality)
135
136 # Check for filtering on info string filters.
137 if options.infoFilters:
138 for filterName, filter in filters.iteritems():
139 value = filterValues[filterName]
140 logic = filterLogic[filterName]
141 if v.infoTags.has_key(filter):
142 if type(value) == int:
143 if logic == "lt" and int(v.infoTags[filter]) < value:
144 filterString = filterString + ";" + filter + str(value) if filterString != "" else filter + str(value)
145 if logic == "gt" and int(v.infoTags[filter]) > value:
146 filterString = filterString + ";" + filter + str(value) if filterString != "" else filter + str(value)
147 elif type(value) == float:
148 if logic == "lt" and float(v.infoTags[filter]) < value:
149 filterString = filterString + ";" + filter + str(value) if filterString != "" else filter + str(value)
150 if logic == "gt" and float(v.infoTags[filter]) > value:
151 filterString = filterString + ";" + filter + str(value) if filterString != "" else filter + str(value)
152
153 filterString = "PASS" if filterString == "" else filterString
154 v.filters = filterString
155 record = v.buildRecord(options.removeGeno)
156 outputFile.write(record)
157
158 # Close the vcf files.
159 v.closeVcf(options.vcfFile)
160
161 # Terminate the program.
162 return 0