일단... 

해당 기능 추가 결과물이고요... 정규식 얘기는 나중에 입 털어드림... 

import re # 정규식용 모듈

정규식은 얘가 있어야 쓸 수 있다. 

elif keyword == "name":
    enzyme_RE = input("효소의 이름이 뭘로 시작하나요? ")
    enzyme_RE_2 = '^' + enzyme_RE

물론 if문에도 관련 코드를 추가했다. (^ 붙으면 그걸로 시작하는 걸 찾아준다)

else: 
    print("Enzyme with start with {0}".format(enzyme_RE))
    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 re.search(enzyme_RE_2,DB_enzyme):
            print("{0} | {1} | {2}".format(DB_enzyme,DB_seq,DB_site))
# 간단 검색(머릿글자)

물론 처리하는 코드도 추가했지... 

FASTA 파일 지원

1) 일단 사용자가 FASTA파일을 읽어오게 되면 따로 시퀀스 정보 입력을 받을 필요가 없어서 그거 관련도니 처리가 필요함. (시퀀스 이름은 Biopython SeqIO로 FASTA파일 ID 가져올 예정)

2) Biopython의 경우 FASTA파일에 꺾쇠가 하나면 read, 여러개면 parse로 가져오는데 parse로 읽을 걸 read로 가져오면 에러뜸. 에러에 대한 처리가 필요하고... parse로 가져오는 파일은 꺾쇠가 여러개인데, 이게 다르게 말하면 ID랑 시퀀스라 여러개라 그거 읽어오면 여러개를 돌리고 저장까지 해야 하는거라 이거는 Jupyter단에서는 힘듭니다... OTL 그래서 단식만 읽을거임... 

3) FASTA 관련된 기능이긴 한데, 본인은 애초에 경로가 고정되어 있으니까 상관 없는데... 그리고 파일이 저장되는 위치가 Jupyter 코드가 있는 경로기도 하고(홈), 리눅스는 경로가 간단간단해요 생각보다... usr 이런거 아니면. 그래서 경로가 고정되어 있지만, 사용자들 입장에서 이러면 상당히 불편하기때문에... FASTA 파일 열 때 뿐 아니라 파일 저장할 때 저장 경로를 선택할 수 있는 창을 띄울 예정입니다... 아니 누가 자기 파일 경로를 다 외워... 

Searcher 기능 관련

1) Wildcard 기능(시퀀스 혹은 효소 이름)

2) restriction site에 N이나 W같은 알파벳이 있을 때 처리법... 이건 finder보다 cutter단에서 처리하는 게 쉽습니다. finder는 찾아서 텍스트를 대체해야 하니... 이거는 둘 다 정규식 찾아서 직접 해 보고 적응해야 합니다. N은 차라리 wildcard니까 적용하기 쉽지... W M 이런거는 or이라;;;; 

사실 언제 추가될지는 모르겠지만 일단 우선적으로 해보겠음. 

기본

파이썬에는 두 가지 반복문이 있는데 한놈은 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는 안에 들어가는 글자가 많을 수록 멋지게 나옵니다. 이 점 유념하시길. 

일반적으로 우리가 쓰는 숫자는 10진수가 맞는데, 컴퓨터는 손가락이 두 개라 이진법을 쓴다. 예, 그겁니다. 일반적으로 10진법을 2진법으로 변환할때는 2로 나누는데, 13을 예시로 들자면 

1. 13/2=6...1
2. 6/2=3...0
3. 3/2=1...1

이렇게 구한 다음 3번 결과의 몫 1부터 시작해 3번 결과의 나머지-2번 결과의 나머지-1번 결과의 나머지 순으로 올라가서 1101이 된다. 1의 보수는 여기서 0을 1로, 1을 0으로 바꾼 0010. 2의 보수는 1의 보수에 1을 더하면 된다. (0011)

참고로 이진수로 변환된 숫자는 8, 16진수와 상호변환이 가능한데 

1. 110010010을 8진수로 바꾸려면 110 010 010으로 세자리씩 끊어서 10진법으로 바꾼다. (622)
2. 110010010을 16진수로 바꾸려면 1 1001 0010으로 네자리씩 끊어서 10진법으로 바꾼다. (192)
3. FF(16)를 8진수로 바꾸려면 2진수로 바꾼 다음 세자리씩 끊어서 10진법으로 바꾼다. (11 111 111->377)
4. 77(8)를 16진수로 바꾸려면 2진수로 바꾼 다음 네자리씩 끊어서 10진법으로 바꾼다. (11 1111->3F)

a = int(input())
bin_list = []
while a >= 1:
    bin_list.append(a%2)
    a = int(a/2)
bin = bin_list[::-1]
bin = ''.join(map(str,bin))
print(bin)
# 입력한 숫자를 이진수로 바꿔준다. (물론 수동 방식이다)

수동변환

for i in range(len(bin_list)):
    if bin_list[i] == 1:
        bin_list[i] = 0
    else:
        bin_list[i] = 1
bin = bin_list[::-1]
bin = ''.join(map(str,bin))
print(bin)
# 1의 보수(2의 보수는 이 방식으로 할 경우 처리가 되게 애매해지는 문제가 있다)

1의 보수

b=int(comp_1,2)
b=b+1
print(b)
# 그래서 2의 보수는 일단 10진수로 바꾼 다음, 1을 더하고 다시 바꿀 예정. 
comp2_list=[]
while b >= 1:
    comp2_list.append(b%2)
    b = int(b/2)
comp_2 = comp2_list[::-1]
comp_2 = ''.join(map(str,comp_2))
print(comp_2)
# 아, 0 빠진건 알아서 붙이세요.

2의 보수

a = int(input())
hex_list = []
while a >= 1:
    if a % 16 == 10: # 16진수에서 10~15까지는 각각 A~F로 표기한다. 즉, 여기에 대한 처리를 따로 진행해야 한다. 
        hex_list.append("A")
    elif a % 16 == 11:
        hex_list.append("B")
    elif a % 16 == 12:
        hex_list.append("C")
    elif a % 16 == 13:
        hex_list.append("D")
    elif a % 16 == 14: 
        hex_list.append("E")
    elif a % 16 == 15:
        hex_list.append("F")
    else: 
        hex_list.append(a % 16)
    a = int(a / 16)
hex = hex_list[::-1]
hex=''.join(map(str,hex))
print(hex)

자매품(16진수)

16진수의 경우 10~15까지를 A~F로 표기하기때문에 그걸 변환해줘야 한다. 

