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

Adding new tool to estimate insertion and deletion rates from 3-way alignments.

Browse Source 
Also, modified 'Fetch Indels' tool.
This commit is contained in:
Guruprasad Anada
2008-03-13 16:21:39 +00:00
parent aa91247bc7
commit 89bb493d3c
6 changed files with 215 additions and 27 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tool_conf.xml.sample
+1
View File
@@ -125,6 +125,7 @@
<tool file="regVariation/maf_cpg_filter.xml" />
<tool file="regVariation/getIndels_2way.xml" />
<tool file="regVariation/getIndels_3way.xml" />
<tool file="regVariation/getIndelRates_3way.xml" />
</section>
<section name="Evolution: HyPhy" id="hyphy">
<tool file="hyphy/hyphy_branch_lengths_wrapper.xml" />
tools/regVariation/getIndelRates_3way.py
Executable
+126
View File
@@ -0,0 +1,126 @@
#!/usr/bin/env python2.4
#Guruprasad Ananda
import sys, os, tempfile, string
fout = open(sys.argv[2],'w')
winsize = int(sys.argv[3])
species_ind = int(sys.argv[4])
def stop_err(msg):
sys.stderr.write(msg)
sys.exit()
def rate_estimator(win, blk_lines, wstart, wend, wspecies):
inserts = 0.0
deletes = 0.0
ilengths = {} #dict containing lengths of blocks(without gaps) having insertion in wspecies
dlengths = {} #dict containing lengths of blocks(without gaps) having deletion in wspecies
prev_bnum = -1
for bline in blk_lines:
items = bline.split('\t')
bnum = int(items[0])
bevent = items[1]
if not(bevent.startswith(wspecies)):
continue
if bevent.endswith('insert'):
inserts += 1
#Add lengths only if the insert belongs to a new alignment block
if not(ilengths.has_key(bnum)):
ilengths[bnum] = int(items[species_ind].split(':')[1])
#prev_bnum = bnum
elif bevent.endswith('delete'):
deletes += 1
#Add lengths only if the delete belongs to a new alignment block
if not(dlengths.has_key(bnum)):
dlengths[bnum] = int(items[species_ind].split(':')[1])
#prev_bnum = bnum
try:
total_ilength = sum(ilengths.values())
irate = inserts/total_ilength
except:
irate = 0
try:
total_dlength = sum(dlengths.values())
drate = deletes/total_dlength
except:
drate = 0
print >>fout, "%s\t%s\t%s\t%s\t%.2e\t%.2e" %(win, wspecies, wstart, wend, irate , drate)
def main():
infile = sys.argv[1]
for i, line in enumerate( file ( infile )):
line = line.rstrip('\r\n')
if len( line )>0 and not line.startswith( '#' ):
elems = line.split( '\t' )
break
if i == 30:
break # Hopefully we'll never get here...
if len( elems ) != 15:
stop_err( "This tool only works on tabular data output by 'Fetch Indels from 3-way alignments' tool. The data in your input dataset is either missing or not formatted properly." )
wspecies = elems[species_ind].split(':')[0].split('.')[0]
fin = open(infile, 'r')
skipped = 0
blk=0
win=0
linestr=""
sorted_infile = tempfile.NamedTemporaryFile()
cmdline = "sort -n -k"+str(species_ind+2)+" -o "+sorted_infile.name+" "+infile
try:
os.system(cmdline)
except:
stop_err("Encountered error while sorting the input file.")
print >>fout, "#Window\tSpecies\tWindow_Start\tWindow_End\tInsertion_Rate\tDeletion_Rate"
for line in sorted_infile.readlines():
line = line.strip("\r\n")
if not(line) or line == "":
continue
elems = line.split('\t')
try:
assert int(elems[0])
assert len(elems) == 15
except Exception, eon:
continue
if not(elems[1].startswith(wspecies)): #Event doesn't belong to the selected species
continue
try:
assert wstart
except NameError:
wstart = int(elems[species_ind+1]) - int(elems[species_ind+1])%winsize + 1
wend = wstart + winsize
lstart = int(elems[species_ind + 1])
if lstart in range(wstart,wend+1):
linestr += line.strip()
linestr += "\n"
else:
try:
win += 1
blk_lines = linestr.strip().split("\n")
rate_estimator(str(win), blk_lines, str(wstart), str(wend), wspecies)
linestr = ""
except:
skipped += 1
pass
linestr=line.strip()+"\n"
wstart = int(elems[species_ind+1]) - int(elems[species_ind+1])%winsize + 1
wend = wstart + winsize
if linestr != "":
try:
win += 1
blk_lines = linestr.strip().split("\n")
rate_estimator(str(win), blk_lines, str(wstart), str(wend), wspecies)
except:
skipped += 1
pass
if skipped:
print "Skipped %s windows as invalid." %(skipped)
if __name__ == "__main__":
main()
tools/regVariation/getIndelRates_3way.xml
+62
View File
@@ -0,0 +1,62 @@
<tool id="getIndelRates_3way" name="Estimate Indel Rates" version="1.0.0">
<description> for 3-way alignments</description>
<command interpreter="python">
getIndelRates_3way.py $input1 $out_file1 $winsize $species
</command>
<inputs>
<page>
<param format="tabular" name="input1" type="data" label="Select data"/>
<param name="winsize" size="10" type="integer" value="1000" label="Estimate rates in windows of size" />
<param name="species" type="select" label="and corresponding to co-ordinates of" multiple="false">
<option value="3">Species 1 (Ingroup 1)</option>
<option value="7">Species 2 (Ingroup 2)</option>
<option value="11">Species 3 (Outgroup)</option>
</param>
<!--
<conditional name="region">
<param name="type" type="select" label="Estimate rates per" multiple="false">
<option value="align">Alignment block</option>
<option value="win">Window</option>
</param>
<when value="win">
<param name="winsize" size="10" type="integer" value="1000" label="of size" />
</when>
<when value="align" />
</conditional>
-->
</page>
</inputs>
<outputs>
<data format="tabular" name="out_file1" metadata_source="input1"/>
</outputs>
<tests>
<test>
<param name="input1" value="indels_3way.tabular"/>
<param name="winsize" value="1000"/>
<param name="species" value="11"/>
<output name="out_file1" file="indelrates_3way.tabular"/>
</test>
</tests>
<help>
.. class:: infomark
**What it does**
This tool estimates the insertion and deletion rates for alignments in a window of specified size.
-----
.. class:: warningmark
**Note**
Any block/s not containing exactly 3 species will be omitted.
</help>
</tool>
tools/regVariation/getIndels_2way.xml
+2 -2
View File
@@ -1,5 +1,5 @@
<tool id="getIndels_2way" name="Get Indels">
<description> for pairwise alignments</description>
<tool id="getIndels_2way" name="Fetch Indels">
<description> from pairwise alignments</description>
<command interpreter="python">
getIndels.py $input1 $out_file1
</command>
tools/regVariation/getIndels_3way.xml
+2 -2
View File
@@ -1,5 +1,5 @@
<tool id="getIndels_3way" name="Get Indel rates">
<description> for 3-way alignments</description>
<tool id="getIndels_3way" name="Fetch Indels" version="1.0.1">
<description> from 3-way alignments</description>
<command interpreter="perl">
parseMAF_smallIndels.pl $input1 $out_file1 $outgroup
</command>
tools/regVariation/parseMAF_smallIndels.pl
+22 -23
View File
@@ -214,7 +214,7 @@ sub get_indels_within_block{
$line1 =~ s/\s+/\t/g;
@line1 = split(/\t/, $line1);
$end1 =($line1[2]+$line1[3]-1);
$seq1 = $line1[1];
$seq1 = $line1[1].":".$line1[3];
$ingroup1 = (split(/\./, $seq1))[0];
$start1 = $line1[2];
$align_length1 = $line1[3];
@@ -231,7 +231,7 @@ sub get_indels_within_block{
$line1 =~ s/\s+/\t/g;
@line1 = split(/\t/, $line1);
$end3 =($line1[2]+$line1[3]-1);
$seq3 = $line1[1];
$seq3 = $line1[1].":".$line1[3];
$start3 = $line1[2];
$align_length3 = $line1[3];
$orient3 = $line1[4];
@@ -248,7 +248,7 @@ sub get_indels_within_block{
$line2 =~ s/\s+/\t/g;
@line2 = split(/\t/, $line2);
$end2 =($line2[2]+$line2[3]-1);
$seq2 = $line2[1];
$seq2 = $line2[1].":".$line2[3];
$ingroup2 = (split(/\./, $seq2))[0];
$start2 = $line2[2];
$align_length2 = $line2[3];
@@ -265,7 +265,7 @@ sub get_indels_within_block{
$line2 =~ s/\s+/\t/g;
@line2 = split(/\t/, $line2);
$end3 =($line2[2]+$line2[3]-1);
$seq3 = $line2[1];
$seq3 = $line2[1].":".$line2[3];
$start3 = $line2[2];
$align_length3 = $line2[3];
$orient3 = $line2[4];
@@ -281,7 +281,7 @@ sub get_indels_within_block{
$line2 =~ s/\s+/\t/g;
@line2 = split(/\t/, $line2);
$end1 =($line2[2]+$line2[3]-1);
$seq1 = $line2[1];
$seq1 = $line2[1].":".$line2[3];
$ingroup1 = (split(/\./, $seq1))[0];
$start1 = $line2[2];
$align_length1 = $line2[3];
@@ -299,7 +299,7 @@ sub get_indels_within_block{
$line3 =~ s/\s+/\t/g;
@line3 = split(/\t/, $line3);
$end2 =($line3[2]+$line3[3]-1);
$seq2 = $line3[1];
$seq2 = $line3[1].":".$line3[3];
$ingroup2 = (split(/\./, $seq2))[0];
$start2 = $line3[2];
$align_length2 = $line3[3];
@@ -316,7 +316,7 @@ sub get_indels_within_block{
$line3 =~ s/\s+/\t/g;
@line3 = split(/\t/, $line3);
$end3 =($line3[2]+$line3[3]-1);
$seq3 = $line3[1];
$seq3 = $line3[1].":".$line3[3];
$start3 = $line3[2];
$align_length3 = $line3[3];
$orient3 = $line3[4];
@@ -339,14 +339,14 @@ sub get_indels_within_block{
$coord1 = $start1_plus;
$coord2 = $start2_plus;
$coord3 = $start3_plus;
for (my $position = 0; $position < $test1; $position++) {
my $indelType = "";
my $indel_line = "";
# seq1 deletes
if ((substr($sequence1,$position,1) eq "-")
&& (substr($sequence2,$position,1) ne "-")
&& (substr($sequence3,$position,1) ne "-")){
&& (substr($sequence2,$position,1) !~ m/[-*\#$?^@]/)
&& (substr($sequence3,$position,1) !~ m/[-*\#$?^@]/)){
$ABC = join("",($ABC,"X"));
$indelType = $seq1."_delete";
@@ -357,9 +357,9 @@ sub get_indels_within_block{
$coord2++; $coord3++;
}
# seq2 deletes
elsif ((substr($sequence1,$position,1) ne "-")
elsif ((substr($sequence1,$position,1) !~ m/[-*\#$?^@]/)
&& (substr($sequence2,$position,1) eq "-")
&& (substr($sequence3,$position,1) ne "-")){
&& (substr($sequence3,$position,1) !~ m/[-*\$?^]/)){
$ABC = join("",($ABC,"Y"));
$indelType = $seq2."_delete";
#print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
@@ -371,7 +371,7 @@ sub get_indels_within_block{
}
# seq1 inserts
elsif ((substr($sequence1,$position,1) ne "-")
elsif ((substr($sequence1,$position,1) !~ m/[-*\#$?^@]/)
&& (substr($sequence2,$position,1) eq "-")
&& (substr($sequence3,$position,1) eq "-")){
$ABC = join("",($ABC,"Z"));
@@ -384,7 +384,7 @@ sub get_indels_within_block{
}
# seq2 inserts
elsif ((substr($sequence1,$position,1) eq "-")
&& (substr($sequence2,$position,1) ne "-")
&& (substr($sequence2,$position,1) !~ m/[-*\#$?^@]/)
&& (substr($sequence3,$position,1) eq "-")){
$ABC = join("",($ABC,"W"));
$indelType = $seq2."_insert";
@@ -395,8 +395,8 @@ sub get_indels_within_block{
$coord2++;
}
# seq3 deletes
elsif ((substr($sequence1,$position,1) ne "-")
&& (substr($sequence2,$position,1) ne "-")
elsif ((substr($sequence1,$position,1) !~ m/[-*\#$?^@]/)
&& (substr($sequence2,$position,1) !~ m/[-*\#$?^@]/)
&& (substr($sequence3,$position,1) eq "-")){
$ABC = join("",($ABC,"S"));
$indelType = $seq3."_delete";
@@ -409,7 +409,7 @@ sub get_indels_within_block{
# seq3 inserts
elsif ((substr($sequence1,$position,1) eq "-")
&& (substr($sequence2,$position,1) eq "-")
&& (substr($sequence3,$position,1) ne "-")){
&& (substr($sequence3,$position,1) !~ m/[-*\#$?^@]/)){
$ABC = join("",($ABC,"T"));
$indelType = $seq3."_insert";
#print OFILE "$count\t$seq1\t$coord1\t$orient1\t$seq2\t$coord2\t$orient2\t$seq3\t$coord3\t$orient3\t$indelType\n";
@@ -426,7 +426,6 @@ sub get_indels_within_block{
@array_return=($seq1,$seq2,$seq3,$ABC);
return (@array_return);
}
# ignore pairwise cases for now, just count the number of blocks
elsif (scalar(@sequences) == 2){
@@ -516,7 +515,7 @@ sub get_starts_only{
chomp($line);
$line =~ s/^\s*//;
$line =~ s/\s+/\t/g;
my @line1 = split(/\t/, $line);
my @line1 = split(/\t/, $line);
$seq1 = $line1[1];
$coord1 = $line1[2];
$seq2 = $line1[4];
@@ -610,12 +609,12 @@ for ($counter6 = 0; $counter6 < @seq1_delete_startOnly; $counter6++){
# # if inserts, increase coords for the sequence inserted, other sequences give coords for 5' and 3' bases flanking the gap
# # for deletes, increase coords for other 2 sequences and the one deleted give coords for 5' and 3' bases flanking the gap
get_final_format(@final1);
get_final_format(@final2);
get_final_format(@final3);
get_final_format(@final4);
get_final_format(@final5);
get_final_format(@final6);
get_final_format(@final3);
get_final_format(@final4);
get_final_format(@final1);
get_final_format(@final2);
sub get_final_format{
my (@final) = @_;