annotate reactivity_cal/react_norm_function.py @ 5:7a8ddf1819b1 draft

Uploaded

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import sys
from Bio import SeqIO
import math
from parse_dis_react import *

def cap(a,value):
    if a>=value:
        return value
    else:
        return a

def react_norm(react_file, result_file, capped_value):
    print("Normalizing.....")
    react1 = parse_dist(react_file)
    react = react1[1]
    h = file(result_file, 'w')

    capped = int(capped_value)

    all_react = []


    for t in react:
        if react[t]!='null':
            for i in range(len(react[t])):
                if react[t][i]!='NA':                   
                    all_react.append(float(react[t][i]))
#                    except:
#                        print(react[t][i])
#                        print(t)
#                        print(i)

    all_react.sort(reverse = True)
    #print((all_react))
    #print(all_react[int(len(all_react)*0.02)])
    #print(all_react[int(len(all_react)*0.03)])
    #print(all_react[int(len(all_react)*0.025)])
    #print(all_react[int(len(all_react)*0.04)])
    #print(all_react[int(len(all_react)*0.05)])
    '''
    mean = sum(all_react)/len(all_react)
    print(mean)
    temp = 0

    for i in range(len(all_react)):
        temp = temp+all_react[i]*all_react[i]
    temp = temp/len(all_react)
    sd = math.sqrt(temp-mean*mean)
    '''
    eight = all_react[int(len(all_react)*0.02):int(len(all_react)*0.1)]
    meight = sum(eight)/len(eight)

    for t in react:
        h.write(t)
        h.write('\n')
        if react[t]!='null':
            if (t.find('AT1G29930')==-1) and (t.find('At1g29930')==-1):
                for i in range((len(react[t])-1)):
                    if react[t][i]!='NA':
                        h.write(str(cap((float(react[t][i])/meight),capped)))
                    else:
                        h.write('NA')
                    h.write('\t')
                if react[t][i+1]!='NA':
                    h.write(str(cap((float(react[t][i+1])/meight),capped)))
                else:
                    h.write('NA')
                h.write('\n')
            else:
                for i in range((len(react[t])-1)):
                    if react[t][i]!='NA':
                        h.write(str(float(react[t][i])*2.6))
                    else:
                        h.write('NA')
                    h.write('\t')
                if react[t][i+1]!='NA':
                    h.write(str(float(react[t][i])*2.6))
                else:
                    h.write('NA')
                h.write('\n')
                
                

    h.close()