기본

파이썬에는 두 가지 반복문이 있는데 한놈은 For고 한놈은 While이다. While은 베이직에서 본 거 같은디... (do while loop) 아무튼... 둘 다 반복문이긴 한데, 둘이 맥락은 좀 다르다. 내가 10페이지의 책을 읽는 것을 for와 while을 이용해 설명하자면 

For: 1쪽부터 10쪽까지 읽어야징
While: 읽은 페이지 수가 10이 될 때까지 읽어야징

이런 차이가 있다. ...사실 이렇게 말하면 모르시겠져? 그래서 가져왔음. 

둘 다 3^1~1^10까지 출력하는 코드인데(사실 저렇게 안하고 프린트문 줘도 됨)

For: 1부터 10까지 3에 제곱해
While: j가 있는데 이게 11보다 작을 동안 3에 제곱하고 하나씩 더해 

이런 식으로 돌아간다. For문은 범위를 주고 반복하는 반복문이고 While은 어떤 조건을 주고 반복하는 조건부 반복문. 그래서 For는 범위가 끝나면 반복문이 끝나고 While은 조건을 만족하면 반복문이 끝난다. 

참고로 베이직의 Do while loop와 비슷한건 While쪽이다. 

For ~ in range

본인은 보통 for i in range: 로 쓴다. 국민 알파벳 뭐 이런건가 바이오파이썬 할 때도 많이 봤던 그거 맞다. 고 때는 for record in records: 형태로 썼다.

for i in range(5):
    print(i)

이런 식으로 쓴다. 이건 0부터 4까지 출력하는 코드로, range는 별 지시가 없으면 0부터 시작해서 ~미만까지 잡는다. 

a=list(range(1,11))
for i in a: 
    print(i*i)

이미 리스트가 존재할때는 이런 식으로 쓰기도 한다. 

a="힘세고 강한 아침!"
for i in a: 
    print(i)

내가 그 얘기를 안해줄 뻔 했는데... 이거 문자열에도 먹힌다. 

While True

그냥 무한루프. 

i=1
while True: 
    i += 1
    if i > 10:
        break
print(i)

이게 While True를 적용한 간단한 코드인데, 밑에 있는 break는 뭐냐... While True는 무한루프라 사용자가 멈추지 않는 이상 계속 돌아간다. 그러면 ctrl+c를 누르거나 프로세스를 죽이지 않는 이상 계에에에에에에에에속 돌아간다 이 얘기. 저기 있는 if문과 break는 i에 1을 계에에에에속 더하다가 i가 10보다 크면 루프문을 빠져나와라 이 얘기다. 즉, while True는 break 없으면 계에에에에에에에에속 돌아간다. 

1. 공통: 이제 어디 자르는지도 나옵니다. 살려줘... 
2. Cutter: 제한효소 누락되던 거 수정했습니다. 근데 된건지 모르겠음. 

def count_func (a,b):
    while a in b:
        global site_count
        global res_loc 
        global res_loc_list
        loc = b.find(a)
        site_count += 1
        b = b[loc+len(a):]
        res_loc = loc + 1
        res_loc_list.append(res_loc)
    return site_count, res_loc
# 이거 통으로 코드에 넣었더니 if 안에 있는데도 시퀀스 없으면 끝내더라... 
# 위치 출력은 되는 것 같은데, 이거 더해야 하는데...

이렇게 하면 위치가 출력되긴 한데 문제가 하나 있다. 저 로직이 해당 시퀀스에서 제한효소 인식 site가 있으면 그 부분 다음 글자부터 slicing을 한다. 그래서 

AfaI [167, 150, 167, 53, 70] [830, 677, 507, 451, 378]
# restriction site 위치와 b(시퀀스에서 restriction site가 끝나는 부분부터 슬라이싱한 결과)의 길이

이렇게 나온다... (왼쪽이 find 결과) 오른쪽은 슬라이싱하고 남은 자리. 

def cut_func (a,b):
    while a in b:
        global res_loc
        global res_loc_list
        global length_seq
        global length_list
        loc = b.find(a)
        b = b[loc+len(a):]
        res_loc = loc + 1
        res_loc_list.append(res_loc)
        length_seq = len(b)
        length_list.append(length_seq)
    return res_loc_list,length_list
# 여기가 위치 관련 함수입니다.
# len(a): 인식 시퀀스 길이, len(b): 슬라이싱 하고 남은 시퀀스 길이

일단 복잡하니까 함수 밖으로 빼고 

AanI [80] [15]
# 결과(AanI의 시퀀스는 TTATAA로 6bp)
# 슬라이싱하고 남은 텍스트의 길이가 15bp이다.

저게 왜 저렇게 나오냐면 

1. AanI이 인식하는 시퀀스(TTATAA)는 6bp이다. (len(AanI)=6)
2. AanI이 인식하는 시퀀스 다음 위치부터 슬라이싱하면 15bp가 남는다. (왼쪽은 find에 +1한 값이다)

length_list.append(seq_length - (length_seq + len(a))) # slicing 후의 길이 목록

그래서 전체 길이-(슬라이싱+인식 시퀀스)+1하면 위치가 같게 나온다. (위 코드의 값에 1을 더하면 된다) 저렇게 하면 AanI의 경우 100-(15+6)+1이 되므로 80. 

def cut_func (a,b):
    while a in b:
        global res_loc # find로 나오는 값
        global res_loc_list
        seq_length = len(sequence)
        loc = b.find(a)
        b = b[loc+len(a):]
        res_loc = len(sequence) - (len(b) + len(a)) + 1
        res_loc_list.append(str(res_loc)) # find로 나오는 위치 목록(slicing에 따른 보정 필요)
    return res_loc_list,length_list
# 여기가 위치 관련 함수입니다.

res_loc_list = ', '.join(res_loc_list)
            f.write("{0}: {1} {2},{3} times cut. Where(bp): {4} \n".format(enzyme,res_find,feature,site_count,res_loc_list))

