석사일기_210927_#2_코드정리

와.............

너무 복잡하다.........

코드가...........

 

정리좀 해야겠다...........

 

22년2월22일 졸업 할 수 있겠지?

 

힘내서 코드 정리하자!

#!/usr/bin/env python
# coding: utf-8
 
# In[1]:
 
 
get_ipython().system('ls')
 
 
# In[2]:
 
 
get_ipython().system('cd')
 
 
# In[3]:
 
 
get_ipython().system('pwd')
 
 
# In[5]:
 
 
import os
path_dir = '/home/bis/KSG/actg/0419/FASTA'
file_list = os.listdir(path_dir)
 
 
# In[16]:
 
 
chromo_str_list= ['_' for i in range(25)]
 
 
# In[17]:
 
 
len(chromo_str_list)
 
 
# In[13]:
 
 
file_list
 
 
# In[15]:
 
 
a = ['1','2','3','4']
b = ''.join(a)
print(b)
 
 
# In[18]:
 
 
def make_one_line(ori):
    chromo_list = []
    for i in ori[1:]:
        chromo_list.append(i.replace("\n",""))
    line = ''.join(chromo_list)
    return line
 
 
# In[20]:
 
 
for i in file_list:
    print(i+" read")
    read_file = open(path_dir+"/"+i)
    chromo_str_ori = read_file.readlines()
    chromo = i.split(".")[0]
    chromo_line = make_one_line(chromo_str_ori)
    if chromo == 'X':
        chromo_str_list[23] = chromo_line
    elif chromo == 'Y':
        chromo_str_list[24] = chromo_line
    elif chromo == 'MT':
        chromo_str_list[0] = chromo_line
    else:
        chromo_str_list[int(chromo)] = chromo_line
 
print(len(chromo_str_list[0]))
print(chromo_str_list[0])
 
 
# In[23]:
 
 
print(chromo_str_list[0][0:5])
 
 
# In[27]:
 
 
print(chromo_str_list[3][52434459:52434462])
 
 
# In[28]:
 
 
get_ipython().system('ls')
 
 
# In[30]:
 
 
utr_filtered = open("utr_filter_0.bed")
utr_filtered_ori = utr_filtered.readlines()
utr_gene_dict = {}
 
for i in utr_filtered_ori[:]:
    gene = i.split("\t")[3].split("_")[2]
#     three_five = i.split("\t")[3].split("_")[3].split("(")[0]
    if utr_gene_dict.get(gene) == None:
        utr_gene_dict[gene] = i
    else:
        utr_gene_dict[gene] = utr_gene_dict.get(gene)+i
        print(utr_gene_dict.get(gene))
        
print(len(utr_gene_dict))
 
 
# In[44]:
 
 
gtf = open('/home/bis/KSG/actg/0419/GTF_whole_chr/Homo_sapiens_only_chro.GRCh37.75.gtf')
gtf_ori = gtf.readlines()
gtf_utr_dict = {}
 
 
# In[45]:
 
 
count = 0
for i in gtf_ori[5:]:
    gene = i.split("\t")[8].split("\"")[1]
    p_type = i.split("\t")[1]
    if utr_gene_dict.get(gene) != None and p_type == 'protein_coding':
        if gtf_utr_dict.get(gene) == None:
            gtf_utr_dict[gene] = i
            count+=1
        else:
            gtf_utr_dict[gene] = gtf_utr_dict.get(gene)+i
            count+=1
 
print(count)
 
 
# In[201]:
 
 
for i in gtf_transcript_utr_dict.keys():
    print('--------')
    print(i)
#     print(gtf_transcript_utr_dict.get(i))
    print('--------')
 
 
# In[203]:
 
 
for i in gtf_gene_transcript_utr_dict.keys():
    print('--------')
    print(i)
    print(gtf_gene_transcript_utr_dict.get(i))
    print('--------')
 
 
