배경

돼지생식기 호흡기 증후군 (Porcine reproductive and respiratory syndrom)은 PRRS 바이러스에 의한 질병으로 자돈이나 육성돈이 이 병에 걸리면 기침, 호흡곤란, 폐렴 등 호흡기 증상을 보이며, 모돈의 경우 임신말기에 유사산 및 조산을 나타내고 허약자돈을 분만하게 되고, 웅돈에서는 정액성상의 이상을 일으키는 등 병명그대로 번식장애와 호흡기 증상으로 인한 성장부진이 특징입니다. 1986년 미국에서 보고되고 1990년 유럽에서 처음 보고된 비교적 새로운 가축질병이며 병원체가 알려지기 전에는 Mystery pig disease라고 불리기도 했습니다.

구제역이나 조류독감처럼 축산산업에는 큰 영향을 주는 질병으로 현재 미국에서는 PADRAP이라는 위해성 평가 프로그램을 구성하여 양돈장이나 종돈장을 체계적으로 관리하기도 합니다.

이번 포스트는 실험에서 얻은 PRRSV 유전자 서열을 기존 서열과 비교해서 PRRSV의 고병원성 여부를 확인하는 Perl 스크립트를 만드는 방법을 살펴보겠습니다. 스크립트를 구현하는 초보자 입장에서 문제를 해결하는 과정을 담아 보려고 했습니다.

고병원성 PRRSV의 특징

Zhou 등에 따르면 고병원성 PRRSV 균주는 Nsp2 유전자 분절의 코딩영역에 고유의 30잔기의 아미노산 결실이 있는 것으로 확인되었습니다. 여기서 Nsp2는 nonstructural protein 2의 약어이며, 중간 부위에 유전자 변이가 매우 많지만, 프로테아제 도메인으로 예상되는 N말단과 트랜스멤브레인으로 예측된 C말단은 보존성이 매우 높은 구조를 가지고 있습니다(Han et al. 2007).


Figure 1. The 30-Amino-Acid Deletion in the Nsp2 of Highly Pathogenic PRRSV. (Zhou et al. 2009)

스크립트 작동조건

실험에서 직접 시퀀싱한 서열을 고병원성인지 확인하는 것이므로 입력 서열은 DNA이며, 이 예제에서의 기준서열은 Figure 1에서 제시된 것처럼 단백질 서열로 가정합니다(VR-2332). 그리고 결실이 있는 위치와 크기의 허용오차를 줄 수 있어야 합니다. 입력인 DNA 서열은 Nsp2의 특정 부위를 증폭한 PCR산물을 직접 시퀀싱한 것입니다.

해결책을 찾아가는...

일단 이 문제는 쌍서열정열(pairwise sequence alignment) 전형적인 예입니다. 서열정열만 하면 거의 해결에 가까이 온 것입니다. 먼저 간단한 퀴즈 하나를 내보겠습니다. 이 문제의 경우 전역정열(global alignment)와 지역정열(local alignment) 알고리즘 중에서 무엇을 선택하는 것이 좋을까요? 이 퀴즈에 답을 하기 위해서는 기준이 되는 단백질 서열의 길이와 입력이 되는 DNA 서열의 길이를 비교해야 합니다. 기준 서열은 PRRSV의 완전유전체 중의 하나이며, 입력 서열은 Nsp2의 PCR 산물이므로 그 길이가 많이 다릅니다. 따라서, 지역정열을 해야 적절한 결과를 얻을 수 있습니다.

많은 지역정열 프로그램이 있지만, 생물정보학의 기본중의 기본이라고 하는 BLAST의 패키지중의 하나인 bl2seq를 이용해서 정열을 구현할 수 있습니다. 여기서 두 번째 퀴즈... 어떤 BLAST 프로그램을 지정해야 할까요? 쿼리가 DNA이고, 서브젝트가 단백질서열이면.... 바로 BLASTX를 이용해야 합니다. 당연히 아시겠지만, bl2seq는 별도의 formatdb없이 바로 두 시퀀스를 지정하면 됩니다.

BLAST를 하고 나면 생물정보 스크립트 단골 메뉴인 파서 구현이 나오겠지만, 여기서는 일단 BioPerl이라는 걸출한 라이브러리를 이용합니다. 이미 잘 짜여 있는데 새로 만들 필요는 없으니까요. 따라서 직접 구현해야하는 스크립트의 핵심 기능은 사용자가 지정하는 오차범위에 일치하는 긴 결실이 있는지를 판단하는 부분입니다.

