[하루에 한 R] pair 형태의 데이터를 matrix로 전환할 때 주의할 점

오늘 쓰려는 주제는 R 활용 기법과 직접적인 관련은 없음을 미리 밝힌다.

한 쌍의 유전체 사이에 어떤 계산을 하여 수치를 얻었다고 가정하자. 기본적으로 얻어지는 결과물은 다음과 같이 pair 형태일 것이다.

친절한 프로그램이라면 matrix 형태로 전환한 파일을 제공해 줄 것이고, 이를 R에서 데이터 프레임으로 읽어들여서 여러 가지의 후속 분석을 실시할 수 있다. 만약 pairwise data만 위 그림과 같이 제공한다면 R에서 reshape 패키지를 이용하여 matrix 형태로 전환하면 된다.

그런데 여기에서 주의할 점이 하나 있다. 바로 자기 자신에 대한 데이터를 포함하고 있느냐를 고려해야 한다. 위에서 보인 그림 중 첫번째 것은 dRep의 결과물이다. FASTA 파일 여러개를 제공하면 프로그램이 알아서 모든 쌍에 대한 ANI를 계산하는데, 이때 자기 자신에 대한 것도 포함한다. 계산값을 보여주는 첫 줄이 그러한 것으로서 alignment coverage와 ani 전부 1,0으로 나타난다.

이것이 전부가 아니다. A->B, B->A 각 경우에 대한 값이 전부 존재하는가? 그렇다면 동일한가? 결과가 없는 것은 어떻게 할 것인가? 의외로 생각할 것이 많다. 그리고 각 경우에 대해서 처리하는 R 코드가 조금씩 달라진다.

위에 보인 그림에서 아래의 것은 POCP.sh를 사용하여 percentage of conserved protein을 계산한 것이다. 이 스크립트에서는 유전체 서열 한 쌍과 스레드 수만 인수로 제공해야 한다. 따라서 query와 subject 서열의 쌍을 별도로 계산하여 넣어 주어야 한다. 나는 Perl의 Algorithm::Combinatorics 모듈을 사용했으므로 당연히 자기 자신에 대한 계산을 하지 않게 만들었다(샘플 코드는 여기를 참조). 계산을 하나마나 POCP = 100%가 나올 것이기 때문이다. 하지만 이를 고려하지 않고 R의 reshape 패키지를 그냥 돌리면 매트릭스 형태가 좀 이상해진다. 전부 1(=100%)로 채워지는 대각선에 대한 배려가 없기 때문이다. 10개의 유전체 서열을 전부 비교하여도 R로 읽어들어 전환한 다음에는 10x10의 매트릭스에서 하나가 부족해진다. query = subject에 해당하는 대각선도 제대로 존재하지 않는다.

이 문제를 해결하려면 pairwise 결과 파일의 뒷부분에 자기 자신에 대한 분석값을 넣는 것이 좋다. A, B, C,...J의 10개 유전체에 대해서 POCP를 계산한다면 다음의 10줄을 넣기만 하면 된다. 자기 자신에 대한 POCP 값은 당연히 100%이므로 이를 세번째 컬럼에 입력해 넣는다.

A A 100.0
B B 100.0
C C 100.0
..
J J 100.0

이렇게 하여 reshape의 cast() 함수를 돌려 처리를 하면 10x10의 매트릭스, 정확히 말하자면 데이터 프레임(df라 하자)이 생길 것이다. A->B, B->A 두 경우에 대해서 한쪽으로만 계산을 하였으니 데이터 프레임으로 전환된 뒤의 구조를 확인해 보면 대각선은 전부 100이고 나머지 절만만 채워진 상태일 것이다.

NA를 전부 0으로 치환한 다음 df와 t(df)를 더하면 대각선 절반에 대한 값이 채워진다. 그러나 대각선은 전부 200이 될 것이다. 이는 또 적절히 처리하면 된다.

오늘의 교훈은 무엇인가? 데이터 프레임을 한번쯤은 엑셀 등에서 열어보고 어떤 상태인지를 눈으로 확인하는 것이 좋다.

2019년 3월 14일 업데이트

Algorithm::Combinatorics의 combination_with_repetition() 서브루틴을 사용하면 자기 자신에 대한 검색 결과도 만들 수 있다. 이것을 사용하여 POPC.sh를 실행하는 Perl 코드는 다음과 같다.

#!/usr/bin/perl

use strict; use warnings;
use Algorithm::Combinatorics qw(combinations_with_repetition);

open LIST, 'list';
my @temp = ();
while (<LIST>) {
    chomp;
    push @temp, $_;
}
my $faa_ref = \@temp;

my $iter = combinations_with_repetition($faa_ref, 2);

while (my $c = $iter->next) {
    my $out = 'pocp' . '_' . $c->[0] . '_vs_' . $c->[1];
    print "Running for $out...\n";
    system("POCP.sh $c->[0] $c->[1] 24 > $out");
}