풀긴 풀었는데 이거 로직이 뭐가 문제인건지 제시한거랑 출력값이 다름. 

일단 문제가 뭐냐... 예를 들어서 아래와 같은 텍스트가 있다면 

츄라이츄라이

출력값은 

000123

이 된다. 네 번째 글자 츄를 기준으로 했을 때 

츄라이에서 츄/이츄/라이츄를 찾는건데 저기서 일치하는 게 츄 하나고 그게 한글자짜리거든. 

다섯번째 글자 라를 기준으로 하면 츄라이츄에서 라/츄라/이츄라/라이츄라 이렇게 찾게 되는거고 그렇게 되면 일치하는 문자열 중 제일 긴 게 츄라라서 2. 

이해하는 데 하루 걸렸다더니 실화냐고요? 네. 잠깐 뇌에 블루스크린 버프 와서요. 

text=input("Insert text: ")
# 전체 텍스트를 slicing하는 for문. 이 안에는 찾을 영역과 찾아야 할 영역이 포함되어 있다.
# 한 글자일때는 찾을 영역이 존재하지만 찾아야 하는 글자가 없으므로 첫 글자는 0이다.
find_list=[]
length=len(text)
for i in range(1,length+1):
    if len(text[0:i]) == 1:
        find_list.append(0)
    else:
        text_sub=text[0:i] # 전체 텍스트
        find=text_sub[0:i-1]# 찾을 범위
        max_list=[]
        for j in range(1,len(find)+1):
            subset=text_sub[j:] # 찾을 텍스트
            find_values=0
            if subset in find:
                max_list.append(len(subset))
            else:
                max_list.append(0)
            max_list=set(max_list)
            max_list=list(max_list)
            max_values=max(max_list) # 리스트에서 최대값을 추출
        find_list.append(max_values) # 했으면 넣어주세요
# 리스트 출력이 뭔가 이상하다.
# set을 적용하게 되면 각 iteration별로 고유값이 나오는 게 아니라 엉뚱한 값이 나오게 된다.
# append의 대상이 되는 리스트는 밖에 있지만, append는 안에 있어서 또 애매하고... append를 밖으로 뺄 수도 없다.
find_text=''.join(map(str,find_list))
# 리스트 형태로 처리한 것을 문자열로 붙여서 출력
print(find_text)

일단 전체 코드는 이거.

​

사실 문제가 정말 궁서체로 개빡세다. 어제 추워서 쉬긴 했지만 그래도 이틀 걸렸다. (주말에는 공부 안 함) 코딩은 컴퓨터 시켜먹으려면 일단 사람이 이해를 해야 하는데 이해하고 착수하는데 꽤 걸림+이게 내 맘대로 안돼서 걸림...

​

위에 설명한 게 문제인데, 그럼 필요한 기능이 뭐냐...

1. 리스트의 인덱싱/슬라이싱(각각 필요한 부분을 단식/범위로 찾는 기능)

2. For(반복문)

3. if(제어문)

4. 특정 문자를 찾는 기능(대충 R의 grep같은 거)

5. 입력(이건 맨 나중이라 우선순위 미뤄도 된다)

크게 이렇게 필요하다.

​

그럼 차근차근 게임코딩빨로 해봅시다. 이사람 대체 근데 솔직히 두뇌풀가동 이후로 엑스트라 밀고 겜코딩 켠 적 없음

Indexing/slicing 일반화 및 찾을 영역 정의하기 

일단 본인은 별찍기때도 그러했듯 small scale(고정값 혹은 실행하는 데 오래 안 걸리는 작은 범위)에서 해보고 일반화 과정을 거친다. 대충 수열로 치자면 어떤 규칙으로 변하는지 몇 개 찍어보고 일반항 공식 도출하는 방식이다. (별찍기 그렇게 한 다음 while로 찍었더니 while이 더 편하더라...)

어떤 n글자짜리 문자열이 있을 때 

-이 문자열의 인덱싱 주소는 0부터 (n-1)이고
-이 문자열의 슬라이싱 주소는 0부터 n까지이다. 이게 왜 그러냐... 아래 텍스트를 예시로 들어보자. 

>>> text="HAPPY DAYS!"
>>> len(text)
11
>>> print(text[0:11])
HAPPY DAYS!