# In[202]:
 
 
gtf_transcript_utr_dict = {}
gtf_gene_transcript_utr_dict = {}
for i in gtf_ori[5:]:
    gene = i.split("\t")[8].split("\"")[1]
    p_type = i.split("\t")[1]
    if utr_gene_dict.get(gene) != None and p_type == 'protein_coding' and '; transcript_id ' in i :
        trans = i.split("\t")[8].split("\"")[3]
        if gtf_gene_transcript_utr_dict.get(gene) == None:
            gtf_gene_transcript_utr_dict[gene] = trans+"\n"
        else:
            if gtf_transcript_utr_dict.get(trans) == None:
                gtf_gene_transcript_utr_dict[gene] = gtf_gene_transcript_utr_dict.get(gene)+trans+"\n"
        if gtf_transcript_utr_dict.get(trans) == None:
            gtf_transcript_utr_dict[trans] = i
        else:
            gtf_transcript_utr_dict[trans] = gtf_transcript_utr_dict.get(trans)+i
            
print(len(gtf_transcript_utr_dict))
print(len(gtf_gene_transcript_utr_dict))
 
 
# In[73]:
 
 
def get_p_utr(UTRs_info,start,stop):
    Utr_3 = []
    Utr_5 = []
    for i in UTRs_info[1:]:
        start_utr = int(i.split("\t")[0])
        end_utr = int(i.split("\t")[1])
        if start_utr < start:
            Utr_5.append(str(start_utr)+":"+str(end_utr))
        elif start_utr >stop:
            Utr_3.append(str(start_utr)+":"+str(end_utr))
    return Utr_5, Utr_3
        
        
def get_m_utr(UTRs_info,start,stop):
    Utr_3 = []
    Utr_5 = []
    for i in UTRs_info[1:]:
        start_utr = int(i.split("\t")[0])
        end_utr = int(i.split("\t")[1])
        if start_utr > start:
            Utr_5.append(str(start_utr)+":"+str(end_utr))
        elif start_utr <stop:
            Utr_3.append(str(start_utr)+":"+str(end_utr))
    return Utr_5, Utr_3
 
 
# In[77]:
 
 
def get_p_start_end(t_info):
    CDSs = t_info.split("\tCDS\t")
    start = int(CDSs[1].split("\t")[0])
    end = int(CDSs[1].split("\t")[1])
    for i in CDSs[1:]:
        temp_start = int(i.split("\t")[0])
        temp_end = int(i.split("\t")[1])
        if temp_start < start:
            start = temp_start
        if temp_end > end:
            end = temp_end
    return start, end
 
def get_m_start_end(t_info):
    CDSs = t_info.split("\tCDS\t")
    start = int(CDSs[1].split("\t")[1])
    end = int(CDSs[1].split("\t")[0])
    for i in CDSs[1:]:
        temp_start = int(i.split("\t")[1])
        temp_end = int(i.split("\t")[0])
        if temp_start > start:
            start = temp_start
        if temp_end < end:
            end = temp_end
    return start, end
 
 
# In[81]:
 
 
a  = '\tstart_codon\t\tstop_codon\t\tUTR\t'
if '\tstart_codon\t'and'\tstop_codon\t'and'\tUTR\t' in trans_info:
    print('please')
 
 
# In[83]:
 
 
print(gtf_transcript_utr_dict.get('ENST00000472925'))
 
 
# In[95]:
 
 
def get_p_stop_codon(t_info):
    stop_codon_start = int(t_info.split("\tstop_codon\t")[1].split("\t")[0])
    stop_codon_end = int(t_info.split("\tstop_codon\t")[1].split("\t")[1])
    chromosome = t_info.split("\t")[0]
    chromo = 333
    if chromosome == 'MT':
        chromo = 0
    elif chromosome == 'X':
        chromo = 23
    elif chromosome == 'Y':
        chromo = 24
    else:
        chromo = int(chromosome)
    stop_codon = chromo_str_list[chromo][stop_codon_start-1:stop_codon_end]
    return stop_codon
 
