diff --git a/lib/galaxy_utils/sequence/vcf.py b/lib/galaxy_utils/sequence/vcf.py
new file mode 100644
index 00000000000..ec56ab151e3
--- /dev/null
+++ b/lib/galaxy_utils/sequence/vcf.py
@@ -0,0 +1,92 @@
+#Dan Blankenberg
+#See: http://1000genomes.org/wiki/doku.php?id=1000_genomes:analysis:vcf3.3
+#See: http://1000genomes.org/wiki/doku.php?id=1000_genomes:analysis:variant_call_format
+
+class VariantCall( object ):
+ version = None
+ header_startswith = None
+ required_header_fields = None
+ required_header_length = None
+
+ @classmethod
+ def get_class_by_format( cls, format ):
+ assert format in VCF_FORMATS, 'Unknown format type specified: %s' % format
+ return VCF_FORMATS[ format ]
+
+ def __init__( self, vcf_line, metadata, sample_names ):
+ raise Exception( 'Abstract Method' )
+
+class VariantCall33( VariantCall ):
+ version = 'VCFv3.3'
+ header_startswith = '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO'
+ required_header_fields = header_startswith.split( '\t' )
+ required_header_length = len( required_header_fields )
+
+ def __init__( self, vcf_line, metadata, sample_names ):
+ self.line = vcf_line.rstrip( '\n\r' )
+ self.metadata = metadata
+ self.sample_names = sample_names
+ self.format = None
+ self.sample_values = []
+
+ #parse line
+ self.fields = self.line.split( '\t' )
+ if sample_names:
+ assert len( self.fields ) == self.required_header_length + len( sample_names ) + 1, 'Provided VCF line (%s) has wrong length (expected: %i)' % ( self.line, self.required_header_length + len( sample_names ) + 1 )
+ else:
+ assert len( self.fields ) == self.required_header_length, 'Provided VCF line (%s) has wrong length (expected: %i)' % ( self.line, self.required_header_length)
+ self.chrom, self.pos, self.id, self.ref, self.alt, self.qual, self.filter, self.info = self.fields[ :self.required_header_length ]
+ self.pos = int( self.pos )
+ self.alt = self.alt.split( ',' )
+ self.qual = float( self.qual )
+ if len( self.fields ) > self.required_header_length:
+ self.format = self.fields[ self.required_header_length ].split( ':' )
+ for sample_value in self.fields[ self.required_header_length + 1: ]:
+ self.sample_values.append( sample_value.split( ':' ) )
+
+#VCF Format version lookup dict
+VCF_FORMATS = {}
+for format in [ VariantCall33 ]:
+ VCF_FORMATS[format.version] = format
+
+class Reader( object ):
+ def __init__( self, fh ):
+ self.vcf_file = fh
+ self.metadata = {}
+ self.header_fields = None
+ self.sample_names = []
+ self.vcf_class = None
+ while True:
+ line = self.vcf_file.readline()
+ assert line, 'Invalid VCF file provided.'
+ line = line.rstrip( '\r\n' )
+ if self.vcf_class and line.startswith( self.vcf_class.header_startswith ):
+ self.header_fields = line.split( '\t' )
+ if len( self.header_fields ) > self.vcf_class.required_header_length:
+ for sample_name in self.header_fields[ self.vcf_class.required_header_length + 1 : ]:
+ self.sample_names.append( sample_name )
+ break
+ assert line.startswith( '##' ), 'Non-metadata line found before header'
+ line = line[2:] #strip ##
+ metadata = line.split( '=', 1 )
+ metadata_name = metadata[ 0 ]
+ if len( metadata ) == 2:
+ metadata_value = metadata[ 1 ]
+ else:
+ metadata_value = None
+ if metadata_name in self.metadata:
+ if not isinstance( self.metadata[ metadata_name ], list ):
+ self.metadata[ metadata_name ] = [ self.metadata[ metadata_name ] ]
+ self.metadata[ metadata_name ].append( metadata_value )
+ else:
+ self.metadata[ metadata_name ] = metadata_value
+ if metadata_name == 'fileformat':
+ self.vcf_class = VariantCall.get_class_by_format( metadata_value )
+ def next( self ):
+ line = self.vcf_file.readline()
+ if not line:
+ raise StopIteration
+ return self.vcf_class( line, self.metadata, self.sample_names )
+ def __iter__( self ):
+ while True:
+ yield self.next()
diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample
index f6e9f86ad91..1ec68813a9f 100644
--- a/tool_conf.xml.sample
+++ b/tool_conf.xml.sample
@@ -144,6 +144,7 @@
+
diff --git a/tools/maf/vcf_to_maf_customtrack.py b/tools/maf/vcf_to_maf_customtrack.py
new file mode 100644
index 00000000000..da193dd46c6
--- /dev/null
+++ b/tools/maf/vcf_to_maf_customtrack.py
@@ -0,0 +1,143 @@
+#Dan Blankenberg
+from optparse import OptionParser
+import galaxy_utils.sequence.vcf
+
+from galaxy import eggs
+import pkg_resources; pkg_resources.require( "bx-python" )
+import bx.align.maf
+
+UNKNOWN_NUCLEOTIDE = '*'
+
+class PopulationVCFParser( object ):
+ def __init__( self, reader, name ):
+ self.reader = reader
+ self.name = name
+ self.counter = 0
+ def next( self ):
+ rval = []
+ vc = self.reader.next()
+ for i, allele in enumerate( vc.alt ):
+ rval.append( ( '%s_%i.%i' % ( self.name, i + 1, self.counter + 1 ), allele ) )
+ self.counter += 1
+ return ( vc, rval )
+ def __iter__( self ):
+ while True:
+ yield self.next()
+
+class SampleVCFParser( object ):
+ def __init__( self, reader ):
+ self.reader = reader
+ self.counter = 0
+ def next( self ):
+ rval = []
+ vc = self.reader.next()
+ alleles = [ vc.ref ] + vc.alt
+
+ if 'GT' in vc.format:
+ gt_index = vc.format.index( 'GT' )
+ for sample_name, sample_value in zip( vc.sample_names, vc.sample_values ):
+ gt_indexes = []
+ for i in sample_value[ gt_index ].replace( '|', '/' ).replace( '\\', '/' ).split( '/' ): #Do we need to consider phase here?
+ try:
+ gt_indexes.append( int( i ) )
+ except:
+ gt_indexes.append( None )
+ for i, allele_i in enumerate( gt_indexes ):
+ if allele_i is not None:
+ rval.append( ( '%s_%i.%i' % ( sample_name, i + 1, self.counter + 1 ), alleles[ allele_i ] ) )
+ self.counter += 1
+ return ( vc, rval )
+ def __iter__( self ):
+ while True:
+ yield self.next()
+
+def main():
+ usage = "usage: %prog [options] output_file dbkey inputfile pop_name"
+ parser = OptionParser( usage=usage )
+ parser.add_option( "-p", "--population", action="store_true", dest="population", default=False, help="Create MAF on a per population basis")
+ parser.add_option( "-s", "--sample", action="store_true", dest="sample", default=False, help="Create MAF on a per sample basis")
+ parser.add_option( "-n", "--name", dest="name", default='Unknown Custom Track', help="Name for Custom Track")
+
+
+ ( options, args ) = parser.parse_args()
+
+ if len ( args ) < 3:
+ parser.error( "Need to specify an output file, a dbkey and at least one input file" )
+
+ if not ( options.population ^ options.sample ):
+ parser.error( 'You must specify either a per population conversion or a per sample conversion, but not both' )
+
+ out = open( args.pop(0), 'wb' )
+ out.write( 'track name="%s" visibility=pack\n' % options.name.replace( "\"", "'" ) )
+
+ maf_writer = bx.align.maf.Writer( out )
+
+ dbkey = args.pop(0)
+
+ vcf_files = []
+ if options.population:
+ i = 0
+ while args:
+ filename = args.pop( 0 )
+ pop_name = args.pop( 0 ).replace( ' ', '_' )
+ if not pop_name:
+ pop_name = 'population_%i' % ( i + 1 )
+ vcf_files.append( PopulationVCFParser( galaxy_utils.sequence.vcf.Reader( open( filename ) ), pop_name ) )
+ i += 1
+ else:
+ while args:
+ filename = args.pop( 0 )
+ vcf_files.append( SampleVCFParser( galaxy_utils.sequence.vcf.Reader( open( filename ) ) ) )
+
+ non_spec_skipped = 0
+ for vcf_file in vcf_files:
+ for vc, variants in vcf_file:
+ num_ins = 0
+ num_dels = 0
+ for variant_name, variant_text in variants:
+ if 'D' in variant_text:
+ num_dels = max( num_dels, int( variant_text[1:] ) )
+ elif 'I' in variant_text:
+ num_ins = max( num_ins, len( variant_text ) - 1 )
+
+ alignment = bx.align.maf.Alignment()
+ ref_text = vc.ref + '-' * num_ins + UNKNOWN_NUCLEOTIDE * ( num_dels - len( vc.ref ) )
+ start_pos = vc.pos - 1
+ if num_dels and start_pos:
+ ref_text = UNKNOWN_NUCLEOTIDE + ref_text
+ start_pos -= 1
+ alignment.add_component( bx.align.maf.Component( src='%s.chr%s' % ( dbkey, vc.chrom ), start = start_pos, size = len( ref_text.replace( '-', '' ) ), strand = '+', src_size = start_pos + len( ref_text ), text = ref_text ) )
+ for variant_name, variant_text in variants:
+ #FIXME:
+ ## skip non-spec. compliant data, see: http://1000genomes.org/wiki/doku.php?id=1000_genomes:analysis:vcf3.3 for format spec
+ ## this check is due to data having indels not represented in the published format spec,
+ ## e.g. 1000 genomes pilot 1 indel data: ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/pilot_data/release/2010_03/pilot1/indels/CEU.SRP000031.2010_03.indels.sites.vcf.gz
+ if variant_text and variant_text[0] in [ '-', '+' ]:
+ non_spec_skipped += 1
+ continue
+
+ #do we need a left padding unknown nucleotide (do we have deletions)?
+ if num_dels and start_pos:
+ var_text = UNKNOWN_NUCLEOTIDE
+ else:
+ var_text = ''
+ if 'D' in variant_text:
+ cur_num_del = int( variant_text[1:] )
+ pre_del = min( len( vc.ref ), cur_num_del )
+ post_del = cur_num_del - pre_del
+ var_text = var_text + '-' * pre_del + '-' * num_ins + '-' * post_del
+ var_text = var_text + UNKNOWN_NUCLEOTIDE * ( len( ref_text ) - len( var_text ) )
+ elif 'I' in variant_text:
+ cur_num_ins = len( variant_text ) - 1
+ var_text = var_text + vc.ref + variant_text[1:] + '-' * ( num_ins - cur_num_ins ) + UNKNOWN_NUCLEOTIDE * max( 0, ( num_dels - 1 ) )
+ else:
+ var_text = var_text + variant_text + '-' * num_ins + UNKNOWN_NUCLEOTIDE * ( num_dels - len( vc.ref ) )
+ alignment.add_component( bx.align.maf.Component( src=variant_name, start = 0, size = len( var_text.replace( '-', '' ) ), strand = '+', src_size = len( var_text.replace( '-', '' ) ), text = var_text ) )
+ maf_writer.write( alignment )
+
+ maf_writer.close()
+
+ if non_spec_skipped:
+ print 'Skipped %i non-specification compliant indels.' % non_spec_skipped
+
+if __name__ == "__main__": main()
diff --git a/tools/maf/vcf_to_maf_customtrack.xml b/tools/maf/vcf_to_maf_customtrack.xml
new file mode 100644
index 00000000000..6ad5309c9ab
--- /dev/null
+++ b/tools/maf/vcf_to_maf_customtrack.xml
@@ -0,0 +1,121 @@
+
+ for display at UCSC
+ vcf_to_maf_customtrack.py $out_file1 ${vcf_source_type.vcf_file[0].vcf_input.dbkey} ${vcf_source_type.vcf_source} -n '$track_name'
+ ##
+ #for $vcf_repeat in $vcf_source_type.vcf_file
+ '${vcf_repeat.vcf_input}'
+ #if $vcf_source_type.vcf_source == '-p'
+ '${vcf_repeat.population_name}'
+ #end if
+ #end for
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+**What it does**
+
+This tool converts a Variant Call Format (VCF) file into a Multiple Alignment Format (MAF) custom track file suitable for display at genome browsers.
+
+This file should be used for display purposes only (e.g as a UCSC Custom Track). Performing an analysis using the output created by this tool as input is not recommended; the source VCF file should be used when performing an analysis.
+
+*Unknown nucleotides* are represented as '*' as required to allow the display to draw properly; these include e.g. reference bases which appear before a deletion and are not available without querying the original reference sequence.
+
+**Example**
+
+Starting with a VCF::
+
+ ##fileformat=VCFv3.3
+ ##fileDate=20090805
+ ##source=myImputationProgramV3.1
+ ##reference=1000GenomesPilot-NCBI36
+ ##phasing=partial
+ ##INFO=NS,1,Integer,"Number of Samples With Data"
+ ##INFO=DP,1,Integer,"Total Depth"
+ ##INFO=AF,-1,Float,"Allele Frequency"
+ ##INFO=AA,1,String,"Ancestral Allele"
+ ##INFO=DB,0,Flag,"dbSNP membership, build 129"
+ ##INFO=H2,0,Flag,"HapMap2 membership"
+ ##FILTER=q10,"Quality below 10"
+ ##FILTER=s50,"Less than 50% of samples have data"
+ ##FORMAT=GT,1,String,"Genotype"
+ ##FORMAT=GQ,1,Integer,"Genotype Quality"
+ ##FORMAT=DP,1,Integer,"Read Depth"
+ ##FORMAT=HQ,2,Integer,"Haplotype Quality"
+ #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA00003
+ 20 14370 rs6054257 G A 29 0 NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ 0|0:48:1:51,51 1|0:48:8:51,51 1/1:43:5:-1,-1
+ 20 17330 . T A 3 q10 NS=3;DP=11;AF=0.017 GT:GQ:DP:HQ 0|0:49:3:58,50 0|1:3:5:65,3 0/0:41:3:-1,-1
+ 20 1110696 rs6040355 A G,T 67 0 NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQ:DP:HQ 1|2:21:6:23,27 2|1:2:0:18,2 2/2:35:4:-1,-1
+ 20 1230237 . T . 47 0 NS=3;DP=13;AA=T GT:GQ:DP:HQ 0|0:54:7:56,60 0|0:48:4:51,51 0/0:61:2:-1,-1
+ 20 1234567 microsat1 G D4,IGA 50 0 NS=3;DP=9;AA=G GT:GQ:DP 0/1:35:4 0/2:17:2 1/1:40:3
+
+
+
+
+Under the following conditions: **VCF Source type:** *Per Population (file)*, **Name for this population:** *CHB+JPT*
+Results in the following MAF custom track::
+
+ track name="Galaxy Custom Track" visibility=pack
+ ##maf version=1
+ a score=0
+ s hg18.chr20 14369 1 + 14370 G
+ s CHB+JPT_1.1 0 1 + 1 A
+
+ a score=0
+ s hg18.chr20 17329 1 + 17330 T
+ s CHB+JPT_1.2 0 1 + 1 A
+
+ a score=0
+ s hg18.chr20 1110695 1 + 1110696 A
+ s CHB+JPT_1.3 0 1 + 1 G
+ s CHB+JPT_2.3 0 1 + 1 T
+
+ a score=0
+ s hg18.chr20 1230236 1 + 1230237 T
+ s CHB+JPT_1.4 0 1 + 1 .
+
+ a score=0
+ s hg18.chr20 1234565 5 + 1234572 *G--***
+ s CHB+JPT_1.5 0 1 + 1 *------
+ s CHB+JPT_2.5 0 7 + 7 *GGA***
+
+
+
+
+