<s> """<STR_LIT>""" <EOL> from __future__ import division <EOL> from __future__ import print_function <EOL> from __future__ import unicode_literals <EOL> import datetime <EOL> import random <EOL> import hashlib <EOL> import re <EOL> import os <EOL> import glob <EOL> import pysam <EOL> import ga4gh . protocol as protocol <EOL> import ga4gh . exceptions as exceptions <EOL> import ga4gh . datamodel as datamodel <EOL> def isUnspecified ( str ) : <EOL> """<STR_LIT>""" <EOL> return str == "<STR_LIT>" or str is None <EOL> def convertVCFPhaseset ( vcfPhaseset ) : <EOL> """<STR_LIT>""" <EOL> if vcfPhaseset is not None and vcfPhaseset != "<STR_LIT:.>" : <EOL> phaseset = vcfPhaseset <EOL> else : <EOL> phaseset = "<STR_LIT:*>" <EOL> return phaseset <EOL> def convertVCFGenotype ( vcfGenotype , vcfPhaseset ) : <EOL> """<STR_LIT>""" <EOL> phaseset = None <EOL> if vcfGenotype is not None : <EOL> delim = "<STR_LIT:/>" <EOL> if "<STR_LIT:|>" in vcfGenotype : <EOL> delim = "<STR_LIT:|>" <EOL> phaseset = convertVCFPhaseset ( vcfPhaseset ) <EOL> if "<STR_LIT:.>" in vcfGenotype : <EOL> genotype = [ - <NUM_LIT:1> ] <EOL> else : <EOL> genotype = map ( int , vcfGenotype . split ( delim ) ) <EOL> else : <EOL> genotype = [ - <NUM_LIT:1> ] <EOL> return genotype , phaseset <EOL> class CallSet ( datamodel . DatamodelObject ) : <EOL> """<STR_LIT>""" <EOL> compoundIdClass = datamodel . CallSetCompoundId <EOL> def toProtocolElement ( self ) : <EOL> """<STR_LIT>""" <EOL> variantSet = self . getParentContainer ( ) <EOL> gaCallSet = protocol . CallSet ( ) <EOL> gaCallSet . created = variantSet . getCreationTime ( ) <EOL> gaCallSet . updated = variantSet . getUpdatedTime ( ) <EOL> gaCallSet . id = self . getId ( ) <EOL> gaCallSet . name = self . getLocalId ( ) <EOL> gaCallSet . sampleId = self . getLocalId ( ) <EOL> gaCallSet . variantSetIds = [ variantSet . getId ( ) ] <EOL> return gaCallSet <EOL> def getSampleName ( self ) : <EOL> """<STR_LIT>""" <EOL> return self . getLocalId ( ) <EOL> class AbstractVariantSet ( datamodel . DatamodelObject ) : <EOL> """<STR_LIT>""" <EOL> compoundIdClass = datamodel . VariantSetCompoundId <EOL> def __init__ ( self , parentContainer , localId ) : <EOL> super ( AbstractVariantSet , self ) . __init__ ( parentContainer , localId ) <EOL> self . _callSetIdMap = { } <EOL> self . _callSetNameMap = { } <EOL> self . _callSetIds = [ ] <EOL> self . _creationTime = None <EOL> self . _updatedTime = None <EOL> self . _referenceSetId = "<STR_LIT>" <EOL> def getCreationTime ( self ) : <EOL> """<STR_LIT>""" <EOL> return self . _creationTime <EOL> def getUpdatedTime ( self ) : <EOL> """<STR_LIT>""" <EOL> return self . _updatedTime <EOL> def addCallSet ( self , sampleName ) : <EOL> """<STR_LIT>""" <EOL> callSet = CallSet ( self , sampleName ) <EOL> callSetId = callSet . getId ( ) <EOL> self . _callSetIdMap [ callSetId ] = callSet <EOL> self . _callSetNameMap [ sampleName ] = callSet <EOL> self . _callSetIds . append ( callSetId ) <EOL> def getCallSets ( self ) : <EOL> """<STR_LIT>""" <EOL> return [ self . _callSetIdMap [ id_ ] for id_ in self . _callSetIds ] <EOL> def getNumCallSets ( self ) : <EOL> """<STR_LIT>""" <EOL> return len ( self . _callSetIds ) <EOL> def getCallSetByName ( self , name ) : <EOL> """<STR_LIT>""" <EOL> if name not in self . _callSetNameMap : <EOL> raise exceptions . CallSetNameNotFoundException ( name ) <EOL> return self . _callSetNameMap [ name ] <EOL> def getCallSetByIndex ( self , index ) : <EOL> """<STR_LIT>""" <EOL> return self . _callSetIdMap [ self . _callSetIds [ index ] ] <EOL> def getCallSet ( self , id_ ) : <EOL> """<STR_LIT>""" <EOL> if id_ not in self . _callSetIdMap : <EOL> raise exceptions . CallSetNotFoundException ( id_ ) <EOL> return self . _callSetIdMap [ id_ ] <EOL> def toProtocolElement ( self ) : <EOL> """<STR_LIT>""" <EOL> protocolElement = protocol . VariantSet ( ) <EOL> protocolElement . id = self . getId ( ) <EOL> protocolElement . datasetId = self . getParentContainer ( ) . getId ( ) <EOL> protocolElement . referenceSetId = self . _referenceSetId <EOL> protocolElement . metadata = self . getMetadata ( ) <EOL> protocolElement . name = self . getLocalId ( ) <EOL> return protocolElement <EOL> def getNumVariants ( self ) : <EOL> """<STR_LIT>""" <EOL> raise NotImplementedError ( ) <EOL> def _createGaVariant ( self ) : <EOL> """<STR_LIT>""" <EOL> ret = protocol . Variant ( ) <EOL> ret . created = self . _creationTime <EOL> ret . updated = self . _updatedTime <EOL> ret . variantSetId = self . getId ( ) <EOL> return ret <EOL> def _createGaVariantAnnotation ( self ) : <EOL> """<STR_LIT>""" <EOL> ret = protocol . VariantAnnotation ( ) <EOL> ret . created = self . _creationTime <EOL> ret . updated = self . _updatedTime <EOL> ret . variantAnnotationSetId = self . getId ( ) <EOL> return ret <EOL> def _createGaTranscriptEffect ( self ) : <EOL> """<STR_LIT>""" <EOL> ret = protocol . TranscriptEffect ( ) <EOL> ret . created = self . _creationTime <EOL> ret . updated = self . _updatedTime <EOL> return ret <EOL> def _createGaOntologyTermSo ( self ) : <EOL> """<STR_LIT>""" <EOL> ret = protocol . OntologyTerm ( ) <EOL> ret . ontologySource = "<STR_LIT>" <EOL> return ret <EOL> def _createGaAlleleLocation ( self ) : <EOL> """<STR_LIT>""" <EOL> ret = protocol . AlleleLocation ( ) <EOL> ret . created = self . _creationTime <EOL> ret . updated = self . _updatedTime <EOL> return ret <EOL> def getVariantId ( self , gaVariant ) : <EOL> """<STR_LIT>""" <EOL> md5 = self . hashVariant ( gaVariant ) <EOL> compoundId = datamodel . VariantCompoundId ( <EOL> self . getCompoundId ( ) , gaVariant . referenceName , <EOL> gaVariant . start , md5 ) <EOL> return str ( compoundId ) <EOL> def getCallSetId ( self , sampleName ) : <EOL> """<STR_LIT>""" <EOL> compoundId = datamodel . CallSetCompoundId ( <EOL> self . getCompoundId ( ) , sampleName ) <EOL> return str ( compoundId ) <EOL> @ classmethod <EOL> def hashVariant ( cls , gaVariant ) : <EOL> """<STR_LIT>""" <EOL> return hashlib . md5 ( <EOL> gaVariant . referenceBases + <EOL> str ( tuple ( gaVariant . alternateBases ) ) ) . hexdigest ( ) <EOL> class SimulatedVariantSet ( AbstractVariantSet ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( <EOL> self , parentContainer , localId , randomSeed = <NUM_LIT:1> , numCalls = <NUM_LIT:1> , <EOL> variantDensity = <NUM_LIT:1> ) : <EOL> super ( SimulatedVariantSet , self ) . __init__ ( parentContainer , localId ) <EOL> self . _randomSeed = randomSeed <EOL> self . _numCalls = numCalls <EOL> for j in range ( numCalls ) : <EOL> self . addCallSet ( "<STR_LIT>" . format ( j ) ) <EOL> self . _variantDensity = variantDensity <EOL> now = protocol . convertDatetime ( datetime . datetime . now ( ) ) <EOL> self . _creationTime = now <EOL> self . _updatedTime = now <EOL> def getNumVariants ( self ) : <EOL> return <NUM_LIT:0> <EOL> def getMetadata ( self ) : <EOL> ret = [ ] <EOL> return ret <EOL> def getVariant ( self , compoundId ) : <EOL> randomNumberGenerator = random . Random ( ) <EOL> start = int ( compoundId . start ) <EOL> randomNumberGenerator . seed ( self . _randomSeed + start ) <EOL> variant = self . generateVariant ( <EOL> compoundId . referenceName , start , randomNumberGenerator ) <EOL> return variant <EOL> def getVariants ( self , referenceName , startPosition , endPosition , <EOL> callSetIds = None ) : <EOL> randomNumberGenerator = random . Random ( ) <EOL> randomNumberGenerator . seed ( self . _randomSeed ) <EOL> i = startPosition <EOL> while i < endPosition : <EOL> if randomNumberGenerator . random ( ) < self . _variantDensity : <EOL> randomNumberGenerator . seed ( self . _randomSeed + i ) <EOL> yield self . generateVariant ( <EOL> referenceName , i , randomNumberGenerator ) <EOL> i += <NUM_LIT:1> <EOL> def generateVariant ( self , referenceName , position , randomNumberGenerator ) : <EOL> """<STR_LIT>""" <EOL> variant = self . _createGaVariant ( ) <EOL> variant . names = [ ] <EOL> variant . referenceName = referenceName <EOL> variant . start = position <EOL> variant . end = position + <NUM_LIT:1> <EOL> bases = [ "<STR_LIT:A>" , "<STR_LIT:C>" , "<STR_LIT>" , "<STR_LIT:T>" ] <EOL> ref = randomNumberGenerator . choice ( bases ) <EOL> variant . referenceBases = ref <EOL> alt = randomNumberGenerator . choice ( <EOL> [ base for base in bases if base != ref ] ) <EOL> variant . alternateBases = [ alt ] <EOL> variant . calls = [ ] <EOL> for callSet in self . getCallSets ( ) : <EOL> call = protocol . Call ( ) <EOL> call . callSetId = callSet . getId ( ) <EOL> randomChoice = randomNumberGenerator . choice ( <EOL> [ [ <NUM_LIT:0> , <NUM_LIT:1> ] , [ <NUM_LIT:1> , <NUM_LIT:0> ] , [ <NUM_LIT:1> , <NUM_LIT:1> ] ] ) <EOL> call . genotype = randomChoice <EOL> call . genotypeLikelihood = [ - <NUM_LIT:100> , - <NUM_LIT:100> , - <NUM_LIT:100> ] <EOL> variant . calls . append ( call ) <EOL> variant . id = self . getVariantId ( variant ) <EOL> return variant <EOL> class SimulatedVariantAnnotationSet ( AbstractVariantSet ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( <EOL> self , parentContainer , localId , variantSet , randomSeed = <NUM_LIT:1> ) : <EOL> super ( SimulatedVariantAnnotationSet , self ) . __init__ ( <EOL> parentContainer , localId ) <EOL> self . _variantSet = variantSet <EOL> self . _variantSetId = str ( variantSet . getCompoundId ( ) ) <EOL> self . _randomSeed = randomSeed <EOL> self . _compoundId = datamodel . VariantAnnotationSetCompoundId ( <EOL> self . getCompoundId ( ) , '<STR_LIT>' ) <EOL> self . _creationTime = datetime . datetime . now ( ) . isoformat ( ) + "<STR_LIT>" <EOL> self . _updatedTime = datetime . datetime . now ( ) . isoformat ( ) + "<STR_LIT>" <EOL> def getAnalysis ( self ) : <EOL> analysis = protocol . Analysis ( ) <EOL> analysis . createDateTime = self . _creationTime <EOL> analysis . updateDateTime = self . _updatedTime <EOL> analysis . software . append ( "<STR_LIT>" ) <EOL> analysis . name = "<STR_LIT:name>" <EOL> analysis . description = "<STR_LIT:description>" <EOL> analysis . id = str ( datamodel . VariantAnnotationSetAnalysisCompoundId ( <EOL> self . _compoundId , "<STR_LIT>" ) ) <EOL> return analysis <EOL> def toProtocolElement ( self ) : <EOL> """<STR_LIT>""" <EOL> protocolElement = protocol . VariantAnnotationSet ( ) <EOL> protocolElement . variantSetId = self . _variantSet . getId ( ) <EOL> protocolElement . id = self . getId ( ) <EOL> protocolElement . name = self . getLocalId ( ) <EOL> protocolElement . analysis = self . getAnalysis ( ) <EOL> return protocolElement <EOL> def getNumVariantAnnotations ( self ) : <EOL> return <NUM_LIT:0> <EOL> def getVariantAnnotation ( self , variant , randomNumberGenerator ) : <EOL> ann = self . generateVariantAnnotation ( <EOL> variant , randomNumberGenerator ) <EOL> return ann <EOL> def getVariantAnnotations ( self , referenceName , start , end ) : <EOL> for variant in self . _variantSet . getVariants ( referenceName , start , end ) : <EOL> yield self . generateVariantAnnotation ( variant ) <EOL> def generateVariantAnnotation ( self , variant ) : <EOL> """<STR_LIT>""" <EOL> seed = self . _randomSeed + variant . start + variant . end <EOL> randomNumberGenerator = random . Random ( ) <EOL> randomNumberGenerator . seed ( seed ) <EOL> ann = protocol . VariantAnnotation ( ) <EOL> ann . variantAnnotationSetId = str ( self . getCompoundId ( ) ) <EOL> ann . variantId = variant . id <EOL> ann . start = variant . start <EOL> ann . end = variant . end <EOL> ann . createDateTime = self . _creationTime <EOL> ann . transcriptEffects = [ ] <EOL> for base in variant . alternateBases * ( <EOL> randomNumberGenerator . randint ( <NUM_LIT:0> , <NUM_LIT:5> ) ) : <EOL> ann . transcriptEffects . append ( self . generateTranscriptEffect ( <EOL> ann , base , randomNumberGenerator ) ) <EOL> ann . id = self . getVariantAnnotationId ( variant , ann ) <EOL> return ann <EOL> def getTranscriptEffectId ( self , gaTranscriptEffect ) : <EOL> effs = [ eff . term for eff in gaTranscriptEffect . effects ] <EOL> return hashlib . md5 ( <EOL> "<STR_LIT>" . format ( <EOL> gaTranscriptEffect . alternateBases , <EOL> gaTranscriptEffect . featureId , <EOL> effs , gaTranscriptEffect . hgvsAnnotation ) <EOL> ) . hexdigest ( ) <EOL> def getVariantAnnotationId ( self , gaVariant , gaAnnotation ) : <EOL> md5 = self . hashVariantAnnotation ( gaVariant , gaAnnotation ) <EOL> compoundId = datamodel . VariantAnnotationCompoundId ( <EOL> self . getCompoundId ( ) , gaVariant . referenceName , <EOL> gaVariant . start , md5 ) <EOL> return str ( compoundId ) <EOL> def _addTranscriptEffectLocations ( self , effect , ann ) : <EOL> effect . hgvsAnnotation = protocol . HGVSAnnotation ( ) <EOL> effect . hgvsAnnotation . genomic = str ( ann . start ) <EOL> effect . hgvsAnnotation . transcript = str ( ann . start ) <EOL> effect . hgvsAnnotation . protein = str ( ann . start ) <EOL> effect . proteinLocation = self . _createGaAlleleLocation ( ) <EOL> effect . proteinLocation . start = ann . start <EOL> effect . CDSLocation = self . _createGaAlleleLocation ( ) <EOL> effect . CDSLocation . start = ann . start <EOL> effect . cDNALocation = self . _createGaAlleleLocation ( ) <EOL> effect . cDNALocation . start = ann . start <EOL> return effect <EOL> def _addTranscriptEffectId ( self , effect ) : <EOL> effect . id = str ( self . getTranscriptEffectId ( effect ) ) <EOL> return effect <EOL> def _getRandomOntologyTerm ( self , randomNumberGenerator ) : <EOL> ontologyTuples = [ ( "<STR_LIT>" , "<STR_LIT>" ) , <EOL> ( "<STR_LIT>" , "<STR_LIT>" ) ] <EOL> term = protocol . OntologyTerm ( ) <EOL> ontologyTuple = randomNumberGenerator . choice ( ontologyTuples ) <EOL> term . term , term . id = ontologyTuple [ <NUM_LIT:0> ] , ontologyTuple [ <NUM_LIT:1> ] <EOL> term . sourceName = "<STR_LIT>" <EOL> term . sourceVersion = "<STR_LIT:0>" <EOL> return term <EOL> def _addTranscriptEffectOntologyTerm ( self , effect , randomNumberGenerator ) : <EOL> effect . effects . append ( <EOL> self . _getRandomOntologyTerm ( randomNumberGenerator ) ) <EOL> return effect <EOL> def _generateAnalysisResult ( self , effect , ann , randomNumberGenerator ) : <EOL> analysisResult = protocol . AnalysisResult ( ) <EOL> analysisResult . analysisId = "<STR_LIT>" <EOL> analysisResult . result = "<STR_LIT>" <EOL> analysisResult . score = randomNumberGenerator . randint ( <NUM_LIT:0> , <NUM_LIT:100> ) <EOL> return analysisResult <EOL> def _addAnalysisResult ( self , effect , ann , randomNumberGenerator ) : <EOL> effect . analysisResults . append ( <EOL> self . _generateAnalysisResult ( <EOL> effect , ann , randomNumberGenerator ) ) <EOL> return effect <EOL> def generateTranscriptEffect ( self , ann , alts , randomNumberGenerator ) : <EOL> effect = self . _createGaTranscriptEffect ( ) <EOL> effect . alternateBases = alts <EOL> effect . effects = [ ] <EOL> effect . analysisResults = [ ] <EOL> effect . featureId = "<STR_LIT>" <EOL> effect = self . _addTranscriptEffectLocations ( effect , ann ) <EOL> effect = self . _addTranscriptEffectOntologyTerm ( <EOL> effect , randomNumberGenerator ) <EOL> effect = self . _addTranscriptEffectOntologyTerm ( <EOL> effect , randomNumberGenerator ) <EOL> effect = self . _addTranscriptEffectId ( effect ) <EOL> effect = self . _addAnalysisResult ( effect , ann , randomNumberGenerator ) <EOL> return effect <EOL> def hashVariantAnnotation ( cls , gaVariant , gaVariantAnnotation ) : <EOL> """<STR_LIT>""" <EOL> treffs = [ treff . id for treff in gaVariantAnnotation . transcriptEffects ] <EOL> return hashlib . md5 ( <EOL> "<STR_LIT>" . format ( <EOL> gaVariant . referenceBases , tuple ( gaVariant . alternateBases ) , <EOL> treffs ) <EOL> ) . hexdigest ( ) <EOL> def _encodeValue ( value ) : <EOL> if isinstance ( value , ( list , tuple ) ) : <EOL> return [ str ( v ) for v in value ] <EOL> else : <EOL> return [ str ( value ) ] <EOL> _nothing = object ( ) <EOL> def isEmptyIter ( it ) : <EOL> """<STR_LIT>""" <EOL> return next ( it , _nothing ) is _nothing <EOL> class HtslibVariantSet ( datamodel . PysamDatamodelMixin , AbstractVariantSet ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , parentContainer , localId , dataDir , dataRepository ) : <EOL> super ( HtslibVariantSet , self ) . __init__ ( parentContainer , localId ) <EOL> self . _dataDir = dataDir <EOL> self . _setAccessTimes ( dataDir ) <EOL> self . _chromFileMap = { } <EOL> self . _metadata = None <EOL> self . _patterns = [ '<STR_LIT>' , '<STR_LIT>' ] <EOL> self . _scanDataFiles ( self . _dataDir , self . _patterns ) <EOL> def _updateMetadata ( self , variantFile ) : <EOL> """<STR_LIT>""" <EOL> metadata = self . _getMetadataFromVcf ( variantFile ) <EOL> if self . _metadata is None : <EOL> self . _metadata = metadata <EOL> def _checkMetadata ( self , variantFile ) : <EOL> """<STR_LIT>""" <EOL> metadata = self . _getMetadataFromVcf ( variantFile ) <EOL> if self . _metadata is not None and self . _metadata != metadata : <EOL> raise exceptions . InconsistentMetaDataException ( <EOL> variantFile . filename ) <EOL> def _checkCallSetIds ( self , variantFile ) : <EOL> """<STR_LIT>""" <EOL> if len ( self . _callSetIdMap ) > <NUM_LIT:0> : <EOL> callSetIds = set ( [ <EOL> self . getCallSetId ( sample ) <EOL> for sample in variantFile . header . samples ] ) <EOL> if callSetIds != set ( self . _callSetIdMap . keys ( ) ) : <EOL> raise exceptions . InconsistentCallSetIdException ( <EOL> variantFile . filename ) <EOL> def checkConsistency ( self ) : <EOL> """<STR_LIT>""" <EOL> filenames = self . _getDataFilenames ( self . _dataDir , self . _patterns ) <EOL> for filename in filenames : <EOL> varFile = self . openFile ( filename ) <EOL> for chrom in varFile . index : <EOL> chrom , _ , _ = self . sanitizeVariantFileFetch ( chrom ) <EOL> if not isEmptyIter ( varFile . fetch ( chrom ) ) : <EOL> self . _checkMetadata ( varFile ) <EOL> self . _checkCallSetIds ( varFile ) <EOL> varFile . close ( ) <EOL> def getNumVariants ( self ) : <EOL> """<STR_LIT>""" <EOL> return <NUM_LIT:0> <EOL> def _updateCallSetIds ( self , variantFile ) : <EOL> """<STR_LIT>""" <EOL> if len ( self . _callSetIdMap ) == <NUM_LIT:0> : <EOL> for sample in variantFile . header . samples : <EOL> self . addCallSet ( sample ) <EOL> def openFile ( self , filename ) : <EOL> return pysam . VariantFile ( filename ) <EOL> def _addDataFile ( self , filename ) : <EOL> varFile = self . openFile ( filename ) <EOL> if varFile . index is None : <EOL> raise exceptions . NotIndexedException ( filename ) <EOL> for chrom in varFile . index : <EOL> chrom , _ , _ = self . sanitizeVariantFileFetch ( chrom ) <EOL> if not isEmptyIter ( varFile . fetch ( chrom ) ) : <EOL> if chrom in self . _chromFileMap : <EOL> raise exceptions . OverlappingVcfException ( filename , chrom ) <EOL> self . _updateMetadata ( varFile ) <EOL> self . _updateCallSetIds ( varFile ) <EOL> self . _chromFileMap [ chrom ] = filename <EOL> varFile . close ( ) <EOL> def _convertGaCall ( self , recordId , name , pysamCall , genotypeData ) : <EOL> compoundId = self . getCallSetId ( name ) <EOL> callSet = self . getCallSet ( compoundId ) <EOL> call = protocol . Call ( ) <EOL> call . callSetId = callSet . getId ( ) <EOL> call . callSetName = callSet . getSampleName ( ) <EOL> call . sampleId = callSet . getSampleName ( ) <EOL> call . phaseset = None <EOL> call . genotype , call . phaseset = convertVCFGenotype ( <EOL> genotypeData , call . phaseset ) <EOL> call . genotypeLikelihood = [ ] <EOL> for key , value in pysamCall . iteritems ( ) : <EOL> if key == '<STR_LIT>' and value is not None : <EOL> call . genotypeLikelihood = list ( value ) <EOL> elif key != '<STR_LIT>' : <EOL> call . info [ key ] = _encodeValue ( value ) <EOL> return call <EOL> def convertVariant ( self , record , callSetIds ) : <EOL> """<STR_LIT>""" <EOL> variant = self . _createGaVariant ( ) <EOL> variant . referenceName = record . contig <EOL> if record . id is not None : <EOL> variant . names = record . id . split ( '<STR_LIT:;>' ) <EOL> variant . start = record . start <EOL> variant . end = record . stop <EOL> variant . referenceBases = record . ref <EOL> if record . alts is not None : <EOL> variant . alternateBases = list ( record . alts ) <EOL> for key , value in record . info . iteritems ( ) : <EOL> if value is not None : <EOL> if isinstance ( value , str ) : <EOL> value = value . split ( '<STR_LIT:U+002C>' ) <EOL> variant . info [ key ] = _encodeValue ( value ) <EOL> sampleData = record . __str__ ( ) . split ( ) [ <NUM_LIT:9> : ] <EOL> variant . calls = [ ] <EOL> sampleIterator = <NUM_LIT:0> <EOL> if callSetIds is not None : <EOL> for name , call in record . samples . iteritems ( ) : <EOL> if self . getCallSetId ( name ) in callSetIds : <EOL> genotypeData = sampleData [ sampleIterator ] . split ( <EOL> "<STR_LIT::>" ) [ <NUM_LIT:0> ] <EOL> variant . calls . append ( self . _convertGaCall ( <EOL> record . id , name , call , genotypeData ) ) <EOL> sampleIterator += <NUM_LIT:1> <EOL> variant . id = self . getVariantId ( variant ) <EOL> return variant <EOL> def getVariant ( self , compoundId ) : <EOL> if compoundId . referenceName in self . _chromFileMap : <EOL> varFileName = self . _chromFileMap [ compoundId . referenceName ] <EOL> else : <EOL> raise exceptions . ObjectNotFoundException ( compoundId ) <EOL> start = int ( compoundId . start ) <EOL> referenceName , startPosition , endPosition = self . sanitizeVariantFileFetch ( <EOL> compoundId . referenceName , start , start + <NUM_LIT:1> ) <EOL> cursor = self . getFileHandle ( varFileName ) . fetch ( <EOL> referenceName , startPosition , endPosition ) <EOL> for record in cursor : <EOL> variant = self . convertVariant ( record , self . _callSetIds ) <EOL> if ( record . start == start and <EOL> compoundId . md5 == self . hashVariant ( variant ) ) : <EOL> return variant <EOL> elif record . start > start : <EOL> raise exceptions . ObjectNotFoundException ( ) <EOL> raise exceptions . ObjectNotFoundException ( compoundId ) <EOL> def getVariants ( self , referenceName , startPosition , endPosition , <EOL> callSetIds = None ) : <EOL> """<STR_LIT>""" <EOL> if callSetIds is None : <EOL> callSetIds = self . _callSetIds <EOL> else : <EOL> for callSetId in callSetIds : <EOL> if callSetId not in self . _callSetIds : <EOL> raise exceptions . CallSetNotInVariantSetException ( <EOL> callSetId , self . getId ( ) ) <EOL> if referenceName in self . _chromFileMap : <EOL> varFileName = self . _chromFileMap [ referenceName ] <EOL> referenceName , startPosition , endPosition = self . sanitizeVariantFileFetch ( <EOL> referenceName , startPosition , endPosition ) <EOL> cursor = self . getFileHandle ( varFileName ) . fetch ( <EOL> referenceName , startPosition , endPosition ) <EOL> for record in cursor : <EOL> yield self . convertVariant ( record , callSetIds ) <EOL> def getMetadata ( self ) : <EOL> return self . _metadata <EOL> def getMetadataId ( self , metadata ) : <EOL> """<STR_LIT>""" <EOL> return str ( datamodel . VariantSetMetadataCompoundId ( <EOL> self . getCompoundId ( ) , '<STR_LIT>' + metadata . key ) ) <EOL> def _getMetadataFromVcf ( self , varFile ) : <EOL> def buildMetadata ( <EOL> key , type_ = "<STR_LIT>" , number = "<STR_LIT:1>" , value = "<STR_LIT>" , id_ = "<STR_LIT>" , <EOL> description = "<STR_LIT>" ) : <EOL> metadata = protocol . VariantSetMetadata ( ) <EOL> metadata . key = key <EOL> metadata . value = value <EOL> metadata . type = type_ <EOL> metadata . number = number <EOL> metadata . description = description <EOL> if id_ == '<STR_LIT>' : <EOL> id_ = self . getMetadataId ( metadata ) <EOL> metadata . id = id_ <EOL> return metadata <EOL> ret = [ ] <EOL> header = varFile . header <EOL> ret . append ( buildMetadata ( key = "<STR_LIT:version>" , value = header . version ) ) <EOL> formats = header . formats . items ( ) <EOL> infos = header . info . items ( ) <EOL> for prefix , content in [ ( "<STR_LIT>" , formats ) , ( "<STR_LIT>" , infos ) ] : <EOL> for contentKey , value in content : <EOL> description = value . description . strip ( '<STR_LIT:">' ) <EOL> key = "<STR_LIT>" . format ( prefix , value . name ) <EOL> if key != "<STR_LIT>" : <EOL> ret . append ( buildMetadata ( <EOL> key = key , type_ = value . type , <EOL> number = "<STR_LIT:{}>" . format ( value . number ) , <EOL> description = description ) ) <EOL> return ret <EOL> def isAnnotated ( self , path ) : <EOL> return self . hasAnnField ( path ) or self . hasConsequenceField ( path ) <EOL> def hasAnnField ( self , path ) : <EOL> return '<STR_LIT>' in self . _getHeaderItems ( path ) <EOL> def hasConsequenceField ( self , path ) : <EOL> return '<STR_LIT>' in self . _getHeaderItems ( path ) <EOL> def _getHeaderItems ( self , path ) : <EOL> pysamreader = self . openFile ( <EOL> glob . glob ( os . path . join ( path , "<STR_LIT>" ) ) [ <NUM_LIT:0> ] ) <EOL> return [ x [ <NUM_LIT:0> ] for x in pysamreader . header . info . items ( ) ] <EOL> class HtslibVariantAnnotationSet ( HtslibVariantSet ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , parentContainer , localId , dataDir , <EOL> backend , variantSet ) : <EOL> super ( HtslibVariantAnnotationSet , self ) . __init__ ( <EOL> parentContainer , localId , dataDir , backend ) <EOL> self . compoundIdClass = datamodel . VariantAnnotationSetCompoundId <EOL> self . _variantSetId = variantSet . getCompoundId ( ) <EOL> self . _variantSet = variantSet <EOL> self . _compoundId = datamodel . VariantAnnotationSetCompoundId ( <EOL> self . getCompoundId ( ) , '<STR_LIT>' ) <EOL> self . _sequenceOntology = backend . getOntologyMap ( '<STR_LIT>' ) <EOL> self . _creationTime = datetime . datetime . now ( ) . isoformat ( ) + "<STR_LIT>" <EOL> self . _updatedTime = datetime . datetime . now ( ) . isoformat ( ) + "<STR_LIT>" <EOL> self . _isVep = "<STR_LIT>" in self . toProtocolElement ( ) . analysis . info <EOL> self . _isCsq = self . hasConsequenceField ( dataDir ) <EOL> self . _annotationCreatedDateTime = self . _creationTime <EOL> for r in self . getMetadata ( ) . records : <EOL> if r . key == "<STR_LIT>" : <EOL> self . _annotationCreatedDateTime = datetime . datetime . strptime ( <EOL> r . value , "<STR_LIT>" ) . isoformat ( ) + "<STR_LIT>" <EOL> def getVariantAnnotations ( self , referenceName , startPosition , endPosition ) : <EOL> """<STR_LIT>""" <EOL> if referenceName in self . _chromFileMap : <EOL> varFileName = self . _chromFileMap [ referenceName ] <EOL> referenceName , startPosition , endPosition = self . sanitizeVariantFileFetch ( <EOL> referenceName , startPosition , endPosition ) <EOL> cursor = self . getFileHandle ( varFileName ) . fetch ( <EOL> referenceName , startPosition , endPosition ) <EOL> transcriptConverter = self . convertTranscriptEffectSnpEff <EOL> if self . _isVep and not self . _isCsq : <EOL> transcriptConverter = self . convertTranscriptEffectVEP <EOL> elif self . _isCsq : <EOL> transcriptConverter = self . convertTranscriptEffectCSQ <EOL> for record in cursor : <EOL> yield self . convertVariantAnnotation ( <EOL> record , transcriptConverter ) <EOL> def convertLocation ( self , pos ) : <EOL> """<STR_LIT>""" <EOL> if isUnspecified ( pos ) : <EOL> return None <EOL> coordLen = pos . split ( '<STR_LIT:/>' ) <EOL> if len ( coordLen ) > <NUM_LIT:1> : <EOL> allLoc = self . _createGaAlleleLocation ( ) <EOL> allLoc . start = int ( coordLen [ <NUM_LIT:0> ] ) - <NUM_LIT:1> <EOL> return allLoc <EOL> return None <EOL> def convertLocationHgvsC ( self , hgvsc ) : <EOL> """<STR_LIT>""" <EOL> if isUnspecified ( hgvsc ) : <EOL> return None <EOL> match = re . match ( "<STR_LIT>" , hgvsc ) <EOL> if match : <EOL> pos = int ( match . group ( <NUM_LIT:1> ) ) <EOL> if pos > <NUM_LIT:0> : <EOL> allLoc = self . _createGaAlleleLocation ( ) <EOL> allLoc . start = pos - <NUM_LIT:1> <EOL> allLoc . referenceSequence = match . group ( <NUM_LIT:2> ) <EOL> allLoc . alternateSequence = match . group ( <NUM_LIT:3> ) <EOL> return allLoc <EOL> return None <EOL> def convertLocationHgvsP ( self , hgvsp ) : <EOL> """<STR_LIT>""" <EOL> if isUnspecified ( hgvsp ) : <EOL> return None <EOL> match = re . match ( "<STR_LIT>" , hgvsp , flags = re . UNICODE ) <EOL> if match is not None : <EOL> allLoc = self . _createGaAlleleLocation ( ) <EOL> allLoc . referenceSequence = match . group ( <NUM_LIT:1> ) <EOL> allLoc . start = int ( match . group ( <NUM_LIT:2> ) ) - <NUM_LIT:1> <EOL> allLoc . alternateSequence = match . group ( <NUM_LIT:3> ) <EOL> return allLoc <EOL> return None <EOL> def addCDSLocation ( self , effect , cdnaPos ) : <EOL> hgvsC = effect . hgvsAnnotation . transcript <EOL> if not isUnspecified ( hgvsC ) : <EOL> effect . CDSLocation = self . convertLocationHgvsC ( hgvsC ) <EOL> if effect . CDSLocation is None : <EOL> effect . CDSLocation = self . convertLocation ( cdnaPos ) <EOL> else : <EOL> effect . CDSLocation . alternateSequence = None <EOL> effect . CDSLocation . referenceSequence = None <EOL> def addProteinLocation ( self , effect , protPos ) : <EOL> hgvsP = effect . hgvsAnnotation . protein <EOL> if not isUnspecified ( hgvsP ) : <EOL> effect . proteinLocation = self . convertLocationHgvsP ( hgvsP ) <EOL> if effect . proteinLocation is None : <EOL> effect . proteinLocation = self . convertLocation ( protPos ) <EOL> def addCDNALocation ( self , effect , cdnaPos ) : <EOL> hgvsC = effect . hgvsAnnotation . transcript <EOL> effect . cDNALocation = self . convertLocation ( cdnaPos ) <EOL> if self . convertLocationHgvsC ( hgvsC ) : <EOL> effect . cDNALocation . alternateSequence = self . convertLocationHgvsC ( hgvsC ) . alternateSequence <EOL> effect . cDNALocation . referenceSequence = self . convertLocationHgvsC ( hgvsC ) . referenceSequence <EOL> def addLocations ( self , effect , protPos , cdnaPos ) : <EOL> """<STR_LIT>""" <EOL> self . addCDSLocation ( effect , cdnaPos ) <EOL> self . addCDNALocation ( effect , cdnaPos ) <EOL> self . addProteinLocation ( effect , protPos ) <EOL> return effect <EOL> def convertTranscriptEffectCSQ ( self , annStr , hgvsG ) : <EOL> """<STR_LIT>""" <EOL> ( alt , gene , featureId , featureType , effects , cdnaPos , <EOL> cdsPos , protPos , aminos , codons , existingVar , distance , <EOL> strand , sift , polyPhen , motifName , motifPos , <EOL> highInfPos , motifScoreChange ) = annStr . split ( '<STR_LIT:|>' ) <EOL> terms = effects . split ( "<STR_LIT:&>" ) <EOL> transcriptEffects = [ ] <EOL> for term in terms : <EOL> transcriptEffects . append ( <EOL> self . _createCsqTranscriptEffect ( <EOL> alt , term , protPos , <EOL> cdnaPos , featureId ) ) <EOL> return transcriptEffects <EOL> def _createCsqTranscriptEffect ( <EOL> self , alt , term , protPos , cdnaPos , featureId ) : <EOL> effect = self . _createGaTranscriptEffect ( ) <EOL> effect . alternateBases = alt <EOL> effect . effects = self . convertSeqOntology ( term ) <EOL> effect . featureId = featureId <EOL> effect . hgvsAnnotation = protocol . HGVSAnnotation ( ) <EOL> effect . hgvsAnnotation . genomic = None <EOL> effect . hgvsAnnotation . transcript = None <EOL> effect . hgvsAnnotation . protein = None <EOL> self . addLocations ( effect , protPos , cdnaPos ) <EOL> effect . id = self . getTranscriptEffectId ( effect ) <EOL> effect . analysisResults = [ ] <EOL> return effect <EOL> def convertTranscriptEffectVEP ( self , annStr , hgvsG ) : <EOL> """<STR_LIT>""" <EOL> effect = self . _createGaTranscriptEffect ( ) <EOL> ( alt , effects , impact , symbol , geneName , featureType , <EOL> featureId , trBiotype , exon , intron , hgvsC , hgvsP , <EOL> cdnaPos , cdsPos , protPos , aminos , codons , <EOL> existingVar , distance , strand , symbolSource , <EOL> hgncId , hgvsOffset ) = annStr . split ( '<STR_LIT:|>' ) <EOL> effect . alternateBases = alt <EOL> effect . effects = self . convertSeqOntology ( effects ) <EOL> effect . featureId = featureId <EOL> effect . hgvsAnnotation = protocol . HGVSAnnotation ( ) <EOL> effect . hgvsAnnotation . genomic = hgvsG <EOL> effect . hgvsAnnotation . transcript = hgvsC <EOL> effect . hgvsAnnotation . protein = hgvsP <EOL> self . addLocations ( effect , protPos , cdnaPos ) <EOL> effect . id = self . getTranscriptEffectId ( effect ) <EOL> effect . analysisResults = [ ] <EOL> return effect <EOL> def convertTranscriptEffectSnpEff ( self , annStr , hgvsG ) : <EOL> """<STR_LIT>""" <EOL> effect = self . _createGaTranscriptEffect ( ) <EOL> ( alt , effects , impact , geneName , geneId , featureType , <EOL> featureId , trBiotype , rank , hgvsC , hgvsP , cdnaPos , <EOL> cdsPos , protPos , distance , errsWarns ) = annStr . split ( '<STR_LIT:|>' ) <EOL> effect . alternateBases = alt <EOL> effect . effects = self . convertSeqOntology ( effects ) <EOL> effect . featureId = featureId <EOL> effect . hgvsAnnotation = protocol . HGVSAnnotation ( ) <EOL> effect . hgvsAnnotation . genomic = hgvsG <EOL> effect . hgvsAnnotation . transcript = hgvsC <EOL> effect . hgvsAnnotation . protein = hgvsP <EOL> self . addLocations ( effect , protPos , cdnaPos ) <EOL> effect . id = self . getTranscriptEffectId ( effect ) <EOL> effect . analysisResults = [ ] <EOL> return effect <EOL> def convertSeqOntology ( self , seqOntStr ) : <EOL> """<STR_LIT>""" <EOL> seqOntTerms = seqOntStr . split ( '<STR_LIT:&>' ) <EOL> soTerms = [ ] <EOL> for soName in seqOntTerms : <EOL> so = self . _createGaOntologyTermSo ( ) <EOL> so . term = soName <EOL> if self . _sequenceOntology is not None : <EOL> so . id = self . _sequenceOntology . getId ( soName , "<STR_LIT>" ) <EOL> soTerms . append ( so ) <EOL> return soTerms <EOL> def convertVariantAnnotation ( self , record , transcriptConverter ) : <EOL> """<STR_LIT>""" <EOL> variant = self . convertVariant ( record , None ) <EOL> annotation = self . _createGaVariantAnnotation ( ) <EOL> annotation . start = variant . start <EOL> annotation . end = variant . end <EOL> annotation . createDateTime = self . _annotationCreatedDateTime <EOL> annotation . variantId = variant . id <EOL> transcriptEffects = [ ] <EOL> hgvsG = record . info . get ( '<STR_LIT>' . encode ( ) ) <EOL> if transcriptConverter != self . convertTranscriptEffectCSQ : <EOL> annotations = record . info . get ( '<STR_LIT>' . encode ( ) ) <EOL> transcriptEffects = self . _convertAnnotations ( <EOL> annotations , variant , hgvsG , transcriptConverter ) <EOL> else : <EOL> annotations = record . info . get ( '<STR_LIT>' . encode ( ) ) <EOL> transcriptEffects = [ ] <EOL> for ann in annotations : <EOL> transcriptEffects . extend ( <EOL> self . convertTranscriptEffectCSQ ( <EOL> ann , hgvsG ) ) <EOL> annotation . transcriptEffects = transcriptEffects <EOL> annotation . id = self . getVariantAnnotationId ( variant , annotation ) <EOL> return annotation <EOL> def _convertAnnotations ( <EOL> self , annotations , variant , hgvsG , transcriptConverter ) : <EOL> transcriptEffects = [ ] <EOL> if annotations is not None : <EOL> for index , ann in enumerate ( annotations ) : <EOL> altshgvsG = "<STR_LIT>" <EOL> if hgvsG is not None : <EOL> altshgvsG = hgvsG [ index % len ( variant . alternateBases ) ] <EOL> transcriptEffects . append ( <EOL> transcriptConverter ( ann , altshgvsG ) ) <EOL> return transcriptEffects <EOL> def getVariantAnnotationId ( self , gaVariant , gaAnnotation ) : <EOL> """<STR_LIT>""" <EOL> md5 = self . hashVariantAnnotation ( gaVariant , gaAnnotation ) <EOL> compoundId = datamodel . VariantAnnotationCompoundId ( <EOL> self . getCompoundId ( ) , gaVariant . referenceName , <EOL> gaVariant . start , md5 ) <EOL> return str ( compoundId ) <EOL> def getVariantId ( self , gaVariant ) : <EOL> """<STR_LIT>""" <EOL> md5 = self . hashVariant ( gaVariant ) <EOL> compoundId = datamodel . VariantCompoundId ( <EOL> self . _variantSetId , gaVariant . referenceName , <EOL> gaVariant . start , md5 ) <EOL> return str ( compoundId ) <EOL> def getAnalysis ( self ) : <EOL> """<STR_LIT>""" <EOL> metadata = self . getMetadata ( ) <EOL> analysis = protocol . Analysis ( ) <EOL> formats = metadata . formats . items ( ) <EOL> infos = metadata . info . items ( ) <EOL> for prefix , content in [ ( "<STR_LIT>" , formats ) , ( "<STR_LIT>" , infos ) ] : <EOL> for contentKey , value in content : <EOL> key = "<STR_LIT>" . format ( prefix , value . name ) <EOL> if key not in analysis . info : <EOL> analysis . info [ key ] = [ ] <EOL> if value . description is not None : <EOL> analysis . info [ key ] . append ( value . description ) <EOL> analysis . createDateTime = self . _creationTime <EOL> analysis . updateDateTime = self . _updatedTime <EOL> for r in metadata . records : <EOL> if r . value is not None : <EOL> if r . key not in analysis . info : <EOL> analysis . info [ r . key ] = [ ] <EOL> analysis . info [ r . key ] . append ( str ( r . value ) ) <EOL> if r . key == "<STR_LIT>" : <EOL> analysis . createDateTime = datetime . datetime . strptime ( <EOL> r . value , "<STR_LIT>" ) . isoformat ( ) + "<STR_LIT>" <EOL> if r . key == "<STR_LIT>" : <EOL> analysis . software . append ( r . value ) <EOL> if r . key == "<STR_LIT:name>" : <EOL> analysis . name = r . value <EOL> if r . key == "<STR_LIT:description>" : <EOL> analysis . description = r . value <EOL> analysis . id = str ( datamodel . VariantAnnotationSetAnalysisCompoundId ( <EOL> self . _compoundId , "<STR_LIT>" ) ) <EOL> return analysis <EOL> def _getMetadataFromVcf ( self , varFile ) : <EOL> header = varFile . header <EOL> return header <EOL> def toProtocolElement ( self ) : <EOL> """<STR_LIT>""" <EOL> protocolElement = protocol . VariantAnnotationSet ( ) <EOL> protocolElement . id = self . getId ( ) <EOL> protocolElement . variantSetId = str ( self . _variantSetId ) <EOL> protocolElement . name = self . getLocalId ( ) <EOL> protocolElement . analysis = self . getAnalysis ( ) <EOL> return protocolElement <EOL> def getTranscriptEffectId ( self , gaTranscriptEffect ) : <EOL> effs = [ eff . term for eff in gaTranscriptEffect . effects ] <EOL> return hashlib . md5 ( <EOL> "<STR_LIT>" . format ( <EOL> gaTranscriptEffect . alternateBases , <EOL> gaTranscriptEffect . featureId , <EOL> effs , gaTranscriptEffect . hgvsAnnotation ) <EOL> ) . hexdigest ( ) <EOL> def hashVariantAnnotation ( cls , gaVariant , gaVariantAnnotation ) : <EOL> """<STR_LIT>""" <EOL> treffs = [ treff . id for treff in gaVariantAnnotation . transcriptEffects ] <EOL> return hashlib . md5 ( <EOL> "<STR_LIT>" . format ( <EOL> gaVariant . referenceBases , tuple ( gaVariant . alternateBases ) , <EOL> treffs ) <EOL> ) . hexdigest ( ) </s>
