Notes: 파이썬 모듈이나 버전은 따라해도 되는데 OS는 본인 편한거 쓰는게 좋음. 개발하시는 분들이 맥이나 리눅스 많이 쓰긴 한데 맥이나 리눅스는 윈도우와 달라서 무턱대고 쓰기 좀 불편하다. 나 고딩때 애들 커맨드 컨트롤 헷갈리는것만 여러번 봤음 진짜. (맥에서는 복사가 Ctrl+C가 아니라 Command+C) 그리고 어지간한 개발툴은 윈도 맥 리눅스 다 지원한다. 카톡 빼고그건 개발툴 아닌데아 그래서 지원 안하나
내가 리눅스 쓸거다 근데 특히 우분투 쓸거다 그러면 우분투에서 여기다 우리꺼 깔면 잘됨 보장함 땅땅땅 있는 노트북 중에서 사는게 좋음.
일단 본인은 노트북 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은 바로가기가 있는데 왜 파이참만...)
아 그리고 바탕화면 그림은 바꿀 수 있는데 윈도처럼 바탕화면에 따라 작업표시줄 색 바꾸는건 안됨. 리누스씨 이 기능 빨리 넣어줘요
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가 편하다보니 그냥 글 쓸 때만 쓴다.
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로 올려야돼서 귀차니즘이 배가됩니다. 간단한거(별찍기)는 이미지 많이 안 올려도 돼서 워드프레스에 올리지만...
Kyoto Encyclopedia of Genes and Genomes의 준말. 그렇다, 이름에 교토가 들어간 걸 보면 아시겠지만 일제 DB다.
여기가 메인페이지
여기가 KEGG brite. 본인이 자주 가는 곳이다. 생각보다 쏠쏠한 정보가 많고 KEGG brite의 경우 golden standard dataset이라고 해서 야먀니시가 인공지능 학습시킬 때 쓴 데이터셋(GPCR, 효소, 핵 리셉터, 트랜스포터) 분류별로 약물 타겟을 알려준다.
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만 보는 코드
다 좋은데 그걸 꼭 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을 찾아준다)
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))
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에서 사용자가 입력한 항목을 찾아준다. (와일드카드는 안 먹는다)
이것도 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()
# 출력
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
주석에는 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)로 수정해줘야 제대로 나온다.
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)
얘는 배열이 아니라 거리 행렬을 입력으로 받는다. 아까 그 괴랄한 형태도 되고, 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)
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)
주성분 분석. 고차원의 데이터를 저차원으로 축소시키는 기법 중 하나라고 한다. 아니 무슨 선형대수학이 왜 나오는거냐 대체... 와 진짜 대환장 파티네 뭐 주성분의 표본을 직교변환한다던가... 이런건 나도 모르겠고... 저는 수학이랑 담 쌓아서요... 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)
이게 한 글에 다 쓰기엔 좀 분량도 분량인데 이게 생각보다 설명이랑 코딩이랑 나뉘어있어서 이론편 실전편 나눕니다.
이건 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이라고요? 그러면 둘이 상관 없는거다.
모티브 찾으면 구조 모티브가 나오는데 이거는 '단백질이나 핵산과 같은 사슬 모양의 생체 분자에서 진화적으로 관련이 없는 다양한 분자에서 나타나는 일반적인 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)
염기별로 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))
# 찾았음!
-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)
# 일단 뽑고 보는 타입
유전자나 단백질 시퀀스 분석(균이라면 16s rRNA라던가)을 통해 얘네가 얼마나 가까운지를 알아내게 되면 그걸 저런 식으로 그려서 나타내는 것. 저렇게 생물 종에 따라 그리는 경우도 있고, 특정 단백질의 homolog나 다른 생물종에서 같은 역할을 하는 단백질에 대해서 저걸 그리기도 한다.
와! 계통수! 그려보자!
from Bio import Phylo
tree = Phylo.read("/home/koreanraichu/Deinococcus.ph", "newick")
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의 끝 노드만, 아래 노드는 중간 노드까지 출력해주는데 왜 위쪽 노드 결과값이 더 긴 지는 묻지 말자.
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
푸룬? 자두? 아니고... 찾아보니까 가지치기라는 뜻이 있다. 근데 뭔 차이인지는 모르겠다.
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
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)
모델은 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 관련된 코드인 듯 하다.
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())
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")
)
from Bio.PDB.PDBList import PDBList
pdbl = PDBList()
pdbl.retrieve_pdb_file("6WO1",file_format="pdb")
# 확장자를 따로 입력하지 않으면 CIF파일로 다운로드 된다.
# file_format="확장자"로 입력하면 특정 파일 형식으로 받을 수 있다.
# pdir="경로"를 입력하면 다운로드 경로도 정할 수 있다.
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)
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)
보통 Biopython을 쓰거나 랜덤, 넘파이, 판다스를 쓸 때는 뭘 모셔와야 하는데, Entrez에 접속하는 모듈도 마찬가지다. 근데 바이오파이썬은 그걸 떠나서 모셔오는 게 너무 핵가족 스케일이여. 아무튼... 그래서 이번에는
from Bio import Entrez
이걸 필두로 뭘 많이 모셔올 예정인데... 아니 아직 아냐 마저 듣고 가.
Entrez에 접속해서 뭘 하려면 저거 말고 필수적으로 입력해야 하는 게 있다. 1. 너님의 API 키 2. 너님의 메일 주소 3. 너님의 매개 변수
셋 중 하나는 반드시 입력해야 하고, 여기서는 이메일을 입력할건데 저거 뭐 이메일 제출한다고 CIA에서 당신 털러 오는 거 아니니까 안심하고 쓰자. 오히려 저거 안 쓰면 에러난다.
Einfo
Entrez 자체에 대한 정보를 볼 수 있다.
from Bio import Entrez
Entrez.email = "blackholekun@gmail.com" # 내가 누구인지 말해주는 과정이 필요하다고...
# 자기 이메일 쓰면 됩니다
handle = Entrez.einfo()
result = handle.read()
handle.close()
print(result)
b'<?xml version="1.0" encoding="UTF-8" ?>\n<!DOCTYPE eInfoResult PUBLIC "-//NLM//DTD einfo 20190110//EN" "https://eutils.ncbi.nlm.nih.gov/eutils/dtd/20190110/einfo.dtd">\n<eInfoResult>\n<DbList>\n\n\t<DbName>pubmed</DbName>\n\t<DbName>protein</DbName>\n\t<DbName>nuccore</DbName>\n\t<DbName>ipg</DbName>\n\t<DbName>nucleotide</DbName>\n\t<DbName>structure</DbName>\n\t<DbName>genome</DbName>\n\t<DbName>annotinfo</DbName>\n\t<DbName>assembly</DbName>\n\t<DbName>bioproject</DbName>\n\t<DbName>biosample</DbName>\n\t<DbName>blastdbinfo</DbName>\n\t<DbName>books</DbName>\n\t<DbName>cdd</DbName>\n\t<DbName>clinvar</DbName>\n\t<DbName>gap</DbName>\n\t<DbName>gapplus</DbName>\n\t<DbName>grasp</DbName>\n\t<DbName>dbvar</DbName>\n\t<DbName>gene</DbName>\n\t<DbName>gds</DbName>\n\t<DbName>geoprofiles</DbName>\n\t<DbName>homologene</DbName>\n\t<DbName>medgen</DbName>\n\t<DbName>mesh</DbName>\n\t<DbName>ncbisearch</DbName>\n\t<DbName>nlmcatalog</DbName>\n\t<DbName>omim</DbName>\n\t<DbName>orgtrack</DbName>\n\t<DbName>pmc</DbName>\n\t<DbName>popset</DbName>\n\t<DbName>proteinclusters</DbName>\n\t<DbName>pcassay</DbName>\n\t<DbName>protfam</DbName>\n\t<DbName>biosystems</DbName>\n\t<DbName>pccompound</DbName>\n\t<DbName>pcsubstance</DbName>\n\t<DbName>seqannot</DbName>\n\t<DbName>snp</DbName>\n\t<DbName>sra</DbName>\n\t<DbName>taxonomy</DbName>\n\t<DbName>biocollections</DbName>\n\t<DbName>gtr</DbName>\n</DbList>\n\n</eInfoResult>\n'
...... 아 개비스콘 아저씨 마렵네 이거.
그래서 짤줍했음.
아니 좀 깔끔하게 못 보여줘요?
from Bio import Entrez
Entrez.email = "blackholekun@gmail.com" # 내가 누구인지 말해주는 과정이 필요하다고...
# 이메일은 자기꺼 그냥 쓰세요
handle = Entrez.einfo()
record = Entrez.read(handle)
handle.close()
print(record)
from Bio import SeqIO
from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.efetch(db="nucleotide", id="NM_180778.4", rettype="gb", retmode="text")
record = SeqIO.read(handle, "genbank") # 레코드 형식으로 불러와서 부분적으로 취사선택할 수 있다.
handle.close() # 이거 꼭 닫아야됨?
print(record.description)
from Bio import Entrez
from Bio import SeqIO
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.efetch(db="nucleotide", id="NM_180778.4", rettype="gb", retmode="text")
record=SeqIO.read(handle,"genbank")
handle.close()
print(record)
SeqIO.write(record,"/home/koreanraichu/NM_180778.4.fasta","fasta")
# 사실 여기까진 쿡북에 없었는데 저장할 수 있는 방법이 없나 해서 해봤음.
사실 genbank 포맷인데 저걸 못열어서 fasta로 저장함...
Elink
뭐 관련된 걸 찾아준다고 한다. DNA 시퀀스로 유전자 찾아주듯이...
from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
pmid = "19304878"
record = Entrez.read(Entrez.elink(dbfrom="pubmed", id=pmid))
print(record[0]['IdList'])
['19304878']
이런 식으로 쓰고
for linksetdb in record[0]["LinkSetDb"]:
print(linksetdb["DbTo"], linksetdb["LinkName"], len(linksetdb["Link"]))
from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.egquery(term="kimchii")
record = Entrez.read(handle)
for row in record["eGQueryResult"]:
print(row["DbName"], row["Count"])
엥? 스펠? 마법사 찍고 올까요? 아니 그거 아냐... 물론 난 작년에 전직했지만 아무튼 그거 아님. 여기서 스펠은 스펠링(철자)의 스펠이지 마법 스펠이 아니다.
from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.espell(term="Conpanilactobacillus")
record = Entrez.read(handle)
print(record["Query"])
print(record["CorrectedQuery"])
Conpanilactobacillus
companilactobacillus
위: 오타 교정 전 아래: 오타 교정 후
사용예
개발자가 또 올려놔서 해봤음
Pubmed와 Medline의 대환장 콜라보레이션
일단 최종 코드가 두 개다.
from Bio import Entrez
from Bio import Medline
import numpy as np
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.egquery(term="kimchii")
record = Entrez.read(handle)
for row in record["eGQueryResult"]:
if row["DbName"] == "pubmed":
print(row["Count"])
# 쿼리를 통해 kimchii가 들어가는 걸 데려온 다음, pubmed에 있는 것만 뽑았다.
handle = Entrez.esearch(db="pubmed", term="kimchii", retmax=50)
record = Entrez.read(handle)
handle.close()
idlist = record["IdList"]
# 이제 리스트가 많으면 뽑아서 확인할 때 알아서 numpy를 불러봅니다.
handle = Entrez.efetch(db="pubmed", id=idlist, rettype="medline",retmode="text")
records = Medline.parse(handle)
# 이렇게 된 이상 medline으로 간다!
handle = Entrez.efetch(db="pubmed", id=idlist, rettype="medline",retmode="text")
records = Medline.parse(handle)
for record in records:
print("title:", record.get("TI", "?"))
print("authors:", record.get("AU", "?"))
print("source:", record.get("SO", "?"))
print("")
# 전체 출력은 여기서(전체 결과 중 제목, 저자, source를 출력)
from Bio import Entrez
from Bio import Medline
import numpy as np
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.egquery(term="kimchii")
record = Entrez.read(handle)
for row in record["eGQueryResult"]:
if row["DbName"] == "pubmed":
print(row["Count"])
# 쿼리를 통해 kimchii가 들어가는 걸 데려온 다음, pubmed에 있는 것만 뽑았다.
handle = Entrez.esearch(db="pubmed", term="kimchii", retmax=50)
record = Entrez.read(handle)
handle.close()
idlist = record["IdList"]
# 이제 리스트가 많으면 뽑아서 확인할 때 알아서 numpy를 불러봅니다.
handle = Entrez.efetch(db="pubmed", id=idlist, rettype="medline",retmode="text")
records = Medline.parse(handle)
# 이렇게 된 이상 medline으로 간다!
search_title = "Leuconostoc"
for record in records:
if not "TI" in record:
continue
if search_title in record["TI"]:
print("Keyword %s found: %s" % (search_title, record["TI"]))
# 제목에 Leuconostoc이 있는 것만 뽑아달라!
위 코드와 아래 코드의 공통 분기는 다음과 같다.
1. 쿼리로 kimchii 찾아서 2. pubmed에 있는 것만 뽑아서 3. ID 리스트업해서 4. 어디가 찾아야지
그리고 갈라지는 분기는 여기서 전체 목록을 출력하는가 or Leuconostoc으로 제목을 한번 더 필터링하는가 여부. 저거 둘 다 쓰면 안먹혀서 분기 갈라진거다...
Nucleotide 코드 받기
from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.egquery(term="Arabidopsis")
record = Entrez.read(handle)
for row in record["eGQueryResult"]:
if row["DbName"] == "nuccore":
print(row["Count"])
# 애기장대로 긁어와서 nuccore에 있는 것만 출력해라
handle = Entrez.esearch(db="nucleotide", term="Arabidopsis", retmax=800, idtype="acc")
record = Entrez.read(handle)
handle.close()
# 올해 얼마나 나왔는지는 모르겠고 일단 800개만 뽑아보자
print(record["IdList"][:5])
idlist = ",".join(record["IdList"][:5])
# 레코드에서 다섯개 뽑아서 리스트업한다.
handle = Entrez.efetch(db="nucleotide", id=idlist, retmode="xml")
records = Entrez.read(handle)
print(records[0].keys())
# 0번째 레코드 키좀 주세요
print(records[0]["GBSeq_length"])
# 시퀀스 길이 줘봐봐
얘도 쿼리 긁는 건 똑같다. 찾는 DB가 다르지... 그리고 예제 코드와 달리 데이터가 방대해서 어쩔 수 없이 위에서 5개만 뽑았다. 안그러면 안끝나 이거...
순순히 Genbank 레코드를 내놓는다면 유혈사태는 면할것이다
from Bio import Entrez
from Bio import SeqIO
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.egquery(term="Leuconostoc AND kimchii") # 뭐야 연산자 이렇게 쓰면 됨?
record = Entrez.read(handle)
for row in record["eGQueryResult"]:
if row["DbName"] == "nuccore":
print(row["Count"])
# nuccore에 수록되어있는 L.kimchii에 대해 찾아보자. (Leuconostoc kimchii)
handle = Entrez.esearch(db="nuccore", term="Leuconostoc AND kimchii")
record = Entrez.read(handle)
gi_list = record["IdList"]
# GI list를 만들어서 출력해보자
gi_str = ",".join(gi_list[0:5])
handle = Entrez.efetch(db="nuccore", id=gi_str, rettype="gb", retmode="text")
text = handle.read()
print(text)
# 이거 저장 안됨?
handle = Entrez.efetch(db="nuccore", id=gi_str, rettype="gb", retmode="text")
records = SeqIO.parse(handle, "gb")
for record in records:
print("%s, length %i, with %i features" % (record.name, len(record), len(record.features)))
# Parsing에는 for문이 국룰이다.
BLAST는 Basic local alignment search tool의 준말로, 보통 블래스트라고 한다. 지나가는 2~30대를 붙잡고 크흐~대 기억이~ 하면 지~난 사랑이~ 가 나오듯, 지나가던 생물학도를 붙잡고 블래스트? 하면 NCBI! 하는 정도로 국민 툴이고 전공수업을 듣다 보면 한 번은 쓰고 넘어가는 툴이기도 하다. 애초에 전공이 생물학인지를 먼저 물어봐야 하는 거 아니냐
이 글에 있는 게 어렵다… 그러면 조용히 구글을 열고 BLAST 검색해서 웹사이트 들어가자. 거기서 할 수 있다.
NCBIWWW
www가 우리말로 ㅋㅋㅋ이긴 한데 그건 일본에서 통하는 말이고... 보통은 월드 와이드 웹이다. BLAST를 따로 깔지 않고도 돌리려면 저게 있어야 한다. 물론 BLAST 깔려있는 서버가 있으면 저거 쓰지 말고 걍 로컬 돌리자. 저거 생각보다 느리다. 어느정도로 느리냐면 돌려놓고 강의실 있는 건물에 있는 카페에서 커피 한 잔 마시고 오면 끝난다. 돌려놓고 설거지 한번 하면 된다
from Bio.Blast import NCBIWWW
아무튼 오늘은 이 친구를 필두로 뭘 또 많이 데려올 예정이다.
인터넷으로 BLAST 돌리기-기본편
참고로 로컬로 돌리는 건 생략. 놋북에 그게 안 깔려있다.
from Bio.Blast import NCBIWWW
result_handle = NCBIWWW.qblast("blastn", "nt", "8332116")
print(result_handle)
기본 코드는 이거다. 이걸 실행하면 더럽게 오래 걸리면서...
<_io.StringIO object at 0x7f2f672ff280>
이러시는 이유가 있으실 것 아니예요.
쟤가 왜 저러는지는 이따 저장할 때 설명하기로 하고, 저 세 개가 뭔지에 대해 알아보자. Biopython에서 BLAST를 돌릴 때는 qblast()를 쓰는데, 거기서 필요한 3대 요소가 뭘로 돌릴래/데이터베이스 어디서 찾을래/GI 번호이다. GI 번호는
이런 식으로 Genbank에서 뭔가를 찾으면 나오는 번호. 즉 저 코드는 GI 번호가 8332116인 유전자의 시퀀스로 blastn을 돌리라는 얘기다. (DB 저거 약어같은데 나도 뭔지 모름) 혹시 여기에 대해 더 궁금한 점이 있다면
from Bio import SeqIO
from Bio.Blast import NCBIWWW
record = SeqIO.read("/home/koreanraichu/CaMV.fasta", format="fasta")
result_handle = NCBIWWW.qblast("blastn", "nt", record.seq)
print(result_handle)
from Bio import SeqIO
from Bio.Blast import NCBIWWW
record = SeqIO.read("m_cold.fasta", format="fasta")
result_handle = NCBIWWW.qblast("blastn", "nt", record.format("fasta"))
print(result_handle)
셋 다 쓸 수 있다. Parse는... 그거 쓸거면 for로 불러와서 시퀀스 픽한다음 몇번쨰 시퀀스로 찾을지 골라야 하지 않나...
결과 파싱하기
from Bio import SeqIO
from Bio.Blast import NCBIWWW
record = SeqIO.read("/home/koreanraichu/CaMV.fasta", format="fasta")
result_handle = NCBIWWW.qblast("blastn", "nt", record.seq)
result_handle = open("/home/koreanraichu/CaMV_BLAST.xml")
from Bio import SeqIO
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML
record = SeqIO.read("/home/koreanraichu/CaMV.fasta", format="fasta")
result_handle = NCBIWWW.qblast("blastn", "nt", record.seq)
blast_record = NCBIXML.read(result_handle)
print(blast_record)
둘 다 해봤는데 전 아래꺼요... (위에꺼 오류남...)
save_file = open("my_blast.xml", "w") # 이거 뭐 하는 구문인데 경초 못 찾는다고 에러뜨냐
save_file.write(result_handle.read())
save_file.close()
result_handle.close()
저장할거면 얘도 추가해야된다. 아까 처음 실행할 때 뭔 괴상한 문자만 땡 하고 끝났던 걸 기억하는가? Biopython에서 BLAST 돌리면 결과를 알아서 저장해주는 게 아니라 BLAST 결과에 대한 핸들을 생성한다. 그리고 그 핸들을 읽는 게 다라서 밑에 올린 저장용 코드를 따로 또 쳐 줘야 읽은 핸들을 xml파일로 저장한다.
from Bio import SeqIO
from Bio.Blast import NCBIWWW
record = SeqIO.read("/home/koreanraichu/CaMV.fasta", format="fasta")
result_handle = NCBIWWW.qblast("blastn", "nt", record.seq)
save_file = open("my_blast.xml", "w") # 이거 뭐 하는 구문인데 경로 못 찾는다고 에러뜨냐
save_file.write(result_handle.read())
save_file.close()
result_handle.close()
# 위에 네 줄까지 추가해야 결과가 저장된다.
from Bio import SeqIO
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML
for seq_record in SeqIO.parse("/home/koreanraichu/coronavirus.gb","genbank"): # 이거 이래 불러야되나
print(seq_record.seq) # Genbank 가져올거면 이렇게 쓰세요
result_handle = NCBIWWW.qblast("blastn", "nt", seq_record.seq)
blast_record = NCBIXML.parse(result_handle)
# 결과가 하나면 read, 여러개면 parse
# 이것도 parse에 for문이 국룰인가?
save_file = open("coronavirus.xml", "w") # 이거 뭐 하는 구문인데 경로 못 찾는다고 에러뜨냐
save_file.write(result_handle.read())
save_file.close()
result_handle.close()
# 위에 네 줄까지 추가해야 결과가 저장된다
근데 Genbank는 레코드 뽑는 것까지 성공했는데 BLAST 결과가 없음으로 나온다. 웹으로 한번 돌려봐야되는데 문제는 .gb파일이 안열린다는 거… (게임보이 뭐시기가 필요해서 파이참으로 바이오파이썬 불러서 열어야 한다)
아무튼 이런 식으로 저장된다. (경로 지정을 안 했더니 파이참 폴더행)
BLAST 레코드
아까 저장한 XML파일! 와! 우리꺼! 결과! 열어보자!
원래 저렇게 생긴 파일이다.
from Bio import SeqIO
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML
with open("my_blast.xml") as bf:
blast_records = NCBIXML.parse(bf)
E_VALUE_THRESH = 0.04 # 일단 이거 높으면 안좋은거다
for blast_record in blast_records:
for alignment in blast_record.alignments:
for hsp in alignment.hsps:
if hsp.expect < E_VALUE_THRESH:
print("****Alignment****")
print("sequence:", alignment.title)
print("length:", alignment.length)
print("e value:", hsp.expect)
print(hsp.query[0:75] + "...")
print(hsp.match[0:75] + "...")
print(hsp.sbjct[0:75] + "...")
# hsp 영역을 추출해 내는 스크립트.
여기서 특정 영역을 추출할 수 있는데, 이 코드에서는 hsp 영역을 추출했다. 웹으로 맨날 보는 그 부분... (결과는 코드블럭에서 길다고 짤렸다...)
응? 벡터요? 시범조교로 벡터 쓰심? 아니 벡터가 아니라 CaMV 시퀀스 썼는데 이게 35S promoter 시퀀스였던 모양이다. CaMV는 콜리플라워 모자이크 바이러스로, 볼티모어 클래스 7(DNA-RT)에 속한다. 엔빌롭이 있고… 근데 저 바이러스가 벡터랑 뭔 상관이냐고? 유전자 과발현 벡터에는 CaMV 바이러스의 35S promoter가 들어간다.
from Bio import SeqIO
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML
record = SeqIO.read("/home/koreanraichu/coronavirus.gb", format="genbank")
result_handle = NCBIWWW.qblast("blastn", "nt", record.seq)
blast_records = NCBIXML.parse(result_handle)
E_VALUE_THRESH = 0.04 # 일단 이거 높으면 안좋은거다
for blast_record in blast_records:
for alignment in blast_record.alignments:
for hsp in alignment.hsps:
if hsp.expect < E_VALUE_THRESH:
print("****Alignment****")
print("sequence:", alignment.title)
print("length:", alignment.length)
print("e value:", hsp.expect)
print(hsp.query[0:75] + "...")
print(hsp.match[0:75] + "...")
print(hsp.sbjct[0:75] + "...")
# hsp 영역을 추출해 내는 스크립트.
아 저장 필요없고 결과나 보여주셈! 하면 다이렉트로 돌려서 결과 보는 것도 된다. 참고로 저기 있는 E-value는 일단 클수록 안 좋다.
SearchIO
뭔 IO가 들어가요? 저거 사실 인풋/아웃풋 약어다.
아까 그 으아악 내눈을 유발했던 XML파일... 거기서 우리가 원하는 정보를 볼 수 있게 해 주는 게 이 친구 역할이다. (대충 개비스콘 아저씨 짤)
from Bio import SearchIO
blast_qresult = SearchIO.read("my_blast.xml", "blast-xml")
print(blast_qresult)
from Bio import SearchIO
blast_qresult = SearchIO.read("my_blast.xml", "blast-xml")
blast_hit = blast_qresult[0]
print(blast_hit)
Query: No
definition line
Hit: gi|1835923315|gb|MN508834.1| (14584)
Binary vector pRATIO3212-SMXL7, complete sequence
HSPs: ---- -------- --------- ------ --------------- ---------------------
# E-value Bit score Span Query range Hit range
---- -------- --------- ------ --------------- ---------------------
0 3.7e-174 623.45 345 [0:345] [791:1136]
색인할 수 있다고 한다. 참고로 저 [0]을 잘 기억해두자.
HSP
힛샥 프로테인...아니고 High scoring pair의 준말이다.
from Bio import SearchIO
blast_qresult = SearchIO.read("my_blast.xml", "blast-xml")
blast_hsp = blast_qresult[0][0]
print(blast_hsp)
Query: No definition line
Hit: gi|1835923315|gb|MN508834.1| Binary vector pRATIO3212-SMXL7, com...
Query range: [0:345] (1)
Hit range: [791:1136] (-1)
Quick stats: evalue 3.7e-174; bitscore 623.45
Fragments: 1 (345 columns)
Query - TCTCTCCAAATGAAATGAACTTCCTTATATAGAGGAAGGGTCTTGCGAAGGATAGTGGG~~~TCAAT
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||~~~|||||
Hit - TCTCTCCAAATGAAATGAACTTCCTTATATAGAGGAAGGGTCTTGCGAAGGATAGTGGG~~~TCAAT
Process finished with exit code 0
얘는 그래서 힛이랑 달리 색인할 때 리스트 두 개가 들어간다. 각각 맨 첫번째 힛에서 맨 첫번째 HSP를 뜻한다. 색인할 때 하나 빼먹으면 힛이 출력된다.
# 얘네는 hsp 위치까지 지정해줘야 출력한다
print(blast_hsp.query_range) # 첫번째 hit의 첫번째 hsp의 query range
print(blast_hsp.evalue) # 첫번째 hit의 첫번째 hsp의 evalue
print(blast_hsp.hit_start) # 첫번째 hit의 첫번째 hsp의 hit start coordinates
print(blast_hsp.query_span) # 첫번째 hit의 첫번째 hsp의 how many residues in the query sequence(residue 몇개냐)
print(blast_hsp.aln_span) # 첫번째 hit의 첫번째 hsp의 alignment 길이
print(blast_hsp.gap_num) # 첫번째 hit의 첫번째 hsp의 gap no.
print(blast_hsp.ident_num) # 첫번째 hit의 첫번째 hsp의 identical residue no.
print(blast_hsp.aln) # 배열 뽑아줘!
(0, 345)
3.74077e-174
791
345
345
0
345
# 밑에는 배열이다
Alignment with 2 rows and 345 columns
TCTCTCCAAATGAAATGAACTTCCTTATATAGAGGAAGGGTCTT...AAT No
TCTCTCCAAATGAAATGAACTTCCTTATATAGAGGAAGGGTCTT...AAT gi|1835923315|gb|MN508834.1|
참고로 힛 상태에서 저거 쓰면 오류난다.
HSPFragment
가장 하위 개념. 얘는 색인할 때 []가 세 개 들어간다.
blast_frag = blast_qresult[0][0][0] # 첫번째 hit의 첫번째 hsp의 첫번째 fragment
print(blast_frag)
Query: No definition line
Hit: gi|1835923315|gb|MN508834.1| Binary vector pRATIO3212-SMXL7, com...
Query range: [0:345] (1)
Hit range: [791:1136] (-1)
Fragments: 1 (345 columns)
Query - TCTCTCCAAATGAAATGAACTTCCTTATATAGAGGAAGGGTCTTGCGAAGGATAGTGGG~~~TCAAT
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||~~~|||||
Hit - TCTCTCCAAATGAAATGAACTTCCTTATATAGAGGAAGGGTCTTGCGAAGGATAGTGGG~~~TCAAT
Process finished with exit code 0
...HSP랑 차이 없는데?
blast_frag = blast_qresult[0][0][0] # 첫번째 hit의 첫번째 hsp의 첫번째 fragment
print(blast_frag.hit) # hit sequence(sequence object)
print(blast_frag.query_start) # 시작 좌표
print(blast_frag.hit_strand) # hi sequence strand 수
ID: gi|1835923315|gb|MN508834.1|
Name: aligned hit sequence
Description: Binary vector pRATIO3212-SMXL7, complete sequence
Number of features: 0
/molecule_type=DNA
Seq('TCTCTCCAAATGAAATGAACTTCCTTATATAGAGGAAGGGTCTTGCGAAGGATA...AAT')
0
-1
코테 본다는 얘기는 못 들은 것 같고... 그냥 면접보다가 어디 실력 좀 볼까? 해서 나온거긴 한데... IDE도 셋업 안 된 상태에서 jupyter 데모로 봤었음... 근데 Python을 안 쓰면 나한테 연락을 안 할것인데 이건 뭔 상황인지 모르겠고...
거기다가 실습할 때 썼던 파일들 올릴랬는데 깃헙 뻑나서 오전중에는 그거랑 씨름함... 해결은 봤죠. 참고로 어제 사용한 FASTAQ파일과 오늘 사용한 FASTAQ파일은 다릅니다. 이거 구하기도 개빡셈.
모듈 모셔오기
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
import numpy as np
import pandas as pd
세번째줄 안불러도 됨 그거 쩌리예요 판다스 넘파이는 4, 5번 풀려고 부른거
1. FastaQ 파일을 Fasta 파일로 변환해라
# FASTAQ file 불러오기
handle="/home/koreanraichu/sra_data_mo.fastq"
for record in SeqIO.parse(handle,"fastq"):
print(record.id)
# FASTA파일로 쓰기
convert_fasta=SeqIO.convert(handle, "fastq", "/home/koreanraichu/sra_data_mo.fasta","fasta")
어제는 내 컴퓨터도 아니고 Jupyter 데모여서 파일 자체 확인을 못 했는데, 이런 상황에서 잘 됐나 확인할거면 SeqIO read로 불러와보자. SeqIO는 뭉텅이로 있으면 parse만 된다. (정확히는 parse+for)
(변환 후)
2. Fasta파일에서 A, T, G, C 세기
handle2="/home/koreanraichu/sra_data_mo.fasta"
for record2 in SeqIO.parse(handle2,"fasta"):
print(Seq(record2.seq).count("A"))
print(type(Seq(record2.seq).count("A")))
# Count로 세는 건 되는데 합계를 안 내준다. int형인데 더해주질 않음.
정수형인데 왜 더하질 못하고 각개출력하니... 내가 저 상태에서 리스트로 만들려고 append 했더니 아 실패요. 그것도 각개 리스트로 나오고 A=A.append 줬더니 None됨...
니네 왜그래...
+그거 해결봤음. 진짜 웃긴 과정이었음...
handle2="/home/koreanraichu/sra_data_mo.fasta"
A = []
for record2 in SeqIO.parse(handle2,"fasta"):
A.append(Seq(record2.seq).count("A"))
print(A)
애초에 만들었던 빈 리스트를 밖으로 빼니까 되는데요? 와 신박하네.
handle2="/home/koreanraichu/sra_data_mo.fasta"
A = [] #아니 이게 이렇게 해결된다고?
T = []
G = []
C = []
for record2 in SeqIO.parse(handle2,"fasta"):
A.append(Seq(record2.seq).count("A"))
T.append(Seq(record2.seq).count("T"))
G.append(Seq(record2.seq).count("G"))
C.append(Seq(record2.seq).count("C"))
print("A: ",sum(A),"T: ",sum(T),"G: ",sum(G),"C: ",sum(C))
코드_최종.py
아니 진짜 저걸로 해결되는거면 판다스 표 개판 오분전 된것도 저걸로 걍 정리해버리자.
3. 전체 Sequence 세기
print(sum(len(r) for r in SeqIO.parse(handle2, "fasta")))
와 이걸 이렇게 해결보네
4. Top 10 length
for record2 in SeqIO.parse(handle2,"fasta"):
record_id=np.array(record2.id)
record_len=np.array(len(record2.seq))
record_table=pd.DataFrame({"ID":record_id, "length":record_len},index=[0])
print(record_table)
# Dataframe이 이상하게 생성된다.
"""
ID length
0 SRR000021.37.2 249
이런 식으로 항목 하나당 한 데이터프레임이 생성되는데 Array까지는 정상적으로 생성됨.
"""
판다스로 표 만들어서 꼽아보려고 했더니 표가 개판 오분전이 되었다... 이쯤되면 array 차원부터 한번 봐야겠는데? array 차원부터 망했는데 이거?
# A4. Top 10 length
record_id=[] # 야 이걸 이렇게 해결보네
record_len=[]
for record2 in SeqIO.parse(handle2,"fasta"):
record_id.append(record2.id)
record_len.append(len(record2.seq))
record_id_array=np.array(record_id)
record_len_array=np.array(record_len)
비효율적으로 보일 지 몰라도 이정도면 본인 수준에서는 장족의 발전임다... 어레이 제대로 된 것도 확인했음.
