Mercurial > repos > devteam > vcf_filter

comparison filter.py @ 0:da1a6f33b504 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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