thippo-Enhancer_Finding/main.py at master · bioinformatics-cnu/thippo-Enhancer_Finding · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
#参数设置

variable_list=['scaffold1','scaffold2','scaffold3','scaffold5']
evalue='20'
species_name='Bacillus.sp.T61'
query_len=80

import re

def find_gene_downstream(position_tuple):
    if '-' in position_tuple:
        n='inside'
        for i in open(variable+'.gff','r').readlines():
            element_list=i.split('	')
            if ((int(element_list[1].strip())-int(position_tuple[1])<0 and int(element_list[2].strip())-int(position_tuple[1])>0) or (int(element_list[1].strip())-int(position_tuple[2])<0 and int(element_list[2].strip())-int(position_tuple[2])>0)):
                return 'inside'
            else:
                if (int(element_list[2].strip())-int(position_tuple[1]))>0:
                    return n
                else:
                    n=element_list[0]
    else:
        for i in open(variable+'.gff','r').readlines():
            element_list=i.split('	')
            if ((int(element_list[1].strip())-int(position_tuple[1])<0 and int(element_list[2].strip())-int(position_tuple[1])>0) or (int(element_list[1].strip())-int(position_tuple[2])<0 and int(element_list[2].strip())-int(position_tuple[2])>0)):
                return 'inside'
            else:
                if (int(element_list[1].strip())-int(position_tuple[2]))>0:
                    return element_list[0]

def find_gene_upstream(position_tuple):
    if '-' in position_tuple:
        n='inside'
        for i in open(variable+'.gff','r').readlines():
            element_list=i.split('	')
            if ((int(element_list[1].strip())-int(position_tuple[1])<0 and int(element_list[2].strip())-int(position_tuple[1])>0) or (int(element_list[1].strip())-int(position_tuple[2])<0 and int(element_list[2].strip())-int(position_tuple[2])>0)):
                return 'inside'
            else:
                if (int(element_list[1].strip())-int(position_tuple[2]))>0:
                    return element_list[0]
    else:
        for i in open(variable+'.gff','r').readlines():
            element_list=i.split('	')
            if ((int(element_list[1].strip())-int(position_tuple[1])<0 and int(element_list[2].strip())-int(position_tuple[1])>0) or (int(element_list[1].strip())-int(position_tuple[2])<0 and int(element_list[2].strip())-int(position_tuple[2])>0)):
                return 'inside'
            else:
                if (int(element_list[2].strip())-int(position_tuple[1]))>0:
                    return n
                else:
                    n=element_list[0]

annotation_dict={}
count_dict={}

for i in open('Bacillus.sp.T61.annotation.table','r').readlines():
    lin=re.findall('(Bacillus.sp.T61GL[0-9]{6})	({.*})',i)
    annotation_dict[lin[0][0]]=lin[0][1]

for variable in variable_list:

    out_down=open('result_e'+evalue+'_'+variable+'_down.th','w')
    out_up=open('result_e'+evalue+'_'+variable+'_up.th','w')

    f=open('e'+evalue+'_'+variable+'.th','r').read()

    range_list=f.split('Range ')
    position_match=[]

    for i in range_list:
        if i:
            n=0
            if 'Minus' in i:
                position_min_list=re.findall(': ([0-9]*?) to',i)
                position_max_list=re.findall(' to ([0-9]*?)GraphicsNext',i)
                number=re.findall('^([0-9]*?): ',i)
                position_match_tuple=re.findall('Query  ([0-9]*?)[^0-9]*?([0-9]*?)\n',i)
                gene_name_down=find_gene_downstream(('-',position_min_list[0],position_max_list[0]))
                gene_name_up=find_gene_upstream(('-',position_min_list[0],position_max_list[0]))
                if 'inside' not in gene_name_down:
                    out_down.write(variable+'	'+'Range '+number[0]+'	'+gene_name_down+'	'+annotation_dict[gene_name_down]+"\n")
                    n+=1
                if 'inside' not in gene_name_up:
                    out_up.write(variable+'	'+'Range '+number[0]+'	'+gene_name_up+'	'+annotation_dict[gene_name_up]+"\n")
                    n+=1
            else:
                position_min_list=re.findall(': ([0-9]*?) to',i)
                position_max_list=re.findall(' to ([0-9]*?)GraphicsNext',i)
                number=re.findall('^([0-9]*?): ',i)
                position_match_tuple=re.findall('Query  ([0-9]*?)[^0-9]*?([0-9]*?)\n',i)
                gene_name_down=find_gene_downstream(('+',position_min_list[0],position_max_list[0]))
                gene_name_up=find_gene_upstream(('+',position_min_list[0],position_max_list[0]))
                if 'inside' not in gene_name_down:
                    out_down.write(variable+'	'+'Range '+number[0]+'	'+gene_name_down+'	'+annotation_dict[gene_name_down]+"\n")
                    n+=1
                if 'inside' not in gene_name_up:
                    out_up.write(variable+'	'+'Range '+number[0]+'	'+gene_name_up+'	'+annotation_dict[gene_name_up]+"\n")
                    n+=1
            if n!=0:
                for p in position_match_tuple:
                    position_match.append(p)

    out_down.close
    out_up.close

    for i in position_match:
        for n in [str(x) for x in range(int(i[0]),int(i[1])+1)]:
            if n in count_dict.keys():
                count_dict[n]+=1
            else:
                count_dict[n]=1

out_R=open('count_plot.R','w')
out_R.write('x=c(')
for k in [x for x in range(0,query_len)]:
    if k==query_len-1:
        out_R.write(str(k+1)+')\n')
    else:
        out_R.write(str(k+1)+',')
out_R.write('y=c(')
for l in [x for x in range(0,query_len)]:
    if l==query_len-1:
        if str(l) in count_dict.keys():
            out_R.write(str(count_dict[str(l)])+')\n')
        else:
            out_R.write('0'+')\n')
    else:
        if str(l) in count_dict.keys():
            out_R.write(str(count_dict[str(l)])+',')
        else:
            out_R.write('0'+',')
out_R.write('plot(x,y,type="h",main="abundance: '+species_name+'",xlab="query_seq",ylab="hit_number")')
out_R.close