misc/varfilter.py

   1 #!/software/bin/python
   2
   3 # Author: lh3, converted to python and modified to add -C option by Aylwyn Scally
   4 #
   5 # About:
   6 #   varfilter.py is a port of Heng's samtools.pl varFilter script into
   7 #   python, with an additional -C INT option. This option sets a minimum
   8 #   consensus score, above which the script will output a pileup line
   9 #   wherever it _could have_ called a variant, even if none is actually
  10 #   called (i.e. hom-ref positions). This is important if you want to
  11 #   subsequently merge the calls with those for another individual to get a
  12 #   synoptic view of calls at each site. Without this option, and in all
  13 #   other respects, it behaves like samtools.pl varFilter.
  14 #
  15 #   Aylwyn Scally as6@sanger.ac.uk
  16
  17
  18 # Filtration code:
  19 #
  20 # C low CNS quality (hom-ref only)
  21 # d low depth
  22 # D high depth
  23 # W too many SNPs in a window (SNP only)
  24 # G close to a high-quality indel (SNP only)
  25 # Q low RMS mapping quality (SNP only)
  26 # g close to another indel with higher quality (indel only)
  27 # s low SNP quality (SNP only)
  28 # i low indel quality (indel only)
  29
  30
  31 import sys
  32 import getopt
  33
  34 def usage():
  35         print '''usage: varfilter.py [options] [cns-pileup]
  36
  37 Options: -Q INT minimum RMS mapping quality for SNPs
  38                  -q INT minimum RMS mapping quality for gaps
  39                  -d INT minimum read depth
  40                  -D INT maximum read depth
  41                  -S INT minimum SNP quality
  42                  -i INT minimum indel quality
  43                  -C INT minimum consensus quality for hom-ref sites
  44
  45                  -G INT min indel score for nearby SNP filtering
  46                  -w INT SNP within INT bp around a gap to be filtered
  47
  48                  -W INT window size for filtering dense SNPs
  49                  -N INT max number of SNPs in a window
  50
  51                  -l INT window size for filtering adjacent gaps
  52
  53                  -p print filtered variants'''
  54
  55 def varFilter_aux(first, is_print):
  56         try:
  57                 if first[1] == 0:
  58                         sys.stdout.write("\t".join(first[4:]) + "\n")
  59                 elif is_print:
  60                         sys.stderr.write("\t".join(["UQdDWGgsiCX"[first[1]]] + first[4:]) + "\n")
  61         except IOError:
  62                 sys.exit()
  63
  64 mindepth = 3
  65 maxdepth = 100
  66 gapgapwin = 30
  67 minsnpmapq = 25
  68 mingapmapq = 10
  69 minindelscore = 25
  70 scorefactor = 100
  71 snpgapwin = 10
  72 densesnpwin = 10
  73 densesnps = 2
  74 printfilt = False
  75 minsnpq = 0
  76 minindelq = 0
  77 mincnsq = 0
  78
  79 try:
  80         options, args = getopt.gnu_getopt(sys.argv[1:], 'pq:d:D:l:Q:w:W:N:G:S:i:C:', [])
  81 except getopt.GetoptError:
  82         usage()
  83         sys.exit(2)
  84 for (oflag, oarg) in options:
  85         if oflag == '-d': mindepth = int(oarg)
  86         if oflag == '-D': maxdepth = int(oarg)
  87         if oflag == '-l': gapgapwin = int(oarg)
  88         if oflag == '-Q': minsnpmapq = int(oarg)
  89         if oflag == '-q': mingapmapq = int(oarg)
  90         if oflag == '-G': minindelscore = int(oarg)
  91         if oflag == '-s': scorefactor = int(oarg)
  92         if oflag == '-w': snpgapwin = int(oarg)
  93         if oflag == '-W': densesnpwin = int(oarg)
  94         if oflag == '-C': mincnsq = int(oarg)
  95         if oflag == '-N': densesnps = int(oarg)
  96         if oflag == '-p': printfilt = True
  97         if oflag == '-S': minsnpq = int(oarg)
  98         if oflag == '-i': minindelq = int(oarg)
  99
 100 if len(args) < 1:
 101         inp = sys.stdin
 102 else:
 103         inp = open(args[0])
 104
 105 # calculate the window size
 106 max_dist = max(gapgapwin, snpgapwin, densesnpwin)
 107
 108 staging = []
 109 for t in (line.strip().split() for line in inp):
 110         (flt, score) = (0, -1)
 111         # non-var sites
 112         if t[3] == '*/*':
 113                 continue
 114         is_snp = t[2].upper() != t[3].upper()
 115         if not (is_snp or mincnsq):
 116                 continue
 117         # clear the out-of-range elements
 118         while staging:
 119                 # Still on the same chromosome and the first element's window still affects this position?
 120                 if staging[0][4] == t[0] and int(staging[0][5]) + staging[0][2] + max_dist >= int(t[1]):
 121                         break
 122                 varFilter_aux(staging.pop(0), printfilt)
 123
 124         # first a simple filter
 125         if int(t[7]) < mindepth:
 126                 flt = 2
 127         elif int(t[7]) > maxdepth:
 128                 flt = 3
 129         if t[2] == '*': # an indel
 130                 if minindelq and minindelq > int(t[5]):
 131                         flt = 8
 132         elif is_snp:
 133                 if minsnpq and minsnpq> int(t[5]):
 134                         flt = 7
 135         else:
 136                 if mincnsq and mincnsq > int(t[4]):
 137                         flt = 9
 138
 139         # site dependent filters
 140         dlen = 0
 141         if flt == 0:
 142                 if t[2] == '*': # an indel
 143                         # If deletion, remember the length of the deletion
 144                         (a,b) = t[3].split('/')
 145                         alen = len(a) - 1
 146                         blen = len(b) - 1
 147                         if alen>blen:
 148                                 if a[0] == '-': dlen=alen
 149                         elif b[0] == '-': dlen=blen
 150
 151                         if int(t[6]) < mingapmapq:
 152                                 flt = 1
 153                         # filtering SNPs
 154                         if int(t[5]) >= minindelscore:
 155                                 for x in (y for y in staging if y[3]):
 156                                         # Is it a SNP and is it outside the SNP filter window?
 157                                         if x[0] >= 0 or int(x[5]) + x[2] + snpgapwin < int(t[1]):
 158                                                 continue
 159                                         if x[1] == 0:
 160                                                 x[1] = 5
 161
 162                         # calculate the filtering score (different from indel quality)
 163                         score = int(t[5])
 164                         if t[8] != '*':
 165                                 score += scorefactor * int(t[10])
 166                         if t[9] != '*':
 167                                 score += scorefactor * int(t[11])
 168                         # check the staging list for indel filtering
 169                         for x in (y for y in staging if y[3]):
 170                           # Is it a SNP and is it outside the gap filter window
 171                                 if x[0] < 0 or int(x[5]) + x[2] + gapgapwin < int(t[1]):
 172                                         continue
 173                                 if x[0] < score:
 174                                         x[1] = 6
 175                                 else:
 176                                         flt = 6
 177                                         break
 178                 else: # a SNP or hom-ref
 179                         if int(t[6]) < minsnpmapq:
 180                                 flt = 1
 181                         # check adjacent SNPs
 182                         k = 1
 183                         for x in (y for y in staging if y[3]):
 184                                 if x[0] < 0 and int(x[5]) + x[2] + densesnpwin >= int(t[1]) and (x[1] == 0 or x[1] == 4 or x[1] == 5):
 185                                         k += 1
 186
 187                         # filtering is necessary
 188                         if k > densesnps:
 189                                 flt = 4
 190                                 for x in (y for y in staging if y[3]):
 191                                         if x[0] < 0 and int(x[5]) + x[2] + densesnpwin >= int(t[1]) and x[1] == 0:
 192                                                 x[1] = 4
 193                         else: # then check gap filter
 194                                 for x in (y for y in staging if y[3]):
 195                                         if x[0] < 0 or int(x[5]) + x[2] + snpgapwin < int(t[1]):
 196                                                 continue
 197                                         if x[0] >= minindelscore:
 198                                                 flt = 5
 199                                                 break
 200
 201         staging.append([score, flt, dlen, is_snp] + t)
 202
 203 # output the last few elements in the staging list
 204 while staging:
 205         varFilter_aux(staging.pop(0), printfilt)