(83)

번외편-코딩테스트 풀이 (3)

풀긴 풀었는데 이거 로직이 뭐가 문제인건지 제시한거랑 출력값이 다름. 일단 문제가 뭐냐... 예를 들어서 아래와 같은 텍스트가 있다면 츄라이츄라이 출력값은 000123 이 된다. 네 번째 글자 츄를 기준으로 했을 때 츄라이에서 츄/이츄/라이츄를 찾는건데 저기서 일치하는 게 츄 하나고 그게 한글자짜리거든. 다섯번째 글자 라를 기준으로 하면 츄라이츄에서 라/츄라/이츄라/라이츄라 이렇게 찾게 되는거고 그렇게 되면 일치하는 문자열 중 제일 긴 게 츄라라서 2. 이해하는 데 하루 걸렸다더니 실화냐고요? 네. 잠깐 뇌에 블루스크린 버프 와서요. text=input("Insert text: ") # 전체 텍스트를 slicing하는 for문. 이 안에는 찾을 영역과 찾아야 할 영역이 포함되어 있다. # 한 글자일때는 ..

번외편-코딩테스트 풀이 (2)

역시나 회사명은 밝히지 않음. 어제 면접보면서 코테 얘기가 나와서 블로그에 올려도 되겠느냐고 여쭤봤고, 그 쪽이 커리어에 도움이 될 테니 올려도 된다+올리면 안 되는 거면 따로 언질을 주겠다는 답변을 받았음. 문제가 이거 말고 하나 더 있는데 그거는 바이오파이썬을 쓸 일이 없고(이것도 안썼지만...), 뇌에 블루스크린이 떠서 이해하느라 시간이 좀 걸려서 아마 구현은 빨라야 내일... 잘못하면 다음주까지도 갈 것 같음. 그래도 두번째껀 실행하는 시간은 빨라서 망정이지 얘는... 저 모듈 Jupyter에서 불러오지도 못함 OTL 전체 코드 import vcf import pandas as pd # 모듈은 항상 위쪽에 부릅니다. vcf_reader = vcf.Reader(open('/home/koreanrai..

Biopython-dbSNP와 Clinvar

이놈들아 이것도 되면 좀 된다고 말좀 해줘... 참고로 이거 어떻게 알았냐면 면접보는 회사에서 발표주제 중 하나가 저 두놈이었는데 찾다보니 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라고 쓰면 안된다. ..

심심해서 써보는 본인 개발환경

Notes: 파이썬 모듈이나 버전은 따라해도 되는데 OS는 본인 편한거 쓰는게 좋음. 개발하시는 분들이 맥이나 리눅스 많이 쓰긴 한데 맥이나 리눅스는 윈도우와 달라서 무턱대고 쓰기 좀 불편하다. 나 고딩때 애들 커맨드 컨트롤 헷갈리는것만 여러번 봤음 진짜. (맥에서는 복사가 Ctrl+C가 아니라 Command+C) 그리고 어지간한 개발툴은 윈도 맥 리눅스 다 지원한다. 카톡 빼고 그건 개발툴 아닌데 아 그래서 지원 안하나 ​ 참고로 카톡이 리눅스에서 안돼서(와인 깔면 된다는데 실패함...) 본인 본의아니게 공부하는동안 연락 안됐음... 리눅스에서는 디코로만 연락 가능합니다. Laptop Lenovo thinkpad E470 (현재 단종) +4GB RAM ​ 내가 리눅스 쓸거다 근데 특히 우분투 쓸거다 그..

Biopython-Q&A

Q&A지만 자문자답이다. 어쨌든 질답은 맞음 Q1. MSA의 그 clustalW랑 MUSCLE은 어찌됐나요? A1. 그거 둘다 깔아야 됩니다. OS 박고 경로 박아서 돌리는 거 나오긴 했는데 트라이 해보려고 했더니 윈도 기준이네... 머슬은 경로 박아서 해봤는데 커맨드만 나와서 MSA는 터미널로 돌리고 있습니다. 리눅스는 clustalw와 MUSCLE 둘 다 일단 설치해두면 터미널에서 돌릴 수 있습니다. (py파일로 아웃풋 설정도 가능) Q2. 일부 건너뛴 챕터들이 있던데...? A2. 실습용 자료 구하기가 빡세거나(얘네도 다 안올려줌...) MSA처럼 할 수 없는 여건인 경우 건너뜁니다. 그래서 16 17(이건 하려고 했는데 wordcloud 하느라 시간 다 잡아먹음) 건너뛰고 케그 했지... Q3. ..

Biopython으로 KEGG 탐방하기

쿡북 분량 개짧음 진짜 이거보다 짧을수가 없음. KEGG? https://www.genome.jp/kegg/ KEGG: Kyoto Encyclopedia of Genes and Genomes www.genome.jp Kyoto Encyclopedia of Genes and Genomes의 준말. 그렇다, 이름에 교토가 들어간 걸 보면 아시겠지만 일제 DB다. 여기가 메인페이지 여기가 KEGG brite. 본인이 자주 가는 곳이다. 생각보다 쏠쏠한 정보가 많고 KEGG brite의 경우 golden standard dataset이라고 해서 야먀니시가 인공지능 학습시킬 때 쓴 데이터셋(GPCR, 효소, 핵 리셉터, 트랜스포터) 분류별로 약물 타겟을 알려준다. Parsing 파싱할거면 일단 파일이 당신 컴퓨터..

텍스트 입력받아서 Wordcloud 만들기

이것도 Project Wordcloud의 초기 코드이다. Text with wordcloud로, 이 당시에는 영문만 됐었다. (한글도 되긴 한데 접속사가 안떼짐) from wordcloud import WordCloud from wordcloud import STOPWORDS import matplotlib.pyplot as plot from PIL import Image import numpy as np # Summon module text = [] input_text= input("wordcloud로 만들 텍스트를 입력해주세요. ") text.append(input_text) yorn = input("더 추가할 텍스트가 있나요?") # 일단 입력을 받는다. (없으면 n, 있으면 n 말고 다른거) wh..

Biopython-Entrez에서 논문 제목 긁어와서 Wordcloud 만들기

Project wordcloud의 극초반 코드. from Bio import Entrez # 논문 긁어올 때 필요한거 from wordcloud import WordCloud from wordcloud import STOPWORDS import matplotlib.pyplot as plot # Wordcloud 그릴 때 필요한거 Entrez.email = "blackholekun@gmail.com" handle = Entrez.esearch(db="pubmed", term="Arabidopsis[title]", retmax="15") record = Entrez.read(handle) IdList=record['IdList'] # 일차적으로 Pubmed에서 논문을 찾을 수단인 PMID를 입수한다. # 날..

Biopython-Clustering 입력 인자

이거 튜토리얼에서 그냥 이런게 있다 하고 넘어가서 적음... Distance matrix Data : 분석할 데이터니까 당연히 있어야 한다. Mask(Default: None) : 결손된 데이터를 0으로 표시하는가? (None이면 모든 데이터가 존재하는 것) Weight(Default: None) : 가중치. 거리를 계산할 때 반영된다고 한다. Transpose(Default: 0) : True면 열간의 거리를, False면 행간의 거리를 계산한다. Dist(Default: 'e') : 거리 함수. 디폴트는 유클리드 거리. (Distance 관련해서는 clustering 이론편을 참조) Cluster distance data : 분석할 데이터 mask : 결손된 데이터를 0으로 표시하는가? (None이면..

Biopython으로 Clusting analysis 하기 (실전편)

와! 드디어 실전인가요? 근데 실전 생각보다 노잼임... Distance matrix 거리 행렬. 두 점간의 거리를 배열해 행렬로 나타낸 것이다. 점이 N개일 때 Distance matrix는 N*N으로 표기할 수 있다. import numpy as np import pandas as pd from Bio.Cluster import distancematrix data=np.array([[0, 1, 2, 3],[4, 5, 6, 7],[8, 9, 10, 11],[1, 2, 3, 4]]) matrix = distancematrix(data) # 뭐야 이거 왜 한영키 안먹어요 distances = distancematrix(data, dist='e') print(distances) [array([], dtype..

Biopython으로 Clusting analysis 하기 (이론편)

분량도 분량인데 spyder에서 한영키가 안돼서 그거땜시 늦었음... 이게 한 글에 다 쓰기엔 좀 분량도 분량인데 이게 생각보다 설명이랑 코딩이랑 나뉘어있어서 이론편 실전편 나눕니다. 이건 clustering 중 하나인 hierarchical clusting. 오늘 할 게 대충 이런거다. Cluster? 비슷한 특성을 가진 데이터 집단을 클러스터라고 한다. 데이터의 특성이 비슷하면 같은 클러스터, 다르면 다른 클러스터에 속한다. 클러스터링 하는 방법이 여러개가 있는데 여기서는 k-mean, k-median, k-medoid랑 hierarchical clustering에 대해 그냥 개 간단하게 설명하고 넘어간다. Hierarchical clustering 앞에 k-들어가는 것과 달리 계층적 클러스터링이라고..

Biopython으로 Sequence motif analysis 하기

모티브 찾으면 구조 모티브가 나오는데 이거는 '단백질이나 핵산과 같은 사슬 모양의 생체 분자에서 진화적으로 관련이 없는 다양한 분자에서 나타나는 일반적인 3차원 구조'로 정의한다. 근데 여기서 다루는 모티브는 그거 말고 시퀀스 모티브... 대충 이런거다. 위 그림은 뭔 시퀀스인지는 모르겠으나 3, 4, 5번째 염기가 GAA가 압도적으로 많은 듯. 오늘은 대충 from Bio import motifs 이런거 부른다. Motif object 객체 생성하기 from Bio import motifs from Bio.Seq import Seq instances=[Seq("TGTCGTATCG"),Seq("GTAAATAGCC"),Seq("GTAAATAACC"),Seq("TCGCGGAGCC"),Seq("ATGTGCCA..

Biopython으로 Phylogenetic tree 그리기

걍 메가 쓰세여... 메가가 편해... 짱이야 메가... 계통수? 계통수... 그러니까 Phylogenetic tree는 이런거다. (...) 유전자나 단백질 시퀀스 분석(균이라면 16s rRNA라던가)을 통해 얘네가 얼마나 가까운지를 알아내게 되면 그걸 저런 식으로 그려서 나타내는 것. 저렇게 생물 종에 따라 그리는 경우도 있고, 특정 단백질의 homolog나 다른 생물종에서 같은 역할을 하는 단백질에 대해서 저걸 그리기도 한다. 와! 계통수! 그려보자! from Bio import Phylo tree = Phylo.read("/home/koreanraichu/Deinococcus.ph", "newick") 이걸 쓰면 그릴 수 있는데... 어디가 일로와 끝까지 듣고 가... 저것만 쓰면 그려는 주는데 ..

Biopython으로 PDB 탐방하기

PDB!! PDB!! 들어가기 전에 PDB는 protein data bank의 약어이기도 하고, 거기서 제공하는 파일 형식이기도 하다. 여기서는 그냥 PDB라고 하면 데이터뱅크, PDB '파일'이라고 하면 PDB 파일이다. 그리고 이새기들 쿡북쓰기 귀찮았는지 모듈 불러와야 하는 거 빼먹더라... 니들도 일하기 싫었구나 파일 읽기 쓰기도 있긴 한데 그건 생략. 읽는것도 하난가 두갠가 오류나서 안됐다. 이 섹션에서 읽을 파일은 1) mmCIF 2) MMTF 3) PDB파일 4) PQR 이다. mmCIF from Bio.PDB.MMCIFParser import MMCIFParser parser = MMCIFParser() structure = parser.get_structure("7f0l", "/home/k..

Biopython으로 Swiss-prot과 ExPASy 데이터베이스 탐방하기

근데 스위스프롯은 긁어와서 저장 안됨? 첫빠따는 파싱이 국룰이지 파싱 방법이 네 가지가 있는데 gZIP은 생략. gZIP 파일을 못 구했다. handle=open('/home/koreanraichu/Q63HQ2.txt') print(handle) import gzip handle = gzip.open("myswissprotfile.dat.gz", "rt") from urllib.request import urlopen url = "https://raw.githubusercontent.com/biopython/biopython/master/Tests/SwissProt/F2CXE6.txt" handle = urlopen(url) print(handle) from Bio import ExPASy handle ..

번외편-코딩테스트 풀이 (3)

Coding/Python 2022. 8. 21. 02:22

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


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

츄라이츄라이

출력값은 

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 # 문제에서 제시한 결과

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

'Coding > Python' 카테고리의 다른 글

오케이 따옴표 떼버렸음  (0) 2022.08.21
10진수->2진수 변환 코드  (0) 2022.08.21
번외편-코딩테스트 풀이 (2)  (0) 2022.08.21
Biopython-dbSNP와 Clinvar  (0) 2022.08.21
심심해서 써보는 본인 개발환경  (0) 2022.08.21
Lv. 35 라이츄

Lv. 35 라이츄

광고 매크로 없는 청정한 블로그를 위해 노력중입니다. 근데 나만 노력하는 것 같음… ㅡㅡ

번외편-코딩테스트 풀이 (2)

Coding/Python 2022. 8. 21. 02:15

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

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

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

VCF file

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만 나온다. 

'Coding > Python' 카테고리의 다른 글

10진수->2진수 변환 코드  (0) 2022.08.21
번외편-코딩테스트 풀이 (3)  (0) 2022.08.21
Biopython-dbSNP와 Clinvar  (0) 2022.08.21
심심해서 써보는 본인 개발환경  (0) 2022.08.21
Biopython-Q&A  (0) 2022.08.21
Lv. 35 라이츄

Lv. 35 라이츄

광고 매크로 없는 청정한 블로그를 위해 노력중입니다. 근데 나만 노력하는 것 같음… ㅡㅡ

Biopython-dbSNP와 Clinvar

Coding/Python 2022. 8. 21. 02:09

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

참고로 이거 어떻게 알았냐면 면접보는 회사에서 발표주제 중 하나가 저 두놈이었는데 찾다보니 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는 알라닌)

'Coding > Python' 카테고리의 다른 글

번외편-코딩테스트 풀이 (3)  (0) 2022.08.21
번외편-코딩테스트 풀이 (2)  (0) 2022.08.21
심심해서 써보는 본인 개발환경  (0) 2022.08.21
Biopython-Q&A  (0) 2022.08.21
Biopython으로 KEGG 탐방하기  (0) 2022.08.21
Lv. 35 라이츄

Lv. 35 라이츄

광고 매크로 없는 청정한 블로그를 위해 노력중입니다. 근데 나만 노력하는 것 같음… ㅡㅡ

심심해서 써보는 본인 개발환경

Coding/Python 2022. 8. 21. 02:05

Notes: 파이썬 모듈이나 버전은 따라해도 되는데 OS는 본인 편한거 쓰는게 좋음. 개발하시는 분들이 맥이나 리눅스 많이 쓰긴 한데 맥이나 리눅스는 윈도우와 달라서 무턱대고 쓰기 좀 불편하다. 나 고딩때 애들 커맨드 컨트롤 헷갈리는것만 여러번 봤음 진짜. (맥에서는 복사가 Ctrl+C가 아니라 Command+C) 그리고 어지간한 개발툴은 윈도 맥 리눅스 다 지원한다. 카톡 빼고 그건 개발툴 아닌데 아 그래서 지원 안하나

참고로 카톡이 리눅스에서 안돼서(와인 깔면 된다는데 실패함...) 본인 본의아니게 공부하는동안 연락 안됐음... 리눅스에서는 디코로만 연락 가능합니다.


Laptop

Lenovo thinkpad E470 (현재 단종) +4GB RAM

내가 리눅스 쓸거다 근데 특히 우분투 쓸거다 그러면 우분투에서 여기다 우리꺼 깔면 잘됨 보장함 땅땅땅 있는 노트북 중에서 사는게 좋음.

일단 본인은 노트북 1호/3호와 달리 얘는 들고 다닐 거 감안해서 사양보다 휴대성에 중점을 뒀는데, 그러면 보통 그램을 삼. 1킬로도 안되는게 되게 가볍다더라고... 얘네 진짜 외계인 고문하나 근데 우분투 측에서 여기다 우리꺼 깔면 잘됨 땅땅땅 한 기종이 아니라 걸렀음... 예산 범위도 오버였고... 100만원 안으로 잡아놨는데 그램은 가격부터 오버여... 그램이 개발할 때 안 좋다는 얘기가 아니라 OS와 예산때문에 거른겁니다. 오해 ㄴㄴ요. 맥북도 예산때문에 걸렀음. 맥이 솔직히 편하긴 한게 고등학생때 써봤거든... 본인 놋북빼고 패드 폰 다 애플기기고. 근데 비싸요...

리눅스는 공짜기때문에 free-DOS 사서 직접 설치하면 됨. 부팅디스크 만드는 법도 찾아보면 나오고 설치법도 나오고... 뭐 우분투의 경우 쓰다보면 터미널에서 돌리는게 좀 많긴 한데 설치까지 터미널로 하는 건 아니고 윈도처럼 GUI로 돌아감. (리눅스도 종류가 많다... 리눅스 민트 센토스 뭐... 우분투는 그 중 하나.)

이 노트북은 원래 내장된 램이 4기가인데, 개발하다보면 생각보다 램이 많이 필요하다고 한다. 당장 노트북 1호가 게이밍이었던 이유도 MATLAB으로 영상 처리를 해야 하기 때문이고(영상처리 생각보다 고사양 필요함)... 그래서 4기가 엑스트라로 박았는데도 파이참 좀 버벅거림... (spyder나 jupyter는 제깍제깍 열린다) 저거 살 때는 파이썬 안했고 자바 했는데 이클립스로 돌아가야 하다 보니 저사양에서 버벅거리는 걸 몇 번 보기도 했고. 놋북 CPU는 i5고... 본인은 아무래도 잘알까지는 아니라서 CPU 사양을 크게 신경쓰는 건 아닌데 일단 i 뒤에 붙는 숫자는 세대가 같으면 큰게 좋고, 본인은 일단 셀러론 거름. 그래서 마지노선이 i3~i5단.

내가 셀러론에 안좋은 추억이 있어... 게임하다 랙걸려서 사리적립만 몇년 했는지 몰라...OTL 무슨 메이플 하는데 랙이 걸리냐고...


OS

Ubuntu 20.04 LTS(2022년 8월 21일 현재는 22.04 LTS)