정리하고 출력 코드 바꾸면 된다. 

DB를 구축할 때 수동으로 구축했는데, CSV 구분자가 ,다보니 .로 들어간게 좀 있어서 그거 수정했음... (마른세수) 
효소 이름 오타난것도 수정했습니다... 

이래서 손끝 다치면 안됨... 

import pandas as pd
enzyme_table = pd.read_csv('/home/koreanraichu/restriction.csv')
# 이쪽은 Finder나 cutter에도 쓰이는 그 DB가 맞습니다. 
enzyme_table2 = pd.read_csv('/home/koreanraichu/restriction_RE.csv')
# 이쪽은 restriction site나 cut site에 A,T,G,C가 아닌 다른 알파벳이 들어가기 때문에 여기서 처음 불러오는 DB입니다. 
enzyme_table = pd.concat([enzyme_table,enzyme_table2])
enzyme_table = enzyme_table.sort_values('Enzyme')
enzyme_table.reset_index(inplace=True)
# 니네는 꼭 합치고 나면 인덱스도 바꿔줘야되더라...

사실 DB가 두갠데(둘다 GitHub에 올라가 있음), 하나는 Finder랑 Cutter에도 적용중이고 다른 하나는 적용하려면 별도의 처리가 필요해서(그게 RE 들어간 거) 여기에 처음 도입했음. 아직까지 통합DB가 없어서 불러온 다음 합치고 인덱싱 다시 하는 절차가 필요합니다. 

enzyme = input("찾고자 하는 효소를 입력해주세요: ")
# 아직 시퀀스로 검색하는 기능은 없습니다.

아직 효소 이름으로만 된다. (시퀀스로 검색하게 되면 밑에 로직도 수정해야 함)

def SeqtoString (a):
    a = enzyme_table.sequence[(enzyme_table['Enzyme'] == enzyme)]
    a = a.to_string(index = False)
    a = str(a)
    return a
def SitetoString (a):
    a = enzyme_table.restriction_site[(enzyme_table['Enzyme'] == enzyme)]
    a = a.to_string(index = False)
    a = str(a)
    return a

문자열화 해 주는 함수가 있고... 

for i in range(len(enzyme_table)):
    find_seq = SeqtoString(enzyme)
    Site_seq = SitetoString(enzyme_table['restriction_site'][i])
    DB_enzyme = enzyme_table['Enzyme'][i]
    DB_seq = enzyme_table['sequence'][i]
    DB_site = enzyme_table['restriction_site'][i]
    if find_seq == str(DB_seq):
        print(DB_enzyme,DB_seq,DB_site)
    else: 
        pass

일단 이렇게만 하고 HaeIII으로 검색해 본 결과. 잘 되는데, 입력한 효소는 빼는 게 맞지 않나... (시퀀스로 검색하게 되면 이 부분도 바뀐다)

print("{0} | {1} | {2} | Input enzyme".format(enzyme,find_seq,Site_seq))
for i in range(len(enzyme_table)):
    find_seq = SeqtoString(enzyme)
    Site_seq = SitetoString(enzyme)
    DB_enzyme = enzyme_table['Enzyme'][i]
    DB_seq = enzyme_table['sequence'][i]
    DB_site = enzyme_table['restriction_site'][i]
    if find_seq == str(DB_seq) and DB_enzyme != enzyme:
        if Site_seq == DB_site:
            print("{0} | {1} | {2} | Isoschizomer".format(DB_enzyme,DB_seq,DB_site))
            # 인식하는 시퀀스와 자르는 방식이 같은 제한효소
        else: 
            print("{0} | {1} | {2} | Neoschizomer".format(DB_enzyme,DB_seq,DB_site))
            # 인식하는 시퀀스는 같으나 자르는 방식이 다른 제한효소
    elif find_seq == str(DB_seq) and DB_enzyme == enzyme:
        pass
    else: 
        pass

Isoschizomer와 Neoschizomer를 구별해준다. (입력한 효소는 물론 뺐다) 시퀀스로 검색하는 기능을 넣게 되면, 이 부분도 로직이 바뀌게 된다. (시퀀스로 검색할 경우 해당 시퀀스를 인식하는 제한효소를 다 뽑아와야 함) 

Iso = ', '.join(Iso)
Neo = ', '.join(Neo)
print("Isoschizomer: {0} \nNeoschizomer: {1}".format(Iso,Neo))
# 실제로 Isoschizomer인데도 Neoscizomer로 표기하는 문제가 있습니다. (BamHI-Nsp29132OO)

물론 본인쟝은 친절하게 텍스트로도 출력해드림. (얘는 따로 데이터 저장할 생각은 없음)

참고

Isoschizomer: 인식하는 시퀀스와 자르는 형태가 같은 제한효소
Neoschizomer: 인식하는 시퀀스는 같은데 자르는 형태가 다른 제한효소(예: GGCC를 인식하지만 한놈은 GG/CC로 자르고 한놈은 /GGCC로 자르는 것)

keyword = input("효소 이름으로 찾으실거면 enzyme을, restriction site sequence로 찾으실거면 sequence를 입력해주세요. ")
if keyword == "enzyme":
    enzyme = input("찾고자 하는 효소를 입력해주세요: ")
elif keyword == "sequence":
    seq = input("찾고자 하는 restriction site sequence를 입력해주세요: ")
else: 
    print("다시 입력해주세요. ")
# 효소 이름으로 찾느냐, 시퀀스로 찾느냐에 따라 검색 결과가 다릅니다.

현재 효소 이름과 인식하는 시퀀스로 찾는 기능을 지원함. (둘 중 하나로만 찾을 수 있습니다)