def get_m_stop_codon(t_info):
    stop_codon_start = int(t_info.split("\tstop_codon\t")[1].split("\t")[1])
    stop_codon_end = int(t_info.split("\tstop_codon\t")[1].split("\t")[0])
    chromosome = t_info.split("\t")[0]
    chromo = 333
    if chromosome == 'MT':
        chromo = 0
    elif chromosome == 'X':
        chromo = 23
    elif chromosome == 'Y':
        chromo = 24
    else:
        chromo = int(chromosome)
    stop_codon = chromo_str_list[chromo][stop_codon_end-1:stop_codon_start]
    return stop_codon
 
 
# In[ ]:
 
 
for 
 
 
# In[107]:
 
 
 
 
 
# In[108]:
 
 
get_utr_info_from_transcript('ENST00000602662')
get_utr_info_from_transcript('ENST00000361705')
get_utr_info_from_transcript('ENST00000402904')
 
 
# In[114]:
 
 
count_up = 0
count_down = 0
for i in gtf_transcript_utr_dict.keys():
    p_m, utr_5, utr_3, stop_codon = get_utr_info_from_transcript(i)
    if p_m == '+':
        start = 0
        check = True
        if len(utr_5) > 1:
            for j in utr_5:
                now = int(j.split(":")[0])
                if now > start:
                    start = now
                else:
                    if check:
                        count_down +=1
                        check = False
        start = 0
        check = True
        if len(utr_3) > 1:
            for j in utr_3:
                now = int(j.split(":")[0])
                if now > start:
                    start = now
                else:
                    if check:
                        count_down +=1
                        check = False
    elif p_m == '-':
        start = 2**31
        check = True
        if len(utr_5) > 1:
            for j in utr_5:
                now = int(j.split(":")[0])
                if now < start:
                    start = now
                else:
                    if check:
                        count_up +=1
                        print(i)
                        print('utr-5')
                        print(utr_5)
                        check = False
        start = 2**31
        check = True
        if len(utr_3) > 1:
            for j in utr_3:
                now = int(j.split(":")[0])
                if now < start:
                    start = now
                else:
                    if check:
                        print(i)
                        print('utr-3')
                        print(utr_3)
                        count_up +=1
                        check = False
 
print(count_up)
print(count_down)
 
 
# In[206]:
 
 
count = 0
filtered_gtf_transcript_utr_dict = {}
def get_utr_info_from_transcript(i):
#     print(i)
    trans_info = gtf_transcript_utr_dict.get(i)
    p_m = trans_info.split("\t")[6]
    utr_3 = []
    utr_5 = []
    stop_codon = ''
    if '\tstart_codon\t'in trans_info and'\tstop_codon\t' in trans_info and'\tUTR\t' in trans_info:
        if p_m == '+':
            start = int(trans_info.split("\tstart_codon\t")[1].split("\t")[0])
            stop = int(trans_info.split("\tstop_codon\t")[1].split("\t")[1])
            UTRs = trans_info.split("\tUTR\t")
            utr_5, utr_3 = get_p_utr(UTRs, start, stop)
            stop_codon = get_p_stop_codon(trans_info)
        elif p_m == '-':
            start = int(trans_info.split("\tstart_codon\t")[1].split("\t")[1])
            stop = int(trans_info.split("\tstop_codon\t")[1].split("\t")[0])
            UTRs = trans_info.split("\tUTR\t")
            utr_5, utr_3 = get_m_utr(UTRs, start, stop)
            stop_codon = get_m_stop_codon(trans_info)
#     print('p_m')
#     print(p_m)
#     print('utr_5')
#     print(utr_5)
#     print('utr_3')
#     print(utr_3)
#     print('stop_codon')
#     print(stop_codon)
#     print("----")
#     print(gtf_transcript_utr_dict.get(i))
    return p_m, utr_5, utr_3, stop_codon
 
 