먼저 bl2seq의 결과 파일을 보면 Figure 2와 같이, 쿼리 요약, HSP (high scoring pair), 파라메터, 통계량 등을 보여줍니다. 

Figure 2. An example BLAST's result of the deletion in Nsp2.

BioPerl의 Bio::SearchIO 모듈을 이용한 BLAST, FASTA, HMMER, Sim4 등의 다양한 정열 프로그램 결과를 읽거나 저장할 수 있습니다. 바로 이 모듈을 이용해서 bl2seq 결과를 파싱(parsing)할 수 있습니다. 구체적인 사용법은 BioPerl 문서의 매뉴얼이나 HOWTO 문서를 읽어보시길 권합니다.

위 의 방법이든 직접 파싱을 하든 일단 HSP에서 긴 결실 부위를 찾는 것이 핵심인데, 연구를 열심히 하는 생물정보 스크립터들이 잘 빠지는 오류중의 하나가 너무 '현학적인'인 방법을 생각하는 것입니다. 이번의 경우에도 입력 서열의 품질이나 정제 상태에 따라서 지역정열의 결과 패턴이 매우 다양하고, 이를 고려하다 보면 각 컬럼별 통계량을 계산하거나, 프레임시프트 여부, 연속 결실을 계산하기 위한 그래프 알고리즘... 이런 상태에 봉착하셨다면 커피한잔이나 산책 후에 "Perl script should be Perlish"를 상기하면 어떨까요?

다음 Perl 코드는 HSP에서 "-"로 연속되는 결실을 포함하고 있는 문자열을 가져오는 것입니다.

my $query_alignment = $hsp->query_string();

이제 Perl의 초강력 정규식을 이용해서 연속하는 결실부위를 찾습니다. 아래 코드는 최소 크기 이상을 가지는 연속되는 결실을 모두 찾는 코드입니다.

my ($ref_start, $ref_end) = $hsp->hit->strand < 0 ?
                                    ( $hsp->hit->end,
                                      $hsp->hit->start ) :
                                    ( $hsp->hit->start,
                                      $hsp->hit->end );

while( $query_alignment =~ /(\-{$threshold_length,})/g ) {
      $indel_start = $-[1] + $ref_start;
      $indel_end = $-[1] + $ref_start + length($1) - 1;
      $indel_length = length($1);
      push  @indels, [$indel_start, $indel_end];
} 

이렇게 @indel 배열에 결실에 대한 정보를 저장하고, threshold % indentification (-p), expected start position of indel (-s), expected end position of indel (-e), tolerant positional error(-d)와 같은 추가적인 필터 조건을 구현하고 실제로 동작시키면 아래와 같습니다.

[bckang@gxs ~]$ ./gvs_find_indel.pl nsp2.blsx -s 921 -e 949 -d 15
#Mutation dectector's parameters: p=30 s=921 e=949 d=15
>AAO13191(NA_prototype)
920     948     29
>AAO13191(NA_prototype)
>AAO13191(NA_prototype)
>AAO13191(NA_prototype)
>AAO13191(NA_prototype)
>AAO13191(NA_prototype)

이 문제를 해결하는 방법은 매우 다양합니다. 위의 방법은 이미 결실 위치와 크기가 알려진 경우에는 간단한 해법이 될 것이라고 생각됩니다. 마지막으로 스크립트를 만드시는 분들께

Perl is supposed to be Perlish and Python looks like Pythonic.

References

1. Zhou, L. et al. (2009) The 30-amino-acid deletion in the Nsp2 of highly pathogenic porcine reproductive and respiratory syndrome virus emerging in China is not related to its virulence. J. Virol. 83(10):5156-67.
2. Han, J. et al. (2007) Identification of Nonessential Regions of the Nsp2 Replicase Protein of PRRSV Strain VR-2332 for Replication in Cell Culture. J. Virol. 81(18):9878-90.
3. Wikipedia: PRRSV.
4. Altscul, S.F. et al. (1999) Basic local alignment search tool. J. Mol. Biol. 215:403-410.
5. Thompson, J.D. et al. (1994) CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 22(22):4673-4680.

6. Perl, http://www.perl.org

7. BioPerl, http://www.bioperl.org
   

   KM사업부장 강병철
   
   
   

Posted by 人Co