if keyword == "enzyme":
    find_seq = SeqtoString(enzyme)
    Site_seq = SitetoString(enzyme)
    Iso = []
    Neo = []
    print("{0} | {1} | {2} | Input enzyme".format(enzyme,find_seq,Site_seq))
    for i in range(len(enzyme_table)):
        DB_enzyme = str(enzyme_table['Enzyme'][i]).strip()
        DB_seq = str(enzyme_table['sequence'][i]).strip().upper()
        DB_site = str(enzyme_table['restriction_site'][i]).strip().upper()
        if find_seq == str(DB_seq) and DB_enzyme != enzyme:
            if Site_seq == DB_site:
                Iso.append(DB_enzyme)
                print("{0} | {1} | {2} | Isoschizomer".format(DB_enzyme,DB_seq,DB_site))
                # 인식하는 시퀀스와 자르는 방식이 같은 제한효소
            elif Site_seq != DB_site: 
                Neo.append(DB_enzyme)
                print("{0} | {1} | {2} | Neoschizomer".format(DB_enzyme,DB_seq,DB_site))
                # 인식하는 시퀀스는 같으나 자르는 방식이 다른 제한효소
        elif find_seq == str(DB_seq) and DB_enzyme == enzyme:
            pass
        else: 
            pass
# 여기까지는 효소 이름으로 검색할 때의 코드

이 쪽은 효소 이름으로 찾을 때 로직이고 

else: 
    find_seq = seq
    print("Searched by: {0}".format(seq))
    for i in range(len(enzyme_table)):
        DB_enzyme = str(enzyme_table['Enzyme'][i]).strip()
        DB_seq = str(enzyme_table['sequence'][i]).strip().upper()
        DB_site = str(enzyme_table['restriction_site'][i]).strip().upper()
        if find_seq == DB_seq:
            print("{0} | {1} | {2}".format(DB_enzyme,DB_seq,DB_site))
        else:
            pass
# 여기까지는 인식 시퀀스로 검색할 때의 코드

이 쪽은 시퀀스로 찾는 로직입니다. 

별개로 sticky만, blunt만 찾는 것도 있으면 좋을 것 같긴 함. 

공통패치

Cut수 세 주는 기능이 추가되었습니다. 

Cutter

Cut수에 따른 효소 리스트 업 기능이 추가되었습니다. 살려줘... (아직 자르는 위치 안했음)
파일 이름 형식이 변경되었습니다. 그래서 이제 시퀀스 이름도 받습니다. 

Finder

Cut수 세 주는 기능에 따른 출력 형식 수정이 있었습니다. 

Cut수 세 주는 기능

사실상 엄청난 노가다의 결과...ㅠㅠ 

(생략)
        if res_find in sequence:
            site_count = 0
            while sequence.find(res_find) != -1:
                loc = sequence.find(res_find)
                site_count += 1
                sequence = sequence[loc+len(res_find):]
                print(enzyme, res_find, sequence.find(res_find))
            f.write("{0}: {1} {2},{3} times cut.\n".format(enzyme,res_find,feature,site_count))
(생략)

이 코드가 단식으로는 되는데 전체 코드에 쌩으로 도입했더니, restriction site가 없는 효소가 나오면 거기서 멈춰버린다. (여기까지가 전날 저녁에 try했던 부분)

def count_site (a,b):
    site_count = 0
    while a in b:
        loc - b.find(a)
        site_count += 1
        b = b[loc+len(a):]
        return site_count

그래서 해당 로직을 아예 함수로 빼버렸다. 

(마른세수)

def count_site (a,b):
    site_count = 0
    while a in b:
        loc = b.find(a)
        site_count += 1
        b = b[loc+len(a):]
        return site_count

근데 생각해보니 네번째 줄에 -가 왜 들어감? (원래 =)

실례지만 0컷이면 안 나와야 정상 아니냐. 

def count_func (a,b):
    site_count = 0
    while a in b:
        loc = b.find(a)
        site_count += 1
        b = b[loc+len(a):]
    return site_count
# 이거 통으로 코드에 넣었더니 if 안에 있는데도 시퀀스 없으면 끝내더라... 
print(count_site("GGCC",sequence))

답답해서 직접 해봤더니 이건 잘 되더라. 

def count_func (a,b):
    while a in b:
        global site_count
        loc = b.find(a)
        site_count += 1
        b = b[loc+len(a):]
    return site_count
# 이거 통으로 코드에 넣었더니 if 안에 있는데도 시퀀스 없으면 끝내더라...

함수 내에서는 그냥 전역변수 선언하고 

(생략)
        if res_find in sequence:
            site_count = 0
            count_func(res_find,sequence)
            count += 1
            f.write("{0}: {1} {2},{3} times cut.\n".format(enzyme,res_find,feature,site_count))
        else: 
            count += 0
(생략)
# 여러분 드디어 저장기능이 추가되었습니다!!!

site_count 변수를 함수 밖으로 빼버렸다. (site_count는 if문 안쪽에 있다)

(흐-뭇)

Cut수별로 나눠주는 기능

(생략)
        if res_find in sequence:
            site_count = 0
            count_func(res_find,sequence)
            count += 1
            count_nocut += 0
            cut_list.append(enzyme)
            f.write("{0}: {1} {2},{3} times cut.\n".format(enzyme,res_find,feature,site_count))
        else: 
            count += 0
            count_nocut += 1
            nocut_list.append(enzyme)
    cut_list = ', '.join(cut_list)
    nocut_list = ', '.join(nocut_list)
    f.write("Total: {0} enzymes cut input sequence, {1} enzymes never cut this sequence. \n".format(count,count_nocut))
    f.write("Enzyme cuts: {0} \nEmzyme no cuts: {1}".format(cut_list,nocut_list))
    f.close()
# 진짜 세주는거 겨우 추가했습니다...ㅠㅠ

사실 길어서 생략했는데 코드 위쪽에 0컷이랑 컷 두 개의 리스트가 있다. 

count = 0
count_nocut = 0
cut_list = []
nocut_list = []

이놈들. (count는 원래 있던 변수고 그 밑에 있는 게 새로 추가한 변수)