커널은 올렸는데 LTS가 새로운게 안 풀린건지 내가 저장소를 따로 추가해야 하는 건지... 리누스형 이러지 말자 진짜

리눅스는 맥과 함께 유닉스 계열 OS고 악마의 명령어 sudo rm -rf/가 먹힌다 우분투의 경우 GUI가 있지만 터미널을 통해 CUI로 해야 하는 게 많다. 대표적인 예가 글꼴 뭉텅이로 설치할 때 복사할 글꼴이 있는 폴더로 가서 sudo cp *.ttf /usr/share/fonts 같은 거(cp=copy). 윈도우의 관리자 권한으로 실행같은 게 sudo인데 생각보다 이거때문에 터미널에서 해야 하는 게 좀 있다. 윈도보다 편한 부분도 있는데 윈도보다 불편한 부분도 존재함. 맥도 그렇잖아요? 대신 우분투 ~설치 이런걸로 찾아보면 정보는 많이 나오는 편.

그리고 얘 종특인지는 모르겠는데 한영키를 한영키가 아니라 오른쪽 알트로 인식해서 그거 셋업하고 한글 언어팩 깔아야 한글 입력 됨. 나는 우분투 깔면 제일 먼저 하는게 한글 깔고-mozc 설정하고(일어 자판)-무선랜 드라이버 추가하는거임. 그 뒤에 부차적으로 필요한 거 깔고 구축함.

대신 터미널에서 하는거에 적응하면 터미널 하나 켜놓고 거기서 다운로드 설치 실행 이동 복사 삭제가 다 가능하다. 본인 깃헙 세팅 터미널에서 하고 깃크라켄 까니까 알아서 셋업 인식하더만. 파이썬도 터미널에서 깔고 필요한 모듈도 일단은 터미널에서 깔았고 아나콘다 베이스도 다 터미널에서 돌리고(conda update 이런거) 파이참도 바로가기 없어서 터미널에서 켰음. (파이썬 IDLE은 바로가기가 있는데 왜 파이참만...)

아 그리고 바탕화면 그림은 바꿀 수 있는데 윈도처럼 바탕화면에 따라 작업표시줄 색 바꾸는건 안됨. 리누스씨 이 기능 빨리 넣어줘요

 


언어

Python 3.8

언어는 가급적 최신꺼 쓰는 편. (R도 깔려는 있음)


Module

-numpy

-sympy(+mpmath)

-scipy

-pandas

-matplotlib(부를때는 matplotlib.pyplot)

-Biopython

-sklearn

-tensorflow

-anaconda

(2022년 8월 21일 현재 request, Flask, pymongo 추가됨)

아 NetworkX도 깔았음. 그래프(이산수학 그래프... 막대그래프는 matplotlib이나 seaborn 쓰세요) 그려준대서...

아나콘다는 파이썬에서 소환하는 모듈은 아니고 뭐라고 해야되나... 저거 베이스로 깔아야 하는 모듈이나 IDE가 또 있음. sympy도 아나콘다 끼고 깔았던 것 같고... 심파이는 방정식 해주는건가 그렇고 scipy는... 뭐지 저거? 이봐요

넘파이는 배열 만드는거고(가끔 바이오파이썬 하다보면 마려운 그거 맞음) 판다스는 데이터프레임, 그러니까 표 만드는겁니다. csv나 엑셀파일 불러와서 표 만드는것도 가능하고 불러와서 계산하고 피벗테이블마냥 그룹별로 묶는것도 됨. 쓰다보면 진짜 편합니다. Biopython은 바이오파이썬이고 matplotlib은 그래프 그려주는거고... (seaborn이 좋다고 하는데 안써봐서 모름)

sklearn, tensorflow는 각각 머신러닝, 딥러닝 관련인데 이쪽은 뭐... 쓸 일이 있을 지 모르겠음. 사용법은 둘째치고 머신러닝이나 딥러닝이 어려워서 내장 데이터도 제대로 못 불렀음. ChEMBL에서 가져온 데이터 처리해서 써보려고 했더니 ChEMBL ID가 문자라서 인식을 못한다는데 ID도 빼줘야되냐고... OTL


IDE

-pycharm

-Jupyter(lab, notebook)

-spyder

-R studio(교육때문에 깔았지만 딱히 쓰지는 않음)

(2022년 8월 21일 현재 VScode도 설치해서 사용중)

spyder는 에리본씨 추천인데 한영키 안먹어서 뭐야 이거 왜이래요 했던거 드디어 해결봤고... ㄹㅇ 디자인 내 취향이었는데 거기서 막혀서 개고생했고... 바이오파이썬 할 때 고정적으로 쓰는건 파이참, 간단한건 IDLE(파이썬 깔면 IDLE도 따라옴), 간단한 코딩인데 IDLE에서 안될 것 같은건(예: 판다스 쿼리검색) jupyter로 돌린다. 저거 랩이랑 놋북 다 깔아야됩니다. (주파이터 랩 놋북 세개 깔아야됨)

리눅스는 바로가기고 뭐고 공통적으로 터미널에서 소환하면 알아서 열리는데(jupyter notebook 치면 실행됨) 파이참은 pycharm 치면 안되고 파이참 설치 경로 찾아가서 쉘 스크립트 실행해줘야 한다(설치도 다른 IDE랑 좀 다름). 즉 다른 IDE는 터미널 키고

spyder

jupyter lab

jupyter notebook

이렇게 치면 바로 열리는데 파이참은

cd /파이참 설치 경로/bin

sh pycharm.sh

해야 열린다. 그나마 탭키 누르면 자동완성 돼서 망정이지... (IDLE은 바로가기도 있고 진짜 간단한거 아니면 IDLE은 안씀) 그래서 편한거는 스파이더랑 주파이터가 편함. 쳐주면 알아서 실행도 되겠다... 근데 spyder는 설치를 pip가 아니라 anaconda 경유해서 해야 하고 한글 입력하려면 셋업을 또 따로 해줘야 한다. 코드는 영어 아닌가 주석도 영어로 다시게요? 파이참은 바로가기 아이콘이 없는 문제가 있지만 일단 설치하고 쓰는 데 크게 지장이 있는 건 아니다. 뭐야 아이콘 돌려줘요

셋업 면에서는 파이참 주파이터가 편하지만 실행 면에서는 스파이더 파이참이 편하다. 그리고 똑같은 다크테마여도 스파이더 다크테마가 내 취향임. (Jupyter는 다크모드 없는거같던디...) 주파이터는 랩보다 노트북을 더 많이 쓰는 편이다.

참고로 파이참은 터미널에서 인스톨하는 게 아니라 파이참 페이지에서 tar파일 받아서 압축풀고 거기 경로 들어가서 쉘 스크립트 실행하는 게 끝이다. 그래서 버전업 할거면 새 파이참 tar파일 받아서 압축 풀고 거기 가서 쉘 스크립트 돌리면 된다.

대부분 Jupyter 주피터로 읽는데 본인은 주파이터로 읽는다. python을 파이썬이라고 하잖음. 개발자는 뭐라고 읽나 궁금하네

여담이지만 처음에 Jupyter notebook이라고 하길래 인공지능 돌릴 용도로 노트북을 따로 준비해야 하나 했음... 나중에 교육 실습에 주파이터 나와서 어? 저거? 하다가 찾아봤더니 Jupyter notebook이래서 충격받았다. 그게 그거였냐고 근데 인공지능은 서버에서 돌아가는데 노트북을 별도로 준비해야됨? 아니 귀하신 분이라 별도로 모셔뒀나 했지

 


기타 툴

-ClustalW

-MUSCLE

-Gitkraken

-ATOM

(2022년 8월 21일 현재 MongoDB가 추가되었음)

노션은 바이오파이썬 Q&A에서도 언급했듯 리눅스에서는 설치가 아니라 웹노션을 쓴다. (윈도우와 맥 버전은 설치 가능) 그래서 쓰고 있지만 설치한 건 아니므로 생략. 코딩한 거 정리할 때는 노션이 짱이다. 코드블럭이 있어...

ClustalW와 MUSCLE은 MSA 툴인데 Biopython에서 돌리려면 깔아야 한대서 깔았더니 안돼서 걍 터미널로 돌리는 중. 역시 spyder나 jupyter 시리즈처럼 터미널에 clustalw, muscle 치면 실행된다. (muscle은 명령어 어떻게 치라고 나온다)

Gitkraken은 본진에도 설치한건데, 깃헙 관련된 건 다 이쪽에서 처리한다. 저장소에는 파일 이동만 하고 풀 커밋 푸시는 깃크라켄 통해서 하는 중. 이것도 okky에서 추천받은 툴이다. 리눅스에서 풀 커밋 푸시는 터미널로 해도 되는데 그거보단 얘가 편함. 이사람 깃헙 있음

ATOM에서 코드 쳐도 된다는데 본인은 오타 교정+자동완성이 되는 IDE가 편하다보니 그냥 글 쓸 때만 쓴다.

'Coding > Python' 카테고리의 다른 글

번외편-코딩테스트 풀이 (2)  (0) 2022.08.21
Biopython-dbSNP와 Clinvar  (0) 2022.08.21
Biopython-Q&A  (0) 2022.08.21
Biopython으로 KEGG 탐방하기  (0) 2022.08.21
텍스트 입력받아서 Wordcloud 만들기  (0) 2022.08.21
Lv. 35 라이츄

Lv. 35 라이츄

광고 매크로 없는 청정한 블로그를 위해 노력중입니다. 근데 나만 노력하는 것 같음… ㅡㅡ

Biopython-Q&A

Coding/Python 2022. 8. 21. 02:00

Q&A지만 자문자답이다. 어쨌든 질답은 맞음 


Q1. MSA의 그 clustalW랑 MUSCLE은 어찌됐나요? 

A1. 그거 둘다 깔아야 됩니다. OS 박고 경로 박아서 돌리는 거 나오긴 했는데 트라이 해보려고 했더니 윈도 기준이네... 머슬은 경로 박아서 해봤는데 커맨드만 나와서 MSA는 터미널로 돌리고 있습니다. 리눅스는 clustalw와 MUSCLE 둘 다 일단 설치해두면 터미널에서 돌릴 수 있습니다. (py파일로 아웃풋 설정도 가능)

 

Q2. 일부 건너뛴 챕터들이 있던데...? 

A2. 실습용 자료 구하기가 빡세거나(얘네도 다 안올려줌...) MSA처럼 할 수 없는 여건인 경우 건너뜁니다. 그래서 16 17(이건 하려고 했는데 wordcloud 하느라 시간 다 잡아먹음) 건너뛰고 케그 했지... 

 

Q3. 하면서 제일 재미없었을 때+제일 재밌었을 때는? 

A3. 실습할 거리가 많으면 많을수록 꿀잼입니다. clustering은 실습할 거리도 없는데다가 초장에 이론 나와서 노잼... 

 

Q4. 쿡북 보면서 할 때 제일 빡쳤던 부분은? 

A4. 처음에는 설명도 막 디테일하게 하고 이거이거이거 다돼! 이러더니 뒤로 갈수록 코드에 오타나고 모듈 빼먹고... 예제대로 따라하면 에러가 반기는... 심지어 그 에러 내가 찾아서 다 수정했어...... 쓰다가 귀찮았는지는 모르겠는데 귀차니즘이 거기서 터지면 안되지 이놈들아... 

 

Q5. 전체적인 개발 환경은? 

A5. 이건 따로 글 올려드리겠음. 여기다 풀면 분량이 길어져서;; 

 

Q6. 가끔 보다보면 응용편이 있던데 그건 대체 어떻게 떠올리신거죠 

A6. 화장실에서 X싸다가요. 이게 측간신 버프인가 그건가 그냥 코드 보면 떠오를 때도 있습니다. for문으로 돌리는거 while로 돌려본다던가... 

 

Q7. 저도 할수 있을까요? 

A7. 저도 파이썬 잘 못하는데 하잖음. 

 

Q8. IDE 배경색 뭐쓰세요?

A8. 원래 밝은 색 썼었는데 지금은 걍 깜장배경 씁니다. 개발자들은 깜장배경이 국룰인가... Jupyter는 흰 바탕입니다. 

 

Q9. 코드블럭에 입력하는 거 귀찮지 않으세요? 

A9. 노션은 코드블럭이나 있으니까 그나마 복붙하면 장땡이지 미디움은 코드블럭도 없어서 스크린샷 찍어서 붙입니다. 네이버도 초창기에 코드블럭 찾기 전에 쓴 건 스크린샷이었는데... 근데 결과 길면 복붙도 귀찮음 진짜... 

 

Q10. 코딩은 어디서 하세요?

A10. 방해받기 싫어서 스터디카페 가서 합니다. 군자역 열공다방 짱짱맨. 근데 밥 먹을데가 없음... 그리고 집중하다보면 끼니도 잊어버려서 음료수로 배 채우다가 집 가서 밥 먹습니다... 그래서 아침먹고 저녁먹고 야식먹고 근데 왜 살이 빠졌지...

 

Q11. 코딩할 때 중요한건 뭐라고 생각하십니까?

A11. 막 한번에 완성하려는 생각은 버리시고, 코드 입력하기 전에 대충 어떻게 돌아가면 좋을지 한번 그려보세요. 기능 구현하는 방법은 모르면 찾으면 되지만 그것도 무슨 기능을 어떤 식으로 구현할 지 알아야 찾아요. 아, 한가지 더. 큰 그림을 그릴 때는 일단 차근차근 해봅시다. 작은 단위(의 고정된 수치)로 먼저 해보고 점차 규모를 키워나가던가 일반화해야지 첫빠따에 일반화하려고 하면 에러가 여러분을 반깁니다. 내가 그래서 에러랑 많이 만났음 

 

Q12. 개발자한테 한마디?

A12. 문의메일 왜 안받냐. 심지어 내가 썬더버드로 보내서 메일 지메일인데. (리눅스 썬더버드에 지메일 연결함) 

 

Q13. 운영체제는 꼭 리눅스여야 하나요? 

A13. 윈도 편하면 그거 쓰세요. 리눅스는 tkinter인가 안됨. 저도 처음엔 그랬는데 이거 터미널 적응하는거랑 키세팅도 빡셉니다. 깔자마자 한글 입력 안돼서 그거 셋업해야되고 저는 놋북 무선랜 드라이버도 따로 넣어줘야되고 개귀찮음 진짜. 그리고 어지간한 개발툴은 윈도 맥 리눅스 지원하니까 편한거 쓰면 됩니다. 개발하시는 분들이 맥이나 리눅스 많이 쓰긴 하지만, 그냥 많이 쓴다고 무작정 그거 쓰는것보단 자기 편한거 쓰는 게 짱입니다. 디코도 지원하는 리눅스를 지원 안 하는 카톡은 뭘까 사랑해요 디코 

 

Q14. Notion 좋아요?

A14. 전직장 개발자님이 전파한건데 이거 ㄹㅇ 와 진짜 개편합니다. 달력 생성도 돼서 스케줄러에 메모장으로도 쓰고 있어요. 실행도 빠릿빠릿하고 (인터넷 연결은 필요하지만) 싱크도 빠르고... 패드에서 그림 올리면 PC에 금방 뜨고, 그림을 옮기거나 글을 쓰면 그것도 빠릿빠릿하게 나타나고... 리눅스에서는 설치형은 없고 웹노션이 있습니다. 노션도 지원하는 웹을 지원 안하는 카톡은 뭘까 

 

Q15. 워드프레스에는 따로 연재 안하시나요?

A15. 거기는 이미지도 HTML로 올려야돼서 귀차니즘이 배가됩니다. 간단한거(별찍기)는 이미지 많이 안 올려도 돼서 워드프레스에 올리지만... 

 

Q16. 책 추천해주세요!

A16. Biopython 책 산 건 있는데 한 번 보고 안읽은듯... 


Q17. 이거 재밌나요?

A17. 저는 재밌습니다. Wordcloud도 재밌었고... 근데 균 학명에 미역김치는 진짜 어디서 나온걸까... 정답: 미역김치 


Q18. 코드는 따로 올려두는 공간이 있나요? 여기에는 py파일이 없는 거 같은데...

A18. 깃헙에 올려두고 있습니다. 그거 덕분에 파이썬이 매트랩 분량 이겼음. 

Lv. 35 라이츄

Lv. 35 라이츄

광고 매크로 없는 청정한 블로그를 위해 노력중입니다. 근데 나만 노력하는 것 같음… ㅡㅡ

Biopython으로 KEGG 탐방하기

Coding/Python 2022. 8. 21. 01:54

쿡북 분량 개짧음 진짜 이거보다 짧을수가 없음. 


KEGG? 

https://www.genome.jp/kegg/

 

KEGG: Kyoto Encyclopedia of Genes and Genomes

 

www.genome.jp

Kyoto Encyclopedia of Genes and Genomes의 준말. 그렇다, 이름에 교토가 들어간 걸 보면 아시겠지만 일제 DB다. 

 

여기가 메인페이지

 

여기가 KEGG brite. 본인이 자주 가는 곳이다. 생각보다 쏠쏠한 정보가 많고 KEGG brite의 경우 golden standard dataset이라고 해서 야먀니시가 인공지능 학습시킬 때 쓴 데이터셋(GPCR, 효소, 핵 리셉터, 트랜스포터) 분류별로 약물 타겟을 알려준다. 


Parsing

파싱할거면 일단 파일이 당신 컴퓨터에 있어야 한다. 

ec_5.4.2.2.txt
0.21MB

급한대로 이거라도 쓰자. (아님 밑에 쿼링에서 갖고오든가...)

 

from Bio.KEGG import Enzyme
records = Enzyme.parse(open("/home/koreanraichu/ec_5.4.2.2.txt"))
record = list(records)[0]
print(record)

전체를 뽑고 싶으면 이렇게 하고 

 

from Bio.KEGG import Enzyme
records = Enzyme.parse(open("/home/koreanraichu/ec_5.4.2.2.txt"))
record = list(records)[0]
print(record.classname)
['Isomerases;', 'Intramolecular transferases;', 'Phosphotransferases (phosphomutases)']

특정 부분만 보고 싶으면 이렇게 한다. 근데 쟤 뭐 하는 효소냐... 

아, parse 친구 read도 있다. 데이터가 단식이면 리드, 뭉텅이면 파스. 

 

Querying

쿼리에 ing면 쿼링인가... 

 

Enzyme