# In[ ]:
 
 
TGA TAA TAG
TCA TTA CTA
 
 
# In[ ]:
 
 
gtf_transcript_utr_dict
filtered_gtf_transcript_utr_dict = {}
 
 
# In[66]:
 
 
##gtf의 시작-1:마지막
print(chromo_str_list[1][2238201:2238204])
print(chromo_str_list[1][3774167:3774170])
 
 
# In[ ]:
 
 
for i in gtf_gene_transcript_utr_dict.keys():
    transcripts = gtf_gene_transcript_utr_dict.get(i).split("\n")
    p_m = gtf_utr_dict.get(i).split("\n")[0].split("\t")[6]
    utr_pep_candidates = double_check(utr_gene_dict.get(i),p_m)
    for j in utr_pep_candidates:
        
 
 
# In[ ]:
 
 
 
 
 
# In[117]:
 
 
##utr 범위 안에 지금 캔디가 있는지 확인
def in_utr_3(utr_info,last_utr_int,p_m):
    safe = False
    check_list = []
    for i in utr_info:
        start = int(i.split(':')[0])
        end = int(i.split(':')[1])
        if start <= last_utr_int and end >=last_utr_int:
            safe = True
            if p_m == '+':
                list_input = str(start)+':'+str(last_utr_int-1)
                check_list.append(list_input)
            else:
                list_input = str(last_utr_int+1)+':'+str(end)
                check_list.append(list_input)
        elif safe == False:
            check_list.append(i)
    return safe, check_list
 
def in_utr_5(utr_info,last_utr_int,p_m):
    safe = False
    check_list = []
    for i in utr_info:
        start = int(i.split(':')[0])
        end = int(i.split(':')[1])
        if start <= last_utr_int and end >=last_utr_int and safe == False:
            safe = True
            if p_m == '+':
                list_input = str(last_utr_int+1)+':'+str(end)
                check_list.append(list_input)
            else:
                list_input = str(start)+':'+str(last_utr_int-1)
                check_list.append(list_input)
        elif safe == True:
            check_list.append(i)
    return safe, check_list
 
 
# In[124]:
 
 
a = ['1','2','3']
for i in range(len(a)-1,-1,-1):
#     print(i)
    print(a[i])
 
 
# In[125]:
 
 
def get_utr_seq(candi_seq_position, transcript,p_m):
    chromosome = gtf_transcript_utr_dict.get(transcript).split("\t")[0]
    chromo = 333
    if chromosome == 'MT':
        chromo = 0
    elif chromosome == 'X':
        chromo = 23
    elif chromosome == 'Y':
        chromo = 24
    else:
        chromo = int(chromosome)
    seq = []
    if p_m == '+':
        for i in candi_seq_position:
            start = int(i.split(":")[0])
            end = int(i.split(":")[1])
            seq.append(chromo_str_list[chromo][start-1:end])
        line = ''.join(seq)
    elif p_m == '-':
        for i in range(len(candi_seq_position)-1,-1,-1):
            start = int(candi_seq_position[i].split(":")[0])
            end = int(candi_seq_position[i].split(":")[1])
            seq.append(chromo_str_list[chromo][start-1:end])
        line = ''.join(seq)
    return line
 
 
# In[187]:
 
 
##3쪽이든 5쪽이든 utr 시작부터 double_check_output까지 다 더하고 3의 배수인지 확인하면서 자른 리스트 안에 stop코돈있는지 확인
##3의 배수이면서 stop 코돈 없으면 true 반환
def concat_utr(transcript, double_check_output):
    p_m, utr_5, utr_3, stop_codon = get_utr_info_from_transcript(transcript)
    three_five = double_check_output.split("_")[0]
    last_utr = int(double_check_output.split("_")[1])
    utr_seq = ''
    check = False
    if three_five == '3':
        check, candi_seq_position = in_utr_3(utr_3, last_utr,p_m)
        if check == False:
            return False
        else:
            if len(utr_3) == 0:
                return False
            else:
                utr_seq = get_utr_seq(candi_seq_position,transcript,p_m)
                if len(utr_seq)%3 != 0:
                    return False
                else:
                    check = check_stop_3(utr_seq,p_m,stop_codon)
    elif three_five == '5':
        check,candi_seq_position = in_utr_5(utr_5, last_utr,p_m)
        if check == False:
            return False
        else:
            if len(utr_5) == 0:
                return False
            else:
                utr_seq = get_utr_seq(candi_seq_position,transcript,p_m)
                if len(utr_seq)%3 != 0:
                    return False
                else:
                    check = check_stop_5(utr_seq, p_m)
    return check
 
 
