mirror of
https://github.com/galaxyproject/galaxy.git
synced 2026-09-24 16:30:27 +08:00
Updates to Human Genome Variation tools
Rename Human Genome Variation to Phenotype Association
This commit is contained in:
Richard Burhans
+17
-13
@@ -192,21 +192,25 @@
|
||||
<tool file="taxonomy/lca.xml" />
|
||||
<tool file="taxonomy/poisson2test.xml" />
|
||||
</section>
|
||||
<section name="Human Genome Variation" id="hgv">
|
||||
<section name="Phenotype Association" id="hgv">
|
||||
<tool file="evolution/codingSnps.xml" />
|
||||
<tool file="evolution/add_scores.xml" />
|
||||
<tool file="human_genome_variation/sift.xml" />
|
||||
<tool file="human_genome_variation/linkToGProfile.xml" />
|
||||
<tool file="human_genome_variation/linkToDavid.xml"/>
|
||||
<tool file="human_genome_variation/ctd.xml" />
|
||||
<tool file="human_genome_variation/funDo.xml" />
|
||||
<tool file="human_genome_variation/snpFreq.xml" />
|
||||
<tool file="human_genome_variation/ldtools.xml" />
|
||||
<tool file="human_genome_variation/pass.xml" />
|
||||
<tool file="human_genome_variation/gpass.xml" />
|
||||
<tool file="human_genome_variation/beam.xml" />
|
||||
<tool file="human_genome_variation/lps.xml" />
|
||||
<tool file="human_genome_variation/hilbertvis.xml" />
|
||||
<tool file="phenotype_association/sift.xml" />
|
||||
<tool file="phenotype_association/linkToGProfile.xml" />
|
||||
<tool file="phenotype_association/linkToDavid.xml"/>
|
||||
<tool file="phenotype_association/ctd.xml" />
|
||||
<tool file="phenotype_association/funDo.xml" />
|
||||
<tool file="phenotype_association/snpFreq.xml" />
|
||||
<tool file="phenotype_association/ldtools.xml" />
|
||||
<tool file="phenotype_association/pass.xml" />
|
||||
<tool file="phenotype_association/gpass.xml" />
|
||||
<tool file="phenotype_association/beam.xml" />
|
||||
<tool file="phenotype_association/lps.xml" />
|
||||
<tool file="phenotype_association/hilbertvis.xml" />
|
||||
<tool file="phenotype_association/freebayes.xml" />
|
||||
<tool file="phenotype_association/master2pg.xml" />
|
||||
<tool file="phenotype_association/vcf2pgSnp.xml" />
|
||||
<tool file="phenotype_association/dividePgSnpAlleles.xml" />
|
||||
</section>
|
||||
<section name="Genome Diversity" id="gd">
|
||||
<tool file="genome_diversity/extract_primers.xml" />
|
||||
|
||||
+17
-15
@@ -461,23 +461,25 @@
|
||||
<tool file="rgenetics/rgGLM.xml"/>
|
||||
<tool file="rgenetics/rgManQQ.xml"/>
|
||||
</section>
|
||||
<section name="Human Genome Variation" id="hgv">
|
||||
<section name="Phenotype Association" id="hgv">
|
||||
<tool file="evolution/codingSnps.xml" />
|
||||
<tool file="evolution/add_scores.xml" />
|
||||
<tool file="human_genome_variation/sift.xml" />
|
||||
<tool file="human_genome_variation/linkToGProfile.xml" />
|
||||
<tool file="human_genome_variation/linkToDavid.xml"/>
|
||||
<tool file="human_genome_variation/ctd.xml" />
|
||||
<tool file="human_genome_variation/funDo.xml" />
|
||||
<tool file="human_genome_variation/snpFreq.xml" />
|
||||
<tool file="human_genome_variation/ldtools.xml" />
|
||||
<tool file="human_genome_variation/pass.xml" />
|
||||
<tool file="human_genome_variation/gpass.xml" />
|
||||
<tool file="human_genome_variation/beam.xml" />
|
||||
<tool file="human_genome_variation/lps.xml" />
|
||||
<tool file="human_genome_variation/hilbertvis.xml" />
|
||||
<tool file="human_genome_variation/freebayes.xml" />
|
||||
<tool file="human_genome_variation/master2pg.xml" />
|
||||
<tool file="phenotype_association/sift.xml" />
|
||||
<tool file="phenotype_association/linkToGProfile.xml" />
|
||||
<tool file="phenotype_association/linkToDavid.xml"/>
|
||||
<tool file="phenotype_association/ctd.xml" />
|
||||
<tool file="phenotype_association/funDo.xml" />
|
||||
<tool file="phenotype_association/snpFreq.xml" />
|
||||
<tool file="phenotype_association/ldtools.xml" />
|
||||
<tool file="phenotype_association/pass.xml" />
|
||||
<tool file="phenotype_association/gpass.xml" />
|
||||
<tool file="phenotype_association/beam.xml" />
|
||||
<tool file="phenotype_association/lps.xml" />
|
||||
<tool file="phenotype_association/hilbertvis.xml" />
|
||||
<tool file="phenotype_association/freebayes.xml" />
|
||||
<tool file="phenotype_association/master2pg.xml" />
|
||||
<tool file="phenotype_association/vcf2pgSnp.xml" />
|
||||
<tool file="phenotype_association/dividePgSnpAlleles.xml" />
|
||||
</section>
|
||||
<section name="Genome Diversity" id="gd">
|
||||
<tool file="genome_diversity/extract_primers.xml" />
|
||||
|
||||
@@ -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 (<FH>) {
|
||||
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 "<html><head><title>g:Profiler link</title></head><body>\n",
|
||||
'<A TARGET=_BLANK HREF="', $link, '">click here to send list of identifiers to g:Profiler</A>', "\n",
|
||||
'</body></html>', "\n";
|
||||
close FH or die "Couldn't close $out, $!\n";
|
||||
|
||||
#also do link that prints text that could be pulled back into Galaxy?
|
||||
exit;
|
||||
+41
@@ -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 (<FH>) {
|
||||
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;
|
||||
|
||||
@@ -0,0 +1,52 @@
|
||||
<tool id="dividePgSnp" name="Separate alleles" hidden="false">
|
||||
<description>in a pgSnp file</description>
|
||||
<command interpreter="perl">
|
||||
dividePgSnpAlleles.pl -ref=$ref_column $input1 > $out_file1
|
||||
</command>
|
||||
<inputs>
|
||||
<param format="interval" name="input1" type="data" label="Personal genome SNP file" />
|
||||
<param name="ref_column" type="data_column" data_ref="input1" label="Column with reference allele if available" />
|
||||
</inputs>
|
||||
<outputs>
|
||||
<data format="interval" name="out_file1" />
|
||||
</outputs>
|
||||
<tests>
|
||||
<test>
|
||||
<param name='input1' value='dividePgSnp_input.pgSnp' ftype='interval' />
|
||||
<param name='ref_column' value='1' />
|
||||
<output name="output" file="dividePgSnp_output.txt" />
|
||||
</test>
|
||||
</tests>
|
||||
|
||||
<help>
|
||||
**What it does**
|
||||
|
||||
This separates the alleles from a pgSnp formated file into separate columns,
|
||||
as well as the frequency and scores that go with the alleles. It will skip
|
||||
any positions with more than 2 alleles. If only a single allele is given "N"
|
||||
will be used for the second with frequency and score of zero. If a column
|
||||
other than the first column is chosen for the reference allele, the value in
|
||||
that column will be used in place of the "N" for single alleles.
|
||||
|
||||
-----
|
||||
|
||||
**Examples**
|
||||
|
||||
- input pgSnp file::
|
||||
|
||||
chr1 256 257 A/C 2 2,4 10,20
|
||||
chr1 56100 56101 A 1 5 30
|
||||
chr1 77052 77053 A/G 2 3,2 40,50
|
||||
chr1 110904 110905 A 1 3 60
|
||||
etc.
|
||||
|
||||
- output::
|
||||
|
||||
chr1 256 257 A 2 10 C 4 20
|
||||
chr1 56100 56101 A 5 30 N 0 0
|
||||
chr1 77052 77053 A 3 40 G 2 50
|
||||
chr1 110904 110905 A 3 60 N 0 0
|
||||
etc.
|
||||
|
||||
</help>
|
||||
</tool>
|
||||
+89
@@ -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 (<FH>) {
|
||||
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 "<html><head><title>g:Profiler link</title></head><body>\n";
|
||||
print FH '<form method="POST" action="http://biit.cs.ut.ee/gprofiler/index.cgi">';
|
||||
foreach my $k (keys %params) {
|
||||
print FH "<input type='hidden' name='$k' value='$params{$k}'>\n";
|
||||
}
|
||||
print FH '<input type="Submit" name="foo" value="Send to g:Profiler">';
|
||||
print FH '</form></body></html>', "\n";
|
||||
close FH or die "Couldn't close $out, $!\n";
|
||||
|
||||
#also do link that prints text that could be pulled back into Galaxy?
|
||||
exit;
|
||||
+8
-4
@@ -2,13 +2,17 @@
|
||||
<description>tools for functional profiling of gene lists</description>
|
||||
|
||||
<command interpreter="perl">
|
||||
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}
|
||||
</command>
|
||||
|
||||
<inputs>
|
||||
<param name="input" type="data" format="tabular" label="Dataset" />
|
||||
<param name="numerical_column" type="data_column" data_ref="input" numerical="True" label="Column with identifiers" />
|
||||
<param name="type" label="Identifier type" type="select">
|
||||
<param name="genes" type="data_column" data_ref="input" label="Column with identifiers" />
|
||||
<param name="region" type="select" label="Or use genomic intervals">
|
||||
<option value="0" selected="true">No</option>
|
||||
<option value="1">Yes</option>
|
||||
</param>
|
||||
<param name="type" label="Identifier type if numeric" type="select">
|
||||
<option value="ENTREZGENE_ACC" selected="true">Entrez Gene Acc</option>
|
||||
<option value="MIM_MORBID">OMIM Morbid Map</option>
|
||||
<option value="MIM_GENE">OMIM Gene ID</option>
|
||||
@@ -31,7 +35,7 @@
|
||||
<tests>
|
||||
<test>
|
||||
<param name="input" ftype="tabular" value="linkToGProfile.tabular" />
|
||||
<param name="numerical_column" value="2" />
|
||||
<param name="genes" value="2" />
|
||||
<param name="type" value="ENTREZGENE_ACC" />
|
||||
<output name="out_file1" file="linkToGProfile_1.out" />
|
||||
</test>
|
||||
Executable
+116
@@ -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 (<FH>) {
|
||||
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;
|
||||
@@ -0,0 +1,78 @@
|
||||
<tool id="vcf2pgSnp" name="Convert VCF" hidden="false">
|
||||
<description>to pgSnp</description>
|
||||
<command interpreter="perl">
|
||||
#if $inType.how == "all" #vcf2pgSnp.pl all $input1 > $out_file1
|
||||
#else if $inType.how == "every" #vcf2pgSnpMult.pl $input1 > $out_file1
|
||||
#else #vcf2pgSnp.pl $inType.ind_column $input1 > $out_file1
|
||||
#end if
|
||||
</command>
|
||||
<inputs>
|
||||
<param format="vcf" name="input1" type="data" label="VCF file" />
|
||||
<conditional name="inType">
|
||||
<param name="how" type="select" label="How to treat individuals">
|
||||
<option value="all">Group all as a population</option>
|
||||
<option value="one">Do just one column</option>
|
||||
<option value="every">Convert each column separately</option>
|
||||
</param>
|
||||
<when value="one">
|
||||
<param name="ind_column" type="data_column" data_ref="input1" label="Column to convert" value="10" />
|
||||
</when>
|
||||
<when value="all">
|
||||
<!-- do nothing -->
|
||||
</when>
|
||||
<when value="every">
|
||||
<!-- do nothing -->
|
||||
</when>
|
||||
</conditional>
|
||||
</inputs>
|
||||
<outputs>
|
||||
<data format="interval" name="out_file1" />
|
||||
</outputs>
|
||||
<tests>
|
||||
<test>
|
||||
<param name="input1" value="vcf2pgSnp_input.vcf" ftype="vcf" />
|
||||
<param name="how" value="all" />
|
||||
<output name="output" file="vcf2pgSnp_output.pgSnp" />
|
||||
</test>
|
||||
<test>
|
||||
<param name="input1" value="vcf2pgSnp_input.vcf" ftype="vcf" />
|
||||
<param name="how" value="every" />
|
||||
<output name="output" file="vcf2pgSnp_output2.pgSnp" />
|
||||
</test>
|
||||
</tests>
|
||||
|
||||
<help>
|
||||
**What it does**
|
||||
|
||||
This converts a VCF file to a pgSnp file with the frequency counts being
|
||||
chromosome counts. If there is more than 1 column of SNP data it will either
|
||||
accumulate all columns as a population, convert the column indicated
|
||||
to pgSnp, or add pgSnp allele data columns for each of the individuals.
|
||||
|
||||
-----
|
||||
|
||||
**Examples**
|
||||
|
||||
- input::
|
||||
|
||||
1 13327 rs144762171 G C 100 PASS VT=SNP;SNPSOURCE=LOWCOV GT:DS:GL 0|0:0.000:-0.03,-1.11,-5.00 0|1:1.000:-1.97,-0.01,-2.51 0|0:0.050:-0.01,-1.69,-5.00 0|0:0.100:-0.48,-0.48,-0.48
|
||||
1 13980 rs151276478 T C 100 PASS VT=SNP;SNPSOURCE=LOWCOV GT:DS:GL 0|0:0.100:-0.48,-0.48,-0.48 0|1:0.950:-0.48,-0.48,-0.48 0|0:0.050:-0.48,-0.48,-0.48 0|0:0.050:-0.48,-0.48,-0.48
|
||||
1 30923 rs140337953 G T 100 PASS VT=SNP;SNPSOURCE=LOWCOV GT:DS:GL 1|1:1.950:-5.00,-0.61,-0.12 0|0:0.450:-0.10,-0.69,-2.81 0|0:0.450:-0.11,-0.64,-3.49 1|1:1.500:-0.48,-0.48,-0.48
|
||||
etc.
|
||||
|
||||
- output as a population::
|
||||
|
||||
chr1 13326 13327 G/C 2 7,1 0,0
|
||||
chr1 13979 13980 T/C 2 7,1 0,0
|
||||
chr1 30922 30923 G/T 2 4,4 0,0
|
||||
etc.
|
||||
|
||||
- output for each column separately::
|
||||
|
||||
chr1 13326 13327 G 1 2 0 G/C 2 1,1 0,0 G 1 2 0 G 1 2 0
|
||||
chr1 13979 13980 T 1 2 0 T/C 2 1,1 0,0 T 1 2 0 T 1 2 0
|
||||
chr1 30922 30923 T 1 2 0 G 1 2 0 G 1 2 0 T 1 2 0
|
||||
etc.
|
||||
|
||||
</help>
|
||||
</tool>
|
||||
Executable
+81
@@ -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 (<FH>) {
|
||||
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;
|
||||
Reference in New Issue
Block a user