<s> """<STR_LIT>""" <EOL> import sys <EOL> import collections <EOL> import pysam <EOL> import CGAT . Experiment as E <EOL> import CGAT . IndexedFasta as IndexedFasta <EOL> import CGAT . IndexedGenome as IndexedGenome <EOL> import CGAT . Genomics as Genomics <EOL> def readJunctions ( filename_junctions ) : <EOL> '''<STR_LIT>''' <EOL> junctions = collections . defaultdict ( dict ) <EOL> infile = open ( filename_junctions , "<STR_LIT:r>" ) <EOL> njunctions = <NUM_LIT:0> <EOL> for line in infile : <EOL> if line . startswith ( "<STR_LIT:#>" ) : <EOL> continue <EOL> data = line [ : - <NUM_LIT:1> ] . split ( "<STR_LIT:\t>" ) <EOL> if data [ <NUM_LIT:0> ] == "<STR_LIT>" : <EOL> continue <EOL> contig , strand , end , start , frame = data [ : <NUM_LIT:5> ] <EOL> start , end , frame = int ( start ) , int ( end ) , int ( frame ) <EOL> j = junctions [ contig ] <EOL> if end in j : <EOL> j [ end ] . add ( ( strand , start , frame ) ) <EOL> else : <EOL> j [ end ] = set ( ( ( strand , start , frame ) , ) ) <EOL> njunctions += <NUM_LIT:1> <EOL> infile . close ( ) <EOL> E . info ( "<STR_LIT>" % ( njunctions , len ( junctions ) ) ) <EOL> return junctions <EOL> class BaseAnnotator ( object ) : <EOL> '''<STR_LIT>''' <EOL> mHeader = ( ) <EOL> def __init__ ( self , fasta = None , * args , ** kwargs ) : <EOL> self . mFasta = fasta <EOL> def __str__ ( self ) : <EOL> return "<STR_LIT>" <EOL> def getHeader ( self ) : <EOL> '''<STR_LIT>''' <EOL> return "<STR_LIT:\t>" . join ( self . mHeader ) <EOL> class BaseAnnotatorSNP ( BaseAnnotator ) : <EOL> '''<STR_LIT>''' <EOL> mHeader = ( <EOL> "<STR_LIT>" , <EOL> "<STR_LIT>" , <EOL> "<STR_LIT>" , <EOL> "<STR_LIT>" , <EOL> "<STR_LIT>" , <EOL> "<STR_LIT>" , <EOL> "<STR_LIT>" , <EOL> "<STR_LIT>" ) <EOL> mAdditionalHeader = ( <EOL> "<STR_LIT>" , <EOL> "<STR_LIT>" ) <EOL> def __init__ ( self , * args , ** kwargs ) : <EOL> BaseAnnotator . __init__ ( self , * args , ** kwargs ) <EOL> def __str__ ( self ) : <EOL> return "<STR_LIT:\t>" . join ( map ( str , self . mSNP ) [ : len ( self . mHeader ) ] ) <EOL> def update ( self , snp ) : <EOL> '''<STR_LIT>''' <EOL> self . mSNP = snp <EOL> class BaseAnnotatorExons ( BaseAnnotator ) : <EOL> '''<STR_LIT>''' <EOL> mHeader = [ "<STR_LIT>" % x for x in ( "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ) ] <EOL> def __init__ ( self , filename_exons , * args , ** kwargs ) : <EOL> BaseAnnotator . __init__ ( self , * args , ** kwargs ) <EOL> exons = IndexedGenome . IndexedGenome ( ) <EOL> nexons = <NUM_LIT:0> <EOL> for g in GTF . iterator ( open ( filename_exons , "<STR_LIT:r>" ) ) : <EOL> exons . add ( g . contig , g . start , g . end , g ) <EOL> nexons += <NUM_LIT:1> <EOL> self . mExons = exons <EOL> E . info ( "<STR_LIT>" % ( nexons , len ( exons ) ) ) <EOL> def update ( self , snp ) : <EOL> '''<STR_LIT>''' <EOL> exons = list ( self . mExons . get ( snp . chromosome , <EOL> snp . pos , <EOL> snp . pos + <NUM_LIT:1> ) ) <EOL> if exons : <EOL> for start , end , gff in exons : <EOL> self . mOTranscripts = "<STR_LIT>" % gff [ "<STR_LIT>" ] <EOL> self . mOUsed = "<STR_LIT>" % gff [ "<STR_LIT>" ] <EOL> self . mOPos = "<STR_LIT:%s>" % gff [ "<STR_LIT>" ] <EOL> else : <EOL> self . mOTranscripts = "<STR_LIT>" <EOL> self . mOUsed = "<STR_LIT>" <EOL> self . mOPos = "<STR_LIT>" <EOL> def __str__ ( self ) : <EOL> return "<STR_LIT:\t>" . join ( ( self . mOTranscripts , self . mOUsed , self . mOPos ) ) <EOL> class BaseAnnotatorSpliceSites ( BaseAnnotator ) : <EOL> '''<STR_LIT>''' <EOL> mHeader = [ "<STR_LIT>" % x <EOL> for x in ( "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT:code>" ) ] <EOL> mIntronTypes = ( ( "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ) , <EOL> ( "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ) , <EOL> ( "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ) ) <EOL> mSize = <NUM_LIT:10> <EOL> def __init__ ( self , filename_junctions , * args , ** kwargs ) : <EOL> BaseAnnotator . __init__ ( self , * args , ** kwargs ) <EOL> junctions = IndexedGenome . IndexedGenome ( ) <EOL> infile = open ( filename_junctions , "<STR_LIT:r>" ) <EOL> njunctions = <NUM_LIT:0> <EOL> for line in infile : <EOL> if line . startswith ( "<STR_LIT:#>" ) : <EOL> continue <EOL> data = line [ : - <NUM_LIT:1> ] . split ( "<STR_LIT:\t>" ) <EOL> if data [ <NUM_LIT:0> ] == "<STR_LIT>" : <EOL> continue <EOL> contig , strand , end , start , frame , gene_id , transcript_id = data <EOL> start , end , frame = int ( start ) , int ( end ) , int ( frame ) <EOL> intron_start , intron_end = end + <NUM_LIT:1> , start <EOL> if strand == "<STR_LIT:->" : <EOL> lcontig = self . mFasta . getLength ( contig ) <EOL> intron_start , intron_end = lcontig - intron_end , lcontig - intron_start <EOL> junctions . add ( contig , intron_start , intron_start + self . mSize , <EOL> ( strand , intron_start , intron_end , gene_id , transcript_id ) ) <EOL> junctions . add ( contig , intron_end - self . mSize , intron_end , <EOL> ( strand , intron_start , intron_end , gene_id , transcript_id ) ) <EOL> njunctions += <NUM_LIT:1> <EOL> infile . close ( ) <EOL> self . mJunctions = junctions <EOL> E . info ( "<STR_LIT>" % <EOL> ( njunctions , len ( junctions ) ) ) <EOL> def hasSpliceMotif ( self , seq5 , seq3 ) : <EOL> """<STR_LIT>""" <EOL> for name , prime5 , prime3 in self . mIntronTypes : <EOL> if seq5 . startswith ( prime5 ) and seq3 . endswith ( prime3 ) : <EOL> return name , prime5 , prime3 <EOL> return None , None , None <EOL> def buildSequenceVariants ( self , seq , strand , pos , snp ) : <EOL> '''<STR_LIT>''' <EOL> is_negative_strand = Genomics . IsNegativeStrand ( strand ) <EOL> reference_base = snp . reference_base <EOL> if reference_base != "<STR_LIT:*>" and is_negative_strand : <EOL> reference_base = Genomics . complement ( reference_base ) <EOL> new_sequences = [ ] <EOL> is_homozygous = True <EOL> if reference_base != "<STR_LIT:*>" : <EOL> if seq [ pos ] . upper ( ) != reference_base . upper ( ) : <EOL> raise ValueError ( "<STR_LIT>" % <EOL> ( snp . pos , reference_base , seq [ pos ] , seq , pos , <EOL> "<STR_LIT:;>" . join ( map ( str , snp ) ) ) ) <EOL> variant_bases = Genomics . resolveAmbiguousNA ( snp . genotype ) <EOL> if len ( variant_bases ) == <NUM_LIT:1> : <EOL> is_homozygous = True <EOL> else : <EOL> is_homozygous = False <EOL> for variant_base in variant_bases : <EOL> if is_negative_strand : <EOL> variant_base = Genomics . complement ( variant_base ) <EOL> s = list ( seq ) <EOL> s [ pos ] = variant_base <EOL> s = "<STR_LIT>" . join ( s ) <EOL> if s != seq : <EOL> new_sequences . append ( s ) <EOL> else : <EOL> variants = snp . genotype . split ( "<STR_LIT:/>" ) <EOL> is_homozygous = False <EOL> for variant in variants : <EOL> s = list ( seq ) <EOL> if variant [ <NUM_LIT:0> ] == "<STR_LIT:+>" : <EOL> toinsert = variant [ <NUM_LIT:1> : ] . upper ( ) <EOL> if is_negative_strand : <EOL> toinsert = Genomics . complement ( toinsert ) <EOL> s . insert ( pos , toinsert ) <EOL> else : <EOL> s . insert ( pos + <NUM_LIT:1> , toinsert ) <EOL> elif variant [ <NUM_LIT:0> ] == "<STR_LIT:->" : <EOL> todelete = variant [ <NUM_LIT:1> : ] . upper ( ) <EOL> l = len ( todelete ) <EOL> if is_negative_strand : <EOL> xstart = max ( <NUM_LIT:0> , pos - l ) <EOL> xend = pos <EOL> todelete = todelete [ : min ( l , pos ) ] <EOL> else : <EOL> xstart = pos + <NUM_LIT:1> <EOL> xend = min ( self . mSize , pos + <NUM_LIT:1> + l ) <EOL> todelete = todelete [ : self . mSize - ( pos + <NUM_LIT:1> ) ] <EOL> deleted = "<STR_LIT>" . join ( s [ xstart : xend ] ) <EOL> if is_negative_strand : <EOL> deleted = Genomics . complement ( deleted ) <EOL> if deleted != todelete : <EOL> raise ValueError ( "<STR_LIT>" % <EOL> ( snp . pos , todelete , deleted , seq , pos , xstart , xend , <EOL> is_negative_strand , <EOL> "<STR_LIT:;>" . join ( map ( str , snp ) ) ) ) <EOL> del s [ xstart : xend ] <EOL> elif variant [ <NUM_LIT:0> ] == "<STR_LIT:*>" : <EOL> is_homozygous = True <EOL> else : <EOL> raise ValueError ( "<STR_LIT>" % variant [ <NUM_LIT:0> ] ) <EOL> s = "<STR_LIT>" . join ( s ) <EOL> if s != seq : <EOL> new_sequences . append ( s ) <EOL> return is_homozygous , new_sequences <EOL> def update ( self , snp ) : <EOL> '''<STR_LIT>''' <EOL> junctions = list ( self . mJunctions . get ( snp . chromosome , <EOL> snp . pos , <EOL> snp . pos + <NUM_LIT:1> ) ) <EOL> self . mResults = [ ] <EOL> reference_base = snp . reference_base <EOL> if junctions : <EOL> for start , end , data in junctions : <EOL> strand , intron_start , intron_end , gene_id , transcript_id = data <EOL> snp_pos = snp . pos <EOL> if strand == "<STR_LIT:->" : <EOL> lcontig = self . mFasta . getLength ( snp . chromosome ) <EOL> intron_start , intron_end = lcontig - intron_end , lcontig - intron_start <EOL> snp_pos = lcontig - snp_pos - <NUM_LIT:1> <EOL> seq5 = self . mFasta . getSequence ( <EOL> snp . chromosome , strand , intron_start , intron_start + self . mSize ) . upper ( ) <EOL> seq3 = self . mFasta . getSequence ( <EOL> snp . chromosome , strand , intron_end - self . mSize , intron_end ) . upper ( ) <EOL> splice_name , prime5 , prime3 = self . hasSpliceMotif ( seq5 , seq3 ) <EOL> if splice_name : <EOL> if intron_start <= snp_pos < intron_start + len ( prime5 ) : <EOL> has_overlap = "<STR_LIT:5>" <EOL> elif intron_end - len ( prime3 ) <= snp_pos < intron_end : <EOL> has_overlap = "<STR_LIT:3>" <EOL> else : <EOL> has_overlap = "<STR_LIT>" <EOL> else : <EOL> has_overlap = "<STR_LIT>" <EOL> E . debug ( "<STR_LIT>" % <EOL> ( snp . chromosome , strand , snp . pos , <EOL> intron_start , intron_end , strand , snp_pos , <EOL> str ( splice_name ) , <EOL> str ( has_overlap ) ) ) <EOL> if intron_start <= snp_pos < intron_start + self . mSize : <EOL> is_homozygous , new_seq5 = self . buildSequenceVariants ( <EOL> seq5 , strand , snp_pos - intron_start , snp ) <EOL> new_seqs = [ ( x , seq3 ) for x in new_seq5 ] <EOL> distance = snp_pos - intron_start <EOL> elif intron_end - self . mSize <= snp_pos < intron_end : <EOL> is_homozygous , new_seq3 = self . buildSequenceVariants ( <EOL> seq3 , strand , self . mSize - ( intron_end - snp_pos ) , snp ) <EOL> new_seqs = [ ( seq5 , x ) for x in new_seq3 ] <EOL> distance = intron_end - snp_pos - <NUM_LIT:1> <EOL> else : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> E . debug ( "<STR_LIT>" % <EOL> ( snp . chromosome , strand , snp . pos , <EOL> intron_start , intron_end , strand , <EOL> str ( splice_name ) , <EOL> distance , <EOL> str ( has_overlap ) , <EOL> str ( is_homozygous ) , <EOL> str ( new_seqs ) ) ) <EOL> if new_seqs == [ ] : <EOL> continue <EOL> for new_seq5 , new_seq3 in new_seqs : <EOL> new_splice_name , new_prime5 , new_prime3 = self . hasSpliceMotif ( <EOL> new_seq5 , new_seq3 ) <EOL> if splice_name is None and new_splice_name is None : <EOL> code = "<STR_LIT>" <EOL> elif new_splice_name is None : <EOL> code = "<STR_LIT:D>" <EOL> elif splice_name is None : <EOL> code = "<STR_LIT:C>" <EOL> elif splice_name == new_splice_name : <EOL> code = "<STR_LIT:S>" <EOL> elif splice_name != new_splice_name : <EOL> code = "<STR_LIT:N>" <EOL> if is_homozygous : <EOL> is_homozygous = "<STR_LIT:1>" <EOL> else : <EOL> is_homozygous = "<STR_LIT:0>" <EOL> if not splice_name : <EOL> splice_name = "<STR_LIT>" <EOL> if not new_splice_name : <EOL> new_splice_name = "<STR_LIT>" <EOL> self . mResults . append ( ( transcript_id , gene_id , has_overlap , is_homozygous , <EOL> splice_name , new_splice_name , str ( distance ) , code ) ) <EOL> def __str__ ( self ) : <EOL> if self . mResults : <EOL> return "<STR_LIT:\t>" . join ( [ "<STR_LIT:;>" . join ( x ) for x in zip ( * self . mResults ) ] ) <EOL> else : <EOL> return "<STR_LIT:\t>" * <NUM_LIT:8> <EOL> class BaseAnnotatorCodon ( BaseAnnotator ) : <EOL> '''<STR_LIT>''' <EOL> mHeader = [ "<STR_LIT:code>" , "<STR_LIT>" , "<STR_LIT>" , <EOL> "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ] <EOL> def __init__ ( self , annotations_file , junctions , * args , ** kwargs ) : <EOL> BaseAnnotator . __init__ ( self , * args , ** kwargs ) <EOL> self . mAnnotations = IndexedFasta . IndexedFasta ( annotations_file ) <EOL> self . mJunctions = junctions <EOL> def updateSNPs ( self , snp , is_negative_strand , pos ) : <EOL> '''<STR_LIT>''' <EOL> contig = snp . chromosome <EOL> lcontig = self . mFasta . getLength ( contig ) <EOL> reference_base = snp . reference_base <EOL> if snp . genotype in '<STR_LIT>' : <EOL> self . mVariantType . append ( "<STR_LIT:O>" ) <EOL> else : <EOL> self . mVariantType . append ( "<STR_LIT:E>" ) <EOL> if reference_base != "<STR_LIT:*>" and is_negative_strand : <EOL> reference_base = Genomics . complement ( reference_base ) <EOL> for reference_codon in self . mReferenceCodons : <EOL> self . mReferenceAAs . append ( Genomics . translate ( reference_codon ) ) <EOL> variant_bases = Genomics . resolveAmbiguousNA ( snp . genotype ) <EOL> if reference_codon [ pos ] != reference_base : <EOL> raise ValueError ( "<STR_LIT>" % <EOL> ( snp . pos , reference_codon , pos , reference_codon [ pos ] , reference_base , "<STR_LIT:;>" . join ( map ( str , snp ) ) ) ) <EOL> for variant_base in variant_bases : <EOL> if is_negative_strand : <EOL> variant_base = Genomics . complement ( variant_base ) <EOL> self . mVariantAAs . extend ( [ Genomics . translate ( x ) <EOL> for x in self . mVariantCodons ] ) <EOL> def updateIndels ( self , snp , is_negative_strand ) : <EOL> contig = snp . chromosome <EOL> lcontig = self . mFasta . getLength ( contig ) <EOL> code = self . mAnnotations . getSequence ( contig , "<STR_LIT:+>" , snp . pos , snp . pos + <NUM_LIT:2> ) <EOL> self . mCode = code <EOL> variants = snp . genotype . split ( "<STR_LIT:/>" ) <EOL> for variant in variants : <EOL> if variant [ <NUM_LIT:0> ] == "<STR_LIT:*>" : <EOL> self . mVariantType . append ( "<STR_LIT>" ) <EOL> elif variant [ <NUM_LIT:0> ] == "<STR_LIT:+>" : <EOL> toinsert = variant [ <NUM_LIT:1> : ] <EOL> self . mVariantType . append ( "<STR_LIT:I>" ) <EOL> elif variant [ <NUM_LIT:0> ] == "<STR_LIT:->" : <EOL> todelete = variant [ <NUM_LIT:1> : ] <EOL> self . mVariantType . append ( "<STR_LIT:D>" ) <EOL> else : <EOL> raise ValueError ( "<STR_LIT>" % variant [ <NUM_LIT:0> ] ) <EOL> if code [ <NUM_LIT:0> ] and code [ <NUM_LIT:1> ] not in '<STR_LIT>' : <EOL> return <EOL> if is_negative_strand : <EOL> variants = [ Genomics . complement ( x ) for x in variants ] <EOL> for reference_codon in self . mReferenceCodons : <EOL> variants = snp . genotype . split ( "<STR_LIT:/>" ) <EOL> variants = [ x [ <NUM_LIT:1> : ] for x in variants ] <EOL> for variant in variants : <EOL> if len ( variant ) % <NUM_LIT:3> != <NUM_LIT:0> : <EOL> self . mVariantCodons . append ( "<STR_LIT:!>" ) <EOL> else : <EOL> self . mVariantCodons . append ( variant ) <EOL> self . mVariantAAs . extend ( <EOL> [ Genomics . translate ( x ) for x in self . mVariantCodons ] ) <EOL> def update ( self , snp ) : <EOL> '''<STR_LIT>''' <EOL> contig = snp . chromosome <EOL> lcontig = self . mFasta . getLength ( contig ) <EOL> self . mReferenceCodons = [ ] <EOL> self . mVariantCodons = [ ] <EOL> self . mReferenceAAs = [ ] <EOL> self . mVariantAAs = [ ] <EOL> self . mVariantType = [ ] <EOL> if snp . pos >= lcontig : <EOL> E . warn ( "<STR_LIT>" % <EOL> ( contig , lcontig , "<STR_LIT:;>" . join ( map ( str , snp ) ) ) ) <EOL> self . mCode = "<STR_LIT:X>" <EOL> return <EOL> reference_base = snp . reference_base <EOL> code = self . mAnnotations . getSequence ( contig , "<STR_LIT:+>" , snp . pos , snp . pos + <NUM_LIT:1> ) <EOL> self . mCode = code <EOL> is_negative_strand = False <EOL> pos = None <EOL> if code in '<STR_LIT>' : <EOL> pos = "<STR_LIT>" . find ( code ) <EOL> is_negative_strand = code in '<STR_LIT>' <EOL> if is_negative_strand : <EOL> pos -= <NUM_LIT:3> <EOL> position = lcontig - snp . pos - <NUM_LIT:1> <EOL> strand = "<STR_LIT:->" <EOL> else : <EOL> position = snp . pos <EOL> strand = "<STR_LIT:+>" <EOL> offset = pos % <NUM_LIT:3> <EOL> start , end = position - offset , position - offset + <NUM_LIT:3> <EOL> self . mReferenceCodons = set ( ) <EOL> if self . mJunctions : <EOL> if contig in self . mJunctions : <EOL> j = self . mJunctions [ contig ] <EOL> for x in xrange ( start , end ) : <EOL> if x in j : <EOL> for junction_strand , next_start , frame in j [ x ] : <EOL> if ( junction_strand == "<STR_LIT:->" ) != is_negative_strand : <EOL> continue <EOL> if frame == <NUM_LIT:0> : <EOL> self . mReferenceCodons . add ( <EOL> self . mFasta . getSequence ( contig , strand , start , end ) . upper ( ) ) <EOL> continue <EOL> codon = self . mFasta . getSequence ( contig , strand , start , x + <NUM_LIT:1> ) + self . mFasta . getSequence ( <EOL> contig , strand , next_start , next_start + end - ( x + <NUM_LIT:1> ) ) <EOL> y = x + <NUM_LIT:1> - start <EOL> E . debug ( "<STR_LIT>" % <EOL> ( contig , strand , snp . pos , position , start , x + <NUM_LIT:1> , next_start , next_start + end - ( x + <NUM_LIT:1> ) , <EOL> codon [ : y ] , codon [ y : ] ) ) <EOL> assert len ( codon ) == <NUM_LIT:3> <EOL> self . mReferenceCodons . add ( codon . upper ( ) ) <EOL> if not self . mReferenceCodons : <EOL> self . mReferenceCodons . add ( <EOL> self . mFasta . getSequence ( contig , strand , start , end ) . upper ( ) ) <EOL> self . mReferenceCodons = tuple ( self . mReferenceCodons ) <EOL> if reference_base == "<STR_LIT:*>" : <EOL> self . updateIndels ( snp , is_negative_strand ) <EOL> else : <EOL> self . updateSNPs ( snp , is_negative_strand , pos ) <EOL> def __str__ ( self ) : <EOL> return "<STR_LIT:\t>" . join ( ( self . mCode , <EOL> "<STR_LIT:U+002C>" . join ( self . mReferenceCodons ) , <EOL> "<STR_LIT:U+002C>" . join ( self . mReferenceAAs ) , <EOL> "<STR_LIT:U+002C>" . join ( self . mVariantType ) , <EOL> "<STR_LIT:U+002C>" . join ( self . mVariantCodons ) , <EOL> "<STR_LIT:U+002C>" . join ( self . mVariantAAs ) , ) ) <EOL> def main ( argv = None ) : <EOL> """<STR_LIT>""" <EOL> if not argv : <EOL> argv = sys . argv <EOL> parser = E . OptionParser ( <EOL> version = "<STR_LIT>" , usage = globals ( ) [ "<STR_LIT>" ] ) <EOL> parser . add_option ( "<STR_LIT>" , "<STR_LIT>" , dest = "<STR_LIT>" , type = "<STR_LIT:string>" , <EOL> help = "<STR_LIT>" ) <EOL> parser . add_option ( "<STR_LIT>" , "<STR_LIT>" , dest = "<STR_LIT>" , type = "<STR_LIT:string>" , <EOL> help = "<STR_LIT>" ) <EOL> parser . add_option ( "<STR_LIT>" , "<STR_LIT>" , dest = "<STR_LIT>" , type = "<STR_LIT:string>" , <EOL> help = "<STR_LIT>" ) <EOL> parser . add_option ( "<STR_LIT>" , "<STR_LIT>" , dest = "<STR_LIT>" , type = "<STR_LIT:string>" , <EOL> help = "<STR_LIT>" ) <EOL> parser . add_option ( "<STR_LIT:-c>" , "<STR_LIT>" , dest = "<STR_LIT>" , type = "<STR_LIT:string>" , <EOL> help = "<STR_LIT>" ) <EOL> parser . add_option ( "<STR_LIT>" , "<STR_LIT>" , dest = "<STR_LIT>" , type = "<STR_LIT>" , <EOL> choices = ( "<STR_LIT>" , "<STR_LIT>" ) , <EOL> help = "<STR_LIT>" ) <EOL> parser . add_option ( "<STR_LIT>" , dest = "<STR_LIT>" , type = "<STR_LIT:string>" , <EOL> help = "<STR_LIT>" ) <EOL> parser . set_defaults ( <EOL> genome_file = None , <EOL> filename_annotations = None , <EOL> filename_exons = None , <EOL> filename_junctions = None , <EOL> input_format = "<STR_LIT>" , <EOL> vcf_sample = None , <EOL> filename_vcf = None , <EOL> ) <EOL> ( options , args ) = E . Start ( parser , argv = argv ) <EOL> ninput , nskipped , noutput = <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> <EOL> if options . genome_file : <EOL> fasta = IndexedFasta . IndexedFasta ( options . genome_file ) <EOL> else : <EOL> fasta = None <EOL> if options . filename_junctions : <EOL> junctions = readJunctions ( options . filename_junctions ) <EOL> else : <EOL> junctions = None <EOL> if options . input_format == "<STR_LIT>" : <EOL> iterator = pysam . Pileup . iterate ( sys . stdin ) <EOL> elif options . input_format == "<STR_LIT>" : <EOL> if not options . vcf_sample : <EOL> raise ValueError ( <EOL> "<STR_LIT>" ) <EOL> if not options . filename_vcf : <EOL> raise ValueError ( <EOL> "<STR_LIT>" ) <EOL> iterator = pysam . Pileup . iterate_from_vcf ( <EOL> options . filename_vcf , options . vcf_sample ) <EOL> modules = [ ] <EOL> modules . append ( BaseAnnotatorSNP ( ) ) <EOL> if options . filename_exons : <EOL> modules . append ( BaseAnnotatorExons ( options . filename_exons , fasta = fasta ) ) <EOL> if options . filename_annotations : <EOL> modules . append ( BaseAnnotatorCodon ( <EOL> options . filename_annotations , fasta = fasta , junctions = junctions ) ) <EOL> if options . filename_junctions : <EOL> modules . append ( <EOL> BaseAnnotatorSpliceSites ( options . filename_junctions , fasta = fasta ) ) <EOL> options . stdout . write ( "<STR_LIT:\t>" . join ( [ x . getHeader ( ) for x in modules ] ) + "<STR_LIT:\n>" ) <EOL> for snp in iterator : <EOL> ninput += <NUM_LIT:1> <EOL> if fasta : <EOL> try : <EOL> snp = snp . _replace ( chromosome = fasta . getToken ( snp . chromosome ) ) <EOL> except KeyError : <EOL> E . warn ( "<STR_LIT>" % <EOL> ( snp . chromosome , str ( snp ) ) ) <EOL> continue <EOL> for module in modules : <EOL> module . update ( snp ) <EOL> options . stdout . write ( "<STR_LIT:\t>" . join ( map ( str , modules ) ) + "<STR_LIT:\n>" ) <EOL> noutput += <NUM_LIT:1> <EOL> E . info ( "<STR_LIT>" % ( ninput , noutput , nskipped ) ) <EOL> E . Stop ( ) <EOL> if __name__ == "<STR_LIT:__main__>" : <EOL> sys . exit ( main ( sys . argv ) ) </s>
