root/galaxy-central/tools/regVariation/getIndelRates_3way.py

#!/usr/bin/env python
#Guruprasad Ananda

from galaxy import eggs
import pkg_resources
pkg_resources.require( "bx-python" )

import sys, os, tempfile
import traceback
import fileinput
from warnings import warn

from galaxy.tools.util.galaxyops import *
from bx.intervals.io import *

from bx.intervals.operations import quicksect

def stop_err(msg):
    sys.stderr.write(msg)
    sys.exit()
    
def counter(node, start, end, sort_col):
    global full, blk_len, blk_list
    if node.start < start:
        if node.right:
            counter(node.right, start, end, sort_col)
    elif start <= node.start <= end and start <= node.end <= end:
        full += 1
        if node.other[0] not in blk_list:
            blk_list.append(node.other[0])
            blk_len += int(node.other[sort_col+2])
        if node.left and node.left.maxend > start:
            counter(node.left, start, end, sort_col)
        if node.right: 
            counter(node.right, start, end, sort_col)
    elif node.start > end:
        if node.left: 
            counter(node.left, start, end, sort_col)
            

infile = sys.argv[1]  
fout = open(sys.argv[2],'w')
int_file = sys.argv[3]
if int_file != "None": #User has specified an interval file
    try:
        fint = open(int_file, 'r')
        dbkey_i = sys.argv[4]
        chr_col_i, start_col_i, end_col_i, strand_col_i = parse_cols_arg( sys.argv[5] )
    except:
        stop_err("Unable to open input Interval file")
        
def main():

    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 ) != 18:
        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." )
    
    for i, line in enumerate( file ( infile )):
        line = line.rstrip('\r\n')
        elems = line.split('\t')
        try:
            assert int(elems[0])
            assert len(elems) == 18
            if int_file != "None":
                if dbkey_i not in elems[3] and  dbkey_i not in elems[8] and dbkey_i not in elems[13]:
                    stop_err("The species build corresponding to your interval file is not present in the Indel file.") 
                if dbkey_i in elems[3]:
                    sort_col = 4
                elif dbkey_i in elems[8]:
                    sort_col = 9
                elif dbkey_i in elems[13]:
                    sort_col = 14
            else:
                species = []
                species.append( elems[3].split('.')[0] )
                species.append( elems[8].split('.')[0] )
                species.append( elems[13].split('.')[0] )
                sort_col = 0    #Based on block numbers
            break
        except:
            continue
        
        
    fin = open(infile, 'r')
    skipped = 0
    
    if int_file == "None":
        sorted_infile = tempfile.NamedTemporaryFile()
        cmdline = "sort -n -k"+str(1)+" -o "+sorted_infile.name+" "+infile
        try:
            os.system(cmdline)
        except:
            stop_err("Encountered error while sorting the input file.")
        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])
        prev_bnum = -1
        sorted_infile.seek(0)
        for line in sorted_infile.readlines():
            line = line.rstrip('\r\n')
            elems = line.split('\t')
            try:
                assert int(elems[0])
                assert len(elems) == 18
                new_bnum = int(elems[0])
                if new_bnum != prev_bnum:
                    if prev_bnum != -1:
                        irate = []
                        drate = []
                        for i,elem in enumerate(inserts):
                            try:
                                irate.append(str("%.2e" %(inserts[i]/blen[i])))
                            except:
                                irate.append('0')
                            try:
                                drate.append(str("%.2e" %(deletes[i]/blen[i])))
                            except:
                                drate.append('0')
                        print >>fout, "%s\t%s\t%s" %(prev_bnum, '\t'.join(irate) , '\t'.join(drate))
                    inserts = [0.0, 0.0, 0.0]
                    deletes = [0.0, 0.0, 0.0]
                    blen = []
                    blen.append( int(elems[6]) )
                    blen.append( int(elems[11]) )
                    blen.append( int(elems[16]) )
                line_sp = elems[1].split('.')[0]
                sp_ind = species.index(line_sp)
                if elems[1].endswith('insert'):
                    inserts[sp_ind] += 1
                elif elems[1].endswith('delete'):
                    deletes[sp_ind] += 1
                prev_bnum = new_bnum 
            except Exception, ei:
                #print >>sys.stderr, ei
                continue
        irate = []
        drate = []
        for i,elem in enumerate(inserts):
            try:
                irate.append(str("%.2e" %(inserts[i]/blen[i])))
            except:
                irate.append('0')
            try:
                drate.append(str("%.2e" %(deletes[i]/blen[i])))
            except:
                drate.append('0')
        print >>fout, "%s\t%s\t%s" %(prev_bnum, '\t'.join(irate) , '\t'.join(drate))
        sys.exit()
    
    
    inf = open(infile, 'r')
    start_met = False
    end_met = False
    sp_file = tempfile.NamedTemporaryFile()
    for n, line in enumerate(inf):
        line = line.rstrip('\r\n')
        elems = line.split('\t')
        try:
            assert int(elems[0])
            assert len(elems) == 18
            if dbkey_i not in elems[1]: 
                if not(start_met):   
                    continue
                else:
                    sp_end = n
                    break
            else:
                print >>sp_file, line
                if not(start_met):
                    start_met = True
                    sp_start = n
        except:
            continue
    
    try:
        assert sp_end
    except:
        sp_end = n+1
    
    sp_file.seek(0)
    win = NiceReaderWrapper( fileinput.FileInput( int_file ),
                                chrom_col=chr_col_i,
                                start_col=start_col_i,
                                end_col=end_col_i,
                                strand_col=strand_col_i,
                                fix_strand=True)
    
    indel = NiceReaderWrapper( fileinput.FileInput( sp_file.name ),
                                chrom_col=1,
                                start_col=sort_col,
                                end_col=sort_col+1,
                                strand_col=-1,
                                fix_strand=True)
    
    indelTree = quicksect.IntervalTree()
    for item in indel:
        if type( item ) is GenomicInterval:
            indelTree.insert( item, indel.linenum, item.fields )
    result=[]
    
    global full, blk_len, blk_list
    for interval in win:
        if type( interval ) is Header:
            pass
        if type( interval ) is Comment:
            pass
        elif type( interval ) == GenomicInterval:
            chrom = interval.chrom
            start = int(interval.start)
            end = int(interval.end)
            if start > end: 
                warn( "Interval start after end!" )
            ins_chr = "%s.%s_insert" %(dbkey_i,chrom)
            del_chr = "%s.%s_delete" %(dbkey_i,chrom)
            irate = 0
            drate = 0
            if ins_chr not in indelTree.chroms and del_chr not in indelTree.chroms:
                pass    
            else:
                if ins_chr in indelTree.chroms:
                    full = 0.0
                    blk_len = 0
                    blk_list = []
                    root = indelTree.chroms[ins_chr]    #root node for the chrom insertion tree
                    counter(root, start, end, sort_col)
                    if blk_len:
                        irate = full/blk_len
                
                if del_chr in indelTree.chroms:
                    full = 0.0
                    blk_len = 0
                    blk_list = []
                    root = indelTree.chroms[del_chr]    #root node for the chrom insertion tree
                    counter(root, start, end, sort_col)
                    if blk_len:
                        drate = full/blk_len
                
            interval.fields.append(str("%.2e" %irate))
            interval.fields.append(str("%.2e" %drate))
            print >>fout, "\t".join(interval.fields)
            fout.flush()

if __name__ == "__main__":
    main()