그래서 이렇게 나오는데... 이게 사람이요... 살다보면 1컷 궁금하다... 

count = 0
count_nocut = 0
once_cut_list = []
multi_cut_list = []
nocut_list = []
(생략)
        if res_find in sequence:
            site_count = 0
            count_func(res_find,sequence)
            count += 1
            count_nocut += 0
            if site_count == 1:
                once_cut_list.append(enzyme)
            else: 
                multi_cut_list.append(enzyme)
            f.write("{0}: {1} {2},{3} times cut.\n".format(enzyme,res_find,feature,site_count))
        else: 
            count += 0
            count_nocut += 1
            nocut_list.append(enzyme)
    once_cut_list = ', '.join(once_cut_list)
    multi_cut_list = ', '.join(multi_cut_list)
    nocut_list = ', '.join(nocut_list)
(생략)

그래서 0컷 1컷 멀티컷 나눠드렸습니다. 

가 이거. 근데... 에픽하이에도 투컷이 있고... (fly랑 one 좋아함) NEB cutter도 투컷까지는 보여줘요... 

count = 0
count_nocut = 0
once_cut_list = []
two_cut_list = []
multi_cut_list = []
nocut_list = []
(생략)
        if res_find in sequence:
            site_count = 0
            count_func(res_find,sequence)
            count += 1
            count_nocut += 0
            if site_count == 1:
                once_cut_list.append(enzyme)
            elif site_count == 2: 
                two_cut_list.append(enzyme)
            else: 
                multi_cut_list.append(enzyme)
            f.write("{0}: {1} {2},{3} times cut.\n".format(enzyme,res_find,feature,site_count))
        else: 
            count += 0
            count_nocut += 1
            nocut_list.append(enzyme)
(생략)

나눠드렸습니다^^ (출력은 1컷 2컷 멀티컷 빵컷 순) 빵컷하니까 빵 자르는 것 같잖아요 

def count_func (a,b):
    while a in b:
        global site_count
        loc = b.find(a)
        site_count += 1
        b = b[loc+len(a):]
    return site_count
# Cutter test하다가 여기에도 추가했음...

Finder에도 이 함수는 똑같이 들어간다. 

with open ('Result_{0}-{1}-{2}_{3}-{4}.txt'.format(year,month,day,enzyme,sequence_name),'w',encoding='utf-8') as f: 
    if sequence.find(res_find) != -1:
        site_count = 0
        cut_count = count_func(res_find,sequence)
        sequence = sequence.replace(res_find,res_site)
        print(enzyme,",",cut_feature)
        print(sequence,cut_count)
        f.write("{0} | {1} | {2} | {3} times cut\n".format(enzyme,res_site,cut_feature,cut_count))
        f.write('Sequence name: {0} \n {1}'.format(sequence_name,sequence))
        f.close()
        # DB에 효소가 있고 일치하는 시퀀스가 있을 때
(하략)

대신 Finder는 DB에 효소가 있고, 효소가 해당 시퀀스를 자를 때만 cut수를 세 준다. 

그래서 이렇게 나온다. (HaeIII이 해당 시퀀스를 자른다)

Finder

with open ('Result_{0}-{1}-{2}_{3}-{4}.txt'.format(year,month,day,enzyme,search_sequence_name),'w',encoding='utf-8') as f:
    if search_sequence.find(res_find) != -1:
        search_sequence = search_sequence.replace(res_find,res_site)
        print(enzyme,",",cut_feature)
        print(search_sequence)
        f.write("{0} | {1} | {2} \n".format(enzyme,res_site,cut_feature))
        f.write('Sequence name: {0} \n'.format(search_sequence_name))
        f.write(search_sequence)
        f.close()
        # DB에 효소가 있고 일치하는 시퀀스가 있을 때
    elif enzyme_table['Enzyme'].isin([enzyme]).any() == True and search_sequence.find(res_find) == -1:
        print("No restriction site in this sequence. ")
        f.write("{0} | {1} | {2} \n".format(enzyme,res_site,cut_feature))
        f.write('Sequence name: {0} \n'.format(search_sequence_name))
        f.write("This restricion enzyme never cut this sequence. ")
        f.close()
        # DB에 효소가 있으나 일치하는 시퀀스가 없을 때
    else:
        print("No data in database. ")
        f.write("{0} \n".format(enzyme))
        f.write("This restriction enzyme not entried in database. ")
        f.close()
        # DB에 효소가 없을 때

저장 형식은 동일하고, 저장용 코드가 간소화 되었습니다. 저게 format이 들을 줄은 몰랐지... 

Cutter

filter = input("sticky로 자르는 제한효소만 보고 싶으면 sticky, blunt로 자르는 제한효소만 보고 싶으면 blunt를 입력해주세요. ")

Filter 기능이 추가되었습니다. 거창한 건 아니고, sticky end나 blunt end만 보고 싶을 때 쓰시면 됩니다. 

if filter == 'sticky':
    enzyme_table = enzyme_table[enzyme_table['cut_feature']== 'sticky']
    enzyme_table.reset_index(inplace=True)
elif filter == 'blunt':
    enzyme_table = enzyme_table[enzyme_table['cut_feature']== 'blunt']
    enzyme_table.reset_index(inplace=True)
else: 
    pass

참고로 filter 입력 여부에 따라 DB에서 해당 조건으로만 테이블을 재구성합니다. (아무것도 안 입력하면 패스)

count = 0
with open('Result.txt','w',encoding='utf-8') as f:
    f.write("Restriction enzyme which cuts this sequence: ")
    for i in range(len(enzyme_table)):
        enzyme = enzyme_table['Enzyme'][i]
        res_find = enzyme_table['sequence'][i]
        res_find = str(res_find)
        if res_find in sequence:
            print(enzyme, res_find, sequence.find(res_find))
            f.write("{0}: {1} \n".format(enzyme,res_find))
            count += 1
        else: 
            count += 0
    print(count)
    f.write("Total: {0} enzymes cut input sequence".format(count))
