/tools/vcf_tools/filter.py

  1#!/usr/bin/python
  2
  3import os.path
  4import sys
  5import optparse
  6
  7import vcfClass
  8from vcfClass import *
  9
 10import tools
 11from tools import *
 12
 13if __name__ == "__main__":
 14  main()
 15
 16def filterFail(text, file):
 17  print >> sys.stderr, text
 18  if file != None: os.remove(file)
 19  exit(1)
 20
 21def 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