annotate tools/regVariation/getIndelRates_3way.py @ 2:c2a356708570

Uploaded

1 #!/usr/bin/env python
2 #Guruprasad Ananda
3 
4 from galaxy import eggs
5 import pkg_resources
6 pkg_resources.require( "bx-python" )
7 
8 import sys, os, tempfile
9 import traceback
10 import fileinput
11 from warnings import warn
12 
13 from galaxy.tools.util.galaxyops import *
14 from bx.intervals.io import *
15 
16 from bx.intervals.operations import quicksect
17 
18 def stop_err(msg):
19     sys.stderr.write(msg)
20     sys.exit()
21     
22 def counter(node, start, end, sort_col):
23     global full, blk_len, blk_list
24     if node.start < start:
25         if node.right:
26             counter(node.right, start, end, sort_col)
27     elif start <= node.start <= end and start <= node.end <= end:
28         full += 1
29         if node.other[0] not in blk_list:
30             blk_list.append(node.other[0])
31             blk_len += int(node.other[sort_col+2])
32         if node.left and node.left.maxend > start:
33             counter(node.left, start, end, sort_col)
34         if node.right: 
35             counter(node.right, start, end, sort_col)
36     elif node.start > end:
37         if node.left: 
38             counter(node.left, start, end, sort_col)
39             
40 
41 infile = sys.argv[1]  
42 fout = open(sys.argv[2],'w')
43 int_file = sys.argv[3]
44 if int_file != "None": #User has specified an interval file
45     try:
46         fint = open(int_file, 'r')
47         dbkey_i = sys.argv[4]
48         chr_col_i, start_col_i, end_col_i, strand_col_i = parse_cols_arg( sys.argv[5] )
49     except:
50         stop_err("Unable to open input Interval file")
51         
52 def main():
53 
54     for i, line in enumerate( file ( infile )):
55         line = line.rstrip('\r\n')
56         if len( line )>0 and not line.startswith( '#' ):
57             elems = line.split( '\t' )
58             break
59         if i == 30:
60             break # Hopefully we'll never get here...
61         
62     if len( elems ) != 18:
63         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." )
64     
65     for i, line in enumerate( file ( infile )):
66         line = line.rstrip('\r\n')
67         elems = line.split('\t')
68         try:
69             assert int(elems[0])
70             assert len(elems) == 18
71             if int_file != "None":
72                 if dbkey_i not in elems[3] and  dbkey_i not in elems[8] and dbkey_i not in elems[13]:
73                     stop_err("The species build corresponding to your interval file is not present in the Indel file.") 
74                 if dbkey_i in elems[3]:
75                     sort_col = 4
76                 elif dbkey_i in elems[8]:
77                     sort_col = 9
78                 elif dbkey_i in elems[13]:
79                     sort_col = 14
80             else:
81                 species = []
82                 species.append( elems[3].split('.')[0] )
83                 species.append( elems[8].split('.')[0] )
84                 species.append( elems[13].split('.')[0] )
85                 sort_col = 0    #Based on block numbers
86             break
87         except:
88             continue
89         
90         
91     fin = open(infile, 'r')
92     skipped = 0
93     
94     if int_file == "None":
95         sorted_infile = tempfile.NamedTemporaryFile()
96         cmdline = "sort -n -k"+str(1)+" -o "+sorted_infile.name+" "+infile
97         try:
98             os.system(cmdline)
99         except:
00             stop_err("Encountered error while sorting the input file.")
01         print >>fout, "#Block\t%s_InsRate\t%s_InsRate\t%s_InsRate\t%s_DelRate\t%s_DelRate\t%s_DelRate" %(species[0],species[1],species[2],species[0],species[1],species[2])
02         prev_bnum = -1
03         sorted_infile.seek(0)
04         for line in sorted_infile.readlines():
05             line = line.rstrip('\r\n')
06             elems = line.split('\t')
07             try:
08                 assert int(elems[0])
09                 assert len(elems) == 18
10                 new_bnum = int(elems[0])
11                 if new_bnum != prev_bnum:
12                     if prev_bnum != -1:
13                         irate = []
14                         drate = []
15                         for i,elem in enumerate(inserts):
16                             try:
17                                 irate.append(str("%.2e" %(inserts[i]/blen[i])))
18                             except:
19                                 irate.append('0')
20                             try:
21                                 drate.append(str("%.2e" %(deletes[i]/blen[i])))
22                             except:
23                                 drate.append('0')
24                         print >>fout, "%s\t%s\t%s" %(prev_bnum, '\t'.join(irate) , '\t'.join(drate))
25                     inserts = [0.0, 0.0, 0.0]
26                     deletes = [0.0, 0.0, 0.0]
27                     blen = []
28                     blen.append( int(elems[6]) )
29                     blen.append( int(elems[11]) )
30                     blen.append( int(elems[16]) )
31                 line_sp = elems[1].split('.')[0]
32                 sp_ind = species.index(line_sp)
33                 if elems[1].endswith('insert'):
34                     inserts[sp_ind] += 1
35                 elif elems[1].endswith('delete'):
36                     deletes[sp_ind] += 1
37                 prev_bnum = new_bnum 
38             except Exception, ei:
39                 #print >>sys.stderr, ei
40                 continue
41         irate = []
42         drate = []
43         for i,elem in enumerate(inserts):
44             try:
45                 irate.append(str("%.2e" %(inserts[i]/blen[i])))
46             except:
47                 irate.append('0')
48             try:
49                 drate.append(str("%.2e" %(deletes[i]/blen[i])))
50             except:
51                 drate.append('0')
52         print >>fout, "%s\t%s\t%s" %(prev_bnum, '\t'.join(irate) , '\t'.join(drate))
53         sys.exit()
54     
55     
56     inf = open(infile, 'r')
57     start_met = False
58     end_met = False
59     sp_file = tempfile.NamedTemporaryFile()
60     for n, line in enumerate(inf):
61         line = line.rstrip('\r\n')
62         elems = line.split('\t')
63         try:
64             assert int(elems[0])
65             assert len(elems) == 18
66             if dbkey_i not in elems[1]: 
67                 if not(start_met):   
68                     continue
69                 else:
70                     sp_end = n
71                     break
72             else:
73                 print >>sp_file, line
74                 if not(start_met):
75                     start_met = True
76                     sp_start = n
77         except:
78             continue
79     
80     try:
81         assert sp_end
82     except:
83         sp_end = n+1
84     
85     sp_file.seek(0)
86     win = NiceReaderWrapper( fileinput.FileInput( int_file ),
87                                 chrom_col=chr_col_i,
88                                 start_col=start_col_i,
89                                 end_col=end_col_i,
90                                 strand_col=strand_col_i,
91                                 fix_strand=True)
92     
93     indel = NiceReaderWrapper( fileinput.FileInput( sp_file.name ),
94                                 chrom_col=1,
95                                 start_col=sort_col,
96                                 end_col=sort_col+1,
97                                 strand_col=-1,
98                                 fix_strand=True)
99     
00     indelTree = quicksect.IntervalTree()
01     for item in indel:
02         if type( item ) is GenomicInterval:
03             indelTree.insert( item, indel.linenum, item.fields )
04     result=[]
05     
06     global full, blk_len, blk_list
07     for interval in win:
08         if type( interval ) is Header:
09             pass
10         if type( interval ) is Comment:
11             pass
12         elif type( interval ) == GenomicInterval:
13             chrom = interval.chrom
14             start = int(interval.start)
15             end = int(interval.end)
16             if start > end: 
17                 warn( "Interval start after end!" )
18             ins_chr = "%s.%s_insert" %(dbkey_i,chrom)
19             del_chr = "%s.%s_delete" %(dbkey_i,chrom)
20             irate = 0
21             drate = 0
22             if ins_chr not in indelTree.chroms and del_chr not in indelTree.chroms:
23                 pass    
24             else:
25                 if ins_chr in indelTree.chroms:
26                     full = 0.0
27                     blk_len = 0
28                     blk_list = []
29                     root = indelTree.chroms[ins_chr]    #root node for the chrom insertion tree
30                     counter(root, start, end, sort_col)
31                     if blk_len:
32                         irate = full/blk_len
33                 
34                 if del_chr in indelTree.chroms:
35                     full = 0.0
36                     blk_len = 0
37                     blk_list = []
38                     root = indelTree.chroms[del_chr]    #root node for the chrom insertion tree
39                     counter(root, start, end, sort_col)
40                     if blk_len:
41                         drate = full/blk_len
42                 
43             interval.fields.append(str("%.2e" %irate))
44             interval.fields.append(str("%.2e" %drate))
45             print >>fout, "\t".join(interval.fields)
46             fout.flush()
47 
48 if __name__ == "__main__":
49     main()