# 아직 저장기능은 없습니다. 지금 출력도 좀 중구난방이라 정리 좀 해야될듯. 
# find로 나오는 위치의 경우 0부터 시작하기떄문에 하나 더해줬습니다. 아울러 해당 메소드가 '가장 처음에 나오는 글자'만 찾아주는거지 전체 검색이 아니기때문에 여러군데를 자르는지 여부는 모릅니다.

아 출력 파일도 만들어준다구! 형식은 (날짜)-filter입니다. (아무것도 선택하지 않을 경우 filter=None)

참고로 여기에 쓰이는 DB와 쓰이지는 않는 DB 둘 다 구축 끝냈습니다. (쓰이지 않는 DB: 인식 혹은 자르는 시퀀스에 N, R, Y같은 게 들어가서 따로 처리해야 함) 

사실 전에 만든 코드 이름은 Finder고, 이놈이 커터임. 

import pandas as pd

이 코드도 판다스가 있어야된다. (DB가 csv)

enzyme_table = pd.read_csv('/home/koreanraichu/restriction.csv')
enzyme_table = enzyme_table.sort_values('Enzyme')
# Finder에도 쓰이는 '그' DB 맞습니다. 현재 수동 구축 중... 
print(enzyme_table)
print(len(enzyme_table))

아직도 갈 길이 멀지만 일단 D까지 추가했음... 아울러 py파일은 print가 빠집니다. 

sequence = input("검색할 시퀀스를 입력해주세요: ")

이건 시퀀스 입력받는 코드(아직 이름은 안 받음)

for i in range(len(enzyme_table)):
    res_find = enzyme_table['sequence'][i]
    res_find = str(res_find)
    if res_find in sequence:
        print(enzyme_table['Enzyme'][i],res_find,sequence.find(res_find))
    else: 
        print(enzyme_table['Enzyme'][i],"Not found")

문자열화가 for문 안에서 이루어지고, 자르는 효소와 자르지 않는 효소의 출력 방식만 다르다. 현재 저장도 지원 안 해줌. 

패치노트

1. 출력 파일명 형식 변경 저기도 format이 먹힐 줄은 몰랐음 
2. 출력 파일의 형식 변경(시퀀스 이름 추가)
3. 깃헙에 해당 코드 py파일도 추가됨

패치 결과

원래 날짜만 추가했었는데 생각해보니 날짜가 같으면 헷갈릴 것 같아서 효소랑 시퀀스 이름도 추가함

시퀀스 이름은 출력 파일에도 저장됩니다. (그래서 이제 세줄임) 

근데 NEB커터 쓰세여 그거 좋음 

일단 이 코드는 이 시퀀스를 자르는 제한효소들을 찾는 게 아님. 이 효소가 이 시퀀스를 자르는가? 를 보는 코드임다. 이 점 유념해주세요. 그리고 이거 올리면서 Jupyter가 매트랩이랑 비중 같아졌다 

코드

Jupyter notebook으로 코딩한거고 나중에 일부 블록은 정리할 예정. 

import pandas as pd

구축한 csv파일을 가져오고 취급하려면 얘가 필요하다. csv파일은 혹시 써보실 분 계시면 말씀주세요. 

참고로 csv파일이 되게 단촐해서 효소 이름, 인식하는 시퀀스, 자르는 시퀀스, 자르는 형태(sticky or blunt)가 들어가 있음. 나중에 여건이 된다면 똑같은 부분을 인식하고 자르는 다른 효소나 처리온도에 대한 정보도 추가할 예정. (물론 출력 파일에도 같이 나오게 해야 하는 게 함정)

enzyme_table = pd.read_csv('/home/koreanraichu/restriction.csv') # 일단 자체적으로 구축했음... 저거 종류 개많습니다 ㅠㅠ 
enzyme_table = enzyme_table.sort_values('Enzyme')
# csv파일의 구성은 크게 효소 이름, 인식하는 시퀀스, 해당 시퀀스를 자르는 형태와 sticky or blunt 여부로 구성되어 있습니다.

위에서 얘기한 csv파일을 불러와서 효소 이름순으로 정렬합니다. 참고로 효소가 인식하는 시퀀스에 N, W, D, S같은 거(ATGC 말고 다른 거) 들어가 있는 효소는 일단 뺐음. 걔들은 따로 처리해야 하는데 그것까지는 아직 무리입니다. 

enzyme=input('시퀀스를 찾을 제한효소를 입력해주세요: ')
search_sequence=input('제한효소 site를 찾을 시퀀스를 입력해주세요: ')

이거는 입력받는 코드. 효소 이름은 토씨 하나 안 틀리고 써야 한다. 

if enzyme_table['Enzyme'].isin([enzyme]).any() == True:
    res_find = enzyme_table.sequence[(enzyme_table['Enzyme'] == enzyme)]
    res_find = res_find.to_string(index=False)
    res_find = str(res_find)
    print(res_find)
    # 효소 이름이 데이터베이스에 있을 경우 검색할 시퀀스 데이터를 가져온다
    res_site = enzyme_table.restriction_site[(enzyme_table['Enzyme'] == enzyme)]
    res_site = res_site.to_string(index=False)
    res_site = str(res_site)
    print(res_site)
    # 효소 이름이 데이터베이스에 있을 경우 검색하고 대체할 시퀀스 데이터를 가져온다
    cut_feature = enzyme_table.cut_feature[(enzyme_table['Enzyme'] == enzyme)]
    cut_feature = cut_feature.to_string(index=False)
    cut_feature = str(cut_feature)
    print(cut_feature)
    # blunt or sticky(나중에 저장 기능 추가할 때 넣을 예정입니다)
else: 
    print("No data in Database")

판다스 종특인지는 모르겠으나... to_string() 없으니까 인식 못하데... 

