mirror of
https://github.com/galaxyproject/galaxy.git
synced 2026-09-24 16:30:27 +08:00
Add a VCF to MAF Custom Track converter tool. 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.
This commit is contained in:
Daniel Blankenberg
@@ -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()
|
||||
@@ -144,6 +144,7 @@
|
||||
<tool file="visualization/GMAJ.xml" />
|
||||
<tool file="visualization/LAJ.xml" />
|
||||
<tool file="visualization/build_ucsc_custom_track.xml" />
|
||||
<tool file="maf/vcf_to_maf_customtrack.xml" />
|
||||
</section>
|
||||
<section name="Regional Variation" id="regVar">
|
||||
<tool file="regVariation/windowSplitter.xml" />
|
||||
|
||||
@@ -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()
|
||||
@@ -0,0 +1,121 @@
|
||||
<tool id="vcf_to_maf_customtrack1" name="VCF to MAF Custom Track">
|
||||
<description>for display at UCSC</description>
|
||||
<command interpreter="python">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
|
||||
</command>
|
||||
<inputs>
|
||||
<param name="track_name" type="text" label="Custom Track Name" value="Galaxy Custom Track" size="30" />
|
||||
<conditional name="vcf_source_type">
|
||||
<param name="vcf_source" type="select" label="VCF Source Source Type">
|
||||
<option value="-p" selected="true">Per Population (file)</option>
|
||||
<option value="-s">Per Sample</option>
|
||||
</param>
|
||||
<when value="-p">
|
||||
<repeat name="vcf_file" title="VCF population file">
|
||||
<param format="tabular" name="vcf_input" type="data" label="VCF file"/>
|
||||
<param name="population_name" type="text" label="Name for this population" value=""/>
|
||||
</repeat>
|
||||
</when>
|
||||
<when value="-s">
|
||||
<repeat name="vcf_file" title="VCF sample file">
|
||||
<param format="tabular" name="vcf_input" type="data" label="VCF file"/>
|
||||
<!-- add column count validator >= 8? -->
|
||||
</repeat>
|
||||
</when>
|
||||
</conditional>
|
||||
</inputs>
|
||||
<outputs>
|
||||
<data format="mafcustomtrack" name="out_file1" />
|
||||
</outputs>
|
||||
<!-- <tests>
|
||||
<test>
|
||||
<param name="track_name" value="Galaxy Custom Track"/>
|
||||
<param name="vcf_source" value="Per Population"/>
|
||||
<param name="vcf_input" value="vcf_to_maf_in.vcf" ftype="tabular"/>
|
||||
<param name="population_name" value=""/>
|
||||
<output name="out_file1" file="vcf_to_maf_population_out.mafcustomtrack"/>
|
||||
</test>
|
||||
<test>
|
||||
<param name="track_name" value="Galaxy Custom Track"/>
|
||||
<param name="vcf_source" value="Per Sample"/>
|
||||
<param name="vcf_input" value="vcf_to_maf_in.vcf" ftype="tabular"/>
|
||||
<output name="out_file1" file="vcf_to_maf_sample_out.mafcustomtrack"/>
|
||||
</test>
|
||||
</tests> -->
|
||||
<help>
|
||||
**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***
|
||||
|
||||
|
||||
</help>
|
||||
</tool>
|
||||
|
||||
Reference in New Issue
Block a user