from Bio.KEGG import REST
from Bio.KEGG import Enzyme
request = REST.kegg_get("ec:7.1.2.2")
# 참고: ATP synthase의 EC번호이다
open("ec_7.1.2.2.txt", "w").write(request.read())
records = Enzyme.parse(open("ec_7.1.2.2.txt"))
record = list(records)[0]
print(record)

참고: EC 7.1.2.2는 ATP synthase다 

 

파일이 생긴 것을 볼 수 있다. (본인은 Pycharm으로 해서 파이참 경로에 생겼다)

 

Pathway

사실상 KEGG의 근본은 pathway data다. 정말 세상 온갖가지 pathway가 다 있음. 진짜로. 

 

from Bio.KEGG import REST
human_pathways = REST.kegg_list("pathway", "hsa").read()
print(human_pathways)

뭔지 모르겠을 땐 일단 다 불러오는 게 국룰이다. 밑에 보여주겠지만 pathway만 써도 불러오는 건 된다. 

 

from Bio.KEGG import REST
human_pathways = REST.kegg_list("pathway", "hsa").read()
# 여기까지 분량이 정말 많은 것을 볼 수 있다.
repair_pathways = []
for line in human_pathways.rstrip().split("\n"):
    entry, description = line.split("\t")
    if "repair" in description:
        repair_pathways.append(entry)
print(repair_pathways)
# 조건부로 Repair pathway만 보는 코드
['path:hsa03410', 'path:hsa03420', 'path:hsa03430']

다 좋은데 그걸 꼭 hsa로 뽑아야겠냐... 알아서 찾아라 이거지 이럴거면 KEGG가지 뭐하러 여기서 찾냐 

 

from Bio.KEGG import REST
human_pathways = REST.kegg_list("pathway", "hsa").read()
# 여기까지 분량이 정말 많은 것을 볼 수 있다.
repair_pathways = []
for line in human_pathways.rstrip().split("\n"):
    entry, description = line.split("\t")
    if "repair" in description:
        repair_pathways.append(description)
print(repair_pathways)
# 조건부로 Repair pathway만 보는 코드
['Base excision repair - Homo sapiens (human)', 'Nucleotide excision repair - Homo sapiens (human)', 'Mismatch repair - Homo sapiens (human)']

아니 그럼 출력을 description으로 바꾸면 되잖아. 와 천잰데? 

 

from Bio.KEGG import REST
human_pathways = REST.kegg_list("pathway", "hsa").read()
# 여기까지 분량이 정말 많은 것을 볼 수 있다.
repair_pathways = []
for line in human_pathways.rstrip().split("\n"):
    entry, description = line.split("\t")
    if "*sis" in description:
        repair_pathways.append(description)
print(repair_pathways)
# 조건부로 Repair pathway만 보는 코드
[]

와일드카드는 안먹거나 *이 아닌 다른 기호인 듯 하다. 

 

Gene

Pathway에서 이어진다. 

 

from Bio.KEGG import REST
human_pathways = REST.kegg_list("pathway", "hsa").read()
# 여기까지 분량이 정말 많은 것을 볼 수 있다.
repair_pathways = []
for line in human_pathways.rstrip().split("\n"):
    entry, description = line.split("\t")
    if "repair" in description:
        repair_pathways.append(entry)
# 조건부로 Repair pathway만 보는 코드
# 여기까지는 파일에 있었으니까 다들 아시리라 믿고...
repair_genes = []
for pathway in repair_pathways:
    pathway_file = REST.kegg_get(pathway).read()  # query and read each pathway

    # iterate through each KEGG pathway file, keeping track of which section
    # of the file we're in, only read the gene in each pathway
    current_section = None
    for line in pathway_file.rstrip().split("\n"):
        section = line[:12].strip()  # section names are within 12 columns
        if not section == "":
            current_section = section

        if current_section == "GENE":
            gene_identifiers, gene_description = line[12:].split("; ")
            gene_id, gene_symbol = gene_identifiers.split()

            if not gene_symbol in repair_genes:
                repair_genes.append(gene_symbol)

print(
    "There are %d repair pathways and %d repair genes. The genes are:"
    % (len(repair_pathways), len(repair_genes))
)
print(", ".join(repair_genes))

일단 코드가 긴데 

1. Pathway 찾아주는 코드(그래서 entry로 쳐야 먹힌다)
2. 에서 Gene 섹션 찾아서 출력해준다. (hsa 번호가 맞는 pathway에서 gene을 찾아준다)

There are 3 repair pathways and 78 repair genes. The genes are:
OGG1, NTHL1, NEIL1, NEIL2, NEIL3, UNG, SMUG1, MUTYH, MPG, MBD4, TDG, APEX1, APEX2, POLB, POLL, HMGB1, XRCC1, PCNA, POLD1, POLD2, POLD3, POLD4, POLE, POLE2, POLE3, POLE4, LIG1, LIG3, PARP1, PARP2, PARP3, PARP4, FEN1, RBX1, CUL4B, CUL4A, DDB1, DDB2, XPC, RAD23B, RAD23A, CETN2, ERCC8, ERCC6, CDK7, MNAT1, CCNH, ERCC3, ERCC2, GTF2H5, GTF2H1, GTF2H2, GTF2H2C_2, GTF2H2C, GTF2H3, GTF2H4, ERCC5, BIVM-ERCC5, XPA, RPA1, RPA2, RPA3, RPA4, ERCC4, ERCC1, RFC1, RFC4, RFC2, RFC5, RFC3, SSBP1, PMS2, MLH1, MSH6, MSH2, MSH3, MLH3, EXO1

그래서 이렇게 나온다. 

근데 솔직히 인간적으로 저것만 찾으면 개노잼 아니냐. 그래서 준비했다. 

 

이건 산화적 인산화(Oxidative phosphorylation)인데, ATP를 얻기 위한 과정 중 하나이다. 여러분들이 일반적으로 알고 있는 ATP 생산 절차는

1. 우리가 먹은 포도당을 해당과정(Glycolysis)을 통해 피루브산 두 개를 만들고

2. TCA cycle(크렙 회로나 시트르산 사이클이라고도 한다) 들어가고(여기서부터 미토콘드리아 영역)

3. Oxidative phosphorylation을 통해 미토콘드리아가 ATP를 만든다(여기서 겁나 나온다)

이건데... Oxidative phosphorylation은 마지막 절차이다. 보통 TCA cycle부터 미토콘드리아 영역이고 피루브산이 TCA cycle로 가면 돌아올 수 없는 강을 건너게 된다. (피루브산 상태에서 포도당신생합성(Gluconeogenesis)을 돌리기도 한다. 번역명 왜이래요. 번역 간지나네)

여담이지만 과당도 해당과정 중간에 낄 수 있다. (Galactose도 낀다...젖당이었나 갈락이었나 아무튼 낌)

 

from Bio.KEGG import REST
human_pathways = REST.kegg_list("pathway", "hsa").read()
# 여기까지 분량이 정말 많은 것을 볼 수 있다.
Oxidative_pathways = []
Oxidative_pathways_name =[]
for line in human_pathways.rstrip().split("\n"):
    entry, description = line.split("\t")
    if "Oxidative" in description:
        Oxidative_pathways.append(entry)
        Oxidative_pathways_name.append(description)
print("Entry no: ",Oxidative_pathways,"description is: ",Oxidative_pathways_name)
# 조건부로 Repair pathway만 보는 코드
# 여기까지는 파일에 있었으니까 다들 아시리라 믿고...
Oxidative_genes = []
for pathway in Oxidative_pathways:
    pathway_file = REST.kegg_get(pathway).read()  # query and read each pathway

    # iterate through each KEGG pathway file, keeping track of which section
    # of the file we're in, only read the gene in each pathway
    current_section = None
    for line in pathway_file.rstrip().split("\n"):
        section = line[:12].strip()  # section names are within 12 columns
        if not section == "":
            current_section = section

        if current_section == "GENE":
            gene_identifiers, gene_description = line[12:].split("; ")
            gene_id, gene_symbol = gene_identifiers.split()

            if not gene_symbol in Oxidative_genes:
                Oxidative_genes.append(gene_symbol)

print(
    "There are %d pathways and %d genes. The genes are:"
    % (len(Oxidative_pathways), len(Oxidative_genes))
)
print(", ".join(Oxidative_genes))
Entry no:  ['path:hsa00190'] description is:  ['Oxidative phosphorylation - Homo sapiens (human)']
There are 1 pathways and 134 genes. The genes are:
ND1, ND2, ND3, ND4, ND4L, ND5, ND6, NDUFS1, NDUFS2, NDUFS3, NDUFS4, NDUFS5, NDUFS6, NDUFS7, NDUFS8, NDUFV1, NDUFV2, NDUFV3, NDUFA1, NDUFA2, NDUFA3, NDUFA4, NDUFA4L2, NDUFA5, NDUFA6, NDUFA7, NDUFA8, NDUFA9, NDUFA10, NDUFAB1, NDUFA11, NDUFA12, NDUFA13, NDUFB1, NDUFB2, NDUFB3, NDUFB4, NDUFB5, NDUFB6, NDUFB7, NDUFB8, NDUFB9, NDUFB10, NDUFB11, NDUFC1, NDUFC2, NDUFC2-KCTD14, SDHA, SDHB, SDHC, SDHD, UQCRFS1, CYTB, CYC1, UQCRC1, UQCRC2, UQCRH, UQCRHL, UQCRB, UQCRQ, UQCR10, UQCR11, COX10, COX3, COX1, COX2, COX4I2, COX4I1, COX5A, COX5B, COX6A1, COX6A2, COX6B1, COX6B2, COX6C, COX7A1, COX7A2, COX7A2L, COX7B, COX7B2, COX7C, COX8C, COX8A, COX11, COX15, COX17, CYCS, ATP5F1A, ATP5F1B, ATP5F1C, ATP5F1D, ATP5F1E, ATP5PO, ATP6, ATP5PB, ATP5MC1, ATP5MC2, ATP5MC3, ATP5PD, ATP5ME, ATP5MF, ATP5MG, ATP5PF, ATP8, ATP6V1A, ATP6V1B1, ATP6V1B2, ATP6V1C2, ATP6V1C1, ATP6V1D, ATP6V1E2, ATP6V1E1, ATP6V1F, ATP6V1G1, ATP6V1G3, ATP6V1G2, ATP6V1H, TCIRG1, ATP6V0A2, ATP6V0A4, ATP6V0A1, ATP6V0C, ATP6V0B, ATP6V0D1, ATP6V0D2, ATP6V0E1, ATP6V0E2, ATP6AP1, ATP4A, ATP4B, ATP12A, PPA2, PPA1, LHPP

물량 쩌네. 


Appendix. 뭐 딴거 없어요? 메인페이지에 개많던데? 

나도 궁금해서 KEGG 메인에 있는거 다 돌려봤음. 

 

from Bio.KEGG import REST
human_pathways = REST.kegg_list('brite').read()
print(human_pathways)
from Bio.KEGG import REST
human_pathways = REST.kegg_list('pathway').read()
print(human_pathways)
from Bio.KEGG import REST
human_pathways = REST.kegg_list('module').read()
print(human_pathways)
from Bio.KEGG import REST
human_pathways = REST.kegg_list('enzyme').read()
print(human_pathways)
from Bio.KEGG import REST
human_pathways = REST.kegg_list('genome').read()
print(human_pathways)
from Bio.KEGG import REST
human_pathways = REST.kegg_list('glycan').read()
print(human_pathways)
from Bio.KEGG import REST
human_pathways = REST.kegg_list('compound').read()
print(human_pathways)
from Bio.KEGG import REST
human_pathways = REST.kegg_list('reaction').read()
print(human_pathways)
from Bio.KEGG import REST
human_pathways = REST.kegg_list('disease').read()
print(human_pathways)
from Bio.KEGG import REST
human_pathways = REST.kegg_list('drug').read()
print(human_pathways)
from Bio.KEGG import REST
human_pathways = REST.kegg_list('network').read()
print(human_pathways)

저기서 되는거 Genes빼고 어지간한건 다 된다. 

1. KEGG Brite
2. KEGG Pathway
3. KEGG Module
4. KEGG Enzyme
5. KEGG Compound
6. KEGG Genome Gene은 안되는데 Genome은 되네? 
7. KEGG Glycan 
8. KEGG Reaction
9. KEGG Network
10. KEGG Drug
11. KEGG Disease


Appendix 2. 응용편

from Bio.KEGG import REST
Disease = REST.kegg_list('disease').read()
# 가져올 수 있는 것: brite, pathway, genome(gene은 안됨), module, enzyme, glycan, compound, reaction, network, drug, disease
disease_leukemia = []
for line in Disease.rstrip().split("\n"):
    entry, description = line.split("\t")
    if "leukemia" in description:
        disease_leukemia.append(description)
print(disease_leukemia)

다 불러올 수 있으면 아까 그 쿼링 코드 써서 검색도 되나 궁금해서 찾아봤음. 

 

from Bio.KEGG import REST
drug = REST.kegg_list('drug').read()
# 가져올 수 있는 것: brite, pathway, genome(gene은 안됨), module, enzyme, glycan, compound, reaction, network, drug, disease
drug_list = []
drug_name=input("찾고 싶은 약과 관련된 것을 입력해주세요: ")
for line in drug.rstrip().split("\n"):
    entry, description = line.split("\t")
    if drug_name in description:
        drug_list.append(description)
print(drug_list)

이건 나름 업그레이드 버전인데 Drug에서 사용자가 입력한 항목을 찾아준다. (와일드카드는 안 먹는다)

Lv. 35 라이츄

Lv. 35 라이츄

광고 매크로 없는 청정한 블로그를 위해 노력중입니다. 근데 나만 노력하는 것 같음… ㅡㅡ

텍스트 입력받아서 Wordcloud 만들기

Coding/Python 2022. 8. 21. 01:45

이것도 Project Wordcloud의 초기 코드이다. Text with wordcloud로, 이 당시에는 영문만 됐었다. (한글도 되긴 한데 접속사가 안떼짐)


from wordcloud import WordCloud
from wordcloud import STOPWORDS
import matplotlib.pyplot as plot
from PIL import Image
import numpy as np  
# Summon module
text = []
input_text= input("wordcloud로 만들 텍스트를 입력해주세요. ")
text.append(input_text)  
yorn = input("더 추가할 텍스트가 있나요?")
# 일단 입력을 받는다. (없으면 n, 있으면 n 말고 다른거)
while yorn != "n":
    yorn = input("더 추가할 텍스트가 있나요? ")
    if yorn == "n":
        text = ''.join(text)
    else:
        input_text= input("wordcloud로 만들 텍스트를 입력해주세요 ")
        text.append(input_text)
# 추가로 입력이 있을 경우 입력이 '없다'고 할 떄까지 입력을 받고 입력이 더 이상 들어오지 않으면 join한다. (입력 받을때마다 join하면 에러남)
image = np.array(Image.open("/home/koreanraichu/600px-793Nihilego.png"))
# 이미지를 부르고
font_path = '/usr/share/fonts/Maplestory Light.ttf'
wordcloud = WordCloud(font_path = font_path,background_color="#ffffff",colormap="PuBu_r",width = 800, height=800, mask=image)
wordcloud = wordcloud.generate_from_text(text)
# 그리고
plot.figure(figsize=(15,15))
plot.axis('off')
plot.imshow(wordcloud)
plot.show()
# 출력

전체 코드(n or N으로 줘봤는데 안먹힘...)


결과(릴리요 도감 설명) 

Lv. 35 라이츄

Lv. 35 라이츄

광고 매크로 없는 청정한 블로그를 위해 노력중입니다. 근데 나만 노력하는 것 같음… ㅡㅡ

Biopython-Entrez에서 논문 제목 긁어와서 Wordcloud 만들기

Coding/Python 2022. 8. 21. 01:43

Project wordcloud의 극초반 코드. 


from Bio import Entrez
# 논문 긁어올 때 필요한거
from wordcloud import WordCloud
from wordcloud import STOPWORDS
import matplotlib.pyplot as plot
# Wordcloud 그릴 때 필요한거
Entrez.email = "blackholekun@gmail.com"
handle = Entrez.esearch(db="pubmed", term="Arabidopsis[title]", retmax="15")
record = Entrez.read(handle)
IdList=record['IdList']
# 일차적으로 Pubmed에서 논문을 찾을 수단인 PMID를 입수한다.
# 날짜검색이 안돼서 키워드에 '애기장대'가 있는 논문 10개를 찾음.
for i in range(len(IdList)):
    handle = Entrez.esummary(db="pubmed", id=IdList[i], retmode="xml")
    records = Entrez.parse(handle)
    title=[]
    for record in records:
        title.append(record['Title'])
# 타이틀 뽑았다.
# wordcloud에서 뺄 단어들(접속사 이런거)
font_path = '/usr/share/fonts/truetype/nanum/BinggraeMelona.ttf'
wordcloud = WordCloud(stopwords=STOPWORDS,   font_path = font_path, width = 800,
    height = 800,background_color="#ffffff")
wordcloud = wordcloud.generate_from_text(record['Title'])
plot.axis('off')
plot.imshow(wordcloud)
plot.show()
# wordcloud generate

코드 구성은 크게 

1. 모듈(wordcloud용 모듈과 Entrez 접속용 모듈)
2. Entrez에서 긁어오기(PMID-title)
3. Wordcloud 그리기


로 나누어진다. 


주석에는 10개로 되어 있었는데, 이게 40개를 가져오다가 서버 에러가 나서 규모를 좀 줄였다. (현재 15개) 코드를 일반화함에 따라 중간에 for문도 바꼈다. 

워드클라우드의 경우 Stopword라고 해서 접속사나 감탄사같은 것들을 다 빼야 하는데(in, at, the같은 거) 영문의 경우 STOPWORDS에 저장되어 있다. 

1차 결과물(10개)

 

2차 결과물(15개, 폰트 변경)

근데 내가 원하는 결과는 대충 위에꺼인 거 같다... 


