From 89bb493d3c0f1691791b323fcf9b5de3be470380 Mon Sep 17 00:00:00 2001 From: Guruprasad Anada Date: Thu, 13 Mar 2008 16:21:39 +0000 Subject: [PATCH] Adding new tool to estimate insertion and deletion rates from 3-way alignments. Also, modified 'Fetch Indels' tool. --- tool_conf.xml.sample | 1 + tools/regVariation/getIndelRates_3way.py | 126 +++++++++++++++++++++ tools/regVariation/getIndelRates_3way.xml | 62 ++++++++++ tools/regVariation/getIndels_2way.xml | 4 +- tools/regVariation/getIndels_3way.xml | 4 +- tools/regVariation/parseMAF_smallIndels.pl | 45 ++++---- 6 files changed, 215 insertions(+), 27 deletions(-) create mode 100755 tools/regVariation/getIndelRates_3way.py create mode 100644 tools/regVariation/getIndelRates_3way.xml diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample index 7ebfbf8516a..dc57f8a62f4 100644 --- a/tool_conf.xml.sample +++ b/tool_conf.xml.sample @@ -125,6 +125,7 @@ +
diff --git a/tools/regVariation/getIndelRates_3way.py b/tools/regVariation/getIndelRates_3way.py new file mode 100755 index 00000000000..021b6612619 --- /dev/null +++ b/tools/regVariation/getIndelRates_3way.py @@ -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() + \ No newline at end of file diff --git a/tools/regVariation/getIndelRates_3way.xml b/tools/regVariation/getIndelRates_3way.xml new file mode 100644 index 00000000000..d8e79d4434d --- /dev/null +++ b/tools/regVariation/getIndelRates_3way.xml @@ -0,0 +1,62 @@ + + for 3-way alignments + + getIndelRates_3way.py $input1 $out_file1 $winsize $species + + + + + + + + + + + + + + + + + + + + + + + + + + + + +.. 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. + + + + + diff --git a/tools/regVariation/getIndels_2way.xml b/tools/regVariation/getIndels_2way.xml index 13e6133b515..9fcb12b34c6 100644 --- a/tools/regVariation/getIndels_2way.xml +++ b/tools/regVariation/getIndels_2way.xml @@ -1,5 +1,5 @@ - - for pairwise alignments + + from pairwise alignments getIndels.py $input1 $out_file1 diff --git a/tools/regVariation/getIndels_3way.xml b/tools/regVariation/getIndels_3way.xml index 54261b5c19a..2f131d715dd 100644 --- a/tools/regVariation/getIndels_3way.xml +++ b/tools/regVariation/getIndels_3way.xml @@ -1,5 +1,5 @@ - - for 3-way alignments + + from 3-way alignments parseMAF_smallIndels.pl $input1 $out_file1 $outgroup diff --git a/tools/regVariation/parseMAF_smallIndels.pl b/tools/regVariation/parseMAF_smallIndels.pl index 75da3d242c9..06a10403c4f 100644 --- a/tools/regVariation/parseMAF_smallIndels.pl +++ b/tools/regVariation/parseMAF_smallIndels.pl @@ -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) = @_;