H,A,P,P,Y, ,D,A,Y,S,!, 해서 11자(중간에 오타가 아니라 공백이 두 개 이고 이 때 인덱싱 주소는 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10이다. 근데 왜 11까지냐... 0에서 10까지니까 슬라이싱을 0:10으로 해 보면 

>>> print(text[0:10])
HAPPY DAYS

뭐야 느낌표 돌려줘요!!! 

이게 왜 이렇게 되는거냐면 파이썬 슬라이싱은 [부터:미만]으로 들어가기 때문이다. 그래서 0:10으로 하면 0부터 10 미만, 즉 9까지 뽑아준다. 

# 전체 텍스트를 slicing하는 for문. 이 안에는 찾을 영역과 찾아야 할 영역이 포함되어 있다. 
# 한 글자일때는 찾을 영역이 존재하지만 찾아야 하는 글자가 없으므로 첫 글자는 0이다. 
for i in range(1,length+1):
    if len(text[0:i]) == 1: 
        print(0)
    else:
        text_sub=text[0:i] # 전체 텍스트
        find=text_sub[0:i-1]# 찾을 범위 
        print(find) # 이게 찾을 범위

이게 찾는 범위를 지정하는 코드. 왜 1부터냐... 저게 [부터:미만]이라고 했는데 일단 첫 빠따는 뭐가 들어가던 0이다. 왜냐, 찾을 영역은 지정할 수 있는데 찾을 글자가 없다. 뭔 소리냐... 

[0,1,2,3,4,5]

위와 같은 리스트가 있을 때, 찾을 범위는 [0], [0,1],...,[0,4]까지이고 찾아야 할 글자는 범바범인데 범위가 [0]일 때 [1], [0,1]일 때 [2],[1,2] 이런 식으로 들어간다. 즉 찾는 영역은 앞에서부터 순차적으로 들어가고 찾아야 하는 영역은 뒤에서부터 순차적으로 들어간다. 

일차적인 for문에서 range가 1,len(text)+1인 이유도 그것때문이다. 저거 없이 그냥 박아버리면 0부터 시작하기때문에 나중에 음수 인덱싱을 하게 될 경우 문제가 생긴다. (-0이나 0이나 그게 그거라...) 그리고 length+1 해야 슬라이싱하면 length값까지 잡아준다. 

찾아야 할 부분집합 슬라이싱하기

# 전체 텍스트를 slicing하는 for문. 이 안에는 찾을 영역과 찾아야 할 영역이 포함되어 있다. 
# 한 글자일때는 찾을 영역이 존재하지만 찾아야 하는 글자가 없으므로 첫 글자는 0이다. 
for i in range(1,length+1):
    if len(text[0:i]) == 1: 
        print(0)
    else:
        text_sub=text[0:i] # 전체 텍스트
        find=text_sub[0:i-1]# 찾을 범위 
        subset_init=text_sub[i-1]
        print(subset_init) # 첫 서브셋이자 시발점

이 부분이 1차 난관이었다. 찾아야 하는 부분집합의 개수는 찾아야 할 영역의 크기에 따라 순차적으로 줄어들고, 찾아야 할 영역의 크기는 문자열의 길이에 따라 달라진다. 위에서 깜빡하고 안 썼는데, 이것은 마치 김밥 한 줄을 썰어서 왼쪽 꼬다리에서부터 오른쪽 꼬다리를 찾는 것과 비슷한 이치. 그러면 범위는 일단 김밥 한 줄로 잡아야 한다. 

이 때 전체 텍스트는 김밥 한 줄, 찾을 부분은 왼쪽 꼬다리라고 보면 된다. 찾는 부분집합은 찾는 영역에 따라 다르지만 김밥 한 줄 기준으로 오른쪽 꼬다리부터 하나씩 시작이다. (찾을 범위에서 오른쪽 꼬다리가 빠지고, 부분집합을 만드는 범위에서 왼쪽 꼬다리가 빠진다) 

그러니까 김밥 한 줄에서 오른쪽 꼬다리부터 순차적으로 지정해야 해서 원래 저 안에 for문이 또 들어간다. 그래서 전체 코드를 보면 for문이 두 개 들어간 것. ...저 시발점을 잘못 잡은 게 오른쪽 꼬다리면 무조건 -1번 인덱싱(맨 뒤)인데... 

제어문

# 전체 텍스트를 slicing하는 for문. 이 안에는 찾을 영역과 찾아야 할 영역이 포함되어 있다. 
# 한 글자일때는 찾을 영역이 존재하지만 찾아야 하는 글자가 없으므로 첫 글자는 0이다. 
find_subset=[]
max_subset=[]
for i in range(1,length+1):
    if len(text[0:i]) == 1: 
        find_subset.append(0)
        print(0)
    else:
        text_sub=text[0:i] # 전체 텍스트
        find=text_sub[0:i-1]# 찾을 범위 
        for j in range(1,len(find)+1):
            subset=text_sub[1:j+1] # 찾을 텍스트
            # 텍스트 유무에 따른 처리는 했는데, 문제는 텍스트가 있을 때 길이가 아니라 리스트 자체가 출력된다... 이거 어쩔겨... 
            if subset in find: 
                print(subset)
                find_subset.append(1)
            else: 
                find_subset.append(0)
print(find_subset)

[0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0]

(마른세수) 

출력이 개판인데 그거는 이따가 또 얘기합시다... 진짜 출력 잡는것도 대환장파티였어요... 

아무튼 저기 안쪽에 for문 하나 더 들어간 게 찾아야 할 부분집합이라는 건 알겠는데, if가 왜 있느냐... 일단 저 코드는 일치하는 게 있으면 1, 없으면 0을 끼워 넣으라는 의미긴 하다. 그러니까 쉽게 말하자면 일치하는 게 있냐 없냐를 먼저 보면서 범위부터 손볼 요량이었다. 

if subset in find: 
    print(subset)
    max_subset.append(1)
    find_subset.append(len(max_subset))
else: 
    find_subset.append(0)

제어문은 일단 이런 식. 해당하는 텍스트가 있을 때 일치하는 텍스트의 '길이'를 넣는다. 

if subset in find:
    max_list.append(len(subset))
else:
    max_list.append(0)
max_list=set(max_list)
max_list=list(max_list)
max_values=max(max_list) # 리스트에서 최대값을 추출

전체 코드에서는 이런 식으로 처리한다. 

출력 처리

여기서 정말 마른세수 여러번 했다... (마른세수)

# 전체 텍스트를 slicing하는 for문. 이 안에는 찾을 영역과 찾아야 할 영역이 포함되어 있다. 
# 한 글자일때는 찾을 영역이 존재하지만 찾아야 하는 글자가 없으므로 첫 글자는 0이다. 
find_subset=[] # 서브셋들의 최대값을 담아서 최종적으로 출력하는 리스트
max_subset=[] # 서브셋들의 최대값을 매기기 위한 리스트
for i in range(1,length+1):
    if len(text[0:i]) == 1: 
        find_subset.append(0)
    else:
        text_sub=text[0:i] # 전체 텍스트
        find=text_sub[0:i-1]# 찾을 범위 
        for j in range(1,len(find)):
            print(j,"th iteration")
            subset=text_sub[j:] # 찾을 텍스트
            print(find,len(find))
            print(subset,"*")
            if subset in find:
                max_subset.append(len(subset))
                max_values=max(max_subset)
                print(max_subset,max_values)
                find_subset.append(max(max_subset))
                print(find_subset)
            else:
                find_subset.append(0)
                print(find_subset)
# 리스트 출력이 뭔가 이상하다. 이거 중복값 처리를 어떻게 해야 하는거지?

이놈은 돌려봤더니 중복값이 고대로 나오고... 

# 전체 텍스트를 slicing하는 for문. 이 안에는 찾을 영역과 찾아야 할 영역이 포함되어 있다. 
# 한 글자일때는 찾을 영역이 존재하지만 찾아야 하는 글자가 없으므로 첫 글자는 0이다. 
find_list=[]
for i in range(1,length+1):
    if len(text[0:i]) == 1: 
        find_list.append(0)
    else:
        text_sub=text[0:i] # 전체 텍스트
        find=text_sub[0:i-1]# 찾을 범위 
        for j in range(1,len(find)):
            subset=text_sub[j:] # 찾을 텍스트
            find_values=0
            if subset in find:
                find_list.append(len(subset))
            else: 
                find_list.append(0)
            find_list=set(find_list)
            find_list=list(find_list)
# 리스트 출력이 뭔가 이상하다.
print(find_list)

이놈은 돌렸더니 중복값이 다 사라져서 글자수랑 안 맞고... 

일단 iteration의 기준이 뭐냐면, 찾을 범위를 하나 지정하고 거기서 문자열 차즌ㄴ 게 하나의 iteration이다. 즉, 첫번째 for문을 돌 때마다 값이 하나씩 저장되어야 한다. 

# 전체 텍스트를 slicing하는 for문. 이 안에는 찾을 영역과 찾아야 할 영역이 포함되어 있다. 
# 한 글자일때는 찾을 영역이 존재하지만 찾아야 하는 글자가 없으므로 첫 글자는 0이다. 
find_list=[]
for i in range(1,length+1):
    if len(text[0:i]) == 1: 
        find_list.append(0)
    else:
        text_sub=text[0:i] # 전체 텍스트
        find=text_sub[0:i-1]# 찾을 범위 
        max_list=[]
        for j in range(1,len(find)+1):
            subset=text_sub[j:] # 찾을 텍스트
            find_values=0
            if subset in find:
                max_list.append(len(subset))
            else: 
                max_list.append(0)
            max_list=set(max_list)
            max_list=list(max_list)
            max_values=max(max_list)
        find_list.append(max_values)
# 리스트 출력이 뭔가 이상하다.
# set을 적용하게 되면 각 iteration별로 고유값이 나오는 게 아니라 엉뚱한 값이 나오게 된다. 
# append의 대상이 되는 리스트는 밖에 있지만, append는 안에 있어서 또 애매하고... append를 밖으로 뺄 수도 없다. 
find_text=''.join(map(str,find_list))
print(find_text)

근데 이게 이거 뺀다고 해결되는거 실화냐... 아, 저기 join은 출력을 리스트 형태로 하기 위해 넣은거다. 

AGTC AGTC CGTG ATCT AGCT AGCT AGTC GCTG CATC AGTA C
0000 1234 1121 1121 1212 3456 7834 2232 2223 3452 1 # 돌려서 나온 결과
0000 1234 0000 1000 1200 1200 1234 0000 0100 1231 0 # 문제에서 제시한 결과

근데 이거 왜 결과가 이렇게 나오지...? 

역시나 회사명은 밝히지 않음. 

어제 면접보면서 코테 얘기가 나와서 블로그에 올려도 되겠느냐고 여쭤봤고, 그 쪽이 커리어에 도움이 될 테니 올려도 된다+올리면 안 되는 거면 따로 언질을 주겠다는 답변을 받았음. 

문제가 이거 말고 하나 더 있는데 그거는 바이오파이썬을 쓸 일이 없고(이것도 안썼지만...), 뇌에 블루스크린이 떠서 이해하느라 시간이 좀 걸려서 아마 구현은 빨라야 내일... 잘못하면 다음주까지도 갈 것 같음. 그래도 두번째껀 실행하는 시간은 빨라서 망정이지 얘는... 저 모듈 Jupyter에서 불러오지도 못함 OTL 

전체 코드

import vcf
import pandas as pd
# 모듈은 항상 위쪽에 부릅니다.
vcf_reader = vcf.Reader(open('/home/koreanraichu/clinvar_20211204.vcf', 'r'))
# gedit으로 여는데 엄청 방대해서 랙걸린다... (리눅스긴 한데 PC 사양이 그렇게 좋은 편은 아님)
Chromosome=[]
CLNSIG=[]
Gene=[]
# 일단 Series를 거쳐서 DataFrame화 할 예정.
for record in vcf_reader:
     INFO_get=record.INFO.get('CLNSIG')
     Gene_get=record.INFO.get('GENEINFO')
     if INFO_get:
          Chromosome.append(record.CHROM)
          CLNSIG.append(str(INFO_get))
     else:
          Chromosome.append(record.CHROM)
          CLNSIG.append("No data")
     # CLNSIG에 키가 없으면 처리하는거
     if Gene_get:
          Gene.append(Gene_get)
     else:
          Gene.append("No data")
     # Gene도 몇개 없더라...
# 이거 자체는 간단하다. 키가 없으면 NaN으로 넣는다는 얘기.
Chr_series=pd.Series(Chromosome)
CLN_series=pd.Series(CLNSIG)
Gene_series=pd.Series(Gene)
# 제발 시리즈 한번에 나오게 해주소서
Gene_df=pd.DataFrame({"Chromosome":Chr_series, "Gene":Gene_series,"CLNSIG":CLN_series})
# 오케이 제발 데이터프레임 이쁘게 나오게 해주소서
Gene_df2=Gene_df.groupby(['Chromosome','CLNSIG']).count()
# Groupby로 시원하게 묶어보자
# 근데 이거 object라 정렬이 1, 2, 3, 4로 안된다...
Gene_df3=Gene_df.groupby('Gene').count()
Gene_df3=Gene_df3.sort_values('CLNSIG',ascending=0)
print(Gene_df3.head(30))
# 정렬하고 30개 뽑아보자

문제

Q1. VCF파일 내 데이터 중 CLNSIC 데이터를 염색체별로 분류할 것
Q2. Pathogenic한 Gene Top 30 꼽기

공통 절차-분석을 위한 Dataframe 만들기

import vcf
vcf_reader = vcf.Reader(open('/home/koreanraichu/clinvar_20211204.vcf', 'r'))
for record in vcf_reader:
     print(record)
# PyVCF로 불렀는데 상당히 방대하다. 와... 이거 판다스 있어야된다...

일단 저걸 하려면 VCF 파일을 열어야 한다. VCF file의 경우 ClinVar에서 받을 수 있고, ftp에 일주일마다 업그레이드가 된다. 아무튼... 저거는 바이오파이썬으로 불러오는 게 아니라 PyVCF라는 모듈이 또 있어서 설치를 해야 하는데 여기서부터 난관이었다. 

1. 설치를 했는데 Jupyter에서 못 불러온다 
2. 파이참에서도 못 불러온다 
3. IDLE은 되네? 
4. 어? 갑자기 파이참에서 되네? (Jupyter는 안됨)

참고로 코드 실행 속도가 정말 어마무시하게 느리다. 데이터가 방대하니 어쩔 수 없지... 

import vcf
import pandas as pd
# 모듈은 항상 위쪽에 부릅니다.
vcf_reader = vcf.Reader(open('/home/koreanraichu/clinvar_20211204.vcf', 'r'))
for record in vcf_reader:
     print(record.INFO['CLNSIG'])
# PyVCF로 불렀는데 상당히 방대하다. 와... 이거 판다스 있어야된다...
     CLNSIG_frame=pd.DataFrame(record.INFO['CLNSIG'])

그리고 그냥 이렇게 Dataframe 만들려고 하면 Key error가 난다. 찾아보니 해당하는 키가 없으면 나는 에러라고 한다. 

import vcf
import pandas as pd
# 모듈은 항상 위쪽에 부릅니다.
vcf_reader = vcf.Reader(open('/home/koreanraichu/clinvar_20211204.vcf', 'r'))
# gedit으로 여는데 엄청 방대해서 랙걸린다... (리눅스긴 한데 PC 사양이 그렇게 좋은 편은 아님)
Chromosome=[]
CLNSIG=[]
# 일단 Series를 거쳐서 DataFrame화 할 예정.
for record in vcf_reader:
     INFO_get=record.INFO.get('CLNSIG')
     if INFO_get:
          Chromosome.append(record.CHROM)
          CLNSIG.append(INFO_get)
          print(record.CHROM,INFO_get)
     else:
          Chromosome.append(record.CHROM)
          CLNSIG.append("NaN")
          print(record.CHROM,"No keys")
# 이거 자체는 간단하다. CLNSIG에서 Key error가 나는데, 없으면 NaN으로 쓴다는 얘기.

그래서 이렇게 if함수를 넣어줍니다. Key가 있으면 해당하는 값을 넣고, 없으면 Nan(현재는 No data)을 넣으라는 얘기. 

Chr_series=pd.Series(Chromosome)
CLN_Series=pd.Series(CLNSIG)
# 심보러시여 제발 시리즈 한번에 나오게 해주소서

에이 솔직히 시리즈는 잘 돼요... 면접관님 보시는데 무슨 심보러 드립을 치고 있냐 

CLN_df=pd.DataFrame({"Chromosome":Chr_series, "CLNSIG":CLN_Series})
print(CLN_df)
# 오케이 제발 데이터프레임 이쁘게 나오게 해주소서

데이터프레임...... (마른세수) 

중간에 개고생하긴 했으나 해결. 

Chromosome                    CLNSIG
0                     1  [Uncertain_significance]
1                     1           [Likely_benign]
2                     1                  [Benign]
3                     1           [Likely_benign]
4                     1  [Uncertain_significance]
...                 ...                       ...
1105836              MT  [Uncertain_significance]
1105837              MT  [Uncertain_significance]
1105838              MT                  [Benign]
1105839              MT  [Uncertain_significance]
1105840  NW_009646201.1                  [Benign]

[1105841 rows x 2 columns]

Process finished with exit code 0

이렇게 데이터프레임이 나왔다. (참고로 최종 Dataframe은 2번 문제도 하나의 데이터프레임으로 풀기 위해 Gene 컬럼도 추가되어있음)

CLNSIG 집계

CLN_df2=CLN_df.groupby('Chromosome').count()
print(CLN_df2)

CLNSIG
Chromosome            
1                89545
10               41010
11               65665
12               49168
13               29113
14               33996
15               38488
16               60126
17               76248
18               18665
19               51709
2               112429
20               20230
21               12514
22               21575
3                56286
4                34468
5                58692
6                48407
7                51289
8                36992
9                50025
MT                2838
NW_009646201.1       1
X                46316
Y                   46

참고로 저게 Object라 정렬이 저렇게 된다... 예전에 scatter plot 그릴때도 저랬는데, 그때는 결측값 지우고 object에서 numeric으로 바꿔서 해결 봤다... 그건 축 값이 순수하게 숫자라 그런거고 저거는 염색체라 지우는 게 안돼요... (주륵) 

CLN_df2=CLN_df.groupby('Chromosome').aggregate(['CLNSIG'])
print(CLN_df2)
# 염색체별로 묶어보자.

이거 뭐 계속 바꿔서 해봤더니 아 리스트라 안된다고!!! 하더라... 그 답이 전체 코드에 있다. 

CLNSIG.append(str(INFO_get))

정확히는 이거. 저 INFO에 CLNSIG 말고 데이터가 많은데 그걸 다 리스트로 불러와서 저런 사태가 일어난 것이라고 한다. 딕셔너리 키값은 튜플이나 문자열같이 절대불변인 것들로 해야 하기 때문. 그래서 if문 돌리고 리스트에 넣을 때 문자열로 변환했다. 

Gene_df2=Gene_df.groupby(['CLNSIG','Chromosome']).count()
print(Gene_df2)
# Groupby(CLNSIG, Chromosome 순)

Gene
CLNSIG          Chromosome      
No data         1             56
                10            29
                11            57
                12            18
                13            10
...                          ...
['risk_factor'] 6             30
                7             16
                8             21
                9             20
                X             10

[527 rows x 1 columns]

니네 근데 진짜 이걸로 해결되는거 실화냐고... (No data: CLNSIG 데이터가 없음)

Gene_df2=Gene_df.groupby(['Chromosome','CLNSIG']).count()
print(Gene_df2)
# Groupby

Gene
Chromosome CLNSIG                                   
1          No data                                56
           ['Affects']                            19
           ['Benign', '_other']                    1
           ['Benign']                          15397
           ['Benign/Likely_benign', '_other']      2
...                                              ...
Y          ['Benign']                              4
           ['Likely_benign']                       5
           ['Likely_pathogenic']                   2
           ['Pathogenic']                         28
           ['Uncertain_significance']              7

[527 rows x 1 columns]

난 저거 말고 염색체별로 볼 건데? 하면 순서를 바꾸시면 됩니다. 

Top 30 pathogenic gene

이거 Top 30 뽑는거는 간단하다. ascending에 0 주고 head(30)으로 뽑으면 된다. 

Gene_df3=Gene_df.groupby('Gene').count()
Gene_df3=Gene_df3.sort_values('CLNSIG',ascending=0)
print(Gene_df3.head(30))
# 정렬하고 30개 뽑아보자

Chromosome  CLNSIG
Gene                                          
BRCA2:675                        13480   13480
BRCA1:672                        11565   11565
TTN:7273|TTN-AS1:100506866        9999    9999
APC:324                           8901    8901
NF1:4763                          7658    7658
TTN:7273                          7601    7601
TSC2:7249                         6392    6392
ATM:472                           6260    6260
MSH6:2956                         5591    5591
POLE:5426                         4939    4939
FBN1:2200                         4796    4796
RYR2:6262                         4358    4358
MSH2:4436                         4262    4262
RYR1:6261                         4063    4063
DMD:1756                          4063    4063
ATM:472|C11orf65:160140           3829    3829
PALB2:79728                       3788    3788
NEB:4703                          3706    3706
SYNE1:23345                       3486    3486
DICER1:23405                      3403    3403
MLH1:4292                         3354    3354
BRIP1:83990                       3350    3350
USH2A:7399                        3300    3300
PLEC:5339                         3267    3267
SMARCA4:6597                      3024    3024
PMS2:5395                         2949    2949
LDLR:3949                         2842    2842
TSC1:7248                         2776    2776
POLD1:5424                        2736    2736
CDH1:999                          2715    2715

이거는 유전자 전체 통틀어서 출력해주는거고 

Gene_df3=Gene_df.groupby(['Gene','CLNSIG']).count()
Gene_df3=Gene_df3.sort_values('Chromosome',ascending=0)
print(Gene_df3.head(30))
# 그냥 Gene으로 정렬하는 코드는 이거

Chromosome
Gene                       CLNSIG                                
BRCA2:675                  ['Uncertain_significance']        5467
APC:324                    ['Uncertain_significance']        4835
TTN:7273|TTN-AS1:100506866 ['Uncertain_significance']        3933
BRCA2:675                  ['Pathogenic']                    3590
TTN:7273                   ['Uncertain_significance']        3238
ATM:472                    ['Uncertain_significance']        3207
NF1:4763                   ['Uncertain_significance']        3101
POLE:5426                  ['Uncertain_significance']        3093
BRCA1:672                  ['Uncertain_significance']        3012
MSH6:2956                  ['Uncertain_significance']        2944
BRCA1:672                  ['Pathogenic']                    2892
TTN:7273|TTN-AS1:100506866 ['Likely_benign']                 2680
BRCA1:672                  ['not_provided']                  2654
BRCA2:675                  ['Likely_benign']                 2445
TTN:7273                   ['Likely_benign']                 2441
TSC2:7249                  ['Uncertain_significance']        2212
RYR2:6262                  ['Uncertain_significance']        2123
MSH2:4436                  ['Uncertain_significance']        2002
ATM:472|C11orf65:160140    ['Uncertain_significance']        1995
PALB2:79728                ['Uncertain_significance']        1956
BRIP1:83990                ['Uncertain_significance']        1917
NF1:4763                   ['Pathogenic']                    1895
RYR1:6261                  ['Uncertain_significance']        1757
DICER1:23405               ['Uncertain_significance']        1745
SYNE1:23345                ['Uncertain_significance']        1712
PLEC:5339                  ['Uncertain_significance']        1670
BRCA1:672                  ['Likely_benign']                 1617
NEB:4703                   ['Likely_benign']                 1617
NF1:4763                   ['Likely_benign']                 1596
RECQL4:9401                ['Uncertain_significance']        1577

Process finished with exit code 0

이렇게 하면 유전자별로 가장 높은 CLNSIG만 나온다. 

이놈들아 이것도 되면 좀 된다고 말좀 해줘... 

참고로 이거 어떻게 알았냐면 면접보는 회사에서 발표주제 중 하나가 저 두놈이었는데 찾다보니 NCBI에서 만든거네? -> Entrez에 있네? -> 비켜봐 시켜볼 게 있어(주섬주섬 파이참을 켠다) 가 된 거임. 

dbSNP

from Bio import Entrez
Entrez.email = "blackholekun@gmail.com" # 내가 누구인지 말해주는 과정이 필요하다고...
# 이메일은 자기꺼 그냥 쓰세요
handle = Entrez.esearch(db="snp", term="EGFR", retmax="40" )
record = Entrez.read(handle)
print(record)

참고로 db에는 snp라고 써야지 dbsnp라고 쓰면 안된다. 

Entrez.email = "blackholekun@gmail.com"
handle = Entrez.esearch(db="snp", term="rs1050171", retmax="40" )
record = Entrez.read(handle)
print(record['IdList'])
IdList=list(record['IdList'])
# dbSNP+esearch
for i in range(len(record['IdList'])):
    handle = Entrez.esummary(db="snp", id=IdList[i], retmode="xml")
    records = Entrez.read((handle))
    print(records)
# 일단 esummary 써서 다 가져왔다.

일단 워드클라우드 코드를 좀 응용해서 있는걸 다 털었는데 분량이 장난없어요... (마른세수) 이와중에 5000자 넘는다고 네이버 코드블록이 짤랐어 또... 

{'DocumentSummarySet': DictElement({'DocumentSummary': [DictElement({'SNP_ID': '1050171', 'ALLELE_ORIGIN': '', 'GLOBAL_MAFS': [{'STUDY': '1000Genomes', 'FREQ': 'A=0.432708/2167'}, {'STUDY': 'ALSPAC', 'FREQ': 'G=0.424235/1635'}, {'STUDY': 'Estonian', 'FREQ': 'G=0.435714/1952'}, {'STUDY': 'ExAC', 'FREQ': 'G=0.477223/57869'}, {'STUDY': 'FINRISK', 'FREQ': 'A=0.486842/148'}, {'STUDY': 'GnomAD', 'FREQ': 'G=0.483832/67753'}, {'STUDY': 'GnomAD_exomes', 'FREQ': 'G=0.476248/119743'}, {'STUDY': 'GoESP', 'FREQ': 'G=0.457558/5951'}, {'STUDY': 'GoNL', 'FREQ': 'G=0.380762/380'}, {'STUDY': 'HapMap', 'FREQ': 'A=0.405488/133'}, {'STUDY': 'KOREAN', 'FREQ': 'A=0.133219/389'}, {'STUDY': 'Korea1K', 'FREQ': 'A=0.126092/231'}, {'STUDY': 'MGP', 'FREQ': 'G=0.440075/235'}, {'STUDY': 'NorthernSweden', 'FREQ': 'G=0.483333/290'}, {'STUDY': 'PRJEB36033', 'FREQ': 'G=0.319149/30'}, {'STUDY': 'PRJEB37584', 'FREQ': 'A=0.153944/121'}, {'STUDY': 'PharmGKB', 'FREQ': 'G=0.344828/40'}, {'STUDY': 'Qatari', 'FREQ': 'G=0.439815/95'}, {'STUDY': 'SGDP_PRJ', 'FREQ': 'G=0.302083/116'}, {'STUDY': 'Siberian', 'FREQ': 'G=0.368421/14'}, {'STUDY': 'TOMMO', 'FREQ': 'A=0.159327/2670'}, {'STUDY': 'TOPMED', 'FREQ': 'G=0.482935/127828'}, {'STUDY': 'TWINSUK', 'FREQ': 'G=0.416667/1545'}, {'STUDY': 'Vietnamese', 'FREQ': 'A=0.285016/175'}, {'STUDY': 'ALFA', 'FREQ': 'G=0.427251/38873'}], 'GLOBAL_POPULATION': '', 'GLOBAL_SAMPLESIZE': '0', 'SUSPECTED': '', 'CLINICAL_SIGNIFICANCE': 'likely-benign,benign', 'GENES': [{'NAME': 'EGFR', 'GENE_ID': '1956'}, {'NAME': 'EGFR-AS1', 'GENE_ID': '100507500'}], 'ACC': 'NC_000007.14', 'CHR': '7', 'HANDLE': 'GENOMED,GRF,BUSHMAN,PJP,EVA_SAMSUNG_MC,SEATTLESEQ,CSHL,ILLUMINA,WUGSC_SSAHASNP,COMPLETE_GENOMICS,CSHL-HAPMAP,EGP_SNPS,KHV_HUMAN_GENOMES,EVA,EVA_FINRISK,BCMHGSC_JDW,EGCUT_WGS,JMKIDD_LAB,JJLAB,CGAP-GAI,CANCER-GENOME,EVA_UK10K_TWINSUK,ACPOP,TOMMO_GENOMICS,BGI,TISHKOFF,SGDP_PRJ,GNOMAD,SSMP,NHLBI-ESP,1000GENOMES,WEILL_CORNELL_DGM,ILLUMINA-UK,HUMANGENOME_JCVI,EVA-GONL,GMI,SI_EXO,EVA_EXAC,PACBIO,ENSEMBL,HAMMER_LAB,EVA_UK10K_ALSPAC,DDI,TOPMED,USC_VALOUEV,ABI,EVA_MGP,PERLEGEN,BIOINF_KMB_FNS_UNIBA,CLINSEQ_SNP,URBANLAB,KRGDB,PHARMGKB_PAAR-UCHI,HUMAN_LONGEVITY,LMM-PCPGM,LEE,OMUKHERJEE_ADBS,HGSV,CORNELL,SWEGEN,KOGIC,FSA-LAB,CPQ_GEN_INCA,EVA_DECODE', 'SPDI': 'NC_000007.14:55181369:G:A,NC_000007.14:55181369:G:C', 'FXN_CLASS': 'non_coding_transcript_variant,coding_sequence_variant,missense_variant,synonymous_variant,genic_downstream_transcript_variant', 'VALIDATED': 'by-frequency,by-alfa,by-cluster', 'DOCSUM': 'HGVS=NC_000007.14:g.55181370G>A,NC_000007.14:g.55181370G>C,NC_000007.13:g.55249063G>A,NC_000007.13:g.55249063G>C,NG_007726.3:g.167339G>A,NG_007726.3:g.167339G>C,NM_005228.5:c.2361G>A,NM_005228.5:c.2361G>C,NM_005228.4:c.2361G>A,NM_005228.4:c.2361G>C,NM_005228.3:c.2361G>A,NM_005228.3:c.2361G>C,NM_001346899.2:c.2226G>A,NM_001346899.2:c.2226G>C,NM_001346899.1:c.2226G>A,NM_001346899.1:c.2226G>C,NM_001346900.2:c.2202G>A,NM_001346900.2:c.2202G>C,NM_001346900.1:c.2202G>A,NM_001346900.1:c.2202G>C,NM_001346941.2:c.1560G>A,NM_001346941.2:c.1560G>C,NM_001346941.1:c.1560G>A,NM_001346941.1:c.1560G>C,NM_001346898.2:c.2361G>A,NM_001346898.2:c.2361G>C,NM_001346898.1:c.2361G>A,NM_001346898.1:c.2361G>C,NM_001346897.2:c.2226G>A,NM_001346897.2:c.2226G>C,NM_001346897.1:c.2226G>A,NM_001346897.1:c.2226G>C,NR_047551.1:n.1201C>T,NR_047551.1:n.1201C>G,NP_005219.2:p.Gln787His,NP_001333828.1:p.Gln742His,NP_001333829.1:p.Gln734His,NP_001333870.1:p.Gln520His,NP_001333827.1:p.Gln787His,NP_001333826.1:p.Gln742His|SEQ=[G/A/C]|LEN=1|GENE=EGFR:1956,EGFR-AS1:100507500', 'TAX_ID': '9606', 'ORIG_BUILD': '86', 'UPD_BUILD': '155', 'CREATEDATE': '2000/10/05 19:10', 'UPDATEDATE': '2021/04/26 11:38', 'SS': '1524562,4415417,14121776,16263880,17933589,23461029,24778961,43094183,52067157,65740167,66861655,74802379,77320037,85017386,86269798,93684223,98280473,104430036,105110075,112044216,116097662,142664329,143002862,159714800,159903775,162355178,166625159,203360709,212046487,223092668,233988791,240940108,279318753,285632487,293873711,342235809,479680970,482283259,485456160,490945938,491906494,534578333,560020602,654383955,779491194,781710098,834961318,836317040,974464379,984303175,1067488412,1074630855,1325217930,1431130856,1584052633,1593883812,1618262516,1661256549,1688745702,1711163743,1805015052,1927546725,1966658513,2024460166,2152655954,2294217561,2463237960,2634609883,2708323052,2736449362,2747825596,2853377601,3001152883,3023062922,3026026429,3347597637,3530770690,3530770691,3629823416,3632517137,3636855046,3646352982,3648637018,3669078603,3719737966,3734653562,3766593686,3785824290,3791125584,3796005564,3809753344,3824277165,3825524129,3825539990,3825720179,3830585782,3838783084,3844235500,3867317995,3914396600,3961507699,3984367631,3984367632,3984588810,3985298625,3986039617,3986382885,4746722701,5183240605,5236855147,5236861046,5237033464,5237196188', 'ALLELE': 'V', 'SNP_CLASS': 'snv', 'CHRPOS': '7:55181370', 'CHRPOS_PREV_ASSM': '7:55249063', 'TEXT': 'MergedRs=1050171', 'SNP_ID_SORT': '0001050171', 'CLINICAL_SORT': '1', 'CITED_SORT': '', 'CHRPOS_SORT': '0055181370', 'MERGED_SORT': '1'}, attributes={'uid': '58115520'})], 'DbBuild': 'Build211004-0955.1'}, attributes={'status': 'OK'})}

(대충 하나 뽑은게 이정도)

from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
handle = Entrez.esearch(db="snp", term="rs1050171", retmax="40" )
record = Entrez.read(handle)
IdList=list(record['IdList'])
# dbSNP+esearch
for i in range(len(record['IdList'])):
    handle = Entrez.esummary(db="snp", id=IdList[i], retmode="xml")
    records = Entrez.read((handle))
    print(records['DocumentSummarySet']['DocumentSummary'][0]['GLOBAL_MAFS'][0])
# 일단 esummary 써서 다 가져왔다. 근데 딕셔너리인데 왜 픽이 안될까...
# 그 전에 저렇게 꼭 4단픽까지 가야겠냐... 아니 진짜 너무한 거 아닙니까...

사실 5단이다. 딕셔너리가 아주 이중삼중이여... 

/home/koreanraichu/PycharmProjects/pythonProject/venv/bin/python "/home/koreanraichu/PycharmProjects/pythonProject/dbSNP in Entrez.py"
{'STUDY': '1000Genomes', 'FREQ': 'A=0.432708/2167'}
{'STUDY': '1000Genomes', 'FREQ': 'A=0.432708/2167'}
{'STUDY': '1000Genomes', 'FREQ': 'A=0.432708/2167'}
{'STUDY': '1000Genomes', 'FREQ': 'A=0.432708/2167'}
{'STUDY': '1000Genomes', 'FREQ': 'A=0.432708/2167'}
{'STUDY': '1000Genomes', 'FREQ': 'A=0.432708/2167'}

다 똑같아보이는데 제대로 가져온 거 맞냐고요? 네, 맞습니다. 

ClinVar

dbSNP는 그 뭐지... single nucleotide polymoorphism, 그러니까 point mutation으로 인한 변이들을 기록해 둔 데이터베이스고 ClinVar는 거기에 대한 보고서에 접근할 수 있게 해 주는 public archive다. 

from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
handle = Entrez.esearch(db="clinvar", term="EGFR", retmax="40" )
record = Entrez.read(handle)
print(record)
# clinVar+esearch

기본적으로는 이렇게 가져오고 

from Bio import Entrez
handle = Entrez.esearch(db="clinvar", term="L858R", retmax="40" )
record = Entrez.read(handle)
print(record['IdList'])
# clinVar+esearch
IdList=list(record['IdList'])
for i in range(len(record['IdList'])):
    handle = Entrez.esummary(db="clinvar", id=IdList[i], retmode="xml")
    records = Entrez.read((handle))
    print(records)
# 일단 esummary 써서 다 가져왔다.

일단 다 쓸어오자... (저 for문은 wordcloud 만들 때 썼던 걸 응용한거다)

{'DocumentSummarySet': DictElement({'DocumentSummary': [DictElement({'obj_type': 'single nucleotide variant', 'accession': 'VCV001314051', 'accession_version': 'VCV001314051.', 'title': 'NM_001243133.2(NLRP3):c.2744T>G (p.Leu915Arg)', 'variation_set': [{'measure_id': '1304312', 'variation_xrefs': [], 'variation_name': 'NM_001243133.2(NLRP3):c.2744T>G (p.Leu915Arg)', 'cdna_change': 'c.2744T>G', 'aliases': [], 'variation_loc': [{'status': 'current', 'assembly_name': 'GRCh38', 'chr': '1', 'band': '1q44', 'start': '247444052', 'stop': '247444052', 'inner_start': '', 'inner_stop': '', 'outer_start': '', 'outer_stop': '', 'display_start': '247444052', 'display_stop': '247444052', 'assembly_acc_ver': 'GCF_000001405.38', 'annotation_release': '', 'alt': '', 'ref': ''}, {'status': 'previous', 'assembly_name': 'GRCh37', 'chr': '1', 'band': '1q44', 'start': '247607354', 'stop': '247607354', 'inner_start': '', 'inner_stop': '', 'outer_start': '', 'outer_stop': '', 'display_start': '247607354', 'display_stop': '247607354', 'assembly_acc_ver': 'GCF_000001405.25', 'annotation_release': '', 'alt': '', 'ref': ''}], 'allele_freq_set': [], 'variant_type': 'single nucleotide variant', 'canonical_spdi': 'NC_000001.11:247444051:T:G'}], 'trait_set': [{'trait_xrefs': [{'db_source': 'MedGen', 'db_id': 'CN517202'}], 'trait_name': 'not provided'}], 'supporting_submissions': {'scv': ['SCV002002472'], 'rcv': ['RCV001771282']}, 'clinical_significance': {'description': 'Uncertain significance', 'last_evaluated': '2021/01/14 00:00', 'review_status': 'criteria provided, single submitter'}, 'record_status': '', 'gene_sort': 'NLRP3', 'chr_sort': '01', 'location_sort': '00000000000247444052', 'variation_set_name': '', 'variation_set_id': '', 'genes': [{'symbol': 'NLRP3', 'GeneID': '114548', 'strand': '+', 'source': 'submitted'}], 'protein_change': 'L801R, L858R, L915R, L917R', 'fda_recognized_database': ''}, attributes={'uid': '1314051'})], 'DbBuild': 'Build211201-1855.1'}, attributes={'status': 'OK'})}

(마른세수)

그래도 얘는 짤리진 않았음... 하지만 딕셔너리가 대체 몇중이지 하면서 쌍욕박는 건 똑같습니다... 

결국 차근차근 겜코딩 하던 머리로 차근차근 딕셔너리 마킹해가면서 찾음... (딕셔너리 다중일때는 ['상위 키']['하위 키'] 이런 식으로 찾으면 된다)

from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
handle = Entrez.esearch(db="clinvar", term="L858R", retmax="40" )
record = Entrez.read(handle)
# clinVar+esearch
IdList=list(record['IdList'])
for i in range(len(record['IdList'])):
    handle = Entrez.esummary(db="clinvar", id=IdList[i], retmode="xml")
    records = Entrez.read((handle))
    print(records['DocumentSummarySet']['DocumentSummary'][0]['title'])
# 일단 esummary 써서 다 가져왔다.
# 야이 근데 무슨 이건 진짜 할 말을 잃었음

NM_001243133.2(NLRP3):c.2744T>G (p.Leu915Arg)
NM_000245.4(MET):c.2573T>G (p.Leu858Arg)
NM_000222.3(KIT):c.2585T>G (p.Leu862Arg)
NM_005228.5(EGFR):c.2573_2574delinsGT (p.Leu858Arg)
NM_005228.5(EGFR):c.2572_2573inv (p.Leu858Arg)
NM_005228.5(EGFR):c.2369C>T (p.Thr790Met)
NM_005228.5(EGFR):c.2573T>G (p.Leu858Arg)

EGFR쪽에 보면 Leu858Arg 있는데 이게 L858R, Thr790Met이 T790M이다. (A는 알라닌)