From ff1b05d283026e23dd536af616ea02574c3c8ebd Mon Sep 17 00:00:00 2001 From: Richard Burhans Date: Thu, 22 Mar 2012 13:03:01 -0400 Subject: [PATCH] Updates to Human Genome Variation tools Rename Human Genome Variation to Phenotype Association --- tool_conf.xml.main | 30 +++-- tool_conf.xml.sample | 32 ++--- .../human_genome_variation/linkToGProfile.pl | 55 --------- .../BEAM2_wrapper.sh | 0 .../beam.xml | 0 .../ctd.pl | 0 .../ctd.xml | 0 .../disease_ontology_gene_fuzzy_selector.pl | 0 .../dividePgSnpAlleles.pl | 41 +++++++ .../dividePgSnpAlleles.xml | 52 ++++++++ .../freebayes.xml | 0 .../funDo.xml | 0 .../gpass.pl | 0 .../gpass.xml | 0 .../hilbertvis.sh | 0 .../hilbertvis.xml | 0 .../ldtools.xml | 0 .../ldtools_wrapper.sh | 0 .../linkToDavid.pl | 0 .../linkToDavid.xml | 0 tools/phenotype_association/linkToGProfile.pl | 89 ++++++++++++++ .../linkToGProfile.xml | 12 +- .../lped_to_geno.pl | 0 .../lps.xml | 0 .../lps_tool_wrapper.sh | 0 .../master2pg.pl | 0 .../master2pg.xml | 0 .../mergeSnps.pl | 0 .../pagetag.py | 0 .../pass.xml | 0 .../pass_wrapper.sh | 0 .../senatag.py | 0 .../sift.xml | 0 .../sift_variants_wrapper.sh | 0 .../snpFreq.xml | 0 .../snpFreq2.pl | 0 tools/phenotype_association/vcf2pgSnp.pl | 116 ++++++++++++++++++ tools/phenotype_association/vcf2pgSnp.xml | 78 ++++++++++++ tools/phenotype_association/vcf2pgSnpMult.pl | 81 ++++++++++++ 39 files changed, 499 insertions(+), 87 deletions(-) delete mode 100755 tools/human_genome_variation/linkToGProfile.pl rename tools/{human_genome_variation => phenotype_association}/BEAM2_wrapper.sh (100%) rename tools/{human_genome_variation => phenotype_association}/beam.xml (100%) rename tools/{human_genome_variation => phenotype_association}/ctd.pl (100%) rename tools/{human_genome_variation => phenotype_association}/ctd.xml (100%) rename tools/{human_genome_variation => phenotype_association}/disease_ontology_gene_fuzzy_selector.pl (100%) create mode 100755 tools/phenotype_association/dividePgSnpAlleles.pl create mode 100644 tools/phenotype_association/dividePgSnpAlleles.xml rename tools/{human_genome_variation => phenotype_association}/freebayes.xml (100%) rename tools/{human_genome_variation => phenotype_association}/funDo.xml (100%) rename tools/{human_genome_variation => phenotype_association}/gpass.pl (100%) rename tools/{human_genome_variation => phenotype_association}/gpass.xml (100%) rename tools/{human_genome_variation => phenotype_association}/hilbertvis.sh (100%) rename tools/{human_genome_variation => phenotype_association}/hilbertvis.xml (100%) rename tools/{human_genome_variation => phenotype_association}/ldtools.xml (100%) rename tools/{human_genome_variation => phenotype_association}/ldtools_wrapper.sh (100%) rename tools/{human_genome_variation => phenotype_association}/linkToDavid.pl (100%) rename tools/{human_genome_variation => phenotype_association}/linkToDavid.xml (100%) create mode 100755 tools/phenotype_association/linkToGProfile.pl rename tools/{human_genome_variation => phenotype_association}/linkToGProfile.xml (84%) rename tools/{human_genome_variation => phenotype_association}/lped_to_geno.pl (100%) rename tools/{human_genome_variation => phenotype_association}/lps.xml (100%) rename tools/{human_genome_variation => phenotype_association}/lps_tool_wrapper.sh (100%) rename tools/{human_genome_variation => phenotype_association}/master2pg.pl (100%) rename tools/{human_genome_variation => phenotype_association}/master2pg.xml (100%) rename tools/{human_genome_variation => phenotype_association}/mergeSnps.pl (100%) rename tools/{human_genome_variation => phenotype_association}/pagetag.py (100%) rename tools/{human_genome_variation => phenotype_association}/pass.xml (100%) rename tools/{human_genome_variation => phenotype_association}/pass_wrapper.sh (100%) rename tools/{human_genome_variation => phenotype_association}/senatag.py (100%) rename tools/{human_genome_variation => phenotype_association}/sift.xml (100%) rename tools/{human_genome_variation => phenotype_association}/sift_variants_wrapper.sh (100%) rename tools/{human_genome_variation => phenotype_association}/snpFreq.xml (100%) rename tools/{human_genome_variation => phenotype_association}/snpFreq2.pl (100%) create mode 100755 tools/phenotype_association/vcf2pgSnp.pl create mode 100644 tools/phenotype_association/vcf2pgSnp.xml create mode 100755 tools/phenotype_association/vcf2pgSnpMult.pl diff --git a/tool_conf.xml.main b/tool_conf.xml.main index 34325ae0d4b..1c08edc6c4a 100644 --- a/tool_conf.xml.main +++ b/tool_conf.xml.main @@ -192,21 +192,25 @@ -
+
- - - - - - - - - - - - + + + + + + + + + + + + + + + +
diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample index 4f05b66e670..11c57896b55 100644 --- a/tool_conf.xml.sample +++ b/tool_conf.xml.sample @@ -461,23 +461,25 @@
-
+
- - - - - - - - - - - - - - + + + + + + + + + + + + + + + +
diff --git a/tools/human_genome_variation/linkToGProfile.pl b/tools/human_genome_variation/linkToGProfile.pl deleted file mode 100755 index 00a12f85b92..00000000000 --- a/tools/human_genome_variation/linkToGProfile.pl +++ /dev/null @@ -1,55 +0,0 @@ -#!/usr/bin/env perl - -use strict; -use warnings; - -################################################### -# linkToGProfile.pl -# Generates a link to gprofile for a list of gene IDs. -# g:Profiler a web-based toolset for functional profiling of gene lists from large-scale experiments (2007) NAR 35 W193-W200 -################################################### - -if (!@ARGV or scalar @ARGV != 4) { - print "usage: linkToGProfile.pl infile.tab 1basedCol idType outfile\n"; - exit 1; -} - -my $in = shift @ARGV; -my $col = shift @ARGV; -my $type = shift @ARGV; -my $out = shift @ARGV; - -if ($col < 1) { - print "ERROR the column number should be 1 based counting\n"; - exit 1; -} -my @gene; -open(FH, $in) or die "Couldn't open $in, $!\n"; -while () { - chomp; - my @f = split(/\t/); - if (scalar @f < $col) { - print "ERROR there is no column $col in $in\n"; - exit 1; - } - if ($f[$col-1]) { push(@gene, $f[$col-1]); } -} -close FH or die "Couldn't close $in, $!\n"; - -my $link = 'http://biit.cs.ut.ee/gprofiler/index.cgi?organism=hsapiens&query=GENELIST&r_chr=1&r_start=start&r_end=end&analytical=1&domain_size_type=annotated&term=&significant=1&sort_by_structure=1&user_thr=1.00&output=png&prefix=TYPE'; -$link =~ s/TYPE/$type/; -my $g = join("+", @gene); -$link =~ s/GENELIST/$g/; -#print output -if (length $link > 2048) { - print "ERROR too many genes to fit in URL, please select a smaller set\n"; - exit; -} -open(FH, ">", $out) or die "Couldn't open $out, $!\n"; -print FH "g:Profiler link\n", - 'click here to send list of identifiers to g:Profiler', "\n", - '', "\n"; -close FH or die "Couldn't close $out, $!\n"; - -#also do link that prints text that could be pulled back into Galaxy? -exit; diff --git a/tools/human_genome_variation/BEAM2_wrapper.sh b/tools/phenotype_association/BEAM2_wrapper.sh similarity index 100% rename from tools/human_genome_variation/BEAM2_wrapper.sh rename to tools/phenotype_association/BEAM2_wrapper.sh diff --git a/tools/human_genome_variation/beam.xml b/tools/phenotype_association/beam.xml similarity index 100% rename from tools/human_genome_variation/beam.xml rename to tools/phenotype_association/beam.xml diff --git a/tools/human_genome_variation/ctd.pl b/tools/phenotype_association/ctd.pl similarity index 100% rename from tools/human_genome_variation/ctd.pl rename to tools/phenotype_association/ctd.pl diff --git a/tools/human_genome_variation/ctd.xml b/tools/phenotype_association/ctd.xml similarity index 100% rename from tools/human_genome_variation/ctd.xml rename to tools/phenotype_association/ctd.xml diff --git a/tools/human_genome_variation/disease_ontology_gene_fuzzy_selector.pl b/tools/phenotype_association/disease_ontology_gene_fuzzy_selector.pl similarity index 100% rename from tools/human_genome_variation/disease_ontology_gene_fuzzy_selector.pl rename to tools/phenotype_association/disease_ontology_gene_fuzzy_selector.pl diff --git a/tools/phenotype_association/dividePgSnpAlleles.pl b/tools/phenotype_association/dividePgSnpAlleles.pl new file mode 100755 index 00000000000..167847d10cc --- /dev/null +++ b/tools/phenotype_association/dividePgSnpAlleles.pl @@ -0,0 +1,41 @@ +#!/usr/bin/perl -w +use strict; + +#divide the alleles and their information into separate columns for pgSnp-like +#files. Keep any additional columns beyond the pgSnp ones. +#reads from stdin, writes to stdout +my $ref; +my $in; +if (@ARGV && $ARGV[0] =~ /-ref=(\d+)/) { + $ref = $1 -1; + if ($ref == -1) { undef $ref; } + shift @ARGV; +} +if (@ARGV) { + $in = shift @ARGV; +} + +open(FH, $in) or die "Couldn't open $in, $!\n"; +while () { + chomp; + my @f = split(/\t/); + my @a = split(/\//, $f[3]); + my @fr = split(/,/, $f[5]); + my @sc = split(/,/, $f[6]); + if ($f[4] == 1) { #homozygous add N, 0, 0 + if ($ref) { push(@a, $f[$ref]); } + else { push(@a, "N"); } + push(@fr, 0); + push(@sc, 0); + } + if ($f[4] > 2) { next; } #skip those with more than 2 alleles + print "$f[0]\t$f[1]\t$f[2]\t$a[0]\t$fr[0]\t$sc[0]\t$a[1]\t$fr[1]\t$sc[1]"; + if (scalar @f > 7) { + splice(@f, 0, 7); #remove first 7 + print "\t", join("\t", @f), "\n"; + }else { print "\n"; } +} +close FH; + +exit; + diff --git a/tools/phenotype_association/dividePgSnpAlleles.xml b/tools/phenotype_association/dividePgSnpAlleles.xml new file mode 100644 index 00000000000..dedbe6c0112 --- /dev/null +++ b/tools/phenotype_association/dividePgSnpAlleles.xml @@ -0,0 +1,52 @@ + diff --git a/tools/human_genome_variation/freebayes.xml b/tools/phenotype_association/freebayes.xml similarity index 100% rename from tools/human_genome_variation/freebayes.xml rename to tools/phenotype_association/freebayes.xml diff --git a/tools/human_genome_variation/funDo.xml b/tools/phenotype_association/funDo.xml similarity index 100% rename from tools/human_genome_variation/funDo.xml rename to tools/phenotype_association/funDo.xml diff --git a/tools/human_genome_variation/gpass.pl b/tools/phenotype_association/gpass.pl similarity index 100% rename from tools/human_genome_variation/gpass.pl rename to tools/phenotype_association/gpass.pl diff --git a/tools/human_genome_variation/gpass.xml b/tools/phenotype_association/gpass.xml similarity index 100% rename from tools/human_genome_variation/gpass.xml rename to tools/phenotype_association/gpass.xml diff --git a/tools/human_genome_variation/hilbertvis.sh b/tools/phenotype_association/hilbertvis.sh similarity index 100% rename from tools/human_genome_variation/hilbertvis.sh rename to tools/phenotype_association/hilbertvis.sh diff --git a/tools/human_genome_variation/hilbertvis.xml b/tools/phenotype_association/hilbertvis.xml similarity index 100% rename from tools/human_genome_variation/hilbertvis.xml rename to tools/phenotype_association/hilbertvis.xml diff --git a/tools/human_genome_variation/ldtools.xml b/tools/phenotype_association/ldtools.xml similarity index 100% rename from tools/human_genome_variation/ldtools.xml rename to tools/phenotype_association/ldtools.xml diff --git a/tools/human_genome_variation/ldtools_wrapper.sh b/tools/phenotype_association/ldtools_wrapper.sh similarity index 100% rename from tools/human_genome_variation/ldtools_wrapper.sh rename to tools/phenotype_association/ldtools_wrapper.sh diff --git a/tools/human_genome_variation/linkToDavid.pl b/tools/phenotype_association/linkToDavid.pl similarity index 100% rename from tools/human_genome_variation/linkToDavid.pl rename to tools/phenotype_association/linkToDavid.pl diff --git a/tools/human_genome_variation/linkToDavid.xml b/tools/phenotype_association/linkToDavid.xml similarity index 100% rename from tools/human_genome_variation/linkToDavid.xml rename to tools/phenotype_association/linkToDavid.xml diff --git a/tools/phenotype_association/linkToGProfile.pl b/tools/phenotype_association/linkToGProfile.pl new file mode 100755 index 00000000000..0f0ee242e09 --- /dev/null +++ b/tools/phenotype_association/linkToGProfile.pl @@ -0,0 +1,89 @@ +#!/usr/bin/env perl + +use strict; +use warnings; + +################################################### +# linkToGProfile.pl +# Generates a link to gprofile for a list of gene IDs. +# g:Profiler a web-based toolset for functional profiling of gene lists from large-scale experiments (2007) NAR 35 W193-W200 +################################################### + +if (!@ARGV or scalar @ARGV < 4) { + print "usage: linkToGProfile.pl infile.tab idType outfile -gene=1basedCol -chr=1basedCol -start=1basedCol -end=1basedCol\n"; + exit 1; +} + +my $in = shift @ARGV; +my $type = shift @ARGV; +my $out = shift @ARGV; + +my $col = 9999; #large unrealistic default +my $chr = 9999; +my $st = 9999; +my $end = 9999; +foreach (@ARGV) { + if (/gene=(\d+)/) { $col = $1; } + elsif (/chr=(\d+)/) { $chr = $1; } + elsif (/start=(\d+)/) { $st = $1; } + elsif (/end=(\d+)/) { $end = $1; } + elsif (/region=1/) { $type = 'region'; } +} + +if ($col < 1 or $chr < 1 or $st < 1 or $end < 1) { + print "ERROR the column number should be 1 based counting\n"; + exit 1; +} +my @gene; +my @pos; +open(FH, $in) or die "Couldn't open $in, $!\n"; +while () { + chomp; + my @f = split(/\t/); + if ($type ne 'region') { + if (scalar @f < $col) { + print "ERROR there is no column $col in $in for type $type\n"; + exit 1; + } + if ($f[$col-1]) { push(@gene, $f[$col-1]); } + }else { + if (scalar @f < $chr or scalar @f < $st or scalar @f < $end) { + print "ERROR there is not enough columns ($chr,$st,$end) in $in\n"; + exit 1; + } + if ($f[$chr-1]) { + $f[$chr-1] =~ s/chr//; + push(@pos, "$f[$chr-1]:$f[$st-1]:$f[$end-1]"); + } + } +} +close FH or die "Couldn't close $in, $!\n"; + +#region_query = 1 for coordinates X:1:10 +#can now do POST method +#http://biit.cs.ut.ee/gprofiler/index.cgi?organism=hsapiens&query=pax6&term=&analytical=1&user_thr=1&sort_by_structure=1&output=txt +my $g = join("+", @gene) if @gene; +if (@pos) { $g = join("+", @pos); } +my %params = ( +"analytical"=>1, +"organism"=>"hsapiens", +"query"=>$g, +"term"=>"", +"output"=>"png", +"prefix"=>$type, +"user_thr"=>"1.00" +); +if (@pos) { $params{"region_query"} = 1; } + +open(FH, ">", $out) or die "Couldn't open $out, $!\n"; +print FH "g:Profiler link\n"; +print FH '
'; +foreach my $k (keys %params) { + print FH "\n"; +} +print FH ''; +print FH '
', "\n"; +close FH or die "Couldn't close $out, $!\n"; + +#also do link that prints text that could be pulled back into Galaxy? +exit; diff --git a/tools/human_genome_variation/linkToGProfile.xml b/tools/phenotype_association/linkToGProfile.xml similarity index 84% rename from tools/human_genome_variation/linkToGProfile.xml rename to tools/phenotype_association/linkToGProfile.xml index 50e72e0d8e1..c9b10e266b5 100644 --- a/tools/human_genome_variation/linkToGProfile.xml +++ b/tools/phenotype_association/linkToGProfile.xml @@ -2,13 +2,17 @@ tools for functional profiling of gene lists - linkToGProfile.pl $input $numerical_column $type $out_file1 + linkToGProfile.pl $input $type $out_file1 -region=$region -gene=$genes -chr=${input.metadata.chromCol} -start=${input.metadata.startCol} -end=${input.metadata.endCol} - - + + + + + + @@ -31,7 +35,7 @@ - + diff --git a/tools/human_genome_variation/lped_to_geno.pl b/tools/phenotype_association/lped_to_geno.pl similarity index 100% rename from tools/human_genome_variation/lped_to_geno.pl rename to tools/phenotype_association/lped_to_geno.pl diff --git a/tools/human_genome_variation/lps.xml b/tools/phenotype_association/lps.xml similarity index 100% rename from tools/human_genome_variation/lps.xml rename to tools/phenotype_association/lps.xml diff --git a/tools/human_genome_variation/lps_tool_wrapper.sh b/tools/phenotype_association/lps_tool_wrapper.sh similarity index 100% rename from tools/human_genome_variation/lps_tool_wrapper.sh rename to tools/phenotype_association/lps_tool_wrapper.sh diff --git a/tools/human_genome_variation/master2pg.pl b/tools/phenotype_association/master2pg.pl similarity index 100% rename from tools/human_genome_variation/master2pg.pl rename to tools/phenotype_association/master2pg.pl diff --git a/tools/human_genome_variation/master2pg.xml b/tools/phenotype_association/master2pg.xml similarity index 100% rename from tools/human_genome_variation/master2pg.xml rename to tools/phenotype_association/master2pg.xml diff --git a/tools/human_genome_variation/mergeSnps.pl b/tools/phenotype_association/mergeSnps.pl similarity index 100% rename from tools/human_genome_variation/mergeSnps.pl rename to tools/phenotype_association/mergeSnps.pl diff --git a/tools/human_genome_variation/pagetag.py b/tools/phenotype_association/pagetag.py similarity index 100% rename from tools/human_genome_variation/pagetag.py rename to tools/phenotype_association/pagetag.py diff --git a/tools/human_genome_variation/pass.xml b/tools/phenotype_association/pass.xml similarity index 100% rename from tools/human_genome_variation/pass.xml rename to tools/phenotype_association/pass.xml diff --git a/tools/human_genome_variation/pass_wrapper.sh b/tools/phenotype_association/pass_wrapper.sh similarity index 100% rename from tools/human_genome_variation/pass_wrapper.sh rename to tools/phenotype_association/pass_wrapper.sh diff --git a/tools/human_genome_variation/senatag.py b/tools/phenotype_association/senatag.py similarity index 100% rename from tools/human_genome_variation/senatag.py rename to tools/phenotype_association/senatag.py diff --git a/tools/human_genome_variation/sift.xml b/tools/phenotype_association/sift.xml similarity index 100% rename from tools/human_genome_variation/sift.xml rename to tools/phenotype_association/sift.xml diff --git a/tools/human_genome_variation/sift_variants_wrapper.sh b/tools/phenotype_association/sift_variants_wrapper.sh similarity index 100% rename from tools/human_genome_variation/sift_variants_wrapper.sh rename to tools/phenotype_association/sift_variants_wrapper.sh diff --git a/tools/human_genome_variation/snpFreq.xml b/tools/phenotype_association/snpFreq.xml similarity index 100% rename from tools/human_genome_variation/snpFreq.xml rename to tools/phenotype_association/snpFreq.xml diff --git a/tools/human_genome_variation/snpFreq2.pl b/tools/phenotype_association/snpFreq2.pl similarity index 100% rename from tools/human_genome_variation/snpFreq2.pl rename to tools/phenotype_association/snpFreq2.pl diff --git a/tools/phenotype_association/vcf2pgSnp.pl b/tools/phenotype_association/vcf2pgSnp.pl new file mode 100755 index 00000000000..8c05e53c1d5 --- /dev/null +++ b/tools/phenotype_association/vcf2pgSnp.pl @@ -0,0 +1,116 @@ +#!/usr/bin/perl -w +use strict; + +#convert from a vcf file to a pgSnp file. +#frequency count = chromosome count +#either a single column/individual +#or all columns as a population + +my $in; +my $stCol = 9; +my $endCol; +if (@ARGV && scalar @ARGV == 2) { + $stCol = shift @ARGV; + $in = shift @ARGV; + if ($stCol eq 'all') { $stCol = 10; } + else { $endCol = $stCol; } + $stCol--; #go from 1 based to zero based column number + if ($stCol < 9) { + print "ERROR genotype fields don't start until column 10\n"; + exit; + } +}elsif (@ARGV && scalar @ARGV == 1) { + $in = shift @ARGV; +}elsif (@ARGV) { + print "usage: vcf2pgSnp.pl [indColNum default=all] file.vcf > file.pgSnp\n"; + exit; +} + +open(FH, $in) or die "Couldn't open $in, $!\n"; +while () { + chomp; + if (/^\s*#/) { next; } #skip comments/headers + if (/^\s*$/) { next; } #skip blank lines + my @f = split(/\t/); + #chr pos1base ID refNt altNt[,|D#|Int] quality filter info format geno1 ... + my $a; + my %nt; + my %all; + my $cnt = 0; + my $var; + if ($f[3] eq 'N') { next; } #ignore ref=N + if ($f[4] =~ /[DI]/ or $f[3] =~ /[DI]/) { next; } #don't do microsatellite + #if ($f[4] =~ /[ACTG],[ACTG]/) { next; } #only do positions with single alternate + if ($f[6] && !($f[6] eq '.' or $f[6] eq 'PASS')) { next; } #filtered for some reason + my $ind = 0; + if ($f[8] ne 'GT') { #more than just genotype + my @t = split(/:/, $f[8]); + foreach (@t) { if ($_ eq 'GT') { last; } $ind++; } + if ($ind == 0 && $f[8] !~ /^GT/) { die "ERROR couldn't find genotype in format $f[8]\n"; } + } + #count 0's, 1's, 2's + if (!$endCol) { $endCol = $#f; } + foreach my $col ($stCol .. $endCol) { + if ($ind > 0) { + my @t = split(/:/, $f[$col]); + $f[$col] = $t[$ind] . ":"; #only keep genotype part + } + if ($f[$col] =~ /^(0|1|2).(0|1|2)/) { + $nt{$1}++; + $nt{$2}++; + }elsif ($f[$col] =~ /^(0|1|2):/) { #chrY or male chrX, single + $nt{$1}++; + } #else ignore + } + if (%nt) { + if ($f[0] !~ /chr/) { $f[0] = "chr$f[0]"; } + print "$f[0]\t", ($f[1]-1), "\t$f[1]\t"; #position info + my $cnt = scalar(keys %nt); + my $fr; + my $sc; + my $all; + if (exists $nt{0}) { + $all = uc($f[3]); + $fr = $nt{0}; + $sc = 0; + } + if (!exists $nt{0} && exists $nt{1}) { + if ($f[4] =~ /([ACTG]),?/) { + $all = $1; + $fr = $nt{1}; + $sc = 0; + }else { die "bad variant nt $f[4] for nt 1"; } + }elsif (exists $nt{1}) { + if ($f[4] =~ /([ACTG]),?/) { + $all .= '/' . $1; + $fr .= ",$nt{1}"; + $sc .= ",0"; + }else { die "bad variant nt $f[4] for nt 1"; } + } + if (exists $nt{2}) { + if ($f[4] =~ /^[ACTG],([ACTG]),?/) { + $all .= '/' . $1; + $fr .= ",$nt{2}"; + $sc .= ",0"; + }else { die "bad variant nt $f[4] for nt 2"; } + } + if (exists $nt{3}) { + if ($f[4] =~ /^[ACTG],[ACTG],([ACTG])/) { + $all .= '/' . $1; + $fr .= ",$nt{3}"; + $sc .= ",0"; + }else { die "bad variant nt $f[4] for nt 3"; } + } + if (exists $nt{4}) { + if ($f[4] =~ /^[ACTG],[ACTG],[ACTG],([ACTG])/) { + $all .= '/' . $1; + $fr .= ",$nt{4}"; + $sc .= ",0"; + }else { die "bad variant nt $f[4] for nt 4"; } + } + print "$all\t$cnt\t$fr\t$sc\n"; + } +} +close FH; + +exit; diff --git a/tools/phenotype_association/vcf2pgSnp.xml b/tools/phenotype_association/vcf2pgSnp.xml new file mode 100644 index 00000000000..1caafbaffb9 --- /dev/null +++ b/tools/phenotype_association/vcf2pgSnp.xml @@ -0,0 +1,78 @@ + diff --git a/tools/phenotype_association/vcf2pgSnpMult.pl b/tools/phenotype_association/vcf2pgSnpMult.pl new file mode 100755 index 00000000000..09e0833da0a --- /dev/null +++ b/tools/phenotype_association/vcf2pgSnpMult.pl @@ -0,0 +1,81 @@ +#!/usr/bin/perl -w +use strict; + +#convert from a vcf file to a pgSnp file with multiple sets of the allele +# specific columns +#frequency count = chromosome count + +my $in; +my $stCol = 9; +my $endCol; +if (@ARGV && scalar @ARGV == 1) { + $in = shift @ARGV; +}else { + print "usage: vcf2pgSnpMult.pl file.vcf > file.pgSnpMult\n"; + exit; +} + +if ($in =~ /.gz$/) { + open(FH, "zcat $in |") or die "Couldn't open $in, $!\n"; +}else { + open(FH, $in) or die "Couldn't open $in, $!\n"; +} +while () { + chomp; + if (/^\s*#/) { next; } #skip comments/headers + if (/^\s*$/) { next; } #skip blank lines + my @f = split(/\t/); + #chr pos1base ID refNt altNt[,|D#|Int] quality filter info format geno1 ... + my $a; + my %nt; + my %all; + my $cnt = 0; + my $var; + if ($f[3] eq 'N') { next; } #ignore ref=N + if ($f[4] =~ /[DI]/ or $f[3] =~ /[DI]/) { next; } #don't do microsatellite + if ($f[6] && !($f[6] eq '.' or $f[6] eq 'PASS')) { next; } #filtered for some reason + my $ind = 0; + if ($f[8] ne 'GT') { #more than just genotype + my @t = split(/:/, $f[8]); + foreach (@t) { if ($_ eq 'GT') { last; } $ind++; } + if ($ind == 0 && $f[8] !~ /^GT/) { die "ERROR couldn't find genotype in format $f[8]\n"; } + } + if (!$endCol) { $endCol = $#f; } + #put f[3] => nt{0} and split f[4] for rest of nt{} + $nt{0} = $f[3]; + my @t = split(/,/, $f[4]); + for (my $i=0; $i<=$#t; $i++) { + my $j = $i + 1; + $nt{$j} = $t[$i]; + } + if ($f[0] !~ /chr/) { $f[0] = "chr$f[0]"; } + print "$f[0]\t", ($f[1]-1), "\t$f[1]"; #position info + foreach my $col ($stCol .. $endCol) { #add each individual (4 columns) + if ($ind > 0) { + my @t = split(/:/, $f[$col]); + $f[$col] = $t[$ind] . ":"; #only keep genotype part + } + print "\t"; + if ($f[$col] =~ /^(\d).(\d)/) { + my $a1 = $1; + my $a2 = $2; + if (!exists $nt{$a1}) { die "ERROR bad allele $a1 in $f[3] $f[4]\n"; } + if (!exists $nt{$a2}) { die "ERROR bad allele $a2 in $f[3] $f[4]\n"; } + if ($a1 eq $a2) { #homozygous + print "$nt{$a1}\t1\t2\t0"; + }else { #heterozygous + print "$nt{$a1}/$nt{$a2}\t2\t1,1\t0,0"; + } + }elsif ($f[$col] =~ /^(\d):/) { #chrY or male chrX, single + my $a1 = $1; + if (!exists $nt{$a1}) { die "ERROR bad allele $a1 in $f[3] $f[4]\n"; } + print "$nt{$a1}\t1\t1\t0"; + }else { #don't know how to parse + die "ERROR unknown genotype $f[$col]\n"; + } + } + print "\n"; #end this SNP +} +close FH; + +exit;