from Bio import Entrez
# 논문 긁어올 때 필요한거
from wordcloud import WordCloud
from wordcloud import STOPWORDS
import matplotlib.pyplot as plot
# Wordcloud 그릴 때 필요한거
import numpy as np
from PIL import Image
# 이미지 마스킹할 때 필요한거
image=np.array(Image.open("/home/koreanraichu/alice.png"))
font_path = '/usr/share/fonts/dovemayo.otf'
# 마스킹할 이미지
Entrez.email = "blackholekun@gmail.com"
handle = Entrez.esearch(db="pubmed", term="EGFR[title]", retmax="20")
record = Entrez.read(handle)
IdList=record['IdList']
# 일차적으로 Pubmed에서 논문을 찾을 수단인 PMID를 입수한다.
# 숫자를 20개로 늘려서 시간이 좀 걸린다.
title=[]
for i in range(len(IdList)):
    handle = Entrez.esummary(db="pubmed", id=IdList[i], retmode="xml")
    records = Entrez.parse(handle)
    for record in records:
        title.append(record['Title'])
title=''.join(title)
# 타이틀 뽑았다.
# 아 이거 근데 묶어줘야되더라...
wordcloud = WordCloud(stopwords=STOPWORDS,font_path = font_path, width = 800,height = 800,
                      background_color="#ffffff",colormap='spring',mask=image)
wordcloud = wordcloud.generate_from_text(record['Title'])
plot.axis('off')
plot.imshow(wordcloud)
plot.show()
# wordcloud generate

PS. 마스크 기능 추가하는 김에 확인해봤는데 하나밖에 안되길래 리스트업한거 join해버림. 

PS2. +코드의 wordcloud = wordcloud.generate_from_text(record['Title'])부분을 wordcloud = wordcloud.generate_from_text(Title)로 수정해줘야 제대로 나온다. 

 

뭐야 왜이렇게 작아요 

 

하는김에 김치로 찾음. 

 

Arabidopsis+2021[Years]로 준 것

Lv. 35 라이츄

Lv. 35 라이츄

광고 매크로 없는 청정한 블로그를 위해 노력중입니다. 근데 나만 노력하는 것 같음… ㅡㅡ

Biopython-Clustering 입력 인자

Coding/Python 2022. 8. 21. 01:40

이거 튜토리얼에서 그냥 이런게 있다 하고 넘어가서 적음... 


Distance matrix

Data : 분석할 데이터니까 당연히 있어야 한다.
Mask(Default: None) : 결손된 데이터를 0으로 표시하는가? (None이면 모든 데이터가 존재하는 것)
Weight(Default: None) : 가중치. 거리를 계산할 때 반영된다고 한다.
Transpose(Default: 0) : True면 열간의 거리를, False면 행간의 거리를 계산한다.
Dist(Default: 'e') : 거리 함수. 디폴트는 유클리드 거리. (Distance 관련해서는 clustering 이론편을 참조)

 

Cluster distance

data : 분석할 데이터
mask : 결손된 데이터를 0으로 표시하는가? (None이면 모든 데이터가 존재하는 것)
weight : 가중치. 거리를 계산할 때 반영된다고 한다.
index1, index2 : 각각 첫번째, 두번째 클러스터에 포함되는 인덱스'들'. 리스트 혹은 정수 형태이다.
method : Cluster간 거리를 정의하는 방법. 
'a': 산술 평균
'm': 중앙값
's': 두 클러스터의 항목 중 가장 짧은 pairwise distance
'x': 두 클러스터의 항목 중 가장 긴 pairwise distance
'v': 두 클러스터의 항목간 평균 pairwise distance
dist : 위에 있는 거리 함수. 디폴트는 유클리드 거리.  관련해서는 위쪽을 참고할 것)
transpose : True면 열간의 거리를, False면 행간의 거리를 계산한다. (True면 컬럼)

 

평균/중앙값 클러스터링

data : 분석할 데이터(필수)
nclusters : 클러스터 개수(기본: 2)
mask : 결손된 데이터를 0으로 표시하는가? (None이면 모든 데이터가 존재하는 것, 기본: 없음)
weight : 가중치(기본: 없음)
transpose : True면 열간의 거리를, False면 행간의 거리를 계산한다. (기본: F)
npass : 알고리즘 가동 횟수(기본: 1)
method : k-mean이면 'a', k-median이면 'm'. 기본값은 a다.
dist : 거리(기본: e)
initialid : initial clustering을 정의한다. (기본: 없음)

 

medoid 클러스터링

distance : 얘는 특이하게도 어레이 자체가 아니라 거리 행렬이 필요하다. 
nclusters : 클러스터 개수(기본: 2)
npass : 알고리즘 가동 횟수(기본: 1)
initialid : initial clustering을 정의한다. (기본: 없음)

 

self organizing maps 

data : 데이터 
mask : 결손된 데이터를 0으로 표시하는가? (None이면 모든 데이터가 존재하는 것, 기본: 없음)
weight : 가중치(기본: 없음)
transpose : True면 열간의 거리를, False면 행간의 거리를 계산한다. (기본: F)
nxgrid, nygrid : 자체 계산되는 직사각형의 셀 수(기본: 2,1)
inittau : initial value(위에 썼던 타우값, 기본: 0.02)
niter : 반복 몇 번 할건디? (기본: 1)
dist : 거리(기본: e)

Lv. 35 라이츄

Lv. 35 라이츄

광고 매크로 없는 청정한 블로그를 위해 노력중입니다. 근데 나만 노력하는 것 같음… ㅡㅡ

Biopython으로 Clusting analysis 하기 (실전편)

Coding/Python 2022. 8. 21. 01:29

와! 드디어 실전인가요? 근데 실전 생각보다 노잼임... 


Distance matrix

거리 행렬. 두 점간의 거리를 배열해 행렬로 나타낸 것이다. 점이 N개일 때 Distance matrix는 N*N으로 표기할 수 있다. 

import numpy as np
import pandas as pd
from Bio.Cluster import distancematrix
data=np.array([[0, 1,  2,  3],[4, 5,  6,  7],[8, 9, 10, 11],[1, 2,  3,  4]])
matrix = distancematrix(data)
# 뭐야 이거 왜 한영키 안먹어요
distances = distancematrix(data, dist='e')
print(distances)
[array([], dtype=float64), array([16.]), array([64., 16.]), array([ 1.,  9., 49.])]

뭐여 행렬 어디갔어요!!! 아 일단 나도 저 반응이어서 이해는 함... 근데 Distance matrix의 성질을 잘 이용하면 저걸 행렬로 만들 수 있다. 

여기서 이용할 Distance matrix의 성질 두 개는 

1. Distance matrix도 인접행렬처럼 주대각선(좌상-우하)이 0이다. 다 0이다. 
2. Distance matrix는 주대각선을 축으로 하는 대칭 행렬(대칭 행렬: 전치행렬이 원 행렬과 같음)이다. 

이거다. 이걸 이용해 저 array를 좌측 하단에 채우고, 대칭되는 방향으로 우측을 채우면 된다. 얘네 이거 안 알려줘서 내가 찾음. 

 

k- 들어가는 군집화

k-mean, k-median

from Bio.Cluster import kcluster
import numpy as np
data=np.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[0,1,2,3]])
clusterid, error, nfound = kcluster(data)
print("Clusterid:",clusterid, "error:",error, "nfound:",nfound)
Clusterid: [0 1 1 0] error: 8.0 nfound: 1

코드가 동일한데, method='a'(기본값)는 평균, method='m'은 중앙값. 위 코드는 평균으로 구한 것. 

 

k-medoid

from Bio.Cluster import kmedoids
from Bio.Cluster import distancematrix
import numpy as np
data=np.array([[0, 1,  2,  3],[4, 5,  6,  7],[8, 9, 10, 11],[1, 2,  3,  4]])
matrix = distancematrix(data)
# 뭐야 이거 왜 한영키 안먹어요
distances = distancematrix(data, dist='e')
clusterid, error, nfound = kmedoids(distances)
print("clusterid:",clusterid,"error:",error,"nfound:",nfound)
clusterid: [0 1 1 0] error: 17.0 nfound: 1

얘는 배열이 아니라 거리 행렬을 입력으로 받는다. 아까 그 괴랄한 형태도 되고, 1차원이나 2차원 배열로도 받는다. 

 

Hierarchical clustering

method

거리를 재는 방식에 따라 나뉜다. 거리 함수랑 별개. 

Pairwise single-linkage clustering: 제일 짧은 거리(s)
Pairwise maximum-linkage clustering: 제일 긴 거리(m)
Pairwise average-linkage clustering: 평균 거리(a)
Pairwise centroid-linkage clustering: 무게중심(centroid, c)

 

노드 생성

from Bio.Cluster import Node
A=Node(2,3) # left, right만 존재(거리=0)하는 노드
B=Node(2,3,0.91) # left, right, distance가 셋 다 존재하는 노드
print(B.left,B.right,B.distance)
2 3 0.91

입력인자는 왼쪽, 오른쪽, 거리. 거리는 option이라 입력 안 해도 된다. 

 

from Bio.Cluster import Node, Tree
nodes = [Node(1, 2, 0.2), Node(0, 3, 0.5), Node(-2, 4, 0.6), Node(-1, -3, 0.9)]
tree = Tree(nodes)
print(tree)
(1, 2): 0.2
(0, 3): 0.5
(-2, 4): 0.6
(-1, -3): 0.9

이런 식으로 트리를 그릴 수도 있다. 물론 Tree 모듈은 모셔와야 한다. 

 

print(tree[0])
# 생성된 tree에서 indexing을 통해 노드에 접근할 수 있다.
print(tree[0].left)
# 노드에 접근해서 노드의 정보를 가져올 수도 있다.
(1, 2): 0.2
1

인덱싱을 할 수도 있다. 트리의 노드를 인덱싱하거나, 아래 코드처럼 트리의 노드에 대한 정보(왼쪽, 오른쪽, 거리)를 가져오는 것도 된다. 

 

indices = tree.sort()
clusterid = tree.cut(nclusters=2)
print("indices: ",indices)
print("clusterid: ",clusterid)
# 누구 저 sort 어떻게 쓰는 지 아는 사람?
indices:  [7 0 1 2 1]
clusterid:  [1 1 1 1 1]

정렬이랑 분할(cut이 분할)도 된다는데 sort 저거 옵션 어떻게 주는건지는 안나왔다. 

 

본격적으로 해보기

from Bio.Cluster import treecluster
import numpy as np
data=np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[0,1,2,3]])
tree = treecluster(data)
print(tree)
(3, 0): 1
(2, 1): 16
(-2, -1): 81

Array로 clustering한 결과. 

 

tree = treecluster(data,dist="b",distancematrix=None)
print(tree)
# 다른 옵션을 줄 수도 있다.
(3, 0): 1
(2, 1): 4
(-2, -1): 9

거리 함수 옵션을 유클리드 거리에서 시티 블록으로 바꿨다. 

 

distances=distancematrix((data))
tree = treecluster(data=None,distancematrix=distances)
print(tree)
# Distance matrix를 미리 계산해 그걸로 그릴 수도 있다. 
# ValueError: use either data or distancematrix; do not use both
# Data와 Distance matrix중 하나는 None이어야 한다. 안그러면 위 에러가 반긴다.
(3, 0): 1
(2, 1): 16
(-2, -1): 81

같은 array로 그린건데 arrya가 아니라 거리 행렬을 구해서 그린 것. Data나 distancematrix는 둘 중 하나만 들어갈 수 있다. 그리고 Distance matrix로 hierarchical clustering을 할 때는 거리 측정 method 중 centroid(‘c’)를 선택할 수 없다. 

 

Self-Organizing Maps

대뇌피질의 시각피질을 모델화했다나… 모델화 한 사람 이름이 Teuvo Kohenen이라 Kohenen map이라고도 부른다. k-mean에 위상 개념이 추가된 것. 

 

Adjustment: 보정치(벡터에서 클러스터 찾을 때 쓰는 값
Radius: 반경

그렇다고 한다. 나도 잘 모름. 

 

from Bio.Cluster import somcluster
import numpy as np
data=np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[0,1,2,3]])
clusterid, celldata = somcluster(data)
print("clusterid:",clusterid,"\n","celldata:",celldata)
clusterid: [[1 0]
 [1 0]
 [1 0]
 [1 0]] 
 celldata: [[[-1.80107153 -0.6130221   0.16768853  0.59348614]]
 [[ 1.6816792  -0.16877007 -0.34863638  1.01090276]]]

 

Principal Component Analysis

주성분 분석. 고차원의 데이터를 저차원으로 축소시키는 기법 중 하나라고 한다. 아니 무슨 선형대수학이 왜 나오는거냐 대체... 와 진짜 대환장 파티네 뭐 주성분의 표본을 직교변환한다던가... 이런건 나도 모르겠고... 저는 수학이랑 담 쌓아서요... I HATE MATH 쌉가능

 

from Bio.Cluster import pca
import numpy as np
data=np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[0,1,2,3]])
columnmean, coordinates, components, eigenvalues = pca(data)
print("columnmean:",columnmean)
print("coordinates:",coordinates)
print("components:",components)
print("eigenvalues:",eigenvalues)
columnmean: [3.75 4.75 5.75 6.75]
coordinates: [[ 5.50000000e+00 -6.00596572e-16 -2.63550106e-49 -0.00000000e+00]
 [-2.50000000e+00 -1.16262943e-16  3.00734525e-47  0.00000000e+00]
 [-1.05000000e+01 -4.49845107e-16 -7.36989119e-48  0.00000000e+00]
 [ 7.50000000e+00 -2.28099977e-16 -1.00093434e-49  0.00000000e+00]]
components: [[-5.00000000e-01 -5.00000000e-01 -5.00000000e-01 -5.00000000e-01]
 [ 8.66025404e-01 -2.88675135e-01 -2.88675135e-01 -2.88675135e-01]
 [ 0.00000000e+00  8.16496581e-01 -4.08248290e-01 -4.08248290e-01]
 [ 0.00000000e+00 -7.19816116e-17 -7.07106781e-01  7.07106781e-01]]
eigenvalues: [1.42478068e+01 7.92857827e-16 3.09646140e-47 0.00000000e+00]

여기 표시된 결과는 대충 이렇다. 

Columnmean: 컬럼 평균(어레이 가로 말고 세로)
Coordinates: 주성분에 대한 데이터의 행 좌표
Components: 주성분
Eigenvalues: 각 주성분의 고유값

 

Cluster/TreeView 파일 핸들링하기

예제파일... 얘네가 깃헙에 올려뒀으니까 받아가세여... 

 

from Bio import Cluster
with open("/home/koreanraichu/cyano.txt") as handle:
    record = Cluster.read(handle)
    print(record)
<Bio.Cluster.Record object at 0x7ff32bb85af0>

불러와서 레코드만 달랑 달라고 하면 저렇게 준다. record.geneid 이런 식으로 세부 요소들을 요청하자.

 

from Bio import Cluster
with open("/home/koreanraichu/cyano.txt") as handle:
    record = Cluster.read(handle)
# 불러왔다
matrix = record.distancematrix()
# Distance matrix 계산
cdata, cmask = record.clustercentroids()
# Cluster 무게중심(Centroid) 계산
distance = record.clusterdistance()
# 클러스터간 거리 계산
tree = record.treecluster()
# hierarchical clustering
# 이거 matplot으로 못뽑나...
clusterid, error, nfound = record.kcluster()
# k-mean clustering
# method='a': k-mean
# method='m': k-median
clusterid, celldata = record.somcluster()
# SOM 계산하기
jobname="cyano_clustering"
record.save(jobname, record.treecluster(), record.treecluster(transpose=1))
# 내 컴퓨터에 저(별)장
# 기본 형식: record.save(jobname, geneclusters, expclusters)
# geneclusters=record.treecluster()
# expclusters=record.treecluster(transpose=1)

본인 깃헙에 있는 전체 소스. 저장의 경우 jobname, geneclusters, expclusters에 입력하는 만큼 파일이 생긴다. 

 

jobname만 입력하면 이거 하나 생긴다. 

 

셋 다 입력하면 세개 생긴다. 

 

from Bio import Cluster
with open("/home/koreanraichu/cyano.txt") as handle:
    record = Cluster.read(handle)
# 불러와서
genetree = record.treecluster(method="s")
genetree.scale()
# Hierarchical clustering
# method는 Pairwise single-linkage clustering
exptree = record.treecluster(dist="u", transpose=1)
# 이것도 Hierarchical clustering(컬럼으로 그림, dist)
# dist는 Uncentered Pearson correlation(어슨이형)
# 그려서
record.save("cyano_result", genetree, exptree)
# 저장

아무튼 이런 식으로 쓴다고... 

Lv. 35 라이츄

Lv. 35 라이츄

광고 매크로 없는 청정한 블로그를 위해 노력중입니다. 근데 나만 노력하는 것 같음… ㅡㅡ

Biopython으로 Clusting analysis 하기 (이론편)

Coding/Python 2022. 8. 21. 01:16

분량도 분량인데 spyder에서 한영키가 안돼서 그거땜시 늦었음... 

이게 한 글에 다 쓰기엔 좀 분량도 분량인데 이게 생각보다 설명이랑 코딩이랑 나뉘어있어서 이론편 실전편 나눕니다. 


이건 clustering 중 하나인 hierarchical clusting. 오늘 할 게 대충 이런거다. 

 

Cluster?

비슷한 특성을 가진 데이터 집단을 클러스터라고 한다. 데이터의 특성이 비슷하면 같은 클러스터, 다르면 다른 클러스터에 속한다. 클러스터링 하는 방법이 여러개가 있는데 여기서는 k-mean, k-median, k-medoid랑 hierarchical clustering에 대해 그냥 개 간단하게 설명하고 넘어간다. 

 

Hierarchical clustering

앞에 k-들어가는 것과 달리 계층적 클러스터링이라고 한다. 계층 나누는 방식에 따라 Top-down과 Bottop-top으로 나뉜다. 

젖산균을 젖산간균(Bacillus, 막대처럼 생긴 균)과 젖산구균(Coccus, 동그란거)으로 나눌 때 

1) Top-down 방식은 젖산균이라는 큰 카테고리를 먼저 잡고 젖산간균과 젖산구균으로 나눈다. 
2) Bottom-top 방식은 젖산간균과 젖산구균이라는 작은 카테고리를 잡고 그걸 젖산균으로 묶는다. 

아, 참고로 말씀드리는거지만 유산균이나 젖산균이나 그게 그거임다. 젖산=유산이라.. 

 

K-들어가는 방법

이건 생각보다 정의 자체가 간단하다. mean은 평균(산술평균), median은 중앙값, medoid는 중앙자(계산하는 게 아니라 실제 데이터 값 중 대표값을 선정한다). 

