Biopython으로 Phylogenetic tree 그리기
걍 메가 쓰세여... 메가가 편해... 짱이야 메가...
계통수?
계통수... 그러니까 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
'Coding > Python' 카테고리의 다른 글
| Biopython으로 Clusting analysis 하기 (이론편) (0) | 2022.08.21 |
|---|---|
| Biopython으로 Sequence motif analysis 하기 (0) | 2022.08.21 |
| Biopython으로 PDB 탐방하기 (0) | 2022.08.21 |
| Biopython으로 Swiss-prot과 ExPASy 데이터베이스 탐방하기 (0) | 2022.08.21 |
| Biopython으로 Entrez database 탐방하기 (0) | 2022.08.21 |