순서대로 '인식하는 시퀀스', '자르는 시퀀스', '자르는 형태' 컬럼. 이 코드 자체는 가져온 DB에 효소 이름이 있어야 진행하므로 DB에 효소가 없을 때는 No data가 뜬다. (참고로 수동으로 구축중이라 시간 무지하게 걸림)

if enzyme_table['Enzyme'].isin([enzyme]).any() == True:
    print(search_sequence.find(str(res_find)))
else: 
    pass
# 여기는 검색결과가 존재하지 않으면 -1로 나옵니다. (윗 블럭이랑 여기는 넘어가도 되는 부분)

이 코드도 마찬가지로 DB에 효소가 있을 때 시퀀스를 찾고, 결과를 출력한다. (find.()는 문자열이 없으면 -1로 뽑음)

with open ('Result.txt','w',encoding='utf-8') as f: 
    if search_sequence.find(res_find) != -1:
        search_sequence = search_sequence.replace(res_find,res_site)
        print(enzyme,",",cut_feature)
        print(search_sequence)
        f.write(enzyme)
        f.write(", ")
        f.write(res_site)
        f.write(", ")
        f.write(cut_feature)
        f.write("\n")
        f.write(search_sequence)
        f.close()
        # DB에 효소가 있고 일치하는 시퀀스가 있을 때
    elif enzyme_table['Enzyme'].isin([enzyme]).any() == True and search_sequence.find(res_find) == -1:  
        print("No restriction site in this sequence. ")
        f.write(enzyme)
        f.write(", ")
        f.write(res_site)
        f.write(", ")
        f.write(cut_feature)
        f.write("\n")
        f.write("This restricion enzyme never cut this sequence. ")
        f.close()
        # DB에 효소가 있으나 일치하는 시퀀스가 없을 때
    else:
        print("No data in database. ")
        f.write(enzyme)
        f.write("\n")
        f.write("This restriction enzyme not entried in database. ")
        f.close()
        # DB에 효소가 없을 때

출력과 저장에 대한 코드. 크게 DB에 효소가 있나/없나로 갈리고 효소가 있을 때 인식 시퀀스가 있나/없나로 갈린다. 

그래서 각각 

이런 식으로 출력된다. 

앞으로 추가할 기능

1. csv 파일 수정에 따른 출력 형식(현재 효소 이름, 자르는 시퀀스, sticky ot blunt 여부만 있는데 처리 온도도 및 같은 시퀀스를 인식하는 다른 효소의 정보도 추가할 예정)

2. 출력 파일명도 형식을 줄 수 있으면 날짜 이런거 들어가면 좋을 것 같음. (아니면 사용자가 입력할 수 있게 하거나) 일단 오늘 날짜 가져오는 게 먼저긴 한데 아니면 하이브리드로 입력 받아서 오늘 날짜랑 같이 묶는 뭐 이런거? 

멀티랑 N, W, K 처리는... 아 그건 제 능력 밖이니까 걍 NEB 커터 쓰세여... 그거 좋음... 

text = []
while True: 
    input_text = input("wordcloud로 만들 텍스트를 입력해주세요. ")
    text.append(input_text) 
    if input_text == "":
        break
text = ' '.join(text)
text=okt.nouns(text)
text = ' '.join(text)

크게는 저 부분을 수정했고, 전체적으로 wordcloud 만들기 위해 입력받는 부분도 간소화했음. (While True 주고 입력 없으면 break 하도록) 

참고로 Wordcloud는 안에 들어가는 글자가 많을 수록 멋지게 나옵니다. 이 점 유념하시길. 

3부작으로 나눌 걸 그랬나... 

라이브러리 다 깔았지들? 

Logistic regression(로지스틱 회귀분석)

> dat=subset(mtcars,select=c(mpg,am,vs))
> dat
                     mpg am vs
Mazda RX4           21.0  1  0
Mazda RX4 Wag       21.0  1  0
Datsun 710          22.8  1  1
Hornet 4 Drive      21.4  0  1
Hornet Sportabout   18.7  0  0
Valiant             18.1  0  1
Duster 360          14.3  0  0
Merc 240D           24.4  0  1
Merc 230            22.8  0  1
Merc 280            19.2  0  1
Merc 280C           17.8  0  1
Merc 450SE          16.4  0  0
Merc 450SL          17.3  0  0
Merc 450SLC         15.2  0  0
Cadillac Fleetwood  10.4  0  0
Lincoln Continental 10.4  0  0
Chrysler Imperial   14.7  0  0
Fiat 128            32.4  1  1
Honda Civic         30.4  1  1
Toyota Corolla      33.9  1  1
Toyota Corona       21.5  0  1
Dodge Challenger    15.5  0  0
AMC Javelin         15.2  0  0
Camaro Z28          13.3  0  0
Pontiac Firebird    19.2  0  0
Fiat X1-9           27.3  1  1
Porsche 914-2       26.0  1  0
Lotus Europa        30.4  1  1
Ford Pantera L      15.8  1  0
Ferrari Dino        19.7  1  0
Maserati Bora       15.0  1  0
Volvo 142E          21.4  1  1

이거 뭐 내장 데이터라고 함... 

예측 변수는 연속적인데 결과가 이분법(얘 아니면 쟤)

> logr_vm=glm(vs~mpg,data=dat,family=binomial)
> logr_vm

Call:  glm(formula = vs ~ mpg, family = binomial, data = dat)

Coefficients:
(Intercept)          mpg  
    -8.8331       0.4304  

Degrees of Freedom: 31 Total (i.e. Null);  30 Residual
Null Deviance:	    43.86 
Residual Deviance: 25.53 	AIC: 29.53

> logr_vm <- glm(vs ~ mpg, data=dat, family=binomial(link="logit"))
> logr_vm

Call:  glm(formula = vs ~ mpg, family = binomial(link = "logit"), data = dat)

Coefficients:
(Intercept)          mpg  
    -8.8331       0.4304  