어떤 클러스터 내의 값이 81, 82, 83, 86일 때 

산술평균: 83
중앙값: 82.5

으로 각각 계산을 통해 값을 도출하게 되지만 중앙자는 저 네 개의 데이터 중에서 값을 하나 골라서 거리를 계산하고 클러스터링하는 방식. 아, 산술평균이 우리가 일반적으로 생각하는 그 평균이다. (다 더해서 갯수로 나누는 그거)


Distance function

초장부터 날 당황하게 만든 그거 맞다. 

 

Euclidean distance

Distance func. option: 'e'

유클리드 공간에서 두 점 사이의 거리를 나타내는 게 유클리드 거리. 따로 거리 함수를 지정해주지 않으면 이걸로 계산한다. 즉 Default값. 수학적인 정의는 뭔 소린지 모르겠지만 대충 피타고라스의 정리같은 거. 

 

City-block distance(Manhattan distance)

Distance func. option: 'b'

식만 보면 이게 뭔 개소리지 싶겠지만 아래 그림을 보면 뭔지 대충 이해가 될 것이다. 

초록색 선이 유클리드 거리이고, 나머지 선들이 맨하탄 거리(시티 블록 거리)이다. 저걸 좌표가 그려진 모눈종이라고 생각해보면, 유클리드 거리는 모눈이고 나바리고 대각선 ㄱㄱ하는거고 맨하탄 거리는 좌표 모눈 따라서 이동하는 것. 그래서 유클리드 거리가 좀 더 짧다. 

대신 도시에서는 맨하탄 거리로 계산하는 게 낫다. 도시의 경우 저 하얀 부분이 대부분 건물인데 길 찾는다고 건물 뚫고 갈거임? 

 

Pearson correlation coefficient

Distance func. option: 'c'

여기서부터 노션으로 정리하다가 사리나왔다. TEX 문법이 정말 (삐-)같거든... 

피어슨 상관 계수는 두 변수간의 선형적인 상관관계를 계량한 것으로, 바 달린건 평균이고 시그마 달린건 표준편차다. 피어슨 상관 계수가 1이면 양의 상관관계, -1이면 음의 상관관계를 갖고 0이면 상관관계가 없다. 이를 이용해 구하는 피어슨 거리의 경우 

이 공식으로 계산하므로 범위는 0~2. 2에 가까울수록 음의 상관관계, 0에 가까울수록 양의 상관관계를 갖고 1이면 상관관계가 없다. 

 

Absolute value of the Pearson correlation coefficient

Distance func. option: 'a'

피어슨 상관 계수에 절대값 끼는거. 위 거리 공식에서 r에 절대값이 들어간다. 즉 0에 가까울수록 어쨌든 상관관계가 존재하는 것이다. 

 

Uncentered correlation (cosine of the angle)

Distance func. option: 'u'

저거 사실 저렇게 검색하면 안나오니까 코사인 유사도(cosine similarity)로 검색하자. 

코사인 유사도는 이렇게 구하고 

저기 들어가는 시그마는 이런 식으로 구한다. 

위키피디아에 나온 공식은 저거. 맨 아래 식을 이용해 거리를 구할 수 있고, 상관계수 자체는 -1에서 1까지의 값을 가지고, 거리는 피어슨때랑 마찬가지로 0에서 2까지의 값을 가진다. 어슨이형꺼랑 범위는 똑같은데 얘는 상관계수가 -1이면 반대인거고 1이면 완전 같은 경우. 0이면 독립정인 경우다. 

 

Absolute uncentered Pearson correlation

Distance func. option: 'x'

위에 거리 공식에서 r에 절대값을 준 것. 

그러니까 이거다. 

 

Spearman’s rank correlation

Distance func. option: 's'

일단 이게 최종 계산 공식. 

피어슨 상관 계수와 공분산을 쓰면 이렇게 계산한다. cov=공분산, rho=피어슨 상관 계수. 

스피어맨 상관 계수는 피어슨 상관 계수와 달리 두 변수의 '순위' 사이의 통계적 의존성을 측정하는 척도이다. 두 변수 사이의 관계가 단조함수(주어진 순서를 보존하는 함수)를 통해 얼마나 잘 설명될 수 있는지를 평가한다. 즉 스피어맨 상관 계수는 변수에 순위가 있을 때의 피어슨 상관 계수. 단, 피어슨 상관 계수가 선형인 것과 달리 얘는 단조함수기 때문에 스피어맨 상관 계수가 1이라고 피어슨도 1인 건 아니다. 

 

Kendall’s τ

Distance func. option: 'k'

켄달 타우라고 부른다. 대충 

이게 타우인데... K랑 L이 뭐냐... 

켄달타우도 스피어맨처럼 순위가 존재하는 변수에 대해 매길 수 있는데, K와 L은 각각 순위가 존재하는 두 변수 중 아무거나 픽했을 때 순위가 일치하는가 여부이다. 이게 무슨 소리냐... 체력장을 예시로 들어보자. 두 사람의 100미터 달리기와 1000미터 달리기의 기록을 비교했을 때 아래와 같은 가짓수가 나오게 된다. 

1. 100m, 1000m 둘 다 빠르다
2. 100m, 1000m 둘 다 느리다
3. 100m에서는 빠르지만 1000m에서는 느리다 
4. 100m에서는 느리지만 1000m에서는 빠르다 

이 때 1, 2번같은 케이스를 concordant pair이라고 하게 되고, 변수 내에 concordant pair가 존재하면 K가 1, concordant pair가 아닌 케이스(3, 4번 케이스)가 존재하면 L이 1. L이 0인 경우가 많으면 타우(켄달타우)가 1에 가까워지고 K가 0인 경우가 많으면 타우가 -1에 가까워진다. 만약 0이라고요? 그러면 둘이 상관 없는거다. 

L 공식에 저 꺾쇠? NOT이다. 

Lv. 35 라이츄

Lv. 35 라이츄

광고 매크로 없는 청정한 블로그를 위해 노력중입니다. 근데 나만 노력하는 것 같음… ㅡㅡ

Biopython으로 Sequence motif analysis 하기

Coding/Python 2022. 8. 21. 01:04

모티브 찾으면 구조 모티브가 나오는데 이거는 '단백질이나 핵산과 같은 사슬 모양의 생체 분자에서 진화적으로 관련이 없는 다양한 분자에서 나타나는 일반적인 3차원 구조'로 정의한다. 근데 여기서 다루는 모티브는 그거 말고 시퀀스 모티브... 

대충 이런거다. 위 그림은 뭔 시퀀스인지는 모르겠으나 3, 4, 5번째 염기가 GAA가 압도적으로 많은 듯. 

오늘은 대충 

from Bio import motifs

이런거 부른다. 


Motif object

객체 생성하기

from Bio import motifs
from Bio.Seq import Seq
instances=[Seq("TGTCGTATCG"),Seq("GTAAATAGCC"),Seq("GTAAATAACC"),Seq("TCGCGGAGCC"),Seq("ATGTGCCATA"),Seq("CTCGTCTGCG")]
m = motifs.create(instances)
print(m)
TGTCGTATCG
GTAAATAGCC
GTAAATAACC
TCGCGGAGCC
ATGTGCCATA
CTCGTCTGCG

이런 식으로 시퀀스 객체를 활용해서 만들 수 있다. 저거는 DNA 랜덤으로 만든거. len() 쓰면 길이도 알려준다. 

 

print(m.counts)
        0      1      2      3      4      5      6      7      8      9
A:   1.00   0.00   2.00   2.00   2.00   0.00   4.00   2.00   0.00   1.00
C:   1.00   1.00   1.00   2.00   0.00   2.00   1.00   0.00   5.00   3.00
G:   2.00   1.00   2.00   1.00   3.00   1.00   0.00   3.00   0.00   2.00
T:   2.00   4.00   1.00   1.00   1.00   3.00   1.00   1.00   1.00   0.00

저 표가 뭔 소리냐면... 일단 저 DNA가 10bp짜리 6개 시퀀스로 만든건데(그래서 0부터 9까지), 6개 시퀀스의 0번 컬럼을 다 봤더니 아데닌과 시토신이 있는 게 하나, 구아닌과 티민이 있는 게 두개라는 얘기. 

 

print(m.counts["A"])
[1, 0, 2, 2, 2, 0, 4, 2, 0, 1]

인덱싱은 넘파이 array처럼 2차원으로 할 수 있다. 저건 염기(인덱스)로 뽑은거. 

 

print(m.counts["A",0])
1

첫번째 염기에 아데닌으로 뽑은 거. 

 

print(m.counts[:,0])
{'A': 1, 'C': 1, 'G': 2, 'T': 2}

컬럼 인덱싱도 된다. 

 

print(m.consensus)
print(m.anticonsensus)
GTAAGTAGCC
AACGCAGCAT

순서대로 

-가장 많은 염기
-가장 적은 염기

를 표시해준다. 

 

print(m.degenerate_consensus)
NTNNRYARCS

...... 어머니... 이게 무슨 개소리죠... 

 

당황하지 말고 이걸 보자. 

이것도 있다. (FASTA 문서에 있음)

그니까 저기 다섯번째에 있는 R은 '아 이거 내가 셔봤는데 대충 퓨린이 많어' 이런 얘기. 

 

r = m.reverse_complement()
print(r)
CGATACGACA
GGCTATTTAC
GGTTATTTAC
GGCTCCGCGA
TATGGCACAT
CGCAGACGAG

reverse complement 만들어서 셀 수도 있다. 

 

Weblogo 만들기

m.weblogo("mymotif.png")

님아 사이즈좀... 

 

Data 읽기

이 와중에 이것도 DB가 있다. JASPAR라고... 

https://jaspar.genereg.net/

 

JASPAR: An open-access database of transcription factor binding profiles

JASPAR is the largest open-access database of curated and non-redundant transcription factor (TF) binding profiles from six different taxonomic groups.

jaspar.genereg.net

여기다. 

 

sites format은 못구했고 jaspar format이랑 pfm 읽어봅시다. meme이랑 tr뭐시기는 에러났음. 

 

from Bio import motifs
fh = open("/home/koreanraichu/MA0004.1.jaspar")
for m in motifs.parse(fh, "jaspar"):
    print(m)
TF name	Arnt
Matrix ID	MA0004.1
Matrix:
        0      1      2      3      4      5
A:   4.00  19.00   0.00   0.00   0.00   0.00
C:  16.00   0.00  20.00   0.00   0.00   0.00
G:   0.00   1.00   0.00  20.00   0.00  20.00
T:   0.00   0.00   0.00   0.00  20.00   0.00

jaspar format

 

from Bio import motifs
with open("Arnt.sites") as handle:
    arnt = motifs.read(handle, "sites")

sites format은 이렇게 부른다. 

 

with open("/home/koreanraichu/MA0007.1.pfm") as handle:
    srf = motifs.read(handle, "pfm")
print(srf)

pfm은 이렇게 부르면 되는데 문제가 하나있다. 

저 위에 있는 ID를 인식을 못 해서 저거 저대로 부르면 ValueError: could not convert string to float: '>MA0007.1' 이 여러분을 반긴다. 

TF name	None
Matrix ID	None
Matrix:
        0      1      2      3      4      5      6      7      8      9     10     11     12     13     14     15     16     17     18     19     20     21
A:   9.00   9.00  11.00  16.00   0.00  12.00  21.00   0.00  15.00   4.00   5.00   6.00   3.00   0.00   4.00  11.00   1.00   3.00   6.00   6.00  10.00   5.00
C:   7.00   2.00   3.00   1.00   0.00   6.00   2.00  24.00   0.00   9.00  11.00   9.00   5.00   0.00   0.00   5.00  22.00  16.00   7.00   5.00  11.00  11.00
G:   2.00   3.00   3.00   7.00  22.00   1.00   0.00   0.00   8.00   2.00   2.00   9.00   1.00  24.00   0.00   1.00   1.00   1.00   5.00   9.00   0.00   6.00
T:   6.00  10.00   7.00   0.00   2.00   5.00   1.00   0.00   1.00   9.00   6.00   0.00  15.00   0.00  20.00   7.00   0.00   4.00   6.00   4.00   3.00   2.00

ID 지우면 이렇게 부를 수는 있는데 그럴거면 불러서 뭐함? 

 

from Bio import motifs
with open("/home/koreanraichu/MA0273.1.meme") as handle:
    record = motifs.parse(handle, "meme")
    print(record)

Meme은 이렇게 부르는데 

이게 이 형식이그덩... 그랬더니 ValueError: Improper MEME XML input file. XML root tag should start with <MEME version= ... 뜸... 

 

from Bio import motifs
with open("/home/koreanraichu/MA0273.1.transfac") as handle:
    record = motifs.parse(handle, "TRANSFAC")
    print(record)

transfac은 이렇게 부르는데 

얘는 또 ValueError: A TRANSFAC motif line should have 2 spaces between key and value columns: "AC MA0273.1"가 반긴다. 니네 뭐하자는... 

 

Data 쓰기

from Bio import motifs
fh = open("/home/koreanraichu/MA0004.1.jaspar")
for m in motifs.parse(fh, "jaspar"):
    print(m.format("pfm"))
  4.00  19.00   0.00   0.00   0.00   0.00
 16.00   0.00  20.00   0.00   0.00   0.00
  0.00   1.00   0.00  20.00   0.00  20.00
  0.00   0.00   0.00   0.00  20.00   0.00

Jaspar->pfm

 

from Bio import motifs
fh = open("/home/koreanraichu/MA0004.1.jaspar")
for m in motifs.parse(fh, "jaspar"):
    print(m.format("transfac"))
P0      A      C      G      T
01      4     16      0      0      C
02     19      0      1      0      A
03      0     20      0      0      C
04      0      0     20      0      G
05      0      0      0     20      T
06      0      0     20      0      G
XX
//

Jaspar->transfac (meme은 지원 안함)

 

instances=[Seq("TGTCGTATCG"),Seq("GTAAATAGCC"),Seq("GTAAATAACC"),Seq("TCGCGGAGCC"),Seq("ATGTGCCATA"),Seq("CTCGTCTGCG")]
m2 = motifs.create(instances)
print(m2.format("jaspar"))
>None 
A [  1.00   0.00   2.00   2.00   2.00   0.00   4.00   2.00   0.00   1.00]
C [  1.00   1.00   1.00   2.00   0.00   2.00   1.00   0.00   5.00   3.00]
G [  2.00   1.00   2.00   1.00   3.00   1.00   0.00   3.00   0.00   2.00]
T [  2.00   4.00   1.00   1.00   1.00   3.00   1.00   1.00   1.00   0.00]

물론 jaspar format으로도 된다. 

 

PWM(position weight matrices)

pwm = m.counts.normalize(pseudocounts=0.5)
print(pwm)
        0      1      2      3      4      5
A:   0.20   0.89   0.02   0.02   0.02   0.02
C:   0.75   0.02   0.93   0.02   0.02   0.02
G:   0.02   0.07   0.02   0.93   0.02   0.93
T:   0.02   0.02   0.02   0.02   0.93   0.02

normalize가 들어간 걸 보면 정규화를 하는 모양. 정규화나 표준화는 원래 인공지능용 Data 전처리 과정에서 진행한다. 

 

pwm = m.counts.normalize(pseudocounts={"A":0.6, "C": 0.4, "G": 0.4, "T": 0.6})
print(pwm)
        0      1      2      3      4      5
A:   0.21   0.89   0.03   0.03   0.03   0.03
C:   0.75   0.02   0.93   0.02   0.02   0.02
G:   0.02   0.06   0.02   0.93   0.02   0.93
T:   0.03   0.03   0.03   0.03   0.94   0.03

염기별로 pseudocount를 지정할 수도 있다. 

 

rpwm = pwm.reverse_complement()
print(rpwm)
# reverse complement로도 된다. (아마 밑에 코드로 계산된 듯)
        0      1      2      3      4      5
A:   0.03   0.94   0.03   0.03   0.03   0.03
C:   0.93   0.02   0.93   0.02   0.06   0.02
G:   0.02   0.02   0.02   0.93   0.02   0.75
T:   0.03   0.03   0.03   0.03   0.89   0.21

reverse complement로도 된다. 

 

print(pwm.consensus)
print(pwm.anticonsensus)
print(pwm.degenerate_consensus)
# 제일 많은거/제일 적은거/Degenerate consensus
CACGTG
GCGCCC
CACGTG

위에 세 개를 저걸로 결론내는 모양. 아 그래서 2:2인데 아데닌 나왔냐 

 

PSSM(position specific scoring matrices)

이거는 PWM이 있어야 구할 수 있는 모양. 

 

pwm = m.counts.normalize(pseudocounts=0.5) # 예제 코드를 보니 얘가 있어야 계산할 수 있는 모양
pssm = pwm.log_odds()
print(pssm)
        0      1      2      3      4      5
A:  -0.29   1.83  -3.46  -3.46  -3.46  -3.46
C:   1.58  -3.46   1.90  -3.46  -3.46  -3.46
G:  -3.46  -1.87  -3.46   1.90  -3.46   1.90
T:  -3.46  -3.46  -3.46  -3.46   1.90  -3.46

저기서 양수면 모티브에서, 음수면 백그라운드에서 많이 보인다. 백그라운드에서 많이 보인다는 건 '주'는 아닌 듯... 0이면? 또이또이 쌤쌤. 

 

background = {"A":0.3,"C":0.2,"G":0.2,"T":0.3}
pssm = pwm.log_odds(background)
print(pssm)
# background를 염기별로 줄 수도 있다.
        0      1      2      3      4      5
A:  -0.55   1.56  -3.72  -3.72  -3.72  -3.72
C:   1.91  -3.14   2.22  -3.14  -3.14  -3.14
G:  -3.14  -1.55  -3.14   2.22  -3.14   2.22
T:  -3.72  -3.72  -3.72  -3.72   1.64  -3.72

염기별로 background setting도 된다. 최대/최소값을 볼 수 있고, background를 임의로 부여한 경우 평균과 표준편차도 구할 수 있다. 

 

print("%4.2f" % pssm.max) # 최대값
print("%4.2f" % pssm.min) # 최소값
mean = pssm.mean(background) # 평균
std = pssm.std(background) # 표준편차
print("mean = %0.2f, standard deviation = %0.2f" % (mean, std))

 

그 이름 검⭐색

시퀀스로 찾기

