opensource/galaxy
3
0
Fork 0
mirror of https://github.com/galaxyproject/galaxy.git synced 2026-09-24 16:30:27 +08:00
Code Issues Packages Projects Releases Wiki Activity

changes to aaChanges tool

Browse Source
This commit is contained in:
Richard Burhans
2012-03-22 13:19:31 -04:00
parent ff1b05d283
commit 77f391c6ae
2 changed files with 91 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tools/evolution/codingSnps.pl
+50 -9
View File
@@ -7,11 +7,15 @@ use strict;
# and a gene bed file with cds start and stop.
# It then checks for changes in coding regions, reporting
# those that cause a frameshift or substitution in the amino acid.
# Output columns:
# chrom, start, end, allele as given (amb code translated)
# Gene ID from genes file, ref amino acid:variant amino acids,
# codon number, (in strand of gene)ref nt, refCodon:variantCodons
#########################################################################
my $seqFlag = "2bit"; #flag to set sequence type 2bit|nib
if (!@ARGV or scalar @ARGV < 3) {
print "Usage: codingSnps.pl snps.bed genes.bed (/dir/*$seqFlag|Galaxy build= loc=) [chr=# start=# end=# snp=# keepColumns=1] > codingSnps.txt\n";
print "Usage: codingSnps.pl snps.bed genes.bed (/dir/*$seqFlag|Galaxy build= loc=) [chr=# start=# end=# snp=# strand=#|-|+ keepColumns=1 synon=1 unique=1] > codingSnps.txt\n";
exit;
}
my $uniq = 0; #flag for whether want uniq positions
@@ -25,6 +29,7 @@ my $col0 = 0; #bed like columns in default positions
my $col1 = 1;
my $col2 = 2;
my $col3 = 3;
my $strand = -1;
#column positions 1 based coming in (for Galaxy)
foreach (@ARGV) {
if (/chr=(\d+)/) { $col0 = $1 -1; }
@@ -32,6 +37,10 @@ foreach (@ARGV) {
elsif (/end=(\d+)/) { $col2 = $1 -1; }
elsif (/snp=(\d+)/) { $col3 = $1 -1; }
elsif (/keepColumns=1/) { $keep = 1; }
elsif (/synon=1/) { $syn = 1; }
elsif (/unique=1/) { $uniq = 1; }
elsif (/strand=(\d+)/) { $strand = $1 -1; } #0 based column
elsif (/strand=-/) { $strand = -99; } #special case of all minus
}
if ($col0 < 0 || $col1 < 0 || $col2 < 0 || $col3 < 0) {
print STDERR "ERROR column numbers are given with origin 1\n";
@@ -42,6 +51,7 @@ my %chrSt; #index in array where each chrom starts
my %codon; #hash of codon amino acid conversions
my $ends = 0; #ends vs sizes in bed 11 position, starts relative to chrom
my $ignoreN = 1; #skip N
my $origAll; #alleles from input file (before changes for strand)
my %amb = (
"R" => "A/G",
@@ -90,6 +100,10 @@ while(<FH>) {
print STDERR "ERROR file has fewer columns than requested, requested columns (0 based) $col0 $col1 $col2 $col3, file has $size\n";
exit 1;
}
if ($strand >= 0 && $strand > $size) {
print STDERR "ERROR file has fewer columns than requested, requested strand in $strand (0 based), file has $size\n";
exit 1;
}
if ($s[$col1] =~ /\D/) {
print STDERR "ERROR the start point must be an integer not $s[$col1]\n";
exit 1;
@@ -100,6 +114,11 @@ while(<FH>) {
}
if ($s[$col3] eq 'N' && $ignoreN) { next; }
if (exists $amb{$s[$col3]}) { $s[$col3] = $amb{$s[$col3]}; }
if (($strand >= 0 && $s[$strand] eq '-') or $strand == -99) {
#reverse complement nts
$origAll = $s[$col3];
$s[$col3] = reverseCompAlleles($s[$col3]);
}else { undef $origAll }
if (!@g && exists $chrSt{$s[$col0]}) { #need to fetch first gene row
$i = $chrSt{$s[$col0]};
@g = split(/\t/, $genes[$i]);
@@ -287,36 +306,45 @@ sub processSnp {
}
#now compute amino acids
my $oldaa = getaa($oldnts);
my $codon = "$oldnts:";
my @newaa;
my $change = 0; #flag for if there is a change
foreach my $v (@newnts) {
my $t = getaa($v);
if ($t ne $oldaa) { $change = 1; }
push(@newaa, $t);
push(@newaa, "$t");
$codon .= "$v/";
}
$codon =~ s/\/$//;
if (!$change && $syn) {
if (!$keep) {
print "$sref->[$col0]\t$sref->[$col1]\t$sref->[$col2]\t$sref->[$col3]";
print "\t$gref->[3]\t$oldaa:", join("/", @newaa), "\n";
print "\t$gref->[3]\t$oldaa:", join("/", @newaa), "\t$cdNum\t$changed\t$codon\n";
}else {
my @s = @{$sref};
print join("\t", @s),
"\t$gref->[3]\t$oldaa:", join("/", @newaa), "\n";
"\t$gref->[3]\t$oldaa:", join("/", @newaa), "\t$cdNum\t$changed\t$codon\n";
}
$done{"$sref->[$col0] $sref->[$col1] $sref->[$col2]"}++;
return;
}elsif ($syn) { return; } #only want synonymous changes
if (!$change) { return; } #no change in amino acids
if (!$keep) {
print "$sref->[$col0]\t$sref->[$col1]\t$sref->[$col2]\t$sref->[$col3]";
if ($gref->[5] eq '-') { $changed = compl($changed); } #use plus for ref
my $a = $sref->[$col3];
if (($strand >= 0 && $origAll) or $strand == -99) { $a = $origAll; }
print "$sref->[$col0]\t$sref->[$col1]\t$sref->[$col2]\t$a";
#my $minus = $changed; #in case minus strand and change back
#if ($gref->[5] eq '-') { $changed = compl($changed); } #use plus for ref
if (!$changed) { return; } #skip this one
print "\t$gref->[3]\t$oldaa:", join("/", @newaa), "\t$cdNum\t$changed\n";
print "\t$gref->[3]\t$oldaa:", join("/", @newaa), "\t$cdNum\t$changed\t$codon\n";
}else {
my @s = @{$sref};
if (($strand >= 0 && $origAll) or $strand == -99) { $s[$col3] = $origAll; }
print join("\t", @s);
if ($gref->[5] eq '-') { $changed = compl($changed); } #use plus for ref
#my $minus = $changed; #in case minus strand and change back
#if ($gref->[5] eq '-') { $changed = compl($changed); } #use plus for ref
if (!$changed) { return; } #skip this one
print "\t$gref->[3]\t$oldaa:", join("/", @newaa), "\t$cdNum\t$changed\n";
print "\t$gref->[3]\t$oldaa:", join("/", @newaa), "\t$cdNum\t$changed\t$codon\n";
}
$done{"$sref->[$col0] $sref->[$col1] $sref->[$col2]"}++;
}
@@ -418,6 +446,19 @@ sub compl {
return $comp;
}
sub reverseCompAlleles {
my $all = shift;
my @nt = split(/\//, $all);
my $rv = '';
foreach my $n (@nt) {
$n = reverse(split(/ */, $n)); #needed for indels
$n = compl($n);
$rv .= "$n/";
}
$rv =~ s/\/$//;
return $rv;
}
sub getaa {
my $nts = shift; #in multiples of 3
my $aa = '';
tools/evolution/codingSnps.xml
+41 -12
View File
@@ -2,7 +2,7 @@
<description>amino-acid changes caused by a set of SNPs</description>
<command interpreter="perl">
codingSnps.pl $input1 $input2 Galaxy build=${input1.metadata.dbkey} loc=${GALAXY_DATA_INDEX_DIR}/codingSnps.loc chr=${input1.metadata.chromCol} start=${input1.metadata.startCol} end=${input1.metadata.endCol} snp=$col1 keepColumns=$keep > $out_file1
codingSnps.pl $input1 $input2 Galaxy build=${input1.metadata.dbkey} loc=${GALAXY_DATA_INDEX_DIR}/codingSnps.loc chr=${input1.metadata.chromCol} start=${input1.metadata.startCol} end=${input1.metadata.endCol} snp=$col1 keepColumns=$keep strand=${strand_source.strand_col} > $out_file1
</command>
<inputs>
@@ -17,6 +17,22 @@
<option value="0" selected="true">No</option>
<option value="1">Yes</option>
</param>
<conditional name="strand_source">
<param name="strand_choice" type="select" label="Strand info">
<option value="data_column">a column in the dataset</option>
<option value="all_pos" selected="true">all on sense/forward/+ strand</option>
<option value="all_neg">all on antisense/reverse/- strand</option>
</param>
<when value="data_column">
<param name="strand_col" type="data_column" data_ref="input1" label="Column with strand"/>
</when>
<when value="all_pos">
<param name="strand_col" type="hidden" value="+"/>
</when>
<when value="all_neg">
<param name="strand_col" type="hidden" value="-"/>
</when>
</conditional>
</inputs>
<outputs>
@@ -36,14 +52,26 @@
<param name="input1" ftype="interval" value="codingSnps_input1.interval" dbkey="hg18" />
<param name="col1" value="6" />
<param name="input2" ftype="interval" value="codingSnps_inputGenes1.bed" dbkey="hg18" />
<param name="strand_choice" value="all_pos" />
<param name="strand_col" value="+" />
<output name="output" file="codingSnps_output1.interval" />
</test>
<test>
<param name="input1" ftype="interval" value="codingSnps_input2.interval" dbkey="hg18" />
<param name="input2" ftype="interval" value="codingSnps_inputGenes2.bed" dbkey="hg18" />
<param name="col1" value="4" />
<param name="strand_choice" value="all_pos" />
<param name="strand_col" value="+" />
<output name="output" file="codingSnps_output2.interval" />
</test>
<test>
<param name="input1" ftype="interval" value="codingSnps_input2.interval" dbkey="hg18" />
<param name="input2" ftype="interval" value="codingSnps_inputGenes2.bed" dbkey="hg18" />
<param name="col1" value="4" />
<param name="strand_choice" value="all_neg" />
<param name="strand_col" value="-" />
<output name="output" file="codingSnps_output3.interval" />
</test>
</tests>
<help>
@@ -79,8 +107,9 @@ to be used for defining the coding regions. This file must be a BED file with
the first 12 columns standard BED columns. The output is the same as the
first input file with
several columns added: the name field from the line of the gene input file
used, the amino acids, the codon number, and the reference nucleotide that
changed in the amino acid.
used, the amino acids, the codon number, the reference nucleotide that
changed in the amino acid (in the same strand as the gene), and the codons
that go with the amino acids.
The amino acids are listed with the reference amino acid first, then a colon,
and then the amino acids for the alleles. If a SNP is not in a coding region
or is synonymous then it is not included in the output file.
@@ -128,15 +157,15 @@ or is synonymous then it is not included in the output file.
- output file, showing non-synonymous substitutions in coding regions::
chr22 15825725 15825726 G/T uc002zlw.1 Gln:Pro/Gln 469 T
chr22 15827035 15827036 G uc002zlw.1 Glu:Asp 414 C
chr22 15827135 15827136 C/G uc002zlw.1 Gly:Gly/Ala 381 C
chr22 15830928 15830929 A/G uc002zlw.1 Ala:Ser/Pro 281 C
chr22 15830951 15830952 G uc002zlw.1 Leu:Pro 273 A
chr22 15830955 15830956 C/T uc002zlw.1 Ser:Gly/Ser 272 T
chr22 15848885 15848886 C/T uc002zlw.1 Ser:Trp/Stop 217 G
chr22 15848885 15848886 C/T uc010gqs.1 Ser:Trp/Stop 200 G
chr22 15849048 15849049 A/C uc002zlw.1 Gly:Stop/Gly 163 C
chr22 15825725 15825726 G/T uc002zlw.1 Gln:Pro/Gln 469 A CAA:CCA/CAA
chr22 15827035 15827036 G uc002zlw.1 Glu:Asp 414 G GAG:GAC
chr22 15827135 15827136 C/G uc002zlw.1 Gly:Gly/Ala 381 G GGT:GGT/GCT
chr22 15830928 15830929 A/G uc002zlw.1 Ala:Ser/Pro 281 G GCA:TCA/CCA
chr22 15830951 15830952 G uc002zlw.1 Leu:Pro 273 T CTT:CCT
chr22 15830955 15830956 C/T uc002zlw.1 Ser:Gly/Ser 272 A AGC:GGC/AGC
chr22 15848885 15848886 C/T uc002zlw.1 Ser:Trp/Stop 217 C TCG:TGG/TAG
chr22 15848885 15848886 C/T uc010gqs.1 Ser:Trp/Stop 200 C TCG:TGG/TAG
chr22 15849048 15849049 A/C uc002zlw.1 Gly:Stop/Gly 163 G GGA:TGA/GGA
etc.
</help>