# In[140]:
 
 
real_utr_bed_candi = {}
 
 
# In[139]:
 
 
a = '123'
print(a[0])
 
 
# In[170]:
 
 
print(utr_gene_dict.get('ENSG00000157933').split("\t")[3].split("_")[1])
 
 
# In[207]:
 
 
count = 0
real_utr_3 = 0
real_utr_5 = 0
fake_utr_3 = 0
fake_utr_5 = 0
 
utr_candi_dict = {}
for i in gtf_gene_transcript_utr_dict.keys():
    transcripts = gtf_gene_transcript_utr_dict.get(i).split("\n")
#     print('--------')
#     print(transcripts)
#     print('--------')
    p_m = gtf_utr_dict.get(i).split("\n")[0].split("\t")[6]
#     print(i)
    utr_pep_candidates, bed_lines = double_check(utr_gene_dict.get(i),p_m)
    count += len(utr_pep_candidates)
    check = False
    for j in range(len(utr_pep_candidates)):
        for transcript in transcripts[:-1]:
            check = concat_utr(transcript, utr_pep_candidates[j])
            if check == True:   
                three_five = utr_pep_candidates[j][0]
                utr_candi_dict[bed_lines[j]] = 1
                if three_five == '3':
                    real_utr_3 +=1
                elif three_five == '5':
                    real_utr_5 +=1
                break
        if check == False:
            three_five = utr_pep_candidates[j][0]
            if three_five == '3':
                fake_utr_3 +=1
            elif three_five == '5':
                fake_utr_5 +=1
    
print(count)
print(real_utr_3+real_utr_5+fake_utr_3+fake_utr_5)
print("r3")
print(real_utr_3)
print("r5")
print(real_utr_5)
print("f3")
print(fake_utr_3)
print("f5")
print(fake_utr_5)
print(len(utr_candi_dict))
print(real_utr_3+real_utr_5)
 
 
# In[127]:
 
 
plus_stop_codon = ['TGA', 'TAA' ,'TAG']
minus_stop_codon = ['TCA', 'TTA', 'CTA']
 
 
# In[132]:
 
 
a = '123456789abcdef'
print(a[3:6])
 
 
# In[135]:
 
 
if '123' in a:
    print('true')
 
 