from Bio import motifs
from Bio.Seq import Seq
test_seq=Seq("CTCGTCTGCG")
m = motifs.create(instances=[Seq("TGTCGTATCG"),Seq("GTAAATAGCC"),Seq("GTAAATAACC"),Seq("TCGCGGAGCC"),Seq("ATGTGCCATA"),Seq("CTCGTCTGCG")])
for pos, seq in m.instances.search(test_seq):
    print(pos,seq)
0 CTCGTCTGCG

test_seq에서 motif를 찾는거기때문에 test_seq보다 motif가 더 길면 에러난다. 

 

PSSM score로 찾기

pwm = m.counts.normalize(pseudocounts=0.5) # 예제 코드를 보니 얘가 있어야 계산할 수 있는 모양
pssm = pwm.log_odds()
print(pssm)
# PSSM으로 찾기(PSSM 도출)
for position, score in pssm.search(test_seq, threshold=3.0):
    print("Position %d: score = %5.3f" % (position, score))
# 찾았음!
Position 0: score = 2.005

threshold요? 그거 밑에 나옴. 아마 3.0이 기준인 듯 하다. 

 

Threshold

background = {"A":0.3,"C":0.2,"G":0.2,"T":0.3}
distribution = pssm.distribution(background=background, precision=10**4)
threshold = distribution.threshold_fpr(0.01)
print("%5.3f" % threshold)
3.784

있으니까 하긴 한다만... 이거 뭐예여... 

 

threshold = distribution.threshold_fnr(0.1)
print("%5.3f" % threshold)
-0.664

일단 위에껀 fpr 밑에껀 fnr이다.

-fpr: false-positive rate (찾을 확률)
-fnr: false-negative rate (못찾을 확률)
-balanced: false-positive rate and the false-negative rate (fnr/fpr≃ t)
-patser: equality between the −log of the false-positive rate and the information content

이렇다고 한다. 저걸로 threshold 계산해서 그거 때려박고 검색할 수도 있다.

 

background = {"A":0.3,"C":0.2,"G":0.2,"T":0.3}
distribution = pssm.distribution(background=background, precision=10**4)
threshold = distribution.threshold_fnr(0.1)
print("%5.3f" % threshold)
for position, score in pssm.search(test_seq, threshold=threshold):
    print("Position %d: score = %5.3f" % (position, score))
Position 0: score = 6.703
Position -22: score = 0.084
Position 5: score = -0.301
Position 9: score = 0.195
Position 12: score = -0.346
Position -12: score = 0.681
Position 15: score = 2.992

이게 threshold 계산해서 때려박은 결과

 

Position 0: score = 6.703
Position -22: score = 0.084
Position 9: score = 0.195
Position -12: score = 0.681
Position 15: score = 2.992

이건 똑같은 모티브에서 고정값threshold로 찾은 결과이다. 

 

PSSM과 motif

print(motif.counts)
print(motif.pwm)
print(motif.pssm)
# 일단 뽑고 보는 타입
0      1      2      3      4      5      6      7      8      9     10     11     12     13     14     15     16     17     18     19     20
A:  64.00 104.00 484.00  78.00 455.00 248.00  98.00  52.00 348.00  12.00 966.00  26.00   9.00  47.00 892.00 487.00 288.00 373.00 178.00 402.00 435.00
C: 217.00 552.00 111.00 401.00  51.00  34.00 724.00 918.00 346.00   3.00   5.00   6.00  10.00 930.00  42.00 123.00 140.00 110.00 367.00 140.00 229.00
G: 425.00 150.00 114.00 186.00 370.00 670.00 143.00   7.00 280.00 977.00  23.00 962.00   4.00   7.00  16.00 320.00 317.00  81.00 184.00 341.00 241.00
T: 292.00 192.00 288.00 333.00 122.00  46.00  33.00  20.00  24.00   6.00   4.00   4.00 975.00  15.00  49.00  68.00 254.00 434.00 268.00 114.00  94.00

        0      1      2      3      4      5      6      7      8      9     10     11     12     13     14     15     16     17     18     19     20
A:   0.06   0.10   0.49   0.08   0.46   0.25   0.10   0.05   0.35   0.01   0.97   0.03   0.01   0.05   0.89   0.49   0.29   0.37   0.18   0.40   0.44
C:   0.22   0.55   0.11   0.40   0.05   0.03   0.73   0.92   0.35   0.00   0.01   0.01   0.01   0.93   0.04   0.12   0.14   0.11   0.37   0.14   0.23
G:   0.43   0.15   0.11   0.19   0.37   0.67   0.14   0.01   0.28   0.98   0.02   0.96   0.00   0.01   0.02   0.32   0.32   0.08   0.18   0.34   0.24
T:   0.29   0.19   0.29   0.33   0.12   0.05   0.03   0.02   0.02   0.01   0.00   0.00   0.98   0.02   0.05   0.07   0.25   0.43   0.27   0.11   0.09

        0      1      2      3      4      5      6      7      8      9     10     11     12     13     14     15     16     17     18     19     20
A:  -1.96  -1.26   0.96  -1.68   0.87  -0.01  -1.35  -2.26   0.48  -4.38   1.95  -3.26  -4.79  -2.41   1.84   0.96   0.21   0.58  -0.49   0.69   0.80
C:  -0.20   1.15  -1.17   0.68  -2.29  -2.88   1.54   1.88   0.47  -6.38  -5.64  -5.38  -4.64   1.90  -2.57  -1.02  -0.84  -1.18   0.56  -0.83  -0.13
G:   0.77  -0.73  -1.13  -0.42   0.57   1.43  -0.80  -5.15   0.17   1.97  -3.44   1.95  -5.96  -5.16  -3.96   0.36   0.34  -1.62  -0.44   0.45  -0.05
T:   0.23  -0.38   0.21   0.42  -1.03  -2.44  -2.92  -3.64  -3.38  -5.38  -5.96  -5.96   1.97  -4.06  -2.35  -1.88   0.02   0.80   0.10  -1.13  -1.41

일단 불러서 뽑아보면 이런 식으로 나온다. (순서대로 count, PWM, PSSM)

 

print(motif.background)
print(motif.pseudocounts)
# 일단 불러온 데이터에서 뽑을 수는 있다.
{'A': 0.25, 'C': 0.25, 'G': 0.25, 'T': 0.25}
{'A': 0.0, 'C': 0.0, 'G': 0.0, 'T': 0.0}

Motif 데이터에 손대지 않은 default값. 

 

motif.pseudocounts = 3.0
motif.background = {"A": 0.2, "C": 0.3, "G": 0.3, "T": 0.2}
# 어? 되네?
print(motif.counts)
print(motif.pwm)
print(motif.pssm)
# 바꾸고 계산해서 보자 한번.
print(motif.background)
print(motif.pseudocounts)
0      1      2      3      4      5      6      7      8      9     10     11     12     13     14     15     16     17     18     19     20
A:   0.07   0.11   0.48   0.08   0.45   0.25   0.10   0.05   0.35   0.01   0.96   0.03   0.01   0.05   0.89   0.49   0.29   0.37   0.18   0.40   0.43
C:   0.22   0.55   0.11   0.40   0.05   0.04   0.72   0.91   0.35   0.01   0.01   0.01   0.01   0.92   0.04   0.12   0.14   0.11   0.37   0.14   0.23
G:   0.42   0.15   0.12   0.19   0.37   0.67   0.14   0.01   0.28   0.97   0.03   0.96   0.01   0.01   0.02   0.32   0.32   0.08   0.19   0.34   0.24
T:   0.29   0.19   0.29   0.33   0.12   0.05   0.04   0.02   0.03   0.01   0.01   0.01   0.97   0.02   0.05   0.07   0.25   0.43   0.27   0.12   0.10

        0      1      2      3      4      5      6      7      8      9     10     11     12     13     14     15     16     17     18     19     20
A:  -1.59  -0.92   1.27  -1.32   1.18   0.31  -1.00  -1.88   0.80  -3.75   2.26  -2.80  -4.07  -2.02   2.15   1.28   0.53   0.90  -0.16   1.00   1.12
C:  -0.46   0.87  -1.41   0.42  -2.49  -3.03   1.26   1.61   0.20  -5.66  -5.24  -5.07  -4.54   1.62  -2.75  -1.27  -1.08  -1.42   0.29  -1.08  -0.39
G:   0.50  -0.99  -1.37  -0.68   0.30   1.15  -1.05  -4.92  -0.10   1.69  -3.54   1.67  -5.44  -4.92  -4.00   0.09   0.08  -1.85  -0.69   0.18  -0.31
T:   0.55  -0.05   0.53   0.73  -0.69  -2.04  -2.49  -3.13  -2.90  -4.49  -4.85  -4.85   2.28  -3.49  -1.96  -1.51   0.35   1.11   0.43  -0.79  -1.06

{'A': 0.2, 'C': 0.3, 'G': 0.3, 'T': 0.2}
{'A': 3.0, 'C': 3.0, 'G': 3.0, 'T': 3.0}

아, 가져온 데이터에서 바꿔줄 수는 있다. 저게 해본거임. 

 

motif.background = None
print(motif.background)
# None: 각 염기별로 0.25씩
motif.background = 0.9
print(motif.background)
# GC 비중보게.
# {'A': 0.04999999999999999, 'C': 0.45, 'G': 0.45, 'T': 0.04999999999999999} 는 심한 거 아니냐 근데.

뭐 이렇게 쓰기도 한다고... 

저 값을 바꿔주게 되면 당연히 평균, 표준편차, Threshold도 바뀐다. 

 

Motif 비교하기

참고로 이거 일단 두 개 불렀다. 

 

motif_1.pseudocounts = {"A":0.6, "C": 0.4, "G": 0.4, "T": 0.6}
motif_1.background = {"A":0.3,"C":0.2,"G":0.2,"T":0.3}
# pseudocount와 background를 설정해준다(깔맞춤)
pssm_motif_1=motif_1.pssm
print(pssm_motif_1)
# pssm을 생성한다
0      1      2      3      4      5      6      7      8      9     10     11     12     13     14
A:   0.78  -0.06  -0.33  -1.50   1.55   1.66   1.66  -1.27  -1.30  -2.21   1.59   1.62   1.59  -0.82  -0.40
C:  -0.57   0.33  -1.28  -2.30  -3.11  -3.63  -4.15   1.65  -1.74  -3.31  -3.39  -3.15  -3.10   0.69   0.39
G:  -0.51   0.30  -0.62   2.02  -1.61  -3.60  -3.70   0.08  -3.51   2.16  -2.01  -2.65  -2.36   0.94  -0.70
T:  -0.52  -0.48   0.87  -2.83  -3.28  -3.81  -3.76  -3.06   1.42  -3.98  -3.47  -3.69  -3.24  -1.28   0.37

귀찮아서 앞쪽 코드를 생략하긴 했는데 jaspar 파일을 두 개 불러서 background와 pesudocount를 하나로 맞추고 PSSM을 도출했다. 그래야 비교가 가능하다. 

 

distance, offset = pssm_motif_2.dist_pearson(pssm_motif_1)
print("distance = %5.3g" % distance)
print("offset: ",offset)
distance = 0.286
offset:  -14

저거 예제 코드에 pssm.dist_pearson으로 되어 있던거 걍 복붙했더니 에러나던데? 그런 이름 없셈! 을 토하더라… 참고로 저 코드에 같은 pssm이 들어가면 당연하게도 distance나 offset이 0이 된다. 

Lv. 35 라이츄

Lv. 35 라이츄

광고 매크로 없는 청정한 블로그를 위해 노력중입니다. 근데 나만 노력하는 것 같음… ㅡㅡ

Biopython으로 Phylogenetic tree 그리기

Coding/Python 2022. 8. 21. 00:44

걍 메가 쓰세여... 메가가 편해... 짱이야 메가... 


계통수?

계통수... 그러니까 Phylogenetic tree는 

이런거다. (...)

유전자나 단백질 시퀀스 분석(균이라면 16s rRNA라던가)을 통해 얘네가 얼마나 가까운지를 알아내게 되면 그걸 저런 식으로 그려서 나타내는 것. 저렇게 생물 종에 따라 그리는 경우도 있고, 특정 단백질의 homolog나 다른 생물종에서 같은 역할을 하는 단백질에 대해서 저걸 그리기도 한다. 

 

와! 계통수! 그려보자! 

from Bio import Phylo
tree = Phylo.read("/home/koreanraichu/Deinococcus.ph", "newick")

이걸 쓰면 그릴 수 있는데... 어디가 일로와 끝까지 듣고 가... 

저것만 쓰면 그려는 주는데 결과물이 없다. 뭐야 트리 돌려줘요! 쌉가능. 

Phylo.draw_ascii(tree)
# Ascii tree

그러니 이걸 끼워넣어보자. 

 

_______________________________________ EU583728.1.1469
  ___|
 |   |_________________________________________ AB195680.1.1544
 |
 |                                           , FJ808613.1.1430
 |                            _______________|
 |                           |               | AY553296.1.1427
 |         __________________|
 |        |                  |  _____________ KJ941155.1.1415
 |________|                  |_|
 |        |                    |_____________ AY667493.1.1404
 |        |
 |        |__________________________________ KF963825.1.1460
 |
 | __________________________________________ AY351389.1.1492
 ||
 ||                                         , CP013862.2043792.2045365
 ||                                         |
 ,|                                         , CP013862.2278622.2280195
 ||                                         |
 ||                                         , CP013862.2582000.2583573
 ||                                         |
 ||                                         | CP013862.2587878.2589451
_||_________________________________________|
 |                                          | CP013862.2511797.2513370
 |                                          |
 |                                          | CP013862.1437203.1438776
 |
 |                                              _ EF612763.1.1522
 |                                          ___|
 |   ______________________________________|   | AB681790.1.1492
 |  |                                      |
 |__|                                      |____ AB016721.1.1502
 |  |
 |  |__________________________________________ AB021192.1.1513
 |  |
 |  |                                       , AGAV01000003.54672.56217
 |  |_______________________________________|
 |                                          | AGAV01000003.371388.372933
 |
 |  __________________________________________ AB116124.1.1501
 |_|
   |__________________________________________ AB116122.1.1490

뭔진 모르겠지만 잘못했다... 아, 내가 미리 말하는 걸 깜빡했는데 본인은 이거 하기 전에 미리 ClustalW로 ph파일 만들었다. (ph파일: Phylogenic tree) 여러분도 미리 만들어놓고 쓰자. 오늘의 도우미 관련해서는 따로 서술하도록 한다. 

 

Phylo.draw(tree)

이 코드를 쓰면 

이런 식으로 matplot으로 출력해준다. (위에 글자는 ASCII)

 

계통수에 색깔 입히기

tree.root.color = (255,0,0)
tree.root.color = "#f7cac9"
tree.root.color = "gray"
# 가지 색을 바꿀 수 있다.

셋 중 하나만 써도 된다. 

 

근데 빨간색은 좀 심했다... 

 

tree.clade[0, 1].color = "blue"
# 부분적으로 색을 바꿔줄 수도 있는 듯.

부분적으로 바꿀거면 clade 인덱싱해서 바꾸면 된다. 

 

이런 식. 

Phylo.draw 윗줄에 있어야 적용된다. 

 

Phyloxml

뭔진 모르겠지만 일단 있으니까 본다... 

 

import sys
Phylo.write(tree, sys.stdout, "phyloxml")
# xml format으로 저장...은 안됨?

저장은 안된다. 

 

IO function

from Bio import Phylo
tree = Phylo.read("/home/koreanraichu/Halobacterium.ph", "newick")
print(tree)
Tree(rooted=False, weight=1.0)
    Clade()
        Clade(branch_length=0.01498)
            Clade(branch_length=0.01948)
                Clade(branch_length=0.36544)
                    Clade(branch_length=0.00448, name='AB021386.1.1495')
                    Clade(branch_length=0.00824, name='AF211861.1.1495')
                Clade(branch_length=0.01186)
                    Clade(branch_length=0.35417, name='AY519654.1.1502')
                    Clade(branch_length=0.36043, name='CP010835.91794.93350')
            Clade(branch_length=0.12384)
                Clade(branch_length=0.0895)
                    Clade(branch_length=0.15987, name='AB041346.1.1480')
                    Clade(branch_length=0.11499, name='AY099168.1.1408')
                Clade(branch_length=0.06337)
                    Clade(branch_length=0.18487, name='AY647169.1.1470')
                    Clade(branch_length=0.20195, name='EU522083.1.1394')
        Clade(branch_length=0.01418)
            Clade(branch_length=0.13573)
                Clade(branch_length=0.13377)
                    Clade(branch_length=0.05661, name='AB074299.1.1473')
                    Clade(branch_length=0.03326)
                        Clade(branch_length=0.02292)
                            Clade(branch_length=0.00218)
                                Clade(branch_length=0.0, name='AB663366.1.1474')
                                Clade(branch_length=0.0, name='AB663368.1.1474')
                            Clade(branch_length=0.002, name='AB663369.1.1474')
                        Clade(branch_length=0.02655, name='AB663367.1.1475')
                Clade(branch_length=0.19268, name='AJ277205.1.1505')
            Clade(branch_length=0.14735)
                Clade(branch_length=0.07756)
                    Clade(branch_length=0.1193, name='AB105159.1.1501')
                    Clade(branch_length=0.11049, name='AY529491.1.1465')
                Clade(branch_length=0.1948, name='CEML01000001.1297981.1299440')
        Clade(branch_length=0.23141)
            Clade(branch_length=0.11327)
                Clade(branch_length=0.006, name='AB663359.1.1473')
                Clade(branch_length=0.00587, name='KC914879.1.1473')
            Clade(branch_length=0.11413, name='U20163.1.1495')

얘도 read랑 parse가 있다. 그리고 parse하면? 그죠 for문이져... 

 

trees = Phylo.parse("/home/koreanraichu/Deinococcus.xml", "phyloxml")
for tree in trees:
    print(tree)
# parse(phyloxml)

예제 코드대로 했더니 뭐 뻑나서 봤는데 변수명이 틀렸더만... 

 

Phylo.convert("tree1.nwk", "newick", "tree1.xml", "nexml")
Phylo.convert("other_trees.xml", "phyloxml", "other_trees.nex", "nexus")

변환도 된단다. 

 

View and export

from Bio import Phylo
tree = Phylo.read("/home/koreanraichu/Halobacterium.ph", "newick")
Phylo.draw(tree, branch_labels=lambda c: c.branch_length)
# 거리 들어가는 건 좋은데 되게 중구난방이시네요.

