|
0
|
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
|