# In[134]:
 
 
for i in range(15//3):
    b = a[int(i*3):int(i*3+3)]
    print(b)
 
 
# In[136]:
 
 
def check_stop_3(utr_seq,p_m,stop_codon):
    seq = []
    temp_plus_stop_codon = plus_stop_codon[:]
    temp_minus_stop_codon = minus_stop_codon[:]
    if p_m == '+':
        if stop_codon in temp_plus_stop_codon:
            temp_plus_stop_codon.remove(stop_codon)
        else:
            return False
        for i in range(len(utr_seq)//3):
            seq.append(utr_seq[int(i*3):int(i*3+3)])
        for i in temp_plus_stop_codon:
            if i in seq:
                return False
    elif p_m == '-':
        if stop_codon in temp_minus_stop_codon:
            temp_minus_stop_codon.remove(stop_codon)
        else:
            return False
        for i in range(len(utr_seq)//3):
            seq.append(utr_seq[int(i*3):int(i*3+3)])
        for i in temp_minus_stop_codon:
            if i in seq:
                return False
    return True
            
def check_stop_5(utr_seq,p_m):
    seq = []
    if p_m == '+':
        for i in range(len(utr_seq)//3):
            seq.append(utr_seq[int(i*3):int(i*3+3)])
        for i in plus_stop_codon:
            if i in seq:
                return False
    elif p_m == '-':
        for i in range(len(utr_seq)//3):
            seq.append(utr_seq[int(i*3):int(i*3+3)])
        for i in minus_stop_codon:
            if i in seq:
                return False
    return True
 
 
# In[60]:
 
 
for i in gtf_utr_dict.keys():
    if gtf_gene_transcript_utr_dict.get(i) == None:
        print*
 
 
# In[61]:
 
 
gtf_transcript_utr_dict.get('ENSG00000261915')
 
 
# In[47]:
 
 
gtf = open('Homo_sapiens_only_chro.GRCh37.75_for_utr_protein_coding_only.gtf','wt')
for i in gtf_utr_dict.keys():
    gtf.write(gtf_utr_dict.get(i))
gtf.close()
 
 
# In[39]:
 
 
len(gtf_utr_dict)
 
 
# In[182]:
 
 
a = 'asdf\nasdf\n'
print(len(a.split("\n")))
bb = a.split("\n")
for i in bb[:]:
    print("fuck!")
    print(i)
 
 
# In[197]:
 
 
## UTR pep 시작점
def double_check(bed, plus_minus):
    utr_num_dict = {}
    utr_num_line_dict = {}
    utr_final_list = []
    utr_bed_final_list = []
    lines = bed.split("\n")
#     print(lines)
    for k in lines[:-1]:
#         print(k)
        num = k.split("\t")[3].split("_")[1]
        three_five = k.split("\t")[3].split("_")[3].split("`")[0]
        if three_five == '3':
            plus_minus_3(plus_minus, k,num,utr_num_dict,utr_num_line_dict)
        elif three_five == '5':
            plus_minus_5(plus_minus, k,num,utr_num_dict,utr_num_line_dict)
    for i in lines[:-1]:
        num = i.split("\t")[3].split("_")[1]
        three_five = i.split("\t")[3].split("_")[3].split("`")[0]
        if three_five == '3':
            put = '3_'+str(utr_num_dict.get(num))
            utr_final_list.append(put)
            utr_bed_final_list.append(utr_num_line_dict.get(num))
        elif three_five == '5':
            put = '5_'+str(utr_num_dict.get(num))
            utr_final_list.append(put)
            utr_bed_final_list.append(utr_num_line_dict.get(num))
    return utr_final_list, utr_bed_final_list
 
 
# In[143]:
 
 
def plus_minus_3(plus_minus, line, num, n_dict, info_dict):
    start = 0
    if plus_minus == '+':
        start = int(line.split("\t")[1])+1
        if n_dict.get(num) == None:
            n_dict[num] = start
            info_dict[num] = line
        else:
            if n_dict.get(num) > start:
                n_dict[num] = start
                info_dict[num] = line
    elif plus_minus == '-':
        start = int(line.split("\t")[2])
        if n_dict.get(num) == None:
            n_dict[num] = start
            info_dict[num] = line
        else:
            if n_dict.get(num) < start:
                n_dict[num] = start
                info_dict[num] = line
 
 
# In[144]:
 
 
def plus_minus_5(plus_minus, line,num,n_dict,info_dict):
    end = 0
    if plus_minus == '+':
        end = int(line.split("\t")[2])
        if n_dict.get(num) == None:
            n_dict[num] = end
            info_dict[num] = line
        else:
            if n_dict.get(num) < end:
                n_dict[num] = end
                info_dict[num] = line
    elif plus_minus == '-':
        end = int(line.split("\t")[1])+1
        if n_dict.get(num) == None:
            n_dict[num] = end
            info_dict[num] = line
        else:
            if n_dict.get(num) > end:
                n_dict[num] = end
                info_dict[num] = line
 
 
# In[210]:
 
 
real_utr = open('real_utr_0927.bed','wt')
for i in utr_candi_dict.keys():
    real_utr.write(i+"\n")
real_utr.close()