익명함수 람다의 힘을 빌리면 계통수에 거리를 넣을 수 있다. 

 

tree.name="Halophile"

그리고 제목 마려우면 이거 쓰고. 

 

Tree and Clade

뭐...뭔지 모르겠다... (기분탓 아님)

 

정보 내놔! 드리겠습니다! 

print(tree.get_terminals()) # tree의 끝 노드를 출력
print(tree.get_nonterminals()) # tree의 중간 노드까지 출력

위 코드는 tree의 끝 노드만, 아래 노드는 중간 노드까지 출력해주는데 왜 위쪽 노드 결과값이 더 긴 지는 묻지 말자. 

 

[Clade(branch_length=0.00448, name='AB021386.1.1495'), Clade(branch_length=0.00824, name='AF211861.1.1495'), Clade(branch_length=0.35417, name='AY519654.1.1502'), Clade(branch_length=0.36043, name='CP010835.91794.93350'), Clade(branch_length=0.15987, name='AB041346.1.1480'), Clade(branch_length=0.11499, name='AY099168.1.1408'), Clade(branch_length=0.18487, name='AY647169.1.1470'), Clade(branch_length=0.20195, name='EU522083.1.1394'), Clade(branch_length=0.05661, name='AB074299.1.1473'), Clade(branch_length=0.0, name='AB663366.1.1474'), Clade(branch_length=0.0, name='AB663368.1.1474'), Clade(branch_length=0.002, name='AB663369.1.1474'), Clade(branch_length=0.02655, name='AB663367.1.1475'), Clade(branch_length=0.19268, name='AJ277205.1.1505'), Clade(branch_length=0.1193, name='AB105159.1.1501'), Clade(branch_length=0.11049, name='AY529491.1.1465'), Clade(branch_length=0.1948, name='CEML01000001.1297981.1299440'), Clade(branch_length=0.006, name='AB663359.1.1473'), Clade(branch_length=0.00587, name='KC914879.1.1473'), Clade(branch_length=0.11413, name='U20163.1.1495')]

(위 코드 결과)

[Clade(), Clade(branch_length=0.01498), Clade(branch_length=0.01948), Clade(branch_length=0.36544), Clade(branch_length=0.01186), Clade(branch_length=0.12384), Clade(branch_length=0.0895), Clade(branch_length=0.06337), Clade(branch_length=0.01418), Clade(branch_length=0.13573), Clade(branch_length=0.13377), Clade(branch_length=0.03326), Clade(branch_length=0.02292), Clade(branch_length=0.00218), Clade(branch_length=0.14735), Clade(branch_length=0.07756), Clade(branch_length=0.23141), Clade(branch_length=0.11327)]

(아래 코드 결과)

 

print(tree.find_clades()) # clade를 찾는다.
print(tree.find_elements()) # element를 찾는다.
print(tree.find_any()) # 다 찾는다.
# find_any()빼면 다 filter object 나온다...
print(tree.common_ancestor("AB021386.1.1495")) # 조상님좀 찾아주세요
print(tree.count_terminals()) # Terminal 갯수좀 세주세요
print(tree.depths()) # depth 얼마임?
print(tree.distance("AB021386.1.1495","U20163.1.1495")) # 얘 거리 얼마냐? (두개 찾아주는 것도 된다)
print(tree.total_branch_length()) # 전체 branch 길이
# 여기까지는 수치로 나오는 정보
print(tree.is_bifurcating()) 
# True if the tree is strictly bifurcating (아마 트리가 둘로만 갈리면 True인 듯)
print(tree.is_monophyletic("CP010835.91794.93350")) 
# Test if all of the given targets comprise a complete subclade(
# 단일 계통 여부)
print(tree.is_parent_of("CP010835.91794.93350")) 
# 얘가 이 트리에서 후손관계야?
print(tree.is_preterminal()) 
# True if all direct descendents are terminal(직계 후손이 terminal이면 참이라는데 뭔 소리지)
# 여기까지는 boolean

설명은 주석을 참고하고... 위 코드들을 입력하면 트리(계통수)에 대한 정보를 출력해준다. 

 

수정하기

아니 근데 실질적인 수정은 안가르쳐줌... 라벨 좀 바꿉시다... 

 

(Before)

 

Collapse

from Bio import Phylo
tree = Phylo.read("/home/koreanraichu/Halobacterium.ph", "newick")
tree.name="Halophile"
tree.collapse("CEML01000001.1297981.1299440")
Phylo.draw(tree)

Collapse를 주면 있었던 게 없어진다. (잘 보면 collapse 코드 안에 있는 항목이 계통수에서 없어졌다) 있었는데요 없어졌습니다 

 

근데 collapse_all을 주면 트리가 빗이 된다. 뭐냐 이건 

 

ladderize

clade를 터미널 노드 수에 따라 정렬해준다. 

 

prune

푸룬? 자두? 아니고... 찾아보니까 가지치기라는 뜻이 있다. 근데 뭔 차이인지는 모르겠다. 

 

root_with_outgroup

tree.root_with_outgroup("CEML01000001.1297981.1299440")

 

tree.root_with_outgroup("AB105159.1.1501")

아웃그룹으로 빼버릴 애를 지정하면 걔를 밖으로 빼준다. 근데 왜 밑으로 빼냐. 

 

root_at_midpoint

제일 중구난방인 걸 가운데로 빼버린 듯... 

 

split

밑에 후손을 더 만들어준다. 기본값은 2. 

 

일단 그려보자

from Bio import Phylo
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord
# FASTA file reading
records = SeqIO.parse("/home/koreanraichu/Halobacterium.fasta","fasta")
save_record=[]
for record in records:
    description = record.description.split(sep=";") # record description 소환
    record.id = description[6] # description 마지막 글자로 변환
# fasta 소환!
    record_save=SeqRecord(record.seq,id=record.id,name=record.id,description=record.description)
    save_record.append(record_save)
print(save_record)
SeqIO.write(save_record,"/home/koreanraichu/Halobacterium_label.fasta","fasta")
# FASTA 파일로 저장은 했는데 description에 세미콜론이 들어가서 안된다...
# 결국 raw file 들어가서 수정함.
# 아 그리고 학명이 띄어쓰기가 되어 있으니까 앞에꺼 빼고 뒤에꺼 넣길래 raw file에서 언더바 추가했음.

tree = Phylo.read("/home/koreanraichu/Halobacterium_label.ph", "newick")
tree.name="Halophile" # Title
tree.as_phyloxml() # phyloxml로 변환
Phylo.draw(tree)

결과

지금까지 그렸던 phylogenic tree를 보면 라벨이 AB 어쩌고 이런 식으로 되어 있다. 이게 Silva에서 rRNA 시퀀스 받아오면 FASTA 형식으로 저장될 때 그게 아이디가 되다 보니 clustalW나 MUSCLE로 트리 그리고 나서도 아이디가 뜬다. 근데 이거 안이쁨. 뭔지 모름. 그래서 학명으로 바꿔줘야 한다… 

그래서 소스코드가 좀 중구난방이 되었으나.. 크게 코드의 역할로 나눠보자면 

1. Silva에서 받은 FASTA파일을 열어서 재구축한 후 저장
2. 아 트리 그려달라고 

이렇게 나뉜다. 

재구축은 ID를 학명으로 대체한 다음(description을 split한 다음 맨 끝에 있는 학명을 ID로 대체) 저장하는 구조이고, 저장한 FASTA file의 경우 ClustalW에서는 안 돌아가고 MUSCLE에서는 돌아가는데 세미콜론때문에 트리 그리다 오류나서 raw file단에서 세미콜론을 다 지웠다. 저장하는것도 애먹었지... 레코드는 다 추출되는데 저장이 하나만 됨...OTL 

Lv. 35 라이츄

Lv. 35 라이츄

광고 매크로 없는 청정한 블로그를 위해 노력중입니다. 근데 나만 노력하는 것 같음… ㅡㅡ

Biopython으로 PDB 탐방하기

Coding/Python 2022. 8. 21. 00:33

PDB!! PDB!! 


들어가기 전에

PDB는 protein data bank의 약어이기도 하고, 거기서 제공하는 파일 형식이기도 하다. 여기서는 그냥 PDB라고 하면 데이터뱅크, PDB '파일'이라고 하면 PDB 파일이다. 

그리고 이새기들 쿡북쓰기 귀찮았는지 모듈 불러와야 하는 거 빼먹더라... 니들도 일하기 싫었구나 

 

파일 읽기

쓰기도 있긴 한데 그건 생략. 읽는것도 하난가 두갠가 오류나서 안됐다. 이 섹션에서 읽을 파일은 

1) mmCIF
2) MMTF
3) PDB파일
4) PQR

이다. 

 

mmCIF

from Bio.PDB.MMCIFParser import MMCIFParser
parser = MMCIFParser()
structure = parser.get_structure("7f0l", "/home/koreanraichu/7f0l.cif")
print(structure)
from Bio.PDB.MMCIF2Dict import MMCIF2Dict
mmcif_dict = MMCIF2Dict("/home/koreanraichu/7f0l.cif")
print(mmcif_dict)
# FileNotFoundError: [Errno 2] No such file or directory: '1FAT.cif'
# 위 에러가 뜨는 경우, 해당하는 확장자의 파일이 내 PC에 있어야 한다.

둘 중 편한 방식으로 불러오면 된다. 참고로 해당 파일이 컴퓨터에 있어야 불러올 수 있다. 다운받는거는 밑에 나옴. 

 

MMTF

from Bio.PDB.mmtf import MMTFParser
structure = MMTFParser.get_structure_from_url("4CUP")
print(structure)
# Bio.MissingPythonDependencyError: Install mmtf to use Bio.PDB.mmtf (e.g. pip install mmtf-python)
# 뭔데 이거
from Bio.PDB.mmtf import MMTFParser
structure = MMTFParser.get_structure_from_url("4CUP")
print(structure)
# Bio.MissingPythonDependencyError: Install mmtf to use Bio.PDB.mmtf (e.g. pip install mmtf-python)

개발자 나와봐 나랑 면담 좀 하자. 저건 뭔 거지같은 오류냐. 

 

PDB파일

from Bio.PDB.PDBParser import PDBParser
parser = PDBParser(PERMISSIVE=1)
structure_id = "2b10"
filename = "/home/koreanraichu/2b10.pdb"
structure = parser.get_structure(structure_id, filename)
print(structure)
from Bio.PDB import parse_pdb_header
with open(filename, "r") as handle:
    header_dict = parse_pdb_header(handle)
    print(header_dict)

일단 둘 다 되긴 되는데 뭔 Warning이 이렇게 뜨냐... 이럴거면 PDB 가서 봤지... 

 

PQR

from Bio.PDB.PDBParser import PDBParser
pqr_parser = PDBParser(PERMISSIVE=1, is_pqr=True)
structure_id = "7lfv"
filename = "/home/koreanraichu/7lfv.pqr"
structure = parser.get_structure(structure_id, filename, is_pqr=True)
print(structure)
# NameError: name 'parser' is not defined

예제 코드는 이렇게 되어 있는데, 저렇게 쓰면 오류난다. 위에는 pqr_parser로 불러놓고 아래에서 parser쓰면 컴퓨터는 이게 뭔데요 하지 찰떡같이 알아듣질 않음. 에이 이건 알파고도 못알아먹어... 

from Bio.PDB.PDBParser import PDBParser
pqr_parser = PDBParser(PERMISSIVE=1, is_pqr=True)
structure_id = "7lfv"
filename = "/home/koreanraichu/7lfv.pqr"
structure = pqr_parser.get_structure(structure_id, filename)
print(structure)

이렇게 하면 경고는 많이 뜨지만 일단 된다. 

 

쓰기

그냥 이렇게 쓰는구나만 보고 가자. 

 

mmCIF

io = MMCIFIO()
io.set_structure(s)
io.save("out.cif")

io = MMCIFIO()
io.set_dict(d)
io.save("out.cif")

 

MMTF

from Bio.PDB.mmtf import MMTFIO
io = MMTFIO()
io.set_structure(s)
io.save("out.mmtf")

 

PDB

io = PDBIO()
io.set_structure(s)
io.save("out.pdb")
class GlySelect(Select):
    def accept_residue(self, residue):
        if residue.get_name()=="GLY":
            return True
        else:
            return False
io = PDBIO()
io.set_structure(s)
io.save("gly_only.pdb", GlySelect())

밑에껀 글라이신만 매긴건가? ㄷㄷ 

 

PQR

io = PDBIO(is_pqr=True)
io.set_structure(s)
io.save("out.pdb")

 

Data representation

PDB에서 제공하는 데이터는 카테고리가 있다. CSS의 클래스 차일드 이런것처럼. 상위 개념부터 순서대로 

 

Structure>Model>Chain>Residue>Atom

 

으로 간다. 보통 model의 경우 NMR로 규명한거는 겁나 많고, X-ray crystallography로 규명한거는 하나. 

 

Model

Structure는 꺼내고 자시고 할 것도 없으니 모델부터 봅시다. 

 

from Bio.PDB.PDBParser import PDBParser
parser = PDBParser(PERMISSIVE=1)
structure_id = "5zi1"
filename = "/home/koreanraichu/5zi1.pdb"
structure = parser.get_structure(structure_id, filename)
# 뭘 불러와야 하죠

print(structure[0])
<Model id=0>

모델은 ID가 int 형태라 저렇게 불러오면 된다. 하위 카테고리고? 그건 또 따로 불러야된다. 

 

Chain

모델의 하위 카테고리. 얘는 ID가 A, B, C, D 이런 식으로 string 형태이다. 

 

model=structure[0]
print(model["A"])
<Chain id=A>

여기서 짐작하신 분들도 계시겠지만, 리스트 전체 뽑을거면 iteration 들어가야 한다. 그것도 밑에서 서술할 예정. 

 

Residue

얘부터 이제 뭔가 겁나 많아질 예정이다. 일단 Residue 자체는 튜플이다. 낙장불입 구성 요소는 

-Hetero-field-W(물분자)/H(다른 HET 영역에 기재된 분자들-예: Heme)/공백(standard molecule)
-Sequence identifier: 해당 residue의 시퀀스 내 위치 정보를 담은 정수
-Insertion code: 잘은 모르겠지만 insertion mutant 관련된 코드인 듯 하다.

이렇게 세 개. 

 

model=structure[0] # Model
chain=model["A"] # Chain
print(chain[(" ", 100, " ")])
<Residue PHE het=  resseq=100 icode= >

보통 저렇게 치기도 하고 숏컷으로 chain[100] 이런 식으로 치기도 한다. 

 

print(residue.get_resname())
print(residue.is_disordered())
print(residue.get_segid())
print(residue.has_id(10))
PHE
0
    
False

Residue에서 뽑아낼 수 있는 속성은 이름/Disordered 여부(1이면 Disordered)/segment ID/얘 아이디 있냐? 정도. 참고로 이 밑에 가면 정신없어서 내가 주석으로 달아놨다. 

 

Atom

# 이제 atom으로 가보자...
atom=residue["CA"]
print(atom.get_name()) # 이름!
print(atom.get_fullname()) # 이름! (풀네임)
print(atom.get_id()) # ID
print(atom.get_coord()) # 좌투더표
print(atom.get_vector()) # 좌표를 벡터 객체로 받음
print(atom.get_bfactor()) # Isotropic B factor
print(atom.get_occupancy()) # Occupancy
print(atom.get_altloc()) # alternative location specifier
print(atom.get_sigatm()) # 원자 매개변수의 표준편차
print(atom.get_siguij()) # Isotropic B factor의 표준편차
print(atom.get_anisou()) # anisotropic B factor
CA
 CA 
CA
[ 20.149 -36.25   32.947]
<Vector 20.15, -36.25, 32.95>
46.05
1.0
 
None
None
None

Process finished with exit code 0

뭔가 많으니 주석을 참고할 것. 

 

# 아 일일이 빼기 귀찮으시죠? 그럼 한줄로 ㄱㄱ하세요.
atom = structure[0]["A"][100]["CA"]
print(atom.get_name())

이런 식으로 한줄로 빼버릴수도 있다. 복잡해서 글치. 

 

Disorder

Point mutation에 대한 정보라는데 조회 안되던데? 

 

Disordered Atoms

from Bio.PDB.PDBParser import PDBParser
parser = PDBParser(PERMISSIVE=1)
structure_id = "4lqm"
filename = "/home/koreanraichu/4lqm.pdb"
structure = parser.get_structure(structure_id, filename)
# EGFR L858R을 불러보자
model=structure[0]
chain=model['A']
print(chain[858])
<Residue ARG het=  resseq=858 icode= >

일단 시범조교 L858R을 불러보자. (아르기닌 한글자 약어가 R이다)

 

atom.disordered_select("A")
print(atom.get_altloc())

근데 이거 써도 오류남... 

 

Disordered Residue

residue = chain[10]
residue.disordered_select("CYS")

왜 Attribute error가 뜨는지 1000자 이내로 서술좀 해주실까. 

 

Iteration 활용하기

다운로드? 그거 맨 밑에요. 

 

전체 보기

from Bio.PDB.PDBParser import PDBParser
parser = PDBParser()
structure = parser.get_structure("2b10", "/home/koreanraichu/2b10.pdb")
model = structure[0]
chain = model["A"]
residue=chain[10]
atom=residue["CB"]
print(atom)
# 이것이 각개전투다 절망편(...)

이것이 기존 방법이다. 근데 이렇게 보면 어때요? 귀찮아요. 그리고 파일을 까 보지 않는 이상 100번째 아미노산이 뭔지 모른다. 

 

from Bio.PDB.PDBParser import PDBParser
p = PDBParser()
structure = p.get_structure("2m3g", "/home/koreanraichu/2m3g.pdb")
for model in structure:
    for chain in structure:
        for residue in structure:
            for atom in structure:
                print(atom)
# iteration 가즈아!!!

일단 가즈아!!! 

 

structure2 = p.get_structure("2b10", "/home/koreanraichu/2b10.pdb")
atoms = structure2.get_atoms()
for atom in atoms:
    print(atom)
# 단축 가즈아!!!

위에 코드도 길고 번거롭다면 이런 방법도 있다. (아예 구조에서 다이렉트로 아톰 뽑아서 가져오기)

 

atoms = chain.get_atoms()
for atom in atoms:
    print(atom)
