종 내에서 두개씩 쌍을 지어 계산한 pairwise similarity를 평균과 오차를 낸 것입니다.
사례로 든 자료를 찾아보니까, 방법이 다음처럼 정리되어 있었습니다.
"Vibrio의 제놈 정보로부터 16S rRNA 유전자 염기서열을 추출 하여 개체 내에 존재하는 16S rRNA 유전자 이질성을 조사하였다. 유전체내 16S rRNA 유전자 염기서열을 비교하기 위하여, 추출한 16S rRNA 염기서열을 CLUSTAL W 1.8 (21)을 이용하 여 하나의 서열로 배열하였다. 이후 BioEdit 5.0.6 (10) 프로그램 을 이용하여, 각각의 염기서열 위치에 대한 엔트로피(entropy) 값 을 계산하였다. 또한 유전체내 16S rRNA 유전자 염기서열 상호 간에 DNA 상동성(similarity)값을 계산하여, 개체 내의 16S rRNA 유전자 상호간의 유사성을 분석하였다" (기장서, 2009).
제가 언듯 느끼기엔 방법에 대한 설명이 많이 부족합니다.
Vc에 8개 샘플이라면 먼저 다중정렬을 합니다. 논문에서는 CLUSTALW 프로그램을 썼습니다. 하지만 요즘은 MUSCLE이나 CLUSTAL-OMEGA를 더 사용하겠죠.
다음으로 pairwise similarity를 계산하는데. 8C2=28쌍의 값이 나오겠죠. Pairwise similarity도 여러가지 식이 있으니까 잘 골라 써야 할텐데(보통은 K2P인 것 같습니다만, 어떤 곳에는 unweighted pairwise similarity를 쓴 곳도 있습니다) 이 논문에는 그 얘기가 없지만, 아마 K2P를 적용했을 겁니다. Blast같은 local alignment 프로그램을 쓰면 안됩니다 (안되는 건 아니지만, 정렬된 구간에 주의해서 사용해야 함. gap을 어떻게 처리하는지도 표준을 따라야 함).
다음으로 평균과 표준편차를 계산합니다.
예를 들어서, 아래의 표는 pairwise similarity를 반 행렬로 정리한 것입니다. 여기서 자기 vs 자기의 100%인 경우를 빼고 나머지를 평균낸 것입니다.
출처: Yoshizawa, Susumu & Wada, Minoru & Yokota, Akira & Kogure, Kazuhiro. (2010). Vibrio sagamiensis sp. nov., luminous marine bacteria isolated from sea water. The Journal of general and applied microbiology. 56. 499-507. 10.2323/jgam.56.499.
제시된 논문은 방법 기술이 부실합니다. 딱 그 해당분야에 제시된 방법론이 있는 것으로 알고 있습니다. 정확히 하기 위해서는 조금 뒤져보셔야 해요.
질문자가 이런 종류의 프로그램 사용이 어렵다면, 통합적인 툴을 사용하는 것이 좋습니다. 무료이면서 기능도 여러가지 잘 갖춰져 있는 MEGA나 DNASP, BioEdit를 접해보세요. 이중에서 서열도 직접 보면서 작업하려면 BioEdit가 나을 겁니다. 여기서 서열 편집, ClustalW 정렬, Distance matrix 계산이 다 가능합니다. 예를 들어서, 거리 행렬의 경우 이런 식으로 나옵니다.
A01_prin_A 0.0000 0.0001 0.0001 0.0001 0.0004 0.0004 0.0005 0.0004
0.0004 0.0011 0.0012 0.0012 0.0013 0.0016 0.0020 0.0021 0.0020
0.0020 0.0022 0.0022 0.0022 0.0018 0.0018 0.0015 0.0018 0.0022
0.0023 0.0025 0.0024 0.0023 0.0018 0.0021 0.0287
A02_argy_A 0.0001 0.0000 0.0000 0.0000 0.0003 0.0003 0.0005 0.0003
0.0004 0.0011 0.0011 0.0012 0.0012 0.0015 0.0019 0.0020 0.0020
0.0019 0.0022 0.0022 0.0022 0.0017 0.0017 0.0015 0.0017 0.0021
0.0023 0.0024 0.0023 0.0022 0.0017 0.0020 0.0286
A03_argy_A 0.0001 0.0000 0.0000 0.0000 0.0003 0.0003 0.0005 0.0003
0.0004 0.0011 0.0011 0.0012 0.0012 0.0015 0.0019 0.0020 0.0020
0.0019 0.0022 0.0022 0.0022 0.0017 0.0017 0.0015 0.0017 0.0021
0.0023 0.0024 0.0023 0.0022 0.0017 0.0020 0.0286
A04_argy_A 0.0001 0.0000 0.0000 0.0000 0.0003 0.0003 0.0005 0.0003
0.0004 0.0011 0.0011 0.0012 0.0012 0.0015 0.0019 0.0020 0.0020
0.0019 0.0022 0.0022 0.0022 0.0017 0.0017 0.0015 0.0017 0.0021
0.0023 0.0024 0.0023 0.0022 0.0017 0.0020 0.0286
A08_stol_A 0.0004 0.0003 0.0003 0.0003 0.0000 0.0004 0.0006 0.0004
0.0004 0.0012 0.0012 0.0012 0.0013 0.0016 0.0020 0.0021 0.0020
0.0020 0.0023 0.0023 0.0023 0.0018 0.0018 0.0015 0.0018 0.0022
0.0023 0.0025 0.0024 0.0023 0.0018 0.0021 0.0287
A05_rubr_A 0.0004 0.0003 0.0003 0.0003 0.0004 0.0000 0.0006 0.0004
0.0005 0.0012 0.0013 0.0013 0.0014 0.0016 0.0021 0.0021 0.0021
0.0021 0.0023 0.0023 0.0023 0.0019 0.0018 0.0016 0.0018 0.0022
0.0023 0.0025 0.0024 0.0023 0.0018 0.0022 0.0287
엑셀에 가져다가 샘플마다 평균과 편차만 구해서 바그래프를 그리면 되겠죠.
-
레벨2
glow216
(대학생)
-
20.07.14 13:25
친절한 답변 감사합니다.
말씀하신대로 pairwise similarity를 반 행렬로 만드는 과정을 보여주신 논문(Vibrio sagamiensis sp. nov., luminous marine bacteria isolated from sea water.)을 참고하여 Vibrio Cholerae의 16S rRNA를 대상으로 진행해보았습니다.
MEGA X를 이용하면
DISTANCE > Compute Pairwise Distances > Model/Method : Kimura 2-parameter model
으로 아래와 같은 결과를 얻을 수 있었습니다.
0.0000000000
0.0012936618 0.0012936618
0.0012928255 0.0012928255 0.0006464125
0.0077923185 0.0077923185 0.0077973819 0.0077923185
0.0000000000 0.0000000000 0.0019357978 0.0012928255 0.0077923185
0.0006459949 0.0006459949 0.0006464125 0.0019404940 0.0084467689 0.0006459949
0.0000000000 0.0000000000 0.0012936618 0.0012928255 0.0077923185 0.0000000000 0.0006459949
그리고
BioEdit를 이용하면 Kimura 2-parameter model를 특별히 설정하는 단계가 없어서 (또는 몰라서)
Alignment > Sequence Identity Matrix
으로 아래와 같은 결과를 얻을 수 있었습니다.
ID 0.999 0.994 0.998 0.990 0.996 0.998 0.999
0.999 ID 0.995 0.998 0.991 0.997 0.999 1.000
0.994 0.995 ID 0.996 0.988 0.997 0.996 0.995
0.998 0.998 0.996 ID 0.991 0.996 0.998 0.998
0.990 0.991 0.988 0.991 ID 0.989 0.990 0.991
0.996 0.997 0.997 0.996 0.989 ID 0.996 0.997
0.998 0.999 0.996 0.998 0.990 0.996 ID 0.999
0.999 1.000 0.995 0.998 0.991 0.997 0.999 ID
두 과정 모두 ClustalW 과정을 거쳤습니다.
위 과정이 올바르게 진행된 것인지,
두 결과가 왜 다른지,
BioEdit를 이용하면 Kimura 2-parameter model 설정이 없는지
궁금합니다.
그리고 엑셀에 샘플마다 평균과 편차를 구한다는 말씀이 행마다 구하는 것인지 잘 이해되지 않습니다.
감사합니다!
Identity(%)는 두개씩의 쌍별 비교에서 같은 염기의 비율로 Similarity와는 비슷하지만 다른 개념입니다. 그 성질을 K2P 계산식에서 생각해 보세요.
MEGA에서는 K2P를 distance로 계산한 값입니다. Similarity는 1-distance로 환산하면 됩니다.
BioEdit에서는 DNADIST ---> Neighborhood phylogenetic tree 에서 나옵니다. 원래 DNADIST 프로그램에서 골라 설정할 수 있는데, K2P가 기본옵션이니 그냥 사용하면 됩니다. 아래는 해당 프로그램의 사용자문서 중 메뉴에 관한 부분입니다.
Nucleic acid sequence Distance Matrix program, version 3.5c
Settings for this run:
D Distance (Kimura, Jin/Nei, ML, J-C)? Kimura 2-parameter
T Transition/transversion ratio? 2.0
C One category of substitution rates? Yes
L Form of distance matrix? Square
M Analyze multiple data sets? No
I Input sequences interleaved? Yes
0 Terminal type (IBM PC, VT52, ANSI)? ANSI
1 Print out the data at start of run No
2 Print indications of progress of run Yes
-
레벨2
glow216
(대학생)
-
20.07.15 21:08
감사합니다.
마지막으로 엑셀에 샘플마다 평균과 편차를 구한다는 말씀이 반행열 값 전체에 대해 구하는 것이 맞는지 궁금합니다.
그리고 Mega X와 BioEdit의 두 결과가 다른 이유가 궁금합니다!
제가 구한 Vc의 결과입니다.
평균 표준편차
0.9974281740 0.00311833381
0.9974392857 0.00312368135
BioEdit와 MEGA에서 결과가 약간 다른 것은 아마 내부적으로 유효숫자 자리수를 처리하는 규칙(변수형)이 달라서 일겁니다. 예시로 든 숫자가 다른가, 같은가는 유효숫자로 판단해야 합니다. 예를 들어서, 100개의 염기에서 15개가 다르면 identity로는 85%가 될 텐데, 이것을 0.85라고 할 지, 0.85000000 이라고 할지 판단할 수 있겠죠? 만약 111개의 염기에서 15개라면, 0.135135135 긴 자릿수의 소수가 되는데, 이걸 어느 자리에서 끊어야 정확한 의미가 되는지 프로그램마다의 기준이 약간씩 다른 듯 합니다. 내부적으로 이러한 부정형의 숫자들을 곱하고 나누고 하다보면 그 오차가 점점 누적되어서 최종의 결과는 상당히 달라보이는 숫자가 나오게 됩니다. 물론 우리는 수학적으로 쉽게 판단할 수 있어요. 두 프로그램의 결과는 유효숫자 범위 안에서는 충분히 정확하기 때문에, 자릿수만 제한시켜서 정리하면 됩니다. 제 생각에, 세자리 수의 정렬 길이라면 소숫점 네자리까지만 반올림으로 표시해야 할 것 같습니다.
-
레벨2
glow216
(대학생)
-
20.07.16 22:14
항상 친절한 답변 감사합니다.
위에서 말씀드린대로,
논문에서는 V. cholerae 99.7% (SD=0.36, N=20)라는 결과가 나왔고,
저는 평균=0.997, 표준편차=0.0031라는 결과가 나왔습니다.
V. cholerae 99.7%는 평균=0.997와 일치하지만, (SD=0.36, N=20)가 어떻게 구해진 것인지 궁금해서 질문드립니다.
SD=0.36가 Standard Deviation이라면 표준편차=0.0031와 다르게 나온 것으로 보이고, 또한, N=20는 무엇을 의미하는지 궁금합니다.
(평균=0.997, 표준편차=0.0031는 좌행렬 값 전체에 대한 평균과 표준편차를 구한 결과입니다.)
학부생 맞아요? 집요하네요!!
해당 논문에서 오류를 여럿 찾았습니다. 그러니 궁금한 부분에 저자가 세심하게 손보지 못한 오류가 있을 수 있다는 걸 감안해야 할 듯 합니다.
N과 관련해서, 16S rRNA gene이 8개인 Vc와 Vs는 이유를 도통 모르겠습니다. 오류가 아닐까 합니다. 11개인 Vh(Table 1)는 N=11*10/2=55개가 되어야 겠지만, 제가 찾아보니까 2번 염색체에 있는 1개가 pseudogene 이었습니다. 그래서 (틀림없이!) 저자가 제외시켜 최종적으로 N=10*9/2=45개가 된 겁니다.
그리고, 본문에 안나와 있고, 또 추측할 수 있는 방법적 가능성이 하나가 아닌 경우 (예를 들면, Similarity 계산방법, 정렬에서 InDel 처리방법, pseudogene처리 등등), 본문에 제시된 값과 내가 손에 들고 있는 값이 왜 다른지 알아내는 일은 고되고, 무의미한 일입니다. 이 분야에서 Similarity 구해서 종을 동정하거나 유전다양성을 분석하거나 하는 방법적 표준을 찾아서 이 전쟁터로 다시 돌아오는 것이 백번 나은 전략이겠죠.
아래 논문을 보니, 깨달아지는게 많았습니다. 비전문가가 괜히 나섰구나 후회도 조금 되구요. 이걸 읽다보니까 오래된 논문이라면 p-distance (uncorrected distance)나 % identity를 썼을 가능성이 있겠구나 생각이 듭니다. 저자가 쓴 다른 논문을 보면, 아마 자세한 방법이 나와 있는 게 있을테니 뒤져 보세요.
https://www.nature.com/articles/npjbiofilms20164
https://help.ezbiocloud.net/how-to-calculate-16s-rrna-sequence-similarity-values-for-bacterial-taxonomy-why-blast-should-be-avoided/
-
레벨2
glow216
(대학생)
-
20.07.23 21:34
친절한 답변 감사합니다.
말씀하신 "N=11*10/2=55"와 "N=10*9/2=45"가
N = 염색체의 개수 * (염색체의 개수-1) / 2
로 구하는지((염색체의 개수-1)가 맞는지?), N이 무엇을 의미하는지 궁금합니다!
또한, 말씀하신대로 방법적 표준을 찾아 결과를 도출하는 것으로
앞서서 진행한
MEGA의 K2P를 1-distance로 계산한 Similarity와
BioEdit의 DNADIST ---> Neighborhood phylogenetic tree를 이용해서 구한
평균 표준편차
0.9974281740 0.00311833381
0.9974392857 0.00312368135
의 과정들이 맞아서, 이대로 진행해도 되는지 궁금합니다!
(평균과 표준편차는 좌행렬 전체에 대해 계산하였습니다.)
그리고 Box Plot에서
박스의 값이 평균이고, 직선으로 표시된 값이 표준편차가 맞는지 궁금합니다!
감사합니다.
###
말씀하신 "N=11*10/2=55"와 "N=10*9/2=45"가
N = 염색체의 개수 * (염색체의 개수-1) / 2
로 구하는지((염색체의 개수-1)가 맞는지?), N이 무엇을 의미하는지 궁금합니다!
###
A: N은 쌍별비교를 하는 pair의 수입니다. '염색체'의 개수가 아니라 '16S 서열'의 갯수이고요. Vh의 16S가 10개 있기 때문에 비교쌍의 수는 10*10/2이지만, 자신과의 비교쌍은 빼야 하기 때문에 10*(10-1)/2=45가 된 겁니다.
###
또한, 말씀하신대로 방법적 표준을 찾아 결과를 도출하는 것으로
앞서서 진행한
MEGA의 K2P를 1-distance로 계산한 Similarity와
BioEdit의 DNADIST ---> Neighborhood phylogenetic tree를 이용해서 구한
평균 표준편차
0.9974281740 0.00311833381
0.9974392857 0.00312368135
의 과정들이 맞아서, 이대로 진행해도 되는지 궁금합니다!
(평균과 표준편차는 좌행렬 전체에 대해 계산하였습니다.)
###
A. 네. 둘 다 K2P로 계산한 값이니까 같은 거죠. 하지만 유효숫자를 따져서 소숫점 이하 자릿수를 끊어야 합니다. 16S 길이가 1kb 정도라면, 유효자리는 4개이기 때문에, 소숫점이하 3자리까지만 넣으면 됩니다.
0.997, 0.003
0.997, 0.003 이 되겠죠.
(그런데, 저도 오기가 생겨서 여러 방법으로 계산을 해 봤는데, 논문의 값과 조금씩 달리 나오네요. K2P, p-distance, J-C 와도 다 달라요. 저도 궁금합니다. 그렇지만, 오차 범위에서는 값이 다 같기 때문에, 오늘은 일단 넘어가도록 합시다).
###
그리고 Box Plot에서
박스의 값이 평균이고, 직선으로 표시된 값이 표준편차가 맞는지 궁금합니다!
###
맞습니다. Box plot은 box-and-whisker plot, box-and-whisker diagram이고, 저것은 막대그래프(bar graph/chart)라고 합니다.