annotate tools/regVariation/microsats_mutability.py @ 0:9071e359b9a3

Uploaded

1 #!/usr/bin/env python
2 #Guruprasad Ananda
3 """
4 This tool computes microsatellite mutability for the orthologous microsatellites fetched from  'Extract Orthologous Microsatellites from pair-wise alignments' tool.
5 """
6 from galaxy import eggs
7 import sys, string, re, commands, tempfile, os, fileinput
8 from galaxy.tools.util.galaxyops import *
9 from bx.intervals.io import *
10 from bx.intervals.operations import quicksect
11 
12 fout = open(sys.argv[2],'w')
13 p_group = int(sys.argv[3])        #primary "group-by" feature
14 p_bin_size = int(sys.argv[4])
15 s_group = int(sys.argv[5])        #sub-group by feature
16 s_bin_size = int(sys.argv[6])
17 mono_threshold = 9
18 non_mono_threshold = 4
19 p_group_cols = [p_group, p_group+7]
20 s_group_cols = [s_group, s_group+7]
21 num_generations = int(sys.argv[7])
22 region = sys.argv[8] 
23 int_file = sys.argv[9]
24 if int_file != "None": #User has specified an interval file
25     try:
26         fint = open(int_file, 'r')
27         dbkey_i = sys.argv[10]
28         chr_col_i, start_col_i, end_col_i, strand_col_i = parse_cols_arg( sys.argv[11] )
29     except:
30         stop_err("Unable to open input Interval file")
31     
32 def stop_err(msg):
33     sys.stderr.write(msg)
34     sys.exit()
35 
36 def reverse_complement(text):
37     DNA_COMP = string.maketrans( "ACGTacgt", "TGCAtgca" )
38     comp = [ch for ch in text.translate(DNA_COMP)]
39     comp.reverse()
40     return "".join(comp)
41 
42 def get_unique_elems(elems):
43     seen=set()
44     return[x for x in elems if x not in seen and not seen.add(x)]
45 
46 def get_binned_lists(uniqlist, binsize):
47     binnedlist=[]
48     uniqlist.sort()
49     start = int(uniqlist[0])
50     bin_ind=0
51     l_ind=0
52     binnedlist.append([])
53     while l_ind < len(uniqlist):
54         elem = int(uniqlist[l_ind])
55         if elem in range(start,start+binsize):
56             binnedlist[bin_ind].append(elem)
57         else:
58             start += binsize
59             bin_ind += 1
60             binnedlist.append([])
61             binnedlist[bin_ind].append(elem)
62         l_ind += 1
63     return binnedlist
64 
65 def fetch_weight(H,C,t):
66     if (H-(C-H)) < t:
67         return 2.0
68     else:
69         return 1.0
70 
71 def mutabilityEstimator(repeats1,repeats2,thresholds):
72     mut_num = 0.0    #Mutability Numerator
73     mut_den = 0.0    #Mutability denominator
74     for ind,H in enumerate(repeats1):
75         C = repeats2[ind]
76         t = thresholds[ind]
77         w = fetch_weight(H,C,t)
78         mut_num += ((H-C)*(H-C)*w)
79         mut_den += w
80     return [mut_num, mut_den]
81 
82 def output_writer(blk, blk_lines):
83     global winspecies, speciesind
84     all_elems_1=[]
85     all_elems_2=[]
86     all_s_elems_1=[]
87     all_s_elems_2=[]
88     for bline in blk_lines:
89         if not(bline):
90             continue
91         items = bline.split('\t')
92         seq1 = items[1]
93         start1 = items[2]
94         end1 = items[3]
95         seq2 = items[8]
96         start2 = items[9]
97         end2 = items[10] 
98         if p_group_cols[0] == 6:
99             items[p_group_cols[0]] = int(items[p_group_cols[0]])
00             items[p_group_cols[1]] = int(items[p_group_cols[1]])
01         if s_group_cols[0] == 6:
02             items[s_group_cols[0]] = int(items[s_group_cols[0]])
03             items[s_group_cols[1]] = int(items[s_group_cols[1]])
04         all_elems_1.append(items[p_group_cols[0]])    #primary col elements for species 1
05         all_elems_2.append(items[p_group_cols[1]])    #primary col elements for species 2
06         if s_group_cols[0] != -1:    #sub-group is not None
07             all_s_elems_1.append(items[s_group_cols[0]])    #secondary col elements for species 1
08             all_s_elems_2.append(items[s_group_cols[1]])    #secondary col elements for species 2
09     uniq_elems_1 = get_unique_elems(all_elems_1)
10     uniq_elems_2 = get_unique_elems(all_elems_2)
11     if s_group_cols[0] != -1:
12         uniq_s_elems_1 = get_unique_elems(all_s_elems_1)
13         uniq_s_elems_2 = get_unique_elems(all_s_elems_2)
14     mut1={}
15     mut2={}
16     count1 = {}
17     count2 = {}
18     """
19     if p_group_cols[0] == 7:    #i.e. the option chosen is group-by unit(AG, GTC, etc)
20         uniq_elems_1 = get_unique_units(j.sort(lambda x, y: len(x)-len(y)))
21     """
22     if p_group_cols[0] == 6:    #i.e. the option chosen is group-by repeat number.
23         uniq_elems_1 = get_binned_lists(uniq_elems_1,p_bin_size)
24         uniq_elems_2 = get_binned_lists(uniq_elems_2,p_bin_size)
25         
26     if s_group_cols[0] == 6:    #i.e. the option chosen is subgroup-by repeat number.
27         uniq_s_elems_1 = get_binned_lists(uniq_s_elems_1,s_bin_size)
28         uniq_s_elems_2 = get_binned_lists(uniq_s_elems_2,s_bin_size)
29 
30     for pitem1 in uniq_elems_1:
31         #repeats1 = []
32         #repeats2 = []
33         thresholds = []
34         if s_group_cols[0] != -1:    #Sub-group by feature is not None
35             for sitem1 in uniq_s_elems_1:
36                 repeats1 = []
37                 repeats2 = []
38                 if type(sitem1) == type(''):
39                     sitem1 = sitem1.strip()
40                 for bline in blk_lines:
41                     belems = bline.split('\t')
42                     if type(pitem1) == list:
43                         if p_group_cols[0] == 6:
44                             belems[p_group_cols[0]] = int(belems[p_group_cols[0]])
45                         if belems[p_group_cols[0]] in pitem1:
46                             if belems[s_group_cols[0]]==sitem1:
47                                 repeats1.append(int(belems[6]))
48                                 repeats2.append(int(belems[13]))
49                                 if belems[4] == 'mononucleotide':
50                                     thresholds.append(mono_threshold)
51                                 else:
52                                     thresholds.append(non_mono_threshold)
53                                 mut1[str(pitem1)+'\t'+str(sitem1)]=mutabilityEstimator(repeats1,repeats2,thresholds)
54                                 if region == 'align':
55                                     count1[str(pitem1)+'\t'+str(sitem1)]=min(sum(repeats1),sum(repeats2))
56                                 else:    
57                                     if winspecies == 1:
58                                         count1["%s\t%s" %(pitem1,sitem1)]=sum(repeats1)
59                                     elif winspecies == 2:
60                                         count1["%s\t%s" %(pitem1,sitem1)]=sum(repeats2)
61                     else:
62                         if type(sitem1) == list:
63                             if s_group_cols[0] == 6:
64                                 belems[s_group_cols[0]] = int(belems[s_group_cols[0]])
65                             if belems[p_group_cols[0]]==pitem1 and belems[s_group_cols[0]] in sitem1:
66                                 repeats1.append(int(belems[6]))
67                                 repeats2.append(int(belems[13]))
68                                 if belems[4] == 'mononucleotide':
69                                     thresholds.append(mono_threshold)
70                                 else:
71                                     thresholds.append(non_mono_threshold)
72                                 mut1["%s\t%s" %(pitem1,sitem1)]=mutabilityEstimator(repeats1,repeats2,thresholds)
73                                 if region == 'align':
74                                     count1[str(pitem1)+'\t'+str(sitem1)]=min(sum(repeats1),sum(repeats2))
75                                 else:    
76                                     if winspecies == 1:
77                                         count1[str(pitem1)+'\t'+str(sitem1)]=sum(repeats1)
78                                     elif winspecies == 2:
79                                         count1[str(pitem1)+'\t'+str(sitem1)]=sum(repeats2)
80                         else:
81                             if belems[p_group_cols[0]]==pitem1 and belems[s_group_cols[0]]==sitem1:
82                                 repeats1.append(int(belems[6]))
83                                 repeats2.append(int(belems[13]))
84                                 if belems[4] == 'mononucleotide':
85                                     thresholds.append(mono_threshold)
86                                 else:
87                                     thresholds.append(non_mono_threshold)
88                                 mut1["%s\t%s" %(pitem1,sitem1)]=mutabilityEstimator(repeats1,repeats2,thresholds)
89                                 if region == 'align':
90                                     count1[str(pitem1)+'\t'+str(sitem1)]=min(sum(repeats1),sum(repeats2))
91                                 else:    
92                                     if winspecies == 1:
93                                         count1["%s\t%s" %(pitem1,sitem1)]=sum(repeats1)
94                                     elif winspecies == 2:
95                                         count1["%s\t%s" %(pitem1,sitem1)]=sum(repeats2)
96         else:   #Sub-group by feature is None
97             for bline in blk_lines:
98                 belems = bline.split('\t')
99                 if type(pitem1) == list:
00                     #print >>sys.stderr, "item: " + str(item1)
01                     if p_group_cols[0] == 6:
02                         belems[p_group_cols[0]] = int(belems[p_group_cols[0]])
03                     if belems[p_group_cols[0]] in pitem1:
04                         repeats1.append(int(belems[6]))
05                         repeats2.append(int(belems[13]))
06                         if belems[4] == 'mononucleotide':
07                             thresholds.append(mono_threshold)
08                         else:
09                             thresholds.append(non_mono_threshold)
10                 else:
11                     if belems[p_group_cols[0]]==pitem1:
12                         repeats1.append(int(belems[6]))
13                         repeats2.append(int(belems[13]))
14                         if belems[4] == 'mononucleotide':
15                             thresholds.append(mono_threshold)
16                         else:
17                             thresholds.append(non_mono_threshold)
18             mut1["%s" %(pitem1)]=mutabilityEstimator(repeats1,repeats2,thresholds)
19             if region == 'align':
20                 count1["%s" %(pitem1)]=min(sum(repeats1),sum(repeats2))
21             else:            
22                 if winspecies == 1:
23                     count1[str(pitem1)]=sum(repeats1)
24                 elif winspecies == 2:
25                     count1[str(pitem1)]=sum(repeats2)
26                 
27     for pitem2 in uniq_elems_2:
28         #repeats1 = []
29         #repeats2 = []
30         thresholds = []
31         if s_group_cols[0] != -1:    #Sub-group by feature is not None
32             for sitem2 in uniq_s_elems_2:
33                 repeats1 = []
34                 repeats2 = []
35                 if type(sitem2)==type(''):
36                     sitem2 = sitem2.strip()
37                 for bline in blk_lines:
38                     belems = bline.split('\t')
39                     if type(pitem2) == list:
40                         if p_group_cols[0] == 6:
41                             belems[p_group_cols[1]] = int(belems[p_group_cols[1]])
42                         if belems[p_group_cols[1]] in pitem2 and belems[s_group_cols[1]]==sitem2:
43                             repeats2.append(int(belems[13]))
44                             repeats1.append(int(belems[6]))
45                             if belems[4] == 'mononucleotide':
46                                 thresholds.append(mono_threshold)
47                             else:
48                                 thresholds.append(non_mono_threshold)
49                             mut2["%s\t%s" %(pitem2,sitem2)]=mutabilityEstimator(repeats2,repeats1,thresholds)
50                             #count2[str(pitem2)+'\t'+str(sitem2)]=len(repeats2)
51                             if region == 'align':
52                                 count2["%s\t%s" %(pitem2,sitem2)]=min(sum(repeats1),sum(repeats2))
53                             else: 
54                                 if winspecies == 1:
55                                     count2["%s\t%s" %(pitem2,sitem2)]=len(repeats2)
56                                 elif winspecies == 2:
57                                     count2["%s\t%s" %(pitem2,sitem2)]=len(repeats1)
58                     else:
59                         if type(sitem2) == list:
60                             if s_group_cols[0] == 6:
61                                 belems[s_group_cols[1]] = int(belems[s_group_cols[1]])
62                             if belems[p_group_cols[1]]==pitem2 and belems[s_group_cols[1]] in sitem2:
63                                 repeats2.append(int(belems[13]))
64                                 repeats1.append(int(belems[6]))
65                                 if belems[4] == 'mononucleotide':
66                                     thresholds.append(mono_threshold)
67                                 else:
68                                     thresholds.append(non_mono_threshold)
69                                 mut2["%s\t%s" %(pitem2,sitem2)]=mutabilityEstimator(repeats2,repeats1,thresholds)
70                                 if region == 'align':
71                                     count2["%s\t%s" %(pitem2,sitem2)]=min(sum(repeats1),sum(repeats2))
72                                 else: 
73                                     if winspecies == 1:
74                                         count2["%s\t%s" %(pitem2,sitem2)]=len(repeats2)
75                                     elif winspecies == 2:
76                                         count2["%s\t%s" %(pitem2,sitem2)]=len(repeats1)
77                         else:
78                             if belems[p_group_cols[1]]==pitem2 and belems[s_group_cols[1]]==sitem2:
79                                 repeats1.append(int(belems[13]))
80                                 repeats2.append(int(belems[6]))
81                                 if belems[4] == 'mononucleotide':
82                                     thresholds.append(mono_threshold)
83                                 else:
84                                     thresholds.append(non_mono_threshold)
85                                 mut2["%s\t%s" %(pitem2,sitem2)]=mutabilityEstimator(repeats2,repeats1,thresholds)
86                                 if region == 'align':
87                                     count2["%s\t%s" %(pitem2,sitem2)]=min(sum(repeats1),sum(repeats2))
88                                 else: 
89                                     if winspecies == 1:
90                                         count2["%s\t%s" %(pitem2,sitem2)]=len(repeats2)
91                                     elif winspecies == 2:
92                                         count2["%s\t%s" %(pitem2,sitem2)]=len(repeats1)
93         else:   #Sub-group by feature is None
94             for bline in blk_lines:
95                 belems = bline.split('\t')
96                 if type(pitem2) == list:
97                     if p_group_cols[0] == 6:
98                         belems[p_group_cols[1]] = int(belems[p_group_cols[1]])
99                     if belems[p_group_cols[1]] in pitem2:
00                         repeats2.append(int(belems[13]))
01                         repeats1.append(int(belems[6]))
02                         if belems[4] == 'mononucleotide':
03                             thresholds.append(mono_threshold)
04                         else:
05                             thresholds.append(non_mono_threshold)
06                 else:
07                     if belems[p_group_cols[1]]==pitem2:
08                         repeats2.append(int(belems[13]))
09                         repeats1.append(int(belems[6]))
10                         if belems[4] == 'mononucleotide':
11                             thresholds.append(mono_threshold)
12                         else:
13                             thresholds.append(non_mono_threshold)
14             mut2["%s" %(pitem2)]=mutabilityEstimator(repeats2,repeats1,thresholds)
15             if region == 'align':
16                 count2["%s" %(pitem2)]=min(sum(repeats1),sum(repeats2))
17             else:
18                 if winspecies == 1:
19                     count2["%s" %(pitem2)]=sum(repeats2)
20                 elif winspecies == 2:
21                     count2["%s" %(pitem2)]=sum(repeats1)
22     for key in mut1.keys():
23         if key in mut2.keys():
24             mut = (mut1[key][0]+mut2[key][0])/(mut1[key][1]+mut2[key][1])
25             count = count1[key]
26             del mut2[key]
27         else:
28             unit_found = False
29             if p_group_cols[0] == 7 or s_group_cols[0] == 7: #if it is Repeat Unit (AG, GCT etc.) check for reverse-complements too
30                 if p_group_cols[0] == 7:
31                     this,other = 0,1
32                 else:
33                     this,other = 1,0
34                 groups1 = key.split('\t')
35                 mutn = mut1[key][0]
36                 mutd = mut1[key][1]
37                 count = 0
38                 for key2 in mut2.keys():
39                     groups2 = key2.split('\t')
40                     if groups1[other] == groups2[other]:
41                         if groups1[this] in groups2[this]*2 or reverse_complement(groups1[this]) in groups2[this]*2:
42                             #mut = (mut1[key][0]+mut2[key2][0])/(mut1[key][1]+mut2[key2][1])
43                             mutn += mut2[key2][0]
44                             mutd += mut2[key2][1]
45                             count += int(count2[key2])
46                             unit_found = True
47                             del mut2[key2]
48                             #break
49             if unit_found:
50                 mut = mutn/mutd
51             else:
52                 mut = mut1[key][0]/mut1[key][1]
53                 count = count1[key]
54         mut = "%.2e" %(mut/num_generations)
55         if region == 'align':
56             print >>fout, str(blk) + '\t'+seq1 + '\t' + seq2 + '\t' +key.strip()+ '\t'+str(mut) + '\t'+ str(count)
57         elif region == 'win':
58             fout.write("%s\t%s\t%s\t%s\n" %(blk,key.strip(),mut,count))
59             fout.flush()
60             
61     #catch any remaining repeats, for instance if the orthologous position contained different repeat units
62     for remaining_key in mut2.keys():
63         mut = mut2[remaining_key][0]/mut2[remaining_key][1]
64         mut = "%.2e" %(mut/num_generations)
65         count = count2[remaining_key]
66         if region == 'align':
67             print >>fout, str(blk) + '\t'+seq1 + '\t'+seq2 + '\t'+remaining_key.strip()+ '\t'+str(mut)+ '\t'+ str(count)
68         elif region == 'win':
69             fout.write("%s\t%s\t%s\t%s\n" %(blk,remaining_key.strip(),mut,count))
70             fout.flush()
71             #print >>fout, blk + '\t'+remaining_key.strip()+ '\t'+str(mut)+ '\t'+ str(count)
72 
73 def counter(node, start, end, report_func):
74     if start <= node.start < end and start < node.end <= end:
75         report_func(node) 
76         if node.right:
77             counter(node.right, start, end, report_func)
78         if node.left:
79             counter(node.left, start, end, report_func)
80     elif node.start < start and node.right:
81         counter(node.right, start, end, report_func)
82     elif node.start >= end and node.left and node.left.maxend > start:
83         counter(node.left, start, end, report_func)
84             
85         
86 def main():
87     infile = sys.argv[1]
88     
89     for i, line in enumerate( file ( infile )):
90         line = line.rstrip('\r\n')
91         if len( line )>0 and not line.startswith( '#' ):
92             elems = line.split( '\t' )
93             break
94         if i == 30:
95             break # Hopefully we'll never get here...
96     
97     if len( elems ) != 15:
98         stop_err( "This tool only works on tabular data output by 'Extract Orthologous Microsatellites from pair-wise alignments' tool. The data in your input dataset is either missing or not formatted properly." )
99     global winspecies, speciesind
00     if region == 'win':
01         if dbkey_i in elems[1]:
02             winspecies = 1
03             speciesind = 1 
04         elif dbkey_i in elems[8]:
05             winspecies = 2
06             speciesind = 8
07         else:
08             stop_err("The species build corresponding to your interval file is not present in the Microsatellite file.") 
09         
10     fin = open(infile, 'r')
11     skipped = 0
12     blk=0
13     win=0
14     linestr=""
15     
16     if region == 'win':
17         
18         msats = NiceReaderWrapper( fileinput.FileInput( infile ),
19                                 chrom_col = speciesind,
20                                 start_col = speciesind+1,
21                                 end_col = speciesind+2,
22                                 strand_col = -1,
23                                 fix_strand = True)
24         msatTree = quicksect.IntervalTree()
25         for item in msats:
26             if type( item ) is GenomicInterval:
27                 msatTree.insert( item, msats.linenum, item.fields )
28         
29         for iline in fint:
30             try:
31                 iline = iline.rstrip('\r\n')
32                 if not(iline) or iline == "":
33                     continue
34                 ielems = iline.strip("\r\n").split('\t')
35                 ichr = ielems[chr_col_i]
36                 istart = int(ielems[start_col_i])
37                 iend = int(ielems[end_col_i])
38                 isrc = "%s.%s" %(dbkey_i,ichr)
39                 if isrc not in msatTree.chroms:
40                     continue
41                 result = []
42                 root = msatTree.chroms[isrc]    #root node for the chrom
43                 counter(root, istart, iend, lambda node: result.append( node ))
44                 if not(result):
45                     continue
46                 tmpfile1 = tempfile.NamedTemporaryFile('wb+')
47                 for node in result:
48                     tmpfile1.write("%s\n" % "\t".join( node.other ))
49                 
50                 tmpfile1.seek(0)
51                 output_writer(iline, tmpfile1.readlines())
52             except:
53                 skipped+=1
54         if skipped:
55             print "Skipped %d intervals as invalid." %(skipped)
56     elif region == 'align':
57         if s_group_cols[0] != -1:
58             print >>fout, "#Window\tSpecies_1\tSpecies_2\tGroupby_Feature\tSubGroupby_Feature\tMutability\tCount"
59         else:
60             print >>fout, "#Window\tSpecies_1\tWindow_Start\tWindow_End\tSpecies_2\tGroupby_Feature\tMutability\tCount"
61         prev_bnum = -1
62         try:
63             for line in fin:
64                 line = line.strip("\r\n")
65                 if not(line) or line == "":
66                     continue
67                 elems = line.split('\t')
68                 try:
69                     assert int(elems[0])
70                     assert len(elems) == 15
71                 except:
72                     continue
73                 new_bnum = int(elems[0])
74                 if new_bnum != prev_bnum:
75                     if prev_bnum != -1:
76                         output_writer(prev_bnum, linestr.strip().replace('\r','\n').split('\n'))
77                     linestr = line + "\n"
78                 else:
79                     linestr += line
80                     linestr += "\n"
81                 prev_bnum = new_bnum
82             output_writer(prev_bnum, linestr.strip().replace('\r','\n').split('\n'))
83         except Exception, ea:
84             print >>sys.stderr, ea
85             skipped += 1
86         if skipped:
87             print "Skipped %d lines as invalid." %(skipped)
88 if __name__ == "__main__":
89     main()