# chain에서 가져올 수도 있는 듯 하다.

이런 식으로 상위 카테고리에서 get_ 걸어도 된다. 

 

residues = model.get_residues()
for residue in residues:
    print(residue)
chains=structure2.get_chains()
for chain in chains:
    print(chain)
# chain도 된다.
models=structure2.get_models()
for model in models:
    print(model)
# 모델은 이걸로도 된다

structure단에서 get_뒷부분 바꿔서 가져와도 된다. 

 

취사선택하기

아 사람이 큰 그림만 볼 수는 없지... 또 조건부로 볼 때가 있어요 또. 

 

from Bio.PDB.PDBParser import PDBParser
p = PDBParser()
structure = p.get_structure("2b10", "/home/koreanraichu/2b10.pdb")
# Cytochrome C peroxidase

일단 시범조교 앞으로. 

 

for model in structure.get_list():
    for chain in model.get_list():
        for residue in chain.get_list():
            if residue.has_id("CA"):
                ca = residue["CA"]
                if ca.get_bfactor() > 50.0:
                    print(ca.get_coord())
# List 내 알파탄소 중 B factor가 50 이상인 것의 '좌표'

Residue 중 알파 탄소(아미노산의 몸통)의 B factor가 50 이상인 것의 좌표를 출력한다. 

 

for residue in chain.get_list():
    residue_id = residue.get_id()
    hetfield = residue_id[0]
    if hetfield[0] == "H":
        print(residue_id)
# HET residue를 전체 출력하는 코드

HET residue를 출력한다. 응? 저 ID? 위에 튜플로 된 그거 말하는거다. 튜플도 수정만 안 될 뿐이지 일단 인덱싱은 된다. 순서대로 W or H or 공백/sequence identifier/insertion code니까 가장 왼쪽(0번 인덱스)에 H가 들어간 것을 찾으면 HET residue. W는 물분자고 공백이면 단백질을 구성하는 아미노산이다. 

 

for model in structure.get_list():
    for chain in model.get_list():
        for residue in chain.get_list():
            if residue.is_disordered():
                resseq = residue.get_id()[1]
                resname = residue.get_resname()
                model_id = model.get_id()
                chain_id = chain.get_id()
                print(model_id, chain_id, resname, resseq)
# Disordered atom을 포함하는 모든 residue들을 출력한다. 참고로 여기에는 없다.

Disordered Atom을 포함하는 모든 residue들을 출력하는 코드. 

 

for model in structure.get_list():
    for chain in model.get_list():
        for residue in chain.get_list():
            if residue.get_resname() == "LEU":
                print(residue.get_id())
# 응용: Residue 중 류신 뽑아줘(Leu)

Residue 중 류신의 ID를 뽑는다. (그니까 몇번 몇번이 류신이냐 이 얘기)

 

polypeptide 관련

이거 하려면 뭐 또 불러와야되는데 이놈들 그거 기재 안했음... 

 

from Bio.PDB.PDBParser import PDBParser
from Bio.PDB.Polypeptide import PPBuilder
ppb=PPBuilder()
# 예제에는 안나왔는데 이거 두 줄 추가해야된다...
p = PDBParser()
structure = p.get_structure("7kz1", "/home/koreanraichu/7kz1.pdb")
model_nr = 1
polypeptide_list = ppb.build_peptides(structure, model_nr)
for polypeptide in polypeptide_list:
    print(polypeptide.get_sequence())
KKWTPPRSPFNLVQETLFHDPWKLLIATIFLNRTSGKMAIPVLWKFLEKYPSAEVARTADWRDVSELLKPLGLYDLRAKTIVKFSDEYLTKQWKYPIELHGIGKYGNDSYRIFCVNEWKQVHPEDHKLNKYHDWLWENHEKLS

이거 하려면 모듈을 불러야 하는데 쿡북에 빠졌다. 이녀석들 보게. 

 

Analyze

이거 항목은 많은데 거리 각도 이면각만 보고 끝냅시다... 뭐 깔라는데 안깔려요 그게... 

아, 그리고 각도랑 이면각은 

from Bio.PDB.vectors import calc_angle
from Bio.PDB.vectors import calc_dihedral

이거 있어야된다. 이새기들 이것도 빼먹음. 

 

거리

from Bio.PDB.PDBParser import PDBParser
p = PDBParser()
structure = p.get_structure("5ngo", "/home/koreanraichu/5ngo.pdb")
# 아무튼 불러와보겠습니다.
residues=structure.get_residues()
for residue in residues:
    print(residue)
# 나도 이거 처음보는거라 목록부터 봐야됨...

일단 불러놓고 생각하는 타입. 

residue1=structure[0]["A"][273]["CA"]
residue2=structure[0]["A"][278]["CA"]
print(residue1-residue2)
# Residue간 거리

거리는 뭐 재고 할 것도 없이 간단하다. 

atom1=structure[0]["A"][423]["CA"]
atom2=structure[0]["A"][423]["CB"]
atom3=structure[0]["A"][423]["N"]
vector1=atom1.get_vector()
vector2=atom2.get_vector()
vector3=atom3.get_vector()
angle = calc_angle(vector1, vector2, vector3)
print(angle)
# 각(별)도
# 0.5872530070961592
atom1=structure[0]["A"][415]["CA"]
atom2=structure[0]["A"][423]["CA"]
atom3=structure[0]["A"][431]["CA"]
atom4=structure[0]["A"][439]["CA"]
vector1 = atom1.get_vector()
vector2 = atom2.get_vector()
vector3 = atom3.get_vector()
vector4 = atom4.get_vector()
torsion = calc_dihedral(vector1,vector2,vector3,vector4)
print(torsion)
# Torsion angle 계산 
# 1.28386400091194

각도와 torsion angle은 벡터로 받아서 계산해야 한다. 

 

Internal coordinates for standard residues

적절한 번역 아시는 분 제보 바랍니다. 

 

model=structure[0]
model.atom_to_internal_coordinates()
for r in model.get_residues():
    if r.internal_coord:
        print(
            r,
            r.internal_coord.get_angle("psi"),
            r.internal_coord.get_angle("phi"),
            r.internal_coord.get_angle("omega"),  # or "omg"
            r.internal_coord.get_angle("chi2"),
            r.internal_coord.get_angle("CB:CA:C"),
            (r.internal_coord.get_length("-1C:0N")  # i-1 to i peptide bond
              if r.internal_coord.rprev else "None")
        )
<Residue SER het=  resseq=451 icode= > None 85.88349606675206 167.36084855239682 None 116.39783860086264 None

 

PDB 접속하기

얘네 이것도 모듈 불러야되는데 안썼더라. 

 

from Bio.PDB.PDBList import PDBList # 이거 불러야된다
pdbl = PDBList()
pdbl.retrieve_pdb_file("6WO1")

이렇게 받으면

이렇게 저장된다. (기본 형식은 mmCIF)

 

from Bio.PDB.PDBList import PDBList
pdbl = PDBList()
pdbl.retrieve_pdb_file("6WO1",file_format="pdb")
# 확장자를 따로 입력하지 않으면 CIF파일로 다운로드 된다.
# file_format="확장자"로 입력하면 특정 파일 형식으로 받을 수 있다.
# pdir="경로"를 입력하면 다운로드 경로도 정할 수 있다.

이렇게 file_format으로 지정해주면 된다. 

 

뭐야 PDB파일 돌려줘요 

혹시 경로지정 되냐고? pdir="경로" 하면 된다. 

Lv. 35 라이츄

Lv. 35 라이츄

광고 매크로 없는 청정한 블로그를 위해 노력중입니다. 근데 나만 노력하는 것 같음… ㅡㅡ

Biopython으로 Swiss-prot과 ExPASy 데이터베이스 탐방하기

Coding/Python 2022. 8. 21. 00:14

근데 스위스프롯은 긁어와서 저장 안됨? 


첫빠따는 파싱이 국룰이지

파싱 방법이 네 가지가 있는데 gZIP은 생략. gZIP 파일을 못 구했다. 

 

handle=open('/home/koreanraichu/Q63HQ2.txt')
print(handle)
import gzip
handle = gzip.open("myswissprotfile.dat.gz", "rt")
from urllib.request import urlopen
url = "https://raw.githubusercontent.com/biopython/biopython/master/Tests/SwissProt/F2CXE6.txt"
handle = urlopen(url)
print(handle)
from Bio import ExPASy
handle = ExPASy.get_sprot_raw("Q63HQ2")
print(handle)

위에서부터 순서대로 

1. 갖고 있는 파일 열기 
2. gZIP 핸들링해서 열기(...)
3. URL 입력해서 갖고오기
4. ExPASy에서 가져오기

이상이다. 

네? 출력이요? 그냥 파싱에 포문 빠져서 그런갑다 하자... 

 

실제로 가져와봤습니다

from Bio import ExPASy
from Bio import SwissProt
# 일단 여기서는 4번, ExPASY에 접속해서 가져와보자.
handle = ExPASy.get_sprot_raw("Q63HQ2")
record = SwissProt.read(handle)
print(record.description)
RecName: Full=Pikachurin; AltName: Full=Agrin-like protein; AltName: Full=EGF-like, fibronectin type-III and laminin G-like domain-containing protein; Flags: Precursor;

참고로 위 방법 중 4번을 선택했다. 네? 결과에 피카츄요? 내가 뭘 잘못 본 거 아니냐고요? 정상인데요? 단백질 이름이 피카츄린이니까 피카츄가 들어가지. 

 

for ref in record.references:
    print("authors:", ref.authors)
    print("title:", ref.title)
    # 아마도 관련 참고문헌 정보?

이런것도 된다. 

 

Prosite 레코드 파싱하기

Prosite는 단백질 도메인, 단백질에서 기능하는 부분과 인식 프로파일 정보, functional family에 대한 정보... 가 들어가 있는데... 이거 CATH 아니냐? CATH와는 다르다 CATH와는! 

 

ftp://ftp.expasy.org/databases/prosite/prosite.dat 

dat파일은 여기서 받자. 리눅스 유저는 터미널에서 wget -i 들어가면 되고... 윈도우는...... 여기다 올려주려고 했는데 10MB 넘어가서 안됨... 리눅스로 갈아타자 뭔소리여 들어가면 텍스트 뜨는데 걍 메모장에 복붙하고 prosite.dat으로 이름 지어주면 된다. 깃헙에요? 저게 올라감? ㄷㄷ 

 

from Bio.ExPASy import Prosite
handle = open("/home/koreanraichu/prosite.dat")
records = Prosite.parse(handle)
print(records)

이렇게 가져오면 된다는데 결과가 또 뭔가 괴랄하다. (dat파일 받는 링크는 위에 올려드림) 

from Bio.ExPASy import Prosite
from urllib.request import urlopen
handle=open("/home/koreanraichu/prosite.dat")
records = Prosite.parse(handle)
record = next(records)
print(record.accession)

이것까지 주면 원하는 정보만 볼 수는 있는데... 어? next? 

 

from Bio.ExPASy import Prosite
from urllib.request import urlopen
handle=open("/home/koreanraichu/prosite.dat")
records = Prosite.parse(handle)
for i in range(0,5):
    record = next(records)
    print(record.accession)

For문 마려우셨죠? 됩니다. 

from Bio.ExPASy import Prosite
from urllib.request import urlopen
handle=open("/home/koreanraichu/prosite.dat")
records = Prosite.parse(handle)
i=0
while i < 5:
    record = next(records)
    print(record.accession)
    i=i+1

아 For문 되면 while도 됨 

저 코드 각개로 해보는 걸 추천한다. 연달아서 하면 1-(2~6)-(7~11) 이렇게 나온다. 

 

n=0
for record in records: n+=1
print(n)

이 코드로 레코드 숫자도 볼 수 있다. 

 

효소 레코드 파싱하기

from Bio.ExPASy import Enzyme
with open("/home/koreanraichu/RuBisCO.txt") as handle:
    record = Enzyme.read(handle)
    print(record['ID'])
4.1.1.39

위 번호는 EC번호로, 효소가 반응을 '어떤 방식으로' 촉매하는가에 따라 번호가 붙는다. 

 

from Bio.ExPASy import Enzyme
with open("/home/koreanraichu/RuBisCO.txt") as handle:
    record = Enzyme.read(handle)
    print(record['ID']) # EC no.
    print(record['DE']) # description
    print(record['AN']) # 대충 synonyms같은건가? 뭐 얘 이렇게도 불러요 이런거
    print(record["CA"]) # 촉매하는 반응(오 이거 식으로 나온다)
    print(record["PR"]) # 이건 모르겠다... 데이터베이스 번호인가... 
    print(record["CC"]) # 아마도 뭐 하는 효소인가에 대한 설명인 듯 
    print(record['DR']) # 뭔진 모르겠지만 일단 잘못했어요... 뭐가 되게 많이떴는데 넘파이 마려웠음

이거 다 띄웠더니 뭐가 많아서 결과는 생략. 작동 하기는 한다. 

 

ExPASy 서버 접속해보기

이제 본격적으로 서버 털러 가 봅시다. 일단 ExPASy 관련해서는 네 가지 쿼리가 있는데 

get_prodoc_entry: To download a Prosite documentation record in HTML format
get_prosite_entry: To download a Prosite record in HTML format
get_prosite_raw: To download a Prosite or Prosite documentation record in raw format
get_sprot_raw: To download a Swiss-Prot record in raw format

이렇게 네 개가 있다. 

다운로드래서 행복회로 태웠는데 저장 안되더만 swissprot은... 

 

Swissprot 검색하기

참고로 이거 Uniprot 번호 먹힌다. 

 

from Bio import ExPASy
from Bio import SwissProt
accessions = ["O23729", "O23730", "O23731"]
records = []
# 저 세개에 대해 찾고싶은갑다.
for accession in accessions:
    handle = ExPASy.get_sprot_raw(accession)
    record = SwissProt.read(handle)
    records.append(record)
    print(records) # IDLE 아니라 이거 있어야 나옴
[<Bio.SwissProt.Record object at 0x7f1abf171f70>, <Bio.SwissProt.Record object at 0x7f1abf9a3d00>, <Bio.SwissProt.Record object at 0x7f1abf117550>]

이러시는 이유가 있으실 것 아니예요. 

 

from Bio import ExPASy
from Bio import SwissProt
accessions = ["C3RT18"] # 단백질 이름은 신경쓰지 맙시다 
records = []
# 저 세개에 대해 찾고싶은갑다.
for accession in accessions:
    handle = ExPASy.get_sprot_raw(accession)
    record = SwissProt.read(handle)
    records.append(record)
print(records) #IDLE 아니라서 이거 써야 나온다

단식 검색도 된다. 

 

>>> for accession in accessions:
...     handle = ExPASy.get_sprot_raw(accession)
...     try:
...         record = SwissProt.read(handle)
...     except ValueException:
...         print("WARNING: Accession %s not found" % accession)
...     records.append(record)

검색했는데 결과가 없을 때 ValueError가 여러분을 반길 것이다. 그럴때는 당황하지 말고 에외 처리를 해 보자. 

 

Prosite 검색하기

from Bio import ExPASy
handle = ExPASy.get_prosite_raw("PS00001")
text = handle.read()
print(text)
# 가장 원시적인 긁는 방법
ID   ASN_GLYCOSYLATION; PATTERN.
AC   PS00001;
DT   01-APR-1990 CREATED; 01-APR-1990 DATA UPDATE; 01-APR-1990 INFO UPDATE.
DE   N-glycosylation site.
PA   N-{P}-[ST]-{P}.
CC   /SITE=1,carbohydrate;
CC   /SKIP-FLAG=TRUE;
CC   /VERSION=1;
PR   PRU00498;
DO   PDOC00001;
//

가장 원시적인 긁어오기 신공. 아 근데 개구린 하려고 했는데 prosite에는 개구린이 없다... (개구린: 단백질 이름) 

 

from Bio import Prosite
handle = ExPASy.get_prosite_raw("PS51036")
record = Prosite.read(handle)
print(record)

이건 된다고 올려놨는데 첫 줄에서 ImportError: cannot import name 'Prosite' from 'Bio' (/home/koreanraichu/PycharmProjects/pythonProject/venv/lib/python3.8/site-packages/Bio/init.py가 나를 반긴다. 넌 또 왜그러냐. 

 

from Bio.ExPASy import Prodoc
handle = ExPASy.get_prosite_raw("PDOC00001")
record = Prodoc.read(handle)
print(record)

Prodoc을 가져올 수도 있다. 

 

handle = ExPASy.get_prosite_entry("PS00001")
html = handle.read()
with open("myprositerecord.html", "w") as out_handle:
    out_handle.write(html)
# HTML format으로 다운로드 받을 수 있다.

handle = ExPASy.get_prodoc_entry("PDOC51036")
html = handle.read()
with open("myprositerecord2.html", "w") as out_handle:
    out_handle.write(html)
# 얘는 prodoc 다운로드 하는 코드

HTML 형식으로 저장도 된다. 

이런 식으로 생성되고(경로 지정을 안 했더니 파이참 폴더에 박아버리네...)

위는 prosite, 아래는 prodoc. 이런 식으로 뜬다. 

 

Prosite database 스캔하기

스캔 아니고 시퀀스 검색 아니냐. 

 

from Bio.ExPASy import ScanProsite
seq="MEHKEVVLLLLLFLKSGQGEPLDDYVNTQGASLFSVTKKQLGAGSIEECAAKCEEDEEFTCRAFQYHSKEQQCVIMAENRKSSIIIRMRDVVLFEKKVYLSECKTGNGKNYRGTMSKTKN"
handle = ScanProsite.scan(seq=seq)
result = ScanProsite.read(handle)
print(result)
[{'sequence_ac': 'USERSEQ1', 'start': 16, 'stop': 98, 'signature_ac': 'PS50948', 'score': '8.873', 'level': '0'}]

어? 딕셔너리네? 

 

print(result.n_seq)
print(result.n_match)
print(result[0])
1
1
{'sequence_ac': 'USERSEQ1', 'start': 16, 'stop': 98, 'signature_ac': 'PS50948', 'score': '8.873', 'level': '0'}

뭐 이렇다는데... 이거 예제꺼 긁어왔는데 왜 결과가 다르냐... 

Lv. 35 라이츄

Lv. 35 라이츄

광고 매크로 없는 청정한 블로그를 위해 노력중입니다. 근데 나만 노력하는 것 같음… ㅡㅡ

방명록