- Biopython-Clustering 입력 인자 2022.08.21
- Biopython으로 Clusting analysis 하기 (실전편) 2022.08.21
- Biopython으로 Clusting analysis 하기 (이론편) 2022.08.21
- Biopython으로 Sequence motif analysis 하기 2022.08.21
- Biopython으로 Phylogenetic tree 그리기 2022.08.21
- Biopython으로 PDB 탐방하기 2022.08.21
- Biopython으로 Swiss-prot과 ExPASy 데이터베이스 탐방하기 2022.08.21
- Biopython으로 Entrez database 탐방하기 2022.08.21
- Biopython으로 BLAST 돌려보기 2022.08.20
- 번외편-코딩테스트 풀이 2022.08.20
이거 튜토리얼에서 그냥 이런게 있다 하고 넘어가서 적음...
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)
'Coding > Python' 카테고리의 다른 글
| 텍스트 입력받아서 Wordcloud 만들기 (0) | 2022.08.21 |
|---|---|
| Biopython-Entrez에서 논문 제목 긁어와서 Wordcloud 만들기 (0) | 2022.08.21 |
| Biopython으로 Clusting analysis 하기 (실전편) (0) | 2022.08.21 |
| Biopython으로 Clusting analysis 하기 (이론편) (0) | 2022.08.21 |
| Biopython으로 Sequence motif analysis 하기 (0) | 2022.08.21 |
와! 드디어 실전인가요? 근데 실전 생각보다 노잼임...
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): 81Array로 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)
# 저장아무튼 이런 식으로 쓴다고...
'Coding > Python' 카테고리의 다른 글
| Biopython-Entrez에서 논문 제목 긁어와서 Wordcloud 만들기 (0) | 2022.08.21 |
|---|---|
| Biopython-Clustering 입력 인자 (0) | 2022.08.21 |
| Biopython으로 Clusting analysis 하기 (이론편) (0) | 2022.08.21 |
| Biopython으로 Sequence motif analysis 하기 (0) | 2022.08.21 |
| Biopython으로 Phylogenetic tree 그리기 (0) | 2022.08.21 |
분량도 분량인데 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이다.
'Coding > Python' 카테고리의 다른 글
| Biopython-Clustering 입력 인자 (0) | 2022.08.21 |
|---|---|
| Biopython으로 Clusting analysis 하기 (실전편) (0) | 2022.08.21 |
| Biopython으로 Sequence motif analysis 하기 (0) | 2022.08.21 |
| Biopython으로 Phylogenetic tree 그리기 (0) | 2022.08.21 |
| Biopython으로 PDB 탐방하기 (0) | 2022.08.21 |
모티브 찾으면 구조 모티브가 나오는데 이거는 '단백질이나 핵산과 같은 사슬 모양의 생체 분자에서 진화적으로 관련이 없는 다양한 분자에서 나타나는 일반적인 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
CGCAGACGAGreverse 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.00jaspar 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.00ID 지우면 이렇게 부를 수는 있는데 그럴거면 불러서 뭐함?
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.00Jaspar->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.02normalize가 들어간 걸 보면 정규화를 하는 모양. 정규화나 표준화는 원래 인공지능용 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.21reverse complement로도 된다.
print(pwm.consensus)
print(pwm.anticonsensus)
print(pwm.degenerate_consensus)
# 제일 많은거/제일 적은거/Degenerate consensusCACGTG
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 CTCGTCTGCGtest_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.005threshold요? 그거 밑에 나옴. 아마 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이 된다.
'Coding > Python' 카테고리의 다른 글
| Biopython으로 Clusting analysis 하기 (실전편) (0) | 2022.08.21 |
|---|---|
| Biopython으로 Clusting analysis 하기 (이론편) (0) | 2022.08.21 |
| Biopython으로 Phylogenetic tree 그리기 (0) | 2022.08.21 |
| Biopython으로 PDB 탐방하기 (0) | 2022.08.21 |
| Biopython으로 Swiss-prot과 ExPASy 데이터베이스 탐방하기 (0) | 2022.08.21 |
걍 메가 쓰세여... 메가가 편해... 짱이야 메가...
계통수?
계통수... 그러니까 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 |
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
FalseResidue에서 뽑아낼 수 있는 속성은 이름/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 factorCA
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.5872530070961592atom1=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="경로" 하면 된다.
'Coding > Python' 카테고리의 다른 글
| Biopython으로 Sequence motif analysis 하기 (0) | 2022.08.21 |
|---|---|
| Biopython으로 Phylogenetic tree 그리기 (0) | 2022.08.21 |
| Biopython으로 Swiss-prot과 ExPASy 데이터베이스 탐방하기 (0) | 2022.08.21 |
| Biopython으로 Entrez database 탐방하기 (0) | 2022.08.21 |
| Biopython으로 BLAST 돌려보기 (0) | 2022.08.20 |
근데 스위스프롯은 긁어와서 저장 안됨?
첫빠따는 파싱이 국룰이지
파싱 방법이 네 가지가 있는데 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'}뭐 이렇다는데... 이거 예제꺼 긁어왔는데 왜 결과가 다르냐...
'Coding > Python' 카테고리의 다른 글
| Biopython으로 Phylogenetic tree 그리기 (0) | 2022.08.21 |
|---|---|
| Biopython으로 PDB 탐방하기 (0) | 2022.08.21 |
| Biopython으로 Entrez database 탐방하기 (0) | 2022.08.21 |
| Biopython으로 BLAST 돌려보기 (0) | 2022.08.20 |
| 번외편-코딩테스트 풀이 (0) | 2022.08.20 |
이것이 그... BLAST 만든 NCBI에 있는 데이터베이스다. 미국답게 스케일 개크다.
들어가기 전에
보통 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){'DbList': ['pubmed', 'protein', 'nuccore', 'ipg', 'nucleotide', 'structure', 'genome', 'annotinfo', 'assembly', 'bioproject', 'biosample', 'blastdbinfo', 'books', 'cdd', 'clinvar', 'gap', 'gapplus', 'grasp', 'dbvar', 'gene', 'gds', 'geoprofiles', 'homologene', 'medgen', 'mesh', 'ncbisearch', 'nlmcatalog', 'omim', 'orgtrack', 'pmc', 'popset', 'proteinclusters', 'pcassay', 'protfam', 'biosystems', 'pccompound', 'pcsubstance', 'seqannot', 'snp', 'sra', 'taxonomy', 'biocollections', 'gtr']}
이러시는 이유가 있으실 것 아니예요... 하씨 넘파이 마렵네 이거...
Esearch
이제 본격적으로 찾아볼 시간이다. 핫핫 데이터내놔!
PubMed에서 논문 찾기
from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.esearch(db="pubmed", term="biopython[title]", retmax="40" )
record = Entrez.read(handle)
print(record['IdList'])['34434786', '22909249', '19304878']바이오파이썬으로 찾았다. (제목에 바이오파이썬이 있는 것) 근데 이 예제 노잼임... 시범조교 불러올거임...
은 바로 이녀석. C. kimchii다. (구 L. kimchii)
김치 related에 당당히 이름을 올렸다. 아니 그야 김치에서 나왔으니까… (자매품으로는 Leuconostoc kimchii도 있다) 근데 얘들아 두부김치도 related 아니냐...
from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.esearch(db="pubmed", term="kimchii[title]", retmax="40" )
record = Entrez.read(handle)
print(record['IdList'])['34320438', '31651376', '31270190', '30361979', '30172442', '29214493', '28920843', '27572507', '27002961', '26370793', '25425317', '25332883', '22140166', '21914872', '21221947', '20494991', '16055277', '15023961', '11931163', '11034505', '11034488'](대충 개비스콘 아저씨 위장 부여잡는 짤)
from Bio import Entrez
import numpy as np
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.esearch(db="pubmed", term="kimchii[title]", retmax="40" )
record = Entrez.read(handle)
result=np.array(record['IdList'])
result=result.reshape(7,3)
print(result)[['34320438' '31651376' '31270190']
['30361979' '30172442' '29214493']
['28920843' '27572507' '27002961']
['26370793' '25425317' '25332883']
['22140166' '21914872' '21221947']
['20494991' '16055277' '15023961']
['11931163' '11034505' '11034488']]
아 편안하네.
from Bio import Entrez
import numpy as np
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.esearch(db="pubmed", term="kimchii[title] and Leuconostoc[title]", retmax="40" )
record = Entrez.read(handle)
# 출력 안이뻐서 배열 만들었음
result=np.array(record['IdList'])
result=result.reshape(2,4)
print(result)
# 바람직한 출력이야![['34320438' '31270190' '30172442' '25425317']
['25332883' '21914872' '20494991' '11034505']]연산자도 먹힌다. 위 코드는 Leuconostoc kimchii를 찾는 연산자. 미역김치 있나요 그건 저거 다 찾아보면 나올듯
Nucleotide DB 찾기
from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.esearch(db="nucleotide", term="Arabidopsis[Orgn] AND LHT1[Gene]", idtype="acc")
record = Entrez.read(handle)
print("%s founded" % record["Count"])
print(record["IdList"])5 founded
['NM_001344354.1', 'NM_180778.4', 'NM_001344353.1', 'NC_003076.8', 'CP002688.1']애기장대(Arabodopsis)의 LHT1에 대해 찾으면 다섯 개가 나온다. 아까 걔는 없음.. ㅋㅋㅋㅋㅋㅋㅋㅋ
저널 찾기
from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.esearch(db="nlmcatalog", term="computational[Journal]", retmax="20")
record = Entrez.read(handle)
print("{} computational journals found".format(record["Count"]))
print("The first 20 are\n{}".format(record["IdList"]))187 computational journals found
The first 20 are
['101775476', '101775136', '468912', '468839', '467370', '101775780', '101774751', '466063', '101768752', '101765300', '101759185', '101752828', '464655', '101768811', '101755127', '101753951', '101753371', '101740904', '101737789', '101736625']저널? 그럼 이거 인공지능 저널 이런것도 찾아주나?
Epost
있으니까 보고 가긴 하는데 이거 뭔지는 모르겠음...
from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
id_list = ["19304878", "18606172", "16403221", "16377612", "14871861", "14630660"]
print(Entrez.epost("pubmed", id=",".join(id_list)).read())b'<?xml version="1.0" encoding="UTF-8" ?>\n<!DOCTYPE ePostResult PUBLIC "-//NLM//DTD epost 20090526//EN" "https://eutils.ncbi.nlm.nih.gov/eutils/dtd/20090526/epost.dtd"><ePostResult>\n\t<QueryKey>1</QueryKey>\n\t<WebEnv>MCID_618c9fbf207466500268a6a2</WebEnv>\n</ePostResult>\n'(대충 개비스콘 아저씨 속쓰린 짤)
Esummary
검색결과의 첫 번째 데이터에 대해 요약본을 보여준다.
from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.esummary(db="pubmed", id="31651376")
record = Entrez.read(handle)
info = record[0]
print("Journal info\nid: {}\nTitle: {}".format(record[0]["Id"], info["Title"]))Journal info
id: 31651376
Title: <i>Lactococcus kimchii</i> sp. nov., a new lactic acid bacterium isolated from kimchi.
Efetch
이거 봤으면 엔트레즈 다 본거다. ...그렇다고 보다말고 가지는 말고...
from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일
handle = Entrez.efetch(db="nucleotide", id="EU490707", rettype="gb", retmode="text")
print(handle.read())LOCUS EU490707 1302 bp DNA linear PLN 26-JUL-2016
DEFINITION Selenipedium aequinoctiale maturase K (matK) gene, partial cds;
chloroplast.
ACCESSION EU490707
VERSION EU490707.1
KEYWORDS .
SOURCE chloroplast Selenipedium aequinoctiale
ORGANISM Selenipedium aequinoctiale
Eukaryota; Viridiplantae; Streptophyta; Embryophyta; Tracheophyta;
Spermatophyta; Magnoliopsida; Liliopsida; Asparagales; Orchidaceae;
Cypripedioideae; Selenipedium.
REFERENCE 1 (bases 1 to 1302)
AUTHORS Neubig,K.M., Whitten,W.M., Carlsward,B.S., Blanco,M.A., Endara,L.,
Williams,N.H. and Moore,M.
TITLE Phylogenetic utility of ycf1 in orchids: a plastid gene more
variable than matK
JOURNAL Plant Syst. Evol. 277 (1-2), 75-84 (2009)
REFERENCE 2 (bases 1 to 1302)
AUTHORS Neubig,K.M., Whitten,W.M., Carlsward,B.S., Blanco,M.A.,
Endara,C.L., Williams,N.H. and Moore,M.J.
TITLE Direct Submission
JOURNAL Submitted (14-FEB-2008) Department of Botany, University of
Florida, 220 Bartram Hall, Gainesville, FL 32611-8526, USA
FEATURES Location/Qualifiers
source 1..1302
/organism="Selenipedium aequinoctiale"
/organelle="plastid:chloroplast"
/mol_type="genomic DNA"
/specimen_voucher="FLAS:Blanco 2475"
/db_xref="taxon:256374"
gene <1..>1302
/gene="matK"
CDS <1..>1302
/gene="matK"
/codon_start=1
/transl_table=11
/product="maturase K"
/protein_id="ACC99456.1"
/translation="IFYEPVEIFGYDNKSSLVLVKRLITRMYQQNFLISSVNDSNQKG
FWGHKHFFSSHFSSQMVSEGFGVILEIPFSSQLVSSLEEKKIPKYQNLRSIHSIFPFL
EDKFLHLNYVSDLLIPHPIHLEILVQILQCRIKDVPSLHLLRLLFHEYHNLNSLITSK
KFIYAFSKRKKRFLWLLYNSYVYECEYLFQFLRKQSSYLRSTSSGVFLERTHLYVKIE
HLLVVCCNSFQRILCFLKDPFMHYVRYQGKAILASKGTLILMKKWKFHLVNFWQSYFH
FWSQPYRIHIKQLSNYSFSFLGYFSSVLENHLVVRNQMLENSFIINLLTKKFDTIAPV
ISLIGSLSKAQFCTVLGHPISKPIWTDFSDSDILDRFCRICRNLCRYHSGSSKKQVLY
RIKYILRLSCARTLARKHKSTVRTFMRRLGSGLLEEFFMEEE"
ORIGIN
1 attttttacg aacctgtgga aatttttggt tatgacaata aatctagttt agtacttgtg
61 aaacgtttaa ttactcgaat gtatcaacag aattttttga tttcttcggt taatgattct
121 aaccaaaaag gattttgggg gcacaagcat tttttttctt ctcatttttc ttctcaaatg
181 gtatcagaag gttttggagt cattctggaa attccattct cgtcgcaatt agtatcttct
241 cttgaagaaa aaaaaatacc aaaatatcag aatttacgat ctattcattc aatatttccc
301 tttttagaag acaaattttt acatttgaat tatgtgtcag atctactaat accccatccc
361 atccatctgg aaatcttggt tcaaatcctt caatgccgga tcaaggatgt tccttctttg
421 catttattgc gattgctttt ccacgaatat cataatttga atagtctcat tacttcaaag
481 aaattcattt acgccttttc aaaaagaaag aaaagattcc tttggttact atataattct
541 tatgtatatg aatgcgaata tctattccag tttcttcgta aacagtcttc ttatttacga
601 tcaacatctt ctggagtctt tcttgagcga acacatttat atgtaaaaat agaacatctt
661 ctagtagtgt gttgtaattc ttttcagagg atcctatgct ttctcaagga tcctttcatg
721 cattatgttc gatatcaagg aaaagcaatt ctggcttcaa agggaactct tattctgatg
781 aagaaatgga aatttcatct tgtgaatttt tggcaatctt attttcactt ttggtctcaa
841 ccgtatagga ttcatataaa gcaattatcc aactattcct tctcttttct ggggtatttt
901 tcaagtgtac tagaaaatca tttggtagta agaaatcaaa tgctagagaa ttcatttata
961 ataaatcttc tgactaagaa attcgatacc atagccccag ttatttctct tattggatca
1021 ttgtcgaaag ctcaattttg tactgtattg ggtcatccta ttagtaaacc gatctggacc
1081 gatttctcgg attctgatat tcttgatcga ttttgccgga tatgtagaaa tctttgtcgt
1141 tatcacagcg gatcctcaaa aaaacaggtt ttgtatcgta taaaatatat acttcgactt
1201 tcgtgtgcta gaactttggc acggaaacat aaaagtacag tacgcacttt tatgcgaaga
1261 ttaggttcgg gattattaga agaattcttt atggaagaag aa
//어... 뭔진 모르겠지만 일단 제가 잘못했어요?
긁어와서 레코드 만들기
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)Arabidopsis thaliana lysine histidine transporter 1 (LHT1), mRNA이런 식으로 취사선택해서 보는 것도 된다.
아 긁어왔으면 저장해야지
참고로 이거 저장은 Entrez 모듈 단독으로 못하고 SeqIO 불러야된다.
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"]))pubmed nuccore_pubmed 1
pubmed nuccore_pubmed_refseq 9
pubmed nuccore_pubmed_weighted 9for문 먹여도 된단다.
EGQuery
뭐야 쿼리가 왜 여기서 나와요? 그러게요.
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"])pubmed 50
pmc 154
mesh 5
books 0
pubmedhealth Error
omim 0
ncbisearch 4
nuccore 928
nucgss 0
nucest 0
protein 33512
genome 2
structure 0
taxonomy 0
snp 0
dbvar 0
gene 10109
sra 3
biosystems 295
unigene 0
cdd 0
clone 0
popset 25
geoprofiles 0
gds 0
homologene 0
pccompound 0
pcsubstance 0
pcassay 0
nlmcatalog 0
probe 0
gap 0
proteinclusters 0
bioproject 11
biosample 15
biocollections 0그냥 얘네 DB에 얼마나 있는지가 나온다. 워우.
ESpell
엥? 스펠? 마법사 찍고 올까요? 아니 그거 아냐... 물론 난 작년에 전직했지만 아무튼 그거 아님. 여기서 스펠은 스펠링(철자)의 스펠이지 마법 스펠이 아니다.
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문이 국룰이다.제목 저거 괜찮냐 이쯤되면 parsing에는 for문이 국룰 아니냐... 아무튼. ...아니 잠깐만 쿼리 연산자 저렇게 쓰는거임?
아무튼 그렇다.
Taxonomy lineage
from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일!!!
handle = Entrez.esearch(db="Taxonomy", term="Lactobacillus")
record = Entrez.read(handle)
print(record['IdList'])
# Taxonomy에서 lactobacillus를 찾아 ID를 내놓아라!
handle = Entrez.efetch(db="Taxonomy", id="1578", retmode="xml")
records = Entrez.read(handle)
print(records[0]['Lineage'])
# Lineage 주세요!from Bio import Entrez
Entrez.email = "blackholekun@gmail.com"
# 아 맞다 메일!!!
handle = Entrez.esearch(db="Taxonomy", term="Enterococcus")
record = Entrez.read(handle)
print(record['IdList'])
# Taxonomy에서 Enterococcus를 찾아 ID를 내놓아라!
handle = Entrez.efetch(db="Taxonomy", id=record['IdList'], retmode="xml")
records = Entrez.read(handle)
print(records[0]['Lineage'])
# Lineage 주세요!위 코드와 아래 코드 둘 다 일단 똑같은데, 아래 코드는 taxonomy 번호 결과를 아예 변수로 떄려박아서 다이렉트로 찾아준다. (위 코드는 직접 입력해야 함)
'Coding > Python' 카테고리의 다른 글
| Biopython으로 PDB 탐방하기 (0) | 2022.08.21 |
|---|---|
| Biopython으로 Swiss-prot과 ExPASy 데이터베이스 탐방하기 (0) | 2022.08.21 |
| Biopython으로 BLAST 돌려보기 (0) | 2022.08.20 |
| 번외편-코딩테스트 풀이 (0) | 2022.08.20 |
| Biopython으로 MSA 해보기 (0) | 2022.08.20 |
쿡북에 있는 것 중 일부 생략했다. (BLAT이랑 8.3~8.5부분)
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 저거 약어같은데 나도 뭔지 모름) 혹시 여기에 대해 더 궁금한 점이 있다면
한번 읽어보자. NCBI에서 제공하는 가이드이다. (물론 영어다)
인터넷으로 BLAST 돌리기-FASTA 파일로 돌리기
Genbank도 돌릴 수는 있는데 뭐가 문젠지 결과를 조졌다... 결과 파일에 암것도 없음.
from Bio.Blast import NCBIWWW
fasta_string = open("/home/koreanraichu/CaMV.fasta").read()
result_handle = NCBIWWW.qblast("blastn", "nt", fasta_string)
print(result_handle)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 영역을 추출했다. 웹으로 맨날 보는 그 부분... (결과는 코드블럭에서 길다고 짤렸다...)
****Alignment****
sequence: gi|1835923315|gb|MN508834.1| Binary vector pRATIO3212-SMXL7, complete sequence
length: 14584
e value: 3.74077e-174
TCTCTCCAAATGAAATGAACTTCCTTATATAGAGGAAGGGTCTTGCGAAGGATAGTGGGATTGTGCGTCATCCCT...
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||...
TCTCTCCAAATGAAATGAACTTCCTTATATAGAGGAAGGGTCTTGCGAAGGATAGTGGGATTGTGCGTCATCCCT...이런 게 계~속 반복된다.
응? 벡터요? 시범조교로 벡터 쓰심? 아니 벡터가 아니라 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)Program: blastn (2.12.0+)
Query: No (345)
definition line
Target: nt
Hits: ---- ----- ----------------------------------------------------------
# # HSP ID + description
---- ----- ----------------------------------------------------------
0 1 gi|1835923315|gb|MN508834.1| Binary vector pRATIO3212-...
1 1 gi|1835923313|gb|MN508833.1| Binary vector pRATIO3212-...
2 1 gi|1835923311|gb|MN508832.1| Binary vector pRATIO1212-...
3 1 gi|1755168873|gb|MN055609.1| Plant binary expression v...
4 1 gi|1755168871|gb|MN055608.1| Plant binary expression v...
5 1 gi|1755168869|gb|MN055607.1| Plant binary expression v...
6 1 gi|1755168867|gb|MN055606.1| Plant binary expression v...
7 1 gi|1755168865|gb|MN055605.1| Plant binary expression v...
8 1 gi|1755168863|gb|MN055604.1| Plant binary expression v...
9 1 gi|1732472841|gb|MK896906.1| Plant binary expression v...
10 1 gi|1732472839|gb|MK896905.1| Plant binary expression v...
11 1 gi|1732472836|gb|MK896904.1| Plant binary expression v...
12 1 gi|1732472834|gb|MK896903.1| Plant binary expression v...
13 3 gi|1732472826|gb|MK896900.1| Plant binary expression v...
14 3 gi|1732472823|gb|MK896899.1| Plant binary expression v...
15 3 gi|1732472819|gb|MK896898.1| Plant binary expression v...
16 3 gi|1732472816|gb|MK896897.1| Plant binary expression v...
17 3 gi|1732472812|gb|MK896896.1| Plant binary expression v...
18 1 gi|1693859709|gb|MK204379.1| Binary vector pBI-Xeg, co...
19 1 gi|1624318504|gb|MK501803.1| Expression vector pHytru,...
20 1 gi|1050047859|dbj|AB830573.1| Gateway vector pB5nYGW D...
21 1 gi|1050047854|dbj|AB830572.1| Gateway vector pB5nRGW D...
22 1 gi|1050047849|dbj|AB830571.1| Gateway vector pB5nGGW D...
23 1 gi|1050047842|dbj|AB830570.1| Gateway vector pB5nCGW D...
24 1 gi|1050047835|dbj|AB830569.1| Gateway vector pB5GWnY D...
25 1 gi|1050047829|dbj|AB830568.1| Gateway vector pB5GWnR D...
26 1 gi|1050047822|dbj|AB830567.1| Gateway vector pB5GWnG D...
27 1 gi|1050047815|dbj|AB830566.1| Gateway vector pB5GWnC D...
28 1 gi|1050047807|dbj|AB830565.1| Gateway vector pB5GWcY D...
29 1 gi|1050047801|dbj|AB830564.1| Gateway vector pB5GWcR D...
~~~
47 3 gi|2087986960|gb|MW026675.1| Plant expression vector p...
48 3 gi|2033742660|gb|MW557527.1| Cloning vector pGWB-cLUC,...
49 3 gi|2033742656|gb|MW557526.1| Cloning vector pGWB-nLUC,...매우 편안한 것을 알 수 있다. 저거 판다스 데이터프레임으로 묶어서 그룹바이 마렵네
SearchIO iteration
for hit in blast_qresult:
print(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]저런 게 반복된다. 뭐... 리스트나 딕셔너리마냥 iteration, 인덱싱에 슬라이싱도 되고, 특정 영역(id, sequence length)만 볼 수도 있다. 물론 슬라이싱과 동시적용도 된다.
sort_key = lambda hit: hit.seq_len
sorted_qresult = blast_qresult.sort(key=sort_key, reverse=True, in_place=False)
for hit in sorted_qresult:
print("%s %i" % (hit.id, hit.seq_len))gi|1835923315|gb|MN508834.1| 14584
gi|1835923313|gb|MN508833.1| 14548
gi|1693859709|gb|MK204379.1| 13667
gi|2087986978|gb|MW026679.1| 13558
gi|2087986968|gb|MW026677.1| 13439
gi|2033742656|gb|MW557526.1| 13430
gi|2087986973|gb|MW026678.1| 13025
gi|1732472816|gb|MK896897.1| 13004
gi|2087986983|gb|MW026680.1| 12906
gi|1835923311|gb|MN508832.1| 12814
gi|1732472819|gb|MK896898.1| 12754
gi|2033742660|gb|MW557527.1| 12656
gi|1050047859|dbj|AB830573.1| 12546
gi|1050047849|dbj|AB830571.1| 12546
gi|1050047842|dbj|AB830570.1| 12546
gi|1050047835|dbj|AB830569.1| 12531
gi|1050047822|dbj|AB830567.1| 12531
gi|1050047815|dbj|AB830566.1| 12531
gi|1050047854|dbj|AB830572.1| 12519
gi|1050047829|dbj|AB830568.1| 12504
gi|2087986964|gb|MW026676.1| 12332
gi|1050047782|dbj|AB830561.1| 12266
gi|1050047801|dbj|AB830564.1| 12249
gi|1050047720|dbj|AB830552.1| 12228
gi|1050047788|dbj|AB830562.1| 12222
gi|1050047776|dbj|AB830560.1| 12222
gi|1050047770|dbj|AB830559.1| 12213
gi|1050047755|dbj|AB830557.1| 12213
gi|1050047748|dbj|AB830556.1| 12213
gi|1050047807|dbj|AB830565.1| 12207
gi|1050047794|dbj|AB830563.1| 12207
gi|1050047762|dbj|AB830558.1| 12186
gi|1050047734|dbj|AB830554.1| 11931
gi|1050047741|dbj|AB830555.1| 11889
gi|1050047727|dbj|AB830553.1| 11889
gi|2087986960|gb|MW026675.1| 11799
gi|1732472826|gb|MK896900.1| 11186
gi|1732472823|gb|MK896899.1| 11185
gi|1732472812|gb|MK896896.1| 10473
gi|1732472836|gb|MK896904.1| 6122
gi|1755168865|gb|MN055605.1| 5507
gi|1732472841|gb|MK896906.1| 5507
gi|1755168863|gb|MN055604.1| 5506
gi|1732472839|gb|MK896905.1| 5506
gi|1755168869|gb|MN055607.1| 5498
gi|1755168867|gb|MN055606.1| 5497
gi|1755168873|gb|MN055609.1| 5465
gi|1755168871|gb|MN055608.1| 5464
gi|1624318504|gb|MK501803.1| 4844
gi|1732472834|gb|MK896903.1| 4794판다스마냥 특정 요소를 기준으로 정렬할 수도 있다. (시퀀스 길이 기준으로 정렬한 것)
filter_func = lambda hit: len(hit.hsps) > 1
filtered_qresult = blast_qresult.hit_filter(filter_func)
print(filtered_qresult)Query: No (345)
definition line
Target: nt
Hits: ---- ----- ----------------------------------------------------------
# # HSP ID + description
---- ----- ----------------------------------------------------------
0 3 gi|1732472826|gb|MK896900.1| Plant binary expression v...
1 3 gi|1732472823|gb|MK896899.1| Plant binary expression v...
2 3 gi|1732472819|gb|MK896898.1| Plant binary expression v...
3 3 gi|1732472816|gb|MK896897.1| Plant binary expression v...
4 3 gi|1732472812|gb|MK896896.1| Plant binary expression v...
5 3 gi|2087986983|gb|MW026680.1| Plant expression vector p...
6 3 gi|2087986978|gb|MW026679.1| Plant expression vector p...
7 3 gi|2087986973|gb|MW026678.1| Plant expression vector p...
8 3 gi|2087986968|gb|MW026677.1| Plant expression vector p...
9 3 gi|2087986964|gb|MW026676.1| Plant expression vector p...
10 3 gi|2087986960|gb|MW026675.1| Plant expression vector p...
11 3 gi|2033742660|gb|MW557527.1| Cloning vector pGWB-cLUC,...
12 3 gi|2033742656|gb|MW557526.1| Cloning vector pGWB-nLUC,...
Process finished with exit code 0함수 걸어서 필터링도 된다는데 람다가 국룰인건가...
Hit, HSP, Fragment
무슨 족보 그리는것도 아니고 하위의 하위 뭔데...
Hit
아까 hit.id 이런 식으로 잠깐 지나쳤던 것. 어쨌거나 상위 개념인 듯 하다.
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아니 strand 수에 -1은 뭐여...
'Coding > Python' 카테고리의 다른 글
| Biopython으로 Swiss-prot과 ExPASy 데이터베이스 탐방하기 (0) | 2022.08.21 |
|---|---|
| Biopython으로 Entrez database 탐방하기 (0) | 2022.08.21 |
| 번외편-코딩테스트 풀이 (0) | 2022.08.20 |
| Biopython으로 MSA 해보기 (0) | 2022.08.20 |
| Biopython SeqIO 써보기 (0) | 2022.08.20 |
코테 본다는 얘기는 못 들은 것 같고... 그냥 면접보다가 어디 실력 좀 볼까? 해서 나온거긴 한데...
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)비효율적으로 보일 지 몰라도 이정도면 본인 수준에서는 장족의 발전임다... 어레이 제대로 된 것도 확인했음.
record_table=pd.DataFrame({"ID":record_id_array,"Length":record_len_array})
print(record_table)아 표 ㄹㅇ 이쁘게 나왔음.
ID Length
0 SRR000021.1.2 176
1 SRR000021.1.4 49
2 SRR000021.2.2 99
3 SRR000021.2.4 106
4 SRR000021.3.2 64
.. ... ...
102 SRR000021.52.2 269
103 SRR000021.52.4 0
104 SRR000021.53.2 276
105 SRR000021.54.2 65
106 SRR000021.54.4 184캬... 이 표 이쁜거 보게. 근데 여기서 끝이 아니고 length를 또 봐야 한다... 살려주세요. 일단 전에 배웠던대로 오름차순 정렬을 해 보자.
record_table2=record_table.sort_values('Length',ascending=0) ID Length
68 SRR000021.35.2 280
104 SRR000021.53.2 276
102 SRR000021.52.2 269
64 SRR000021.33.2 269
70 SRR000021.36.2 268
.. ... ...
63 SRR000021.32.4 0
73 SRR000021.37.4 0
65 SRR000021.33.4 0
71 SRR000021.36.4 0
53 SRR000021.27.4 0엥? 중복값이 있구나 이거.
record_table2=record_table.groupby('Length').count()
record_table2=record_table2.sort_values('Length',ascending=0)
print(record_table2.head(10)) ID
Length
280 1
276 1
269 2
268 1
259 1
257 1
249 1
247 1
244 1
240 1ID가 왜 저렇게 나오냐면 저거 그룹바이로 같은 값 묶어서 세달라고 카운트 걸어서 그렇다. 즉, 269에 2 있는 건 269bp가 두개라는 얘기. 그럼 아이디는 어떻게 뽑나요 쿼리 걸어 쿼리
근데 생각해보니까 이거 빈도수 Top10이지 길이가 아니다.
record_table2=record_table2.sort_values('ID',ascending=0) ID
Length
0 15
107 3
63 3
99 3
44 2
105 2
111 2
178 2
64 2
124 2이거 뭐 하는 데이터인데 0개가 이래 많음?
5. 조건부 합계(여기서는 100bp 기준으로 나눴다)
for record2 in SeqIO.parse(handle2,"fasta"):
if len(record2.seq) > 100:
i=0
i=i+len(record2.seq)
else:
j=0
j=j+len(record2.seq)
print(i,j)얘도 각개로 나온다. 일단 2번부터 해결해야 할 듯. 참고로 4번 코드를 보시면 아시겠지만, 넘파이 판다스 둘 다 불렀다.
굳이 if 하지 말고 저거 해결 봤으니까 걍 위에 표 불러서 판다스 소환하자.
print(record_table[record_table.Length >= 100].sum())
print(record_table[record_table.Length < 100].sum())ID SRR000021.1.2SRR000021.2.4SRR000021.3.4SRR0000...
Length 9984
dtype: object
ID SRR000021.1.4SRR000021.2.2SRR000021.3.2SRR0000...
Length 1955
dtype: object이거 셀렉터 못줌?
print("length >= 100bp:",record_table['Length'][record_table.Length >= 100].sum())
print("length < 100bp:",record_table['Length'][record_table.Length < 100].sum())length >= 100bp: 9984
length < 100bp: 1955되는데요?
'Coding > Python' 카테고리의 다른 글
| Biopython으로 Entrez database 탐방하기 (0) | 2022.08.21 |
|---|---|
| Biopython으로 BLAST 돌려보기 (0) | 2022.08.20 |
| Biopython으로 MSA 해보기 (0) | 2022.08.20 |
| Biopython SeqIO 써보기 (0) | 2022.08.20 |
| Biopython으로 시퀀스 레코드 생성하고 만져보기 (0) | 2022.08.20 |