Degrees of Freedom: 31 Total (i.e. Null);  30 Residual
Null Deviance:	    43.86 
Residual Deviance: 25.53 	AIC: 29.53

일단 둘 중 하나 적용하면 된다. 어차피 결과는 같아서... 

> summary(logr_vm)

Call:
glm(formula = vs ~ mpg, family = binomial(link = "logit"), data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.2127  -0.5121  -0.2276   0.6402   1.6980  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -8.8331     3.1623  -2.793  0.00522 **
mpg           0.4304     0.1584   2.717  0.00659 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 43.860  on 31  degrees of freedom
Residual deviance: 25.533  on 30  degrees of freedom
AIC: 29.533

Number of Fisher Scoring iterations: 6

이쪽은 회귀분석과 마찬가지로 summary()로 더보기가 된다. 같은 회귀분석이라 그런가 쟤는 선형인데?

> library(ggplot2)
> ggplot(dat, aes(x=mpg, y=vs)) + geom_point() + 
+   stat_smooth(method="glm", method.args=list(family="binomial"), se=FALSE)
`geom_smooth()` using formula 'y ~ x'

처음에 geom 보고 점이라 저렇게 쓴건가 했음... ㅋㅋㅋㅋㅋㅋ 아 이 그래프는 그리시려면 ggplot2 깔아야됨. 어차피 다음 챕터도 저거니까 걍 지금 까세요. 

> par(mar=c(4,4,1,1))
> plot(dat$mpg,dat$vs)
> curve(predict(logr_vm,data.frame(mpg=x),type="response"),add=TRUE)

이거는 ggplot2 없어도 될건데... 

예측변수와 결과가 둘 다 이분법(둘다 얘 아니면 쟤) 무슨 시그모이드냐 

> logr_va=glm(vs~am,data=dat,family=binomial)
> logr_va

Call:  glm(formula = vs ~ am, family = binomial, data = dat)

Coefficients:
(Intercept)           am  
    -0.5390       0.6931  

Degrees of Freedom: 31 Total (i.e. Null);  30 Residual
Null Deviance:	    43.86 
Residual Deviance: 42.95 	AIC: 46.95

근데 am은 뭐임? am pm은 아닐것인데... 각도도 아니고. 

> ggplot(dat, aes(x=am, y=vs)) + 
+   geom_point(shape=1, position=position_jitter(width=.05,height=.05)) + 
+   stat_smooth(method="glm", method.args=list(family="binomial"), se=FALSE)
`geom_smooth()` using formula 'y ~ x'

뭐야 내 인테그랄 돌려줘요 로지스틱 의문의 적분기호행 

> par(mar=c(4,4,1,1))
> plot(jitter(dat$am,.2),jitter(dat$vs,.2))
> curve(predict(logr_va, data.frame(am=x), type="response"), add=TRUE)

(마른세수)

예측 데이터가 이분법+연속형 짬뽕이고 결과는 이분법인 예

> logr_vma=glm(vs~mpg+am,data=dat,family=binomial)
> logr_vma

Call:  glm(formula = vs ~ mpg + am, family = binomial, data = dat)

Coefficients:
(Intercept)          mpg           am  
   -12.7051       0.6809      -3.0073  

Degrees of Freedom: 31 Total (i.e. Null);  29 Residual
Null Deviance:	    43.86 
Residual Deviance: 20.65 	AIC: 26.65

(대충 이렇게 분석하면 된다는 얘기) 

예측 데이터 여러개와 상호작용(Interaction)

> logr_vmai=glm(vs~mpg*am,data=dat,family=binomial)
> logr_vmai

Call:  glm(formula = vs ~ mpg * am, family = binomial, data = dat)

Coefficients:
(Intercept)          mpg           am       mpg:am  
   -20.4784       1.1084      10.1055      -0.6637  

Degrees of Freedom: 31 Total (i.e. Null);  28 Residual
Null Deviance:	    43.86 
Residual Deviance: 19.12 	AIC: 27.12

> logr_vmai=glm(vs~mpg+am+mpg:am,data=dat,family=binomial)
> logr_vmai

Call:  glm(formula = vs ~ mpg + am + mpg:am, family = binomial, data = dat)

Coefficients:
(Intercept)          mpg           am       mpg:am  
   -20.4784       1.1084      10.1055      -0.6637  

Degrees of Freedom: 31 Total (i.e. Null);  28 Residual
Null Deviance:	    43.86 
Residual Deviance: 19.12 	AIC: 27.12

둘 중 아무거나 쓰면 된다. 근데 얘는 왜 그래프 없음? 

등분산 검정

> head(InsectSprays)
  count spray
1    10     A
2     7     A
3    20     A
4    14     A
5    14     A
6    12     A

이거랑 

> tg=ToothGrowth
> tg$dose=factor(tg$dose)
> head(tg)
   len supp dose
1  4.2   VC  0.5
2 11.5   VC  0.5
3  7.3   VC  0.5
4  5.8   VC  0.5
5  6.4   VC  0.5
6 10.0   VC  0.5

이거 쓸 건데 독립변수가 하나냐 멀티냐 차이다. 둘 다 내장 데이터인듯. 

> plot(count~spray,data=InsectSprays)

> plot(len ~ interaction(dose,supp), data=ToothGrowth)

각 데이터별 박스 그래프는 대충 이런 모양이다. 

Bartlett’s test

> bartlett.test(count~spray,data=InsectSprays)

	Bartlett test of homogeneity of variances

data:  count by spray
Bartlett's K-squared = 25.96, df = 5, p-value = 9.085e-05

독립 변수가 하나일 때

> bartlett.test(len~interaction(supp,dose),data=ToothGrowth)

	Bartlett test of homogeneity of variances

data:  len by interaction(supp, dose)
Bartlett's K-squared = 6.9273, df = 5, p-value = 0.2261

독립 변수가 여러개일 때(interaction()으로 합쳐야 한다)

Levene’s test

이거 쓰려면 뭐 또 깔아야됨... car 깔고 오세여. 

> leveneTest(count~spray,data=InsectSprays)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value   Pr(>F)   
group  5  3.8214 0.004223 **
      66                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

독립 변수가 하나일 때

> leveneTest(len~supp*dose,data=tg)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  5  1.7086 0.1484
      54

독립 변수가 여러개일 때

Fligner-Killeen test

> fligner.test(count~spray,data=InsectSprays)

	Fligner-Killeen test of homogeneity of variances

data:  count by spray
Fligner-Killeen:med chi-squared = 14.483, df = 5, p-value = 0.01282

독립 변수가 하나일 때

> fligner.test(len~interaction(supp,dose),data=ToothGrowth)

	Fligner-Killeen test of homogeneity of variances

data:  len by interaction(supp, dose)
Fligner-Killeen:med chi-squared = 7.7488, df = 5, p-value = 0.1706

독립 변수가 여러개일 때(얘도 interaction()으로 묶는다)

Inter-rater reliability

네 irr 깔고 오시면 되겠습니다. 

> library(irr)
필요한 패키지를 로딩중입니다: lpSolve
> data(diagnoses)
> dat=diagnoses[,1:3]
> head(dat)
                   rater1                  rater2                  rater3
1             4. Neurosis             4. Neurosis             4. Neurosis
2 2. Personality Disorder 2. Personality Disorder 2. Personality Disorder
3 2. Personality Disorder        3. Schizophrenia        3. Schizophrenia
4                5. Other                5. Other                5. Other
5 2. Personality Disorder 2. Personality Disorder 2. Personality Disorder
6           1. Depression           1. Depression        3. Schizophrenia

그리고 이걸 또 데려옴... 

Kohen's kappa

> kappa2(dat[,c(1,2)],"unweighted")
 Cohen's Kappa for 2 Raters (Weights: unweighted)

 Subjects = 30 
   Raters = 2 
    Kappa = 0.651 

        z = 7 
  p-value = 2.63e-12

응? 가중치? 이거는 나중에 설명드림. 일단 이렇게 하시면 됩니다. 

N raters: Fleiss’s Kappa, Conger’s Kappa

> kappam.fleiss(dat)
 Fleiss' Kappa for m Raters

 Subjects = 30 
   Raters = 3 
    Kappa = 0.534 

        z = 9.89 
  p-value = 0

Fleiss's kappa

> kappam.fleiss(dat,exact=TRUE)
 Fleiss' Kappa for m Raters (exact value)

 Subjects = 30 
   Raters = 3 
    Kappa = 0.55

Conger's exact kappa

Kohen's kappa(weighted)

데이터가 순서가 있는 데이터일 경우 가중치가 추가된다고... 

> data(anxiety)
> dfa=anxiety[,c(1,2)]
> head(dfa)
  rater1 rater2
1      3      3
2      3      6
3      3      4
4      4      6
5      5      2
6      5      4

이 데이터를 쓸 예정이다. ...여기 순서가 있었음? 

> kappa2(dfa,"squared")
 Cohen's Kappa for 2 Raters (Weights: squared)

 Subjects = 20 
   Raters = 2 
    Kappa = 0.297 

        z = 1.34 
  p-value = 0.18

Weight=squared

> kappa2(dfa,"equal")
 Cohen's Kappa for 2 Raters (Weights: equal)

 Subjects = 20 
   Raters = 2 
    Kappa = 0.189 

        z = 1.42 
  p-value = 0.157

Weight=equal

> kappa2(dfa,"unweighted")
 Cohen's Kappa for 2 Raters (Weights: unweighted)

 Subjects = 20 
   Raters = 2 
    Kappa = 0.119 

        z = 1.16 
  p-value = 0.245

가중치 없음

Kohen's kappa(weighted)+factor

팩터: 레벨이 있는 그거 맞음

> dfa2=dfa
> dfa2$rater1=factor(dfa2$rater1,levels=1:6,labels=LETTERS[1:6])
> dfa2$rater2=factor(dfa2$rater2,levels=1:6,labels=LETTERS[1:6])
> head(dfa2)
  rater1 rater2
1      C      C
2      C      F
3      C      D
4      D      F
5      E      B
6      E      D

위에서 불러 온 데이터들을 팩터화한 다음, 레벨을 매겼다. 

> levels(dfa2$rater1)
[1] "A" "B" "C" "D" "E" "F"
> levels(dfa2$rater2)
[1] "A" "B" "C" "D" "E" "F"

그래서 이렇게 됐다. (아마 label이 letter 1~6이라 알파벳으로 표기된 듯)

> kappa2(dfa2,"squared")
 Cohen's Kappa for 2 Raters (Weights: squared)

 Subjects = 20 
   Raters = 2 
    Kappa = 0.297 

        z = 1.34 
  p-value = 0.18

Weight=squared

> kappa2(dfa2,"equal")
 Cohen's Kappa for 2 Raters (Weights: equal)

 Subjects = 20 
   Raters = 2 
    Kappa = 0.189 

        z = 1.42 
  p-value = 0.157

Weight=equal

> kappa2(dfa2,"unweighted")
 Cohen's Kappa for 2 Raters (Weights: unweighted)

 Subjects = 20 
   Raters = 2 
    Kappa = 0.119 

        z = 1.16 
  p-value = 0.245

가중치 없음

ICC(Intraclass correlation coefficient)

이건 또 무슨 저세상 용어냐...OTL 

> icc(anxiety,model="twoway",type="agreement")
 Single Score Intraclass Correlation

   Model: twoway 
   Type : agreement 

   Subjects = 20 
     Raters = 3 
   ICC(A,1) = 0.198

 F-Test, H0: r0 = 0 ; H1: r0 > 0 
 F(19,39.7) = 1.83 , p = 0.0543 

 95%-Confidence Interval for ICC Population Values:
  -0.039 < ICC < 0.494

정말 이건 무슨 저세상 용어인가요...