Mercurial > repos > rliterman > csp2

annotate CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/lib/python3.8/site-packages/pysam/libcvcf.pyx @ 69:33d812a61356

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 2e9511a184a1ca667c7be0c6321a36dc4e3d116d
author jpayne
date Tue, 18 Mar 2025 17:55:14 -0400
parents
children
rev   line source
jpayne@69 1 # cython: embedsignature=True
jpayne@69 2 #
jpayne@69 3 # Code to read, write and edit VCF files
jpayne@69 4 #
jpayne@69 5 # VCF lines are encoded as a dictionary with these keys (note: all lowercase):
jpayne@69 6 # 'chrom': string
jpayne@69 7 # 'pos': integer
jpayne@69 8 # 'id': string
jpayne@69 9 # 'ref': string
jpayne@69 10 # 'alt': list of strings
jpayne@69 11 # 'qual': integer
jpayne@69 12 # 'filter': None (missing value), or list of keys (strings); empty list parsed as ["PASS"]
jpayne@69 13 # 'info': dictionary of values (see below)
jpayne@69 14 # 'format': list of keys (strings)
jpayne@69 15 # sample keys: dictionary of values (see below)
jpayne@69 16 #
jpayne@69 17 # The sample keys are accessible through vcf.getsamples()
jpayne@69 18 #
jpayne@69 19 # A dictionary of values contains value keys (defined in ##INFO or
jpayne@69 20 # ##FORMAT lines) which map to a list, containing integers, floats,
jpayne@69 21 # strings, or characters. Missing values are replaced by a particular
jpayne@69 22 # value, often -1 or .
jpayne@69 23 #
jpayne@69 24 # Genotypes are not stored as a string, but as a list of 1 or 3
jpayne@69 25 # elements (for haploid and diploid samples), the first (and last) the
jpayne@69 26 # integer representing an allele, and the second the separation
jpayne@69 27 # character. Note that there is just one genotype per sample, but for
jpayne@69 28 # consistency the single element is stored in a list.
jpayne@69 29 #
jpayne@69 30 # Header lines other than ##INFO, ##FORMAT and ##FILTER are stored as
jpayne@69 31 # (key, value) pairs and are accessible through getheader()
jpayne@69 32 #
jpayne@69 33 # The VCF class can be instantiated with a 'regions' variable
jpayne@69 34 # consisting of tuples (chrom,start,end) encoding 0-based half-open
jpayne@69 35 # segments. Only variants with a position inside the segment will be
jpayne@69 36 # parsed. A regions parser is available under parse_regions.
jpayne@69 37 #
jpayne@69 38 # When instantiated, a reference can be passed to the VCF class. This
jpayne@69 39 # may be any class that supports a fetch(chrom, start, end) method.
jpayne@69 40 #
jpayne@69 41 # NOTE: the position that is returned to Python is 0-based, NOT
jpayne@69 42 # 1-based as in the VCF file.
jpayne@69 43 # NOTE: There is also preliminary VCF functionality in the VariantFile class.
jpayne@69 44 #
jpayne@69 45 # TODO:
jpayne@69 46 # only v4.0 writing is complete; alleles are not converted to v3.3 format
jpayne@69 47 #
jpayne@69 48
jpayne@69 49 from collections import namedtuple, defaultdict
jpayne@69 50 from operator import itemgetter
jpayne@69 51 import sys, re, copy, bisect
jpayne@69 52
jpayne@69 53 from libc.stdlib cimport atoi
jpayne@69 54 from libc.stdint cimport int8_t, int16_t, int32_t, int64_t
jpayne@69 55 from libc.stdint cimport uint8_t, uint16_t, uint32_t, uint64_t
jpayne@69 56
jpayne@69 57 cimport pysam.libctabix as libctabix
jpayne@69 58 cimport pysam.libctabixproxies as libctabixproxies
jpayne@69 59
jpayne@69 60 from pysam.libcutils cimport force_str
jpayne@69 61
jpayne@69 62 import pysam
jpayne@69 63
jpayne@69 64 gtsRegEx = re.compile("[|/\\\\]")
jpayne@69 65 alleleRegEx = re.compile('^[ACGTN]+$')
jpayne@69 66
jpayne@69 67 # Utility function. Uses 0-based coordinates
jpayne@69 68 def get_sequence(chrom, start, end, fa):
jpayne@69 69 # obtain sequence from .fa file, without truncation
jpayne@69 70 if end<=start: return ""
jpayne@69 71 if not fa: return "N"*(end-start)
jpayne@69 72 if start<0: return "N"*(-start) + get_sequence(chrom, 0, end, fa).upper()
jpayne@69 73 sequence = fa.fetch(chrom, start, end).upper()
jpayne@69 74 if len(sequence) < end-start: sequence += "N"*(end-start-len(sequence))
jpayne@69 75 return sequence
jpayne@69 76
jpayne@69 77 # Utility function. Parses a region string
jpayne@69 78 def parse_regions( string ):
jpayne@69 79 result = []
jpayne@69 80 for r in string.split(','):
jpayne@69 81 elts = r.split(':')
jpayne@69 82 chrom, start, end = elts[0], 0, 3000000000
jpayne@69 83 if len(elts)==1: pass
jpayne@69 84 elif len(elts)==2:
jpayne@69 85 if len(elts[1])>0:
jpayne@69 86 ielts = elts[1].split('-')
jpayne@69 87 if len(ielts) != 2: ValueError("Don't understand region string '%s'" % r)
jpayne@69 88 try: start, end = int(ielts[0])-1, int(ielts[1])
jpayne@69 89 except: raise ValueError("Don't understand region string '%s'" % r)
jpayne@69 90 else:
jpayne@69 91 raise ValueError("Don't understand region string '%s'" % r)
jpayne@69 92 result.append( (chrom,start,end) )
jpayne@69 93 return result
jpayne@69 94
jpayne@69 95
jpayne@69 96 FORMAT = namedtuple('FORMAT','id numbertype number type description missingvalue')
jpayne@69 97
jpayne@69 98 ###########################################################################################################
jpayne@69 99 #
jpayne@69 100 # New class
jpayne@69 101 #
jpayne@69 102 ###########################################################################################################
jpayne@69 103
jpayne@69 104 cdef class VCFRecord(libctabixproxies.TupleProxy):
jpayne@69 105 '''vcf record.
jpayne@69 106
jpayne@69 107 initialized from data and vcf meta
jpayne@69 108 '''
jpayne@69 109
jpayne@69 110 cdef vcf
jpayne@69 111 cdef char * contig
jpayne@69 112 cdef uint32_t pos
jpayne@69 113
jpayne@69 114 def __init__(self, vcf):
jpayne@69 115 self.vcf = vcf
jpayne@69 116 self.encoding = vcf.encoding
jpayne@69 117
jpayne@69 118 # if len(data) != len(self.vcf._samples):
jpayne@69 119 # self.vcf.error(str(data),
jpayne@69 120 # self.BAD_NUMBER_OF_COLUMNS,
jpayne@69 121 # "expected %s for %s samples (%s), got %s" % \
jpayne@69 122 # (len(self.vcf._samples),
jpayne@69 123 # len(self.vcf._samples),
jpayne@69 124 # self.vcf._samples,
jpayne@69 125 # len(data)))
jpayne@69 126
jpayne@69 127 def __cinit__(self, vcf):
jpayne@69 128 # start indexed access at genotypes
jpayne@69 129 self.offset = 9
jpayne@69 130
jpayne@69 131 self.vcf = vcf
jpayne@69 132 self.encoding = vcf.encoding
jpayne@69 133
jpayne@69 134 def error(self, line, error, opt=None):
jpayne@69 135 '''raise error.'''
jpayne@69 136 # pass to vcf file for error handling
jpayne@69 137 return self.vcf.error(line, error, opt)
jpayne@69 138
jpayne@69 139 cdef update(self, char * buffer, size_t nbytes):
jpayne@69 140 '''update internal data.
jpayne@69 141
jpayne@69 142 nbytes does not include the terminal '\0'.
jpayne@69 143 '''
jpayne@69 144 libctabixproxies.TupleProxy.update(self, buffer, nbytes)
jpayne@69 145
jpayne@69 146 self.contig = self.fields[0]
jpayne@69 147 # vcf counts from 1 - correct here
jpayne@69 148 self.pos = atoi(self.fields[1]) - 1
jpayne@69 149
jpayne@69 150 def __len__(self):
jpayne@69 151 return max(0, self.nfields - 9)
jpayne@69 152
jpayne@69 153 property contig:
jpayne@69 154 def __get__(self): return self.contig
jpayne@69 155
jpayne@69 156 property pos:
jpayne@69 157 def __get__(self): return self.pos
jpayne@69 158
jpayne@69 159 property id:
jpayne@69 160 def __get__(self): return self.fields[2]
jpayne@69 161
jpayne@69 162 property ref:
jpayne@69 163 def __get__(self):
jpayne@69 164 return self.fields[3]
jpayne@69 165
jpayne@69 166 property alt:
jpayne@69 167 def __get__(self):
jpayne@69 168 # convert v3.3 to v4.0 alleles below
jpayne@69 169 alt = self.fields[4]
jpayne@69 170 if alt == ".": alt = []
jpayne@69 171 else: alt = alt.upper().split(',')
jpayne@69 172 return alt
jpayne@69 173
jpayne@69 174 property qual:
jpayne@69 175 def __get__(self):
jpayne@69 176 qual = self.fields[5]
jpayne@69 177 if qual == b".": qual = -1
jpayne@69 178 else:
jpayne@69 179 try: qual = float(qual)
jpayne@69 180 except: self.vcf.error(str(self),self.QUAL_NOT_NUMERICAL)
jpayne@69 181 return qual
jpayne@69 182
jpayne@69 183 property filter:
jpayne@69 184 def __get__(self):
jpayne@69 185 f = self.fields[6]
jpayne@69 186 # postpone checking that filters exist. Encode missing filter or no filtering as empty list
jpayne@69 187 if f == b"." or f == b"PASS" or f == b"0": return []
jpayne@69 188 else: return f.split(';')
jpayne@69 189
jpayne@69 190 property info:
jpayne@69 191 def __get__(self):
jpayne@69 192 col = self.fields[7]
jpayne@69 193 # dictionary of keys, and list of values
jpayne@69 194 info = {}
jpayne@69 195 if col != b".":
jpayne@69 196 for blurp in col.split(';'):
jpayne@69 197 elts = blurp.split('=')
jpayne@69 198 if len(elts) == 1: v = None
jpayne@69 199 elif len(elts) == 2: v = elts[1]
jpayne@69 200 else: self.vcf.error(str(self),self.ERROR_INFO_STRING)
jpayne@69 201 info[elts[0]] = self.vcf.parse_formatdata(elts[0], v, self.vcf._info, str(self.vcf))
jpayne@69 202 return info
jpayne@69 203
jpayne@69 204 property format:
jpayne@69 205 def __get__(self):
jpayne@69 206 return self.fields[8].split(':')
jpayne@69 207
jpayne@69 208 property samples:
jpayne@69 209 def __get__(self):
jpayne@69 210 return self.vcf._samples
jpayne@69 211
jpayne@69 212 def __getitem__(self, key):
jpayne@69 213
jpayne@69 214 # parse sample columns
jpayne@69 215 values = self.fields[self.vcf._sample2column[key]].split(':')
jpayne@69 216 alt = self.alt
jpayne@69 217 format = self.format
jpayne@69 218
jpayne@69 219 if len(values) > len(format):
jpayne@69 220 self.vcf.error(str(self.line),self.BAD_NUMBER_OF_VALUES,"(found %s values in element %s; expected %s)" %\
jpayne@69 221 (len(values),key,len(format)))
jpayne@69 222
jpayne@69 223 result = {}
jpayne@69 224 for idx in range(len(format)):
jpayne@69 225 expected = self.vcf.get_expected(format[idx], self.vcf._format, alt)
jpayne@69 226 if idx < len(values): value = values[idx]
jpayne@69 227 else:
jpayne@69 228 if expected == -1: value = "."
jpayne@69 229 else: value = ",".join(["."]*expected)
jpayne@69 230
jpayne@69 231 result[format[idx]] = self.vcf.parse_formatdata(format[idx], value, self.vcf._format, str(self.data))
jpayne@69 232 if expected != -1 and len(result[format[idx]]) != expected:
jpayne@69 233 self.vcf.error(str(self.data),self.BAD_NUMBER_OF_PARAMETERS,
jpayne@69 234 "id=%s, expected %s parameters, got %s" % (format[idx],expected,result[format[idx]]))
jpayne@69 235 if len(result[format[idx]] ) < expected: result[format[idx]] += [result[format[idx]][-1]]*(expected-len(result[format[idx]]))
jpayne@69 236 result[format[idx]] = result[format[idx]][:expected]
jpayne@69 237
jpayne@69 238 return result
jpayne@69 239
jpayne@69 240
jpayne@69 241 cdef class asVCFRecord(libctabix.Parser):
jpayne@69 242 '''converts a :term:`tabix row` into a VCF record.'''
jpayne@69 243 cdef vcffile
jpayne@69 244 def __init__(self, vcffile):
jpayne@69 245 self.vcffile = vcffile
jpayne@69 246
jpayne@69 247 cdef parse(self, char * buffer, int len):
jpayne@69 248 cdef VCFRecord r
jpayne@69 249 r = VCFRecord(self.vcffile)
jpayne@69 250 r.copy(buffer, len)
jpayne@69 251 return r
jpayne@69 252
jpayne@69 253 class VCF(object):
jpayne@69 254
jpayne@69 255 # types
jpayne@69 256 NT_UNKNOWN = 0
jpayne@69 257 NT_NUMBER = 1
jpayne@69 258 NT_ALLELES = 2
jpayne@69 259 NT_NR_ALLELES = 3
jpayne@69 260 NT_GENOTYPES = 4
jpayne@69 261 NT_PHASED_GENOTYPES = 5
jpayne@69 262
jpayne@69 263 _errors = { 0:"UNKNOWN_FORMAT_STRING:Unknown file format identifier",
jpayne@69 264 1:"BADLY_FORMATTED_FORMAT_STRING:Formatting error in the format string",
jpayne@69 265 2:"BADLY_FORMATTED_HEADING:Did not find 9 required headings (CHROM, POS, ..., FORMAT) %s",
jpayne@69 266 3:"BAD_NUMBER_OF_COLUMNS:Wrong number of columns found (%s)",
jpayne@69 267 4:"POS_NOT_NUMERICAL:Position column is not numerical",
jpayne@69 268 5:"UNKNOWN_CHAR_IN_REF:Unknown character in reference field",
jpayne@69 269 6:"V33_BAD_REF:Reference should be single-character in v3.3 VCF",
jpayne@69 270 7:"V33_BAD_ALLELE:Cannot interpret allele for v3.3 VCF",
jpayne@69 271 8:"POS_NOT_POSITIVE:Position field must be >0",
jpayne@69 272 9:"QUAL_NOT_NUMERICAL:Quality field must be numerical, or '.'",
jpayne@69 273 10:"ERROR_INFO_STRING:Error while parsing info field",
jpayne@69 274 11:"ERROR_UNKNOWN_KEY:Unknown key (%s) found in formatted field (info; format; or filter)",
jpayne@69 275 12:"ERROR_FORMAT_NOT_NUMERICAL:Expected integer or float in formatted field; got %s",
jpayne@69 276 13:"ERROR_FORMAT_NOT_CHAR:Eexpected character in formatted field; got string",
jpayne@69 277 14:"FILTER_NOT_DEFINED:Identifier (%s) in filter found which was not defined in header",
jpayne@69 278 15:"FORMAT_NOT_DEFINED:Identifier (%s) in format found which was not defined in header",
jpayne@69 279 16:"BAD_NUMBER_OF_VALUES:Found too many of values in sample column (%s)",
jpayne@69 280 17:"BAD_NUMBER_OF_PARAMETERS:Found unexpected number of parameters (%s)",
jpayne@69 281 18:"BAD_GENOTYPE:Cannot parse genotype (%s)",
jpayne@69 282 19:"V40_BAD_ALLELE:Bad allele found for v4.0 VCF (%s)",
jpayne@69 283 20:"MISSING_REF:Reference allele missing",
jpayne@69 284 21:"V33_UNMATCHED_DELETION:Deleted sequence does not match reference (%s)",
jpayne@69 285 22:"V40_MISSING_ANGLE_BRACKETS:Format definition is not deliminted by angular brackets",
jpayne@69 286 23:"FORMAT_MISSING_QUOTES:Description field in format definition is not surrounded by quotes",
jpayne@69 287 24:"V40_FORMAT_MUST_HAVE_NAMED_FIELDS:Fields in v4.0 VCF format definition must have named fields",
jpayne@69 288 25:"HEADING_NOT_SEPARATED_BY_TABS:Heading line appears separated by spaces, not tabs",
jpayne@69 289 26:"WRONG_REF:Wrong reference %s",
jpayne@69 290 27:"ERROR_TRAILING_DATA:Numerical field ('%s') has semicolon-separated trailing data",
jpayne@69 291 28:"BAD_CHR_TAG:Error calculating chr tag for %s",
jpayne@69 292 29:"ZERO_LENGTH_ALLELE:Found zero-length allele",
jpayne@69 293 30:"MISSING_INDEL_ALLELE_REF_BASE:Indel alleles must begin with single reference base",
jpayne@69 294 31:"ZERO_FOR_NON_FLAG_FIELD: number set to 0, but type is not 'FLAG'",
jpayne@69 295 32:"ERROR_FORMAT_NOT_INTEGER:Expected integer in formatted field; got %s",
jpayne@69 296 33:"ERROR_FLAG_HAS_VALUE:Flag fields should not have a value",
jpayne@69 297 }
jpayne@69 298
jpayne@69 299 # tag-value pairs; tags are not unique; does not include fileformat, INFO, FILTER or FORMAT fields
jpayne@69 300 _header = []
jpayne@69 301
jpayne@69 302 # version number; 33=v3.3; 40=v4.0
jpayne@69 303 _version = 40
jpayne@69 304
jpayne@69 305 # info, filter and format data
jpayne@69 306 _info = {}
jpayne@69 307 _filter = {}
jpayne@69 308 _format = {}
jpayne@69 309
jpayne@69 310 # header; and required columns
jpayne@69 311 _required = ["CHROM","POS","ID","REF","ALT","QUAL","FILTER","INFO","FORMAT"]
jpayne@69 312 _samples = []
jpayne@69 313
jpayne@69 314 # control behaviour
jpayne@69 315 _ignored_errors = set([11,31]) # ERROR_UNKNOWN_KEY, ERROR_ZERO_FOR_NON_FLAG_FIELD
jpayne@69 316 _warn_errors = set([])
jpayne@69 317 _leftalign = False
jpayne@69 318
jpayne@69 319 # reference sequence
jpayne@69 320 _reference = None
jpayne@69 321
jpayne@69 322 # regions to include; None includes everything
jpayne@69 323 _regions = None
jpayne@69 324
jpayne@69 325 # statefull stuff
jpayne@69 326 _lineno = -1
jpayne@69 327 _line = None
jpayne@69 328 _lines = None
jpayne@69 329
jpayne@69 330 def __init__(self, _copy=None, reference=None, regions=None,
jpayne@69 331 lines=None, leftalign=False):
jpayne@69 332 # make error identifiers accessible by name
jpayne@69 333 for id in self._errors.keys():
jpayne@69 334 self.__dict__[self._errors[id].split(':')[0]] = id
jpayne@69 335 if _copy != None:
jpayne@69 336 self._leftalign = _copy._leftalign
jpayne@69 337 self._header = _copy._header[:]
jpayne@69 338 self._version = _copy._version
jpayne@69 339 self._info = copy.deepcopy(_copy._info)
jpayne@69 340 self._filter = copy.deepcopy(_copy._filter)
jpayne@69 341 self._format = copy.deepcopy(_copy._format)
jpayne@69 342 self._samples = _copy._samples[:]
jpayne@69 343 self._sample2column = copy.deepcopy(_copy._sample2column)
jpayne@69 344 self._ignored_errors = copy.deepcopy(_copy._ignored_errors)
jpayne@69 345 self._warn_errors = copy.deepcopy(_copy._warn_errors)
jpayne@69 346 self._reference = _copy._reference
jpayne@69 347 self._regions = _copy._regions
jpayne@69 348 if reference: self._reference = reference
jpayne@69 349 if regions: self._regions = regions
jpayne@69 350 if leftalign: self._leftalign = leftalign
jpayne@69 351 self._lines = lines
jpayne@69 352 self.encoding = "ascii"
jpayne@69 353 self.tabixfile = None
jpayne@69 354
jpayne@69 355 def error(self,line,error,opt=None):
jpayne@69 356 if error in self._ignored_errors: return
jpayne@69 357 errorlabel, errorstring = self._errors[error].split(':')
jpayne@69 358 if opt: errorstring = errorstring % opt
jpayne@69 359 errwarn = ["Error","Warning"][error in self._warn_errors]
jpayne@69 360 errorstring += " in line %s: '%s'\n%s %s: %s\n" % (self._lineno,line,errwarn,errorlabel,errorstring)
jpayne@69 361 if error in self._warn_errors: return
jpayne@69 362 raise ValueError(errorstring)
jpayne@69 363
jpayne@69 364 def parse_format(self,line,format,filter=False):
jpayne@69 365 if self._version == 40:
jpayne@69 366 if not format.startswith('<'):
jpayne@69 367 self.error(line,self.V40_MISSING_ANGLE_BRACKETS)
jpayne@69 368 format = "<"+format
jpayne@69 369 if not format.endswith('>'):
jpayne@69 370 self.error(line,self.V40_MISSING_ANGLE_BRACKETS)
jpayne@69 371 format += ">"
jpayne@69 372 format = format[1:-1]
jpayne@69 373 data = {'id':None,'number':None,'type':None,'descr':None}
jpayne@69 374 idx = 0
jpayne@69 375 while len(format.strip())>0:
jpayne@69 376 elts = format.strip().split(',')
jpayne@69 377 first, rest = elts[0], ','.join(elts[1:])
jpayne@69 378 if first.find('=') == -1 or (first.find('"')>=0 and first.find('=') > first.find('"')):
jpayne@69 379 if self._version == 40: self.error(line,self.V40_FORMAT_MUST_HAVE_NAMED_FIELDS)
jpayne@69 380 if idx == 4: self.error(line,self.BADLY_FORMATTED_FORMAT_STRING)
jpayne@69 381 first = ["ID=","Number=","Type=","Description="][idx] + first
jpayne@69 382 if first.startswith('ID='): data['id'] = first.split('=')[1]
jpayne@69 383 elif first.startswith('Number='): data['number'] = first.split('=')[1]
jpayne@69 384 elif first.startswith('Type='): data['type'] = first.split('=')[1]
jpayne@69 385 elif first.startswith('Description='):
jpayne@69 386 elts = format.split('"')
jpayne@69 387 if len(elts)<3:
jpayne@69 388 self.error(line,self.FORMAT_MISSING_QUOTES)
jpayne@69 389 elts = first.split('=') + [rest]
jpayne@69 390 data['descr'] = elts[1]
jpayne@69 391 rest = '"'.join(elts[2:])
jpayne@69 392 if rest.startswith(','): rest = rest[1:]
jpayne@69 393 else:
jpayne@69 394 self.error(line,self.BADLY_FORMATTED_FORMAT_STRING)
jpayne@69 395 format = rest
jpayne@69 396 idx += 1
jpayne@69 397 if filter and idx==1: idx=3 # skip number and type fields for FILTER format strings
jpayne@69 398 if not data['id']: self.error(line,self.BADLY_FORMATTED_FORMAT_STRING)
jpayne@69 399 if 'descr' not in data:
jpayne@69 400 # missing description
jpayne@69 401 self.error(line,self.BADLY_FORMATTED_FORMAT_STRING)
jpayne@69 402 data['descr'] = ""
jpayne@69 403 if not data['type'] and not data['number']:
jpayne@69 404 # fine, ##filter format
jpayne@69 405 return FORMAT(data['id'],self.NT_NUMBER,0,"Flag",data['descr'],'.')
jpayne@69 406 if not data['type'] in ["Integer","Float","Character","String","Flag"]:
jpayne@69 407 self.error(line,self.BADLY_FORMATTED_FORMAT_STRING)
jpayne@69 408 # I would like a missing-value field, but it isn't there
jpayne@69 409 if data['type'] in ['Integer','Float']: data['missing'] = None # Do NOT use arbitrary int/float as missing value
jpayne@69 410 else: data['missing'] = '.'
jpayne@69 411 if not data['number']: self.error(line,self.BADLY_FORMATTED_FORMAT_STRING)
jpayne@69 412 try:
jpayne@69 413 n = int(data['number'])
jpayne@69 414 t = self.NT_NUMBER
jpayne@69 415 except ValueError:
jpayne@69 416 n = -1
jpayne@69 417 if data['number'] == '.': t = self.NT_UNKNOWN
jpayne@69 418 elif data['number'] == '#alleles': t = self.NT_ALLELES
jpayne@69 419 elif data['number'] == '#nonref_alleles': t = self.NT_NR_ALLELES
jpayne@69 420 elif data['number'] == '#genotypes': t = self.NT_GENOTYPES
jpayne@69 421 elif data['number'] == '#phased_genotypes': t = self.NT_PHASED_GENOTYPES
jpayne@69 422 elif data['number'] == '#phased_genotypes': t = self.NT_PHASED_GENOTYPES
jpayne@69 423 # abbreviations added in VCF version v4.1
jpayne@69 424 elif data['number'] == 'A': t = self.NT_ALLELES
jpayne@69 425 elif data['number'] == 'G': t = self.NT_GENOTYPES
jpayne@69 426 else:
jpayne@69 427 self.error(line,self.BADLY_FORMATTED_FORMAT_STRING)
jpayne@69 428 # if number is 0 - type must be Flag
jpayne@69 429 if n == 0 and data['type'] != 'Flag':
jpayne@69 430 self.error( line, self.ZERO_FOR_NON_FLAG_FIELD)
jpayne@69 431 # force type 'Flag' if no number
jpayne@69 432 data['type'] = 'Flag'
jpayne@69 433
jpayne@69 434 return FORMAT(data['id'],t,n,data['type'],data['descr'],data['missing'])
jpayne@69 435
jpayne@69 436 def format_format( self, fmt, filter=False ):
jpayne@69 437 values = [('ID',fmt.id)]
jpayne@69 438 if fmt.number != None and not filter:
jpayne@69 439 if fmt.numbertype == self.NT_UNKNOWN: nmb = "."
jpayne@69 440 elif fmt.numbertype == self.NT_NUMBER: nmb = str(fmt.number)
jpayne@69 441 elif fmt.numbertype == self.NT_ALLELES: nmb = "#alleles"
jpayne@69 442 elif fmt.numbertype == self.NT_NR_ALLELES: nmb = "#nonref_alleles"
jpayne@69 443 elif fmt.numbertype == self.NT_GENOTYPES: nmb = "#genotypes"
jpayne@69 444 elif fmt.numbertype == self.NT_PHASED_GENOTYPES: nmb = "#phased_genotypes"
jpayne@69 445 else:
jpayne@69 446 raise ValueError("Unknown number type encountered: %s" % fmt.numbertype)
jpayne@69 447 values.append( ('Number',nmb) )
jpayne@69 448 values.append( ('Type', fmt.type) )
jpayne@69 449 values.append( ('Description', '"' + fmt.description + '"') )
jpayne@69 450 if self._version == 33:
jpayne@69 451 format = ",".join([v for k,v in values])
jpayne@69 452 else:
jpayne@69 453 format = "<" + (",".join( ["%s=%s" % (k,v) for (k,v) in values] )) + ">"
jpayne@69 454 return format
jpayne@69 455
jpayne@69 456 def get_expected(self, format, formatdict, alt):
jpayne@69 457 fmt = formatdict[format]
jpayne@69 458 if fmt.numbertype == self.NT_UNKNOWN: return -1
jpayne@69 459 if fmt.numbertype == self.NT_NUMBER: return fmt.number
jpayne@69 460 if fmt.numbertype == self.NT_ALLELES: return len(alt)+1
jpayne@69 461 if fmt.numbertype == self.NT_NR_ALLELES: return len(alt)
jpayne@69 462 if fmt.numbertype == self.NT_GENOTYPES: return ((len(alt)+1)*(len(alt)+2)) // 2
jpayne@69 463 if fmt.numbertype == self.NT_PHASED_GENOTYPES: return (len(alt)+1)*(len(alt)+1)
jpayne@69 464 return 0
jpayne@69 465
jpayne@69 466
jpayne@69 467 def _add_definition(self, formatdict, key, data, line ):
jpayne@69 468 if key in formatdict: return
jpayne@69 469 self.error(line,self.ERROR_UNKNOWN_KEY,key)
jpayne@69 470 if data == None:
jpayne@69 471 formatdict[key] = FORMAT(key,self.NT_NUMBER,0,"Flag","(Undefined tag)",".")
jpayne@69 472 return
jpayne@69 473 if data == []: data = [""] # unsure what type -- say string
jpayne@69 474 if type(data[0]) == type(0.0):
jpayne@69 475 formatdict[key] = FORMAT(key,self.NT_UNKNOWN,-1,"Float","(Undefined tag)",None)
jpayne@69 476 return
jpayne@69 477 if type(data[0]) == type(0):
jpayne@69 478 formatdict[key] = FORMAT(key,self.NT_UNKNOWN,-1,"Integer","(Undefined tag)",None)
jpayne@69 479 return
jpayne@69 480 formatdict[key] = FORMAT(key,self.NT_UNKNOWN,-1,"String","(Undefined tag)",".")
jpayne@69 481
jpayne@69 482
jpayne@69 483 # todo: trim trailing missing values
jpayne@69 484 def format_formatdata( self, data, format, key=True, value=True, separator=":" ):
jpayne@69 485 output, sdata = [], []
jpayne@69 486 if type(data) == type([]): # for FORMAT field, make data with dummy values
jpayne@69 487 d = {}
jpayne@69 488 for k in data: d[k] = []
jpayne@69 489 data = d
jpayne@69 490 # convert missing values; and silently add definitions if required
jpayne@69 491 for k in data:
jpayne@69 492 self._add_definition( format, k, data[k], "(output)" )
jpayne@69 493 for idx,v in enumerate(data[k]):
jpayne@69 494 if v == format[k].missingvalue: data[k][idx] = "."
jpayne@69 495 # make sure GT comes first; and ensure fixed ordering; also convert GT data back to string
jpayne@69 496 for k in data:
jpayne@69 497 if k != 'GT': sdata.append( (k,data[k]) )
jpayne@69 498 sdata.sort()
jpayne@69 499 if 'GT' in data:
jpayne@69 500 sdata = [('GT',map(self.convertGTback,data['GT']))] + sdata
jpayne@69 501 for k,v in sdata:
jpayne@69 502 if v == []: v = None
jpayne@69 503 if key and value:
jpayne@69 504 if v != None: output.append( k+"="+','.join(map(str,v)) )
jpayne@69 505 else: output.append( k )
jpayne@69 506 elif key: output.append(k)
jpayne@69 507 elif value:
jpayne@69 508 if v != None: output.append( ','.join(map(str,v)) )
jpayne@69 509 else: output.append( "." ) # should not happen
jpayne@69 510 # snip off trailing missing data
jpayne@69 511 while len(output) > 1:
jpayne@69 512 last = output[-1].replace(',','').replace('.','')
jpayne@69 513 if len(last)>0: break
jpayne@69 514 output = output[:-1]
jpayne@69 515 return separator.join(output)
jpayne@69 516
jpayne@69 517
jpayne@69 518 def enter_default_format(self):
jpayne@69 519 for f in [FORMAT('GT',self.NT_NUMBER,1,'String','Genotype','.'),
jpayne@69 520 FORMAT('DP',self.NT_NUMBER,1,'Integer','Read depth at this position for this sample',-1),
jpayne@69 521 FORMAT('FT',self.NT_NUMBER,1,'String','Sample Genotype Filter','.'),
jpayne@69 522 FORMAT('GL',self.NT_UNKNOWN,-1,'Float','Genotype likelihoods','.'),
jpayne@69 523 FORMAT('GLE',self.NT_UNKNOWN,-1,'Float','Genotype likelihoods','.'),
jpayne@69 524 FORMAT('GQ',self.NT_NUMBER,1,'Integer','Genotype Quality',-1),
jpayne@69 525 FORMAT('PL',self.NT_GENOTYPES,-1,'Integer','Phred-scaled genotype likelihoods', '.'),
jpayne@69 526 FORMAT('GP',self.NT_GENOTYPES,-1,'Float','Genotype posterior probabilities','.'),
jpayne@69 527 FORMAT('GQ',self.NT_GENOTYPES,-1,'Integer','Conditional genotype quality','.'),
jpayne@69 528 FORMAT('HQ',self.NT_UNKNOWN,-1,'Integer','Haplotype Quality',-1), # unknown number, since may be haploid
jpayne@69 529 FORMAT('PS',self.NT_UNKNOWN,-1,'Integer','Phase set','.'),
jpayne@69 530 FORMAT('PQ',self.NT_NUMBER,1,'Integer','Phasing quality',-1),
jpayne@69 531 FORMAT('EC',self.NT_ALLELES,1,'Integer','Expected alternate allel counts',-1),
jpayne@69 532 FORMAT('MQ',self.NT_NUMBER,1,'Integer','RMS mapping quality',-1),
jpayne@69 533 ]:
jpayne@69 534 if f.id not in self._format:
jpayne@69 535 self._format[f.id] = f
jpayne@69 536
jpayne@69 537 def parse_header(self, line):
jpayne@69 538
jpayne@69 539 assert line.startswith('##')
jpayne@69 540 elts = line[2:].split('=')
jpayne@69 541 key = elts[0].strip()
jpayne@69 542 value = '='.join(elts[1:]).strip()
jpayne@69 543 if key == "fileformat":
jpayne@69 544 if value == "VCFv3.3":
jpayne@69 545 self._version = 33
jpayne@69 546 elif value == "VCFv4.0":
jpayne@69 547 self._version = 40
jpayne@69 548 elif value == "VCFv4.1":
jpayne@69 549 # AH - for testing
jpayne@69 550 self._version = 40
jpayne@69 551 elif value == "VCFv4.2":
jpayne@69 552 # AH - for testing
jpayne@69 553 self._version = 40
jpayne@69 554 else:
jpayne@69 555 self.error(line,self.UNKNOWN_FORMAT_STRING)
jpayne@69 556 elif key == "INFO":
jpayne@69 557 f = self.parse_format(line, value)
jpayne@69 558 self._info[ f.id ] = f
jpayne@69 559 elif key == "FILTER":
jpayne@69 560 f = self.parse_format(line, value, filter=True)
jpayne@69 561 self._filter[ f.id ] = f
jpayne@69 562 elif key == "FORMAT":
jpayne@69 563 f = self.parse_format(line, value)
jpayne@69 564 self._format[ f.id ] = f
jpayne@69 565 else:
jpayne@69 566 # keep other keys in the header field
jpayne@69 567 self._header.append( (key,value) )
jpayne@69 568
jpayne@69 569
jpayne@69 570 def write_header( self, stream ):
jpayne@69 571 stream.write("##fileformat=VCFv%s.%s\n" % (self._version // 10, self._version % 10))
jpayne@69 572 for key,value in self._header: stream.write("##%s=%s\n" % (key,value))
jpayne@69 573 for var,label in [(self._info,"INFO"),(self._filter,"FILTER"),(self._format,"FORMAT")]:
jpayne@69 574 for f in var.itervalues(): stream.write("##%s=%s\n" % (label,self.format_format(f,filter=(label=="FILTER"))))
jpayne@69 575
jpayne@69 576
jpayne@69 577 def parse_heading( self, line ):
jpayne@69 578 assert line.startswith('#')
jpayne@69 579 assert not line.startswith('##')
jpayne@69 580 headings = line[1:].split('\t')
jpayne@69 581 # test for 8, as FORMAT field might be missing
jpayne@69 582 if len(headings)==1 and len(line[1:].split()) >= 8:
jpayne@69 583 self.error(line,self.HEADING_NOT_SEPARATED_BY_TABS)
jpayne@69 584 headings = line[1:].split()
jpayne@69 585
jpayne@69 586 for i,s in enumerate(self._required):
jpayne@69 587
jpayne@69 588 if len(headings)<=i or headings[i] != s:
jpayne@69 589
jpayne@69 590 if len(headings) <= i:
jpayne@69 591 err = "(%sth entry not found)" % (i+1)
jpayne@69 592 else:
jpayne@69 593 err = "(found %s, expected %s)" % (headings[i],s)
jpayne@69 594
jpayne@69 595 #self.error(line,self.BADLY_FORMATTED_HEADING,err)
jpayne@69 596 # allow FORMAT column to be absent
jpayne@69 597 if len(headings) == 8:
jpayne@69 598 headings.append("FORMAT")
jpayne@69 599 else:
jpayne@69 600 self.error(line,self.BADLY_FORMATTED_HEADING,err)
jpayne@69 601
jpayne@69 602 self._samples = headings[9:]
jpayne@69 603 self._sample2column = dict( [(y,x+9) for x,y in enumerate( self._samples ) ] )
jpayne@69 604
jpayne@69 605 def write_heading( self, stream ):
jpayne@69 606 stream.write("#" + "\t".join(self._required + self._samples) + "\n")
jpayne@69 607
jpayne@69 608 def convertGT(self, GTstring):
jpayne@69 609 if GTstring == ".": return ["."]
jpayne@69 610 try:
jpayne@69 611 gts = gtsRegEx.split(GTstring)
jpayne@69 612 if len(gts) == 1: return [int(gts[0])]
jpayne@69 613 if len(gts) != 2: raise ValueError()
jpayne@69 614 if gts[0] == "." and gts[1] == ".": return [gts[0],GTstring[len(gts[0]):-len(gts[1])],gts[1]]
jpayne@69 615 return [int(gts[0]),GTstring[len(gts[0]):-len(gts[1])],int(gts[1])]
jpayne@69 616 except ValueError:
jpayne@69 617 self.error(self._line,self.BAD_GENOTYPE,GTstring)
jpayne@69 618 return [".","|","."]
jpayne@69 619
jpayne@69 620 def convertGTback(self, GTdata):
jpayne@69 621 return ''.join(map(str,GTdata))
jpayne@69 622
jpayne@69 623 def parse_formatdata( self, key, value, formatdict, line ):
jpayne@69 624 # To do: check that the right number of values is present
jpayne@69 625 f = formatdict.get(key,None)
jpayne@69 626 if f == None:
jpayne@69 627 self._add_definition(formatdict, key, value, line )
jpayne@69 628 f = formatdict[key]
jpayne@69 629 if f.type == "Flag":
jpayne@69 630 if value is not None: self.error(line,self.ERROR_FLAG_HAS_VALUE)
jpayne@69 631 return []
jpayne@69 632 values = value.split(',')
jpayne@69 633 # deal with trailing data in some early VCF files
jpayne@69 634 if f.type in ["Float","Integer"] and len(values)>0 and values[-1].find(';') > -1:
jpayne@69 635 self.error(line,self.ERROR_TRAILING_DATA,values[-1])
jpayne@69 636 values[-1] = values[-1].split(';')[0]
jpayne@69 637 if f.type == "Integer":
jpayne@69 638 for idx,v in enumerate(values):
jpayne@69 639 try:
jpayne@69 640 if v == ".": values[idx] = f.missingvalue
jpayne@69 641 else: values[idx] = int(v)
jpayne@69 642 except:
jpayne@69 643 self.error(line,self.ERROR_FORMAT_NOT_INTEGER,"%s=%s" % (key, str(values)))
jpayne@69 644 return [0] * len(values)
jpayne@69 645 return values
jpayne@69 646 elif f.type == "String":
jpayne@69 647 self._line = line
jpayne@69 648 if f.id == "GT": values = list(map( self.convertGT, values ))
jpayne@69 649 return values
jpayne@69 650 elif f.type == "Character":
jpayne@69 651 for v in values:
jpayne@69 652 if len(v) != 1: self.error(line,self.ERROR_FORMAT_NOT_CHAR)
jpayne@69 653 return values
jpayne@69 654 elif f.type == "Float":
jpayne@69 655 for idx,v in enumerate(values):
jpayne@69 656 if v == ".": values[idx] = f.missingvalue
jpayne@69 657 try: return list(map(float,values))
jpayne@69 658 except:
jpayne@69 659 self.error(line,self.ERROR_FORMAT_NOT_NUMERICAL,"%s=%s" % (key, str(values)))
jpayne@69 660 return [0.0] * len(values)
jpayne@69 661 else:
jpayne@69 662 # can't happen
jpayne@69 663 self.error(line,self.ERROR_INFO_STRING)
jpayne@69 664
jpayne@69 665 def inregion(self, chrom, pos):
jpayne@69 666 if not self._regions: return True
jpayne@69 667 for r in self._regions:
jpayne@69 668 if r[0] == chrom and r[1] <= pos < r[2]: return True
jpayne@69 669 return False
jpayne@69 670
jpayne@69 671 def parse_data( self, line, lineparse=False ):
jpayne@69 672 cols = line.split('\t')
jpayne@69 673 if len(cols) != len(self._samples)+9:
jpayne@69 674 # gracefully deal with absent FORMAT column
jpayne@69 675 # and those missing samples
jpayne@69 676 if len(cols) == 8:
jpayne@69 677 cols.append("")
jpayne@69 678 else:
jpayne@69 679 self.error(line,
jpayne@69 680 self.BAD_NUMBER_OF_COLUMNS,
jpayne@69 681 "expected %s for %s samples (%s), got %s" % (len(self._samples)+9, len(self._samples), self._samples, len(cols)))
jpayne@69 682
jpayne@69 683 chrom = cols[0]
jpayne@69 684
jpayne@69 685 # get 0-based position
jpayne@69 686 try: pos = int(cols[1])-1
jpayne@69 687 except: self.error(line,self.POS_NOT_NUMERICAL)
jpayne@69 688 if pos < 0: self.error(line,self.POS_NOT_POSITIVE)
jpayne@69 689
jpayne@69 690 # implement filtering
jpayne@69 691 if not self.inregion(chrom,pos): return None
jpayne@69 692
jpayne@69 693 # end of first-pass parse for sortedVCF
jpayne@69 694 if lineparse: return chrom, pos, line
jpayne@69 695
jpayne@69 696 id = cols[2]
jpayne@69 697
jpayne@69 698 ref = cols[3].upper()
jpayne@69 699 if ref == ".":
jpayne@69 700 self.error(line,self.MISSING_REF)
jpayne@69 701 if self._version == 33: ref = get_sequence(chrom,pos,pos+1,self._reference)
jpayne@69 702 else: ref = ""
jpayne@69 703 else:
jpayne@69 704 for c in ref:
jpayne@69 705 if c not in "ACGTN": self.error(line,self.UNKNOWN_CHAR_IN_REF)
jpayne@69 706 if "N" in ref: ref = get_sequence(chrom,pos,pos+len(ref),self._reference)
jpayne@69 707
jpayne@69 708 # make sure reference is sane
jpayne@69 709 if self._reference:
jpayne@69 710 left = max(0,pos-100)
jpayne@69 711 faref_leftflank = get_sequence(chrom,left,pos+len(ref),self._reference)
jpayne@69 712 faref = faref_leftflank[pos-left:]
jpayne@69 713 if faref != ref: self.error(line,self.WRONG_REF,"(reference is %s, VCF says %s)" % (faref,ref))
jpayne@69 714 ref = faref
jpayne@69 715
jpayne@69 716 # convert v3.3 to v4.0 alleles below
jpayne@69 717 if cols[4] == ".": alt = []
jpayne@69 718 else: alt = cols[4].upper().split(',')
jpayne@69 719
jpayne@69 720 if cols[5] == ".": qual = -1
jpayne@69 721 else:
jpayne@69 722 try: qual = float(cols[5])
jpayne@69 723 except: self.error(line,self.QUAL_NOT_NUMERICAL)
jpayne@69 724
jpayne@69 725 # postpone checking that filters exist. Encode missing filter or no filtering as empty list
jpayne@69 726 if cols[6] == "." or cols[6] == "PASS" or cols[6] == "0": filter = []
jpayne@69 727 else: filter = cols[6].split(';')
jpayne@69 728
jpayne@69 729 # dictionary of keys, and list of values
jpayne@69 730 info = {}
jpayne@69 731 if cols[7] != ".":
jpayne@69 732 for blurp in cols[7].split(';'):
jpayne@69 733 elts = blurp.split('=')
jpayne@69 734 if len(elts) == 1: v = None
jpayne@69 735 elif len(elts) == 2: v = elts[1]
jpayne@69 736 else: self.error(line,self.ERROR_INFO_STRING)
jpayne@69 737 info[elts[0]] = self.parse_formatdata(elts[0],
jpayne@69 738 v,
jpayne@69 739 self._info,
jpayne@69 740 line)
jpayne@69 741
jpayne@69 742 # Gracefully deal with absent FORMAT column
jpayne@69 743 if cols[8] == "": format = []
jpayne@69 744 else: format = cols[8].split(':')
jpayne@69 745
jpayne@69 746 # check: all filters are defined
jpayne@69 747 for f in filter:
jpayne@69 748 if f not in self._filter: self.error(line,self.FILTER_NOT_DEFINED, f)
jpayne@69 749
jpayne@69 750 # check: format fields are defined
jpayne@69 751 if self._format:
jpayne@69 752 for f in format:
jpayne@69 753 if f not in self._format: self.error(line,self.FORMAT_NOT_DEFINED, f)
jpayne@69 754
jpayne@69 755 # convert v3.3 alleles
jpayne@69 756 if self._version == 33:
jpayne@69 757 if len(ref) != 1: self.error(line,self.V33_BAD_REF)
jpayne@69 758 newalts = []
jpayne@69 759 have_deletions = False
jpayne@69 760 for a in alt:
jpayne@69 761 if len(a) == 1: a = a + ref[1:] # SNP; add trailing reference
jpayne@69 762 elif a.startswith('I'): a = ref[0] + a[1:] + ref[1:] # insertion just beyond pos; add first and trailing reference
jpayne@69 763 elif a.startswith('D'): # allow D<seq> and D<num>
jpayne@69 764 have_deletions = True
jpayne@69 765 try:
jpayne@69 766 l = int(a[1:]) # throws ValueError if sequence
jpayne@69 767 if len(ref) < l: # add to reference if necessary
jpayne@69 768 addns = get_sequence(chrom,pos+len(ref),pos+l,self._reference)
jpayne@69 769 ref += addns
jpayne@69 770 for i,na in enumerate(newalts): newalts[i] = na+addns
jpayne@69 771 a = ref[l:] # new deletion, deleting pos...pos+l
jpayne@69 772 except ValueError:
jpayne@69 773 s = a[1:]
jpayne@69 774 if len(ref) < len(s): # add Ns to reference if necessary
jpayne@69 775 addns = get_sequence(chrom,pos+len(ref),pos+len(s),self._reference)
jpayne@69 776 if not s.endswith(addns) and addns != 'N'*len(addns):
jpayne@69 777 self.error(line,self.V33_UNMATCHED_DELETION,
jpayne@69 778 "(deletion is %s, reference is %s)" % (a,get_sequence(chrom,pos,pos+len(s),self._reference)))
jpayne@69 779 ref += addns
jpayne@69 780 for i,na in enumerate(newalts): newalts[i] = na+addns
jpayne@69 781 a = ref[len(s):] # new deletion, deleting from pos
jpayne@69 782 else:
jpayne@69 783 self.error(line,self.V33_BAD_ALLELE)
jpayne@69 784 newalts.append(a)
jpayne@69 785 alt = newalts
jpayne@69 786 # deletion alleles exist, add dummy 1st reference allele, and account for leading base
jpayne@69 787 if have_deletions:
jpayne@69 788 if pos == 0:
jpayne@69 789 # Petr Danacek's: we can't have a leading nucleotide at (1-based) position 1
jpayne@69 790 addn = get_sequence(chrom,pos+len(ref),pos+len(ref)+1,self._reference)
jpayne@69 791 ref += addn
jpayne@69 792 alt = [allele+addn for allele in alt]
jpayne@69 793 else:
jpayne@69 794 addn = get_sequence(chrom,pos-1,pos,self._reference)
jpayne@69 795 ref = addn + ref
jpayne@69 796 alt = [addn + allele for allele in alt]
jpayne@69 797 pos -= 1
jpayne@69 798 else:
jpayne@69 799 # format v4.0 -- just check for nucleotides
jpayne@69 800 for allele in alt:
jpayne@69 801 if not alleleRegEx.match(allele):
jpayne@69 802 self.error(line,self.V40_BAD_ALLELE,allele)
jpayne@69 803
jpayne@69 804 # check for leading nucleotide in indel calls
jpayne@69 805 for allele in alt:
jpayne@69 806 if len(allele) != len(ref):
jpayne@69 807 if len(allele) == 0: self.error(line,self.ZERO_LENGTH_ALLELE)
jpayne@69 808 if ref[0].upper() != allele[0].upper() and "N" not in (ref[0]+allele[0]).upper():
jpayne@69 809 self.error(line,self.MISSING_INDEL_ALLELE_REF_BASE)
jpayne@69 810
jpayne@69 811 # trim trailing bases in alleles
jpayne@69 812 # AH: not certain why trimming this needs to be added
jpayne@69 813 # disabled now for unit testing
jpayne@69 814 # if alt:
jpayne@69 815 # for i in range(1,min(len(ref),min(map(len,alt)))):
jpayne@69 816 # if len(set(allele[-1].upper() for allele in alt)) > 1 or ref[-1].upper() != alt[0][-1].upper():
jpayne@69 817 # break
jpayne@69 818 # ref, alt = ref[:-1], [allele[:-1] for allele in alt]
jpayne@69 819
jpayne@69 820 # left-align alleles, if a reference is available
jpayne@69 821 if self._leftalign and self._reference:
jpayne@69 822 while left < pos:
jpayne@69 823 movable = True
jpayne@69 824 for allele in alt:
jpayne@69 825 if len(allele) > len(ref):
jpayne@69 826 longest, shortest = allele, ref
jpayne@69 827 else:
jpayne@69 828 longest, shortest = ref, allele
jpayne@69 829 if len(longest) == len(shortest) or longest[:len(shortest)].upper() != shortest.upper():
jpayne@69 830 movable = False
jpayne@69 831 if longest[-1].upper() != longest[len(shortest)-1].upper():
jpayne@69 832 movable = False
jpayne@69 833 if not movable:
jpayne@69 834 break
jpayne@69 835 ref = ref[:-1]
jpayne@69 836 alt = [allele[:-1] for allele in alt]
jpayne@69 837 if min([len(allele) for allele in alt]) == 0 or len(ref) == 0:
jpayne@69 838 ref = faref_leftflank[pos-left-1] + ref
jpayne@69 839 alt = [faref_leftflank[pos-left-1] + allele for allele in alt]
jpayne@69 840 pos -= 1
jpayne@69 841
jpayne@69 842 # parse sample columns
jpayne@69 843 samples = []
jpayne@69 844 for sample in cols[9:]:
jpayne@69 845 dict = {}
jpayne@69 846 values = sample.split(':')
jpayne@69 847 if len(values) > len(format):
jpayne@69 848 self.error(line,self.BAD_NUMBER_OF_VALUES,"(found %s values in element %s; expected %s)" % (len(values),sample,len(format)))
jpayne@69 849 for idx in range(len(format)):
jpayne@69 850 expected = self.get_expected(format[idx], self._format, alt)
jpayne@69 851 if idx < len(values): value = values[idx]
jpayne@69 852 else:
jpayne@69 853 if expected == -1: value = "."
jpayne@69 854 else: value = ",".join(["."]*expected)
jpayne@69 855
jpayne@69 856 dict[format[idx]] = self.parse_formatdata(format[idx],
jpayne@69 857 value,
jpayne@69 858 self._format,
jpayne@69 859 line)
jpayne@69 860 if expected != -1 and len(dict[format[idx]]) != expected:
jpayne@69 861 self.error(line,self.BAD_NUMBER_OF_PARAMETERS,
jpayne@69 862 "id=%s, expected %s parameters, got %s" % (format[idx],expected,dict[format[idx]]))
jpayne@69 863 if len(dict[format[idx]] ) < expected: dict[format[idx]] += [dict[format[idx]][-1]]*(expected-len(dict[format[idx]]))
jpayne@69 864 dict[format[idx]] = dict[format[idx]][:expected]
jpayne@69 865 samples.append( dict )
jpayne@69 866
jpayne@69 867 # done
jpayne@69 868 d = {'chrom':chrom,
jpayne@69 869 'pos':pos, # return 0-based position
jpayne@69 870 'id':id,
jpayne@69 871 'ref':ref,
jpayne@69 872 'alt':alt,
jpayne@69 873 'qual':qual,
jpayne@69 874 'filter':filter,
jpayne@69 875 'info':info,
jpayne@69 876 'format':format}
jpayne@69 877 for key,value in zip(self._samples,samples):
jpayne@69 878 d[key] = value
jpayne@69 879
jpayne@69 880 return d
jpayne@69 881
jpayne@69 882
jpayne@69 883 def write_data(self, stream, data):
jpayne@69 884 required = ['chrom','pos','id','ref','alt','qual','filter','info','format'] + self._samples
jpayne@69 885 for k in required:
jpayne@69 886 if k not in data: raise ValueError("Required key %s not found in data" % str(k))
jpayne@69 887 if data['alt'] == []: alt = "."
jpayne@69 888 else: alt = ",".join(data['alt'])
jpayne@69 889 if data['filter'] == None: filter = "."
jpayne@69 890 elif data['filter'] == []:
jpayne@69 891 if self._version == 33: filter = "0"
jpayne@69 892 else: filter = "PASS"
jpayne@69 893 else: filter = ';'.join(data['filter'])
jpayne@69 894 if data['qual'] == -1: qual = "."
jpayne@69 895 else: qual = str(data['qual'])
jpayne@69 896
jpayne@69 897 output = [data['chrom'],
jpayne@69 898 str(data['pos']+1), # change to 1-based position
jpayne@69 899 data['id'],
jpayne@69 900 data['ref'],
jpayne@69 901 alt,
jpayne@69 902 qual,
jpayne@69 903 filter,
jpayne@69 904 self.format_formatdata(
jpayne@69 905 data['info'], self._info, separator=";"),
jpayne@69 906 self.format_formatdata(
jpayne@69 907 data['format'], self._format, value=False)]
jpayne@69 908
jpayne@69 909 for s in self._samples:
jpayne@69 910 output.append(self.format_formatdata(
jpayne@69 911 data[s], self._format, key=False))
jpayne@69 912
jpayne@69 913 stream.write( "\t".join(output) + "\n" )
jpayne@69 914
jpayne@69 915 def _parse_header(self, stream):
jpayne@69 916 self._lineno = 0
jpayne@69 917 for line in stream:
jpayne@69 918 line = force_str(line, self.encoding)
jpayne@69 919 self._lineno += 1
jpayne@69 920 if line.startswith('##'):
jpayne@69 921 self.parse_header(line.strip())
jpayne@69 922 elif line.startswith('#'):
jpayne@69 923 self.parse_heading(line.strip())
jpayne@69 924 self.enter_default_format()
jpayne@69 925 else:
jpayne@69 926 break
jpayne@69 927 return line
jpayne@69 928
jpayne@69 929 def _parse(self, line, stream):
jpayne@69 930 # deal with files with header only
jpayne@69 931 if line.startswith("##"): return
jpayne@69 932 if len(line.strip()) > 0:
jpayne@69 933 d = self.parse_data( line.strip() )
jpayne@69 934 if d: yield d
jpayne@69 935 for line in stream:
jpayne@69 936 self._lineno += 1
jpayne@69 937 if self._lines and self._lineno > self._lines: raise StopIteration
jpayne@69 938 d = self.parse_data( line.strip() )
jpayne@69 939 if d: yield d
jpayne@69 940
jpayne@69 941 ######################################################################################################
jpayne@69 942 #
jpayne@69 943 # API follows
jpayne@69 944 #
jpayne@69 945 ######################################################################################################
jpayne@69 946
jpayne@69 947 def getsamples(self):
jpayne@69 948 """ List of samples in VCF file """
jpayne@69 949 return self._samples
jpayne@69 950
jpayne@69 951 def setsamples(self,samples):
jpayne@69 952 """ List of samples in VCF file """
jpayne@69 953 self._samples = samples
jpayne@69 954
jpayne@69 955 def getheader(self):
jpayne@69 956 """ List of header key-value pairs (strings) """
jpayne@69 957 return self._header
jpayne@69 958
jpayne@69 959 def setheader(self,header):
jpayne@69 960 """ List of header key-value pairs (strings) """
jpayne@69 961 self._header = header
jpayne@69 962
jpayne@69 963 def getinfo(self):
jpayne@69 964 """ Dictionary of ##INFO tags, as VCF.FORMAT values """
jpayne@69 965 return self._info
jpayne@69 966
jpayne@69 967 def setinfo(self,info):
jpayne@69 968 """ Dictionary of ##INFO tags, as VCF.FORMAT values """
jpayne@69 969 self._info = info
jpayne@69 970
jpayne@69 971 def getformat(self):
jpayne@69 972 """ Dictionary of ##FORMAT tags, as VCF.FORMAT values """
jpayne@69 973 return self._format
jpayne@69 974
jpayne@69 975 def setformat(self,format):
jpayne@69 976 """ Dictionary of ##FORMAT tags, as VCF.FORMAT values """
jpayne@69 977 self._format = format
jpayne@69 978
jpayne@69 979 def getfilter(self):
jpayne@69 980 """ Dictionary of ##FILTER tags, as VCF.FORMAT values """
jpayne@69 981 return self._filter
jpayne@69 982
jpayne@69 983 def setfilter(self,filter):
jpayne@69 984 """ Dictionary of ##FILTER tags, as VCF.FORMAT values """
jpayne@69 985 self._filter = filter
jpayne@69 986
jpayne@69 987 def setversion(self, version):
jpayne@69 988 if version != 33 and version != 40: raise ValueError("Can only handle v3.3 and v4.0 VCF files")
jpayne@69 989 self._version = version
jpayne@69 990
jpayne@69 991 def setregions(self, regions):
jpayne@69 992 self._regions = regions
jpayne@69 993
jpayne@69 994 def setreference(self, ref):
jpayne@69 995 """ Provide a reference sequence; a Python class supporting a fetch(chromosome, start, end) method, e.g. PySam.FastaFile """
jpayne@69 996 self._reference = ref
jpayne@69 997
jpayne@69 998 def ignoreerror(self, errorstring):
jpayne@69 999 try: self._ignored_errors.add(self.__dict__[errorstring])
jpayne@69 1000 except KeyError: raise ValueError("Invalid error string: %s" % errorstring)
jpayne@69 1001
jpayne@69 1002 def warnerror(self, errorstring):
jpayne@69 1003 try: self._warn_errors.add(self.__dict__[errorstring])
jpayne@69 1004 except KeyError: raise ValueError("Invalid error string: %s" % errorstring)
jpayne@69 1005
jpayne@69 1006 def parse(self, stream):
jpayne@69 1007 """ Parse a stream of VCF-formatted lines. Initializes class instance and return generator """
jpayne@69 1008 last_line = self._parse_header(stream)
jpayne@69 1009 # now return a generator that does the actual work. In this way the pre-processing is done
jpayne@69 1010 # before the first piece of data is yielded
jpayne@69 1011 return self._parse(last_line, stream)
jpayne@69 1012
jpayne@69 1013 def write(self, stream, datagenerator):
jpayne@69 1014 """ Writes a VCF file to a stream, using a data generator (or list) """
jpayne@69 1015 self.write_header(stream)
jpayne@69 1016 self.write_heading(stream)
jpayne@69 1017 for data in datagenerator: self.write_data(stream,data)
jpayne@69 1018
jpayne@69 1019 def writeheader(self, stream):
jpayne@69 1020 """ Writes a VCF header """
jpayne@69 1021 self.write_header(stream)
jpayne@69 1022 self.write_heading(stream)
jpayne@69 1023
jpayne@69 1024 def compare_calls(self, pos1, ref1, alt1, pos2, ref2, alt2):
jpayne@69 1025 """ Utility function: compares two calls for equality """
jpayne@69 1026 # a variant should always be assigned to a unique position, one base before
jpayne@69 1027 # the leftmost position of the alignment gap. If this rule is implemented
jpayne@69 1028 # correctly, the two positions must be equal for the calls to be identical.
jpayne@69 1029 if pos1 != pos2: return False
jpayne@69 1030 # from both calls, trim rightmost bases when identical. Do this safely, i.e.
jpayne@69 1031 # only when the reference bases are not Ns
jpayne@69 1032 while len(ref1)>0 and len(alt1)>0 and ref1[-1] == alt1[-1]:
jpayne@69 1033 ref1 = ref1[:-1]
jpayne@69 1034 alt1 = alt1[:-1]
jpayne@69 1035 while len(ref2)>0 and len(alt2)>0 and ref2[-1] == alt2[-1]:
jpayne@69 1036 ref2 = ref2[:-1]
jpayne@69 1037 alt2 = alt2[:-1]
jpayne@69 1038 # now, the alternative alleles must be identical
jpayne@69 1039 return alt1 == alt2
jpayne@69 1040
jpayne@69 1041 ###########################################################################################################
jpayne@69 1042 ###########################################################################################################
jpayne@69 1043 ## API functions added by Andreas
jpayne@69 1044 ###########################################################################################################
jpayne@69 1045
jpayne@69 1046 def connect(self, filename, encoding="ascii"):
jpayne@69 1047 '''connect to tabix file.'''
jpayne@69 1048 self.encoding=encoding
jpayne@69 1049 self.tabixfile = pysam.Tabixfile(filename, encoding=encoding)
jpayne@69 1050 self._parse_header(self.tabixfile.header)
jpayne@69 1051
jpayne@69 1052 def __del__(self):
jpayne@69 1053 self.close()
jpayne@69 1054 self.tabixfile = None
jpayne@69 1055
jpayne@69 1056 def close(self):
jpayne@69 1057 if self.tabixfile:
jpayne@69 1058 self.tabixfile.close()
jpayne@69 1059 self.tabixfile = None
jpayne@69 1060
jpayne@69 1061 def fetch(self,
jpayne@69 1062 reference=None,
jpayne@69 1063 start=None,
jpayne@69 1064 end=None,
jpayne@69 1065 region=None ):
jpayne@69 1066 """ Parse a stream of VCF-formatted lines.
jpayne@69 1067 Initializes class instance and return generator """
jpayne@69 1068 return self.tabixfile.fetch(
jpayne@69 1069 reference,
jpayne@69 1070 start,
jpayne@69 1071 end,
jpayne@69 1072 region,
jpayne@69 1073 parser = asVCFRecord(self))
jpayne@69 1074
jpayne@69 1075 def validate(self, record):
jpayne@69 1076 '''validate vcf record.
jpayne@69 1077
jpayne@69 1078 returns a validated record.
jpayne@69 1079 '''
jpayne@69 1080
jpayne@69 1081 raise NotImplementedError("needs to be checked")
jpayne@69 1082
jpayne@69 1083 chrom, pos = record.chrom, record.pos
jpayne@69 1084
jpayne@69 1085 # check reference
jpayne@69 1086 ref = record.ref
jpayne@69 1087 if ref == ".":
jpayne@69 1088 self.error(str(record),self.MISSING_REF)
jpayne@69 1089 if self._version == 33: ref = get_sequence(chrom,pos,pos+1,self._reference)
jpayne@69 1090 else: ref = ""
jpayne@69 1091 else:
jpayne@69 1092 for c in ref:
jpayne@69 1093 if c not in "ACGTN": self.error(str(record),self.UNKNOWN_CHAR_IN_REF)
jpayne@69 1094 if "N" in ref: ref = get_sequence(chrom,
jpayne@69 1095 pos,
jpayne@69 1096 pos+len(ref),
jpayne@69 1097 self._reference)
jpayne@69 1098
jpayne@69 1099 # make sure reference is sane
jpayne@69 1100 if self._reference:
jpayne@69 1101 left = max(0,self.pos-100)
jpayne@69 1102 faref_leftflank = get_sequence(chrom,left,self.pos+len(ref),self._reference)
jpayne@69 1103 faref = faref_leftflank[pos-left:]
jpayne@69 1104 if faref != ref: self.error(str(record),self.WRONG_REF,"(reference is %s, VCF says %s)" % (faref,ref))
jpayne@69 1105 ref = faref
jpayne@69 1106
jpayne@69 1107 # check: format fields are defined
jpayne@69 1108 for f in record.format:
jpayne@69 1109 if f not in self._format: self.error(str(record),self.FORMAT_NOT_DEFINED, f)
jpayne@69 1110
jpayne@69 1111 # check: all filters are defined
jpayne@69 1112 for f in record.filter:
jpayne@69 1113 if f not in self._filter: self.error(str(record),self.FILTER_NOT_DEFINED, f)
jpayne@69 1114
jpayne@69 1115 # convert v3.3 alleles
jpayne@69 1116 if self._version == 33:
jpayne@69 1117 if len(ref) != 1: self.error(str(record),self.V33_BAD_REF)
jpayne@69 1118 newalts = []
jpayne@69 1119 have_deletions = False
jpayne@69 1120 for a in alt:
jpayne@69 1121 if len(a) == 1: a = a + ref[1:] # SNP; add trailing reference
jpayne@69 1122 elif a.startswith('I'): a = ref[0] + a[1:] + ref[1:] # insertion just beyond pos; add first and trailing reference
jpayne@69 1123 elif a.startswith('D'): # allow D<seq> and D<num>
jpayne@69 1124 have_deletions = True
jpayne@69 1125 try:
jpayne@69 1126 l = int(a[1:]) # throws ValueError if sequence
jpayne@69 1127 if len(ref) < l: # add to reference if necessary
jpayne@69 1128 addns = get_sequence(chrom,pos+len(ref),pos+l,self._reference)
jpayne@69 1129 ref += addns
jpayne@69 1130 for i,na in enumerate(newalts): newalts[i] = na+addns
jpayne@69 1131 a = ref[l:] # new deletion, deleting pos...pos+l
jpayne@69 1132 except ValueError:
jpayne@69 1133 s = a[1:]
jpayne@69 1134 if len(ref) < len(s): # add Ns to reference if necessary
jpayne@69 1135 addns = get_sequence(chrom,pos+len(ref),pos+len(s),self._reference)
jpayne@69 1136 if not s.endswith(addns) and addns != 'N'*len(addns):
jpayne@69 1137 self.error(str(record),self.V33_UNMATCHED_DELETION,
jpayne@69 1138 "(deletion is %s, reference is %s)" % (a,get_sequence(chrom,pos,pos+len(s),self._reference)))
jpayne@69 1139 ref += addns
jpayne@69 1140 for i,na in enumerate(newalts): newalts[i] = na+addns
jpayne@69 1141 a = ref[len(s):] # new deletion, deleting from pos
jpayne@69 1142 else:
jpayne@69 1143 self.error(str(record),self.V33_BAD_ALLELE)
jpayne@69 1144 newalts.append(a)
jpayne@69 1145 alt = newalts
jpayne@69 1146 # deletion alleles exist, add dummy 1st reference allele, and account for leading base
jpayne@69 1147 if have_deletions:
jpayne@69 1148 if pos == 0:
jpayne@69 1149 # Petr Danacek's: we can't have a leading nucleotide at (1-based) position 1
jpayne@69 1150 addn = get_sequence(chrom,pos+len(ref),pos+len(ref)+1,self._reference)
jpayne@69 1151 ref += addn
jpayne@69 1152 alt = [allele+addn for allele in alt]
jpayne@69 1153 else:
jpayne@69 1154 addn = get_sequence(chrom,pos-1,pos,self._reference)
jpayne@69 1155 ref = addn + ref
jpayne@69 1156 alt = [addn + allele for allele in alt]
jpayne@69 1157 pos -= 1
jpayne@69 1158 else:
jpayne@69 1159 # format v4.0 -- just check for nucleotides
jpayne@69 1160 for allele in alt:
jpayne@69 1161 if not alleleRegEx.match(allele):
jpayne@69 1162 self.error(str(record),self.V40_BAD_ALLELE,allele)
jpayne@69 1163
jpayne@69 1164
jpayne@69 1165 # check for leading nucleotide in indel calls
jpayne@69 1166 for allele in alt:
jpayne@69 1167 if len(allele) != len(ref):
jpayne@69 1168 if len(allele) == 0: self.error(str(record),self.ZERO_LENGTH_ALLELE)
jpayne@69 1169 if ref[0].upper() != allele[0].upper() and "N" not in (ref[0]+allele[0]).upper():
jpayne@69 1170 self.error(str(record),self.MISSING_INDEL_ALLELE_REF_BASE)
jpayne@69 1171
jpayne@69 1172 # trim trailing bases in alleles
jpayne@69 1173 # AH: not certain why trimming this needs to be added
jpayne@69 1174 # disabled now for unit testing
jpayne@69 1175 # for i in range(1,min(len(ref),min(map(len,alt)))):
jpayne@69 1176 # if len(set(allele[-1].upper() for allele in alt)) > 1 or ref[-1].upper() != alt[0][-1].upper():
jpayne@69 1177 # break
jpayne@69 1178 # ref, alt = ref[:-1], [allele[:-1] for allele in alt]
jpayne@69 1179
jpayne@69 1180 # left-align alleles, if a reference is available
jpayne@69 1181 if self._leftalign and self._reference:
jpayne@69 1182 while left < pos:
jpayne@69 1183 movable = True
jpayne@69 1184 for allele in alt:
jpayne@69 1185 if len(allele) > len(ref):
jpayne@69 1186 longest, shortest = allele, ref
jpayne@69 1187 else:
jpayne@69 1188 longest, shortest = ref, allele
jpayne@69 1189 if len(longest) == len(shortest) or longest[:len(shortest)].upper() != shortest.upper():
jpayne@69 1190 movable = False
jpayne@69 1191 if longest[-1].upper() != longest[len(shortest)-1].upper():
jpayne@69 1192 movable = False
jpayne@69 1193 if not movable:
jpayne@69 1194 break
jpayne@69 1195 ref = ref[:-1]
jpayne@69 1196 alt = [allele[:-1] for allele in alt]
jpayne@69 1197 if min([len(allele) for allele in alt]) == 0 or len(ref) == 0:
jpayne@69 1198 ref = faref_leftflank[pos-left-1] + ref
jpayne@69 1199 alt = [faref_leftflank[pos-left-1] + allele for allele in alt]
jpayne@69 1200 pos -= 1
jpayne@69 1201
jpayne@69 1202 __all__ = [
jpayne@69 1203 "VCF", "VCFRecord", ]